Fix low-coverage tiles counting toward min coverage threshold.
[lightning.git] / export.go
1 // Copyright (C) The Lightning Authors. All rights reserved.
2 //
3 // SPDX-License-Identifier: AGPL-3.0
4
5 package lightning
6
7 import (
8         "bufio"
9         "bytes"
10         "context"
11         "errors"
12         "flag"
13         "fmt"
14         "io"
15         "net/http"
16         _ "net/http/pprof"
17         "os"
18         "path/filepath"
19         "runtime"
20         "sort"
21         "strconv"
22         "strings"
23         "sync"
24         "time"
25
26         "git.arvados.org/arvados.git/sdk/go/arvados"
27         "github.com/arvados/lightning/hgvs"
28         "github.com/klauspost/pgzip"
29         "github.com/kshedden/gonpy"
30         "github.com/sirupsen/logrus"
31         log "github.com/sirupsen/logrus"
32 )
33
34 type tvVariant struct {
35         hgvs.Variant
36         librefs map[tileLibRef]bool
37 }
38
39 type outputFormat interface {
40         Filename() string
41         PadLeft() bool
42         Head(out io.Writer, cgs []CompactGenome, cases []bool, p float64) error
43         Print(out io.Writer, seqname string, varslice []tvVariant) error
44         Finish(outdir string, out io.Writer, seqname string) error
45         MaxGoroutines() int
46 }
47
48 var outputFormats = map[string]func() outputFormat{
49         "hgvs-numpy": func() outputFormat {
50                 return &formatHGVSNumpy{alleles: map[string][][]int8{}}
51         },
52         "hgvs-onehot": func() outputFormat { return formatHGVSOneHot{} },
53         "hgvs":        func() outputFormat { return formatHGVS{} },
54         "pvcf":        func() outputFormat { return formatPVCF{} },
55         "vcf":         func() outputFormat { return formatVCF{} },
56 }
57
58 type exporter struct {
59         outputFormat   outputFormat
60         outputPerChrom bool
61         compress       bool
62         maxTileSize    int
63         filter         filter
64         maxPValue      float64
65         cases          []bool
66 }
67
68 func (cmd *exporter) RunCommand(prog string, args []string, stdin io.Reader, stdout, stderr io.Writer) int {
69         var err error
70         defer func() {
71                 if err != nil {
72                         fmt.Fprintf(stderr, "%s\n", err)
73                 }
74         }()
75         flags := flag.NewFlagSet("", flag.ContinueOnError)
76         flags.SetOutput(stderr)
77         pprof := flags.String("pprof", "", "serve Go profile data at http://`[addr]:port`")
78         pprofdir := flags.String("pprof-dir", "", "write Go profile data to `directory` periodically")
79         runlocal := flags.Bool("local", false, "run on local host (default: run in an arvados container)")
80         projectUUID := flags.String("project", "", "project `UUID` for output data")
81         priority := flags.Int("priority", 500, "container request priority")
82         refname := flags.String("ref", "", "reference genome `name`")
83         inputDir := flags.String("input-dir", ".", "input `directory`")
84         cases := flags.String("cases", "", "file indicating which genomes are positive cases (for computing p-values)")
85         flags.Float64Var(&cmd.maxPValue, "p-value", 1, "do chi square test and omit columns with p-value above this threshold")
86         outputDir := flags.String("output-dir", ".", "output `directory`")
87         outputFormatStr := flags.String("output-format", "hgvs", "output `format`: hgvs, pvcf, or vcf")
88         outputBed := flags.String("output-bed", "", "also output bed `file`")
89         flags.BoolVar(&cmd.outputPerChrom, "output-per-chromosome", true, "output one file per chromosome")
90         flags.BoolVar(&cmd.compress, "z", false, "write gzip-compressed output files")
91         labelsFilename := flags.String("output-labels", "", "also output genome labels csv `file`")
92         flags.IntVar(&cmd.maxTileSize, "max-tile-size", 50000, "don't try to make annotations for tiles bigger than given `size`")
93         cmd.filter.Flags(flags)
94         err = flags.Parse(args)
95         if err == flag.ErrHelp {
96                 err = nil
97                 return 0
98         } else if err != nil {
99                 return 2
100         }
101         if flag.NArg() > 0 {
102                 err = fmt.Errorf("extra unparsed command line arguments: %q", flag.Args())
103                 return 2
104         }
105
106         if f, ok := outputFormats[*outputFormatStr]; !ok {
107                 err = fmt.Errorf("invalid output format %q", *outputFormatStr)
108                 return 2
109         } else {
110                 cmd.outputFormat = f()
111         }
112
113         if *pprof != "" {
114                 go func() {
115                         log.Println(http.ListenAndServe(*pprof, nil))
116                 }()
117         }
118         if *pprofdir != "" {
119                 go writeProfilesPeriodically(*pprofdir)
120         }
121
122         if !*runlocal {
123                 if *outputDir != "." {
124                         err = errors.New("cannot specify output directory in container mode: not implemented")
125                         return 1
126                 }
127                 runner := arvadosContainerRunner{
128                         Name:        "lightning export",
129                         Client:      arvados.NewClientFromEnv(),
130                         ProjectUUID: *projectUUID,
131                         RAM:         750000000000,
132                         VCPUs:       96,
133                         Priority:    *priority,
134                         APIAccess:   true,
135                 }
136                 err = runner.TranslatePaths(inputDir, cases)
137                 if err != nil {
138                         return 1
139                 }
140                 if *outputBed != "" {
141                         if strings.Contains(*outputBed, "/") {
142                                 err = fmt.Errorf("cannot use -output-bed filename %q containing '/' char", *outputBed)
143                                 return 1
144                         }
145                         *outputBed = "/mnt/output/" + *outputBed
146                 }
147                 runner.Args = []string{"export", "-local=true",
148                         "-pprof", ":6000",
149                         "-pprof-dir", "/mnt/output",
150                         "-ref", *refname,
151                         "-cases", *cases,
152                         "-p-value", fmt.Sprintf("%f", cmd.maxPValue),
153                         "-output-format", *outputFormatStr,
154                         "-output-bed", *outputBed,
155                         "-output-labels", "/mnt/output/labels.csv",
156                         "-output-per-chromosome=" + fmt.Sprintf("%v", cmd.outputPerChrom),
157                         "-max-tile-size", fmt.Sprintf("%d", cmd.maxTileSize),
158                         "-input-dir", *inputDir,
159                         "-output-dir", "/mnt/output",
160                         "-z=" + fmt.Sprintf("%v", cmd.compress),
161                 }
162                 runner.Args = append(runner.Args, cmd.filter.Args()...)
163                 var output string
164                 output, err = runner.Run()
165                 if err != nil {
166                         return 1
167                 }
168                 fmt.Fprintln(stdout, output)
169                 return 0
170         }
171
172         var cgs []CompactGenome
173         tilelib := &tileLibrary{
174                 retainNoCalls:       true,
175                 retainTileSequences: true,
176                 compactGenomes:      map[string][]tileVariantID{},
177         }
178         err = tilelib.LoadDir(context.Background(), *inputDir)
179         if err != nil {
180                 return 1
181         }
182
183         refseq, ok := tilelib.refseqs[*refname]
184         if !ok {
185                 err = fmt.Errorf("reference name %q not found in input; have %v", *refname, func() (names []string) {
186                         for name := range tilelib.refseqs {
187                                 names = append(names, name)
188                         }
189                         return
190                 }())
191                 return 1
192         }
193
194         log.Infof("filtering: %+v", cmd.filter)
195         cmd.filter.Apply(tilelib)
196
197         names := cgnames(tilelib)
198         for _, name := range names {
199                 cgs = append(cgs, CompactGenome{Name: name, Variants: tilelib.compactGenomes[name]})
200         }
201         if *labelsFilename != "" {
202                 log.Infof("writing labels to %s", *labelsFilename)
203                 var f *os.File
204                 f, err = os.OpenFile(*labelsFilename, os.O_CREATE|os.O_WRONLY, 0777)
205                 if err != nil {
206                         return 1
207                 }
208                 defer f.Close()
209                 for i, name := range names {
210                         _, err = fmt.Fprintf(f, "%d,%q,%q\n", i, trimFilenameForLabel(name), cmd.outputFormat.Filename())
211                         if err != nil {
212                                 err = fmt.Errorf("write %s: %w", *labelsFilename, err)
213                                 return 1
214                         }
215                 }
216                 err = f.Close()
217                 if err != nil {
218                         err = fmt.Errorf("close %s: %w", *labelsFilename, err)
219                         return 1
220                 }
221         }
222
223         cmd.cases = make([]bool, len(names))
224         if *cases != "" {
225                 log.Infof("reading cases file: %s", *cases)
226                 var f io.ReadCloser
227                 f, err = open(*cases)
228                 if err != nil {
229                         return 1
230                 }
231                 defer f.Close()
232                 var buf []byte
233                 buf, err = io.ReadAll(f)
234                 if err != nil {
235                         return 1
236                 }
237                 for _, pattern := range bytes.Split(buf, []byte("\n")) {
238                         if len(pattern) == 0 {
239                                 continue
240                         }
241                         pattern := string(pattern)
242                         idx := -1
243                         for i, name := range names {
244                                 if !strings.Contains(name, pattern) {
245                                         continue
246                                 } else if idx >= 0 {
247                                         err = fmt.Errorf("pattern %q in cases file matches multiple genome IDs: %q, %q", pattern, names[idx], name)
248                                         return 1
249                                 } else {
250                                         idx = i
251                                 }
252                         }
253                         if idx < 0 {
254                                 log.Warnf("pattern %q in cases file does not match any genome IDs", pattern)
255                                 continue
256                         }
257                         cmd.cases[idx] = true
258                 }
259         }
260
261         var bedout io.Writer
262         var bedfile *os.File
263         var bedbufw *bufio.Writer
264         if *outputBed != "" {
265                 bedfile, err = os.OpenFile(*outputBed, os.O_CREATE|os.O_WRONLY, 0666)
266                 if err != nil {
267                         return 1
268                 }
269                 defer bedfile.Close()
270                 bedbufw = bufio.NewWriterSize(bedfile, 16*1024*1024)
271                 bedout = bedbufw
272         }
273
274         err = cmd.export(*outputDir, bedout, tilelib, refseq, cgs)
275         if err != nil {
276                 return 1
277         }
278         if bedout != nil {
279                 err = bedbufw.Flush()
280                 if err != nil {
281                         return 1
282                 }
283                 err = bedfile.Close()
284                 if err != nil {
285                         return 1
286                 }
287         }
288         return 0
289 }
290
291 func (cmd *exporter) export(outdir string, bedout io.Writer, tilelib *tileLibrary, refseq map[string][]tileLibRef, cgs []CompactGenome) error {
292         var seqnames []string
293         var missing []tileLibRef
294         for seqname, librefs := range refseq {
295                 seqnames = append(seqnames, seqname)
296                 for _, libref := range librefs {
297                         if libref.Variant != 0 && tilelib.TileVariantSequence(libref) == nil {
298                                 missing = append(missing, libref)
299                         }
300                 }
301         }
302         sort.Strings(seqnames)
303
304         if len(missing) > 0 {
305                 if limit := 100; len(missing) > limit {
306                         log.Warnf("first %d missing tiles: %v", limit, missing[:limit])
307                 } else {
308                         log.Warnf("missing tiles: %v", missing)
309                 }
310                 return fmt.Errorf("%d needed tiles are missing from library", len(missing))
311         }
312
313         outw := make([]io.WriteCloser, len(seqnames))
314         bedw := make([]io.WriteCloser, len(seqnames))
315
316         var merges sync.WaitGroup
317         merge := func(dst io.Writer, src []io.WriteCloser, label string) {
318                 var mtx sync.Mutex
319                 for i, seqname := range seqnames {
320                         pr, pw := io.Pipe()
321                         src[i] = pw
322                         merges.Add(1)
323                         seqname := seqname
324                         go func() {
325                                 defer merges.Done()
326                                 log.Infof("writing %s %s", seqname, label)
327                                 scanner := bufio.NewScanner(pr)
328                                 for scanner.Scan() {
329                                         mtx.Lock()
330                                         dst.Write(scanner.Bytes())
331                                         dst.Write([]byte{'\n'})
332                                         mtx.Unlock()
333                                 }
334                                 log.Infof("writing %s %s done", seqname, label)
335                         }()
336                 }
337         }
338         if cmd.outputPerChrom {
339                 for i, seqname := range seqnames {
340                         fnm := filepath.Join(outdir, strings.Replace(cmd.outputFormat.Filename(), ".", "."+seqname+".", 1))
341                         if cmd.compress {
342                                 fnm += ".gz"
343                         }
344                         f, err := os.OpenFile(fnm, os.O_CREATE|os.O_WRONLY|os.O_TRUNC, 0666)
345                         if err != nil {
346                                 return err
347                         }
348                         defer f.Close()
349                         log.Infof("writing %q", f.Name())
350                         outw[i] = f
351                         if cmd.compress {
352                                 z := pgzip.NewWriter(f)
353                                 defer z.Close()
354                                 outw[i] = z
355                         }
356                         err = cmd.outputFormat.Head(outw[i], cgs, cmd.cases, cmd.maxPValue)
357                         if err != nil {
358                                 return err
359                         }
360                 }
361         } else {
362                 fnm := filepath.Join(outdir, cmd.outputFormat.Filename())
363                 if cmd.compress {
364                         fnm += ".gz"
365                 }
366                 f, err := os.OpenFile(fnm, os.O_CREATE|os.O_WRONLY|os.O_TRUNC, 0666)
367                 if err != nil {
368                         return err
369                 }
370                 defer f.Close()
371                 log.Infof("writing %q", fnm)
372                 var out io.Writer = f
373                 if cmd.compress {
374                         z := pgzip.NewWriter(out)
375                         defer z.Close()
376                         out = z
377                 }
378                 cmd.outputFormat.Head(out, cgs, cmd.cases, cmd.maxPValue)
379                 merge(out, outw, "output")
380         }
381         if bedout != nil {
382                 merge(bedout, bedw, "bed")
383         }
384
385         throttle := throttle{Max: runtime.NumCPU()}
386         if max := cmd.outputFormat.MaxGoroutines(); max > 0 {
387                 throttle.Max = max
388         }
389         log.Infof("assembling %d sequences in %d goroutines", len(seqnames), throttle.Max)
390         for seqidx, seqname := range seqnames {
391                 seqidx, seqname := seqidx, seqname
392                 outw := outw[seqidx]
393                 bedw := bedw[seqidx]
394                 throttle.Acquire()
395                 go func() {
396                         defer throttle.Release()
397                         if bedw != nil {
398                                 defer bedw.Close()
399                         }
400                         outwb := bufio.NewWriterSize(outw, 8*1024*1024)
401                         eachVariant(bedw, tilelib.taglib.keylen, seqname, refseq[seqname], tilelib, cgs, cmd.outputFormat.PadLeft(), cmd.maxTileSize, func(varslice []tvVariant) {
402                                 err := cmd.outputFormat.Print(outwb, seqname, varslice)
403                                 throttle.Report(err)
404                         })
405                         err := cmd.outputFormat.Finish(outdir, outwb, seqname)
406                         throttle.Report(err)
407                         err = outwb.Flush()
408                         throttle.Report(err)
409                         err = outw.Close()
410                         throttle.Report(err)
411                 }()
412         }
413
414         merges.Wait()
415         throttle.Wait()
416         return throttle.Err()
417 }
418
419 // Align genome tiles to reference tiles, call callback func on each
420 // variant, and (if bedw is not nil) write tile coverage to bedw.
421 func eachVariant(bedw io.Writer, taglen int, seqname string, reftiles []tileLibRef, tilelib *tileLibrary, cgs []CompactGenome, padLeft bool, maxTileSize int, callback func(varslice []tvVariant)) {
422         t0 := time.Now()
423         progressbar := time.NewTicker(time.Minute)
424         defer progressbar.Stop()
425         var outmtx sync.Mutex
426         defer outmtx.Lock()
427         refpos := 0
428         variantAt := map[int][]tvVariant{} // variantAt[chromOffset][genomeIndex*2+phase]
429         for refstep, libref := range reftiles {
430                 select {
431                 case <-progressbar.C:
432                         var eta interface{}
433                         if refstep > 0 {
434                                 fin := t0.Add(time.Duration(float64(time.Now().Sub(t0)) * float64(len(reftiles)) / float64(refstep)))
435                                 eta = fmt.Sprintf("%v (%v)", fin.Format(time.RFC3339), fin.Sub(time.Now()))
436                         } else {
437                                 eta = "N/A"
438                         }
439                         log.Printf("exportSeq: %s: refstep %d of %d, %.0f/s, ETA %v", seqname, refstep, len(reftiles), float64(refstep)/time.Now().Sub(t0).Seconds(), eta)
440                 default:
441                 }
442                 diffs := map[tileLibRef][]hgvs.Variant{}
443                 refseq := tilelib.TileVariantSequence(libref)
444                 tagcoverage := 0 // number of times the start tag was found in genomes -- max is len(cgs)*2
445                 for cgidx, cg := range cgs {
446                         for phase := 0; phase < 2; phase++ {
447                                 var variant tileVariantID
448                                 if i := int(libref.Tag)*2 + phase; len(cg.Variants) > i {
449                                         variant = cg.Variants[i]
450                                 }
451                                 if variant > 0 {
452                                         tagcoverage++
453                                 }
454                                 if variant == libref.Variant || variant == 0 {
455                                         continue
456                                 }
457                                 glibref := tileLibRef{Tag: libref.Tag, Variant: variant}
458                                 vars, ok := diffs[glibref]
459                                 if !ok {
460                                         genomeseq := tilelib.TileVariantSequence(glibref)
461                                         if len(genomeseq) == 0 {
462                                                 // Hash is known but sequence
463                                                 // is not, e.g., retainNoCalls
464                                                 // was false during import
465                                                 continue
466                                         }
467                                         if len(genomeseq) > maxTileSize {
468                                                 continue
469                                         }
470                                         refSequence := refseq
471                                         // If needed, extend the
472                                         // reference sequence up to
473                                         // the tag at the end of the
474                                         // genomeseq sequence.
475                                         refstepend := refstep + 1
476                                         for refstepend < len(reftiles) && len(refSequence) >= taglen && !bytes.EqualFold(refSequence[len(refSequence)-taglen:], genomeseq[len(genomeseq)-taglen:]) && len(refSequence) <= maxTileSize {
477                                                 if &refSequence[0] == &refseq[0] {
478                                                         refSequence = append([]byte(nil), refSequence...)
479                                                 }
480                                                 refSequence = append(refSequence, tilelib.TileVariantSequence(reftiles[refstepend])...)
481                                                 refstepend++
482                                         }
483                                         // (TODO: handle no-calls)
484                                         if len(refSequence) <= maxTileSize {
485                                                 refstr := strings.ToUpper(string(refSequence))
486                                                 genomestr := strings.ToUpper(string(genomeseq))
487                                                 vars, _ = hgvs.Diff(refstr, genomestr, time.Second)
488                                         }
489                                         diffs[glibref] = vars
490                                 }
491                                 for _, v := range vars {
492                                         if padLeft {
493                                                 v = v.PadLeft()
494                                         }
495                                         v.Position += refpos
496                                         varslice := variantAt[v.Position]
497                                         if varslice == nil {
498                                                 varslice = make([]tvVariant, len(cgs)*2)
499                                                 variantAt[v.Position] = varslice
500                                         }
501                                         varslice[cgidx*2+phase].Variant = v
502                                         if varslice[cgidx*2+phase].librefs == nil {
503                                                 varslice[cgidx*2+phase].librefs = map[tileLibRef]bool{glibref: true}
504                                         } else {
505                                                 varslice[cgidx*2+phase].librefs[glibref] = true
506                                         }
507                                 }
508                         }
509                 }
510                 refpos += len(refseq) - taglen
511
512                 // Flush entries from variantAt that are behind
513                 // refpos. Flush all entries if this is the last
514                 // reftile of the path/chromosome.
515                 flushpos := make([]int, 0, len(variantAt))
516                 lastrefstep := refstep == len(reftiles)-1
517                 for pos := range variantAt {
518                         if lastrefstep || pos <= refpos {
519                                 flushpos = append(flushpos, pos)
520                         }
521                 }
522                 sort.Slice(flushpos, func(i, j int) bool { return flushpos[i] < flushpos[j] })
523                 flushvariants := make([][]tvVariant, len(flushpos))
524                 for i, pos := range flushpos {
525                         varslice := variantAt[pos]
526                         delete(variantAt, pos)
527                         // Check for uninitialized (zero-value)
528                         // elements in varslice
529                         for i := range varslice {
530                                 if varslice[i].Position != 0 {
531                                         // Not a zero-value element
532                                         continue
533                                 }
534                                 // Set the position so
535                                 // varslice[*].Position are all equal
536                                 varslice[i].Position = pos
537                                 // This could be either =ref or a
538                                 // missing/low-quality tile. Figure
539                                 // out which.
540                                 vidx := int(libref.Tag)*2 + i%2
541                                 if vidx >= len(cgs[i/2].Variants) {
542                                         // Missing tile.
543                                         varslice[i].New = "-"
544                                         continue
545                                 }
546                                 v := cgs[i/2].Variants[vidx]
547                                 if v < 1 || len(tilelib.TileVariantSequence(tileLibRef{Tag: libref.Tag, Variant: v})) == 0 {
548                                         // Missing/low-quality tile.
549                                         varslice[i].New = "-" // fasta "gap of indeterminate length"
550                                 }
551                         }
552                         flushvariants[i] = varslice
553                 }
554                 outmtx.Lock()
555                 go func() {
556                         defer outmtx.Unlock()
557                         for _, varslice := range flushvariants {
558                                 callback(varslice)
559                         }
560                 }()
561                 if bedw != nil && len(refseq) > 0 {
562                         tilestart := refpos - len(refseq) + taglen
563                         tileend := refpos
564                         if !lastrefstep {
565                                 tileend += taglen
566                         }
567                         thickstart := tilestart + taglen
568                         if refstep == 0 {
569                                 thickstart = 0
570                         }
571                         thickend := refpos
572
573                         // coverage score, 0 to 1000
574                         score := 1000
575                         if len(cgs) > 0 {
576                                 score = 1000 * tagcoverage / len(cgs) / 2
577                         }
578
579                         fmt.Fprintf(bedw, "%s %d %d %d %d . %d %d\n",
580                                 seqname, tilestart, tileend,
581                                 libref.Tag,
582                                 score,
583                                 thickstart, thickend)
584                 }
585         }
586 }
587
588 func bucketVarsliceByRef(varslice []tvVariant) map[string]map[string]int {
589         byref := map[string]map[string]int{}
590         for _, v := range varslice {
591                 if v.Ref == "" && v.New == "" {
592                         // =ref
593                         continue
594                 }
595                 if v.New == "-" {
596                         // no-call
597                         continue
598                 }
599                 alts := byref[v.Ref]
600                 if alts == nil {
601                         alts = map[string]int{}
602                         byref[v.Ref] = alts
603                 }
604                 alts[v.New]++
605         }
606         return byref
607 }
608
609 type formatVCF struct{}
610
611 func (formatVCF) MaxGoroutines() int                     { return 0 }
612 func (formatVCF) Filename() string                       { return "out.vcf" }
613 func (formatVCF) PadLeft() bool                          { return true }
614 func (formatVCF) Finish(string, io.Writer, string) error { return nil }
615 func (formatVCF) Head(out io.Writer, cgs []CompactGenome, cases []bool, p float64) error {
616         _, err := fmt.Fprint(out, "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n")
617         return err
618 }
619 func (formatVCF) Print(out io.Writer, seqname string, varslice []tvVariant) error {
620         for ref, alts := range bucketVarsliceByRef(varslice) {
621                 altslice := make([]string, 0, len(alts))
622                 for alt := range alts {
623                         altslice = append(altslice, alt)
624                 }
625                 sort.Strings(altslice)
626
627                 info := "AC="
628                 for i, a := range altslice {
629                         if i > 0 {
630                                 info += ","
631                         }
632                         info += strconv.Itoa(alts[a])
633                 }
634                 _, err := fmt.Fprintf(out, "%s\t%d\t.\t%s\t%s\t.\t.\t%s\n", seqname, varslice[0].Position, ref, strings.Join(altslice, ","), info)
635                 if err != nil {
636                         return err
637                 }
638         }
639         return nil
640 }
641
642 type formatPVCF struct{}
643
644 func (formatPVCF) MaxGoroutines() int                     { return 0 }
645 func (formatPVCF) Filename() string                       { return "out.vcf" }
646 func (formatPVCF) PadLeft() bool                          { return true }
647 func (formatPVCF) Finish(string, io.Writer, string) error { return nil }
648 func (formatPVCF) Head(out io.Writer, cgs []CompactGenome, cases []bool, p float64) error {
649         fmt.Fprintln(out, `##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">`)
650         fmt.Fprintf(out, "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT")
651         for _, cg := range cgs {
652                 fmt.Fprintf(out, "\t%s", cg.Name)
653         }
654         _, err := fmt.Fprintf(out, "\n")
655         return err
656 }
657
658 func (formatPVCF) Print(out io.Writer, seqname string, varslice []tvVariant) error {
659         for ref, alts := range bucketVarsliceByRef(varslice) {
660                 altslice := make([]string, 0, len(alts))
661                 for alt := range alts {
662                         altslice = append(altslice, alt)
663                 }
664                 sort.Strings(altslice)
665                 for i, a := range altslice {
666                         alts[a] = i + 1
667                 }
668                 _, err := fmt.Fprintf(out, "%s\t%d\t.\t%s\t%s\t.\t.\t.\tGT", seqname, varslice[0].Position, ref, strings.Join(altslice, ","))
669                 if err != nil {
670                         return err
671                 }
672                 for i := 0; i < len(varslice); i += 2 {
673                         v1, v2 := varslice[i], varslice[i+1]
674                         a1, a2 := alts[v1.New], alts[v2.New]
675                         if v1.Ref != ref {
676                                 // variant on allele 0 belongs on a
677                                 // different output line -- same
678                                 // chr,pos but different "ref" length
679                                 a1 = 0
680                         }
681                         if v2.Ref != ref {
682                                 a2 = 0
683                         }
684                         _, err := fmt.Fprintf(out, "\t%d/%d", a1, a2)
685                         if err != nil {
686                                 return err
687                         }
688                 }
689                 _, err = out.Write([]byte{'\n'})
690                 if err != nil {
691                         return err
692                 }
693         }
694         return nil
695 }
696
697 type formatHGVS struct{}
698
699 func (formatHGVS) MaxGoroutines() int                                                     { return 0 }
700 func (formatHGVS) Filename() string                                                       { return "out.tsv" }
701 func (formatHGVS) PadLeft() bool                                                          { return false }
702 func (formatHGVS) Head(out io.Writer, cgs []CompactGenome, cases []bool, p float64) error { return nil }
703 func (formatHGVS) Finish(string, io.Writer, string) error                                 { return nil }
704 func (formatHGVS) Print(out io.Writer, seqname string, varslice []tvVariant) error {
705         for i := 0; i < len(varslice)/2; i++ {
706                 if i > 0 {
707                         out.Write([]byte{'\t'})
708                 }
709                 var1, var2 := varslice[i*2], varslice[i*2+1]
710                 if var1.New == "-" || var2.New == "-" {
711                         _, err := out.Write([]byte{'N'})
712                         if err != nil {
713                                 return err
714                         }
715                         continue
716                 }
717                 if var1.Variant == var2.Variant {
718                         if var1.Ref == var1.New {
719                                 _, err := out.Write([]byte{'.'})
720                                 if err != nil {
721                                         return err
722                                 }
723                         } else {
724                                 _, err := fmt.Fprintf(out, "%s:g.%s", seqname, var1.String())
725                                 if err != nil {
726                                         return err
727                                 }
728                         }
729                 } else {
730                         _, err := fmt.Fprintf(out, "%s:g.[%s];[%s]", seqname, var1.String(), var2.String())
731                         if err != nil {
732                                 return err
733                         }
734                 }
735         }
736         _, err := out.Write([]byte{'\n'})
737         return err
738 }
739
740 type formatHGVSOneHot struct{}
741
742 func (formatHGVSOneHot) MaxGoroutines() int { return 0 }
743 func (formatHGVSOneHot) Filename() string   { return "out.tsv" }
744 func (formatHGVSOneHot) PadLeft() bool      { return false }
745 func (formatHGVSOneHot) Head(out io.Writer, cgs []CompactGenome, cases []bool, p float64) error {
746         return nil
747 }
748 func (formatHGVSOneHot) Finish(string, io.Writer, string) error { return nil }
749 func (formatHGVSOneHot) Print(out io.Writer, seqname string, varslice []tvVariant) error {
750         vars := map[hgvs.Variant]bool{}
751         for _, v := range varslice {
752                 if v.Ref != v.New {
753                         vars[v.Variant] = true
754                 }
755         }
756
757         // sort variants to ensure output is deterministic
758         sorted := make([]hgvs.Variant, 0, len(vars))
759         for v := range vars {
760                 sorted = append(sorted, v)
761         }
762         sort.Slice(sorted, func(a, b int) bool { return hgvs.Less(sorted[a], sorted[b]) })
763
764         for _, v := range sorted {
765                 if v.New == "-" {
766                         continue
767                 }
768                 fmt.Fprintf(out, "%s.%s", seqname, v.String())
769                 for i := 0; i < len(varslice); i += 2 {
770                         if varslice[i].Variant == v || varslice[i+1].Variant == v {
771                                 out.Write([]byte("\t1"))
772                         } else {
773                                 out.Write([]byte("\t0"))
774                         }
775                 }
776                 _, err := out.Write([]byte{'\n'})
777                 if err != nil {
778                         return err
779                 }
780         }
781         return nil
782 }
783
784 type formatHGVSNumpy struct {
785         sync.Mutex
786         writelock sync.Mutex
787         alleles   map[string][][]int8 // alleles[seqname][variantidx][genomeidx*2+phase]
788         cases     []bool
789         maxPValue float64
790 }
791
792 func (*formatHGVSNumpy) MaxGoroutines() int { return 4 }
793 func (*formatHGVSNumpy) Filename() string   { return "annotations.csv" }
794 func (*formatHGVSNumpy) PadLeft() bool      { return false }
795 func (f *formatHGVSNumpy) Head(out io.Writer, cgs []CompactGenome, cases []bool, p float64) error {
796         f.cases = cases
797         f.maxPValue = p
798         return nil
799 }
800 func (f *formatHGVSNumpy) Print(outw io.Writer, seqname string, varslice []tvVariant) error {
801         // sort variants to ensure output is deterministic
802         sorted := make([]hgvs.Variant, 0, len(varslice))
803         for _, v := range varslice {
804                 sorted = append(sorted, v.Variant)
805         }
806         sort.Slice(sorted, func(a, b int) bool { return hgvs.Less(sorted[a], sorted[b]) })
807
808         f.Lock()
809         seqalleles := f.alleles[seqname]
810         f.Unlock()
811
812         chi2x := make([]bool, 0, len(varslice))
813         chi2y := make([]bool, 0, len(varslice))
814
815         // append a row to seqalleles for each unique non-ref variant
816         // in varslice.
817         var previous hgvs.Variant
818         for _, v := range sorted {
819                 if previous == v || v.Ref == v.New || v.New == "-" {
820                         continue
821                 }
822                 previous = v
823                 chi2x, chi2y := chi2x, chi2y
824                 newrow := make([]int8, len(varslice))
825                 for i, allele := range varslice {
826                         if allele.Variant == v {
827                                 newrow[i] = 1
828                                 chi2x = append(chi2x, true)
829                                 chi2y = append(chi2y, f.cases[i/2])
830                         } else if allele.Variant.New == "-" {
831                                 newrow[i] = -1
832                         } else {
833                                 chi2x = append(chi2x, false)
834                                 chi2y = append(chi2y, f.cases[i/2])
835                         }
836                 }
837                 if f.maxPValue < 1 && pvalue(chi2x, chi2y) > f.maxPValue {
838                         continue
839                 }
840                 seqalleles = append(seqalleles, newrow)
841                 _, err := fmt.Fprintf(outw, "%d,%q\n", len(seqalleles)-1, seqname+"."+v.String())
842                 if err != nil {
843                         return err
844                 }
845         }
846
847         f.Lock()
848         f.alleles[seqname] = seqalleles
849         f.Unlock()
850         return nil
851 }
852 func (f *formatHGVSNumpy) Finish(outdir string, _ io.Writer, seqname string) error {
853         // Write seqname's data to a .npy matrix with one row per
854         // genome and 2 columns per variant.
855         f.Lock()
856         seqalleles := f.alleles[seqname]
857         delete(f.alleles, seqname)
858         f.Unlock()
859         if len(seqalleles) == 0 {
860                 return nil
861         }
862         out := make([]int8, len(seqalleles)*len(seqalleles[0]))
863         rows := len(seqalleles[0]) / 2
864         cols := len(seqalleles) * 2
865         // copy seqalleles[varidx][genome*2+phase] to
866         // out[genome*nvars*2 + varidx*2 + phase]
867         for varidx, alleles := range seqalleles {
868                 for g := 0; g < len(alleles)/2; g++ {
869                         aa, ab := alleles[g*2], alleles[g*2+1]
870                         if aa < 0 || ab < 0 {
871                                 // no-call
872                                 out[g*cols+varidx*2] = -1
873                                 out[g*cols+varidx*2+1] = -1
874                         } else if aa > 0 && ab > 0 {
875                                 // hom
876                                 out[g*cols+varidx*2] = 1
877                         } else if aa > 0 || ab > 0 {
878                                 // het
879                                 out[g*cols+varidx*2+1] = 1
880                         }
881                 }
882         }
883         outf, err := os.OpenFile(outdir+"/matrix."+seqname+".npy", os.O_CREATE|os.O_EXCL|os.O_WRONLY, 0777)
884         if err != nil {
885                 return err
886         }
887         defer outf.Close()
888         bufw := bufio.NewWriter(outf)
889         npw, err := gonpy.NewWriter(nopCloser{bufw})
890         if err != nil {
891                 return err
892         }
893         log.WithFields(logrus.Fields{
894                 "seqname": seqname,
895                 "rows":    rows,
896                 "cols":    cols,
897         }).Info("writing numpy")
898         npw.Shape = []int{rows, cols}
899         f.writelock.Lock() // serialize because WriteInt8 uses lots of memory
900         npw.WriteInt8(out)
901         f.writelock.Unlock()
902         err = bufw.Flush()
903         if err != nil {
904                 return err
905         }
906         err = outf.Close()
907         if err != nil {
908                 return err
909         }
910         return nil
911 }