given a FASTQ file and a BAM file of aligned reads, is there an efficient way to get all FASTQ reads that are in the original FASTQ but not in the BAM? Perhaps using bedtools. i.e.:
unmapped_script original.fastq aligned.bam > unmapped.fastq
should create an unmapped.fastq file, which is a subset of original.fastq containing only those entries that do not appear in aligned.bam
thank you.
4 answers
If you sort the mapped BAM names and the read-names, you can use linux join to print out unpaired--e.g. unmapped reads.
FQ=$1
BAM=$2
# print out mapped headers
samtools view -F 4 $BAM | awk '{ print "@"$1"/1" }' | sort -u > mapped.txt
# grap only the fq headers.
awk '(NR % 4 == 1)' $FQ | sort -u > reads.txt
# join and print only un-paired (-v)
join -v 1 reads.txt mapped.txt > unmapped.txt
call like::
bash unmapped.sh input.fq mapped.bam
and the unmapped reads will be in unmapped.txt, from there, you can use grep -f to get back a fastq file. Note, you may have to change the awk command (currently '{ print "@"$1"/1" }') to either not add the trailing "/1" or to add "/2" depending on the read end (or depending on the samtools flags.
Get the query names from the BAM file:
samtools view <input file> | cut -f1 | sort > BAM_headers.txt
Get the query names from your FASTQ file (assuming your read length is 101bp):
awk '$0 ~ /^@/ && length($0) < 101' <fastq_file> | sed 's/^@//' | sort > FASTQ_headers.txt
Use diff to compare the two files.
diff BAM_headers.txt FASTQ_headers.txt
Since all reads in the FASTQ headers file will also be in the BAM headers file, diff should only show query names that are in the FASTQ file that are not in the BAM file. The output should show something like this:
7a8,9
> D74RYQN1_0189:1:2208:17093:200915#0
> D74RYQN1_0189:1:2208:17337:200823#0
8a11,12
> D74RYQN1_0189:1:2208:19283:200878#0
> D74RYQN1_0189:1:2208:19604:200869#0
That's your list of missing read names. You can parse the diff output to get the read names without leading '>' symbols, then use grep to get the actual reads from the FASTQ file if you'd like. If you want to use this approach and need help figuring out the last two parts, let me know.
•
written
by
Mitch Bekritsky
I do like your approach by the way, but it will have to take into account the trailing / or # signs that are used for paired end reads, since many alignment programs drop these. for example tophat will drop the /1 or /2 of the FASTQ read mate 1 and read mate 2 and only report the name that comes before the '/'. it's quite annoying
For rlong's approach to work, you would need to sort the FASTQ headers as well. Otherwise, you'd have to store the BAM headers as a dictionary or a hash, then as you read the FASTQ file, check to see if it occurs in the dictionary. This would be extremely memory inefficient, since you'd have to hold all the BAM files in memory. Sadly, my python expertise is limited. All the parsing is definitely doable in python, but I don't know if there's an implementation of diff in python. If the output of TopHat chops off the '/' or '#', you can use sed to get rid of them. See next comment for that
sed 's/[#/].*$//' FASTQ_headers.txt (or a python regular expression equivalent) would get rid of any # or / symbols in the read. Be careful, because if they occur anywhere else in the read (they shouldn't), then it will truncate from there. Also, in my previous comment, I meant that rlong's approach would have to hold all the BAM query names in memory, not the whole file.
One more thing...will tophat output reads where only one paired end maps?
I like this solution. You only traverse each input file once, limiting IO. Then if you UNKNOWN cares about the contents of the reads, and not just the read names, you have final traversal of the fastq to drop out the reads that match your list of those in the fastq but not the bam.
Thanks! It comes from long hours trying to find BASH commands that do things I need without me needing to do lots of coding.
Great to know @google. If you were mapping paired-end reads together, there's another bash trick that would tell you if one or both reads didn't map, but it's not necessary here.
The sort for reads ending in a '/1' isn't the problem, from my view. It's the string comparison with the BAM file, which (from what I've seen) won't contain the trailing '/1', and so diff will not be able to match any read headers between the FASTQ and BAM files
Unfortunately, yes; for large files, sorting can take a bit of time. Depending on the sorting algorithm used and how many entries you're trying to sort, 10 minutes isn't terrible, even if it seems like it is ;)
My pleasure! Happy to help...
Are you sorting in python or on the command line in bash? It sounds like something is out of space. Either it's your hard drive (check with df if you're on unix), or the sort is running out of RAM when it loads all the read headers into memory. If your hard drive still has lots of space, your best bet is to go with command line sort or samtools sort. When they encounter a file that's too large, they'll split the file into smaller chunks, sort the smaller chunks, then merge them together. It would be a pain to code that yourself if you're trying to get it done on a Sunday night.
that error actually looks like it's when the command line sort is running out of space. How big are your header files from the bam or fastq?
From the samtools man page, sort -n sorts by read header, which should give the same output. When you say it doesn't work, does the file not look sorted? Or is there some other issue?
Example:
$ cat test.fq | perl mergelines.pl | sort --stable -k1,1 | perl splitlines.pl > sorted
$ head sorted
@42EBKAAXX090828:6:100:1699:328/2
@42EBKAAXX090828:6:10:1077:1883/2
as you can see it thinks that "@42EBKAAXX090828:6:100:1699:328/2" comes before "@42EBKAAXX090828:6:10:1077:1883/2"
It's an alphanumeric sort, so it is treating the query names as ASCII values. The sort you're seeing is because 42EBKAAXX090828:6:100:1699:328/2 has a '0' in the same position that 42EBKAAXX090828:6:10:1077:1883/2 has a ':', and ':' comes after '0'. Samtools appears to split the query name by ':' and sort, which will be different than unix sort output. Use this command: sort -t ':' -k 1,1 -k2,2n -k3,3n -k4,4n -k 5,5 (split by ':', numeric sort on fields 2,3,4, alphanumeric sort on 1 and 5) if you want to use samtools sort, or you can just use unix sort on both files.
I looked in my local unix sort and in samtools sort and I don't see -V anywhere. What sort binary are you using that has that option?
Follow this "seqanswers" thread to extract unmapped fastqs using original Fastq and bam file (especially from output from tophat or other programms which does not include unmapped reads in bam file)
to extract unmapped reads from bam file (from newer versions of bwa etc...) which includes unmapped reads follow this.
One approach would be to sort your bam by read name:
samtools sort -n your.bam
[EDIT: thanks to unknown] The fastq will also need to be sorted by read_names for this to work.
Open a samtools view of the bam in a pipe, and step through the fastq file record by record. If the read-name of your current record sorts to before the current read_name in your samtools view pipe, then that fastq line is not in your bam. If they match, you have a read that exists in both files. If the the read_name from your fastq file sorts to after the current read_name in your samtools view pipe, you have reads in your bam and not in your fastq.
I'm not sure about unix sort. Indeed it would have to mimic the sorting of samtools sort. Regarding the repetitions in the bam, that shouldn't matter. Let's say you have a read name from your fastq, for simplicity we'll say it is "aab". So if you have two "aab" reads in your bam, which you certainly will if this is paired end data, then when you hit "aab" in your fastq, three things can happen. If the last read name from the bam is "aaa", grab another read from the bam pipe. If this is also "aaa", draw another. If you get "aac", there are no "aab" reads in the bam. --- continued ---
if you get "aab" from the bam, you grab another read from the fastq. Suppose it is "aac". Then you will compare "aab" with "aac" and draw another read from the bam pipe. If the next read is "aab" as well, no problem, you'll keep drawing them until you get "aac" or something beyond it.
Log in to answer this question.
-
Sam To Bam - Loss Of Data Or Just Great Compression?
written by Clare
-
Conversion From Sam To Bed With Max Information
written by roll
-
Samtools View Of A Sam File Entire Chromosome
written by dfernan
-
Samtools Merge/Reheader Of Bam Files With Different @Sq Lines
written by Nathan S. Watson-Haigh
-
Bedtools Multicov Need A Bam Index File Specification Option
written by ewre
-
How To Rearrange Paired End Bam File?
written by deepthithomaskannan
-
Adding Custom Bed And Bam Tracks To UCSC Browser By URL
written by user
-
Filtering Stream Bam Output Idealy From Pysam/Python
written by User 9996
-
Reporting The Bam Reads Overlapping A Set Of Intervals With Bedtools
written by User 9996
-
Fetching Read By Its Id From A Bam File With Pysam In Python
written by User 9996
Every BAM file I've seen contains all reads, including unmapped ones. Are you sure that unmapped reads aren't present in your BAM file? Or am I misundersatnding
bwa definately included unmapped reads, I'm not sure bowtie does by default.
These BAMs are produced by Tophat. My understanding is that Tophat only reports aligned reads from FASTQ, isn't that right?
I've mainly used BWA, so that's been my experience. I just checked bowtie, and if you run it as /path/to/bowtie --sam <genome> <input_file>, all reads are still reported, including unmapped reads (FLAG is 4, RNAME is *, POS and MAPQ are 0)
I've never used TopHat, but a quick peek online makes it look like they only report mapped reads. First time I've seen that...
TopHat only puts in the mapped reads - it'll put them in multiple times if they map to many locations as well - so you may need to pipe column 1 into sort/uniq UNIX pipes to get unique mapped read names.