Unconfigured Ad

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

    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
            Proteomic Platforms: How to Choose the Right Analytical Strategy to Improve Detection and Clinical Applications
            by SEQadmin2


            Proteomics platforms are evolving rapidly, with advances in mass spectrometry and affinity-based approaches expanding what researchers can detect and at what scale. As the field moves toward deeper proteome coverage and clinical applications, scientists face an increasingly complex landscape of tools. This article will explore how researchers are navigating these choices to find the right platform for their work.

            The systematic characterization of the human proteome has
            ...
            Yesterday, 11:48 AM
          • SEQadmin2
            Advanced Sequencing Platforms Tackle Neuroscience’s Toughest Genomics Problems
            by SEQadmin2



            Genomics studies in neuroscience face a special challenge due to the brain’s complexity and scarcity of samples. Mapping changes in cell type and state using conventional next-generation sequencing methods remains challenging. Advances in technologies like single-cell sequencing, spatial transcriptomics, and long-read sequencing have opened the door to deeper studies of the brain and diseases like Alzheimer’s, amyotrophic lateral sclerosis (ALS), and schizophrenia.
            ...
            07-09-2026, 11:10 AM
          • SEQadmin2
            Cancer Drug Resistance: The Lingering Barrier to Rising Survival
            by SEQadmin2



            Cancer survival rates have significantly increased in the last few decades in the United States, reaching a combined 70% 5-year survival rate by 2021. Behind this number, there are years of research to find new therapies, drug targets, and early detection methods. But there is one core challenge that keeps slowing down these advances, and it’s about drug resistance.

            There is no single reason why many patients don’t respond to treatment as expected. Cancer is...
            07-08-2026, 05:17 AM

          ad_right_rmr

          Collapse

          News

          Collapse

          Topics Statistics Last Post
          Started by SEQadmin2, Yesterday, 11:10 AM
          0 responses
          9 views
          0 reactions
          Last Post SEQadmin2  
          Started by SEQadmin2, 07-13-2026, 10:26 AM
          0 responses
          30 views
          0 reactions
          Last Post SEQadmin2  
          Started by SEQadmin2, 07-09-2026, 10:04 AM
          0 responses
          40 views
          0 reactions
          Last Post SEQadmin2  
          Started by SEQadmin2, 07-08-2026, 10:08 AM
          0 responses
          25 views
          0 reactions
          Last Post SEQadmin2  
          Working...