Hello! I am trying to quantify expression of NMD mutations in paired end WGS and RNAseq data. I am trying to use a custom genome file because when I use the GATK.GRCh38 I lose the insertions that are present in my WGS. I first ran the Sarek and then I used the gatk FastaAlternateReferenceMaker with the deepvariant vcf file (from Sarek) and the Ensembl: Homo_sapiens.GRCh38.dna.primary_assembly.fa. I have tried to use the igenomes Enseml gtf: /fdb/igenomes_nf/Homo_sapiens/Ensembl/pub/release-110/gtf/Homo_sapiens.GRCh38.110.gtf, I have tried using: gffread /fdb/igenomes_nf/Homo_sapiens/Ensembl/pub/release-110/gtf/Homo_sapiens.GRCh38.110.gtf -g custom_fasta.fasta -o sample.custom.gtf, and g2gtools vcf2chain and g2gtools transform to attempt to liftover to a custom GTF file. It seems no matter the gtf file used I get an error:

Parsing gtf File is done!

"rsem-extract-reference-transcripts rsem/genome 0 Homo_sapiens.GRCh38.110.gtf None 0 rsem/custom_fasta.fasta" failed! Plase check if you provide correct parameters/options for the pipeline!

Here is the code sbatch scripts I have tried all with the three different gtf files stated above. If anyone has any insight on how to properly create the gtf file for a custom genome please let me! Thanks so much for your help!

#! /bin/bash
#SBATCH --job-name=nextflow-main
#SBATCH --cpus-per-task=16
#SBATCH --mem=200G
#SBATCH --gres=lscratch:200
#SBATCH --time=40:00:00


module load nextflow
export NXF_SINGULARITY_CACHEDIR=/data/$USER/nxf_singularity_cache;
export SINGULARITY_CACHEDIR=/data/$USER/.singularity;
export TMPDIR=/lscratch/$SLURM_JOB_ID
export NXF_JVM_ARGS="-Xms2g -Xmx4g"

nextflow run nf-core/rnaseq \
-profile biowulflocal \
--input samplesheet_RNA.csv \
--outdir rnaseq_out \
--fasta /data/LBD/1000genomes/SRA_downloads/RNA_Test/Ref/Sample_HG00176_.nfcore.fasta \
--gtf /fdb/igenomes_nf/Homo_sapiens/Ensembl/pub/release-110/gtf/Homo_sapiens.GRCh38.110.gtf \
-resume

#! /bin/bash
#SBATCH --job-name=nextflow-main
#SBATCH --cpus-per-task=16
#SBATCH --mem=200G
#SBATCH --gres=lscratch:200
#SBATCH --time=40:00:00


module load nextflow
export NXF_SINGULARITY_CACHEDIR=/data/$USER/nxf_singularity_cache;
export SINGULARITY_CACHEDIR=/data/$USER/.singularity;
export TMPDIR=/lscratch/$SLURM_JOB_ID
export NXF_JVM_ARGS="-Xms2g -Xmx4g"

nextflow run nf-core/rnaseq \
-profile biowulflocal \
--input samplesheet_RNA.csv \
--outdir rnaseq_out \
--fasta /data/LBD/1000genomes/SRA_downloads/RNA_Test/Ref/Sample_HG00176_.nfcore.fasta \
--gtf /fdb/igenomes_nf/Homo_sapiens/Ensembl/pub/release-110/gtf/Homo_sapiens.GRCh38.110.gtf \
--aligner hisat2 \
--skip_gtf_filter \
--skip_gtf_transcript_filter \
-resume

#! /bin/bash
#SBATCH --job-name=nextflow-main
#SBATCH --cpus-per-task=16
#SBATCH --mem=200G
#SBATCH --gres=lscratch:200
#SBATCH --time=40:00:00


