SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
BWA MEM - which reads were mapped by mate-pair rescue? Flo1000 Bioinformatics 4 09-09-2014 08:01 AM
BWA MEM mate pair incorrect insert sizes? hartmaier Bioinformatics 7 07-24-2014 11:38 AM
how to set parameters for mapping reads to a bacterial genome using "bwa-mem/bowtie"? joyjane88 RNA Sequencing 1 11-13-2013 02:54 AM
Samtools Flag to identify outward facing (mate-pair not paired end) reads lesander Bioinformatics 5 10-28-2011 01:24 PM

Reply
 
Thread Tools
Old 07-05-2015, 06:42 PM   #1
gada
Junior Member
 
Location: china

Join Date: Dec 2014
Posts: 3
Default 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.
gada is offline   Reply With Quote
Old 07-06-2015, 05:25 AM   #2
dariober
Senior Member
 
Location: Cambridge, UK

Join Date: May 2010
Posts: 310
Default

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...
dariober is offline   Reply With Quote
Old 07-06-2015, 05:32 AM   #3
gada
Junior Member
 
Location: china

Join Date: Dec 2014
Posts: 3
Default

thanks. got it. Maybe bwa mem just like paired reads overlapped? I'll have a try.
gada is offline   Reply With Quote
Old 07-06-2015, 06:02 AM   #4
dariober
Senior Member
 
Location: Cambridge, UK

Join Date: May 2010
Posts: 310
Default

Quote:
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.
dariober is offline   Reply With Quote
Old 07-06-2015, 06:19 AM   #5
gada
Junior Member
 
Location: china

Join Date: Dec 2014
Posts: 3
Default

$ 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.
gada is offline   Reply With Quote
Reply

Tags
bwa, samtools

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 02:13 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2017, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO