variantAt := map[int][]hgvs.Variant{} // variantAt[chromOffset][genomeIndex*2+phase]
for refstep, libref := range reftiles {
reftile := tileVariant[libref]
- coverage := int64(0) // number of ref bases that are called in genomes -- max is len(reftile.Sequence)*len(cgs)*2
+ tagcoverage := 0 // number of times the start tag was found in genomes -- max is len(cgs)*2
for cgidx, cg := range cgs {
for phase := 0; phase < 2; phase++ {
if len(cg.Variants) <= int(libref.Tag)*2+phase {
if variant == 0 {
continue
}
+ tagcoverage++
genometile := tileVariant[tileLibRef{Tag: libref.Tag, Variant: variant}]
if variant == libref.Variant {
continue
}
varslice[cgidx*2+phase] = v
}
- coverage += int64(len(reftile.Sequence))
}
}
refpos += len(reftile.Sequence) - taglen
}
thickend := refpos
// coverage score, 0 to 1000
- score := 1000 * coverage / int64(len(reftile.Sequence)) / int64(len(cgs)) / 2
+ score := 1000 * tagcoverage / len(cgs) / 2
fmt.Fprintf(bedw, "%s %d %d %d %d . %d %d\n",
seqname, tilestart, tileend,
libref.Tag,
c.Check(err, check.IsNil)
c.Logf("%s", string(bedout))
c.Check(string(bedout), check.Equals, `chr1 0 248 0 500 . 0 224
-chr1 224 372 1 250 . 248 348
-chr1 348 496 2 0 . 372 472
-chr1 472 572 3 0 . 496 572
-chr2 0 248 4 750 . 0 224
-chr2 224 372 5 0 . 248 348
-chr2 348 496 6 500 . 372 472
-chr2 472 572 7 0 . 496 572
+chr1 224 372 1 1000 . 248 348
+chr1 348 496 2 1000 . 372 472
+chr1 472 572 3 1000 . 496 572
+chr2 0 248 4 1000 . 0 224
+chr2 224 372 5 750 . 248 348
+chr2 348 496 6 1000 . 372 472
+chr2 472 572 7 1000 . 496 572
`)
}