I have fasta files with multiple sequences. They are reads that mapped to the genome outside of probe-targeted regions. From a quick perusal, they appear to be repetitive and have low complexity. Is there a way I can quantify their complexity and use that as an explanation as to why those regions were mapped to?
3 Answers
"Complexity" is also known as entropy. A quick and dirty way is to use a file compression utility - zip will get good compression for a repetitive sequence.
You can also look at the k-mer distribution for some suitable k: in a repetitive sequence it will be biased. Here is a tool that give a local plot: https://rdrr.io/github/vsbuffalo/qrqc/man/kmerEntropyPlot.html
-
$\begingroup$ I looked into compressibility as a measure of complexity. I found this article very interesting. They do not recommend it for DNA sequences. sci-hub.do/10.1109/iwobi.2014.6913941 $\endgroup$ Mar 2, 2021 at 8:51
Use Picard EstimateLibraryComplexity.
From the docs:
Estimates the numbers of unique molecules in a sequencing library.
java -jar picard.jar EstimateLibraryComplexity \
I=input.bam \
O=est_lib_complex_metrics.txt
-
1$\begingroup$ Although I mentioned I had fasta files, I was able to take the original bam files and pipe only those regions into picard with ´samtools view -b -L regions.bed sample.bed | picard....´ $\endgroup$ Mar 2, 2021 at 9:24
How long are your reads, and how many? Some colleagues of mine wrote https://github.com/eclarke/komplexity to quickly score complexity in large numbers of short reads. The score is derived from the number of unique k-mers normalized by the sequence length, so it's simple and fast, though the score inherently becomes less informative for longer and longer sequences (so, best for short reads). We've used it for filtering repetitive or other low-complexity sequences in metagenomics projects.
-
$\begingroup$ This looks promising! I don't want to run it on the unmapped per se as you can see in my question. But since all of the non-targeted regions are less than 300 bps, I think komplixity will be helpful. I'll let you know how it goes, and possibly award this as the accepted answer. $\endgroup$ Mar 5, 2021 at 8:54
-
$\begingroup$ do you think you could ask your colleagues to take a look at their issues on github? I'm not sure it's being maintained. $\endgroup$ Mar 9, 2021 at 13:14
-
$\begingroup$ I see what you mean, those issues aren't getting any attention. I think the short answer for most of those is that it'd need new features to address them, like handling paired-end reads automatically. (But for that use case you can at least filter forward and reverse separately, and take the intersection of what remains.) If I can carve out the time I'll send eclarke some PRs to help with those though. $\endgroup$– JesseMar 9, 2021 at 16:31
-