Different results when running Tama on subset vs full BAM file

[skudashev]

Hello,

Thank you for a great tool. I've recently tested Tama on a small region BAM file that was generated from the big BAM using samtools view. Tama managed to correctly assign strand to the transcript models, however after running tama_mapped_sam_splitter.py on the big BAM and then running Tama with the same parameters, it created a different set of isoforms on the wrong strand for the same region. Tama was run with the following parameters: -x no_cap -c 90 -i 80 -rm low_mem -icm ident_cov -a 200 -m 20 -e common_ends. I feel like I am missing something here and would really appreciate some help.

Kind regards,

Sofia

[GenomeRIK]

Hello Sofia,

Thank you for using TAMA!

This is quite strange. Could you show me the sam file lines from the original big bam file and then in the split bam file for this region where the strand seemed to be reversed?

Thank you, Richard

[skudashev]

Hello Richard,

I have tried looking for a reason for the difference and I think it might be that I removed all the secondary and supplementary (all non-primary) alignments from the bam file before running tama on it the second time. For some reason after running samtools view -@ 20 -h -F 2048 -F 2064 -F 256 -F 272 samtools also removed some reads with sam flag 16 (read reverse strand) which led to most of the reads that aligned to the region to be removed from the file. Example reads that all aligned to the same transcript model according to "tama_trans_read.bed" that got removed after the samtools command:

