func ReadFai(path string, chrom string) []*sam.Reference {
	f, err := os.Open(path)
	if err != nil {
		fmt.Fprintf(os.Stderr, "error opening fai: %s. Is is present?\n", path)
		panic(err)
	}
	idx, err := fai.ReadFrom(f)
	if err != nil {
		fmt.Fprintf(os.Stderr, "error opening fai: %s. are you sure this i  s a valid fasta index?\n", path)
		panic(err)
	}
	recs := make([]fai.Record, 0, len(idx))
	for _, rec := range idx {
		recs = append(recs, rec)
	}
	sort.Slice(recs, func(i, j int) bool { return recs[i].Start < recs[j].Start })
	refs := make([]*sam.Reference, 0, len(idx))
	for _, rec := range recs {
		if chrom != "" && rec.Name != chrom {
			continue
		}
		ref, err := sam.NewReference(rec.Name, "", "", rec.Length, nil, nil)
		if err != nil {
			panic(err)
		}
		refs = append(refs, ref)
	}
	if len(refs) == 0 {
		if chrom != "" {
			log.Printf("ERROR: didn't find %s in %s", chrom, path)
		}
		panic(fmt.Sprintf("ERROR: didn't find any usable chromosomes in %s", path))
	}
	h, err := sam.NewHeader(nil, refs)
	if err != nil {
		panic(err)
	}
	return h.Refs()
}
ReadFai returns a slit of references from the fai path.
If chrom is "" all chromosomes are returned.
package indexcov
import (
	"bufio"
	"compress/gzip"
	"fmt"
	"html/template"
	"io"
	"log"
	"os"
	"os/exec"
	"path/filepath"
	"reflect"
	"regexp"
	"sort"
	"strconv"
	"strings"
	"sync"
	"time"
	"gonum.org/v1/gonum/floats"
	"gonum.org/v1/gonum/mat"
	"gonum.org/v1/gonum/stat"
	arg "github.com/alexflint/go-arg"
	"github.com/biogo/biogo/io/seqio/fai"
	"github.com/biogo/hts/bam"
	"github.com/biogo/hts/bgzf"
	"github.com/biogo/hts/sam"
	chartjs "github.com/brentp/go-chartjs"
	"github.com/brentp/go-chartjs/types"
	"github.com/brentp/goleft"
	"github.com/brentp/goleft/indexcov/crai"
	"github.com/brentp/xopen"
)
var Ploidy = 2
var cli = &struct {
	Directory      string         `arg:"-d,required,help:directory for output files"`
	IncludeGL      bool           `arg:"-e,help:plot GL chromosomes like: GL000201.1 which are not plotted by default"`
	ExcludePatt    string         `arg:-p,help:regular expression of chromosome names to exclude"`
	Sex            string         `arg:"-X,help:comma delimited names of the sex chromosome(s) used to infer sex. Set to '' if no sex chromosomes are present."`
	Chrom          string         `arg:"-c,help:optional chromosome to extract depth. default is entire genome."`
	Fai            string         `arg:"-f,help:fasta index file. Required when crais are used."`
	ExtraNormalize bool           `arg:"-n,help:normalize across samples and do local smoothign within sample. this is recommended for CRAI"`
	Bam            []string       `arg:"positional,required,help:bam(s) or crais for which to estimate coverage"`
	sex            []string       `arg:"-"`
	exclude        *regexp.Regexp `arg:"-"`
}{Sex: "X,Y", ExcludePatt: `^chrEBV$|^NC|_random$|Un_|^HLA\-|_alt$|hap\d$`}
// MaxCN is the maximum normalized value.
var MaxCN = float32(8)
// Index wraps a bam.Index to cache calculated values.
type Index struct {
	*bam.Index
	crai *crai.Index
	path string
	//mu                *sync.RWMutex
	medianSizePerTile float64
	sizes             [][]int64
	mapped            uint64
	unmapped          uint64
}
func (i *Index) Path() string {
	return i.path
}
// Sizes returns the size of each block in slices of chromosomes.
func (i *Index) Sizes() [][]int64 {
	return i.sizes
}
func vOffset(o bgzf.Offset) int64 {
	return o.File<<16 | int64(o.Block)
}
// init sets the medianSizePerTile
func (x *Index) init() {
	if x.Index != nil {
		x.sizes, x.mapped, x.unmapped = getSizes(x.Index)
		x.Index = nil
	} else if x.crai != nil {
		x.sizes = x.crai.Sizes()
		if x.sizes == nil {
			log.Fatal("bad index:", x.path)
		}
		x.crai = nil
	}
	// sizes is used to get the median.
	sizes := make([]int64, 0, 16384)
	for k := 0; k < len(x.sizes); k++ {
		sizes = append(sizes, x.sizes[k]...)
	}
	if len(sizes) < 1 {
		log.Fatalf("indexcov: no usable chromsomes in bam: %s", x.path)
	}
	sort.Slice(sizes, func(i, j int) bool { return sizes[i] < sizes[j] })
	total := int64(0)
	cumsum := make([]int64, len(sizes))
	// some bins with very high coverage can screw the normalization.
	// so we cap at the 98th percentile and then take the bin that passed then
	// median.
	n98 := sizes[int(0.98*float64(len(sizes)))]
	for i, s := range sizes {
		if s > n98 {
			s = n98
		}
		total += s
		cumsum[i] = total
	}
	var idx = sort.Search(len(cumsum), func(i int) bool { return (cumsum[i]) > (total / 2) })
	for idx >= len(sizes) {
		idx--
	}
	x.medianSizePerTile = float64(sizes[idx])
}
// NormalizedDepth returns a list of numbers for the normalized depth of the given region.
// Values are scaled to have a mean of 1. If end is 0, the full chromosome is returned.
func (x *Index) NormalizedDepth(refID int) []float32 {
	if x.medianSizePerTile == 0.0 {
		x.init()
	}
	if refID >= len(x.sizes) {
		return make([]float32, 0)
	}
	ref := x.sizes[refID]
	depths := make([]float32, 0, len(ref))
	if x.medianSizePerTile == 0 {
		return depths
	}
	for i, o := range ref {
		depths = append(depths, float32(float64(o)/x.medianSizePerTile))
		if depths[i] > 50000 {
			depths[i] = 50000
		}
	}
	return depths
}
const slots = 70
// with 0.5, we'll get centered at 1 and max of 2.
// so the max is 1/slotsMid
const slotsMid = float64(2) / float64(3)
func tint(f float32) int {
	if v := int(f); v < slots {
		if v < 0 {
			return 0
		}
		return v
	}
	return slots - 1
}
// CountsAtDepth calculates the count of items in depths that are at 100 * d
func CountsAtDepth(depths []float32, counts []int) {
	if len(counts) != slots {
		panic(fmt.Sprintf("indexcov: expecting counts to be length %d, was: %d", slots, len(counts)))
	}
	for _, d := range depths {
		counts[tint(d*(slots*float32(slotsMid))+0.5)]++
	}
}
// CountsROC returns a slice that indicates the cumulative proportion of
// 16KB chunks that were at least (normalized) depth given by their index.
func CountsROC(counts []int) []float32 {
	totals := make([]int, len(counts))
	totals[len(totals)-1] = counts[len(totals)-1]
	for i := len(totals) - 2; i >= 0; i-- {
		totals[i] = totals[i+1] + counts[i]
	}
	max := float32(totals[0])
	roc := make([]float32, len(counts))
	for i := 0; i < len(roc); i++ {
		roc[i] = float32(totals[i]) / max
	}
	return roc
}
func getRef(b *bam.Reader, chrom string) *sam.Reference {
	refs := b.Header().Refs()
	if strings.HasPrefix(chrom, "chr") {
		chrom = chrom[3:]
	}
	for _, ref := range refs {
		if chrom == ref.Name() {
			return ref
		}
		if strings.HasPrefix(ref.Name(), "chr") {
			if chrom == ref.Name()[3:] {
				return ref
			}
		}
	}
	return nil
}
func GetShortName(b string, isCrai bool) (string, error) {
	// TODO: replace this with samplename.Names()
	if !isCrai {
		fh, err := os.Open(b)
		if err != nil {
			return "", err
		}
		defer fh.Close()
		br, err := bam.NewReader(fh, 1)
		if err != nil {
			return "", err
		}
		defer br.Close()
		m := make(map[string]bool)
		for _, rg := range br.Header().RGs() {
			m[rg.Get(sam.Tag([2]byte{'S', 'M'}))] = true
		}
		if len(m) > 1 {
			return "", fmt.Errorf("bam reagroup: more than one RG for %s", b)
		}
		for sm := range m {
			return sm, nil
		}
	}
	vs := strings.Split(b, "/")
	v := vs[len(vs)-1]
	vs = strings.Split(v, ".")
	if len(vs) <= 2 {
		return vs[0], nil
	}
	v = strings.Join(vs[0:len(vs)-1], "-")
	return v, nil
}
func getWriter(base string) (*bgzf.Writer, error) {
	fh, err := os.Create(fmt.Sprintf("%s.bed.gz", base))
	if err != nil {
		return nil, err
	}
	w := bgzf.NewWriter(fh, 1)
	w.ModTime = time.Unix(0, 0)
	w.OS = 0xff
	return w, nil
}
func zero(ints []int) {
	for i := range ints {
		ints[i] = 0
	}
}
func getDirectory(path string) (bool, error) {
	fi, err := os.Stat(path)
	if err != nil {
		if err := os.MkdirAll(path, 0755); err != nil {
			return false, err
		}
		return true, nil
	}
	return fi.IsDir(), nil
}
// ReadFai returns a slit of references from the fai path.
// If chrom is "" all chromosomes are returned.