18 "github.com/kshedden/gonpy"
21 type gvcf2numpy struct {
27 func (cmd *gvcf2numpy) RunCommand(prog string, args []string, stdin io.Reader, stdout, stderr io.Writer) int {
31 fmt.Fprintf(stderr, "%s\n", err)
34 flags := flag.NewFlagSet("", flag.ContinueOnError)
35 flags.SetOutput(stderr)
36 flags.StringVar(&cmd.tagLibraryFile, "tag-library", "", "tag library fasta `file`")
37 flags.StringVar(&cmd.refFile, "ref", "", "reference fasta `file`")
38 err = flags.Parse(args)
39 if err == flag.ErrHelp {
42 } else if err != nil {
44 } else if cmd.refFile == "" || cmd.tagLibraryFile == "" {
45 fmt.Fprintln(os.Stderr, "cannot run without -tag-library and -ref arguments")
47 } else if flags.NArg() == 0 {
53 infiles, err := listVCFFiles(flags.Args())
58 log.Printf("tag library %s load starting", cmd.tagLibraryFile)
59 f, err := os.Open(cmd.tagLibraryFile)
63 var rdr io.ReadCloser = f
64 if strings.HasSuffix(cmd.tagLibraryFile, ".gz") {
65 rdr, err = gzip.NewReader(f)
67 err = fmt.Errorf("%s: gzip: %s", cmd.tagLibraryFile, err)
72 err = taglib.Load(rdr)
77 err = fmt.Errorf("cannot tile: tag library is empty")
80 log.Printf("tag library %s load done", cmd.tagLibraryFile)
82 tilelib := tileLibrary{taglib: &taglib}
83 tseqs, err := cmd.tileGVCFs(&tilelib, infiles)
87 err = cmd.printVariants(tseqs)
94 func listVCFFiles(paths []string) (files []string, err error) {
95 for _, path := range paths {
96 if fi, err := os.Stat(path); err != nil {
97 return nil, fmt.Errorf("%s: stat failed: %s", path, err)
98 } else if !fi.IsDir() {
99 files = append(files, path)
102 d, err := os.Open(path)
104 return nil, fmt.Errorf("%s: open failed: %s", path, err)
107 names, err := d.Readdirnames(0)
109 return nil, fmt.Errorf("%s: readdir failed: %s", path, err)
112 for _, name := range names {
113 if strings.HasSuffix(name, ".vcf") || strings.HasSuffix(name, ".vcf.gz") {
114 files = append(files, filepath.Join(path, name))
119 for _, file := range files {
120 if _, err := os.Stat(file + ".csi"); err == nil {
122 } else if _, err = os.Stat(file + ".tbi"); err == nil {
125 return nil, fmt.Errorf("%s: cannot read without .tbi or .csi index file", file)
131 func (cmd *gvcf2numpy) tileGVCFs(tilelib *tileLibrary, infiles []string) ([]tileSeq, error) {
132 limit := make(chan bool, runtime.NumCPU())
133 errs := make(chan error, 1)
134 tseqs := make([]tileSeq, len(infiles)*2)
135 var wg sync.WaitGroup
136 for i, infile := range infiles {
137 for phase := 0; phase < 2; phase++ {
139 go func(i int, infile string, phase int) {
142 defer func() { <-limit }()
143 log.Printf("%s phase %d starting", infile, phase+1)
144 defer log.Printf("%s phase %d done", infile, phase+1)
146 tseqs[i*2+phase], err = cmd.tileGVCF(tilelib, infile, phase)
161 if err := <-errs; err != nil {
167 func (cmd *gvcf2numpy) printVariants(tseqs []tileSeq) error {
169 for _, tseq := range tseqs {
170 for _, path := range tseq {
171 for _, tvar := range path {
172 if maxtag < tvar.tag {
178 out := make([]uint16, len(tseqs)*int(maxtag+1))
179 for i := 0; i < len(tseqs)/2; i++ {
180 for phase := 0; phase < 2; phase++ {
181 for _, path := range tseqs[i*2+phase] {
182 for _, tvar := range path {
183 out[2*int(tvar.tag)+phase] = uint16(tvar.variant)
188 w := bufio.NewWriter(cmd.output)
189 npw, err := gonpy.NewWriter(nopCloser{w})
193 npw.Shape = []int{len(tseqs) / 2, 2 * (int(maxtag) + 1)}
198 type nopCloser struct {
202 func (nopCloser) Close() error { return nil }
204 func (cmd *gvcf2numpy) tileGVCF(tilelib *tileLibrary, infile string, phase int) (tileseq tileSeq, err error) {
205 args := []string{"bcftools", "consensus", "--fasta-ref", cmd.refFile, "-H", fmt.Sprint(phase + 1), infile}
206 indexsuffix := ".tbi"
207 if _, err := os.Stat(infile + ".csi"); err == nil {
210 if out, err := exec.Command("docker", "image", "ls", "-q", "lightning-runtime").Output(); err == nil && len(out) > 0 {
211 args = append([]string{
212 "docker", "run", "--rm",
214 "--volume=" + infile + ":" + infile + ":ro",
215 "--volume=" + infile + indexsuffix + ":" + infile + indexsuffix + ":ro",
216 "--volume=" + cmd.refFile + ":" + cmd.refFile + ":ro",
220 consensus := exec.Command(args[0], args[1:]...)
221 consensus.Stderr = os.Stderr
222 stdout, err := consensus.StdoutPipe()
226 err = consensus.Start()
230 tileseq, err = tilelib.TileFasta(fmt.Sprintf("%s phase %d", infile, phase+1), stdout)
238 err = consensus.Wait()
240 err = fmt.Errorf("%s phase %d: bcftools: %s", infile, phase, err)