[Brian Bushnell]

Yes, there is a reference that is pretty close, <2% difference at the base level, with a few rearrangements likely. I use Btrim for the trimming.

Sounds close enough to use Quast to get a very good idea of assembly quality. You can run it with the assembly and reference as the inputs and it will tell you the quality of the assembly with respect to the reference. Even with 2% difference and some rearrangements, the assembly closest to the reference should be the highest quality assembly.

[Genomics101]

Yes, there is a reference that is pretty close, <2% difference at the base level, with a few rearrangements likely. I use Btrim for the trimming.

[Brian Bushnell]

Well... if there is a reference, I would compare them to the reference. With or without a reference, I think Quast is fairly useful. It tells you the number of predicted genes and so forth in your assembly; generally, the more longer genes you have and the fewer shorter genes, the higher the assembly quality. Most useful if you have more than one assembly. If you have a reference, it will also tell you the number of misassemblies.

Similarly, if you map all of the source reads to the assembly, then the number of non-match symbols (insertions, deletions, substitutions) correlates with assembly quality, as does % of reads mapped, % paired, and % unambiguous mappings. We use BBMap for that purpose here, though again, it's better for determining the relative quality of two assemblies than the absolute quality of one assembly. You can also try ALE, which attempts to judge the probability of an assembly being correct given the reads. This is based on mismatches, deviation of paired insert size from normal, and coverage depth.

If you have EST data, then calculating the EST capture rate is also very useful in determining genome completeness and misassemblies.

Incidentally, I think Q30 is much too high for trimming. Also, what trimming tool are you using? In my tests, the top performing ones use the phred algorithm (this includes seqtk and my own); all other algorithms were consistently inferior. The attached graph shows error rate of reads versus bases remaining from the initial 150Mbp dataset after quality-trimming; each point represents a different quality cutoff (so the leftmost point is trimming to Q40, then Q39, etc). The black line is the best, but that uses error-correction AND quality-trimming; all the others just use quality-trimming. seqtk and bbtrim appear to be identical. The error rates were calculated by mapping back to the Pedobacter reference.

[Genomics101]

Trimming to ~180bp seems like overkill for data of that quality; I would probably try trimming to something very conservative like Q10. .

Thanks. I didn't trim by length, but by quality (Q30 cut off) and the reads just came out with most of them at around. But doing a less strict trimming may be the answer.

Since I have your very helpful attention here, since I am doing the assmebly with the untrimmed reads, do you have a suggestion for a good way to assess how the sequencing errors are affecting the accuracy of the contigs? Should I just BLAST a few regions I have done with Sanger sequencing?

[Brian Bushnell]

Read length uniformity shouldn't matter to Velvet. Trimming to ~180bp seems like overkill for data of that quality; I would probably try trimming to something very conservative like Q10. Excessive trimming can also cause biases.

[Genomics101]

I also have the option of using the longer and more uniform untrimmed reads, but the quality is pretty questionable:



I tried doing an assembly with these and got better N50s at very high kmers (~99-135) but the kmer depth has a weird bell curve relationship with kmer size rather than a direct one. The lowest kmer I tried (45) had a coverage of only 1,2 and it the kmer coverage peaked at 41.5X at kmer =93 (also a relative good N50 at ~20kb). Also, I am very wary of using data with so many errors.

[Brian Bushnell]

Ah - what I actually mean is, a graph for a fixed kmer length (of, say, 31) where the X axis is depth and Y axis is number of kmers found at that depth, both log-scale. Ideally, you should have a sharp peak at some depth (maybe 40) and it should drop dramatically on the left and right.

My attachment shows the 31-mer frequency histogram for e.coli synthetic reads. You can see a main peak at about 200 and a few repeat peaks after that. If the data was real and had uneven coverage, there would be a broad peak rather than a sharp one.

FYI, I generated this with the 'khist.sh' script in the BBMap package and plotted it in Excel.

[Genomics101]

Hmm, well in that case, I am surprised that the N50 is best at such a low K. Unless the coverage is highly non-uniform, as can happen if data is over-amplified. Do you have a kmer-depth histogram?

[AdrianP]

Hmm, well in that case, I am surprised that the N50 is best at such a low K. Unless the coverage is highly non-uniform, as can happen if data is over-amplified. Do you have a kmer-depth histogram?

To obtain that, I can recommend:

http://kmergenie.bx.psu.edu/

CMD:
./kmergenie --diploid -k <higher_kmer> -e 1 -l <lower_kmer> -t <cpu_threads> -o <output_name> <read_location>

Start with higher 101, and lower 41, see what graphs.

[Brian Bushnell]

@Brian Bushnell The insert size mean valus is 407 with an SD of 130, and the 23X is kmer coverage, the base coverage is about 43X. The per sequence quality is between 36 and 38 for almost all of them. There are no adaptors as they are removed by the Illumina 1.9 pipeline.

Hmm, well in that case, I am surprised that the N50 is best at such a low K. Unless the coverage is highly non-uniform, as can happen if data is over-amplified. Do you have a kmer-depth histogram?

[AdrianP]

@Brian Bushnell The insert size mean valus is 407 with an SD of 130, and the 23X is kmer coverage, the base coverage is about 43X. The per sequence quality is between 36 and 38 for almost all of them. There are no adaptors as they are removed by the Illumina 1.9 pipeline.

Okay, with a base coverage of 43X you do not want Velvet. DBG needs high coverage for repeat resolution.

My advice, is to use SeqPrep to merge your reads. You should have 3 files, forward, reverse, and merged after using it. Feed those to the MIRA assembler, which is an OLC assembler, and I expect you to get better results.

[Genomics101]

@AdrianP Thanks for your reply, but I actually did kmers (with a larger gap between them ) all the way up to 191 as the initial analysis:

[Genomics101]

@Brian Bushnell The insert size mean valus is 407 with an SD of 130, and the 23X is kmer coverage, the base coverage is about 43X. The per sequence quality is between 36 and 38 for almost all of them. There are no adaptors as they are removed by the Illumina 1.9 pipeline.

[AdrianP]

Kmer coverage

Okay, than what the person in post #2 said doesn't apply, because they assumed nucleotide coverage. You need to go higher kmers. kmer coverage of higher than 20 is a waste, you need to aim between 10 and 20.

Use VelvetOptimiser, and try kmers to 160-180, you might see a second peak in N50, this is common.

[Genomics101]

Kmer coverage