voila modulize error on .psicov

[Maoting Chen]

Hi San,

I tried to run voila modulize on .psicov to get the AS event types on different conditions. Below is how I got the .psicov and .sgc file and ran voila modulize: 

majiq psi-coverage ../build_4rep_group_majiq-files/splicegraph.zarr ./psi_coverage_4rep_sep_group/OP.psicov ../build_4rep_group_majiq-files/OP1.Aligned.out.sj ../build_4rep_group_majiq-files/OP2.Aligned.out.sj ../build_4rep_group_majiq-files/OP7.Aligned.out.sj ../build_4rep_group_majiq-files/OP8.Aligned.out.sj --prefixes OP1 OP2 OP7 OP8 

majiq-v3 sg-coverage ./build_4rep_group_majiq-files/splicegraph.zarr ./sg-coverage_4rep_group/OP.sgc ./build_4rep_group_majiq-files/OP1.Aligned.out.sj ./build_4rep_group_majiq-files/OP2.Aligned.out.sj ./build_4rep_group_majiq-files/OP7.Aligned.out.sj ./build_4rep_group_majiq-files/OP8.Aligned.out.sj --prefixes OP1 OP2 OP7 OP8  (Note: no sg-coverage in majiq but only in majiq-v3)

voila modulize ../build_4rep_group_majiq-files/splicegraph.zarr ./psi_coverage_4rep_sep_group/OP.psicov ../sg-coverage_4rep_group/OP.sgc -d OP_modulize/ --show-all 

The error was reported: trying to get bins and more than one group returned, this should not happen Traceback (most recent call last): File "/n/home09/mtchen/.conda/envs/majiqV3/lib/python3.12/site-packages/rna_voila/classify.py", line 104, in classify_gene writer.alt5prime() File "/n/home09/mtchen/.conda/envs/majiqV3/lib/python3.12/site-packages/rna_voila/classifier/tsv_writer.py", line 783, in alt5prime quants = [event_non_changing, event_changing] + cases_non_changing + cases_changing + self.quantifications(module, 's', edge=event['Proximal']) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "/n/home09/mtchen/.conda/envs/majiqV3/lib/python3.12/site-packages/rna_voila/classifier/quantification_finder.py", line 75, in quantifications quants = self.quantifications_int[field][0](*self.quantifications_int[field][1:])(lsv_id, edge) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "/n/home09/mtchen/.conda/envs/majiqV3/lib/python3.12/site-packages/rna_voila/classifier/quantification_finder.py", line 188, in f return self._inner_edge_aggregate(lsv, generate_variances([lsv.bins][0]), edge) ^^^^^^^^ File "/n/home09/mtchen/.conda/envs/majiqV3/lib/python3.12/site-packages/rna_voila/api/view_matrix_zarr.py", line 453, in bins assert False ^^^^^ AssertionError Processing Genes and Modules [0/1]^M

However, if I ran voila modulize on .voila from majiq deltapsi and the two corresponding .sgc files, there was no problem.

I am using Voila v3.0.17.dev1+g11dc0f550. Could you please kindly help to look into this error? Happy to provide any other information needed.

Thank you,

Maoting

[San Jewell]

Hi Maoting, 

I believe that this was an issue with an untested combination of file/experiment grouping in voila, which was fixed in a recent update. (following the same pip command you used before should show version "voila 3.0.18.dev1+g69f22130d") ; if not, please let me know and I will check deeper into the cause. 

Thanks, 

-San

[Maoting Chen]

Hi San,

The issue was indeed solved after I installed the newest majiq and voila. Thank you for your clarification. 

I checked my gene of interest in the output of modulize result and want to ask how modulize define each AS event. One LSV of interest is defined as"other" in the modulize output, while I personally will think it is alternative first exon. Please see the attached splicegraph and the LSV of interest is gene:WBGene00018514:t:5793996-5794065. Could you please explain the inconsistency of modulize output and human judgement?

Thank youm

maoting

[San Jewell]

Hi Maoting, 

Each of the various modulizer event detection functions runs a specific algorithm to try to figure if something matches any of the very specific patterns we have defined, if it ends up not matching anything, then we call it 'other'. 

In this case the alternate first exon specifically looks for an exon which is not the first exon (numbered 1), but which lacks junctions or introns extending to it's 5' direction. Exon 4 doesn't match this definition because of the de-novo junction, and the half-exon detected at the end of the de-novo junction also doesn't fit because the definition only looks for full exons as alternate first. (However, it would be detected as putative alternate first, which is specified in a separate tsv file) 

Let me know if this makes sense, or if you think this definition is less useful than it could be, I can put your suggestion to review in my lab as a potential improvement. 

Thanks, 

-San

[Maoting Chen]

Hi San, 

Thank you for your clarification.

In this case, how should I evaluate the effect of this de-novo junction on the splicing outcome? Does this de-novo junction mean exon3 has an additional sequence at the 3' end and then is spliced to exon4 directly? Also, this de-novo junction doesn't belong to any LSV. So how is the effect of this de-novo junction demonstrated in the PSI quantification or voila modulize output?

Thank you,

Maoting

[San Jewell]

Hi Maoting, 

I will preface by mentioning I am a little more on the computer science side of this software than the genetics side, so, if you would like, I can ask others in my lab to contribute to this answer if it is not sufficient. That being said, the 'half exon' in this case, was too far from exon 3 for our software to consider it a de-novo exon extension, it doesn't preclude the possibility, but from our informed defaults it has the potential to be not an extension, but rather the 3' end of an additional exon between 3 and 4, but, no evidence was found for the 5' end of that potential exon, so it's left as a half exon. (it would be helpful to look at the distances in coordinates to get an idea of whether this is more or less viable). For your second question, in Majiq we define all quantifications in terms of LSV (local splicing variation) as there is no split in the splicegraph towards the half exon, we can not quantify for or against that path with our model. I would find it interesting if there is a proxy for the inclusion of that junction based on the PSI of the 317-read junction, weighted to the number of reads of the de-novo junction, but currently our model doesn't quantify this. 

Thanks, 

-San