This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Forum: Careful!! Velvet Generates Different Assemblies From The Same Input Data And Same Input Parameters With Openmp Enabled.

Hi all,

I am using four different libraries of read size 76bps, insert sizes are 300bps, 1kb, 8kb and 12kb. Expected genome size is 80MB.

I am running velvet using these four libraries. Actually I tried velvet first for different k-mers and then for the best k-mer I tried different cov_cutoffs. In all these assemblies from same k-mer (69) with different cov_cutoffs, I used the same Roadmap and Sequences files from the initial velveth run (K-mer 69 and cov_cutoff default).

Surprisingly I got 10MB of N50 and 23MB of largest scaffold size using cov_cutoff of 12 (median coverage is 30.76 in Log file) on the previously generated Roadmap and Sequences files. Then later I tried all new assembly from the same reads, K-mer 69 and cov_cutoff 12, now my N50 is 2MB and largest scaffold size is 6.78 MB.

Later I tried the same input files and same velvet parameters and I figured out that velveth is generating different Roadmap file for the same k-mer 69, for all three runs. What could be the reason behind this? In this case it is not possible to regenerate the results.

I would really appreciate your comments on this.

Best regards,
Rahul

velvet genome
• 11,206 views • link
updated by Ram Ram • written by Rahul Sharma Rahul Sharma

4 answers

The solution is to set the OMP thread environment in your shell script:

#!/bin/bash

export OMP_THREAD_LIMIT=1
export OMP_NUM_THREADS=1

velveth ....
velvetg .

This will ensure the same results, but it will mean that velveth will likely take longer to execute. For finding the optimal kmer and cutoff, I recommend using VelvetOptimser. It is especially important to use these OMP settings if you are using VelvetOptimiser, as every single thread will try to use all the processors on that node if you don't have these variables set appropriately.

• 145 views • link
updated by Ram Ram • written by SES SES

Thanks for this, but then velvet takes ages to generate assemblies on a single processor. Regarding VelvetOptimiser, I cannot use it in my analysis. Actually I have huge dataset and it consumes ~95% of RAM of our whole group's machine. I need to submit jobs accordingly, in this case manual way of optimizing velvet is the right choice.

• 0 views • link
written by Rahul Sharma Rahul Sharma

Hi SES,

I am currently trying to setup my submission script for VelvetOptimizer on our cluster and was a bit confused as to how the various parameters had to be setup. The plan was to spread the instances over 24 threads, so I initially had the following parameters:

  • In the SLURM submission script: OMP_NUM_THREADS=24, OPENBLAS_NUM_THREADS=24, --cpus=24, --cpus-per-task=24, mem=256Gb
  • In the VelvetOptimizer command line: -t 24

