[anchor] Add getL50()

anchor.go

@@ -91,13 +91,12 @@ type PathSet map[*Path]bool
 // Run kicks off the merging algorithm
 func (r *Anchorer) Run() {
 	// Prepare the paths to run
-	flanksize := int64(1000000)
 	nIterations := 1
 	r.ExtractInterContigLinks()
+	flanksize := int64(1000000)
 	paths := r.makeTrivialPaths(r.contigs, flanksize)
 	for i := 0; i < nIterations; i++ {
 		r.iterativeGraphMerge(paths, flanksize)
-
 		// Attempt to split bad joins
 		r.makeContigStarts()
 		piler := r.inspectGaps(LinkDist)
@@ -118,9 +117,11 @@ func (r *Anchorer) iterativeGraphMerge(paths PathSet, flanksize int64) {
 	prevPaths := len(paths)
 	graphRemake := true
 	for prevPaths > 1 {
+		// flanksize = getL50(paths) / 4
 		if graphRemake {
 			i++
-			log.Noticef("Starting iteration %d with %d paths", i, len(paths))
+			log.Noticef("Starting iteration %d with %d paths (L50=%d)",
+				i, len(paths), getL50(paths))
 			G = r.makeGraph(paths)
 		}
 		CG := r.makeConfidenceGraph(G)
@@ -596,3 +597,13 @@ func (r *Anchorer) serialize(res int64, jsonfile, npyfile string) {
 	w.WriteInt32(C)
 	log.Noticef("Matrix (resolution=%d) written to `%s`", res, npyfile)
 }
+
+// getL50 computes the L50 of all component contigs within a path
+func getL50(paths PathSet) int64 {
+	pathLengths := []int64{}
+	for path := range paths {
+		pathLengths = append(pathLengths, path.length)
+	}
+
+	return L50(pathLengths)
+}

base.go

@@ -239,6 +239,26 @@ func median(s []float64) float64 {
 	return result
 }
 
+// L50 returns the sequence length L where half of the genome is covered in contigs
+// of length >= L50
+func L50(lengths []int64) int64 {
+	sortInt64s(lengths)
+	total := int64(0)
+	for _, length := range lengths {
+		total += length
+	}
+	halfTotal := total / 2
+	cumsize := int64(0)
+	i := len(lengths) - 1
+	for ; i >= 0; i-- {
+		cumsize += lengths[i]
+		if cumsize > halfTotal {
+			break
+		}
+	}
+	return lengths[i]
+}
+
 // OutlierCutoff implements Iglewicz and Hoaglin's robust, returns the cutoff values -
 // lower bound and upper bound.
 func OutlierCutoff(a []float64) (float64, float64) {

graph.go

@@ -191,7 +191,7 @@ func (r *Anchorer) getUniquePaths() PathSet {
 		}
 	}
 
-	log.Noticef("%d paths (nComplex=%d nSingleton=%d), %d contigs (nComplex=%d nSingleton=%d)",
+	log.Noticef("%d paths (nComplex=%d nSingle=%d), %d contigs (nComplex=%d nSingle=%d)",
 		nComplex+nSingleton, nComplex, nSingleton,
 		nComplexContigs+nSingletonContigs, nComplexContigs, nSingletonContigs)