Calculating an N50 from Velvet output
[This article was first published on Jermdemo Raised to the Law, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
In sequencing circles the N50 length is a useful heuristic for judging the quality of an assembly. Here is my definition of N50 length, which you may or may not find intuitive:Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
N50 length is the length of the shortest contig such that the sum of contigs of equal length or longer is at least 50% of the total length of all contigsFor example’s sake imagine an assembler has created contigs of the following length (in descending order):
91 77 70 69 62 56 45 29 16 4The sum of these is 519bp, so the sum of all contigs equal to or greater than N50 must be equal to or greater than 519/2 or 259.5
We can see by brute force that
91+77=168
91+77+70=238
91+77+70+69=307 (that’ll do)
so the N50 for this assembly is 69bp
Another way to look at this:
at least half the nucleotides in this assembly belong to contigs of size 69bp or longer.
N50 vs N50 length
Technically N50, as opposed to N50 length, refers to the ordinal of that last contig that pushes it over the brink – in this example 4 (since 69bp is the 4th largest contig). Unfortunately, a higher N50 implies the opposite of a longer N50 length. Some papers refer to N50 length as L50, while most have simply followed the lazy convention of dropping “length” off of “N50 length”. I think it is important to include units with your N50 to minimize confusion.Contig N50 vs Scaffold N50
Another distinction is often made between contig N50 and scaffold N50. Contigs are “contiguous segments”, while scaffolds (aka supercontigs) consist of contigs separated by gaps. Scaffolds are constructed using paired-end information at the read level and, in major sequencing projects, paired BAC ends. Because the scaffolds sequences are filled with varying quantities of empty N’s, the scaffold N50 should not solely be used as a comparative score of assembly quality.Velvet, when used with sane expCov settings, is very conservative with regard to scaffolding – so much that the contigs.fa N50 can be virtually considered a contig N50, as opposed to a scaffold N50. More aggressive programs, such as SOAPdenovo, produce separate contig and scaffold files.
Velvet N50 from stats.txt
Velvet is a popular assembler for short sequences that uses DeBruijn graphs and Eulerian graph theory instead of a repetitive align-consensus-align approach. Although it returns an N50 in the course of assembling, I wanted to derive it from the contigs themselves. These contigs are summarized in a table called stats.txtUsing R and its cumulative summation function we can easily compute N50.
Strategy: Order the contigs by decreasing size and find the first value for which the cumulative summation is at least half the total sum.
Thanks to Barry Rowlingson for providing this solution
myStatsTable<-read.table("stats.txt",header=TRUE) contigs<-rev(sort(myStatsTable$lgth+myKmerValue-1)) n50<-contigs[cumsum(contigs) >= sum(contigs)/2][1]
Beware the kmer
Note: The length in the stats.txt file is given as length=lgth+kmer-1, where kmer is the kmer value chosen for that assembly. The N50 length given in the Log file also appears to be in kmers. You cannot convert an N50 in kmers to bp by adding kmer-1. The math doesn’t work like that – you need to convert each contig to bp before recalculating N50.Finally, you can calculate N50 from sequences in the contigs.fa, but this file only contains contigs longer than 2-kmers by default. The contigs.fa bp-N50 will sometimes approximate the kmer-N50 in the Log file, but that is not a rule to depend on.
To leave a comment for the author, please follow the link and comment on their blog: Jermdemo Raised to the Law.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.