1 // Copyright (C) The Lightning Authors. All rights reserved.
3 // SPDX-License-Identifier: AGPL-3.0
23 "git.arvados.org/arvados.git/sdk/go/arvados"
24 log "github.com/sirupsen/logrus"
27 type anno2vcf struct {
30 func (cmd *anno2vcf) RunCommand(prog string, args []string, stdin io.Reader, stdout, stderr io.Writer) int {
34 fmt.Fprintf(stderr, "%s\n", err)
37 flags := flag.NewFlagSet("", flag.ContinueOnError)
38 flags.SetOutput(stderr)
39 pprof := flags.String("pprof", "", "serve Go profile data at http://`[addr]:port`")
40 runlocal := flags.Bool("local", false, "run on local host (default: run in an arvados container)")
41 projectUUID := flags.String("project", "", "project `UUID` for output data")
42 priority := flags.Int("priority", 500, "container request priority")
43 inputDir := flags.String("input-dir", "./in", "input `directory`")
44 outputDir := flags.String("output-dir", "./out", "output `directory`")
45 err = flags.Parse(args)
46 if err == flag.ErrHelp {
49 } else if err != nil {
55 log.Println(http.ListenAndServe(*pprof, nil))
60 runner := arvadosContainerRunner{
61 Name: "lightning anno2vcf",
62 Client: arvados.NewClientFromEnv(),
63 ProjectUUID: *projectUUID,
70 err = runner.TranslatePaths(inputDir)
74 runner.Args = []string{"anno2vcf", "-local=true",
76 "-input-dir", *inputDir,
77 "-output-dir", "/mnt/output",
80 output, err = runner.Run()
84 fmt.Fprintln(stdout, output)
88 d, err := open(*inputDir)
94 fis, err := d.Readdir(-1)
100 sort.Slice(fis, func(i, j int) bool { return fis[i].Name() < fis[j].Name() })
110 allcalls := map[string][]*call{}
112 thr := throttle{Max: runtime.GOMAXPROCS(0)}
113 log.Print("reading input files")
114 for _, fi := range fis {
115 if !strings.HasSuffix(fi.Name(), "annotations.csv") {
118 filename := *inputDir + "/" + fi.Name()
119 thr.Go(func() error {
120 log.Printf("reading %s", filename)
121 buf, err := ioutil.ReadFile(filename)
123 return fmt.Errorf("%s: %s", filename, err)
125 lines := bytes.Split(buf, []byte{'\n'})
126 calls := map[string][]*call{}
127 for lineIdx, line := range lines {
131 if lineIdx&0xff == 0 && thr.Err() != nil {
134 fields := bytes.Split(line, []byte{','})
136 return fmt.Errorf("%s line %d: wrong number of fields (%d < %d): %q", fi.Name(), lineIdx+1, len(fields), 8, line)
140 // "=" reference or ""
141 // non-diffable tile variant
144 tile, _ := strconv.ParseInt(string(fields[0]), 10, 64)
145 variant, _ := strconv.ParseInt(string(fields[2]), 10, 64)
146 position, _ := strconv.ParseInt(string(fields[5]), 10, 64)
147 seq := string(fields[4])
148 if calls[seq] == nil {
149 calls[seq] = make([]*call, 0, len(lines)/50)
153 if (len(del) == 0 || len(ins) == 0) && len(fields) >= 9 {
154 // "123,,AA,T" means 123insAA
155 // preceded by T. We record it
156 // here as "122 T TAA" to
157 // avoid writing an empty
158 // "ref" field in our
159 // VCF. Similarly, we record
160 // deletions as "122 TAA T"
161 // rather than "123 AA .".
162 del = append(append(make([]byte, 0, len(fields[8])+len(del)), fields[8]...), del...)
163 ins = append(append(make([]byte, 0, len(fields[8])+len(ins)), fields[8]...), ins...)
164 position -= int64(len(fields[8]))
166 del = append([]byte(nil), del...)
167 ins = append([]byte(nil), ins...)
169 calls[seq] = append(calls[seq], &call{
171 variant: int(variant),
172 position: int(position),
179 for seq, seqcalls := range calls {
180 allcalls[seq] = append(allcalls[seq], seqcalls...)
190 thr = throttle{Max: len(allcalls)}
191 for seq, seqcalls := range allcalls {
192 seq, seqcalls := seq, seqcalls
193 thr.Go(func() error {
194 log.Printf("%s: sorting", seq)
195 sort.Slice(seqcalls, func(i, j int) bool {
196 ii, jj := seqcalls[i], seqcalls[j]
197 if cmp := ii.position - jj.position; cmp != 0 {
200 if cmp := len(ii.deletion) - len(jj.deletion); cmp != 0 {
203 if cmp := bytes.Compare(ii.insertion, jj.insertion); cmp != 0 {
206 if cmp := ii.tile - jj.tile; cmp != 0 {
209 return ii.variant < jj.variant
212 vcfFilename := fmt.Sprintf("%s/annotations.%s.vcf", *outputDir, seq)
213 log.Printf("%s: writing %s", seq, vcfFilename)
215 f, err := os.Create(vcfFilename)
220 bufw := bufio.NewWriterSize(f, 1<<20)
221 _, err = fmt.Fprintf(bufw, `##fileformat=VCFv4.0
222 ##INFO=<ID=TV,Number=.,Type=String,Description="tile-variant">
223 #CHROM POS ID REF ALT QUAL FILTER INFO
228 placeholder := []byte{'.'}
229 for i := 0; i < len(seqcalls); {
232 info := fmt.Sprintf("TV=,%d-%d,", call.tile, call.variant)
233 for i < len(seqcalls) &&
234 call.position == seqcalls[i].position &&
235 len(call.deletion) == len(seqcalls[i].deletion) &&
236 bytes.Equal(call.insertion, seqcalls[i].insertion) {
239 info += fmt.Sprintf("%d-%d,", call.tile, call.variant)
241 deletion := call.deletion
242 if len(deletion) == 0 {
243 deletion = placeholder
245 insertion := call.insertion
246 if len(insertion) == 0 {
247 insertion = placeholder
249 _, err = fmt.Fprintf(bufw, "%s\t%d\t%s\t%s\t%s\t.\t.\t%s\n", seq, call.position, call.hgvsID, deletion, insertion, info)
262 log.Printf("%s: done", seq)