Tweet
Hello,
I am working with RNAseq data, produced on the Illumina Hiseq platform. My reads are paired-end and 76bp long. I map them against the mouse genome with Tophat.
I have a question regarding how Tophat scores and reports the alignments. As far as I understand, the alignment with the higher score will be reported; in cases where several alignments have the same score, it will report as many as allowed by the --max-multihits option. However, the alignment score is not reported in the output BAM file, since all unique alignments get 255 (in the mapping quality field).
I have reads aligned to the mt-Co2 (ENSMUSG00000064354) gene in the mitochondrial chromosome, and they are reported as unique hits. However, this gene has an identical copy in chromosome 1, so all the reads also map perfectly there (as confirmed by BLAST), but those alignments are not reported. I have checked that the alignments in chr1 are in the correct orientation and have similar (or even identical) insert sizes.
Any ideas on why Tophat is not reporting them? Am I missing something here?
Here are some examples from the SAM output file:
IL13_5283:1:100:10005:4660 147 MT 7602 255 76M = 7363 -315 CTGAAATTTGTGGATCTAACCATAGCTTTATGCCCATTGTCCTAGAAATGGTTCCACTAAAATATTTCGAAAACTG EEE
ECDBEEEEEEEAEEEDEEE?EEECEECEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEECEE NM:i:0 NH:i:1 XF:Z:ENSMUSG00000064354
IL13_5283:1:100:10005:4660 99 MT 7363 255 76M = 7602 315 CTTTGATTCATATATAATCCCAACAAACGACCTAAAACCTGGTGAACTACGACTGCTAGAAGTTGATAACCGAGTC HHH
HHFHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH@HHHHHHHHFHHHHHEH@BH>HHFCEHCBCCF=D NM:i:0 NH:i:1 XF:Z:ENSMUSG00000064354
IL13_5283:1:100:10007:10895 147 MT 7282 255 76M = 7040 -318 CAACAACCCCGTATTAACCGTTAAAACCATAGGGCACCAATGATACTGAAGCTACGAATATACTGACTATGAAGAC ??=
??7B??=<<A>AD@D?B?ADAD=DDDDDDBDD?DDDDDDDD?DADDDDDDDDDDDDDDDDDDDDBDDDDDDDD NM:i:0 NH:i:1 XF:Z:ENSMUSG00000064354
IL13_5283:1:100:10007:10895 99 MT 7040 255 76M = 7282 318 CAAGACGCCACATCCCCTATTATAGAAGAGCTAATAAATTTCCATGATCACACACTAATAATTGTTTTCCTAATTA HHH
HHHHHHHHHHHHHHHHHHHHHGFHHDHHHGHHHHHHHHHGGHHAGHEFFEDEHAGHCBFEH?DEDBAFBAFF@ NM:i:0 NH:i:1 XF:Z:ENSMUSG00000064354
IL13_5283:1:100:10015:17002 163 MT 7026 255 76M = 7318 368 TCCAACTTGGTCTACAAGACGCCACATCCCCTATTATAGAAGAGCTAATAACTTTCCATGATAACACACTAATAAT 969
56(4959@<505:3086=87::31286>>962@:,5.@:2569;2>4;.579>+2<3:0((3(+:<)5;=5*A NM:i:2 NH:i:1 XF:Z:ENSMUSG00000064354
IL13_5283:1:100:10015:17002 83 MT 7318 255 76M = 7026 -368 CCAATGATACTGAAGCTACGAATATACTGACTATGAAGACCTATGCTTTGATTCATATATAATCCCAACAAACGAC 8:;
=;=:87::8<:::=7;8:=:6:?<>=:<5?5::<=6.<8<=>D<789853886@+>CD@*9>57498:6B9;B NM:i:0 NH:i:1 XF:Z:ENSMUSG00000064354
All suggestions are very much appreciated!
I am working with RNAseq data, produced on the Illumina Hiseq platform. My reads are paired-end and 76bp long. I map them against the mouse genome with Tophat.
I have a question regarding how Tophat scores and reports the alignments. As far as I understand, the alignment with the higher score will be reported; in cases where several alignments have the same score, it will report as many as allowed by the --max-multihits option. However, the alignment score is not reported in the output BAM file, since all unique alignments get 255 (in the mapping quality field).
I have reads aligned to the mt-Co2 (ENSMUSG00000064354) gene in the mitochondrial chromosome, and they are reported as unique hits. However, this gene has an identical copy in chromosome 1, so all the reads also map perfectly there (as confirmed by BLAST), but those alignments are not reported. I have checked that the alignments in chr1 are in the correct orientation and have similar (or even identical) insert sizes.
Any ideas on why Tophat is not reporting them? Am I missing something here?
Here are some examples from the SAM output file:
IL13_5283:1:100:10005:4660 147 MT 7602 255 76M = 7363 -315 CTGAAATTTGTGGATCTAACCATAGCTTTATGCCCATTGTCCTAGAAATGGTTCCACTAAAATATTTCGAAAACTG EEE
ECDBEEEEEEEAEEEDEEE?EEECEECEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEECEE NM:i:0 NH:i:1 XF:Z:ENSMUSG00000064354
IL13_5283:1:100:10005:4660 99 MT 7363 255 76M = 7602 315 CTTTGATTCATATATAATCCCAACAAACGACCTAAAACCTGGTGAACTACGACTGCTAGAAGTTGATAACCGAGTC HHH
HHFHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH@HHHHHHHHFHHHHHEH@BH>HHFCEHCBCCF=D NM:i:0 NH:i:1 XF:Z:ENSMUSG00000064354
IL13_5283:1:100:10007:10895 147 MT 7282 255 76M = 7040 -318 CAACAACCCCGTATTAACCGTTAAAACCATAGGGCACCAATGATACTGAAGCTACGAATATACTGACTATGAAGAC ??=
??7B??=<<A>AD@D?B?ADAD=DDDDDDBDD?DDDDDDDD?DADDDDDDDDDDDDDDDDDDDDBDDDDDDDD NM:i:0 NH:i:1 XF:Z:ENSMUSG00000064354
IL13_5283:1:100:10007:10895 99 MT 7040 255 76M = 7282 318 CAAGACGCCACATCCCCTATTATAGAAGAGCTAATAAATTTCCATGATCACACACTAATAATTGTTTTCCTAATTA HHH
HHHHHHHHHHHHHHHHHHHHHGFHHDHHHGHHHHHHHHHGGHHAGHEFFEDEHAGHCBFEH?DEDBAFBAFF@ NM:i:0 NH:i:1 XF:Z:ENSMUSG00000064354
IL13_5283:1:100:10015:17002 163 MT 7026 255 76M = 7318 368 TCCAACTTGGTCTACAAGACGCCACATCCCCTATTATAGAAGAGCTAATAACTTTCCATGATAACACACTAATAAT 969
56(4959@<505:3086=87::31286>>962@:,5.@:2569;2>4;.579>+2<3:0((3(+:<)5;=5*A NM:i:2 NH:i:1 XF:Z:ENSMUSG00000064354
IL13_5283:1:100:10015:17002 83 MT 7318 255 76M = 7026 -368 CCAATGATACTGAAGCTACGAATATACTGACTATGAAGACCTATGCTTTGATTCATATATAATCCCAACAAACGAC 8:;
=;=:87::8<:::=7;8:=:6:?<>=:<5?5::<=6.<8<=>D<789853886@+>CD@*9>57498:6B9;B NM:i:0 NH:i:1 XF:Z:ENSMUSG00000064354
All suggestions are very much appreciated!
Comment