However it appears that this generated thread allocation issues (it maybe tried to allocate 24*24 threads...?), and the program eventually crashed with the following error messages (not sure if the thread allocation actually caused the crash though; if you don't think so, I would appreciate any recommendation :-) ):

Use of uninitialized value in numeric ne (!=) at /home/umons/orgmb/csheridan/software/BioPerl/BioPerl-1.6.1/VelvetOptimiser.pl line 289.
Use of uninitialized value $maxScore in numeric gt (>) at /home/umons/orgmb/csheridan/software/BioPerl/BioPerl-1.6.1/VelvetOptimiser.pl line 290.
Use of uninitialized value in numeric gt (>) at /home/umons/orgmb/csheridan/software/BioPerl/BioPerl-1.6.1/VelvetOptimiser.pl line 290.
Aug 11 16:12:51 Hash value of best assembly by assembly score: 79
Aug 11 16:12:51 Optimisation routine chosen for best assembly: shortPaired
Aug 11 16:12:51 Looking for the expected coverage
Unable to open /gpfsuser/home/users/c/s/csheridan/data/auto_data_79/stats.txt for exp_cov determination.
 at /home/umons/orgmb/csheridan/software/BioPerl/BioPerl-1.6.1/VelvetOptimiser.pl line 838

Anyway, I then tried to correct this by setting the threads parameters as follows, but things appeared much slower with Velvet calculating only 3 hash values at a time:

  • In the SLURM submission script: OMP_NUM_THREADS=8, OPENBLAS_NUM_THREADS=8, --cpus=24, --cpus-per-task=24, mem=256Gb
  • In the VelvetOptimizer command line: -t 3

So my question is, how do you exactly set the threads parameters in order to gain from the use of OPENMP, and how can Velvet/Oases be optimally parallelised?

Thank you!

• 2 views • link
updated by Ram Ram • written by sheridan.christopher • sheridan.christopher

It is not a good idea to post a new question in the comment section. On Biostar we like to keep a single topic in thread. In addition very few people see your post so it is also inefficient.

• 2 views • link
updated by Ram Ram • written by Istvan Albert Istvan Albert

I ran a test and my files are identical:

velveth a31 31 -shortPaired -fastq -separate r1.fq r2.fq
velveth b31 31 -shortPaired -fastq -separate r1.fq r2.fq
cmp a31/Roadmaps b31/Roadmaps

No output returned from cmp. Running the velvetg on each also produces identical results:

velvetg a31 -exp_cov auto
Final graph has 1055 nodes and n50 of 15415, max 65035, total 1446385, using 159599/200002 reads

velvetg b31 -exp_cov auto
Final graph has 1055 nodes and n50 of 15415, max 65035, total 1446385, using 159599/200002 reads
• 171 views • link
updated by Ram Ram • written by Istvan Albert Istvan Albert

Thanks for your answer, but please try 3-4 different libraries. The result will be same with the single library. Here is what I did with my test dataset from four libraries of different insert:

nohup velveth test1 69 -shortPaired -fastq demo_300b.fastq -shortPaired2 -fastq demo_1kb.fastq -shortPaired3 -fastq demo_8kb.fastq -shortPaired4 -fastq demo_20kb.fastq &

nohup velveth test2 69 -shortPaired -fastq demo_300b.fastq -shortPaired2 -fastq demo_1kb.fastq -shortPaired3 -fastq demo_8kb.fastq -shortPaired4 -fastq demo_20kb.fastq &

$ wc -l ./test1/Roadmaps #Lines in Roadmap from test1
4416 Roadmaps

$ wc -l ./test2/Roadmaps #Lines in Roadmap from test2
4436 Roadmaps

Assemblies are also different, but I am wondering how? Daniel also replied me that velvet generates different assemblies from same input reads and input parameters. Then how can one regenerate the results?

Best regards,
Rahul

• 0 views • link
updated by Ram Ram • written by Rahul Sharma Rahul Sharma

I don't fully understand why results would be same with a single library and different for multiple libraries. One could just put all data into one library, in your case all data is short paired, no?

That being said I have never had the situation where two consecutive assemblies were radically different. As you point out it makes no sense because one could not regenerate results.

That is the reason I gave it a shot, I would not be surprised in seeing small differences but radical changes would be very strange.

I am not saying your problem is not valid, but if true it would have big repercussions, and it may have to do perhaps with the order of combining reads with very large size differences, 300bp to 20Kb.

• 0 views • link
updated by Ram Ram • written by Istvan Albert Istvan Albert

For what it's worth, this link to a Velvet-users thread was mentioned on the equivalent post on SeqAnswers. So things should be deterministic in simple cases (of course, whether they actually are is a different question).

• 145 views • link
written by Devon Ryan Devon Ryan

Thanks for the link, I checked it and found that velvet gives different results with OPENMP and reordering the reads. In my case I have compiled velvet using OPENMP=1 and using multithreading. Thats why it is giving different Roadmaps and assemblies. Yes my data is paired-end, as I have explained in my question.

• 0 views • link
written by Rahul Sharma Rahul Sharma

what is shortPaired3 and shortPaired4 there my velvet build does not have those, I can just list the -shortPaired again.

plus shouldn't the 20Kb reads be assembled as long?

• 0 views • link
written by Istvan Albert Istvan Albert

The number after the read type specifies the different read sets, and the max number used cannot be larger than whatever the CATEGORIES parameter was set to during compilation. The 20kb is likely an insert size (short reads, large insert). I don't think it's even possible to generate reads of that length except maybe from an assembly.

• 151 views • link
written by SES SES

Because you have installed velvet with default Makefile setting. I have five different libraries, thats why I installed velvet with CATEGORIES =5. Regarding long, it refers to the read size, not the insert size. All of the libraries are Illumina, has read length <100bps so I will not use long option. For 454 or longer reads >400bps one should use long option.

• 0 views • link
written by Rahul Sharma Rahul Sharma

Yes, the long option is for use with Sanger and 454 data. Five seems like a very low number to set the CATEGORIES to, you may want to set that to something arbitrarily high so you don't have to debug, then recompile in the future. For example, my CATEGORIES is set to 62.

• 0 views • link
written by SES SES

Yes thats true, I have set it to 10. But I just mentioned 5 for that example dataset. But setting MAXKMERLENGTH=250 consumes more RAM than MAXKMERLENGTH=100, if someone wants to use k-mer upto k-mer=70. So one should take care of these things while compilation.

• 0 views • link
written by Rahul Sharma Rahul Sharma 
Hello
i run velvet for transcriptome paired data by using cov_cutoff 4.0 min_contig_lgth 200
Final graph has 251014 nodes and n50 of 562, max 8255, total 32301310, using 0/32799252 reads

what are the nodes 251014 and maxand and toatal numbers indicates in velvet output. i could not able get it.

Please answer me on this.

• 1 views • link
written by reddy.renuka5 • reddy.renuka5

Now it has been clear that velvet generates different assemblies because of OPENMP, now what should be considered. Should I go for the 10MB N50 size assemblies? How much this will be realistic? For the publication point of view, our results should be re-generated using same methods. But when I run velvet 5 different times, with same data, k-mer 69 and cov_cutoff 12 it is giving me five different assemblies. Which assembly one should trust then?

I would appreciate your suggestions on this!

Thanks in advance! Rahul

• 147 views • link
written by Rahul Sharma Rahul Sharma

Well as long as we accept that velvet works and each assembly is valid just that one is more complete than the other I would pick the best assembly that works.

After all this is what people do anyhow when sweeping over the parameter space.

Then describe in the supplementary materials the issue with reproducibility.

• 238 views • link
written by Istvan Albert Istvan Albert

Thanks for your answer, but please try 3-4 different libraries. The result will be same with the single library. Here is what I did with my test dataset from four libraries of different insert:

nohup velveth test1 69 -shortPaired -fastq demo_300b.fastq -shortPaired2 -fastq demo_1kb.fastq -shortPaired3 -fastq demo_8kb.fastq -shortPaired4 -fastq demo_20kb.fastq &

nohup velveth test2 69 -shortPaired -fastq demo_300b.fastq -shortPaired2 -fastq demo_1kb.fastq -shortPaired3 -fastq demo_8kb.fastq -shortPaired4 -fastq demo_20kb.fastq &


$ wc -l ./test1/Roadmaps #Lines in Roadmap from test1
4416 Roadmaps


$ wc -l ./test2/Roadmaps #Lines in Roadmap from test2
4436 Roadmaps

Assemblies are also different, but I am wondering how? Daniel also replied me that velvet generates different assemblies from same input reads and input parameters. Then how can one regenerate the results?

Best regards,

Rahul

• 0 views • link
written by Rahul Sharma Rahul Sharma

Please move this to a comment on Istvan's answer.

• 0 views • link
written by Devon Ryan Devon Ryan

Log in to answer this question.

More posts like this