tilelib.go

   1 package main
   2
   3 import (
   4         "bufio"
   5         "bytes"
   6         "io"
   7         "strings"
   8         "sync"
   9
  10         log "github.com/sirupsen/logrus"
  11         "golang.org/x/crypto/blake2b"
  12 )
  13
  14 type tileVariantID uint16 // 1-based
  15
  16 type tileLibRef struct {
  17         tag     tagID
  18         variant tileVariantID
  19 }
  20
  21 type tileSeq map[string][]tileLibRef
  22
  23 func (tseq tileSeq) Variants() []tileVariantID {
  24         maxtag := 0
  25         for _, refs := range tseq {
  26                 for _, ref := range refs {
  27                         if maxtag < int(ref.tag) {
  28                                 maxtag = int(ref.tag)
  29                         }
  30                 }
  31         }
  32         vars := make([]tileVariantID, maxtag+1)
  33         for _, refs := range tseq {
  34                 for _, ref := range refs {
  35                         vars[int(ref.tag)] = ref.variant
  36                 }
  37         }
  38         return vars
  39 }
  40
  41 type tileLibrary struct {
  42         skipOOO bool
  43         taglib  *tagLibrary
  44         variant [][][blake2b.Size256]byte
  45         // count [][]int
  46         // seq map[[blake2b.Size]byte][]byte
  47         variants int
  48
  49         mtx sync.Mutex
  50 }
  51
  52 func (tilelib *tileLibrary) TileFasta(filelabel string, rdr io.Reader) (tileSeq, error) {
  53         ret := tileSeq{}
  54         type jobT struct {
  55                 label string
  56                 fasta []byte
  57         }
  58         todo := make(chan jobT)
  59         scanner := bufio.NewScanner(rdr)
  60         go func() {
  61                 defer close(todo)
  62                 var fasta []byte
  63                 var seqlabel string
  64                 for scanner.Scan() {
  65                         buf := scanner.Bytes()
  66                         if len(buf) == 0 || buf[0] == '>' {
  67                                 todo <- jobT{seqlabel, fasta}
  68                                 seqlabel, fasta = string(buf[1:]), nil
  69                                 log.Debugf("%s %s reading fasta", filelabel, seqlabel)
  70                         } else {
  71                                 fasta = append(fasta, bytes.ToLower(buf)...)
  72                         }
  73                 }
  74                 todo <- jobT{seqlabel, fasta}
  75         }()
  76         type foundtag struct {
  77                 pos    int
  78                 tagid  tagID
  79                 taglen int
  80         }
  81         found := make([]foundtag, 2000000)
  82         path := make([]tileLibRef, 2000000)
  83         totalFoundTags := 0
  84         totalPathLen := 0
  85         skippedSequences := 0
  86         for job := range todo {
  87                 if len(job.fasta) == 0 {
  88                         continue
  89                 } else if strings.Contains(job.label, "_") {
  90                         skippedSequences++
  91                         continue
  92                 }
  93                 log.Debugf("%s %s tiling", filelabel, job.label)
  94
  95                 found = found[:0]
  96                 tilelib.taglib.FindAll(job.fasta, func(tagid tagID, pos, taglen int) {
  97                         found = append(found, foundtag{pos: pos, tagid: tagid, taglen: taglen})
  98                 })
  99                 totalFoundTags += len(found)
 100
 101                 path = path[:0]
 102                 last := foundtag{tagid: -1}
 103                 for i, f := range found {
 104                         log.Tracef("%s %s found[%d] == %#v", filelabel, job.label, i, f)
 105                         if tilelib.skipOOO {
 106                                 if f.tagid < last.tagid+1 {
 107                                         log.Debugf("%s %s skipped out-of-order tag %d (found at %d) because it appears after tag %d (found at %d)", filelabel, job.label, f.tagid, f.pos, last.tagid, last.pos)
 108                                         continue
 109                                 }
 110                                 if f.tagid > last.tagid+1 && // accepting this tag would mean skipping some tags
 111                                         i+1 < len(found) && // there is a "next" found tag after this one
 112                                         found[i+1].tagid > last.tagid && // next found tag is usable (we haven't already passed it in accepted sequence)
 113                                         found[i+1].tagid <= f.tagid { // next found tag is expected before this one (so we can't use both)
 114                                         log.Debugf("%s %s skipped out-of-order tag %d (found at %d) because it appears between tag %d (found at %d) and %d (found at %d)", filelabel, job.label, f.tagid, f.pos, last.tagid, last.pos, found[i+1].tagid, found[i+1].pos)
 115                                         continue
 116                                 }
 117                         }
 118                         if last.taglen > 0 {
 119                                 path = append(path, tilelib.getRef(last.tagid, job.fasta[last.pos:f.pos+f.taglen]))
 120                         }
 121                         last = f
 122                 }
 123                 if last.taglen > 0 {
 124                         path = append(path, tilelib.getRef(last.tagid, job.fasta[last.pos:]))
 125                 }
 126
 127                 pathcopy := make([]tileLibRef, len(path))
 128                 copy(pathcopy, path)
 129                 ret[job.label] = pathcopy
 130                 log.Debugf("%s %s tiled with path len %d, skipped %d", filelabel, job.label, len(path), len(found)-len(path))
 131                 totalPathLen += len(path)
 132         }
 133         log.Printf("%s tiled with total path len %d in %d sequences (skipped %d sequences with '_' in name, skipped %d out-of-order tags)", filelabel, totalPathLen, len(ret), skippedSequences, totalFoundTags-totalPathLen)
 134         return ret, scanner.Err()
 135 }
 136
 137 func (tilelib *tileLibrary) Len() int {
 138         tilelib.mtx.Lock()
 139         defer tilelib.mtx.Unlock()
 140         return tilelib.variants
 141 }
 142
 143 // Return a tileLibRef for a tile with the given tag and sequence,
 144 // adding the sequence to the library if needed.
 145 func (tilelib *tileLibrary) getRef(tag tagID, seq []byte) tileLibRef {
 146         for _, b := range seq {
 147                 if b != 'a' && b != 'c' && b != 'g' && b != 't' {
 148                         // return "tile not found" if seq has any
 149                         // no-calls
 150                         return tileLibRef{tag: tag}
 151                 }
 152         }
 153         tilelib.mtx.Lock()
 154         defer tilelib.mtx.Unlock()
 155         // if tilelib.seq == nil {
 156         //      tilelib.seq = map[[blake2b.Size]byte][]byte{}
 157         // }
 158         if tilelib.variant == nil {
 159                 tilelib.variant = make([][][blake2b.Size256]byte, tilelib.taglib.Len())
 160         }
 161         seqhash := blake2b.Sum256(seq)
 162         for i, varhash := range tilelib.variant[tag] {
 163                 if varhash == seqhash {
 164                         return tileLibRef{tag: tag, variant: tileVariantID(i + 1)}
 165                 }
 166         }
 167         tilelib.variants++
 168         tilelib.variant[tag] = append(tilelib.variant[tag], seqhash)
 169         // tilelib.seq[seqhash] = append([]byte(nil), seq...)
 170         return tileLibRef{tag: tag, variant: tileVariantID(len(tilelib.variant[tag]))}
 171 }