This is a test version of Biostars. For the public version, visit https://www.biostars.org.
remove a percent of reads from a BAM file

I would like to remove x% of reads from a BAM file at a specific region . For instance : remove 50% of reads in chr4:1-1000000

What's the most efficient approach ?

bam remove
• 4,564 views • link
written by sacha sacha

3 answers

This should work:

 java -jar GenomeAnalysisTK.jar \
   -T PrintReads \
   -R reference.fasta \
   -I input.bam \
   -o output.bam \
   -dfrac 0.5 \
   -L chr4:1-1000000

https://software.broadinstitute.org/gatk/documentation/tooldocs/org_broadinstitute_gatk_tools_walkers_readutils_PrintReads.php

• 1 views • link
written by Zaag Zaag

I prefer something using samtools / bedtools ! But thanks, I will try it

• 1 views • link
written by sacha sacha

Maybe something like this:

samtools view -sb 0.5 in.bam chr4:1-1000000 > out.bam
• 1 views • link
written by Zaag Zaag

Those commands keeps reads . I want to exclude them. I though to remove all reads in chr4:1-1000000, and then merge with your output. Unfortunally, I didn't find how to remove easily reads in the region chr4:1-1000000

• 1 views • link
written by sacha sacha

This removes all the reads in the specified region (-XL option) to include use L, don't know a samtools equivalent:

 java -jar GenomeAnalysisTK.jar \
   -T PrintReads \
   -R reference.fasta \
   -I input.bam \
   -o output.bam \
   -XL chr4:1-1000000
• 1 views • link
written by Zaag Zaag

Thanks ! it seems working. But I wonder why the reference is required.

• 1 views • link
written by sacha sacha

Similarly, you can do this with Picard Tools 

java -jar ~/tools/picard.jar DownsampleSam I=Input.bam o=Output.bam PROBABILITY=0.5

Or use use samtools view

samtools view -s 0.5 Input.bam > Output.bam

Picardtools has the advantage that it keeps or discards read pairs as a whole. E.g. if one mate is removed/kept, the other is as well.

• 1 views • link
written by stam • stam

Seems those commands keeps reads. I want to exclude reads from a region chr4:1-1000000.

• 1 views • link
written by sacha sacha

Randomly drawing 50% of the reads from a region is the same as randomly discarding 50% of the reads, isn't it?

Both samtools and picard allow you to specify the region as well.

• 1 views • link
written by stam • stam

Attached a picture of what I want :

https://snag.gy/Hyb4Cu.jpg

• 1 views • link
written by sacha sacha

To remove reads from a bam in a particular region, I only found the python way using pysam. Here is the pseudo-code

for read in bamfile:
     if read.pos > region_start and (read.pos + read.length) < region_end:
          do_nothing()
     else:
        output.write(read)
• 1 views • link
written by sacha sacha

My version of pysam doesn't have read.length - what is it? If you meant len(read.seq), strictly speaking you should parse the CIGAR string as that takes into consideration insertions/deletions/skips/softclipped bases/etc. Also, obviously, read pairs will be all screwed up as stam says, but perhaps thats OK for your application.

• 1 views • link
written by John John

I m using : read.query_length

• 1 views • link
written by sacha sacha

Yeah, I mean, you can do that by all means, but it won't be as accurate as parsing the CIGAR strings. To be honest, for your tests it'll probably be fine. But this is what i'm talking about:

>>> b
<pysam.calignedsegment.AlignedSegment object at 0x109a18f58>
>>> print b
SOLEXA1_0001:4:87:13822:5416#0  16  0   3299480 0   8M16308N28M -1  -1  36  GAAACTAACACTCTCTCTCTCTCTCTCTCTCTCTCT    array('B', [34, 34, 34, 34, 34, 34, 34, 34, 34, 33, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34])    [('NM', 2), ('XS', '-'), ('NH', 5), ('CC', 'chr10'), ('CP', 85523284)]
>>> print b.query_length
36

So this is a read where, according to the CIGAR, the first 8bp map, then theres a 16,308bp gap, then 28bp map. To get the end position you want the POS + 8 + 16308 + 28, and not POS + 36. However, since reads like this can be rare, you might want to just add code that checks if the read has a weird CIGAR and ask you manually to accept/reject the read. Might be faster in the long run.

I also don't know how pysam handles soft-clipped reads (and for the record, this is exactly why i dislike pysam - it causes more ambiguity than it gives back in abstraction). Perhaps query_length, whatever that is, takes into account soft-clipped reads. Perhaps not. Who knows. But if it does not and you have soft-clipped adapters, you'll have to parse CIGAR strings for sure.

• 1 views • link
written by John John

Log in to answer this question.

More posts like this