SEQanswers

Go Back   SEQanswers > Applications Forums > De novo discovery



Similar Threads
Thread Thread Starter Forum Replies Last Post
inchworm and paired ends ians Bioinformatics 4 10-05-2011 06:55 AM
bfast bgeneratereads for paired ends sdvie Bioinformatics 2 03-23-2011 10:04 AM
paired ends with cuffdiff Greg Bioinformatics 1 07-05-2010 10:12 AM
MAQ paired ends prm36 Bioinformatics 0 04-15-2010 07:29 AM
consed and illumina paired ends ggloor Bioinformatics 1 06-22-2009 08:46 PM

Reply
 
Thread Tools
Old 04-23-2011, 10:27 AM   #1
AdrianP
Senior Member
 
Location: Ottawa

Join Date: Apr 2011
Posts: 130
Default Velvet & paired-ends

Greetings to you all,

First of all thank you for creating this forum, it seems like a great way to share knowledge.

I am a undergraduate student starting on a job at my university involving assembly of sequenced data. I am an expirienced linux user so I will try to use linux programs.

Currently I have familiarized myself with velvet a bit and I wish to try to assemble some data that has been previously well assembled and I wanted to test velvets capabilities on it as well. After making sure I know how to use velvet I will also try Ray and SOAP.

So I have 2 files:
100611_s_4_1_seq_GDR-7.txt (1.6 GB)
100611_s_4_2_seq_GDR-7.txt (1.6 GB)

I have used this as a refrence for my work:
http://wiki.bioinformatics.ucdavis.e...y_using_velvet

From what I understand I need to merge the two files into 1 file with shuffleSequences_fastq.pl. Is this correct?
Code:
shuffleSequences_fastq.pl 100611_s_4_1_seq_GDR-7.txt100611_s_4_2_seq_GDR-7.txt 100611_s_4_both_seq_GDR-7.txt
I have trouble understanding what does subsetting mean from that page. If we look at:
Code:
17. Do the subsetting. Soon we will compare the single ended assembly to the paired-end assembly. In order for the comparison to be fair, we must use the same total number of reads. Therefore each paired end file will contain 1/4 of the reads:
It seems like they are trying to compare single ended with paired-end. I am only doing paired-end, do I need to do subsetting?
Code:
Final graph has 302039 nodes and n50 of 175, max 1779, total 5984104, using 13024530/17609332 reads
Now one more thing, after running velveth_de and velvetg_de, at the end of velvetg_de it tells me how many nodes have been created. Is that the number of contigs? How do I interpret that last line?

I am using the -shortPaired option for velveth.

My last question is, has anyone here used consed? I just wanted to ask since I have some problems setting up that program.

Thank you.

Last edited by AdrianP; 04-23-2011 at 10:42 AM.
AdrianP is offline   Reply With Quote
Old 04-26-2011, 01:17 AM   #2
Thorondor
Member
 
Location: Heidelberg

Join Date: Feb 2011
Posts: 69
Default

if you don't want to compare paired end to single end you don't need to do "subsetting".

the last line just tells you your N50, how many reads were used (you can also request velvet to output a UnusedReads.fa file) etc, i am not sure if the number of nodes is the number of your contigs. You can check that by grep ">" -c contigs.fa.
Thorondor is offline   Reply With Quote
Old 04-26-2011, 06:35 AM   #3
AdrianP
Senior Member
 
Location: Ottawa

Join Date: Apr 2011
Posts: 130
Default

Thank you for your previous reply.
Is there any way to see which contig is largest? Somehow sort contigs by their size?
AdrianP is offline   Reply With Quote
Old 04-26-2011, 07:25 AM   #4
Thorondor
Member
 
Location: Heidelberg

Join Date: Feb 2011
Posts: 69
Default

sure there are a lot of ways. ;-) Use a script, you even have the (length + kmer -1) in the id of the contig so it is really easy.

here are some perl scripts that might help:
http://wiki.bioinformatics.ucdavis.e.../Data_Analysis
Thorondor is offline   Reply With Quote
Old 04-27-2011, 06:34 PM   #5
AdrianP
Senior Member
 
Location: Ottawa

Join Date: Apr 2011
Posts: 130
Default

If Velvet assembles poor (small) contigs when other programs with same settings (coverage and insertsize) do much much better, what can be my conclusions?

By the way most of those scripts are for fasta, and i got fastq, is there a script that converts?

