ORF error?

[moschmi]

Hello,
I'm getting consistently this error (the KeyError sometimes changes):
Traceback (most recent call last): File "/home/mos6/venv/bin/jcast", line 8, in <module> sys.exit(main()) File "/home/mos6/venv/lib/python3.8/site-packages/jcast/main.py", line 438, in main args.func(args) File "/home/mos6/venv/lib/python3.8/site-packages/jcast/main.py", line 166, in runjcast main_log.info(translate_one_partial(jx)) File "/home/mos6/venv/lib/python3.8/site-packages/jcast/main.py", line 208, in _translate_one sequence.translate(use_phase=True) File "/home/mos6/venv/lib/python3.8/site-packages/jcast/sequences.py", line 145, in translate self.slice1_aa = h.make_pep(self.slice1_nt.seq, self.j.strand, i, terminate=True) File "/home/mos6/venv/lib/python3.8/site-packages/jcast/helpers.py", line 227, in make_pep aa = code[nt[i:i + 3]] KeyError: 'AGX'

I'm using this command with reference files from gencode.
jcast /home/mos6/RNAseq/yumm_RNAseq/rMats_clones/rmats_output/rmats_UV3R-D9_novel /home/references/gencode_ref_mm/gencode.vM25.annotation.gtf /home/references/gencode_ref_mm/GRCm38.primary_assembly.genome.fa -o jcast_D9_novel #error

As recommended I looked for NA in the GTF and Ns in the fasta. None in the GTF and Ns are only in the beginning of each chromosome.

I tried this also with references from NCBI, Ensembl, UCSC.
It looks like there is something wrong with the ORF? What do I need to fix?
Best
Marcel

[ed-lau]

It looks like one of the retrieved codons reads AGX which is not in the codon table. I assume the genome/annotation files are what you used for RNA-seq alignment and rMATS? HAlso have you tried the GRCm38.p6 genome fa file on Gencode?

[moschmi]

I used gencode release 25 which is a GRCm38.p6.
I used the same references for hisat2, and rMATS.
The error is not consistent. If I repeat it with exactly the same parameters the error stops the process at a different gene.

[ed-lau]

I would like to look into this. Would you mind sharing your rMATS output?

[moschmi]

Sure, I uploaded some examples here:
https://github.com/moschmi/files
The A5_G4_ worked once. The others never always giving the ORF error.

[ed-lau]

I was able to produce at least one error consistently from the A5_G4 sample. This particular rMATS result line in the A3SS file gave an error:

13271	"ENSMUSG00000053483"	"Usp21"	chr1	-	171282496	171282611	171282496	171282594	171282694	171283167	13271	119	1	220	0	166	149	1	1.0	0.991	1.0-0.009

Apparently the retrieved nucleotide sequence has a hard mask.

INFO:jcast.seq:Retrieved nucleotide for A3SS 13271 Usp21 slice 1: CGGGAGTCATTGTAGACATGCCAACCAGTCTGGCACCGGCACAGGGCTGTGTAGTGGCCATAGTGGACACTACCCGAGTGGTTGCAAAGGGCATACAGCTGGTAGACAGGGCTTCCCGCTTTGTCGCTGGCAAAGTCCCCTAGGCTTAGTCTCTGCAGTGGGAAGTCTACACCTACTGAACTTTTCTTGATGGAGCCTCGGGAGGTGGAAAATCGGTTCAGATCTAAGCCCTGGTTAAGGAGCACTGGGGAGGCTAAGGGTCTCTAGGGGGGACTGCAGTAATAAATTCAGTACCCTCCATTCCCCCAACACATCCTGAAAGGTGGCAAAGAAAGTCAGGGTGGCAGGTTTTCCCTAAGGGAAACTCAAACGAAGGCGAGGAGACACTGTATGGGAGATACCCTAAGAGAAATAAGAACGTGAAGATTAGGATATGGAGCACGAGGATTCGGGGGAATCTTTGTACTGTCAACTTTTTGGTACTTCGGGTTTTCTGCCGACATCGGTCACACACCTAAATAAACAGCAAATAGTGTTAGGTCAGGGTGCCCTCACCTAAGTTGCTTTTTC**N**TGACACCCCCAACTTAC

This is a known issue that JCAST does not deal with masked nucleotides. I assume hisat2 could align to masked regions? Which I wasn't aware of.

For the time being it could be a quick temporary fix to simply not translate any slices that contain masked nucleotides ('NnXx'). I will push a new version to GitHub today or tomorrow after doing a bit more testing that should circumvent this issue.

[moschmi]

Thanks a lot.
I wasn't aware of that as well. I'll go over the sequences and exclude those. I I also thought to just replace all the Ns with As. That should do the trick.
Can't wait for the new version...

[ed-lau]

Can you try to see if v0.3.5a patched the issue?

[moschmi]

