Unconfigured Ad

Collapse
X
 
  • Time
  • Show
Clear All
new posts
  • gada
    Junior Member
    • Dec 2014
    • 3

    #1

    bwa mem don't set PROPER_PAIR flag in pair end align

    $ cat r1.fq
    @1
    GGGCTAAGAGACAGCAAATACACATGAACAGAAAGAAGAGGGCAAAGAAAAGGCTGACGG
    +
    CCCFFFFFHHFHHJIJIJIJJJJJJJJDIJJJJJIJJIJJIIGIJJJJHIHJJJJIHIHI

    $ cat r2.fq
    @1
    TAGTTTAAACGAGATTGCCAGGACCGGGTATCATTCACCATTTTTCTTTTCGTTAACTTG
    +
    CCCFFFFFHHHHHJJJIIIHIIJGGIGIJEIJIFGGGIJJJGFIIIJIFFEF?CFFEDCB

    $ bwa mem ../ref/ref.fa r1.fq r2.fq > mem.sam;cat mem.sam
    @SQ SN:1 LN:100000
    @PG ID:bwa PN:bwa VN:0.7.10-r1017-dirty CL:bwa mem ../ref/ref.fa r1.fq r2.fq
    1 81 1 1741 60 60M = 1681 -120 CCGTCAGCCTTTTCTTTGCCCTCTTCTTTCTGTTCATGTGTATTTGCTGTCTCTTAGCCC IHIHIJJJJHIHJJJJIGIIJJIJJIJJJJJIDJJJJJJJJIJIJIJHHFHHFFFFFCCC NM:i:1 MD:Z:18A41 AS:i:55 XS:i:0
    1 161 1 1681 60 60M = 1741 120 TAGTTTAAACGAGATTGCCAGGACCGGGTATCATTCACCATTTTTCTTTTCGTTAACTTG CCCFFFFFHHHHHJJJIIIHIIJGGIGIJEIJIFGGGIJJJGFIIIJIFFEF?CFFEDCB NM:i:1 MD:Z:21C38 AS:i:55 XS:i:0

    $ bwa aln ../ref/ref.fa r1.fq >r1.sai
    $ bwa aln ../ref/ref.fa r2.fq >r2.sai
    $ bwa sampe -r "@RG\tID:foo\tSM:bar" ../ref/ref.fa r1.sai r2.sai r1.fq r2.fq |\
    samtools view -bS - > aln.bam
    $ samtools view aln.bam
    1 83 1 1741 60 60M = 1681 -120 CCGTCAGCCTTTTCTTTGCCCTCTTCTTTCTGTTCATGTGTATTTGCTGTCTCTTAGCCC IHIHIJJJJHIHJJJJIGIIJJIJJIJJJJJIDJJJJJJJJIJIJIJHHFHHFFFFFCCC RG:Z:foo XT:A:U NM:i:1 SM:i:37 AM:i:37 X0:i:1 X1:i:0 XM:i:1 XO:i:0 XG:i:0 MD:Z:18A41
    1 163 1 1681 60 60M = 1741 120 TAGTTTAAACGAGATTGCCAGGACCGGGTATCATTCACCATTTTTCTTTTCGTTAACTTG CCCFFFFFHHHHHJJJIIIHIIJGGIGIJEIJIFGGGIJJJGFIIIJIFFEF?CFFEDCB RG:Z:foo XT:A:U NM:i:1 SM:i:37 AM:i:37 X0:i:1 X1:i:0 XM:i:1 XO:i:0 XG:i:0 MD:Z:21C38

    the problem is:
    samtools mpileup can't output valid results when the PROPER_PAIR flag not set.
  • dariober
    Senior Member
    • May 2010
    • 311

    #2
    If you want mpileup to include these cases use the -A option (-A: count anomalous read pairs). However, the question of way bwa mem doesn't like this pair remains...

    Comment

    • gada
      Junior Member
      • Dec 2014
      • 3

      #3
      thanks. got it. Maybe bwa mem just like paired reads overlapped? I'll have a try.

      Comment

      • dariober
        Senior Member
        • May 2010
        • 311

        #4
        Originally posted by gada View Post
        Maybe bwa mem just like paired reads overlapped?
        Mmm.. I doubt it... Non overlapping pairs make perfect sense and I don't think such a bug could have sneaked into bwa mem undetected.

        Comment

        • gada
          Junior Member
          • Dec 2014
          • 3

          #5
          $ bwa mem ../ref/ref.fa r1.fq r2.fq
          [M::bwa_idx_load_from_disk] read 0 ALT contigs
          @SQ SN:1 LN:100000
          @PG ID:bwa PN:bwa VN:0.7.10-r1017-dirty CL:bwa mem ../ref/ref.fa r1.fq r2.fq
          [M:rocess] read 2 sequences (182 bp)...
          [M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (0, 1, 0, 0)
          [M::mem_pestat] skip orientation FF as there are not enough pairs
          [M::mem_pestat] skip orientation FR as there are not enough pairs
          [M::mem_pestat] skip orientation RF as there are not enough pairs
          [M::mem_pestat] skip orientation RR as there are not enough pairs
          [M::mem_process_seqs] Processed 2 reads in 0.000 CPU sec, 0.000 real sec
          1 97 1 1861 60 120M = 1979 180 CTTCATCTGCAGGTGTCTGACTTCCAGCAACTGCTGGCCTGTGCCAGGGTGCAAGCTGAGCACTGGAGTGGAGTTTTCCTGTGGAGAGGAGCCATGCCTAGAGTGGGATGGGCCATTGTT CCCCFFFFFHHFHHJIJIJIJJJJJJJJDIJJJJJIJJIJJIIGIJJJJHIHJJJJIHIHICCFFFFFHHFHHJIJIJIJJJJJJJJDIJJJJJIJJIJJIIGIJJJJHIHJJJJIHIHI NM:i:0 MD:Z:120 AS:i:120 XS:i:0
          1 145 1 1979 60 62M = 1861 -180 TTCATCTTCTGGCCCCTGTTGTCTGCATGTAACTTAATACCACAACCAGGCATAGGGGAAAG BBBCDEFFC?FEFFIJIIIFGJJJIGGGFIJIEJIGIGGJIIHIIIJJJHHHHHFFFFFCCC NM:i:0 MD:Z:62 AS:i:62 XS:i:0
          [main] Version: 0.7.10-r1017-dirty
          [main] CMD: bwa mem ../ref/ref.fa r1.fq r2.fq
          [main] Real time: 0.019 sec; CPU: 0.004 sec

          $ ruby -e 'puts 97.to_s(2)'
          1100001
          ruby -e 'puts 145.to_s(2)'
          10010001

          nothing changed.

          Comment

          Latest Articles

          Collapse

          • SEQadmin2
            Multiomic Liquid Biopsy: What DNA/RNA Panels Reveal That DNA Alone Can’t
            by SEQadmin2



            Tumors contain complex mixtures of cell populations that evolve in response to treatment and host-specific factors. Tissue biopsies provide only static snapshots of a tumor’s state. But what if we could follow that evolution over time using a minimally invasive technique?

            Liquid biopsies can detect tumor-derived material circulating in the bloodstream, providing a way to study treatment response, detect emerging resistance, and discover new biomarkers. Circulating
            ...
            10-14-2026, 09:55 AM
          • SEQadmin2
            New Genomics Technologies Take Aim at Long-Standing Limits
            by SEQadmin2


            Researchers using sequencing and genomics tools often have to make trade-offs. They can choose between speed or scale, short reads or long-range information, or targeted panels or a view of the whole transcriptome. New technologies that have been released this year are built to address those tough choices.

            We asked six companies the same four questions to learn about their latest products. The new technologies bring a lot to the table, including rethinking sequencing
            ...
            09-28-2026, 10:25 AM

          ad_right_rmr

          Collapse

          News

          Collapse

          Topics Statistics Last Post
          Started by SEQadmin2, Yesterday, 10:52 AM
          0 responses
          12 views
          0 reactions
          Last Post SEQadmin2  
          Started by SEQadmin2, 10-14-2026, 05:39 AM
          0 responses
          19 views
          0 reactions
          Last Post SEQadmin2  
          Started by SEQadmin2, 10-13-2026, 09:56 AM
          0 responses
          19 views
          0 reactions
          Last Post SEQadmin2  
          Started by SEQadmin2, 10-09-2026, 11:14 AM
          0 responses
          27 views
          0 reactions
          Last Post SEQadmin2  
          Working...