maxFrequency float64
pcaComponents int
minCoverage int
+ minCoverageAll bool
includeVariant1 bool
debugTag tagID
flags.IntVar(&cmd.pcaComponents, "pca-components", 4, "number of PCA components to compute / use in logistic regression")
maxPCATiles := flags.Int("max-pca-tiles", 0, "maximum tiles to use as PCA input (filter, then drop every 2nd colum pair until below max)")
debugTag := flags.Int("debug-tag", -1, "log debugging details about specified tag")
+ flags.BoolVar(&cmd.minCoverageAll, "min-coverage-all", false, "apply -min-coverage filter based on all samples, not just training set")
flags.IntVar(&cmd.threads, "threads", 16, "number of memory-hungry assembly threads, and number of VCPUs to request for arvados container")
flags.Float64Var(&cmd.chi2PValue, "chi2-p-value", 1, "do Χ² test (or logistic regression if -samples file has PCA components) and omit columns with p-value above this threshold")
flags.Float64Var(&cmd.pvalueMinFrequency, "pvalue-min-frequency", 0.01, "skip p-value calculation on tile variants below this frequency in the training set")
"-chunked-onehot=" + fmt.Sprintf("%v", *onehotChunked),
"-samples=" + *samplesFilename,
"-case-control-only=" + fmt.Sprintf("%v", *caseControlOnly),
+ "-min-coverage-all=" + fmt.Sprintf("%v", cmd.minCoverageAll),
"-pca=" + fmt.Sprintf("%v", *onlyPCA),
"-pca-components=" + fmt.Sprintf("%d", cmd.pcaComponents),
"-max-pca-tiles=" + fmt.Sprintf("%d", *maxPCATiles),
}
}
}
- if cmd.filter.MinCoverage == 1 {
- // In the generic formula below, floating point
- // arithmetic can effectively push the coverage
- // threshold above 1.0, which is impossible/useless.
- // 1.0 needs to mean exactly 100% coverage.
+
+ if cmd.minCoverageAll {
cmd.minCoverage = len(cmd.cgnames)
} else {
- cmd.minCoverage = int(math.Ceil(cmd.filter.MinCoverage * float64(len(cmd.cgnames))))
+ cmd.minCoverage = cmd.trainingSetSize
+ }
+ if cmd.filter.MinCoverage < 1 {
+ cmd.minCoverage = int(math.Ceil(cmd.filter.MinCoverage * float64(cmd.minCoverage)))
}
if len(cmd.samples[0].pcaComponents) > 0 {
return err
}
foundthistag := false
- taglib.FindAll(tiledata[:len(tiledata)-1], func(tagid tagID, offset, _ int) {
+ taglib.FindAll(bufio.NewReader(bytes.NewReader(tiledata[:len(tiledata)-1])), nil, func(tagid tagID, offset, _ int) {
if !foundthistag && tagid == libref.Tag {
foundthistag = true
return
if err == errSkip {
return nil
} else if err != nil {
- return fmt.Errorf("%04d: DecodeLibrary(%s): err", infileIdx, infile)
+ return fmt.Errorf("%04d: DecodeLibrary(%s): %w", infileIdx, infile, err)
}
tagstart := cgs[cmd.cgnames[0]].StartTag
tagend := cgs[cmd.cgnames[0]].EndTag
count[blake2b.Sum256(rt.tiledata)] = 0
}
- for cgname, cg := range cgs {
+ for cgidx, cgname := range cmd.cgnames {
+ if !cmd.minCoverageAll && !cmd.samples[cgidx].isTraining {
+ continue
+ }
+ cg := cgs[cgname]
idx := int(tag-tagstart) * 2
for allele := 0; allele < 2; allele++ {
v := cg.Variants[idx+allele]
if cmd.filter.MaxTag >= 0 && tag > tagID(cmd.filter.MaxTag) {
break
}
- if rt := reftile[tag]; rt == nil || rt.excluded {
+ if rt := reftile[tag]; mask != nil && (rt == nil || rt.excluded) {
continue
}
if v == 0 {
for i := range cmd.samples {
cmd.samples[i].pcaComponents = make([]float64, outcols)
for c := 0; c < outcols; c++ {
- cmd.samples[i].pcaComponents[i] = pca.At(i, c)
+ cmd.samples[i].pcaComponents[c] = pca.At(i, c)
}
}
log.Print("done")
}
tagoffset := tag - chunkstarttag
coverage := 0
- for _, cg := range cgs {
+ for cgidx, cgname := range cmd.cgnames {
+ if !cmd.minCoverageAll && !cmd.samples[cgidx].isTraining {
+ continue
+ }
+ cg := cgs[cgname]
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 {