tophat : sam-flag 115 = properly-paired + read.reverse + mate.reverse ?

[Pierre Lindenbaum]

cross posted on biostars: http://www.biostars.org/p/87071 and http://seqanswers.com/forums/showthread.php?t=36055


I ran tophat2 using the standard options.

Code:

    $ tophat2 -v                                                                    
    TopHat v2.0.10

Code:

    $ samtools view -H TOPHAT/accepted_hits.bam   | grep PG
    @PG     ID:TopHat       VN:2.0.10       CL:/commun/data/packages/tophat-2.0.10.Linux_x86_64/tophat -p 10 -G genes.gtf -o  TOPHAT --rg-id g24 --rg-library 6VGWT3 --rg-sample 6VGWT3 --rg-description 6VGWT3  34 fastqs --rg-platform-unit 1 2 3 4 --rg-center Nantes --rg-platform Illumina  mm10 6VGWT3_ATGTCA_L002_R1_002.fastq.gz,6VGWT3_ATGTCA_L004_R1_002.fastq.gz,6VGWT3_ATGTCA_L003_R1_002.fastq.gz,6VGWT3_ATGTCA_L002_R1_003.fastq.gz,6VGWT3_ATGTCA_L004_R1_003.fastq.gz,6VGWT3_ATGTCA_L003_R1_003.fastq.gz,6VGWT3_ATGTCA_L002_R1_004.fastq.gz,6VGWT3_ATGTCA_L004_R1_004.fastq.gz,6VGWT3_ATGTCA_L003_R1_004.fastq.gz,6VGWT3_ATGTCA_L004_R1_001.fastq.gz,6VGWT3_ATGTCA_L002_R1_001.fastq.gz,6VGWT3_ATGTCA_L003_R1_001.fastq.gz,6VGWT3_ATGTCA_L001_R1_002.fastq.gz,6VGWT3_ATGTCA_L001_R1_003.fastq.gz,6VGWT3_ATGTCA_L001_R1_004.fastq.gz,6VGWT3_ATGTCA_L001_R1_001.fastq.gz 6VGWT3_ATGTCA_L002_R2_002.fastq.gz,6VGWT3_ATGTCA_L004_R2_002.fastq.gz,6VGWT3_ATGTCA_L003_R2_002.fastq.gz,6VGWT3_ATGTCA_L002_R2_003.fastq.gz,6VGWT3_ATGTCA_L004_R2_003.fastq.gz,6VGWT3_ATGTCA_L003_R2_003.fastq.gz,6VGWT3_ATGTCA_L002_R2_004.fastq.gz,6VGWT3_ATGTCA_L004_R2_004.fastq.gz,6VGWT3_ATGTCA_L003_R2_004.fastq.gz,6VGWT3_ATGTCA_L004_R2_001.fastq.gz,6VGWT3_ATGTCA_L002_R2_001.fastq.gz,6VGWT3_ATGTCA_L003_R2_001.fastq.gz,6VGWT3_ATGTCA_L001_R2_002.fastq.gz,6VGWT3_ATGTCA_L001_R2_003.fastq.gz,6VGWT3_ATGTCA_L001_R2_004.fastq.gz,6VGWT3_ATGTCA_L001_R2_001.fastq.gz

I found some sam flags= 115 !

115=

• read paired
• read mapped in proper pair
• read reverse strand
• mate reverse strand
• first in pair

Code:

    $ samtools view  -f 115 -F 256  dir/accepted_hits.bam | head -n 2
    HWI-1KL149:61:D2C11ACXX:4:2204:6848:94129       115     chr1    24611547        1       70M2D31M        chrM    10906   0       GTAGGCGATTAGTGATTTTAAATCTGTTTGGCGTAAGCAGATTGAGCTAGTTATAATTATTCCTCATAGGGAGAAGGATGAAGGGGTATGCTATATATTTT      DDDDDBDDDDDEEEEEEFFFFFFFHHHHJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJIJJJJIJJJJJJJJJJJJJJHHHHHFFFFFCCC      AS:i:-11        XN:i:0  XM:i:0  XO:i:1  XG:i:2  NM:i:2  MD:Z:70^GA31    YT:Z:UU NH:i:4  CC:Z:chrM CP:i:10928       HI:i:2  RG:Z:g24
    HWI-1KL149:61:D2C11ACXX:2:2109:12004:4228       115     chr1    24611549        1       101M    chrM    10815   0       TGGCGATTAGTGATTTTAAATCTGTTTGGCGTAAGCAGATTGAGCTAGTTATAATTATTCCTCATAGGGAGAGAAGGATGAAGGGGTATGCTATATATTTT      DDDDDDDDDDDEEEEEEEFFFFFHHHJJJJJJJJJJJJJJJJJJIJJJJIJJJJJJJJIJJJJIGJJJJJJJJJJIJJJJJJJJJJJJHHHHHFFFFFCCC      AS:i:-5 XN:i:0  XM:i:1  XO:i:0  XG:i:0  NM:i:1  MD:Z:0A100      YT:Z:UU NH:i:4  CC:Z:chrM       CP:i:10928HI:i:2   RG:Z:g24

how can a read be "mapped in proper pair" and read reverse strand+ mate reverse strand ? what is the consequence for a tool like htseqcount ? Does it only count the reads in proper pair ?

Thank you.


Pierre

[bert hubert]

Hi Pierre,

In Antonie I frequently see pairs that would align both in the same direction, and close by too. I think this is some kind of Illumina artifact. In Antonie, I do map both reads, but I don't pair them up. Perhaps Tophat decides that if they are close enough together they should be mapped as a pair?

     Bert

[Pierre Lindenbaum]

I frequently see pairs that would align both in the same direction

but then the sam flag "properly paired" shouldn't be set.