8
$\begingroup$

Working on a fastQ file, I would like to get the occurrences of repeated sequences for all possible kmers of a given length that cover at least 90% of the read's length for the whole data set.

example :

for a length 6 with the kmer "ATTGGG" and a data set containing reads of 300 bases each

I want to obtain the occurrences of the following kmers :

ATTGGGATTGGGATTGGGATTGGGATTGGGATTGGG......  (90 % of the read's length without  gaps) 
ATTGGGATTGGGATTGGGATTGGGATTGGGATTGGGA.....  (91%)
ATTGGGATTGGGATTGGGATTGGGATTGGGATTGGGAT....  (92%)

this for all possible combinations of kmers of length = 6.

I have done some research and tested many softwares such as (KMC2,Jellyfish,Khmer,Dsk,Scturtle) but none matched what I would like to do, because I want to obtain the occurrences of the repeated kmer without any possible gaps in-between with a certain % of read coverage.

I've done a naive algorithm that works but the optimization isn't the best obviously. Is there a tool that does the trick ?

$\endgroup$
0

1 Answer 1

6
$\begingroup$

The following javascript does what you want. You need node.js to run it. It should be easy to translate the code to Python. I have not carefully tested it. Use with caution.

EDIT (response to new comments): the program has been changed to compute the sum of lengths. Note that only a tandem repeat of length k*2 or longer is counted. For example in sequence ATTGGGATTGGGATTcGGGATTGGG, the script returns 15, because the second is not long enough. Counting length smaller than k*2 will be slower and more complex, but is often not the right thing to do. For example, in ATTGGGATTGGGcccccccgaaatcgatagcatcgaGGGATTGcgatc, we would not count the second GGGATTG as a tandem repeat.

On performance, for an input string of length $l$, the time complexity of this script is $O(l)$ associated with a small constant. It will be efficient and adequate to practical uses.

If you use TRF, you should run it on your sequences and then test if any tandem repeats it finds match GGGAAT – the procedure will be similar to my script, except that TRF will also find impure repeats.

// given the repeat unit length _k_, find all tandem repeats of 2*k or longer
function trf_k(k, str)
{
    var streak = 0, a = [];
    for (var i = k; i <= str.length; ++i) {
        if (i < str.length && str[i] == str[i - k]) {
            ++streak;
        } else {
            if (streak >= k)
                a.push([i - streak - k, i]);
            streak = 0;
        }
    }
    return a;
}
// find the sum of tandem repeat length with _kmer_ being the repeat unit
function trf_kmer(kmer, str)
{
    var a = trf_k(kmer.length, str), sum = 0;
    for (var i = 0; i < a.length; ++i)
        if (str.substring(a[i][0], a[i][1]).indexOf(kmer) >= 0)
            sum += a[i][1] - a[i][0];
    return sum;
}

var str = 'xATTGGGATTGGGATTGxGATTGGGATTGGGATTGGGATTGxx';
console.log(trf_kmer('GGGATT', str));
$\endgroup$
3
  • $\begingroup$ Sorry i mean that i want the perfect match of each possible repeated sequences, and in the case of your algorithm, when there is a gap it's still count the kmer. Indeed the task is highly customized. I guess I'll just start using parallelism on my program. I'll take a look into TRF, thanks. $\endgroup$
    – hilta007
    Mar 15, 2018 at 10:00
  • 1
    $\begingroup$ @hilta007 Thanks for the explanation (I am deleting my old comment BTW). If you want to find the total length, just change "max" to "sum", but with a caveat (see edits for detail). My implementation should be fairly efficient. I wouldn't worry too much about performance. A good single-threaded script can sometimes be faster than a bad multi-threaded C/C++ program. $\endgroup$
    – user172818
    Mar 15, 2018 at 15:13
  • 1
    $\begingroup$ Thanks for the edits, waiting the results of TRF and my own algorithm to compare quality/speed for now. I'll keep your algorithm in mind. $\endgroup$
    – hilta007
    Mar 15, 2018 at 16:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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