findGtNodes.nim

## Genrate GtNode per variant (SNP)

import tables
import hts
import utils

# minSnpDepth = 2, this SNP should at least be covered by two cells
proc findGtNodes*(rec:Variant, variantIndex: int, ibam:Bam, maxTotalReads:int, minTotalReads:int, mapq: int,
                  barcodeTable:TableRef,
                  minbsq:int,
                  barcodeTag:string ): Table[string,GtNode] =

  var barcodedNodesTable = initTable[string,GtNode]()
  var rec_alt:char
  var total_reads = 0
  rec_alt = rec.ALT[0][0]
  for aln in ibam.query(chrom = $rec.CHROM,start = rec.POS.cint-1, stop = rec.POS.cint):
    var cbt = tag[string](aln, barcodeTag)
    var currentCB: string
    var base_off: int
    var base: char
    if cbt.isNone:
      continue  
    else:
      currentCB = cbt.get      
    if aln.flag.unmapped or aln.mapping_quality.cint < mapq or aln.flag.dup:
#      echo "read skipped, bad mapping"
      continue
    ## skip unmapped, duplicated, mapq low reads or aln.flag.secondary or aln.flag.supplementary
    ## if not aln.flag.proper_pair: continue
    if not barcodeTable.hasKey(currentCB): 
      continue
    base_off = get_base_offset(position = rec.POS.cint, align = aln) 
    base = aln.base_at(base_off)   
    if aln.base_quality_at(base_off).cint < minbsq: 
      continue
    total_reads+=1
    if barcodedNodesTable.hasKey(currentCB):
      if base == rec.REF[0]: 
    #     echo "adding ref geno to cell " & currentCB
        barcodedNodesTable[currentCB].alleles = add_allele(barcodedNodesTable[currentCB],alleleBinGeno = 0)
        continue
      if base == rec_alt: 
        barcodedNodesTable[currentCB].alleles = add_allele(barcodedNodesTable[currentCB],alleleBinGeno = 1)
        continue
    else:
    #  echo "creating new node for cell " & currentCB
      var newGtNode = GtNode(chrom: $rec.CHROM, variantIndex: variantIndex)
      if base == rec.REF[0]: 
        newGtNode.alleles = "0"
        barcodedNodesTable[currentCB] = newGtNode
        continue
      if base == rec_alt: 
        newGtNode.alleles = "1"
        barcodedNodesTable[currentCB] = newGtNode
        continue
  if total_reads > maxTotalReads or total_reads < minTotalReads:
    return initTable[string,GtNode]()
  return barcodedNodesTable