Hi all. I am facing an issue with RNAsplice. My f...
# rnasplice
p
Hi all. I am facing an issue with RNAsplice. My feeling is that’s related to the genome I use. Is it best to use gencode? Thanks so much for your help. Here is my code,
Copy code
nextflow run nf-core/rnasplice   --input /shared/projects/rna_methylation_m6a/test_samplesheet.csv   --contrasts /shared/projects/rna_methylation_m6a/contrastsheet.csv   --fasta /shared/bank/homo_sapiens/GRCh38/fasta/Homo_sapiens.GRCh38.dna.primary_assembly.fa   --gtf /shared/bank/homo_sapiens/GRCh38/gtf/Homo_sapiens.GRCh38.101.gtf   --rmats true   --paired true   --sashimi_plots true   -profile ifb_core   --outdir /shared/projects/rna_methylation_m6a/nfcore   --min_samps_gene_expr 2   -resume and I get the following error: aused by:
Process
NFCORE_RNASPLICE:RNASPLICE:ISOFORMSWITCHANALYZER (1)
terminated with an error exit status (1) Command executed: mkdir -p results run_isoformswitchanalyzer.R Homo_sapiens.GRCh38.101.gtf genome.transcripts.fa test_samplesheet.csv contrastsheet.csv 0.05 0.1 cat <-END_VERSIONS versions.yml “NFCORE_RNASPLICERNASPLICEISOFORMSWITCHANALYZER”: r-base: $(echo $(R --version 2>&1) | sed ‘s/^.*R version //; s/ .*$//’) bioconductor-stager: $(Rscript -e “library(IsoformSwitchAnalyzeR); cat(as.character(packageVersion(‘IsoformSwitchAnalyzeR’)))“) END_VERSIONS Command exit status: 1 Command output: (empty) Command error: INFO: Environment variable SINGULARITYENV_TMPDIR is set, but APPTAINERENV_TMPDIR is preferred INFO: Environment variable SINGULARITYENV_NXF_TASK_WORKDIR is set, but APPTAINERENV_NXF_TASK_WORKDIR is preferred INFO: Environment variable SINGULARITYENV_NXF_DEBUG is set, but APPTAINERENV_NXF_DEBUG is preferred WARNING: While bind mounting ‘/shared/shared’ destination is already in the mount point list Warning message: In read.table(file = file, header = header, sep = sep, quote = quote, : incomplete final line found by readTableHeader on ‘contrastsheet.csv’ Error in importRdata(isoformCountMatrix = salmonQuant$counts, isoformRepExpression = salmonQuant$abundance, : The annotation and quantification (count/abundance matrix and isoform annotation) seems to be different (Jaccard similarity < 0.925). Either isforoms found in the annotation are not quantifed or vise versa. Specifically: 228801 isoforms were quantified. 229487 isoforms are annotated. Only 0 overlap. 228801 isoforms quantifed had no corresponding annoation This combination cannot be analyzed since it will cause discrepencies between quantification and annotation thereby skewing all analysis. If there is no overlap (as in zero or close) there are two options: 1) The files do not fit together (e.g. different databases, versions, etc) (no fix except using propperly paired files). 2) It is somthing to do with how the isoform ids are stored in the different files. This problem might be solvable using some of the ‘ignoreAfterBar’, ‘ignoreAfterSpace’ or ‘ignoreAfterPeriod’ arguments. Examples from expression matrix are : ENST00000613533, ENST00000470736, ENST00000557538 Examples of annoation are : ENST00000559930.5, ENST00000648220.1, ENST00000598390.5 Examples of isoforms which were only found im the quantification are : ENST00000480948, ENST00000371658, ENST00000452305 If there is a large overlap but still far from complete there are 3 possibilites: 1) The files do not fit together (e.g different databases versions etc.) (no fix except using propperly paired files). 2) If you are using Ensembl data you have supplied the GTF without phaplotyps. You need to supply the Ensembl_version.chr_patch_hapl_scaff.gtf file - NOT the Ensembl_version.chr.gtf 3) One file could contain non-chanonical chromosomes while the other do not (might be solved using the ‘removeNonConvensionalChr’ argument.) 4) It is somthing to do with how a subset of the isoform ids are stored in the different files. This problem might be solvable using some of the ‘ignoreAfterBar’, ‘ignoreAfterSpace’ or ‘ignoreAfterPeriod’ arguments. For more info see the FAQ Execution halted
t
Hmm, I have run rnaspliuce successfully with reference files from Ensembl, including a "regular" GTF, but that was a while ago, and the pipeline did not use the
ISOFORMSWITCHANALYZER
process. I am not clear where the files it uses come from, but it looks like they are not compatible:
Copy code
Examples from expression matrix are : ENST00000613533, ENST00000470736, ENST00000557538
    Examples of annoation are : ENST00000559930.5, ENST00000648220.1, ENST00000598390.5
    Examples of isoforms which were only found im the quantification are : ENST00000480948, ENST00000371658, ENST00000452305
If these are examples from your files, the problem is that one uses the versions after the transcript ID (e.g.
.1
), while the other does not. One of those files is almost certainly the GTF, but I am not clear where the other comes from. There is another tip, though:
Copy code
2) If you are using Ensembl data you have supplied the GTF without phaplotyps. You need to supply the <Ensembl_version>.chr_patch_hapl_scaff.gtf file - NOT the <Ensembl_version>.chr.gtf
Since I have not used this, I can't say much about it, but since the developers seem to have strong feelings about it, I would take it seriously. 🤞
p
Thank you!!
👍 1
s
Hi! I have the same issue running dev version of the pipeline with Ensembl annotation.
Copy code
nextflow run nf-core/rnasplice -r dev \
 -profile singularity \
 --input "$PWD"/samplesheet.csv \
 --contrasts "$PWD"/contrasts.csv \
 --leafcutter true \
 --gtf "$GTF" \
 --sashimi_plot false \
 --fasta "$FASTA" \
 --save_reference \
 --outdir $PWD \
 -resume

