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