diff.go

   1 package main
   2
   3 import (
   4         "bufio"
   5         "bytes"
   6         "flag"
   7         "fmt"
   8         "io"
   9         "os"
  10
  11         "github.com/arvados/lightning/hgvs"
  12 )
  13
  14 type diffFasta struct{}
  15
  16 func (cmd *diffFasta) RunCommand(prog string, args []string, stdin io.Reader, stdout, stderr io.Writer) int {
  17         var err error
  18         defer func() {
  19                 if err != nil {
  20                         fmt.Fprintf(stderr, "%s\n", err)
  21                 }
  22         }()
  23         flags := flag.NewFlagSet("", flag.ContinueOnError)
  24         flags.SetOutput(stderr)
  25         offset := flags.Int("offset", 0, "coordinate offset")
  26         sequence := flags.String("sequence", "chr1", "sequence label")
  27         err = flags.Parse(args)
  28         if err == flag.ErrHelp {
  29                 err = nil
  30                 return 0
  31         } else if err != nil {
  32                 return 2
  33         }
  34         if len(flags.Args()) != 2 {
  35                 err = fmt.Errorf("usage: %s [options] a.fasta b.fasta", prog)
  36                 return 2
  37         }
  38
  39         var fasta [2][]byte
  40         errs := make(chan error, 2)
  41         for idx, fnm := range flags.Args() {
  42                 idx, fnm := idx, fnm
  43                 go func() {
  44                         f, err := os.Open(fnm)
  45                         if err != nil {
  46                                 errs <- err
  47                                 return
  48                         }
  49                         defer f.Close()
  50                         scanner := bufio.NewScanner(f)
  51                         for scanner.Scan() {
  52                                 buf := scanner.Bytes()
  53                                 if len(buf) > 0 && buf[0] != '>' {
  54                                         fasta[idx] = append(fasta[idx], bytes.ToUpper(buf)...)
  55                                 }
  56                         }
  57                         errs <- scanner.Err()
  58                 }()
  59         }
  60         for range flags.Args() {
  61                 if err = <-errs; err != nil {
  62                         return 1
  63                 }
  64         }
  65
  66         variants := hgvs.Diff(string(fasta[0]), string(fasta[1]))
  67         if *offset != 0 {
  68                 for i := range variants {
  69                         variants[i].Position += *offset
  70                 }
  71         }
  72         for _, v := range variants {
  73                 fmt.Fprintf(stdout, "%s:g.%s\t%s\t%d\t%s\t%s\n", *sequence, v.String(), *sequence, v.Position, v.Ref, v.New)
  74         }
  75         return 0
  76 }