Hello, I was wondering if the nf-core/rnaseq pipel...
# rnaseq
m
Hello, I was wondering if the nf-core/rnaseq pipeline has a problem with .fna files for the reference genome.
m
Can you elaborate? The schema allows for that ending and I am unaware that there are any differences in the file contents:
Copy code
"fasta": {
                    "type": "string",
                    "format": "file-path",
                    "exists": true,
                    "mimetype": "text/plain",
                    "pattern": "^\\S+\\.fn?a(sta)?(\\.gz)?$",
                    "description": "Path to FASTA genome file.",
                    "help_text": "This parameter is mandatory if `--genome` is not specified. If you don't have the appropriate alignment index, it will be generated automatically. Use with `--save_reference` to store the index for future runs.",
                    "fa_icon": "far fa-file-code",
                    "errorMessage": "The FASTA file path must end with .fa, .fna, .fasta optionally with .gz, must not contain spaces, and must exist."
                },
What kind of error are you getting?
m
ERROR ~ Error executing process > 'NFCORE_RNASEQPREPARE GENOMEMAKE_TRANSCRIPTS_FASTA (rsem/GCF_000001215.4_Release_6_plus_ISO1_MT_genomic.fna)' Caused by: Process
NFCORE_RNASEQ:PREPARE_GENOME:MAKE_TRANSCRIPTS_FASTA (rsem/GCF_000001215.4_Release_6_plus_ISO1_MT_genomic.fna)
terminated with an error exit status (255) Command executed: rsem-prepare-reference \ --gtf GCF_000001215.4_Release_6_plus_ISO1_MT_genomic.filtered.filtered.gtf \ --num-threads 12 \ \ rsem/GCF_000001215.4_Release_6_plus_ISO1_MT_genomic.fna \ rsem/genome cp rsem/genome.transcripts.fa . cat <<-END_VERSIONS > versions.yml "NFCORE_RNASEQPREPARE GENOMEMAKE_TRANSCRIPTS_FASTA": rsem: $(rsem-calculate-expression --version | sed -e "s/Current version: RSEM v//g") star: $(STAR --version | sed -e "s/STAR_//g") END_VERSIONS Command exit status: 255 Command output: rsem-extract-reference-transcripts rsem/genome 0 GCF_000001215.4_Release_6_plus_ISO1_MT_genomic.filtered.filtered.gtf None 0 rsem/GCF_000001215.4_Release_6_plus_ISO1_MT_genomic.fna Parsed 200000 lines Parsed 400000 lines "rsem-extract-reference-transcripts rsem/genome 0 GCF_000001215.4_Release_6_plus_ISO1_MT_genomic.filtered.filtered.gtf None 0 rsem/GCF_000001215.4_Release_6_plus_ISO1_MT_genomic.fna" failed! Plase check if you provide correct parameters/options for the pipeline! Command error: rsem-extract-reference-transcripts rsem/genome 0 GCF_000001215.4_Release_6_plus_ISO1_MT_genomic.filtered.filtered.gtf None 0 rsem/GCF_000001215.4_Release_6_plus_ISO1_MT_genomic.fna According to the GTF file given, transcript NM_001260299.1 has exons from different orientations! Parsed 200000 lines Parsed 400000 lines "rsem-extract-reference-transcripts rsem/genome 0 GCF_000001215.4_Release_6_plus_ISO1_MT_genomic.filtered.filtered.gtf None 0 rsem/GCF_000001215.4_Release_6_plus_ISO1_MT_genomic.fna" failed! Plase check if you provide correct parameters/options for the pipeline!
So I get the error above, but I am not sure why?
My pipeline is very simple as I am in just try now. nextflow run nf-core/rnaseq \ --input samplesheet.csv \ -profile conda \ --save_align_intermeds \ --outdir results_initial \ --skipBiotypeQC \ --fasta /home/nspmpego/Dm6_plus_ISO1_MT/GCF/GCF_000001215.4_Release_6_plus_ISO1_MT_genomic.fna\ --gtf /home/nspmpego/Dm6_plus_ISO1_MT/GCF/GCF_000001215.4_Release_6_plus_ISO1_MT_genomic.filtered.gtf \ --save_reference
m
The relevant part of the error message is:
Copy code
Command error:
  rsem-extract-reference-transcripts rsem/genome 0 GCF_000001215.4_Release_6_plus_ISO1_MT_genomic.filtered.filtered.gtf None 0 rsem/GCF_000001215.4_Release_6_plus_ISO1_MT_genomic.fna
  According to the GTF file given, transcript NM_001260299.1 has exons from different orientations!
