Use tilelib to load for export and pca.

[lightning.git] / stats.go

diff --git a/stats.go b/stats.go
index 0fcb897afae156aa1952b56ec4909100e5be85c7..8da2a69ad75afb3ebfe72e28aa2f3d59c80bd792 100644 (file)
--- a/stats.go
+++ b/stats.go
@@ -17,11 +17,11 @@ import (
        log "github.com/sirupsen/logrus"
 )
 
-type stats struct {
+type statscmd struct {
        debugUnplaced bool
 }
 
-func (cmd *stats) RunCommand(prog string, args []string, stdin io.Reader, stdout, stderr io.Writer) int {
+func (cmd *statscmd) RunCommand(prog string, args []string, stdin io.Reader, stdout, stderr io.Writer) int {
        var err error
        defer func() {
                if err != nil {
@@ -101,7 +101,10 @@ func (cmd *stats) RunCommand(prog string, args []string, stdin io.Reader, stdout
        }
 
        bufw := bufio.NewWriter(output)
-       cmd.doStats(input, bufw)
+       err = cmd.doStats(input, bufw)
+       if err != nil {
+               return 1
+       }
        err = bufw.Flush()
        if err != nil {
                return 1
@@ -113,9 +116,10 @@ func (cmd *stats) RunCommand(prog string, args []string, stdin io.Reader, stdout
        return 0
 }
 
-func (cmd *stats) doStats(input io.Reader, output io.Writer) error {
+func (cmd *statscmd) doStats(input io.Reader, output io.Writer) error {
        var ret struct {
                Genomes          int
+               CalledBases      []int64
                Tags             int
                TagsPlacedNTimes []int // a[x]==y means there were y tags that placed x times
                TileVariants     int
@@ -126,6 +130,7 @@ func (cmd *stats) doStats(input io.Reader, output io.Writer) error {
 
        var tagSet [][]byte
        var tagPlacements []int
+       tileVariantCalls := map[tileLibRef]int{}
        dec := gob.NewDecoder(bufio.NewReaderSize(input, 1<<26))
        for {
                var ent LibraryEntry
@@ -133,7 +138,7 @@ func (cmd *stats) doStats(input io.Reader, output io.Writer) error {
                if err == io.EOF {
                        break
                } else if err != nil {
-                       return err
+                       return fmt.Errorf("gob decode: %w", err)
                }
                ret.Genomes += len(ent.CompactGenomes)
                ret.TileVariants += len(ent.TileVariants)
@@ -144,25 +149,18 @@ func (cmd *stats) doStats(input io.Reader, output io.Writer) error {
                        ret.Tags = len(ent.TagSet)
                        tagSet = ent.TagSet
                }
-               for _, g := range ent.CompactGenomes {
-                       if need := (len(g.Variants)+1)/2 - len(tagPlacements); need > 0 {
-                               tagPlacements = append(tagPlacements, make([]int, need)...)
-                       }
-                       for idx, v := range g.Variants {
-                               if v > 0 {
-                                       tagPlacements[idx/2]++
-                               }
-                       }
-               }
                for _, tv := range ent.TileVariants {
                        if need := 1 + len(tv.Sequence) - len(ret.VariantsBySize); need > 0 {
                                ret.VariantsBySize = append(ret.VariantsBySize, make([]int, need)...)
                                ret.NCVariantsBySize = append(ret.NCVariantsBySize, make([]int, need)...)
                        }
 
+                       calls := 0
                        hasNoCalls := false
                        for _, b := range tv.Sequence {
-                               if b != 'a' && b != 'c' && b != 'g' && b != 't' {
+                               if b == 'a' || b == 'c' || b == 'g' || b == 't' {
+                                       calls++
+                               } else {
                                        hasNoCalls = true
                                }
                        }
@@ -172,6 +170,21 @@ func (cmd *stats) doStats(input io.Reader, output io.Writer) error {
                        } else {
                                ret.VariantsBySize[len(tv.Sequence)]++
                        }
+
+                       tileVariantCalls[tileLibRef{Tag: tv.Tag, Variant: tv.Variant}] = calls
+               }
+               for _, g := range ent.CompactGenomes {
+                       if need := (len(g.Variants)+1)/2 - len(tagPlacements); need > 0 {
+                               tagPlacements = append(tagPlacements, make([]int, need)...)
+                       }
+                       calledBases := int64(0)
+                       for idx, v := range g.Variants {
+                               if v > 0 {
+                                       tagPlacements[idx/2]++
+                                       calledBases += int64(tileVariantCalls[tileLibRef{Tag: tagID(idx / 2), Variant: v}])
+                               }
+                       }
+                       ret.CalledBases = append(ret.CalledBases, calledBases)
                }
        }
        for id, p := range tagPlacements {