anno2vcf.go

   1 // Copyright (C) The Lightning Authors. All rights reserved.
   2 //
   3 // SPDX-License-Identifier: AGPL-3.0
   4
   5 package lightning
   6
   7 import (
   8         "bufio"
   9         "bytes"
  10         "flag"
  11         "fmt"
  12         "io"
  13         "io/ioutil"
  14         "net/http"
  15         _ "net/http/pprof"
  16         "os"
  17         "runtime"
  18         "sort"
  19         "strconv"
  20         "strings"
  21         "sync"
  22
  23         "git.arvados.org/arvados.git/sdk/go/arvados"
  24         log "github.com/sirupsen/logrus"
  25 )
  26
  27 type anno2vcf struct {
  28 }
  29
  30 func (cmd *anno2vcf) RunCommand(prog string, args []string, stdin io.Reader, stdout, stderr io.Writer) int {
  31         var err error
  32         defer func() {
  33                 if err != nil {
  34                         fmt.Fprintf(stderr, "%s\n", err)
  35                 }
  36         }()
  37         flags := flag.NewFlagSet("", flag.ContinueOnError)
  38         flags.SetOutput(stderr)
  39         pprof := flags.String("pprof", "", "serve Go profile data at http://`[addr]:port`")
  40         runlocal := flags.Bool("local", false, "run on local host (default: run in an arvados container)")
  41         projectUUID := flags.String("project", "", "project `UUID` for output data")
  42         priority := flags.Int("priority", 500, "container request priority")
  43         inputDir := flags.String("input-dir", "./in", "input `directory`")
  44         outputDir := flags.String("output-dir", "./out", "output `directory`")
  45         err = flags.Parse(args)
  46         if err == flag.ErrHelp {
  47                 err = nil
  48                 return 0
  49         } else if err != nil {
  50                 return 2
  51         }
  52
  53         if *pprof != "" {
  54                 go func() {
  55                         log.Println(http.ListenAndServe(*pprof, nil))
  56                 }()
  57         }
  58
  59         if !*runlocal {
  60                 runner := arvadosContainerRunner{
  61                         Name:        "lightning anno2vcf",
  62                         Client:      arvados.NewClientFromEnv(),
  63                         ProjectUUID: *projectUUID,
  64                         RAM:         500000000000,
  65                         VCPUs:       64,
  66                         Priority:    *priority,
  67                         KeepCache:   2,
  68                         APIAccess:   true,
  69                 }
  70                 err = runner.TranslatePaths(inputDir)
  71                 if err != nil {
  72                         return 1
  73                 }
  74                 runner.Args = []string{"anno2vcf", "-local=true",
  75                         "-pprof", ":6060",
  76                         "-input-dir", *inputDir,
  77                         "-output-dir", "/mnt/output",
  78                 }
  79                 var output string
  80                 output, err = runner.Run()
  81                 if err != nil {
  82                         return 1
  83                 }
  84                 fmt.Fprintln(stdout, output)
  85                 return 0
  86         }
  87
  88         d, err := open(*inputDir)
  89         if err != nil {
  90                 log.Print(err)
  91                 return 1
  92         }
  93         defer d.Close()
  94         fis, err := d.Readdir(-1)
  95         if err != nil {
  96                 log.Print(err)
  97                 return 1
  98         }
  99         d.Close()
 100         sort.Slice(fis, func(i, j int) { return fis[i].Name() < fis[j].Name() })
 101
 102         type call struct {
 103                 tile      int
 104                 variant   int
 105                 sequence  []byte
 106                 position  int
 107                 deletion  []byte
 108                 insertion []byte
 109         }
 110         var allcalls []*call
 111         var mtx sync.Mutex
 112         throttle := throttle{Max: runtime.GOMAXPROCS(0)}
 113         log.Print("reading input files")
 114         for _, fi := range fis {
 115                 if !strings.HasSuffix(fi.Name(), "annotations.csv") {
 116                         continue
 117                 }
 118                 filename := *inputDir + "/" + fi.Name()
 119                 throttle.Acquire()
 120                 go func() {
 121                         defer throttle.Release()
 122                         log.Printf("reading %s", filename)
 123                         buf, err := ioutil.ReadFile(filename)
 124                         if err != nil {
 125                                 throttle.Report(fmt.Errorf("%s: %s", filename, err))
 126                                 return
 127                         }
 128                         lines := bytes.Split(buf, []byte{'\n'})
 129                         calls := make([]*call, 0, len(lines))
 130                         for lineIdx, line := range lines {
 131                                 if len(line) == 0 {
 132                                         continue
 133                                 }
 134                                 if lineIdx & ^0xfff == 0 && throttle.Err() != nil {
 135                                         return
 136                                 }
 137                                 fields := bytes.Split(line, []byte{','})
 138                                 if len(fields) != 8 {
 139                                         throttle.Report(fmt.Errorf("%s line %d: wrong number of fields (%d != %d): %q", fi.Name(), lineIdx+1, len(fields), 8, line))
 140                                         return
 141                                 }
 142                                 tile, _ := strconv.ParseInt(string(fields[0]), 10, 64)
 143                                 variant, _ := strconv.ParseInt(string(fields[2]), 10, 64)
 144                                 position, _ := strconv.ParseInt(string(fields[5]), 10, 64)
 145                                 calls = append(calls, &call{
 146                                         tile:      int(tile),
 147                                         variant:   int(variant),
 148                                         sequence:  append([]byte(nil), fields[4]...),
 149                                         position:  int(position),
 150                                         deletion:  append([]byte(nil), fields[6]...),
 151                                         insertion: append([]byte(nil), fields[7]...),
 152                                 })
 153                         }
 154                         mtx.Lock()
 155                         allcalls = append(allcalls, calls...)
 156                         mtx.Unlock()
 157                 }()
 158         }
 159         throttle.Wait()
 160         if throttle.Err() != nil {
 161                 log.Print(throttle.Err())
 162                 return 1
 163         }
 164         log.Print("sorting")
 165         sort.Slice(allcalls, func(i, j int) bool {
 166                 ii, jj := allcalls[i], allcalls[j]
 167                 if cmp := bytes.Compare(ii.sequence, jj.sequence); cmp != 0 {
 168                         return cmp < 0
 169                 }
 170                 if cmp := ii.position - jj.position; cmp != 0 {
 171                         return cmp < 0
 172                 }
 173                 if cmp := len(ii.deletion) - len(jj.deletion); cmp != 0 {
 174                         return cmp < 0
 175                 }
 176                 if cmp := bytes.Compare(ii.insertion, jj.insertion); cmp != 0 {
 177                         return cmp < 0
 178                 }
 179                 if cmp := ii.tile - jj.tile; cmp != 0 {
 180                         return cmp < 0
 181                 }
 182                 return ii.variant < jj.variant
 183         })
 184
 185         vcfFilename := *outputDir + "/annotations.vcf"
 186         log.Printf("writing %s", vcfFilename)
 187         f, err := os.Create(vcfFilename)
 188         if err != nil {
 189                 return 1
 190         }
 191         defer f.Close()
 192         bufw := bufio.NewWriterSize(f, 1<<20)
 193         _, err = fmt.Fprintf(bufw, `##fileformat=VCFv4.0
 194 ##INFO=<ID=TV,Number=.,Type=String,Description="tile-variant">
 195 #CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO
 196 `)
 197         if err != nil {
 198                 return 1
 199         }
 200         placeholder := []byte{'.'}
 201         for i := 0; i < len(allcalls); {
 202                 call := allcalls[i]
 203                 i++
 204                 info := fmt.Sprintf("TV=,%d-%d,", call.tile, call.variant)
 205                 for i < len(allcalls) &&
 206                         bytes.Equal(call.sequence, allcalls[i].sequence) &&
 207                         call.position == allcalls[i].position &&
 208                         len(call.deletion) == len(allcalls[i].deletion) &&
 209                         bytes.Equal(call.insertion, allcalls[i].insertion) {
 210                         call = allcalls[i]
 211                         i++
 212                         info += fmt.Sprintf("%d-%d,", call.tile, call.variant)
 213                 }
 214                 deletion := call.deletion
 215                 if len(deletion) == 0 {
 216                         deletion = placeholder
 217                 }
 218                 insertion := call.insertion
 219                 if len(insertion) == 0 {
 220                         insertion = placeholder
 221                 }
 222                 _, err = fmt.Fprintf(bufw, "%s\t%d\t.\t%s\t%s\t.\t.\t%s\n", call.sequence, call.position, deletion, insertion, info)
 223                 if err != nil {
 224                         return 1
 225                 }
 226         }
 227         err = bufw.Flush()
 228         if err != nil {
 229                 return 1
 230         }
 231         err = f.Close()
 232         if err != nil {
 233                 return 1
 234         }
 235         return 0
 236 }