25 "git.arvados.org/arvados.git/sdk/go/arvados"
26 log "github.com/sirupsen/logrus"
29 type vcf2fasta struct {
34 gvcfRegionsPyData []byte
45 func (cmd *vcf2fasta) RunCommand(prog string, args []string, stdin io.Reader, stdout, stderr io.Writer) int {
49 fmt.Fprintf(stderr, "%s\n", err)
52 flags := flag.NewFlagSet("", flag.ContinueOnError)
53 flags.SetOutput(stderr)
54 flags.StringVar(&cmd.refFile, "ref", "", "reference fasta `file`")
55 flags.StringVar(&cmd.genomeFile, "genome", "", "reference genome `file`")
56 flags.BoolVar(&cmd.mask, "mask", false, "mask uncalled regions (default: output hom ref)")
57 flags.StringVar(&cmd.gvcfRegionsPy, "gvcf-regions.py", "https://raw.githubusercontent.com/lijiayong/gvcf_regions/master/gvcf_regions.py", "source of gvcf_regions.py")
58 flags.StringVar(&cmd.gvcfType, "gvcf-type", "gatk", "gvcf_type argument to gvcf_regions.py: gatk, complete_genomics, freebayes")
59 flags.StringVar(&cmd.projectUUID, "project", "", "project `UUID` for containers and output data")
60 flags.StringVar(&cmd.outputDir, "output-dir", "", "output directory")
61 flags.IntVar(&cmd.vcpus, "vcpus", 0, "number of VCPUs to request for arvados container (default: 2*number of input files, max 32)")
62 flags.BoolVar(&cmd.runLocal, "local", false, "run on local host (default: run in an arvados container)")
63 cmd.batchArgs.Flags(flags)
64 priority := flags.Int("priority", 500, "container request priority")
65 pprof := flags.String("pprof", "", "serve Go profile data at http://`[addr]:port`")
66 err = flags.Parse(args)
67 if err == flag.ErrHelp {
70 } else if err != nil {
72 } else if cmd.refFile == "" {
73 err = errors.New("reference data (-ref) not specified")
75 } else if flags.NArg() == 0 {
83 log.Println(http.ListenAndServe(*pprof, nil))
88 err = cmd.loadRegionsPy()
95 if cmd.outputDir != "" {
96 err = errors.New("cannot specify output dir in non-local mode")
101 infiles, err = listInputFiles(flags.Args())
105 batchsize := (len(infiles) + cmd.batchArgs.batches - 1) / cmd.batchArgs.batches
106 if cmd.vcpus = batchsize * 2; cmd.vcpus > 32 {
110 client := arvados.NewClientFromEnv()
111 runner := arvadosContainerRunner{
112 Name: "lightning vcf2fasta",
114 ProjectUUID: cmd.projectUUID,
115 RAM: 2<<30 + int64(cmd.vcpus)<<28,
118 Mounts: map[string]map[string]interface{}{
119 "/gvcf_regions.py": map[string]interface{}{
121 "content": string(cmd.gvcfRegionsPyData),
125 err = runner.TranslatePaths(&cmd.refFile, &cmd.genomeFile)
129 inputs := flags.Args()
130 for i := range inputs {
131 err = runner.TranslatePaths(&inputs[i])
137 outputs, err = cmd.batchArgs.RunBatches(context.Background(), func(ctx context.Context, batch int) (string, error) {
140 runner.RAM += int64(cmd.vcpus) << 31
142 runner.Args = []string{"vcf2fasta",
144 "-ref", cmd.refFile, fmt.Sprintf("-mask=%v", cmd.mask),
145 "-genome", cmd.genomeFile,
146 "-gvcf-regions.py", "/gvcf_regions.py",
147 "-gvcf-type", cmd.gvcfType,
148 "-output-dir", "/mnt/output",
150 runner.Args = append(runner.Args, cmd.batchArgs.Args(batch)...)
151 runner.Args = append(runner.Args, inputs...)
152 log.Printf("batch %d: %v", batch, runner.Args)
153 return runner.RunContext(ctx)
158 fmt.Fprintln(stdout, strings.Join(outputs, " "))
162 infiles, err := listInputFiles(flags.Args())
166 infiles = cmd.batchArgs.Slice(infiles)
172 todo := make(chan job)
174 for _, infile := range infiles {
175 for phase := 1; phase <= 2; phase++ {
176 todo <- job{vcffile: infile, phase: phase}
182 done := make(chan error, runtime.NumCPU()*2)
183 var wg sync.WaitGroup
184 for i := 0; i < runtime.NumCPU(); i++ {
188 for job := range todo {
190 // a different worker encountered an error
193 err := cmd.vcf2fasta(job.vcffile, job.phase)
195 done <- fmt.Errorf("%s phase %d: %s", job.vcffile, job.phase, err)
213 func maybeInDocker(args, mountfiles []string) []string {
214 if out, err := exec.Command("docker", "image", "ls", "-q", "lightning-runtime").Output(); err != nil || len(out) == 0 {
217 dockerrun := []string{
218 "docker", "run", "--rm", "-i",
221 for _, f := range mountfiles {
222 dockerrun = append(dockerrun, "--volume="+f+":"+f+":ro")
224 dockerrun = append(dockerrun, "lightning-runtime")
225 dockerrun = append(dockerrun, args...)
229 func (cmd *vcf2fasta) vcf2fasta(infile string, phase int) error {
230 ctx, cancel := context.WithCancel(context.Background())
233 _, basename := filepath.Split(infile)
234 outfile := filepath.Join(cmd.outputDir, fmt.Sprintf("%s.%d.fasta.gz", basename, phase))
235 outf, err := os.OpenFile(outfile, os.O_CREATE|os.O_WRONLY, 0777)
237 return fmt.Errorf("error opening output file: %s", err)
240 bufw := bufio.NewWriterSize(outf, 8*1024*1024)
241 gzipw := gzip.NewWriter(bufw)
244 var maskfifo string // filename of mask fifo if we're running bedtools, otherwise ""
246 var wg sync.WaitGroup
247 errs := make(chan error, 2)
249 chrSize := map[string]int{}
251 vcffile, err := os.Open(infile)
255 defer vcffile.Close()
256 var rdr io.Reader = vcffile
257 if strings.HasSuffix(infile, ".gz") {
258 rdr, err = gzip.NewReader(vcffile)
263 contigre := regexp.MustCompile(`([^=,]*)=([^>,]*)`)
264 scanner := bufio.NewScanner(rdr)
266 if s := scanner.Text(); !strings.HasPrefix(s, "##") {
268 } else if !strings.HasPrefix(s, "##contig=<") {
271 kv := map[string]string{}
272 for _, m := range contigre.FindAllStringSubmatch(s[10:], -1) {
275 if kv["ID"] != "" && kv["length"] != "" {
276 chrSize[kv["ID"]], _ = strconv.Atoi(kv["length"])
280 if err = scanner.Err(); err != nil {
281 return fmt.Errorf("error scanning input file %q: %s", infile, err)
283 var regions bytes.Buffer
284 bedargs := []string{"python2", "-", "--gvcf_type", cmd.gvcfType, infile}
285 bed := exec.CommandContext(ctx, bedargs[0], bedargs[1:]...)
286 bed.Stdin = bytes.NewBuffer(cmd.gvcfRegionsPyData)
287 bed.Stdout = ®ions
288 bed.Stderr = cmd.stderr
289 log.Printf("running %v", bed.Args)
291 log.Printf("exited %v", bed.Args)
293 return fmt.Errorf("gvcf_regions: %s", err)
296 if cmd.genomeFile != "" {
297 // Read chromosome sizes from genome file in
298 // case any weren't specified in the VCF
300 genomeFile, err := os.Open(cmd.genomeFile)
302 return fmt.Errorf("error opening genome file %q: %s", cmd.genomeFile, err)
304 scanner := bufio.NewScanner(genomeFile)
308 _, err := fmt.Sscanf(scanner.Text(), "%s\t%d", &chr, &size)
310 return fmt.Errorf("error parsing genome file %q: %s", cmd.genomeFile, err)
312 if chrSize[chr] == 0 {
316 if err = scanner.Err(); err != nil {
317 return fmt.Errorf("error scanning genome file %q: %s", cmd.genomeFile, err)
321 // "bedtools complement" expects the chromosome sizes
322 // ("genome file") to appear in the same order as the
323 // chromosomes in the input vcf, so we need to sort
325 scanner = bufio.NewScanner(bytes.NewBuffer(append([]byte(nil), regions.Bytes()...)))
326 var sortedGenomeFile bytes.Buffer
330 _, err := fmt.Sscanf(scanner.Text(), "%s\t%d", &chr, &size)
332 return fmt.Errorf("error parsing gvcf_regions output: %s", err)
334 if size, ok := chrSize[chr]; ok {
335 fmt.Fprintf(&sortedGenomeFile, "%s\t%d\n", chr, size)
340 // The bcftools --mask argument needs to end in ".bed"
341 // in order to be parsed as a BED file, so we need to
342 // use a named pipe instead of stdin.
343 tempdir, err := ioutil.TempDir("", "")
345 return fmt.Errorf("TempDir: %s", err)
347 defer os.RemoveAll(tempdir)
348 maskfifo = filepath.Join(tempdir, "fifo.bed")
349 err = syscall.Mkfifo(maskfifo, 0600)
351 return fmt.Errorf("mkfifo: %s", err)
354 // bedtools complement can't seem to read from a pipe
355 // reliably -- "Error: line number 1 of file
356 // /dev/stdin has 1 fields, but 3 were expected." --
357 // so we stage to a temp file.
358 regionsFile := filepath.Join(tempdir, "gvcf_regions.bed")
359 err = ioutil.WriteFile(regionsFile, regions.Bytes(), 0644)
368 maskfifow, err := os.OpenFile(maskfifo, os.O_WRONLY, 0)
373 defer maskfifow.Close()
375 bedcompargs := []string{"bedtools", "complement", "-i", regionsFile, "-g", "/dev/stdin"}
376 bedcompargs = maybeInDocker(bedcompargs, []string{cmd.genomeFile})
377 bedcomp := exec.CommandContext(ctx, bedcompargs[0], bedcompargs[1:]...)
378 bedcomp.Stdin = &sortedGenomeFile
379 bedcomp.Stdout = maskfifow
380 bedcomp.Stderr = cmd.stderr
381 log.Printf("running %v", bedcomp.Args)
383 log.Printf("exited %v", bedcomp.Args)
385 errs <- fmt.Errorf("bedtools complement: %s", err)
388 err = maskfifow.Close()
399 consargs := []string{"bcftools", "consensus", "--fasta-ref", cmd.refFile, "-H", fmt.Sprint(phase)}
401 consargs = append(consargs, "--mask", maskfifo)
403 consargs = append(consargs, infile)
404 indexsuffix := ".tbi"
405 if _, err := os.Stat(infile + ".csi"); err == nil {
408 mounts := []string{infile, infile + indexsuffix, cmd.refFile}
410 mounts = append(mounts, maskfifo)
412 consargs = maybeInDocker(consargs, mounts)
414 consensus := exec.CommandContext(ctx, consargs[0], consargs[1:]...)
415 consensus.Stderr = os.Stderr
416 consensus.Stdout = gzipw
417 consensus.Stderr = cmd.stderr
418 log.Printf("running %v", consensus.Args)
419 err = consensus.Run()
421 errs <- fmt.Errorf("bcftools consensus: %s", err)
424 log.Printf("exited %v", consensus.Args)
443 for err := range errs {
453 func (cmd *vcf2fasta) loadRegionsPy() error {
454 if strings.HasPrefix(cmd.gvcfRegionsPy, "http") {
455 resp, err := http.Get(cmd.gvcfRegionsPy)
459 if resp.StatusCode != http.StatusOK {
460 return fmt.Errorf("get %q: http status %d", cmd.gvcfRegionsPy, resp.StatusCode)
462 buf, err := ioutil.ReadAll(resp.Body)
464 return fmt.Errorf("get %q: read body: %s", cmd.gvcfRegionsPy, err)
466 cmd.gvcfRegionsPyData = buf
469 buf, err := ioutil.ReadFile(cmd.gvcfRegionsPy)
473 cmd.gvcfRegionsPyData = buf