24 "git.arvados.org/arvados.git/sdk/go/arvados"
25 log "github.com/sirupsen/logrus"
28 type importer struct {
40 func (cmd *importer) RunCommand(prog string, args []string, stdin io.Reader, stdout, stderr io.Writer) int {
44 fmt.Fprintf(stderr, "%s\n", err)
47 flags := flag.NewFlagSet("", flag.ContinueOnError)
48 flags.SetOutput(stderr)
49 flags.StringVar(&cmd.tagLibraryFile, "tag-library", "", "tag library fasta `file`")
50 flags.StringVar(&cmd.refFile, "ref", "", "reference fasta `file`")
51 flags.StringVar(&cmd.outputFile, "o", "-", "output `file`")
52 flags.StringVar(&cmd.projectUUID, "project", "", "project `UUID` for output data")
53 flags.BoolVar(&cmd.runLocal, "local", false, "run on local host (default: run in an arvados container)")
54 flags.BoolVar(&cmd.skipOOO, "skip-ooo", false, "skip out-of-order tags")
55 flags.BoolVar(&cmd.outputTiles, "output-tiles", false, "include tile variant sequences in output file")
56 flags.BoolVar(&cmd.includeNoCalls, "include-no-calls", false, "treat tiles with no-calls as regular tiles")
57 priority := flags.Int("priority", 500, "container request priority")
58 pprof := flags.String("pprof", "", "serve Go profile data at http://`[addr]:port`")
59 loglevel := flags.String("loglevel", "info", "logging threshold (trace, debug, info, warn, error, fatal, or panic)")
60 err = flags.Parse(args)
61 if err == flag.ErrHelp {
64 } else if err != nil {
66 } else if cmd.tagLibraryFile == "" {
67 fmt.Fprintln(os.Stderr, "cannot import without -tag-library argument")
69 } else if flags.NArg() == 0 {
76 log.Println(http.ListenAndServe(*pprof, nil))
80 lvl, err := log.ParseLevel(*loglevel)
87 runner := arvadosContainerRunner{
88 Name: "lightning import",
89 Client: arvados.NewClientFromEnv(),
90 ProjectUUID: cmd.projectUUID,
95 err = runner.TranslatePaths(&cmd.tagLibraryFile, &cmd.refFile, &cmd.outputFile)
99 inputs := flags.Args()
100 for i := range inputs {
101 err = runner.TranslatePaths(&inputs[i])
106 if cmd.outputFile == "-" {
107 cmd.outputFile = "/mnt/output/library.gob"
109 // Not yet implemented, but this should write
110 // the collection to an existing collection,
111 // possibly even an in-place update.
112 err = errors.New("cannot specify output file in container mode: not implemented")
115 runner.Args = append([]string{"import",
117 "-loglevel=" + *loglevel,
118 fmt.Sprintf("-skip-ooo=%v", cmd.skipOOO),
119 fmt.Sprintf("-output-tiles=%v", cmd.outputTiles),
120 fmt.Sprintf("-include-no-calls=%v", cmd.includeNoCalls),
121 "-tag-library", cmd.tagLibraryFile,
123 "-o", cmd.outputFile,
126 output, err = runner.Run()
130 fmt.Fprintln(stdout, output+"/library.gob")
134 infiles, err := listInputFiles(flags.Args())
139 taglib, err := cmd.loadTagLibrary()
144 var output io.WriteCloser
145 if cmd.outputFile == "-" {
146 output = nopCloser{stdout}
148 output, err = os.OpenFile(cmd.outputFile, os.O_CREATE|os.O_WRONLY, 0777)
154 bufw := bufio.NewWriter(output)
155 cmd.encoder = gob.NewEncoder(bufw)
157 tilelib := &tileLibrary{taglib: taglib, includeNoCalls: cmd.includeNoCalls, skipOOO: cmd.skipOOO}
159 tilelib.encoder = cmd.encoder
162 for range time.Tick(10 * time.Minute) {
163 log.Printf("tilelib.Len() == %d", tilelib.Len())
167 err = cmd.tileInputs(tilelib, infiles)
182 func (cmd *importer) tileFasta(tilelib *tileLibrary, infile string) (tileSeq, error) {
183 var input io.ReadCloser
184 input, err := os.Open(infile)
189 if strings.HasSuffix(infile, ".gz") {
190 input, err = gzip.NewReader(input)
196 return tilelib.TileFasta(infile, input)
199 func (cmd *importer) loadTagLibrary() (*tagLibrary, error) {
200 log.Printf("tag library %s load starting", cmd.tagLibraryFile)
201 f, err := os.Open(cmd.tagLibraryFile)
206 var rdr io.ReadCloser = f
207 if strings.HasSuffix(cmd.tagLibraryFile, ".gz") {
208 rdr, err = gzip.NewReader(f)
210 return nil, fmt.Errorf("%s: gzip: %s", cmd.tagLibraryFile, err)
214 var taglib tagLibrary
215 err = taglib.Load(rdr)
219 if taglib.Len() < 1 {
220 return nil, fmt.Errorf("cannot tile: tag library is empty")
222 log.Printf("tag library %s load done", cmd.tagLibraryFile)
226 func listInputFiles(paths []string) (files []string, err error) {
227 for _, path := range paths {
228 if fi, err := os.Stat(path); err != nil {
229 return nil, fmt.Errorf("%s: stat failed: %s", path, err)
230 } else if !fi.IsDir() {
231 if !strings.HasSuffix(path, ".2.fasta") || strings.HasSuffix(path, ".2.fasta.gz") {
232 files = append(files, path)
236 d, err := os.Open(path)
238 return nil, fmt.Errorf("%s: open failed: %s", path, err)
241 names, err := d.Readdirnames(0)
243 return nil, fmt.Errorf("%s: readdir failed: %s", path, err)
246 for _, name := range names {
247 if strings.HasSuffix(name, ".vcf") || strings.HasSuffix(name, ".vcf.gz") {
248 files = append(files, filepath.Join(path, name))
249 } else if strings.HasSuffix(name, ".1.fasta") || strings.HasSuffix(name, ".1.fasta.gz") {
250 files = append(files, filepath.Join(path, name))
255 for _, file := range files {
256 if strings.HasSuffix(file, ".1.fasta") || strings.HasSuffix(file, ".1.fasta.gz") {
258 } else if _, err := os.Stat(file + ".csi"); err == nil {
260 } else if _, err = os.Stat(file + ".tbi"); err == nil {
263 return nil, fmt.Errorf("%s: cannot read without .tbi or .csi index file", file)
269 func (cmd *importer) tileInputs(tilelib *tileLibrary, infiles []string) error {
270 starttime := time.Now()
271 errs := make(chan error, 1)
272 todo := make(chan func() error, len(infiles)*2)
273 var encodeJobs sync.WaitGroup
274 for _, infile := range infiles {
276 var phases sync.WaitGroup
278 variants := make([][]tileVariantID, 2)
279 if strings.HasSuffix(infile, ".1.fasta") || strings.HasSuffix(infile, ".1.fasta.gz") {
280 todo <- func() error {
282 log.Printf("%s starting", infile)
283 defer log.Printf("%s done", infile)
284 tseqs, err := cmd.tileFasta(tilelib, infile)
285 var kept, dropped int
286 variants[0], kept, dropped = tseqs.Variants()
287 log.Printf("%s found %d unique tags plus %d repeats", infile, kept, dropped)
290 infile2 := regexp.MustCompile(`\.1\.fasta(\.gz)?$`).ReplaceAllString(infile, `.2.fasta$1`)
291 todo <- func() error {
293 log.Printf("%s starting", infile2)
294 defer log.Printf("%s done", infile2)
295 tseqs, err := cmd.tileFasta(tilelib, infile2)
296 var kept, dropped int
297 variants[1], kept, dropped = tseqs.Variants()
298 log.Printf("%s found %d unique tags plus %d repeats", infile, kept, dropped)
302 for phase := 0; phase < 2; phase++ {
304 todo <- func() error {
306 log.Printf("%s phase %d starting", infile, phase+1)
307 defer log.Printf("%s phase %d done", infile, phase+1)
308 tseqs, err := cmd.tileGVCF(tilelib, infile, phase)
309 var kept, dropped int
310 variants[phase], kept, dropped = tseqs.Variants()
311 log.Printf("%s phase %d found %d unique tags plus %d repeats", infile, phase+1, kept, dropped)
318 defer encodeJobs.Done()
323 ntags := len(variants[0])
324 if ntags < len(variants[1]) {
325 ntags = len(variants[1])
327 flat := make([]tileVariantID, ntags*2)
328 for i := 0; i < ntags; i++ {
329 for hap := 0; hap < 2; hap++ {
330 if i < len(variants[hap]) {
331 flat[i*2+hap] = variants[hap][i]
335 err := cmd.encoder.Encode(LibraryEntry{
336 CompactGenomes: []CompactGenome{{Name: infile, Variants: flat}},
347 var tileJobs sync.WaitGroup
349 for i := 0; i < runtime.NumCPU()*9/8+1; i++ {
351 atomic.AddInt64(&running, 1)
353 defer tileJobs.Done()
354 defer atomic.AddInt64(&running, -1)
355 for fn := range todo {
366 remain := len(todo) + int(atomic.LoadInt64(&running)) - 1
367 ttl := time.Now().Sub(starttime) * time.Duration(remain) / time.Duration(cap(todo)-remain)
368 eta := time.Now().Add(ttl)
369 log.Printf("progress %d/%d, eta %v (%v)", cap(todo)-remain, cap(todo), eta, ttl)
379 func (cmd *importer) tileGVCF(tilelib *tileLibrary, infile string, phase int) (tileseq tileSeq, err error) {
380 if cmd.refFile == "" {
381 err = errors.New("cannot import vcf: reference data (-ref) not specified")
384 args := []string{"bcftools", "consensus", "--fasta-ref", cmd.refFile, "-H", fmt.Sprint(phase + 1), infile}
385 indexsuffix := ".tbi"
386 if _, err := os.Stat(infile + ".csi"); err == nil {
389 if out, err := exec.Command("docker", "image", "ls", "-q", "lightning-runtime").Output(); err == nil && len(out) > 0 {
390 args = append([]string{
391 "docker", "run", "--rm",
393 "--volume=" + infile + ":" + infile + ":ro",
394 "--volume=" + infile + indexsuffix + ":" + infile + indexsuffix + ":ro",
395 "--volume=" + cmd.refFile + ":" + cmd.refFile + ":ro",
399 consensus := exec.Command(args[0], args[1:]...)
400 consensus.Stderr = os.Stderr
401 stdout, err := consensus.StdoutPipe()
406 err = consensus.Start()
410 defer consensus.Wait()
411 tileseq, err = tilelib.TileFasta(fmt.Sprintf("%s phase %d", infile, phase+1), stdout)
419 err = consensus.Wait()
421 err = fmt.Errorf("%s phase %d: bcftools: %s", infile, phase, err)