"golang.org/x/crypto/blake2b"
)
+const annotationMaxTileSpan = 100
+
type sliceNumpy struct {
filter filter
threads int
cmd.cgnames = nil
var tagset [][]byte
- DecodeLibrary(in0, strings.HasSuffix(infiles[0], ".gz"), func(ent *LibraryEntry) error {
+ err = DecodeLibrary(in0, strings.HasSuffix(infiles[0], ".gz"), func(ent *LibraryEntry) error {
if len(ent.TagSet) > 0 {
tagset = ent.TagSet
}
if err != nil {
return 1
}
+ if len(cmd.cgnames) == 0 {
+ err = fmt.Errorf("fatal: 0 cases, 0 controls, nothing to do")
+ return 1
+ }
cmd.minCoverage = int(math.Ceil(cmd.filter.MinCoverage * float64(len(cmd.cgnames))))
{
seqname string // chr1
pos int // distance from start of chromosome to starttag
tiledata []byte // acgtggcaa...
+ excluded bool // true if excluded by regions file
+ nexttag tagID // tagID of following tile (-1 for last tag of chromosome)
}
isdup := map[tagID]bool{}
reftile := map[tagID]*reftileinfo{}
for seqname, cseq := range refseq {
pos := 0
+ lastreftag := tagID(-1)
for _, libref := range cseq {
if cmd.filter.MaxTag >= 0 && libref.Tag > tagID(cmd.filter.MaxTag) {
continue
variant: libref.Variant,
tiledata: tiledata,
pos: pos,
+ nexttag: -1,
}
+ if lastreftag >= 0 {
+ reftile[lastreftag].nexttag = libref.Tag
+ }
+ lastreftag = libref.Tag
}
pos += len(tiledata) - taglen
}
}
log.Printf("before applying mask, len(reftile) == %d", len(reftile))
log.Printf("deleting reftile entries for regions outside %d intervals", mask.Len())
- for tag, rt := range reftile {
+ for _, rt := range reftile {
if !mask.Check(strings.TrimPrefix(rt.seqname, "chr"), rt.pos, rt.pos+len(rt.tiledata)) {
- delete(reftile, tag)
+ rt.excluded = true
}
}
log.Printf("after applying mask, len(reftile) == %d", len(reftile))
if tv.Ref {
continue
}
+ // Skip tile with no
+ // corresponding ref tile, if
+ // mask is in play (we can't
+ // determine coordinates for
+ // these)
if mask != nil && reftile[tv.Tag] == nil {
- // Don't waste
- // time/memory on
- // masked-out tiles.
+ continue
+ }
+ // Skip tile whose
+ // corresponding ref tile is
+ // outside target regions --
+ // unless it's a potential
+ // spanning tile.
+ if mask != nil && reftile[tv.Tag].excluded &&
+ (int(tv.Tag+1) >= len(tagset) ||
+ (bytes.HasSuffix(tv.Sequence, tagset[tv.Tag+1]) && reftile[tv.Tag+1] != nil && !reftile[tv.Tag+1].excluded)) {
continue
}
if tv.Tag == cmd.debugTag {
for tag := tagstart; tag < tagend; tag++ {
rt := reftile[tag]
if rt == nil && mask != nil {
- // Excluded by specified regions
+ // With no ref tile, we don't
+ // have coordinates to say
+ // this is in the desired
+ // regions -- so it's not.
+ // TODO: handle ref spanning
+ // tile case.
+ continue
+ }
+ if rt != nil && rt.excluded {
+ // TODO: don't skip yet --
+ // first check for spanning
+ // tile variants that
+ // intersect non-excluded ref
+ // tiles.
continue
}
if cmd.filter.MaxTag >= 0 && tag > tagID(cmd.filter.MaxTag) {
}
}
if *onehotChunked || *onehotSingle {
- onehot, xrefs := cmd.tv2homhet(cgs, maxv, remap, tag, tagstart)
+ onehot, xrefs := cmd.tv2homhet(cgs, maxv, remap, tag, tagstart, seq)
if tag == cmd.debugTag {
log.WithFields(logrus.Fields{
"onehot": onehot,
if rt == nil {
// Reference does not use any
// variant of this tile
+ //
+ // TODO: diff against the
+ // relevant portion of the
+ // ref's spanning tile
outcol++
continue
}
} else {
done[v] = true
}
- if len(tv.Sequence) < taglen || !bytes.HasSuffix(rt.tiledata, tv.Sequence[len(tv.Sequence)-taglen:]) {
+ if len(tv.Sequence) < taglen {
+ continue
+ }
+ // if reftilestr doesn't end
+ // in the same tag as tv,
+ // extend reftilestr with
+ // following ref tiles until
+ // it does (up to an arbitrary
+ // sanity-check limit)
+ reftilestr := reftilestr
+ endtagstr := strings.ToUpper(string(tv.Sequence[len(tv.Sequence)-taglen:]))
+ for i, rt := 0, rt; i < annotationMaxTileSpan && !strings.HasSuffix(reftilestr, endtagstr) && rt.nexttag >= 0; i++ {
+ rt = reftile[rt.nexttag]
+ reftilestr += strings.ToUpper(string(rt.tiledata[taglen:]))
+ }
+ if mask != nil && !mask.Check(strings.TrimPrefix(rt.seqname, "chr"), rt.pos, rt.pos+len(reftilestr)) {
+ continue
+ }
+ if !strings.HasSuffix(reftilestr, endtagstr) {
fmt.Fprintf(annow, "%d,%d,%d,,%s,%d,,,\n", tag, outcol, v, rt.seqname, rt.pos)
continue
}
- if lendiff := len(rt.tiledata) - len(tv.Sequence); lendiff < -1000 || lendiff > 1000 {
+ if lendiff := len(reftilestr) - len(tv.Sequence); lendiff < -1000 || lendiff > 1000 {
fmt.Fprintf(annow, "%d,%d,%d,,%s,%d,,,\n", tag, outcol, v, rt.seqname, rt.pos)
continue
}
if cmd.filter.MaxTag >= 0 && tag > tagID(cmd.filter.MaxTag) {
break
}
- if mask != nil && reftile[tag] == nil {
+ if rt := reftile[tag]; rt == nil || rt.excluded {
continue
}
if v == 0 {
return 1
}
fnm = fmt.Sprintf("%s/onehot-columns.npy", *outputDir)
- err = writeNumpyInt32(fnm, onehotXref2int32(xrefs), 4, len(xrefs))
+ err = writeNumpyInt32(fnm, onehotXref2int32(xrefs), 5, len(xrefs))
if err != nil {
return 1
}
// variants of a single tile/tag#.
//
// Return nil if no tile variant passes Χ² filter.
-func (cmd *sliceNumpy) tv2homhet(cgs map[string]CompactGenome, maxv tileVariantID, remap []tileVariantID, tag, chunkstarttag tagID) ([][]int8, []onehotXref) {
+func (cmd *sliceNumpy) tv2homhet(cgs map[string]CompactGenome, maxv tileVariantID, remap []tileVariantID, tag, chunkstarttag tagID, seq map[tagID][]TileVariant) ([][]int8, []onehotXref) {
if tag == cmd.debugTag {
tv := make([]tileVariantID, len(cmd.cgnames)*2)
for i, name := range cmd.cgnames {
tagoffset := tag - chunkstarttag
coverage := 0
for _, cg := range cgs {
- if cg.Variants[tagoffset*2] > 0 && cg.Variants[tagoffset*2+1] > 0 {
+ alleles := 0
+ for _, v := range cg.Variants[tagoffset*2 : tagoffset*2+2] {
+ if v > 0 && int(v) < len(seq[tag]) && len(seq[tag][v].Sequence) > 0 {
+ alleles++
+ }
+ }
+ if alleles == 2 {
coverage++
}
}
// P-value row contains 1000000x actual p-value.
func onehotXref2int32(xrefs []onehotXref) []int32 {
xcols := len(xrefs)
- xdata := make([]int32, 4*xcols)
+ xdata := make([]int32, 5*xcols)
for i, xref := range xrefs {
xdata[i] = int32(xref.tag)
xdata[xcols+i] = int32(xref.variant)
xdata[xcols*2+i] = 1
}
xdata[xcols*3+i] = int32(xref.pvalue * 1000000)
+ xdata[xcols*4+i] = int32(-math.Log10(xref.pvalue) * 1000000)
}
return xdata
}
projects / lightning.git / blobdiff