86:477|cce6a542-1e5b-4b1e-a2b3-df63abcc3f70 16 chr17 79089344 60 82M1D15M1310N18M1D29M1D7M1D16M3D6M1D26M2D21M3567N24M2I50M2I5M983N59M5S 0 0 GTCCTGGTTTGAAGAAAGAAATCTTTATTAATAATCTTAATAATACTGTTAGTTTGAATGGTCACACCTTGGAAGCATACAGATGGTAAGAAAGGAAGCCTCTTGGTTTGGTTGGTTTTTTTTTGTTGCTTGGATCTTAACATCTTTTTTGTTTTTTTGTTTTGTGATTTTTGTTTTTTGTTTTTGCACTTCATGGTCCGAAGGAAACGGTGGAAGGTTTCATGGTTCCAATGCTGTAGGTCGCCACGCGGGCCCGATGGTGTGATGGTACGGGTCGGCAGCTGCGTAGACTCTGCTGCGTAACTGTCGCTGTAGGCTGCCGCCGCTGCAGCGGGCTGAGCGTATCTGTAGGCTGCGTAGCCCCCAA NM:i:15 ms:i:325 AS:i:227 nn:i:0 ts:A:- tp:A:P cm:i:80 s1:i:296 s2:i:0 de:f:0.0272 rl:i:42 109:524|282b1a4b-b2ea-4010-a8d5-5ab600a6d728 16 chr17 79089344 60 16M1I77M1311N15M4I7M12I47M1I8M1D33M1D6M1I20M3565N79M983N59M2S 0 0 GTCCTGGTTTGAAGAAGAAAAATCTTTATTAATAATCTTAATAATACTGTTAGTTTGAATGGTCACACCTTGGAAGCATACAGAATGGTAAGAAGGAAGCCTCTTGGTTTTGGTGGTTGGTTTTTTTTTTTTTTTTTTTTTTGTTGCTTGGATCTTAACATCTTTTTTTGTTTTTTTTGTTTTTGTGATTTTTTTGTTTTTTGTTTTTTGCCCTTCATGGTCGAGAAGGGAAACGGTGGAAGGTTTCACATGGTTCCAATGCTGTAGGTCGCCGCGGGCCCGATGGTGTGATGGTACGGGTCGGCAGCTGCGTAGACTCTGCCGTAACTGTCGCCGTAGGCTGCCGCCGCTGCAGCGGGCTGAGCGTATCTGTAGGCTGCGTAGCCCC NM:i:26 ms:i:327 AS:i:221 nn:i:0 ts:A:- tp:A:P cm:i:76 s1:i:247 s2:i:0 de:f:0.0321 rl:i:70 97:494|d5f8497d-0542-4b0a-8998-bee05f242486 16 chr17 79089344 60 1S18M1I12M3I68M1310N5M1D12M1I15M1D6M2I16M2D14M1I62M3565N17M3I2M2I53M5I7M983N10M2I10M2I38M 0 0 GGTCCTGGTTTGAAGAAAGTAAATCTTTATTACTGATAATCTTAATAATACTGTTAGTTTGAATGGTCACACCTTGGAAGCATACAGAATGGTAAGAAAGGAAGCCTCTGGTTTGGTTGGTTTTTTTTTTTGTTGCTGGATCTTTTAACATCTTTTTTTGTTTTTCGTTTTGTGCATTTTTTTTGCTTCTTTGTTTTTGCCCTTCATGGTCCGAGAAGGAAACGGTGGAAGGTTTCACATGGTTCCAATGCTGTGCGTGCCGTCGCCGCGGGCCCGATGGTGTGATGGTACGGGTCGGCAGCTGCGTAGACTCTGAACGGCCGTAACTGTCGCTGTGGAGGCTGCCGCTGTGCTGCAGCGGGCTGAGCGTATCTGTAGGCTGCGTAGC NM:i:31 ms:i:304 AS:i:202 nn:i:0 ts:A:- tp:A:P cm:i:67 s1:i:231 s2:i:0 de:f:0.0476 rl:i:18 72:1715|c7f3a309-9c9f-41c6-b8db-642bc239fc93 16 chr17 79089346 60 19M1I5M1I73M1301N26M2I17M1D20M5I53M2I11M1D2M1D10M3565N9M1D51M2I18M983N38M1I1M1I22M1078N8M5I32M599N18M5D2M2D35M1D48M1D22M267N53M3838N60M1517N15M1I53M4I6M1D20M817N49M1I8M2523N10M1D85M1D29M1I11M8705N62M1D88M2D20M2I14M1D13M2D11M3D30M1I3M1D2M120018N43M71916N11M1D12M2I16M1D6M1D1M1D51M174629N12M1D6M3D16M1D6M4D4M5D4M1D5M2D4M2D7M2I13M2I3M1I4M1D39M128227N14M1I17M1D19M108S 0 0 CCTGGTTTGAAGAAAGAAACTTTTTAAGTAATAATCTTAATAATACTGTTAGTTTGAATGGTCACACCTTGGAAGCATACAGAATGGTAAGAAAGGAAGAGATGGATGCCTCTTGGTTTGGTTGGTTTTTTTTTTTTGTTGCTTGATCTTAACATCTTTTTTTGTTTTCTTTTTTTTGTTTTGTGATTTTTTTTGTTTTTTTGTTTTTGCCCTTCATGGTCCGAGAGAAGGAAACGTGAAGGTTTCACATGGTTCCACGTTGCAGGTCGCTGCGGGCCCGATGGTGTGATGGTACGGGTCGGCAGCTGCGCGTAGACTCTGCCGTAACTGTCGCTGTAGGCTGCCGCCGCTGCAGCGGGCTGAGCTGCAATCTGTAGGCTGCGTAGCCTCCATAAATCTCTTTTCAGCACCATAAAATCCATCCTGATACATGACCGCTCCGTAAGTCGGGATGTAGGTGGCGCAGCCCGAAATGGGCTATACACGGCCCACCCCGGCCCCGAAGATGTGCGCCCTGGTAGGCAACGGCTGTGCCGGTGAGGGGTAGGGGAACCCCGTCACTGCATAGAATTCAGGCCCGTAGACTGCGCCGACCACTGGATTTAGCTTCCAGCCGTTGGTGTAGGGGTTCCCCGTCTTCTTGTTGGTCATCACTCGGGCTGTGGCATTATTGACCTCAATTTTCCGTCCCCTCTACGATCGTCCCATTCAGCTTCTCCCGCCTTCGGTCAGCATCTGAGCTTAGCAGTTTCAAAGCTACAAACCCAAAAACCCTTGGAGCCCTGCTCGTTAAAAATGATCTCCACGTCTAAAATTTTTCCGGAATTGCTCCGAACATTTACGCAAGTCGGGGTCCCTGAACCGGAAGGGGATGTTGGAGACGTGTAGCCGTCAGGGCTGCCGCTTCTCTGTAGGGTCGGAGGTGGGAGCGGCTGGCTGTCCGTCTGTGCCGCCTTCGTCTGTCTGCGGCACTGTCTGGGTCCCGGCGATGGGCTGTGTGCTGGCCCCCGAGCCTGGCTGCTCGGGGTGGTCTGTGCTGGTGTGTACAGGGTCATGCCATGCTCTGTGGGGACCGGGGTCTGGCCGGAGTAGTCCTGCGTGGGGTGCGGTGGGGGCAGCGTACTCGGCAGGGATGCAGCGTTCTGTGGCGGAGGGGGTACTGGGCGGGGGTAGGGCGCCATCGCTTCAGGCGGAGCCGTGGCGTCCTTGACGCCCTTCATTCAGCCCATTCTCCCCGAGGCACTCTGGGCTCTCTCTGTTTGCAGAGGATGGCTCCAGGAGAGCAAGCCCCGTGCTGAGGGGTAAGACCTGGGAAAATGGTGGCAGAATTTCCTCTGGCTGCTGTGACTTCAAGCTGGTAGTGGGAGTGCGGGAATTGCCGAAAGGGGTGCCGCGGGCCTCGCCCTCGGGCGGTGGGGCCGCGCCGGGGGAGGAAAGTAATCGCCGCAGGAAGACTGTTAGCTTCTCAGAGGCAGAATTTCAACGGGCTCAGCGCTCCCCCGCGGCAGGTGACTGGGGCCGGCGGCCATGCCCGTTGAAATTCTTGCCTCTGAGAAGCTAACAGTCTTCCCAGGGTCGCCACCTCAGCAGCCCCCTCCTTGCCAAGGACGGTCCGGTGAAGGAGCTCACCTGGGCCTCCCC NM:i:150 ms:i:1141 AS:i:641 nn:i:0 ts:A:- tp:A:P cm:i:240 s1:i:897 s2:i:0 de:f:0.0724 rl:i:54

I understand that this isn't really a Tama issue but if you have seen an issue like this before, I'd be very grateful if you shared your advice.

Kind regards, Sofia

[skudashev]

Upon reading the samtools documentation, I realised that the -F command excludes all reads with any bits set in the command, and both the 2064 and 272 include the reverse strand bit (0x10). Apologies for taking up your time.