"bufio"
"compress/gzip"
"encoding/gob"
+ "errors"
"flag"
"fmt"
"io"
- "log"
"net/http"
_ "net/http/pprof"
"os"
"sync"
"sync/atomic"
"time"
+
+ "git.arvados.org/arvados.git/sdk/go/arvados"
+ log "github.com/sirupsen/logrus"
)
type importer struct {
tagLibraryFile string
refFile string
+ outputFile string
+ projectUUID string
+ runLocal bool
+ skipOOO bool
+ outputTiles bool
+ includeNoCalls bool
encoder *gob.Encoder
}
flags.SetOutput(stderr)
flags.StringVar(&cmd.tagLibraryFile, "tag-library", "", "tag library fasta `file`")
flags.StringVar(&cmd.refFile, "ref", "", "reference fasta `file`")
+ flags.StringVar(&cmd.outputFile, "o", "-", "output `file`")
+ flags.StringVar(&cmd.projectUUID, "project", "", "project `UUID` for output data")
+ flags.BoolVar(&cmd.runLocal, "local", false, "run on local host (default: run in an arvados container)")
+ flags.BoolVar(&cmd.skipOOO, "skip-ooo", false, "skip out-of-order tags")
+ flags.BoolVar(&cmd.outputTiles, "output-tiles", false, "include tile variant sequences in output file")
+ flags.BoolVar(&cmd.includeNoCalls, "include-no-calls", false, "treat tiles with no-calls as regular tiles")
+ priority := flags.Int("priority", 500, "container request priority")
pprof := flags.String("pprof", "", "serve Go profile data at http://`[addr]:port`")
+ loglevel := flags.String("loglevel", "info", "logging threshold (trace, debug, info, warn, error, fatal, or panic)")
err = flags.Parse(args)
if err == flag.ErrHelp {
err = nil
return 0
} else if err != nil {
return 2
- } else if cmd.refFile == "" || cmd.tagLibraryFile == "" {
- fmt.Fprintln(os.Stderr, "cannot run without -tag-library and -ref arguments")
+ } else if cmd.tagLibraryFile == "" {
+ fmt.Fprintln(os.Stderr, "cannot import without -tag-library argument")
return 2
} else if flags.NArg() == 0 {
flags.Usage()
}()
}
+ lvl, err := log.ParseLevel(*loglevel)
+ if err != nil {
+ return 2
+ }
+ log.SetLevel(lvl)
+
+ if !cmd.runLocal {
+ runner := arvadosContainerRunner{
+ Name: "lightning import",
+ Client: arvados.NewClientFromEnv(),
+ ProjectUUID: cmd.projectUUID,
+ RAM: 60000000000,
+ VCPUs: 16,
+ Priority: *priority,
+ }
+ err = runner.TranslatePaths(&cmd.tagLibraryFile, &cmd.refFile, &cmd.outputFile)
+ if err != nil {
+ return 1
+ }
+ inputs := flags.Args()
+ for i := range inputs {
+ err = runner.TranslatePaths(&inputs[i])
+ if err != nil {
+ return 1
+ }
+ }
+ if cmd.outputFile == "-" {
+ cmd.outputFile = "/mnt/output/library.gob"
+ } else {
+ // Not yet implemented, but this should write
+ // the collection to an existing collection,
+ // possibly even an in-place update.
+ err = errors.New("cannot specify output file in container mode: not implemented")
+ return 1
+ }
+ runner.Args = append([]string{"import",
+ "-local=true",
+ "-loglevel=" + *loglevel,
+ fmt.Sprintf("-skip-ooo=%v", cmd.skipOOO),
+ fmt.Sprintf("-output-tiles=%v", cmd.outputTiles),
+ fmt.Sprintf("-include-no-calls=%v", cmd.includeNoCalls),
+ "-tag-library", cmd.tagLibraryFile,
+ "-ref", cmd.refFile,
+ "-o", cmd.outputFile,
+ }, inputs...)
+ var output string
+ output, err = runner.Run()
+ if err != nil {
+ return 1
+ }
+ fmt.Fprintln(stdout, output+"/library.gob")
+ return 0
+ }
+
infiles, err := listInputFiles(flags.Args())
if err != nil {
return 1
}
- tilelib, err := cmd.loadTileLibrary()
+ taglib, err := cmd.loadTagLibrary()
if err != nil {
return 1
}
+
+ var output io.WriteCloser
+ if cmd.outputFile == "-" {
+ output = nopCloser{stdout}
+ } else {
+ output, err = os.OpenFile(cmd.outputFile, os.O_CREATE|os.O_WRONLY, 0777)
+ if err != nil {
+ return 1
+ }
+ defer output.Close()
+ }
+ bufw := bufio.NewWriter(output)
+ cmd.encoder = gob.NewEncoder(bufw)
+
+ tilelib := &tileLibrary{taglib: taglib, includeNoCalls: cmd.includeNoCalls, skipOOO: cmd.skipOOO}
+ if cmd.outputTiles {
+ cmd.encoder.Encode(LibraryEntry{TagSet: taglib.Tags()})
+ tilelib.encoder = cmd.encoder
+ }
go func() {
- for range time.Tick(10 * time.Second) {
+ for range time.Tick(10 * time.Minute) {
log.Printf("tilelib.Len() == %d", tilelib.Len())
}
}()
- w := bufio.NewWriter(stdout)
- cmd.encoder = gob.NewEncoder(w)
- err = cmd.tileGVCFs(tilelib, infiles)
+
+ err = cmd.tileInputs(tilelib, infiles)
if err != nil {
return 1
}
- err = w.Flush()
+ err = bufw.Flush()
+ if err != nil {
+ return 1
+ }
+ err = output.Close()
if err != nil {
return 1
}
return tilelib.TileFasta(infile, input)
}
-func (cmd *importer) loadTileLibrary() (*tileLibrary, error) {
+func (cmd *importer) loadTagLibrary() (*tagLibrary, error) {
log.Printf("tag library %s load starting", cmd.tagLibraryFile)
f, err := os.Open(cmd.tagLibraryFile)
if err != nil {
return nil, fmt.Errorf("cannot tile: tag library is empty")
}
log.Printf("tag library %s load done", cmd.tagLibraryFile)
- return &tileLibrary{taglib: &taglib}, nil
+ return &taglib, nil
}
+var (
+ vcfFilenameRe = regexp.MustCompile(`\.vcf(\.gz)?$`)
+ fasta1FilenameRe = regexp.MustCompile(`\.1\.fa(sta)?(\.gz)?$`)
+ fasta2FilenameRe = regexp.MustCompile(`\.2\.fa(sta)?(\.gz)?$`)
+ fastaFilenameRe = regexp.MustCompile(`\.fa(sta)?(\.gz)?$`)
+)
+
func listInputFiles(paths []string) (files []string, err error) {
for _, path := range paths {
if fi, err := os.Stat(path); err != nil {
return nil, fmt.Errorf("%s: stat failed: %s", path, err)
} else if !fi.IsDir() {
- if !strings.HasSuffix(path, ".2.fasta") || strings.HasSuffix(path, ".2.fasta.gz") {
+ if !fasta2FilenameRe.MatchString(path) {
files = append(files, path)
}
continue
}
sort.Strings(names)
for _, name := range names {
- if strings.HasSuffix(name, ".vcf") || strings.HasSuffix(name, ".vcf.gz") {
+ if vcfFilenameRe.MatchString(name) {
files = append(files, filepath.Join(path, name))
- } else if strings.HasSuffix(name, ".1.fasta") || strings.HasSuffix(name, ".1.fasta.gz") {
+ } else if fastaFilenameRe.MatchString(name) && !fasta2FilenameRe.MatchString(name) {
files = append(files, filepath.Join(path, name))
}
}
d.Close()
}
for _, file := range files {
- if strings.HasSuffix(file, ".1.fasta") || strings.HasSuffix(file, ".1.fasta.gz") {
- continue
- } else if _, err := os.Stat(file + ".csi"); err == nil {
- continue
- } else if _, err = os.Stat(file + ".tbi"); err == nil {
+ if fastaFilenameRe.MatchString(file) {
continue
+ } else if vcfFilenameRe.MatchString(file) {
+ if _, err := os.Stat(file + ".csi"); err == nil {
+ continue
+ } else if _, err = os.Stat(file + ".tbi"); err == nil {
+ continue
+ } else {
+ return nil, fmt.Errorf("%s: cannot read without .tbi or .csi index file", file)
+ }
} else {
- return nil, fmt.Errorf("%s: cannot read without .tbi or .csi index file", file)
+ return nil, fmt.Errorf("don't know how to handle filename %s", file)
}
}
return
}
-func (cmd *importer) tileGVCFs(tilelib *tileLibrary, infiles []string) error {
+func (cmd *importer) tileInputs(tilelib *tileLibrary, infiles []string) error {
starttime := time.Now()
errs := make(chan error, 1)
todo := make(chan func() error, len(infiles)*2)
var phases sync.WaitGroup
phases.Add(2)
variants := make([][]tileVariantID, 2)
- if strings.HasSuffix(infile, ".1.fasta") || strings.HasSuffix(infile, ".1.fasta.gz") {
+ if fasta1FilenameRe.MatchString(infile) {
todo <- func() error {
defer phases.Done()
log.Printf("%s starting", infile)
defer log.Printf("%s done", infile)
tseqs, err := cmd.tileFasta(tilelib, infile)
- variants[0] = tseqs.Variants()
+ var kept, dropped int
+ variants[0], kept, dropped = tseqs.Variants()
+ log.Printf("%s found %d unique tags plus %d repeats", infile, kept, dropped)
return err
}
- infile2 := regexp.MustCompile(`\.1\.fasta(\.gz)?$`).ReplaceAllString(infile, `.2.fasta$1`)
+ infile2 := fasta1FilenameRe.ReplaceAllString(infile, `.2.fa$1$2`)
todo <- func() error {
defer phases.Done()
log.Printf("%s starting", infile2)
defer log.Printf("%s done", infile2)
tseqs, err := cmd.tileFasta(tilelib, infile2)
- variants[1] = tseqs.Variants()
+ var kept, dropped int
+ variants[1], kept, dropped = tseqs.Variants()
+ log.Printf("%s found %d unique tags plus %d repeats", infile2, kept, dropped)
return err
}
- } else {
+ } else if fastaFilenameRe.MatchString(infile) {
+ todo <- func() error {
+ defer phases.Done()
+ defer phases.Done()
+ log.Printf("%s starting", infile)
+ defer log.Printf("%s done", infile)
+ tseqs, err := cmd.tileFasta(tilelib, infile)
+ if err != nil {
+ return err
+ }
+ totlen := 0
+ for _, tseq := range tseqs {
+ totlen += len(tseq)
+ }
+ log.Printf("%s tiled %d seqs, total len %d", infile, len(tseqs), totlen)
+ return cmd.encoder.Encode(LibraryEntry{
+ CompactSequences: []CompactSequence{{Name: infile, TileSequences: tseqs}},
+ })
+ }
+ // Don't write out a CompactGenomes entry
+ continue
+ } else if vcfFilenameRe.MatchString(infile) {
for phase := 0; phase < 2; phase++ {
phase := phase
todo <- func() error {
log.Printf("%s phase %d starting", infile, phase+1)
defer log.Printf("%s phase %d done", infile, phase+1)
tseqs, err := cmd.tileGVCF(tilelib, infile, phase)
- variants[phase] = tseqs.Variants()
+ var kept, dropped int
+ variants[phase], kept, dropped = tseqs.Variants()
+ log.Printf("%s phase %d found %d unique tags plus %d repeats", infile, phase+1, kept, dropped)
return err
}
}
+ } else {
+ panic(fmt.Sprintf("bug: unhandled filename %q", infile))
}
encodeJobs.Add(1)
go func() {
if len(errs) > 0 {
return
}
- ntags := len(variants[0])
- if ntags < len(variants[1]) {
- ntags = len(variants[1])
- }
- flat := make([]tileVariantID, ntags*2)
- for i := 0; i < ntags; i++ {
- flat[i*2] = variants[0][i]
- flat[i*2+1] = variants[1][i]
- }
err := cmd.encoder.Encode(LibraryEntry{
- CompactGenomes: []CompactGenome{{Name: infile, Variants: flat}},
+ CompactGenomes: []CompactGenome{{Name: infile, Variants: flatten(variants)}},
})
if err != nil {
select {
}
go close(todo)
var tileJobs sync.WaitGroup
- running := int64(runtime.NumCPU())
- for i := 0; i < runtime.NumCPU(); i++ {
+ var running int64
+ for i := 0; i < runtime.NumCPU()*9/8+1; i++ {
tileJobs.Add(1)
+ atomic.AddInt64(&running, 1)
go func() {
defer tileJobs.Done()
defer atomic.AddInt64(&running, -1)
}
func (cmd *importer) tileGVCF(tilelib *tileLibrary, infile string, phase int) (tileseq tileSeq, err error) {
+ if cmd.refFile == "" {
+ err = errors.New("cannot import vcf: reference data (-ref) not specified")
+ return
+ }
args := []string{"bcftools", "consensus", "--fasta-ref", cmd.refFile, "-H", fmt.Sprint(phase + 1), infile}
indexsuffix := ".tbi"
if _, err := os.Stat(infile + ".csi"); err == nil {
}
return
}
+
+func flatten(variants [][]tileVariantID) []tileVariantID {
+ ntags := 0
+ for _, v := range variants {
+ if ntags < len(v) {
+ ntags = len(v)
+ }
+ }
+ flat := make([]tileVariantID, ntags*2)
+ for i := 0; i < ntags; i++ {
+ for hap := 0; hap < 2; hap++ {
+ if i < len(variants[hap]) {
+ flat[i*2+hap] = variants[hap][i]
+ }
+ }
+ }
+ return flat
+}
projects / lightning.git / blobdiff