Adrian
AdrianP is offline   Reply With Quote
Old 04-27-2011, 11:17 PM   #6
Thorondor
Member
 
Location: Heidelberg

Join Date: Feb 2011
Posts: 69
Default

http://brianknaus.com/software/srtoolbox/fastq2fasta.pl

first hit in google. ;-) Also normally trimming "programs" takes fastq as input and output a fasta.

velvet needs a good coverage to do well because it's de brujin based. Since I don't know on what data you run velvet I and what you expect there is no help. Try smaller kmers, try different parameters, try multiple kmer....
Thorondor is offline   Reply With Quote
Old 04-28-2011, 02:25 AM   #7
AdrianP
Senior Member
 
Location: Ottawa

Join Date: Apr 2011
Posts: 130
Default

One thing that is puzzeling me is that Genegenious takes a few days to assemble the data that velvet assembles in 10-15 mins. Is it normal that velvet runs so fast?
AdrianP is offline   Reply With Quote
Old 04-28-2011, 02:58 AM   #8
Thorondor
Member
 
Location: Heidelberg

Join Date: Feb 2011
Posts: 69
Default

i have no clue what genegenious is and on what algorithm it is based. So if it is a ovelap-based method, yes it is possible and depends on kmer, amount of reads you have, expected coverage, read length .....
Thorondor is offline   Reply With Quote
Old 05-03-2011, 10:14 AM   #9
AdrianP
Senior Member
 
Location: Ottawa

Join Date: Apr 2011
Posts: 130
Default

I was wondering a bit more about Velvet's last line output.

