Assemble short reads based on k-mers

Question

[venu]

Hello all,

I am completely new for this kind of tasks. I have data like this,

>in0
GATCCTCGAAGTTACACGGG
>in1
TACGTCGACGTCAATCCGGG
>in2
TACACGGGCCGCTCCTGGGC
>in3
ACGGGGTACTACGAGACGCG
>in4
AGGGGGAATGTGGTCCACAT
>in5
TCCACATGGCTTGCTCCTGA
>in6
CTTGACGTTATGAATTTCGC

and so on..I need to assemble these short reads. I want to use perl for this. I just need a pseudo code on how to do this or direct me to a good resource. At the end I need a single string containing consensus sequence.

[Devon Ryan]

Is their a reason you want to reinvent the assembly wheel? There are a good number of assemblers already written, why bother writing yet another one without a good reason?

[venu]

Can you give me some examples, so that I can find them directly on the internet.

[dylan.storey]

You could try google

[Devon Ryan]

As orange said, SOAPdenovo is one option. Others would include Trinity and Minia. There are quite a few of these if you just search pubmed for "DNA assembler" or "DNA assemble".

Answer 1

[orange]

why not try SOAPdenovo instead of perl ? Did you not aseembly genome sequence before ?

Answer 2

[thackl]

Assuming you want to stick to Perl for educational purpose:

Here is some code to quite efficiently create kmers with Perl: https://github.com/thackl/perl5lib-Kmer/blob/master/lib/Kmer.pm.

Perl is not really made for handling graph structures, but there is one module that you could use to set up a De-Bruijn structure: http://search.cpan.org/~jhi/Graph-0.96/lib/Graph.pod. I played around with it some time ago but did not follow through.

[venu]

I just need a single string of consensus sequence from the above shown file

[thackl]

But definitely in Perl?

[venu]

Not exactly, but any simple program that receives the above file as input and outputs the consensus string. I can understand the perl code easily, so the tag.

[thackl]

But your data sets are small - you want to do some form of microassembly?

[venu]

Yes. I've 200 such reads in a file and I want an output like

 GGCATTTAACCGAAGCCGGTGGGTTAGACTATGATCCTCGAAGTTACACGGGCCGCTCCTGGGCGTGGCTGCTCCCAGCCCTAGCCCCAATGTAATATAAAGGTCGTGCCCAGTTAGCGTTAAGCAAGAGGTGTTACAAATATCTTGGAGAGTCATGTCGCAATTCTTGACGTTATGAATTTCGCGGTGAACAATGTCGCCCAGAATGGCAGGTCATGAAAAGCTTCAGCGGGAACCAGCAC....

[thackl]

What is the coverage of the reads and the length?

[venu]

I am doing this kind of work for the first time. What I know is each read has different lengths of k-mers.

[nterhoeven]

Take a look at this: http://www.homolog.us/Tutorials/index.php?p=1.1&s=2

This should give you a basic overview about assemblies

Answer 3

[thackl]

You can use SPAdes:

spades.py --only-assembler --sc -k 33 -s in.fq -out asm
  # results are in asm/contigs.fa

This should also work for lowish coverage (5-10X) but assumes little to no errors in your data. Also, there can be multiple contigs and regions with very low coverage (just a 1-3 reads at a certain position, e.g. the ends of you target sequence) will be missing.