Hi Guys,
I have a BAM file, and a big read list. What I want to do is to remove the reads in the read list from the BAM file. I can transform Bam to Sam file and then use a Python script to remove unwanted reads. And then transform Sam to Bam again. But I am wondering if there is a more efficient way, which I mean faster, easier, and memory-efficient, to achieve this goal?
Any advice is appreciated!
Tao
•
written
by
Tao
5 answers
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
samtools view -h sample1.bam | grep -vf read_ids_to_remove.txt | samtools view -bS -o sample1_filter.bam -
http://genometoolbox.blogspot.com/2013/06/remove-list-of-reads-from-bam-file.html
Hi Vivek,
Thanks for your nice reply. I also found this method, but unfortunately, it's kind of some inefficient. I think may be 'grep' is always doing global search. Here, given a read ID, grep has to confirm all the bam context doesn't contain this read ID(pattern). So it's much slower when a big list of unwanted reads occur. While in my script, I only need to confirm the first column(reads ID) of the Bam file doesn't contain those unwanted reads. But, anyway, your method is a good choice to handle small read list. Thanks!
Tao
The fastest and easiest solution is probably to use BBMap + samtools:
filterbyname.sh in=mapped.bam out=filtered.bam names=names.txt include=false
Samtools needs to be in the path. The memory usage depends on the number of names; the speed doesn't (well, not much).
•
written
by
Brian Bushnell
Indeed, this works:
samtools view -h sample1.bam | grep -vf read_ids_to_remove.txt | samtools view -bS -o sample1_filter.bam -
If possible, adding --threads to samtools view is advisable. \ Using awk instead of grep should be faster.
samtools view ... | awk 'FNR==NR{a[$1];next} {if(!($1 in a)) {print $0 }}' read_ids_to_remove.txt - | samtools view ...
not
grep -vf read_ids_to_remove.txt
but:
LC_ALL=C grep -vFwf read_ids_to_remove.txt
Since version 1.12, samtools has an option -N / --qname-file for filtering on the read name field:
samtools view -N read_ids_to_remove.txt -U filtered.bam -o /dev/null input.bam
(Add threading -@ / --threads options to taste, as described in other answers.)
Log in to answer this question.
-
Specific case of "bedtools intersect"
written by predeus
-
Efficiently merge two BAM files while retaining reads from only one file in overlapping regions
written by kiranchari •
-
Can I do eqtl mapping combining RNA-seq and Microarray data?
written by Tao
-
How to remove a list of reads from fastq file?
written by Tao
-
Should I use reads with good quality but failed-vendor flag?
written by Tao
-
How to count files/folders on NCBI ftp server?
written by Tao
-
Extract data using awk/sed and output to different files
written by Tao
-
Filtering nth length from mapped SAM/BAM file.
written by unique379
-
Converting Sam File To Hdf5
written by Abhi
-
Fetching All Alignments From A Sam/Bam By Read Header In Perl
written by Abhi
