mhalushka
commented
1 year ago Thank you for this insightful work. This is very interesting. I am not at all surprised about the tRNA data having multi-mapping reads. That is expected due to the tRNA families scattered throughout the genome. We had spent a considerable amount of time coming up with the best strategy to report collapsed tRNA results and report them as combined "Comb" in the tRF.Counts and tRF.RP100K. It would be interesting to learn how these final numbers change due to multimapping or if they are relatively stable. Do you have any sense of the multimapping impact on miRNAs? I again suspect that since we collapse similar miRNAs together, this is negligible, but perhaps I am mistaken. In building miRge, we were most interested in counting the species in the RNA-seq data. Thus we worked hard to capture isomiRs. If one is more interested in knowing from which genetic loci the read came from, miRge would not be the best approach and an alignment tool that does not allow mismatches and that aligns directly to the genome would be a better approach (but at the cost of a lot of isomiR data).
Glfrey
arunhpatil
Hi @arunhpatil,
I don't expect this particular pull request to be merged but figure this is best way to discuss what I've implemented so far, and why.
To handle the ambiguous alignment problem I've had to make one major change to miRges operations, which is using dictionary of dictionaries to store the output of the bowtie alignment, before processing into the data frame. This dictionary has the following format: {sequence:{aligned_mol_name_0: full sam alignment information, ..., aligned_mol_name_n: full sam alignment information}}.
The reason for this is because iterating through the alignment line by line and storing the results in the data frame directly means we loose the alignment information for the aligned first molecule unless we want to iterate through the output multiple times (which is an option). I thought carefully about whether a dictionary was the best option given that they can become quite hefty to store in RAM, which adds to the computational burden for users, but for now I think it's the most straightforward for me to implement. We can swap this strategy out for another, more RAM-friendly, alternative if preferred bearing in mind the unavoidable trade-off of time and/or I/O/diskspace costs.
Once the dictionary is generated, I then iterate through the keys (sequences) and for single mappers, simply populate the data frame as was being done previously. For multi mappers I check which alignment has the highest scoring MAPQ and enter that entry into the data frame. For instances where there's a tie for the high scoring alignment, the first alignment that was iterated through is recorded in the data frame with the addition of the "" character I mentioned. I know you said you wanted to hold off that, I included it only so I could double check everything was running as I intended (check that sequences with "" have an entry in the secondary alignment file) and we can remove this as you wish. Multi-mapping alignments are stored in an iteration-specific file e.g. mature_tRNA_secondary_matches.csv and the number of multi-mappers, and the portion of the unique reads this makes up is reported.
On that note, I have downloaded some public data from the ENA, study accession ERP112400 and ran samples ERR2938231, ERR2938232, ERR2938234, ERR2938235 through this new code. I have found that similarly to my data, only tRNAs are reporting multiple equally valid alignments and that proportion of multi-mapping unique reads is very high for mature tRNAs (from 62% (public)-78% (my data)) and lower for primary tRNAs (4.8% (public) - 17.6 (my data)). I'm not highly experienced with tRNAs but I wonder if this impacts the files output when using the
'trfflag (i.e. tRF.Counts.csv, tRF.RP100K.csv, tRFs.aligned.report.tsv)? Incidently, I noticed when using this flag the directorytRFs.samples.tmpisn't being deleted, which I guess it should given the ".tmp" extension?This doesn't help our quest to identify how a multi-mapper might be handled for miRNAs. My guess is that multi mapping miRNAs are most likely to map to sequences from the same family given their highly similar sequence which means the reported counts in miR.Counts.csv and miR.RPM.csv would still be valid. If we wanted to be extra cautious we could implement a subroutine to make sure all the reported alignments for a unique read do, indeed, fall into set of miRNAs with merged names?