Error in importRdata(isoformCountMatrix = salmonQuant$counts, isoformRepExpression = salmonQuant$abundance,  :                                 
    The annotation and quantification (count/abundance matrix and isoform annotation) seems to be different (Jaccard similarity < 0.925).        
  Either isforoms found in the annotation are not quantifed or vise versa.                                                                       
  Specifically:                                                                                                                                     253220 isoforms were quantified.                                                                                                              
   254129 isoforms are annotated.                                                                                                                
   Only 0 overlap.                                                                                                                               
   253220 isoforms quantifed had no corresponding annoation                                                                                      
                                                                                                                                                 
  This combination cannot be analyzed since it will cause discrepencies between quantification and annotation thereby skewing all analysis.      
                                                                                                                                                 
  If there is no overlap (as in zero or close) there are two options:                                                                            
   1) The files do not fit together (e.g. different databases, versions, etc) (no fix except using propperly paired files).                      
   2) It is somthing to do with how the isoform ids are stored in the different files. This problem might be solvable using some of the 'ignoreAf
terBar', 'ignoreAfterSpace' or 'ignoreAfterPeriod' arguments.                                                                                    
       Examples from expression matrix are : ENST00000583894, ENST00000544125, ENST00000496497                                                   
       Examples of annoation are : ENST00000358894.11, ENST00000548136.1, ENST00000577166.5                                                      
       Examples of isoforms which were only found im the quantification are  : ENST00000528501, ENST00000365067, ENST00000392028
I went to workdir and all annotation files as well as salmon quant files seem to have non-dotted transcript id.
Copy code
grep ENST00000358894 Homo_sapiens.GRCh38.112.gtf
20      ensembl_havana  three_prime_utr 63330720        63334806        .       +       .       gene_id "ENSG00000101203"; gene_version "17"; transcript_id "ENST00000358894"; transcript_version "11"; gene_name "COL20A1"; gene_source "ensembl_havana"; gene_biotype "protein_coding"; transcript_name "COL20A1-201"; transcript_source "ensembl_havana"; transcript_biotype "protein_coding"; tag "CCDS"; ccds_id "CCDS46628"; tag "basic"; tag "Ensembl_canonical"; tag "GENCODE Primary"; tag "MANE_Select"; transcript_support_level "1 (assigned to previous version 10)";

grep ENST00000358894 genome.transcripts.fa 
>ENST00000358894

grep ENST00000358894 IR_2p5Gy_48h_1/quant.sf 
ENST00000358894 8045    7807.825        0.027318        16.647
There seems to be no input file where those dotted versions may have come from 🤷 Would be cool to know how you solved the problem 🙂 Thanks!!
t
I have no further insights beside what I wrote above. Unless @Pierre Klein has a solution to share, I would suggest doing as much detective work as you can to find the source of the gene ID and starting a new thread with what you found out ...
s
It appears that adding 'ignoreAfterPeriod' option does circumvent this specific problem. I am not ready to start a new thread now since I have more problems with the tool, but if it might help someone I did this:
Copy code
SwitchList <- importRdata(
    isoformCountMatrix   = salmonQuant$counts,
    isoformRepExpression = salmonQuant$abundance,
    designMatrix         = design,
    isoformExonAnnoation = gtf,
    isoformNtFasta       = transcript_sequences,
    comparisonsToMake    = contrasts,
    showProgress = FALSE,
    quiet = TRUE,
    ignoreAfterPeriod=TRUE #for Ensembl annotations only
)
t
Well,
ignoreAfterPeriod
will just ignore the version number of the gene name, but just because it doesn't crash doesn't mean the analysis is correct. But I would want to figure out where those names come from. You don't want to mix different gene annotations (e.g. Ensembl and Gencode), because even though most of the IDs may the same (if you ignore versions), some will not be (because they were newly introduced in one of the versions, or deprecated, etc.).
s
The version numbers are added explicitly from the very same gtf - see line 1364 in the source code (https://github.com/kvittingseerup/IsoformSwitchAnalyzeR/blob/master/R/import_data.R) . So I am absolutely sure it is the same annotation.
👍 1
t
If this is the case, that appears to be a bug in the pipeline. Maybe you can put together a reproducible test case and open an issue on GitHub ...