module load nextflow
export NXF_SINGULARITY_CACHEDIR=/data/$USER/nxf_singularity_cache;
export SINGULARITY_CACHEDIR=/data/$USER/.singularity;
export TMPDIR=/lscratch/$SLURM_JOB_ID
export NXF_JVM_ARGS="-Xms2g -Xmx4g"

nextflow run nf-core/rnaseq \
-profile biowulflocal \
--input samplesheet_RNA.csv \
--outdir rnaseq_out \
--fasta /data/LBD/1000genomes/SRA_downloads/RNA_Test/Ref/Sample_HG00176_.nfcore.fasta \
--gtf /fdb/igenomes_nf/Homo_sapiens/Ensembl/pub/release-110/gtf/Homo_sapiens.GRCh38.110.gtf \
--skip_gtf_filter \
--skip_gtf_transcript_filter \
-resume

RSEM error is probably due to the GTF coordinates pointing to incorrect positions in your custom fasta. You'll need to figure out a way to correct the coordinates in your GTF file to match your fasta file. Also, I think gffread runs successfully since it only checks that chromosome names match not if the coords are within sequence length bounds. Hopefully somebody else can chime in to help you further. Good luck

I have no experience with this, but I would double-check that the contig names are compatible between all of your reference files (both downloaded and created by you). You could look at the differences between your custom FASTA and GTF and original versions to see if there were any unintentional changes. And if you have indels in your VCF (which would shift coordinates), maybe spot-check a few affected intervals (exons, etc.) in the custom GTF (compared to the original) to see if that lift-over occurred correctly. I am also not convinced that you really need all these custom files for simple gene expression analysis -- the mapper has quite a bit of tolerance for mismatches, so unless you have major structural changes (I don't think deepvariant would show those), any differences should be minimal.

Hi @Thomas Danhorn I have double checked that both the FASTA and GTF file have the same contig names and I still get the same error. The reason I believe I need to use a custom fasta is because after I run the paired WGS and RNAseq through NF-core Sarek and NF-core RNASeq I use phASER for allele specific expression and I lose the insertions that are present in the WGS. So that is why I think I need the custom genome so NF-core RNAseq so I can properly get expression of these insertions and then the phased expression from phASER. If. you think there is an alternative route or easier router than the custom fasta please let me know!

OK, that makes sense. Sorry, I work with somatic variants and don't have experience with phasing or suggestions for that. I have looked into the source code of RSEM and to me it seems that the parsing of the GTF was successful, given the message

Parsing gtf File is done!

is at the end of the function. It is not clear to me what the actual problem was, since there is no specific error message, we just know the program

rsem-extract-reference-transcripts

returned a non-zero exit code. What I would do is make a minimal GTF (maybe one gene) from the one you have and provide the FASTA (maybe just the contig on which that gene is) and run

rsem-extract-reference-transcripts

manually. (You should have everything needed in the work directory of the failed process.) If that succeeds, you can add more genes slowly until you find one that fails (maybe start with the genes that have been modified from the original). (If it fails, you are already at that point.) Then you can take a hard look at that gene in the GTF and the FASTA and see if there are any obvious mistakes. You can change things to see if you can make it work to diagnose the issue. Random idea: Look at your reference files with

less -U

to see if there are any inappropriate control characters in there, e.g.

^M

from Windows line endings. Probably unrelated to your immediate issue, but I noticed

#SBATCH --cpus-per-task=16
#SBATCH --mem=200G

in your script(s). If your config is set up in a sane, efficient way, this makes no sense -- the script you are firing off just starts the

nextflow

process itself, which needs one CPU and maybe 4-8 GB of memory. This process then submits the processes doing actual work (parsing GTFs, mapping reads to the genome, etc.) as SLURM jobs to your cluster (I assume that is the case for you, since you have a SLURM script). I suggest investing a bit of time getting a reasonable config together, it will pay off multiple times over. Here are some tips. 🤞

Thank you for your reply! I am going to start testing

rsem-extract-reference-transcripts

and looking into my config file! Thanks again for the suggestions!