I have created a denovo reference "genome" of RAD contigs and have mapped back my reads using BWA. Using the the samtools flagstat option I queried my read alignments. An example:
3232117 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
3226784 + 0 mapped (99.84% : N/A)
3232117 + 0 paired in sequencing
1595145 + 0 read1
1636972 + 0 read2
2662532 + 0 properly paired (82.38% : N/A)
2871947 + 0 with itself and mate mapped
354837 + 0 singletons (10.98% : N/A)
213393 + 0 with mate mapped to a different chr
213393 + 0 with mate mapped to a different chr (mapQ>=5)
Two things are obvious from the alignment: 1) singletons must arise because a mate fails the quality check during the mapping procedure, and 2) in some cases mates map to different "chromosomes" (RAD contigs).
So the questions.
Firstly, I feel like eliminating singletons would essentially be throwing out information. Should these be kept or is there too much uncertainty in their origin to be reliable because their mate hasn't mapped?
Secondly (and more importantly), how do I remove reads in which mate pairs have mapped to different chromosomes? I have played with the samtools view -f 3 option (read paired + read mapped in proper pair) and this reduces the split mates but they are still present; it also removes singletons. From what I gather, there doesn't seem to be a specific option combination to do this (https://broadinstitute.github.io/picard/explain-flags.html). NOTE: I have also tried using -f 11 (read paired + read mapped in proper pair + mate unmapped), but this removes all reads (I think because only reads that are properly paired are allowed make it through, and thus cannot have an unmapped mate).
A follow on question is should such pairs be removed, given that they might be indicative of repetitive sequences?
Thoughts and comments?
These pairs might help you identify mis-assemblies [or even sites for recent transposon mobilization] thus there is no reason to remove them.
Singletons should be retained as single reads only.
It might help if you specify 'why' you want to remove the reads.
@David,
Thanks for the answer. :)
Do you have any tips as to how to treat the singletons as unpaired reads in analysis that also includes paired reads in the same sample?
I guess reads whose mates map to diffrent chromosomes make me nervous in that they don't have an undispurted origin. Also, because they are relatively low in frequency I reasoned that I wouldn't lose too much data by throwing these out.
Still learning about "best practises" for genomic data, so any insights are really useful.
First of all, BWA is antiquated. Use this software, it is literally a step ahead of the rest: https://ccb.jhu.edu/software/hisat2/index.shtml
But to answer your question, most read alignment softwares have this option (to map reads as singletons). For BWA, apparently you have to use something like this: bwa samse
Again, I would recommend to use Hisat2. See the help menu for options to do all the things you asked about and more.
@David.
Thanks for suggestion. Will definitely have a play with this program.
- Badges
- Science method