What is n50? (I understand that it is a measurement of quality, the higher the better???)
What is max?
What is total?
What are nodes? (this isn't as important since i believe it is related to the graph that velvet builds, and I do not use the graph, just the final contigs.fa file. Should I use the graph?)

Also, what is the diffrence between shortpaired and shortpaired2 ? Something to do with inert libraries...

Thanks a lot!
AdrianP is offline   Reply With Quote
Old 05-04-2011, 12:47 AM   #10
Thorondor
Member
 
Location: Heidelberg

Join Date: Feb 2011
Posts: 69
Default

come on, do a bit more research on your own. :P

read the velvet paper:
http://genome.cshlp.org/content/18/5/821.short

You don't want to use the graph and if you know what a graph is you should also know what nodes are. Nodes are the vertices of a graph.
shortpaired2 ist the same as shortpaired but for a separate insert size library (also stated in the manual).
Thorondor is offline   Reply With Quote
Old 05-04-2011, 09:33 AM   #11
AdrianP
Senior Member
 
Location: Ottawa

Join Date: Apr 2011
Posts: 130
Default

The Manual I have read it a few times but I guess what I was asking is what does that mean "separate insert size library" ? Separate from what?

As for the research paper, i had a look at it before, I can't find a defenition for n50, it jumps straightly to using that term, eve in the abstract. I would not have asked if I did not do the research myself.

I found the answer here after googling n50 http://seqanswers.com/forums/showthread.php?t=2332
AdrianP is offline   Reply With Quote
Old 05-04-2011, 10:19 AM   #12
Thorondor
Member
 
Location: Heidelberg

Join Date: Feb 2011
Posts: 69
Default

yes, it is the first hit when you google "definition: N50"

separate to "shortpaired" I assume. So you can use 2 different PE libraries with different insert sizes, but maybe I am wrong.

Total is their calculated base pairs total, but keep in mind that the bp length of eacht transcript is the "real" bp length minus kmer plus 1, as mentioned somewhere.^^

max might be the longest contig, but I don't remember it exactly, since you need to do more statistics anyway. ;-)
Thorondor is offline   Reply With Quote
Old 05-05-2011, 01:48 AM   #13
tonybolger
Senior Member
 
Location: berlin

Join Date: Feb 2010
Posts: 156
Default

Quote:
Originally Posted by AdrianP View Post
If Velvet assembles poor (small) contigs when other programs with same settings (coverage and insertsize) do much much better, what can be my conclusions?
Could be many things: appropriate vs inappropriate settings, more/less sensitive to data error vs coverage, or simply wrong/right tool for this particular job.

You also could have one tool making a decent size but completely incorrect assembly, with another making a cautious but correct assembly. N50/size isn't everything. Can you validate your results somehow, e.g. against another closely-related known genome?

Either way, getting the best out of a dataset may require months of trial and error, tweaking etc. Even getting a particular tool do run properly and produce decent output can take weeks and can make a massive difference - the DBG assemblers all seem to have glass jaws. Pre-filtering the data seems to make a massive difference to most though.

Incidentally, I don't have a lot of experience with velvet in particular, since it's simply too heavy for my project (>1GBase genome)
tonybolger is offline   Reply With Quote
Old 05-05-2011, 03:31 AM   #14
AdrianP
Senior Member
 
Location: Ottawa

Join Date: Apr 2011
Posts: 130
Default

Yeah actually my next genome to work with is a mitGenome and is about 70k, pretty cool.

I will start working with consed, not an easy program to work with but as I understand incredibly useful.
AdrianP is offline   Reply With Quote
Old 05-31-2011, 10:16 AM   #15
seb567
Senior Member
 
Location: Québec, Canada

Join Date: Jul 2008
Posts: 260
Default

Quote:
Originally Posted by AdrianP View Post
Greetings to you all,

After making sure I know how to use velvet I will also try Ray and SOAP.

Hi !

I am the author of Ray so if you have any question, ask away.

Basically, with Ray, you will need to convert your two files to fasta or fastq format.

There is a script in maq for that.
http://maq.sourceforge.net/

maq-0.7.1/scripts/fq_all2std.pl export2std 100611_s_4_1_seq_GDR-7.txt > 100611_s_4_1_seq_GDR-7.txt.fastq
maq-0.7.1/scripts/fq_all2std.pl export2std 100611_s_4_2_seq_GDR-7.txt > 100611_s_4_2_seq_GDR-7.txt.fastq


Ray is available at http://denovoassembler.sf.net

Then, using Ray, you assemble these reads:

mpirun -np 8 Ray -k 31 -p 100611_s_4_1_seq_GDR-7.txt.fastq 100611_s_4_2_seq_GR-7.txt.fastq -o Ray-test-1.4.0

I encourage you to explore the files written by Ray:

ls Ray-test-1.4.0.*

see http://denovoassembler.sourceforge.n...00000000000000
seb567 is offline   Reply With Quote
Old 05-31-2011, 10:18 AM   #16
AdrianP
Senior Member
 
Location: Ottawa

Join Date: Apr 2011
Posts: 130
Default

Quick Question, does Ray support hybrid assembly? 454 and illumina? Can ray do a reference assembly? Is there GUI?

I almost got to using RAY but had huge problems with openmpi.

Adrian
AdrianP is offline   Reply With Quote
Old 05-31-2011, 10:56 AM   #17
seb567
Senior Member
 
Location: Québec, Canada

Join Date: Jul 2008
Posts: 260
Default

Quote:
Originally Posted by AdrianP View Post
Quick Question, does Ray support hybrid assembly?
Yes.

Ray: simultaneous assembly of reads from a mix of high-throughput sequencing technologies.
Sébastien Boisvert, François Laviolette, and Jacques Corbeil.
Journal of Computational Biology (Mary Ann Liebert, Inc. publishers).
November 2010, 17(11): 1519-1533.
doi:10.1089/cmb.2009.0238

See others here too: http://denovoassembler.sourceforge.n...lications.html




Quote:
Originally Posted by AdrianP View Post
454 and illumina?
Yes. See paper above.

Quote:
Originally Posted by AdrianP View Post
Can ray do a reference assembly?
No. Ray only does de novo stuff.

Quote:
Originally Posted by AdrianP View Post
Is there GUI?
Not in the vanilla, but I know that Applied Maths Inc. is working on one.

See http://www.applied-maths.com/opensource/opensource.htm

Quote:
Originally Posted by AdrianP View Post
I almost got to using RAY but had huge problems with openmpi.

Adrian
Open-MPI is very easy to install. If you have a decent GNU/Linux distribution, then you can install it right away as a package.

Otherwise:

wget http://www.open-mpi.org/software/omp...-1.4.3.tar.bz2
tar xjf openmpi-1.4.3.tar.bz2
cd openmpi-1.4.3
./configure --prefix=$(pwd)/build
make
make install

You now have Open-MPI in $(pwd)/build

In particular, $(pwd)/build/bin/mpic++ and $(pwd)/build/bin/mpirun.

Now add $(pwd)/build/bin to your path.

For example in your .bashrc:

export PATH=/home/boiseb01/software/openmpi-1.4.1/build/bin:$PATH


Anyway, you should use a decent GNU/Linux distribution -- it makes things easier.
seb567 is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 03:23 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2019, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO