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