Pierre posted an efficient method for step 2 earlier:
I have paired-end data as a BAM file and I am trying to extract reads from a specific region. This sounds pretty straightforward, however, the catch is that some of the reads may have their mates mapped in a completely different location, maybe due to large structural variation.
I was going through some of samtools' options however, I did not find any particular parameter that can get me what I want.
Any suggestions would be much appreciated.
3 answers
I think it could be done using a two-step approach.
First extract all the reads that are mapped within your region of interest using samtools view command.
samtools view Input.bam chrA:x-y | cut -f1 > IDs.txtNow you can use these read ids and extract both forward and reverse reads corresponding to that read id from the original bam file. You can use a single grep command with multiple strings or read ids. See
egrepand|. Of course , this is not an efficient way but it should work.
egrep "^#|Id1|Id2|....|Idn" Input.sam (sam format) > Output.sam or you can use grep -f IDs.txt Input.sam and add header afterwards.
I wrote something today: http://lindenb.github.io/jvarkit/SamViewWithMate.html
$ java -jar dist/samviewwithmate.jar -r "9:137230721-137230796" ./src/test/resources/HG02260.transloc.chr9.14.bam | cut -f 1-9 | tail
ERR251239.10989793 83 9 137230747 60 30S70M = 137230326 -490
ERR251239.3385449 147 9 137230754 60 1S99M = 137230352 -500
ERR251240.17111373 99 9 137230764 60 100M = 137231150 475
ERR251240.46859433 147 9 137230777 60 65S35M = 137230342 -469
ERR251240.74563730 147 9 137230787 60 1S99M = 137230407 -478
ERR251240.1291708 83 9 137230789 60 100M = 137230411 -477
ERR251240.11887757 97 9 137230795 37 100M 14 79839451 0
ERR251239.34016218 81 14 79839349 37 100M 9 137230679 0
ERR251240.10196873 81 14 79839368 37 100M 9 137230721 0
ERR251240.11887757 145 14 79839451 37 100M 9 137230795 0
Log in to answer this question.