Mention undiffable variants in annotations, write -2 in hgvs matrix.

[lightning.git] / slicenumpy.go

diff --git a/slicenumpy.go b/slicenumpy.go
index fdaa02e658de2a59dc6b1bbbbdaad0ecb250f1ef..a8fdeb064eefa682437b0558958785b003ed81fc 100644 (file)
--- a/slicenumpy.go
+++ b/slicenumpy.go
@@ -10,6 +10,7 @@ import (
        "flag"
        "fmt"
        "io"
+       "io/ioutil"
        "net/http"
        _ "net/http/pprof"
        "os"
@@ -206,11 +207,11 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s
                                return 1
                        }
                        if isdup[libref.Tag] {
-                               log.Printf("dropping reference tile %+v from %s, tag not unique", libref, seqname)
+                               log.Printf("dropping reference tile %+v from %s @ %d, tag not unique", libref, seqname, pos)
                        } else if reftile[libref.Tag] != nil {
-                               log.Printf("dropping reference tile %+v from %s, tag not unique", tileLibRef{Tag: libref.Tag, Variant: reftile[libref.Tag].variant}, reftile[libref.Tag].seqname)
+                               log.Printf("dropping reference tile %+v from %s @ %d, tag not unique", tileLibRef{Tag: libref.Tag, Variant: reftile[libref.Tag].variant}, reftile[libref.Tag].seqname, reftile[libref.Tag].pos)
                                delete(reftile, libref.Tag)
-                               log.Printf("dropping reference tile %+v from %s, tag not unique", libref, seqname)
+                               log.Printf("dropping reference tile %+v from %s @ %d, tag not unique", libref, seqname, pos)
                                isdup[libref.Tag] = true
                        } else {
                                reftile[libref.Tag] = &reftileinfo{
@@ -375,8 +376,8 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s
                                }
                                variants, ok := seq[tag]
                                if !ok {
-                                       outcol++
-                                       continue
+                                       // how could we even have a reftile if there is no sequence data??
+                                       return fmt.Errorf("bug: have no variants for tag %d but reftile is %+v", tag, rt)
                                }
                                reftilestr := strings.ToUpper(string(rt.tiledata))
                                remap := variantRemap[tag-tagstart]
@@ -389,9 +390,11 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s
                                                done[v] = true
                                        }
                                        if len(tv.Sequence) < taglen || !bytes.HasSuffix(rt.tiledata, tv.Sequence[len(tv.Sequence)-taglen:]) {
+                                               fmt.Fprintf(annow, "%d,%d,%d,.,%s,%d,.,,\n", tag, outcol, v, rt.seqname, rt.pos)
                                                continue
                                        }
                                        if lendiff := len(rt.tiledata) - len(tv.Sequence); lendiff < -1000 || lendiff > 1000 {
+                                               fmt.Fprintf(annow, "%d,%d,%d,,%s,%d,,,\n", tag, outcol, v, rt.seqname, rt.pos)
                                                continue
                                        }
                                        diffs, _ := hgvs.Diff(reftilestr, strings.ToUpper(string(tv.Sequence)), 0)
@@ -468,6 +471,7 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s
                        cols += len(chunk) / rows
                }
                out := make([]int16, rows*cols)
+               hgvs := map[string][2][]int16{} // hgvs -> [[g0,g1,g2,...], [g0,g1,g2,...]] (slice of genomes for each phase)
                startcol := 0
                for outIdx, chunk := range toMerge {
                        chunkcols := len(chunk) / rows
@@ -490,9 +494,41 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s
                                if len(line) == 0 {
                                        continue
                                }
-                               fields := bytes.SplitN(line, []byte{','}, 3)
+                               fields := bytes.SplitN(line, []byte{','}, 8)
                                incol, _ := strconv.Atoi(string(fields[1]))
-                               fmt.Fprintf(annow, "%s,%d,%s\n", fields[0], incol+startcol/2, fields[2])
+                               tileVariant, _ := strconv.Atoi(string(fields[2]))
+                               hgvsID := string(fields[3])
+                               seqname := string(fields[4])
+                               pos, _ := strconv.Atoi(string(fields[5]))
+                               refseq := fields[6]
+                               if mask != nil && !mask.Check(strings.TrimPrefix(seqname, "chr"), pos, pos+len(refseq)) {
+                                       // The tile intersects one of
+                                       // the selected regions, but
+                                       // this particular HGVS
+                                       // variant does not.
+                                       continue
+                               }
+                               fmt.Fprintf(annow, "%s,%d,%d,%s,%s,%d,%s,%s\n", fields[0], incol+startcol/2, tileVariant, hgvsID, seqname, pos, refseq, fields[7])
+                               hgvscols := hgvs[hgvsID]
+                               if hgvscols[0] == nil {
+                                       hgvscols = [2][]int16{make([]int16, len(cgnames)), make([]int16, len(cgnames))}
+                                       hgvs[hgvsID] = hgvscols
+                               }
+                               for ph := 0; ph < 2; ph++ {
+                                       for row := 0; row < rows; row++ {
+                                               v := chunk[row*chunkcols+incol*2+ph]
+                                               if int(v) == tileVariant {
+                                                       if len(hgvsID) == 0 {
+                                                               // we have the tile variant sequence, but the diff against ref didn't work out (see lendiff above)
+                                                               hgvscols[ph][row] = -2
+                                                       } else {
+                                                               hgvscols[ph][row] = 1
+                                                       }
+                                               } else if v < 0 {
+                                                       hgvscols[ph][row] = -1
+                                               }
+                                       }
+                               }
                        }
 
                        startcol += chunkcols
@@ -509,6 +545,37 @@ func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, s
                if err != nil {
                        return 1
                }
+               out = nil
+
+               cols = len(hgvs) * 2
+               log.Printf("building hgvs-based matrix: %d rows x %d cols", rows, cols)
+               out = make([]int16, rows*cols)
+               hgvsIDs := make([]string, 0, len(hgvs))
+               for hgvsID := range hgvs {
+                       hgvsIDs = append(hgvsIDs, hgvsID)
+               }
+               sort.Strings(hgvsIDs)
+               var hgvsLabels bytes.Buffer
+               for idx, hgvsID := range hgvsIDs {
+                       fmt.Fprintf(&hgvsLabels, "%d,%s\n", idx, hgvsID)
+                       for ph := 0; ph < 2; ph++ {
+                               hgvscol := hgvs[hgvsID][ph]
+                               for row, val := range hgvscol {
+                                       out[row*cols+idx*2+ph] = val
+                               }
+                       }
+               }
+               err = writeNumpyInt16(fmt.Sprintf("%s/hgvs.npy", *outputDir), out, rows, cols)
+               if err != nil {
+                       return 1
+               }
+
+               fnm := fmt.Sprintf("%s/hgvs.annotations.csv", *outputDir)
+               log.Printf("writing hgvs labels: %s", fnm)
+               err = ioutil.WriteFile(fnm, hgvsLabels.Bytes(), 0777)
+               if err != nil {
+                       return 1
+               }
        }
        return 0
 }