I am unfamiliar with Drosophila as a model organism, so I do not know whether that annotation has a bug or whether this is indeed some fancy splicing event that breaks common tool assumptions? I recall that the 1bp exon in Arabidopsis also broke the pipeline…
m
I can try to download the files from somewhere else and see how it goes. Is there anything wrong with the pipeline?
m
No the pipeline is fine. But that one gene is apparently breaking the tool
rsem-extract-reference-transcripts
, because trans-splicing of from different DNA strands is not something commonly encountered in other organisms. So unfortunately, you will have to remove that transcript from your GTF and need to manually quantify it if it is of relevance for your research project.
m
Thank you. I can try to remove the mdg4. It was the one giving me trouble with "?" for the strand column in the GTF. I'll try tomorrow because my brain hurts now...
t
Maybe have a look how Ensembl deals with this gene/transcript ...
🙌 1
m
Hello, another day, another problem. I managed to get the pipeline to move forward quite a bit, but I am now presented with another failure. my pipeline: nextflow run nf-core/rnaseq \ --input samplesheet.csv \ -profile conda \ --save_align_intermeds \ --outdir results_initial \ --skipBiotypeQC \ --fasta /home/nspmpego/Dm6_plus_ISO1_MT/GCF/GCF_000001215.4_Release_6_plus_ISO1_MT_genomic.fna\ --gtf /home/nspmpego/Dm6_plus_ISO1_MT/GCF/GCF_000001215.4_Release_6_plus_ISO1_MT_genomic.filtered.2.gtf \ --save_reference
The Error message: ERROR ~ Error executing process > 'NFCORE_RNASEQRNASEQBAM_RSEQC:RSEQC_READDISTRIBUTION (Axe_SD_rep2)' Caused by: Process
NFCORE_RNASEQ:RNASEQ:BAM_RSEQC:RSEQC_READDISTRIBUTION (Axe_SD_rep2)
terminated with an error exit status (1) Command executed: read_distribution.py \ -i Axe_SD_rep2.markdup.sorted.bam \ -r GCF_000001215.4_Release_6_plus_ISO1_MT_genomic.filtered.2.filtered.bed \ > Axe_SD_rep2.read_distribution.txt cat <<-END_VERSIONS > versions.yml "NFCORE_RNASEQRNASEQBAM_RSEQC:RSEQC_READDISTRIBUTION": rseqc: $(read_distribution.py --version | sed -e "s/read_distribution.py //g") END_VERSIONS Command exit status: 1 Command output: (empty) Command error: processing GCF_000001215.4_Release_6_plus_ISO1_MT_genomic.filtered.2.filtered.bed ... Traceback (most recent call last): File "/home/nspmpego/microbiota_photoperiodism/gut_RNA_seq/ld_sd_rnaseq_analysis/work/conda/env-0e90878ad6b4a0932ac4d3282b759d1e/bin/read_distribution.py", line 302, in <module> main() File "/home/nspmpego/microbiota_photoperiodism/gut_RNA_seq/ld_sd_rnaseq_analysis/work/conda/env-0e90878ad6b4a0932ac4d3282b759d1e/bin/read_distribution.py", line 191, in main intergenic_down1kb_base,intergenic_down5kb_base,intergenic_down10kb_base) = process_gene_model(options.ref_gene_model) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "/home/nspmpego/microbiota_photoperiodism/gut_RNA_seq/ld_sd_rnaseq_analysis/work/conda/env-0e90878ad6b4a0932ac4d3282b759d1e/bin/read_distribution.py", line 81, in process_gene_model cds_exon=BED.unionBed3(cds_exon) ^^^^^^^^^^^^^^^^^^^^^^^ File "/home/nspmpego/microbiota_photoperiodism/gut_RNA_seq/ld_sd_rnaseq_analysis/work/conda/env-0e90878ad6b4a0932ac4d3282b759d1e/lib/python3.12/site-packages/qcmodule/BED.py", line 2504, in unionBed3 bitsets = binned_bitsets_from_list(lst) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "/home/nspmpego/microbiota_photoperiodism/gut_RNA_seq/ld_sd_rnaseq_analysis/work/conda/env-0e90878ad6b4a0932ac4d3282b759d1e/lib/python3.12/site-packages/bx/bitset_builders.py", line 155, in binned_bitsets_from_list last_bitset.set_range(start, end - start) File "lib/bx/bitset.pyx", line 217, in bx.bitset.BinnedBitSet.set_range File "lib/bx/bitset.pyx", line 187, in bx.bitset.bb_check_range_count IndexError: Count (-456) must be non-negative. Work dir: /home/nspmpego/microbiota_photoperiodism/gut_RNA_seq/ld_sd_rnaseq_analysis/work/8f/17d3d39315079c440343745150bbef
Do you have any suggestions on how to fix the problem? or what the problem may be?
t
This looks like some issue with a binary file to me, which may or may not be from some random corruption ... Is this error reproducible?
m
I can try and run it again. I'll post here when I have done so.
🤞 1
Hello, I have attached the execution report of my last run. The pipeline failed again but perhaps with a different error compared to yesterday.... I am not totally sure. Any suggestions? Do you think there may be an issue with my installation?
Hello, Just wanted to let you know that I managed to run the pipeline successfully. I just downloaded a new GTF and .fa files from ensembl (ensembl_BDGP6.54_GCA_000001215.4). I cat together the primary assembly .fa files. and I run the following pipeline nextflow run nf-core/rnaseq \ --input samplesheet.csv \ -profile conda \ --save_align_intermeds \ --outdir results_initial \ --skipBiotypeQC \ --fasta /home/nspmpego/ensembl_BDGP6.54_GCA_000001215.4/drosophila_genome_assembly.fa\ --gtf /home/nspmpego/ensembl_BDGP6.54_GCA_000001215.4/Drosophila_melanogaster.BDGP6.54.115.gtf\ --save_reference
🎉 2
t
Glad you got it to work. I suspect your original GTF had an issue that caused the error.
👍 1