2
$\begingroup$

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?

$\endgroup$
1
$\begingroup$

"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

$\endgroup$
1
  • $\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$
    – Glubbdrubb
    Mar 2 at 8:51
1
$\begingroup$

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
$\endgroup$
1
  • 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$
    – Glubbdrubb
    Mar 2 at 9:24
0
$\begingroup$

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.

$\endgroup$
4
  • $\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$
    – Glubbdrubb
    Mar 5 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$
    – Glubbdrubb
    Mar 9 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$
    – Jesse
    Mar 9 at 16:31
  • $\begingroup$ Thanks Jesse. That would be helpful. $\endgroup$
    – Glubbdrubb
    Mar 10 at 8:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.