Hi all, Hope someone can help.

We can pass a region to samtools view such as 1:1000-1010.

Is it possible to pass a list of regions to samtools view and have the output into a single sam file i.e. (something like)

$ samtools view -h .bam regionlist.txt > mysam.sam

what would the syntax for this be?

also, in what format does the region list file should be in?

thanks, a

Samtools have the -L option

samtools view -b -L ROI.bed  file.bam > ROI_file.bam

BUT -L does not use the samtools index so the search is slooow depending on how large is your bam file. In my benchmarks querying for 10, 100 or 1000 sequences take exactly the same time so seem that bed size does not matter too much. You need to do a bit of benchmarking if you prefer to do a -L or to do a loop (for a bam with 23 million reads take the same time the -L that writting all the stuff for looping ;-)) YMMV

$ time samtools view -b -L 100_ROI.bed file.bam > ROI.bam
real    0m38.831s
user    0m37.666s
sys    0m0.556s


$ time (samtools view -H file.bam > roi_xargs.sam; \
      cat 100_ROI.bed | perl -lane 'print "$F[0]:$F[1]-$F[2]"' | xargs -n1 -t -I{} samtools view file.bam {} >> roi_xargs.sam; \
      samtools view -bSh roi_xargs.sam > roi_xargs.bam \
  )
real   0m7.188s
user   0m5.080s
sys    0m1.304s

And if you feel adventurous and your disks are using isilon or lustre you can use gnu-parallel and forget about same file concurrency and do it in parallel in a breeze. The only issue is that you would need to do the bam merge afterwards.

From the samtools documentation:

 -L FILE  output alignments overlapping the input BED FILE [null]

samtools view -L regionlist.bed -h in.bam > out.sam

The region list file would be in BED format.

You can also use tabix, which is an index and search tool developed by the Samtools author. Here is an example to fit your needs:

Your SAM file should be coordinate sorted. If it is not, then use gnu sort.

sort -k 3,3 -k 4,4n mysam.sam > mysam.sorted.sam

Then, tabix needs a bgzip compressed file.

bgzip mysam.sorted.sam

Then you need tabix to index your file. This will create a file ending in .tbi

tabix -p sam mysam.sorted.sam.gz

Lastly, you may retrieve your regions in this format:

tabix mysam.sorted.sam.gz chr1:1000-1010 chr2:2000-2020 ... > mysam.sorted.regions.sam

You could use a for loop in a shell script to do this.

# ouput header only, once
samtools view -H aln.bam > mysam.sam
while read line;do
    samtools view aln.bam $line >> mysam.sam
done < regions.txt

where regions.txt is line-based, with each line being:

ref1:start-end

or xargs:

samtools view -H aln.bam > mysam.sam
xargs -a regions.txt -I {} samtools view aln.bam {} >> mysam.sam

Sorry for the bump (it's for future Google visitors), I used cl.parallel's answer as slight guidance but I wanted to get the same region across a set of reads, so I wrote a script. Uses a BED file as input.

for i in $( ls *.bam ); do
    samtools view $i -H > OUTPUTPREFIX${i%.bam}.sam 
    while read -r -a myArray; do 
        samtools view $i ${myArray[0]}:${myArray[1]}-${myArray[2]} >> OUTPUTPREFIX${i%.bam}.sam
     done < INPUTFILE.bed 
done

This is my solution so far and hope it helps. The key to the question is: samtools view accept multiple position. This is useful for small files.

samtools view -O BAM -o OUTPUT.bam INPUT.bam `perl -ane '{print "$F[0]:$F[1]-$F[2] "}' INPUT.bed `

where INPUT.bam is your BAM file

OUTPUT.bam is your output.

INPUT.bed is a file with the format of chr, start and end.