csmatch.py 文件源码

def crossGenotypeWindows(commonSNPsCHR, commonSNPsPOS, snpsP1, snpsP2, inFile, binLen, outFile, logDebug = True):
    ## inFile are the SNPs of the sample
    (snpCHR, snpPOS, snpGT, snpWEI, DPmean) = snpmatch.parseInput(inFile = inFile, logDebug = logDebug)
    # identifying the segregating SNPs between the accessions
    # only selecting 0 or 1
    segSNPsind = np.where((snpsP1 != snpsP2) & (snpsP1 >= 0) & (snpsP2 >= 0) & (snpsP1 < 2) & (snpsP2 < 2))[0]
    log.info("number of segregating snps between parents: %s", len(segSNPsind))
    (ChrBins, PosBins) = getBinsSNPs(commonSNPsCHR, commonSNPsPOS, binLen)
    log.info("number of bins: %s", len(ChrBins))
    outfile = open(outFile, 'w')
    for i in range(len(PosBins)):
      start = np.sum(PosBins[0:i])
      end = start + PosBins[i]
      # first snp positions which are segregating and are in this window
      reqPOSind = segSNPsind[np.where((segSNPsind < end) & (segSNPsind >= start))[0]]
      reqPOS = commonSNPsPOS[reqPOSind]
      perchrTarPosind = np.where(snpCHR == ChrBins[i])[0]
      perchrTarPos = snpPOS[perchrTarPosind]
      matchedAccInd = reqPOSind[np.where(np.in1d(reqPOS, perchrTarPos))[0]]
      matchedTarInd = perchrTarPosind[np.where(np.in1d(perchrTarPos, reqPOS))[0]]
      matchedTarGTs = snpGT[matchedTarInd]
      try:
        TarGTBinary = snpmatch.parseGT(matchedTarGTs)
        TarGTBinary[np.where(TarGTBinary == 2)[0]] = 4
        genP1 = np.subtract(TarGTBinary, snpsP1[matchedAccInd])
        genP1no = len(np.where(genP1 == 0)[0])
        (geno, pval) = getWindowGenotype(genP1no, len(genP1))
        outfile.write("%s\t%s\t%s\t%s\t%s\n" % (i+1, genP1no, len(genP1), geno, pval))
      except:
        outfile.write("%s\tNA\tNA\tNA\tNA\n" % (i+1))
      if i % 40 == 0:
        log.info("progress: %s windows", i+10)
    log.info("done!")
    outfile.close()