func (pm *pileupMutable) addReadPair(reads []readSNP, isMinus PosType, pCtx *pileupContext) (err error) {
	for i, r := range reads {
		if err = alignRelevantBases(&(pm.alignedBaseBufs[i]), r, &pCtx.bedPart); err != nil {
			return
		}
		clipQuals(r.samr, pCtx.clip)
	}
	abb0 := pm.alignedBaseBufs[0]
	abb1 := pm.alignedBaseBufs[1]
	minBaseQual := pCtx.minBaseQual
	perReadNeeded := pCtx.perReadNeeded
	if (len(reads) == 1) || (len(abb1) == 0) {
		if len(abb0) != 0 {
			pm.addUnstitchedSegment(&(reads[0]), isMinus, abb0, minBaseQual, perReadNeeded)
			curEndMax := abb0[len(abb0)-1].posInRef + 1
			if pm.endMax < curEndMax {
				pm.endMax = curEndMax
			}
		}
		return
	}
	if len(abb0) == 0 {
		pm.addUnstitchedSegment(&(reads[1]), isMinus, abb1, minBaseQual, perReadNeeded)
		curEndMax := abb1[len(abb1)-1].posInRef + 1
		if pm.endMax < curEndMax {
			pm.endMax = curEndMax
		}
		return
	}
	mask := pm.nCirc() - 1
	seq0 := reads[0].seq8
	seq1 := reads[1].seq8
	qual0 := reads[0].samr.Qual
	qual1 := reads[1].samr.Qual
	idx0 := 0
	idx1 := 0
	for (len(abb0) != idx0) && (len(abb1) != idx1) {
		posInRef0 := abb0[idx0].posInRef
		posInRef1 := abb1[idx1].posInRef
		if posInRef0 == posInRef1 {
			row := &pm.resultRingBuffer[posInRef0&mask]
			row.Depth++
			posInRead0 := abb0[idx0].posInRead
			curSeq0 := seq0[posInRead0]  
			posInRead1 := abb1[idx1].posInRead
			curSeq1 := seq1[posInRead1]
			if curSeq0 == curSeq1 {
				base := pileup.Seq8ToEnumTable[curSeq0]
				if !perReadNeeded {
					if pCtx.qpt.lookup2(qual0[posInRead0], qual1[posInRead1]) || (base == pileup.BaseX) {
						row.Counts[base][isMinus]++
					}
				} else {
					panic("stitched per-read features not yet supported")
				}
			} else {
				row.Counts[pileup.BaseX][isMinus]++
			}
			idx0++
			idx1++
		} else if posInRef0 < posInRef1 {
			if !perReadNeeded {
				pm.addBase(posInRef0&mask, abb0[idx0].posInRead, isMinus, seq0, qual0, minBaseQual)
			} else {
				panic("stitched per-read features not yet supported")
			}
			idx0++
		} else {
			if !perReadNeeded {
				pm.addBase(posInRef1&mask, abb1[idx1].posInRead, isMinus, seq1, qual1, minBaseQual)
			} else {
				panic("stitched per-read features not yet supported")
			}
			idx1++
		}
	}
	var curEndMax PosType
	if len(abb0) != idx0 {
		pm.addUnstitchedSegment(&(reads[0]), isMinus, abb0[idx0:], minBaseQual, perReadNeeded)
		curEndMax = abb0[len(abb0)-1].posInRef + 1
	} else if len(abb1) != idx1 {
		pm.addUnstitchedSegment(&(reads[1]), isMinus, abb1[idx1:], minBaseQual, perReadNeeded)
		curEndMax = abb1[len(abb1)-1].posInRef + 1
	} else {
		curEndMax = abb0[len(abb0)-1].posInRef + 1
	}
	if pm.endMax < curEndMax {
		pm.endMax = curEndMax
	}
	return
}
addReadPair adds a single read or read-pair to the pileup.
package snp
import (
	"context"
	"fmt"
	"io/ioutil"
	"os"
	"runtime"
	"strconv"
	"github.com/grailbio/base/log"
	"github.com/grailbio/base/recordio"
	"github.com/grailbio/base/traverse"
	"github.com/grailbio/bio/circular"
	gbam "github.com/grailbio/bio/encoding/bam"
	"github.com/grailbio/bio/encoding/bamprovider"
	"github.com/grailbio/bio/encoding/fasta"
	"github.com/grailbio/bio/interval"
	"github.com/grailbio/bio/pileup"
	"github.com/grailbio/hts/sam"
)
type Opts struct {
	BedPath      string
	Region       string
	BamIndexPath string
	Clip         int
	Cols         string
	FlagExclude  int
	Mapq         int
	MaxReadLen   int
	MaxReadSpan  int
	MinBagDepth  int
	MinBaseQual  int
	Parallelism  int
	PerStrand    bool
	RemoveSq     bool
	Stitch       bool
	TempDir      string
}
var DefaultOpts = Opts{
	Clip:        0,
	FlagExclude: 0xf00,
	Mapq:        60,
	MaxReadLen:  500,
	MaxReadSpan: 511,
	MinBagDepth: 0,
	MinBaseQual: 0,
	Parallelism: 0,
	PerStrand:   false,
	RemoveSq:    false,
	Stitch:      false,
}
const (
	colBitDpRef = 1 << iota
	colBitDpAlt
	colBitEndDists
	colBitQuals
	colBitFraglens
	colBitStrands
	colBitHighQ
	colBitLowQ
)
const colPerReadMask = (colBitEndDists | colBitQuals | colBitFraglens | colBitStrands)
var colNameMap = map[string]int{
	"dpref":    colBitDpRef,
	"dpalt":    colBitDpAlt,
	"enddists": colBitEndDists,
	"quals":    colBitQuals,
	"fraglens": colBitFraglens,
	"strands":  colBitStrands,
	"highq":    colBitHighQ,
	"lowq":     colBitLowQ,
}
type refContext struct {
	refSeq8 string
	refID int
	refName string
}
type alignedPos struct {
	posInRef  PosType
	posInRead PosType
}
type pileupMutable struct {
	resultRingBuffer []PileupPayload  
	seq8Buf          []byte           
	alignedBaseBufs  [2][]alignedPos  
	firstReads       firstreadSNPTable
	endMax           PosType  
	w                recordio.Writer
	writePosScanner  interval.UnionScanner
}
func newPileupMutable(nCirc PosType, maxReadLen int, stitch bool, w *os.File) (pm pileupMutable) {
	pm = pileupMutable{
		resultRingBuffer: make([]PileupPayload, nCirc),
		seq8Buf:          make([]byte, 0, maxReadLen),
		alignedBaseBufs: [2][]alignedPos{
			make([]alignedPos, 0, maxReadLen),
			nil,
		},
	}
	if stitch {
		pm.firstReads = newFirstreadSNPTable(nCirc)
		pm.alignedBaseBufs[1] = make([]alignedPos, 0, maxReadLen)
	}
	if w != nil {
		pm.w = recordio.NewWriter(w, recordio.WriterOpts{
			Marshal:      MarshalPileupRow,
			Transformers: []string{"zstd 1"},
		})
	}
	return
}
type pileupContext struct {
	bedPart       interval.BEDUnion  
	clip          int                
	ignoreStrand  bool               
	minBaseQual   byte
	perReadNeeded bool            
	qpt           *qualPassTable  
	stitch        bool
}
func (pm *pileupMutable) addBase(circPos, posInRead, isMinus PosType, seq, qual []byte, minBaseQual byte) {
	row := &pm.resultRingBuffer[circPos]
	row.Depth++
	base := pileup.Seq8ToEnumTable[seq[posInRead]]
	if (qual[posInRead] >= minBaseQual) || (base == pileup.BaseX) {
		row.Counts[base][isMinus]++
	}
}
func (pm *pileupMutable) appendBase(circPos, posInRead, isMinus PosType, seq, qual []byte, minBaseQual, strandByte byte) {
	row := &pm.resultRingBuffer[circPos]
	row.Depth++
	base := pileup.Seq8ToEnumTable[seq[posInRead]]
	if base == pileup.BaseX {
		row.Counts[base][isMinus]++
	} else if qual[posInRead] >= minBaseQual {
		row.Counts[base][isMinus]++
		row.PerRead[base] = append(row.PerRead[base], PerReadFeatures{
			Dist5p:  uint16(posInRead),
			Fraglen: uint16(len(qual)),
			Qual:    qual[posInRead],
			Strand:  strandByte,
		})
	}
}
func (pm *pileupMutable) addUnstitchedSegment(read *readSNP, isMinus PosType, alignedBases []alignedPos, minBaseQual byte, perReadNeeded bool) {
	mask := pm.nCirc() - 1
	qual := read.samr.Qual
	if !perReadNeeded {
		for _, ab := range alignedBases {
			pm.addBase(ab.posInRef&mask, ab.posInRead, isMinus, read.seq8, qual, minBaseQual)
		}
	} else {
		strandByte := byte(pileup.GetStrand(read.samr))
		for _, ab := range alignedBases {
			pm.appendBase(ab.posInRef&mask, ab.posInRead, isMinus, read.seq8, qual, minBaseQual, strandByte)
		}
	}
}
func alignRelevantBases(result *[]alignedPos, read readSNP, bedPart *interval.BEDUnion) (err error) {
	*result = (*result)[:0]
	refID := read.samr.Ref.ID()
	posInRef := PosType(read.samr.Pos)
	endPos := read.mapEnd
	relevantIntervals := bedPart.OverlapByID(refID, posInRef, endPos)
	riIdx := 0
	nextIntervalStart := relevantIntervals[riIdx]
	nextIntervalEnd := relevantIntervals[riIdx+1]
	cigar := read.samr.Cigar
	posInRead := PosType(0)
	for _, co := range cigar {
		cLen := PosType(co.Len())
		switch co.Type() {
		case sam.CigarMatch:
			nextPosInRef := posInRef + cLen
			if nextPosInRef > nextIntervalStart {
				cigarMatchOffset := PosType(0)
				for {
					if posInRef+cigarMatchOffset < nextIntervalStart {
						cigarMatchOffset = nextIntervalStart - posInRef
					}
					if nextIntervalEnd > nextPosInRef {
						for ; cigarMatchOffset < cLen; cigarMatchOffset++ {
							(*result) = append(*result, alignedPos{
								posInRef:  posInRef + cigarMatchOffset,
								posInRead: posInRead + cigarMatchOffset,
							})
						}
						break
					}
					offsetStop := nextIntervalEnd - posInRef
					for ; cigarMatchOffset < offsetStop; cigarMatchOffset++ {
						(*result) = append(*result, alignedPos{
							posInRef:  posInRef + cigarMatchOffset,
							posInRead: posInRead + cigarMatchOffset,
						})
					}
					riIdx += 2
					if riIdx == len(relevantIntervals) {
						return
					}
					nextIntervalStart = relevantIntervals[riIdx]
					nextIntervalEnd = relevantIntervals[riIdx+1]
				}
			}
			posInRef = nextPosInRef
			posInRead += cLen
		case sam.CigarInsertion:
			posInRead += cLen
		case sam.CigarSkipped:
			fallthrough
		case sam.CigarDeletion:
			posInRef += cLen
			for posInRef >= nextIntervalEnd {
				riIdx += 2
				if riIdx == len(relevantIntervals) {
					return
				}
				nextIntervalEnd = relevantIntervals[riIdx+1]
			}
			nextIntervalStart = relevantIntervals[riIdx]
		case sam.CigarSoftClipped:
			posInRead += cLen
		case sam.CigarHardClipped:
		default:
			return fmt.Errorf("alignRelevantBases: unexpected CIGAR code %v", co)
		}
	}
	return
}