I ran jcast with A5 but got the same error, trying now with A5_G4
error below:
INFO:root:Extracted GTF attributes: ['gene_id', 'gene_type', 'gene_name', 'level', 'mgi_id', 'havana_gene', 'transcript_id', 'transcript_type', 'transcript_name', 'transcript_support_level', 'tag', 'havana_transcript', 'exon_number', 'exon_id', 'protein_id', 'ccdsid', 'ont'] Processing MXE Junctions: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████| 1844/1844 [21:24<00:00, 1.44it/s] Processing SE Junctions: 100%|██████████████████████████████████████████████████████████████████████████████████████████████| 30831/30831 [4:25:23<00:00, 1.94it/s] Processing RI Junctions: 100%|██████████████████████████████████████████████████████████████████████████████████████████████████| 8423/8423 [54:25<00:00, 2.58it/s] Processing A5SS Junctions: 100%|████████████████████████████████████████████████████████████████████████████████████████████| 10318/10318 [1:03:06<00:00, 2.73it/s] Processing A3SS Junctions: 75%|█████████████████████████████████████████████████████████████████████ | 11214/14948 [1:05:50<21:55, 2.84it/s] Traceback (most recent call last): File "/etc/samba/mos6/user_homes/mos6/anaconda3/bin/jcast", line 8, in <module> sys.exit(main()) File "/etc/samba/mos6/user_homes/mos6/anaconda3/lib/python3.7/site-packages/jcast/main.py", line 439, in main args.func(args) File "/etc/samba/mos6/user_homes/mos6/anaconda3/lib/python3.7/site-packages/jcast/main.py", line 166, in runjcast main_log.info(translate_one_partial(jx)) File "/etc/samba/mos6/user_homes/mos6/anaconda3/lib/python3.7/site-packages/jcast/main.py", line 209, in _translate_one sequence.translate(use_phase=True) File "/etc/samba/mos6/user_homes/mos6/anaconda3/lib/python3.7/site-packages/jcast/sequences.py", line 145, in translate self.slice1_aa = h.make_pep(self.slice1_nt.seq, self.j.strand, i, terminate=True) File "/etc/samba/mos6/user_homes/mos6/anaconda3/lib/python3.7/site-packages/jcast/helpers.py", line 227, in make_pep aa = code[nt[i:i + 3]] KeyError: 'CAX'

[ed-lau]

Just making sure that v.0.3.5a is being used by your conda environment? Can you upload the log file?

[moschmi]

I think I ran v.0.3.5:
2022-11-01 00:23:06,510 - jcast - INFO - Namespace(canonical=False, func=<function runjcast at 0x7f00eb519710>, g_or_ln=0, genome='/etc/samba/mos6/user_homes/mos6/genes_mm/GRCm38.primary_assembly.genome.fa', gtf_file='/etc/samba/mos6/user_homes/mos6/genes_mm/gencode.vM25.annotation.gtf', model=False, out='jcast_A5_novel', qvalue=[0, 1], read=1, rmats_folder='/etc/samba/mos6/user_homes/mos6/mel_LN_Attac/rmats_output/rmats_UV2R-A5_novel') 2022-11-01 00:23:06,510 - jcast - INFO - 0.3.5 2022-11-01 00:23:06,826 - jcast.gtf - INFO - Reading GTF file. This could take a minute. 2022-11-01 00:25:39,208 - jcast.genome - INFO - Read from genome.

Sorry I need to do a clean install

[ed-lau]

Got it - v0.3.5a is not on PyPI yet, but you can pip install the zip file inside the v0.3.5a link, or install from the repo directly.

[moschmi]

Yay, looks like it's working! I'm running my samples and so far no error.
Thanks a lot for the super quick fix

[moschmi]

Unfortunately, while there is no error, the pipeline just stops after 5-6 h and the only output is the log file.
It stops during the A3SS step at various genes.
Here's the tail:
tail jcast_A5_novel5a/jcast_20221103040257/jcast_main.log 2022-11-03 08:54:46,154 - jcast.junction - INFO - Anchor exon start 48929681 Anchor exon end 48929729 2022-11-03 08:54:46,161 - jcast.junction - INFO - Transcription start: None Transcript end: -1 2022-11-03 08:54:46,161 - jcast.junction - INFO - Retrieved phase: None 2022-11-03 08:54:46,162 - jcast.seq - INFO - Retrieved nucleotide for A3SS 14947 Orc4 slice 1: CAATCTACATGTAAGTCCAATAACTGCTACTGGTGTCTGTGCAGATTGAGAAATGTCAAAAAGATTATAAAGGAGTGTTTGATTTTTCTGATGAGCAAAAATATCAAATTCATCCAGTATGAAGATCACTGGGCAACTGCTAGTCCGGTCACCTTTTTGCAAAGCTTCCAGAAGAAATGAAAGGTTTTCAGCAAAGCTTCC 2022-11-03 08:54:46,162 - jcast.seq - INFO - Retrieved nucleotide for A3SS 14947 Orc4 slice 2: CAATCTACATGTAAGTCCAATAACTGCTACTGGTGTCTTTTTGCAAAGCTTCCAGAAGAAATGAAAGGTTTTCAGCAAAGCTTCC 2022-11-03 08:54:46,172 - jcast.seq - INFO - Locally cached sequence retrieved for ENSMUSG00000026761. 2022-11-03 08:54:46,172 - jcast.seq - INFO - Using non-phase translation for 14947 Orc4 2022-11-03 08:54:46,173 - jcast - INFO - SKIPPED. Sequence discarded due to low coverage.

[ed-lau]

I was able to complete the A5_novel sample from your rMATS output after ~6-7 hours with the -q 0 1 and -c options. Happy to send you the output files if that helps. The SE took a while with 30000+ junctions because JCAST currently has to make a call to UniProt to retrieve the canonical protein sequence for each gene. If you suspect this is due to a time out on your computer I might suggest using the read count filter -r to remove low-read junctions, or using -s MXE RI etc. to only run some splice types at a time.