e6b86afd81ca22038c025c669e39a6fd3a0c3b88
[lightning.git] / import.go
1 package main
2
3 import (
4         "bufio"
5         "compress/gzip"
6         "encoding/gob"
7         "errors"
8         "flag"
9         "fmt"
10         "io"
11         "net/http"
12         _ "net/http/pprof"
13         "os"
14         "os/exec"
15         "path/filepath"
16         "regexp"
17         "runtime"
18         "sort"
19         "strings"
20         "sync"
21         "sync/atomic"
22         "time"
23
24         "git.arvados.org/arvados.git/sdk/go/arvados"
25         log "github.com/sirupsen/logrus"
26 )
27
28 type importer struct {
29         tagLibraryFile string
30         refFile        string
31         outputFile     string
32         projectUUID    string
33         runLocal       bool
34         encoder        *gob.Encoder
35 }
36
37 func (cmd *importer) RunCommand(prog string, args []string, stdin io.Reader, stdout, stderr io.Writer) int {
38         var err error
39         defer func() {
40                 if err != nil {
41                         fmt.Fprintf(stderr, "%s\n", err)
42                 }
43         }()
44         flags := flag.NewFlagSet("", flag.ContinueOnError)
45         flags.SetOutput(stderr)
46         flags.StringVar(&cmd.tagLibraryFile, "tag-library", "", "tag library fasta `file`")
47         flags.StringVar(&cmd.refFile, "ref", "", "reference fasta `file`")
48         flags.StringVar(&cmd.outputFile, "o", "-", "output `file`")
49         flags.StringVar(&cmd.projectUUID, "project", "", "project `UUID` for output data")
50         flags.BoolVar(&cmd.runLocal, "local", false, "run on local host (default: run in an arvados container)")
51         priority := flags.Int("priority", 500, "container request priority")
52         pprof := flags.String("pprof", "", "serve Go profile data at http://`[addr]:port`")
53         err = flags.Parse(args)
54         if err == flag.ErrHelp {
55                 err = nil
56                 return 0
57         } else if err != nil {
58                 return 2
59         } else if cmd.tagLibraryFile == "" {
60                 fmt.Fprintln(os.Stderr, "cannot import without -tag-library argument")
61                 return 2
62         } else if flags.NArg() == 0 {
63                 flags.Usage()
64                 return 2
65         }
66
67         if *pprof != "" {
68                 go func() {
69                         log.Println(http.ListenAndServe(*pprof, nil))
70                 }()
71         }
72
73         if !cmd.runLocal {
74                 runner := arvadosContainerRunner{
75                         Name:        "lightning import",
76                         Client:      arvados.NewClientFromEnv(),
77                         ProjectUUID: cmd.projectUUID,
78                         RAM:         30000000000,
79                         VCPUs:       16,
80                         Priority:    *priority,
81                 }
82                 err = runner.TranslatePaths(&cmd.tagLibraryFile, &cmd.refFile, &cmd.outputFile)
83                 if err != nil {
84                         return 1
85                 }
86                 inputs := flags.Args()
87                 for i := range inputs {
88                         err = runner.TranslatePaths(&inputs[i])
89                         if err != nil {
90                                 return 1
91                         }
92                 }
93                 if cmd.outputFile == "-" {
94                         cmd.outputFile = "/mnt/output/library.gob"
95                 } else {
96                         // Not yet implemented, but this should write
97                         // the collection to an existing collection,
98                         // possibly even an in-place update.
99                         err = errors.New("cannot specify output file in container mode: not implemented")
100                         return 1
101                 }
102                 runner.Args = append([]string{"import", "-local=true", "-tag-library", cmd.tagLibraryFile, "-ref", cmd.refFile, "-o", cmd.outputFile}, inputs...)
103                 var output string
104                 output, err = runner.Run()
105                 if err != nil {
106                         return 1
107                 }
108                 fmt.Fprintln(stdout, output+"/library.gob")
109                 return 0
110         }
111
112         infiles, err := listInputFiles(flags.Args())
113         if err != nil {
114                 return 1
115         }
116
117         tilelib, err := cmd.loadTileLibrary()
118         if err != nil {
119                 return 1
120         }
121         go func() {
122                 for range time.Tick(10 * time.Minute) {
123                         log.Printf("tilelib.Len() == %d", tilelib.Len())
124                 }
125         }()
126
127         var output io.WriteCloser
128         if cmd.outputFile == "-" {
129                 output = nopCloser{stdout}
130         } else {
131                 output, err = os.OpenFile(cmd.outputFile, os.O_CREATE|os.O_WRONLY, 0777)
132                 if err != nil {
133                         return 1
134                 }
135                 defer output.Close()
136         }
137         bufw := bufio.NewWriter(output)
138         cmd.encoder = gob.NewEncoder(bufw)
139
140         err = cmd.tileInputs(tilelib, infiles)
141         if err != nil {
142                 return 1
143         }
144         err = bufw.Flush()
145         if err != nil {
146                 return 1
147         }
148         err = output.Close()
149         if err != nil {
150                 return 1
151         }
152         return 0
153 }
154
155 func (cmd *importer) tileFasta(tilelib *tileLibrary, infile string) (tileSeq, error) {
156         var input io.ReadCloser
157         input, err := os.Open(infile)
158         if err != nil {
159                 return nil, err
160         }
161         defer input.Close()
162         if strings.HasSuffix(infile, ".gz") {
163                 input, err = gzip.NewReader(input)
164                 if err != nil {
165                         return nil, err
166                 }
167                 defer input.Close()
168         }
169         return tilelib.TileFasta(infile, input)
170 }
171
172 func (cmd *importer) loadTileLibrary() (*tileLibrary, error) {
173         log.Printf("tag library %s load starting", cmd.tagLibraryFile)
174         f, err := os.Open(cmd.tagLibraryFile)
175         if err != nil {
176                 return nil, err
177         }
178         defer f.Close()
179         var rdr io.ReadCloser = f
180         if strings.HasSuffix(cmd.tagLibraryFile, ".gz") {
181                 rdr, err = gzip.NewReader(f)
182                 if err != nil {
183                         return nil, fmt.Errorf("%s: gzip: %s", cmd.tagLibraryFile, err)
184                 }
185                 defer rdr.Close()
186         }
187         var taglib tagLibrary
188         err = taglib.Load(rdr)
189         if err != nil {
190                 return nil, err
191         }
192         if taglib.Len() < 1 {
193                 return nil, fmt.Errorf("cannot tile: tag library is empty")
194         }
195         log.Printf("tag library %s load done", cmd.tagLibraryFile)
196         return &tileLibrary{taglib: &taglib}, nil
197 }
198
199 func listInputFiles(paths []string) (files []string, err error) {
200         for _, path := range paths {
201                 if fi, err := os.Stat(path); err != nil {
202                         return nil, fmt.Errorf("%s: stat failed: %s", path, err)
203                 } else if !fi.IsDir() {
204                         if !strings.HasSuffix(path, ".2.fasta") || strings.HasSuffix(path, ".2.fasta.gz") {
205                                 files = append(files, path)
206                         }
207                         continue
208                 }
209                 d, err := os.Open(path)
210                 if err != nil {
211                         return nil, fmt.Errorf("%s: open failed: %s", path, err)
212                 }
213                 defer d.Close()
214                 names, err := d.Readdirnames(0)
215                 if err != nil {
216                         return nil, fmt.Errorf("%s: readdir failed: %s", path, err)
217                 }
218                 sort.Strings(names)
219                 for _, name := range names {
220                         if strings.HasSuffix(name, ".vcf") || strings.HasSuffix(name, ".vcf.gz") {
221                                 files = append(files, filepath.Join(path, name))
222                         } else if strings.HasSuffix(name, ".1.fasta") || strings.HasSuffix(name, ".1.fasta.gz") {
223                                 files = append(files, filepath.Join(path, name))
224                         }
225                 }
226                 d.Close()
227         }
228         for _, file := range files {
229                 if strings.HasSuffix(file, ".1.fasta") || strings.HasSuffix(file, ".1.fasta.gz") {
230                         continue
231                 } else if _, err := os.Stat(file + ".csi"); err == nil {
232                         continue
233                 } else if _, err = os.Stat(file + ".tbi"); err == nil {
234                         continue
235                 } else {
236                         return nil, fmt.Errorf("%s: cannot read without .tbi or .csi index file", file)
237                 }
238         }
239         return
240 }
241
242 func (cmd *importer) tileInputs(tilelib *tileLibrary, infiles []string) error {
243         starttime := time.Now()
244         errs := make(chan error, 1)
245         todo := make(chan func() error, len(infiles)*2)
246         var encodeJobs sync.WaitGroup
247         for _, infile := range infiles {
248                 infile := infile
249                 var phases sync.WaitGroup
250                 phases.Add(2)
251                 variants := make([][]tileVariantID, 2)
252                 if strings.HasSuffix(infile, ".1.fasta") || strings.HasSuffix(infile, ".1.fasta.gz") {
253                         todo <- func() error {
254                                 defer phases.Done()
255                                 log.Printf("%s starting", infile)
256                                 defer log.Printf("%s done", infile)
257                                 tseqs, err := cmd.tileFasta(tilelib, infile)
258                                 variants[0] = tseqs.Variants()
259                                 return err
260                         }
261                         infile2 := regexp.MustCompile(`\.1\.fasta(\.gz)?$`).ReplaceAllString(infile, `.2.fasta$1`)
262                         todo <- func() error {
263                                 defer phases.Done()
264                                 log.Printf("%s starting", infile2)
265                                 defer log.Printf("%s done", infile2)
266                                 tseqs, err := cmd.tileFasta(tilelib, infile2)
267                                 variants[1] = tseqs.Variants()
268                                 return err
269                         }
270                 } else {
271                         for phase := 0; phase < 2; phase++ {
272                                 phase := phase
273                                 todo <- func() error {
274                                         defer phases.Done()
275                                         log.Printf("%s phase %d starting", infile, phase+1)
276                                         defer log.Printf("%s phase %d done", infile, phase+1)
277                                         tseqs, err := cmd.tileGVCF(tilelib, infile, phase)
278                                         variants[phase] = tseqs.Variants()
279                                         return err
280                                 }
281                         }
282                 }
283                 encodeJobs.Add(1)
284                 go func() {
285                         defer encodeJobs.Done()
286                         phases.Wait()
287                         if len(errs) > 0 {
288                                 return
289                         }
290                         ntags := len(variants[0])
291                         if ntags < len(variants[1]) {
292                                 ntags = len(variants[1])
293                         }
294                         flat := make([]tileVariantID, ntags*2)
295                         for i := 0; i < ntags; i++ {
296                                 flat[i*2] = variants[0][i]
297                                 flat[i*2+1] = variants[1][i]
298                         }
299                         err := cmd.encoder.Encode(LibraryEntry{
300                                 CompactGenomes: []CompactGenome{{Name: infile, Variants: flat}},
301                         })
302                         if err != nil {
303                                 select {
304                                 case errs <- err:
305                                 default:
306                                 }
307                         }
308                 }()
309         }
310         go close(todo)
311         var tileJobs sync.WaitGroup
312         var running int64
313         for i := 0; i < runtime.NumCPU()*9/8+1; i++ {
314                 tileJobs.Add(1)
315                 atomic.AddInt64(&running, 1)
316                 go func() {
317                         defer tileJobs.Done()
318                         defer atomic.AddInt64(&running, -1)
319                         for fn := range todo {
320                                 if len(errs) > 0 {
321                                         return
322                                 }
323                                 err := fn()
324                                 if err != nil {
325                                         select {
326                                         case errs <- err:
327                                         default:
328                                         }
329                                 }
330                                 remain := len(todo) + int(atomic.LoadInt64(&running)) - 1
331                                 ttl := time.Now().Sub(starttime) * time.Duration(remain) / time.Duration(cap(todo)-remain)
332                                 eta := time.Now().Add(ttl)
333                                 log.Printf("progress %d/%d, eta %v (%v)", cap(todo)-remain, cap(todo), eta, ttl)
334                         }
335                 }()
336         }
337         tileJobs.Wait()
338         encodeJobs.Wait()
339         go close(errs)
340         return <-errs
341 }
342
343 func (cmd *importer) tileGVCF(tilelib *tileLibrary, infile string, phase int) (tileseq tileSeq, err error) {
344         if cmd.refFile == "" {
345                 err = errors.New("cannot import vcf: reference data (-ref) not specified")
346                 return
347         }
348         args := []string{"bcftools", "consensus", "--fasta-ref", cmd.refFile, "-H", fmt.Sprint(phase + 1), infile}
349         indexsuffix := ".tbi"
350         if _, err := os.Stat(infile + ".csi"); err == nil {
351                 indexsuffix = ".csi"
352         }
353         if out, err := exec.Command("docker", "image", "ls", "-q", "lightning-runtime").Output(); err == nil && len(out) > 0 {
354                 args = append([]string{
355                         "docker", "run", "--rm",
356                         "--log-driver=none",
357                         "--volume=" + infile + ":" + infile + ":ro",
358                         "--volume=" + infile + indexsuffix + ":" + infile + indexsuffix + ":ro",
359                         "--volume=" + cmd.refFile + ":" + cmd.refFile + ":ro",
360                         "lightning-runtime",
361                 }, args...)
362         }
363         consensus := exec.Command(args[0], args[1:]...)
364         consensus.Stderr = os.Stderr
365         stdout, err := consensus.StdoutPipe()
366         defer stdout.Close()
367         if err != nil {
368                 return
369         }
370         err = consensus.Start()
371         if err != nil {
372                 return
373         }
374         defer consensus.Wait()
375         tileseq, err = tilelib.TileFasta(fmt.Sprintf("%s phase %d", infile, phase+1), stdout)
376         if err != nil {
377                 return
378         }
379         err = stdout.Close()
380         if err != nil {
381                 return
382         }
383         err = consensus.Wait()
384         if err != nil {
385                 err = fmt.Errorf("%s phase %d: bcftools: %s", infile, phase, err)
386                 return
387         }
388         return
389 }