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
•
written
by
Chirag Parsania
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.
•
written
by
Pierre Lindenbaum
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.
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
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()
Thanks, this helped me!
For future users, it should be
for aln in infile: #in place of samfile
Log in to answer this question.
