20 "git.arvados.org/arvados.git/sdk/go/arvados"
21 "github.com/arvados/lightning/hgvs"
22 log "github.com/sirupsen/logrus"
25 type annotatecmd struct {
29 tag2tagid map[string]tagID
32 func (cmd *annotatecmd) RunCommand(prog string, args []string, stdin io.Reader, stdout, stderr io.Writer) int {
36 fmt.Fprintf(stderr, "%s\n", err)
39 flags := flag.NewFlagSet("", flag.ContinueOnError)
40 flags.SetOutput(stderr)
41 pprof := flags.String("pprof", "", "serve Go profile data at http://`[addr]:port`")
42 runlocal := flags.Bool("local", false, "run on local host (default: run in an arvados container)")
43 projectUUID := flags.String("project", "", "project `UUID` for output data")
44 priority := flags.Int("priority", 500, "container request priority")
45 inputFilename := flags.String("i", "-", "input `file` (library)")
46 outputFilename := flags.String("o", "-", "output `file`")
47 flags.BoolVar(&cmd.variantHash, "variant-hash", false, "output variant hash instead of index")
48 flags.IntVar(&cmd.maxTileSize, "max-tile-size", 50000, "don't try to make annotations for tiles bigger than given `size`")
49 err = flags.Parse(args)
50 if err == flag.ErrHelp {
53 } else if err != nil {
59 log.Println(http.ListenAndServe(*pprof, nil))
63 if *outputFilename != "-" {
64 err = errors.New("cannot specify output file in container mode: not implemented")
67 runner := arvadosContainerRunner{
68 Name: "lightning annotate",
69 Client: arvados.NewClientFromEnv(),
70 ProjectUUID: *projectUUID,
75 err = runner.TranslatePaths(inputFilename)
79 runner.Args = []string{"annotate", "-local=true", fmt.Sprintf("-variant-hash=%v", cmd.variantHash), "-max-tile-size", strconv.Itoa(cmd.maxTileSize), "-i", *inputFilename, "-o", "/mnt/output/tilevariants.csv"}
81 output, err = runner.Run()
85 fmt.Fprintln(stdout, output+"/tilevariants.csv")
89 var input io.ReadCloser
90 if *inputFilename == "-" {
91 input = ioutil.NopCloser(stdin)
93 input, err = os.Open(*inputFilename)
100 var output io.WriteCloser
101 if *outputFilename == "-" {
102 output = nopCloser{stdout}
104 output, err = os.OpenFile(*outputFilename, os.O_CREATE|os.O_WRONLY, 0666)
110 bufw := bufio.NewWriterSize(output, 4*1024*1024)
112 tilelib := &tileLibrary{
114 retainTileSequences: true,
116 err = tilelib.LoadGob(context.Background(), input, strings.HasSuffix(*inputFilename, ".gz"), nil)
120 err = cmd.exportTileDiffs(bufw, tilelib)
139 func (cmd *annotatecmd) exportTileDiffs(outw io.Writer, tilelib *tileLibrary) error {
140 tagset := tilelib.taglib.Tags()
141 if len(tagset) == 0 {
142 return errors.New("cannot annotate library without tags")
144 taglen := len(tagset[0])
146 for name := range tilelib.refseqs {
147 refs = append(refs, name)
149 cmd.tag2tagid = make(map[string]tagID, len(tagset))
150 for tagid, tagseq := range tagset {
151 cmd.tag2tagid[string(tagseq)] = tagID(tagid)
154 log.Infof("len(refs) %d", len(refs))
156 outch := make(chan string, runtime.NumCPU()*2)
157 var outwg sync.WaitGroup
162 for s := range outch {
163 io.WriteString(outw, s)
169 for _, refcs := range tilelib.refseqs {
173 throttle := &throttle{Max: runtime.NumCPU()*2 + nseqs*2 + 1}
174 defer throttle.Wait()
176 for _, refname := range refs {
178 refcs := tilelib.refseqs[refname]
179 var seqnames []string
180 for seqname := range refcs {
181 seqnames = append(seqnames, seqname)
183 sort.Strings(seqnames)
184 for _, seqname := range seqnames {
187 if throttle.Err() != nil {
191 defer throttle.Release()
192 throttle.Report(cmd.annotateSequence(throttle, outch, tilelib, taglen, refname, seqname, refcs[seqname], len(refs) > 1))
197 return throttle.Err()
200 func (cmd *annotatecmd) annotateSequence(throttle *throttle, outch chan<- string, tilelib *tileLibrary, taglen int, refname, seqname string, reftiles []tileLibRef, refnamecol bool) error {
203 refnamefield = "," + trimFilenameForLabel(refname)
206 // tilestart[123] is the index into refseq
207 // where the tile for tag 123 was placed.
208 tilestart := map[tagID]int{}
209 tileend := map[tagID]int{}
210 for _, libref := range reftiles {
211 if libref.Variant < 1 {
212 return fmt.Errorf("reference %q seq %q uses variant zero at tag %d", refname, seqname, libref.Tag)
214 seq := tilelib.TileVariantSequence(libref)
215 if len(seq) < taglen {
216 return fmt.Errorf("reference %q seq %q uses tile %d variant %d with sequence len %d < taglen %d", refname, seqname, libref.Tag, libref.Variant, len(seq), taglen)
219 if len(refseq) == 0 {
222 tilestart[libref.Tag] = len(refseq) - overlap
223 refseq = append(refseq, seq[overlap:]...)
224 tileend[libref.Tag] = len(refseq)
226 log.Infof("seq %s len(refseq) %d len(tilestart) %d", seqname, len(refseq), len(tilestart))
227 // outtag is tag's index in the subset of tags that aren't
228 // dropped. If there are 10M tags and half are dropped by
229 // dropTiles, tag ranges from 0 to 10M-1 and outtag ranges
232 // IOW, in the matrix built by cgs2array(), {tag} is
233 // represented by columns {outtag}*2 and {outtag}*2+1.
235 for tag, tvs := range tilelib.variant {
236 if len(cmd.dropTiles) > tag && cmd.dropTiles[tag] {
241 // Must shadow outcol var to use safely in goroutine below.
243 refstart, ok := tilestart[tag]
245 // Tag didn't place on this
246 // reference sequence. (It
247 // might place on the same
248 // chromosome in a genome
249 // anyway, but we don't output
250 // the annotations that would
252 // outch <- fmt.Sprintf("%d,%d,-1%s\n", tag, outcol, refnamefield)
255 for variant := 1; variant <= len(tvs); variant++ {
256 variant, hash := tileVariantID(variant), tvs[variant-1]
257 tileseq := tilelib.TileVariantSequence(tileLibRef{Tag: tag, Variant: variant})
258 if len(tileseq) == 0 {
260 } else if len(tileseq) < taglen {
261 return fmt.Errorf("tilevar %d,%d has sequence len %d < taglen %d", tag, variant, len(tileseq), taglen)
264 endtag := string(tileseq[len(tileseq)-taglen:])
265 if endtagid, ok := cmd.tag2tagid[endtag]; !ok {
266 // Tile variant doesn't end on a tag, so it can only place at the end of a chromosome.
267 refpart = refseq[refstart:]
268 log.Warnf("%x tilevar %d,%d endtag not in ref: %s", hash[:13], tag, variant, endtag)
269 } else if refendtagstart, ok := tilestart[endtagid]; !ok {
270 // Ref ends a chromsome with a (possibly very large) variant of this tile, but genomes with this tile don't.
271 // Give up. (TODO: something smarter)
272 log.Debugf("%x not annotating tilevar %d,%d because end tag %d is not in ref", hash[:13], tag, variant, endtagid)
275 // Non-terminal tile vs. non-terminal reference.
276 refpart = refseq[refstart : refendtagstart+taglen]
277 log.Tracef("\n%x tilevar %d,%d endtag %s endtagid %d refendtagstart %d", hash[:13], tag, variant, endtag, endtagid, refendtagstart)
279 if len(refpart) > cmd.maxTileSize {
280 log.Warnf("%x tilevar %d,%d skipping long diff, ref %s seq %s pos %d ref len %d", hash[:13], tag, variant, refname, seqname, refstart, len(refpart))
283 if len(tileseq) > cmd.maxTileSize {
284 log.Warnf("%x tilevar %d,%d skipping long diff, ref %s seq %s variant len %d", hash[:13], tag, variant, refname, seqname, len(tileseq))
287 // log.Printf("\n%x @ refstart %d \n< %s\n> %s\n", tv.Blake2b, refstart, refpart, tileseq)
291 defer throttle.Release()
292 diffs, _ := hgvs.Diff(strings.ToUpper(string(refpart)), strings.ToUpper(string(tileseq)), 0)
293 for _, diff := range diffs {
294 diff.Position += refstart
297 varid = fmt.Sprintf("%x", hash)[:13]
299 varid = fmt.Sprintf("%d", variant)
301 outch <- fmt.Sprintf("%d,%d,%s%s,%s:g.%s\n", tag, outcol, varid, refnamefield, seqname, diff.String())