Extract reads from bam/sam files using read id

Question

[Chirag Parsania]

Hi,

I have fastq header of 1000s of reads. I want to check presence all of those reads in another bam file. I have written following script

while IFS='' read -r line || [[ -n "$line" ]]; do samtools view accepted_hits.bam | grep $line ;done < fastqHeaders.txt ;

But this is taking too much time as my .bam file is too big. does anyone know faster way to do this.

Kindly help

Chirag

Answer 1

[Pierre Lindenbaum]

use picard FilterSamReads: https://broadinstitute.github.io/picard/command-line-overview.html

READ_LIST_FILE (File) Read List File containing reads that will be included or excluded from the OUTPUT SAM or BAM file. Default value: null.

[Carlo Yague]

After trying to play with grep and a sam file (it worked but was so sloooow), I found your answer and its much faster, thx.

[max_19]

any idea why FilterSamReads would result in

"ERROR: Invalid argument ' '.

code:

> java -jar $EBROOTPICARD/picard.jar FilterSamReads \ I=usritt.bam \ 
> O=output.bam \ READ_LIST_FILE=keep.txt \ FILTER=includeReadList

[Carlo Yague]

try to remove the "\" except those at the end of the lines. For instance:

java -jar $EBROOTPICARD/picard.jar FilterSamReads I=usritt.bam O=output.bam READ_LIST_FILE=keep.txt FILTER=includeReadList

Answer 2

[dariober]

That's how I would do it. I don't know how fast it is, but for a bam file of few 10s of millions of reads should be acceptable. It depends on python's pysam library. File fastqHeader.txt should have one read name per line.

!#/usr/bin/end python

import pysam

fq= open('fastqHeader.txt').readlines()
fq= [x.strip() for x in fq]
fq= set(fq)

infile= pysam.AlignmentFile('in.bam')
outfile= pysam.AlignmentFile('out.bam', template= infile, mode= 'wb')
for aln in samfile:
    if aln.query_name in fq:
        outfile.write(aln)

infile.close()
outfile.close()

[sidwekhande]

Thanks, this helped me!

For future users, it should be

for aln in infile: #in place of samfile

[GA]

Hi dariober! Can you post a sample of what the fastqHeader.txt is supposed to look like? Thank you so much!!