export.go

   1 package lightning
   2
   3 import (
   4         "bufio"
   5         "bytes"
   6         "context"
   7         "errors"
   8         "flag"
   9         "fmt"
  10         "io"
  11         "net/http"
  12         _ "net/http/pprof"
  13         "os"
  14         "path/filepath"
  15         "runtime"
  16         "sort"
  17         "strconv"
  18         "strings"
  19         "sync"
  20         "time"
  21
  22         "git.arvados.org/arvados.git/sdk/go/arvados"
  23         "github.com/arvados/lightning/hgvs"
  24         "github.com/klauspost/pgzip"
  25         log "github.com/sirupsen/logrus"
  26 )
  27
  28 type tvVariant struct {
  29         hgvs.Variant
  30         librefs map[tileLibRef]bool
  31 }
  32
  33 type outputFormat struct {
  34         Filename string
  35         Head     func(out io.Writer, cgs []CompactGenome)
  36         Print    func(out io.Writer, seqname string, varslice []tvVariant)
  37         PadLeft  bool
  38 }
  39
  40 var (
  41         outputFormats = map[string]outputFormat{
  42                 "hgvs-onehot": outputFormatHGVSOneHot,
  43                 "hgvs":        outputFormatHGVS,
  44                 "pvcf":        outputFormatPVCF,
  45                 "vcf":         outputFormatVCF,
  46         }
  47         outputFormatHGVS       = outputFormat{Filename: "out.tsv", Head: headNone, Print: printHGVS}
  48         outputFormatHGVSOneHot = outputFormat{Filename: "out.tsv", Head: headNone, Print: printHGVSOneHot}
  49         outputFormatPVCF       = outputFormat{Filename: "out.vcf", Head: headPVCF, Print: printPVCF, PadLeft: true}
  50         outputFormatVCF        = outputFormat{Filename: "out.vcf", Head: headVCF, Print: printVCF, PadLeft: true}
  51         headNone               = func(io.Writer, []CompactGenome) {}
  52 )
  53
  54 type exporter struct {
  55         outputFormat   outputFormat
  56         outputPerChrom bool
  57         compress       bool
  58         maxTileSize    int
  59 }
  60
  61 func (cmd *exporter) RunCommand(prog string, args []string, stdin io.Reader, stdout, stderr io.Writer) int {
  62         var err error
  63         defer func() {
  64                 if err != nil {
  65                         fmt.Fprintf(stderr, "%s\n", err)
  66                 }
  67         }()
  68         flags := flag.NewFlagSet("", flag.ContinueOnError)
  69         flags.SetOutput(stderr)
  70         pprof := flags.String("pprof", "", "serve Go profile data at http://`[addr]:port`")
  71         pprofdir := flags.String("pprof-dir", "", "write Go profile data to `directory` periodically")
  72         runlocal := flags.Bool("local", false, "run on local host (default: run in an arvados container)")
  73         projectUUID := flags.String("project", "", "project `UUID` for output data")
  74         priority := flags.Int("priority", 500, "container request priority")
  75         refname := flags.String("ref", "", "reference genome `name`")
  76         inputDir := flags.String("input-dir", ".", "input `directory`")
  77         outputDir := flags.String("output-dir", ".", "output `directory`")
  78         outputFormatStr := flags.String("output-format", "hgvs", "output `format`: hgvs, pvcf, or vcf")
  79         outputBed := flags.String("output-bed", "", "also output bed `file`")
  80         flags.BoolVar(&cmd.outputPerChrom, "output-per-chromosome", true, "output one file per chromosome")
  81         flags.BoolVar(&cmd.compress, "z", false, "write gzip-compressed output files")
  82         labelsFilename := flags.String("output-labels", "", "also output genome labels csv `file`")
  83         flags.IntVar(&cmd.maxTileSize, "max-tile-size", 50000, "don't try to make annotations for tiles bigger than given `size`")
  84         err = flags.Parse(args)
  85         if err == flag.ErrHelp {
  86                 err = nil
  87                 return 0
  88         } else if err != nil {
  89                 return 2
  90         }
  91         if flag.NArg() > 0 {
  92                 err = fmt.Errorf("extra unparsed command line arguments: %q", flag.Args())
  93                 return 2
  94         }
  95
  96         if f, ok := outputFormats[*outputFormatStr]; !ok {
  97                 err = fmt.Errorf("invalid output format %q", *outputFormatStr)
  98                 return 2
  99         } else {
 100                 cmd.outputFormat = f
 101         }
 102
 103         if *pprof != "" {
 104                 go func() {
 105                         log.Println(http.ListenAndServe(*pprof, nil))
 106                 }()
 107         }
 108         if *pprofdir != "" {
 109                 go writeProfilesPeriodically(*pprofdir)
 110         }
 111
 112         if !*runlocal {
 113                 if *outputDir != "." {
 114                         err = errors.New("cannot specify output directory in container mode: not implemented")
 115                         return 1
 116                 }
 117                 runner := arvadosContainerRunner{
 118                         Name:        "lightning export",
 119                         Client:      arvados.NewClientFromEnv(),
 120                         ProjectUUID: *projectUUID,
 121                         RAM:         700000000000,
 122                         VCPUs:       96,
 123                         Priority:    *priority,
 124                         APIAccess:   true,
 125                 }
 126                 err = runner.TranslatePaths(inputDir)
 127                 if err != nil {
 128                         return 1
 129                 }
 130                 if *outputBed != "" {
 131                         if strings.Contains(*outputBed, "/") {
 132                                 err = fmt.Errorf("cannot use -output-bed filename %q containing '/' char", *outputBed)
 133                                 return 1
 134                         }
 135                         *outputBed = "/mnt/output/" + *outputBed
 136                 }
 137                 runner.Args = []string{"export", "-local=true",
 138                         "-pprof", ":6000",
 139                         "-pprof-dir", "/mnt/output",
 140                         "-ref", *refname,
 141                         "-output-format", *outputFormatStr,
 142                         "-output-bed", *outputBed,
 143                         "-output-labels", "/mnt/output/labels.csv",
 144                         "-output-per-chromosome=" + fmt.Sprintf("%v", cmd.outputPerChrom),
 145                         "-max-tile-size", fmt.Sprintf("%d", cmd.maxTileSize),
 146                         "-input-dir", *inputDir,
 147                         "-output-dir", "/mnt/output",
 148                         "-z=" + fmt.Sprintf("%v", cmd.compress),
 149                 }
 150                 var output string
 151                 output, err = runner.Run()
 152                 if err != nil {
 153                         return 1
 154                 }
 155                 fmt.Fprintln(stdout, output)
 156                 return 0
 157         }
 158
 159         var cgs []CompactGenome
 160         tilelib := &tileLibrary{
 161                 retainNoCalls:       true,
 162                 retainTileSequences: true,
 163                 compactGenomes:      map[string][]tileVariantID{},
 164         }
 165         err = tilelib.LoadDir(context.Background(), *inputDir, nil)
 166         if err != nil {
 167                 return 1
 168         }
 169
 170         refseq, ok := tilelib.refseqs[*refname]
 171         if !ok {
 172                 err = fmt.Errorf("reference name %q not found in input; have %v", *refname, func() (names []string) {
 173                         for name := range tilelib.refseqs {
 174                                 names = append(names, name)
 175                         }
 176                         return
 177                 }())
 178                 return 1
 179         }
 180
 181         names := cgnames(tilelib)
 182         for _, name := range names {
 183                 cgs = append(cgs, CompactGenome{Name: name, Variants: tilelib.compactGenomes[name]})
 184         }
 185         if *labelsFilename != "" {
 186                 log.Infof("writing labels to %s", *labelsFilename)
 187                 var f *os.File
 188                 f, err = os.OpenFile(*labelsFilename, os.O_CREATE|os.O_WRONLY, 0777)
 189                 if err != nil {
 190                         return 1
 191                 }
 192                 defer f.Close()
 193                 for i, name := range names {
 194                         _, err = fmt.Fprintf(f, "%d,%q,%q\n", i, trimFilenameForLabel(name), cmd.outputFormat.Filename)
 195                         if err != nil {
 196                                 err = fmt.Errorf("write %s: %w", *labelsFilename, err)
 197                                 return 1
 198                         }
 199                 }
 200                 err = f.Close()
 201                 if err != nil {
 202                         err = fmt.Errorf("close %s: %w", *labelsFilename, err)
 203                         return 1
 204                 }
 205         }
 206
 207         var bedout io.Writer
 208         var bedfile *os.File
 209         var bedbufw *bufio.Writer
 210         if *outputBed != "" {
 211                 bedfile, err = os.OpenFile(*outputBed, os.O_CREATE|os.O_WRONLY, 0666)
 212                 if err != nil {
 213                         return 1
 214                 }
 215                 defer bedfile.Close()
 216                 bedbufw = bufio.NewWriterSize(bedfile, 16*1024*1024)
 217                 bedout = bedbufw
 218         }
 219
 220         err = cmd.export(*outputDir, bedout, tilelib, refseq, cgs)
 221         if err != nil {
 222                 return 1
 223         }
 224         if bedout != nil {
 225                 err = bedbufw.Flush()
 226                 if err != nil {
 227                         return 1
 228                 }
 229                 err = bedfile.Close()
 230                 if err != nil {
 231                         return 1
 232                 }
 233         }
 234         return 0
 235 }
 236
 237 func (cmd *exporter) export(outdir string, bedout io.Writer, tilelib *tileLibrary, refseq map[string][]tileLibRef, cgs []CompactGenome) error {
 238         var seqnames []string
 239         var missing []tileLibRef
 240         for seqname, librefs := range refseq {
 241                 seqnames = append(seqnames, seqname)
 242                 for _, libref := range librefs {
 243                         if libref.Variant != 0 && tilelib.TileVariantSequence(libref) == nil {
 244                                 missing = append(missing, libref)
 245                         }
 246                 }
 247         }
 248         sort.Strings(seqnames)
 249
 250         if len(missing) > 0 {
 251                 if limit := 100; len(missing) > limit {
 252                         log.Warnf("first %d missing tiles: %v", limit, missing[:limit])
 253                 } else {
 254                         log.Warnf("missing tiles: %v", missing)
 255                 }
 256                 return fmt.Errorf("%d needed tiles are missing from library", len(missing))
 257         }
 258
 259         outw := make([]io.WriteCloser, len(seqnames))
 260         bedw := make([]io.WriteCloser, len(seqnames))
 261
 262         var merges sync.WaitGroup
 263         merge := func(dst io.Writer, src []io.WriteCloser, label string) {
 264                 var mtx sync.Mutex
 265                 for i, seqname := range seqnames {
 266                         pr, pw := io.Pipe()
 267                         src[i] = pw
 268                         merges.Add(1)
 269                         seqname := seqname
 270                         go func() {
 271                                 defer merges.Done()
 272                                 log.Infof("writing %s %s", seqname, label)
 273                                 scanner := bufio.NewScanner(pr)
 274                                 for scanner.Scan() {
 275                                         mtx.Lock()
 276                                         dst.Write(scanner.Bytes())
 277                                         dst.Write([]byte{'\n'})
 278                                         mtx.Unlock()
 279                                 }
 280                                 log.Infof("writing %s %s done", seqname, label)
 281                         }()
 282                 }
 283         }
 284         if cmd.outputPerChrom {
 285                 for i, seqname := range seqnames {
 286                         fnm := filepath.Join(outdir, strings.Replace(cmd.outputFormat.Filename, ".", "."+seqname+".", 1))
 287                         if cmd.compress {
 288                                 fnm += ".gz"
 289                         }
 290                         f, err := os.OpenFile(fnm, os.O_CREATE|os.O_WRONLY|os.O_TRUNC, 0666)
 291                         if err != nil {
 292                                 return err
 293                         }
 294                         defer f.Close()
 295                         log.Infof("writing %q", f.Name())
 296                         outw[i] = f
 297                         if cmd.compress {
 298                                 z := pgzip.NewWriter(f)
 299                                 defer z.Close()
 300                                 outw[i] = z
 301                         }
 302                         cmd.outputFormat.Head(outw[i], cgs)
 303                 }
 304         } else {
 305                 fnm := filepath.Join(outdir, cmd.outputFormat.Filename)
 306                 if cmd.compress {
 307                         fnm += ".gz"
 308                 }
 309                 f, err := os.OpenFile(fnm, os.O_CREATE|os.O_WRONLY|os.O_TRUNC, 0666)
 310                 if err != nil {
 311                         return err
 312                 }
 313                 defer f.Close()
 314                 log.Infof("writing %q", fnm)
 315                 var out io.Writer = f
 316                 if cmd.compress {
 317                         z := pgzip.NewWriter(out)
 318                         defer z.Close()
 319                         out = z
 320                 }
 321                 cmd.outputFormat.Head(out, cgs)
 322                 merge(out, outw, "output")
 323         }
 324         if bedout != nil {
 325                 merge(bedout, bedw, "bed")
 326         }
 327
 328         throttle := throttle{Max: runtime.NumCPU()}
 329         log.Infof("assembling %d sequences in %d goroutines", len(seqnames), throttle.Max)
 330         for seqidx, seqname := range seqnames {
 331                 seqidx, seqname := seqidx, seqname
 332                 outw := outw[seqidx]
 333                 bedw := bedw[seqidx]
 334                 throttle.Acquire()
 335                 go func() {
 336                         defer throttle.Release()
 337                         if bedw != nil {
 338                                 defer bedw.Close()
 339                         }
 340                         outwb := bufio.NewWriterSize(outw, 8*1024*1024)
 341                         cmd.exportSeq(outwb, bedw, tilelib.taglib.keylen, seqname, refseq[seqname], tilelib, cgs)
 342                         err := outwb.Flush()
 343                         throttle.Report(err)
 344                         err = outw.Close()
 345                         throttle.Report(err)
 346                 }()
 347         }
 348
 349         merges.Wait()
 350         throttle.Wait()
 351         return throttle.Err()
 352 }
 353
 354 // Align genome tiles to reference tiles, write diffs to outw, and (if
 355 // bedw is not nil) write tile coverage to bedw.
 356 func (cmd *exporter) exportSeq(outw, bedw io.Writer, taglen int, seqname string, reftiles []tileLibRef, tilelib *tileLibrary, cgs []CompactGenome) {
 357         t0 := time.Now()
 358         progressbar := time.NewTicker(time.Minute)
 359         defer progressbar.Stop()
 360         var outmtx sync.Mutex
 361         defer outmtx.Lock()
 362         refpos := 0
 363         variantAt := map[int][]tvVariant{} // variantAt[chromOffset][genomeIndex*2+phase]
 364         for refstep, libref := range reftiles {
 365                 select {
 366                 case <-progressbar.C:
 367                         var eta interface{}
 368                         if refstep > 0 {
 369                                 fin := t0.Add(time.Duration(float64(time.Now().Sub(t0)) * float64(len(reftiles)) / float64(refstep)))
 370                                 eta = fmt.Sprintf("%v (%v)", fin.Format(time.RFC3339), fin.Sub(time.Now()))
 371                         } else {
 372                                 eta = "N/A"
 373                         }
 374                         log.Printf("exportSeq: %s: refstep %d of %d, %.0f/s, ETA %v", seqname, refstep, len(reftiles), float64(refstep)/time.Now().Sub(t0).Seconds(), eta)
 375                 default:
 376                 }
 377                 diffs := map[tileLibRef][]hgvs.Variant{}
 378                 refseq := tilelib.TileVariantSequence(libref)
 379                 tagcoverage := 0 // number of times the start tag was found in genomes -- max is len(cgs)*2
 380                 for cgidx, cg := range cgs {
 381                         for phase := 0; phase < 2; phase++ {
 382                                 if len(cg.Variants) <= int(libref.Tag)*2+phase {
 383                                         continue
 384                                 }
 385                                 variant := cg.Variants[int(libref.Tag)*2+phase]
 386                                 if variant == 0 {
 387                                         continue
 388                                 }
 389                                 tagcoverage++
 390                                 if variant == libref.Variant {
 391                                         continue
 392                                 }
 393                                 glibref := tileLibRef{Tag: libref.Tag, Variant: variant}
 394                                 vars, ok := diffs[glibref]
 395                                 if !ok {
 396                                         genomeseq := tilelib.TileVariantSequence(glibref)
 397                                         if len(genomeseq) == 0 {
 398                                                 // Hash is known but sequence
 399                                                 // is not, e.g., retainNoCalls
 400                                                 // was false during import
 401                                                 continue
 402                                         }
 403                                         if len(genomeseq) > cmd.maxTileSize {
 404                                                 continue
 405                                         }
 406                                         refSequence := refseq
 407                                         // If needed, extend the
 408                                         // reference sequence up to
 409                                         // the tag at the end of the
 410                                         // genomeseq sequence.
 411                                         refstepend := refstep + 1
 412                                         for refstepend < len(reftiles) && len(refSequence) >= taglen && !bytes.EqualFold(refSequence[len(refSequence)-taglen:], genomeseq[len(genomeseq)-taglen:]) && len(refSequence) <= cmd.maxTileSize {
 413                                                 if &refSequence[0] == &refseq[0] {
 414                                                         refSequence = append([]byte(nil), refSequence...)
 415                                                 }
 416                                                 refSequence = append(refSequence, tilelib.TileVariantSequence(reftiles[refstepend])...)
 417                                                 refstepend++
 418                                         }
 419                                         // (TODO: handle no-calls)
 420                                         if len(refSequence) <= cmd.maxTileSize {
 421                                                 refstr := strings.ToUpper(string(refSequence))
 422                                                 genomestr := strings.ToUpper(string(genomeseq))
 423                                                 vars, _ = hgvs.Diff(refstr, genomestr, time.Second)
 424                                         }
 425                                         diffs[glibref] = vars
 426                                 }
 427                                 for _, v := range vars {
 428                                         if cmd.outputFormat.PadLeft {
 429                                                 v = v.PadLeft()
 430                                         }
 431                                         v.Position += refpos
 432                                         varslice := variantAt[v.Position]
 433                                         if varslice == nil {
 434                                                 varslice = make([]tvVariant, len(cgs)*2)
 435                                                 variantAt[v.Position] = varslice
 436                                         }
 437                                         varslice[cgidx*2+phase].Variant = v
 438                                         if varslice[cgidx*2+phase].librefs == nil {
 439                                                 varslice[cgidx*2+phase].librefs = map[tileLibRef]bool{glibref: true}
 440                                         } else {
 441                                                 varslice[cgidx*2+phase].librefs[glibref] = true
 442                                         }
 443                                 }
 444                         }
 445                 }
 446                 refpos += len(refseq) - taglen
 447
 448                 // Flush entries from variantAt that are behind
 449                 // refpos. Flush all entries if this is the last
 450                 // reftile of the path/chromosome.
 451                 flushpos := make([]int, 0, len(variantAt))
 452                 lastrefstep := refstep == len(reftiles)-1
 453                 for pos := range variantAt {
 454                         if lastrefstep || pos <= refpos {
 455                                 flushpos = append(flushpos, pos)
 456                         }
 457                 }
 458                 sort.Slice(flushpos, func(i, j int) bool { return flushpos[i] < flushpos[j] })
 459                 flushvariants := make([][]tvVariant, len(flushpos))
 460                 for i, pos := range flushpos {
 461                         varslice := variantAt[pos]
 462                         delete(variantAt, pos)
 463                         for i := range varslice {
 464                                 if varslice[i].Position == 0 {
 465                                         varslice[i].Position = pos
 466                                 }
 467                         }
 468                         flushvariants[i] = varslice
 469                 }
 470                 outmtx.Lock()
 471                 go func() {
 472                         defer outmtx.Unlock()
 473                         for _, varslice := range flushvariants {
 474                                 cmd.outputFormat.Print(outw, seqname, varslice)
 475                         }
 476                 }()
 477                 if bedw != nil && len(refseq) > 0 {
 478                         tilestart := refpos - len(refseq) + taglen
 479                         tileend := refpos
 480                         if !lastrefstep {
 481                                 tileend += taglen
 482                         }
 483                         thickstart := tilestart + taglen
 484                         if refstep == 0 {
 485                                 thickstart = 0
 486                         }
 487                         thickend := refpos
 488
 489                         // coverage score, 0 to 1000
 490                         score := 1000
 491                         if len(cgs) > 0 {
 492                                 score = 1000 * tagcoverage / len(cgs) / 2
 493                         }
 494
 495                         fmt.Fprintf(bedw, "%s %d %d %d %d . %d %d\n",
 496                                 seqname, tilestart, tileend,
 497                                 libref.Tag,
 498                                 score,
 499                                 thickstart, thickend)
 500                 }
 501         }
 502 }
 503
 504 func bucketVarsliceByRef(varslice []tvVariant) map[string]map[string]int {
 505         byref := map[string]map[string]int{}
 506         for _, v := range varslice {
 507                 if v.Ref == "" && v.New == "" {
 508                         continue
 509                 }
 510                 alts := byref[v.Ref]
 511                 if alts == nil {
 512                         alts = map[string]int{}
 513                         byref[v.Ref] = alts
 514                 }
 515                 alts[v.New]++
 516         }
 517         return byref
 518 }
 519
 520 func headVCF(out io.Writer, cgs []CompactGenome) {
 521         fmt.Fprint(out, "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n")
 522 }
 523
 524 func printVCF(out io.Writer, seqname string, varslice []tvVariant) {
 525         for ref, alts := range bucketVarsliceByRef(varslice) {
 526                 altslice := make([]string, 0, len(alts))
 527                 for alt := range alts {
 528                         altslice = append(altslice, alt)
 529                 }
 530                 sort.Strings(altslice)
 531
 532                 info := "AC="
 533                 for i, a := range altslice {
 534                         if i > 0 {
 535                                 info += ","
 536                         }
 537                         info += strconv.Itoa(alts[a])
 538                 }
 539                 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)
 540         }
 541 }
 542
 543 func headPVCF(out io.Writer, cgs []CompactGenome) {
 544         fmt.Fprintln(out, `##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">`)
 545         fmt.Fprintf(out, "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT")
 546         for _, cg := range cgs {
 547                 fmt.Fprintf(out, "\t%s", cg.Name)
 548         }
 549         fmt.Fprintf(out, "\n")
 550 }
 551
 552 func printPVCF(out io.Writer, seqname string, varslice []tvVariant) {
 553         for ref, alts := range bucketVarsliceByRef(varslice) {
 554                 altslice := make([]string, 0, len(alts))
 555                 for alt := range alts {
 556                         altslice = append(altslice, alt)
 557                 }
 558                 sort.Strings(altslice)
 559                 for i, a := range altslice {
 560                         alts[a] = i + 1
 561                 }
 562                 fmt.Fprintf(out, "%s\t%d\t.\t%s\t%s\t.\t.\t.\tGT", seqname, varslice[0].Position, ref, strings.Join(altslice, ","))
 563                 for i := 0; i < len(varslice); i += 2 {
 564                         v1, v2 := varslice[i], varslice[i+1]
 565                         a1, a2 := alts[v1.New], alts[v2.New]
 566                         if v1.Ref != ref {
 567                                 // variant on allele 0 belongs on a
 568                                 // different output line -- same
 569                                 // chr,pos but different "ref" length
 570                                 a1 = 0
 571                         }
 572                         if v2.Ref != ref {
 573                                 a2 = 0
 574                         }
 575                         fmt.Fprintf(out, "\t%d/%d", a1, a2)
 576                 }
 577                 out.Write([]byte{'\n'})
 578         }
 579 }
 580
 581 func printHGVS(out io.Writer, seqname string, varslice []tvVariant) {
 582         for i := 0; i < len(varslice)/2; i++ {
 583                 if i > 0 {
 584                         out.Write([]byte{'\t'})
 585                 }
 586                 var1, var2 := varslice[i*2], varslice[i*2+1]
 587                 if var1.Variant == var2.Variant {
 588                         if var1.Ref == var1.New {
 589                                 out.Write([]byte{'.'})
 590                         } else {
 591                                 fmt.Fprintf(out, "%s:g.%s", seqname, var1.String())
 592                         }
 593                 } else {
 594                         fmt.Fprintf(out, "%s:g.[%s];[%s]", seqname, var1.String(), var2.String())
 595                 }
 596         }
 597         out.Write([]byte{'\n'})
 598 }
 599
 600 func printHGVSOneHot(out io.Writer, seqname string, varslice []tvVariant) {
 601         vars := map[hgvs.Variant]bool{}
 602         for _, v := range varslice {
 603                 if v.Ref != v.New {
 604                         vars[v.Variant] = true
 605                 }
 606         }
 607
 608         // sort variants to ensure output is deterministic
 609         sorted := make([]hgvs.Variant, 0, len(vars))
 610         for v := range vars {
 611                 sorted = append(sorted, v)
 612         }
 613         sort.Slice(sorted, func(a, b int) bool { return hgvs.Less(sorted[a], sorted[b]) })
 614
 615         for _, v := range sorted {
 616                 fmt.Fprintf(out, "%s.%s", seqname, v.String())
 617                 for i := 0; i < len(varslice); i += 2 {
 618                         if varslice[i].Variant == v || varslice[i+1].Variant == v {
 619                                 out.Write([]byte("\t1"))
 620                         } else {
 621                                 out.Write([]byte("\t0"))
 622                         }
 623                 }
 624                 out.Write([]byte{'\n'})
 625         }
 626 }