Python Samfile.pileup方法代码示例

# 需要导入模块: from pysam import Samfile [as 别名]
# 或者: from pysam.Samfile import pileup [as 别名]
def callbase(bamfile, snpsites, out):
    BF = Samfile(bamfile, 'rb') #open your bam file
    SF = open(snpsites, 'r')    #the file contain snp sites info
    RF = open(out, 'w')         #resulte file
    RF.write('ref_name\tpos\tRbase\tAbase\tA\tT\tC\tG\tN\tothers\n')
    for i in SF:
        if i.startswith('#'):
            continue
        else:
            line = ParseSNPsitesLine(i)
            vcf_pos = line.pos-1 #change 1-base to 0-based
            vcf_refname = line.chrom
            print 'processing: %s %s...'%(vcf_refname, str(vcf_pos))
            At, Tt, Ct, Gt, Nt, othert = 0, 0, 0, 0, 0, 0
            for i in BF.pileup(vcf_refname, vcf_pos, vcf_pos+1):
                if i.pos == vcf_pos:
                    vcf_Rbase = line.Rbase
                    vcf_Abase = line.Abase
                    for j in i.pileups:
                        yourbase = j.alignment.seq[j.qpos]
                        if yourbase == 'A': At += 1
                        elif yourbase == 'T': Tt += 1
                        elif yourbase == 'C': Ct += 1
                        elif yourbase == 'G': Gt += 1
                        elif yourbase == 'N': Nt += 1
                        else: othert += 1
        RF.write('%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n'%(vcf_refname, \
str(vcf_pos+1), vcf_Rbase, vcf_Abase, str(At), str(Tt), str(Ct), str(Gt), \
str(Nt), str(othert)))
    BF.close()

# 需要导入模块: from pysam import Samfile [as 别名]
# 或者: from pysam.Samfile import pileup [as 别名]
 def removeEdgeMismatches(self,bamFile,minDistance, minBaseQual):
     startTime=Helper.getTime()
     minDistance=int(minDistance)
     counter=0;j=0  
     num_lines = len(self.variantDict)
     Helper.info(" [%s] remove Missmatches from the first %s bp from read edges" % (startTime.strftime("%c"),str(minDistance)),self.logFile,self.textField)
     
     bamFile = Samfile(bamFile, "rb")
     
     for varKey in self.variantDict.keys():
         variant = self.variantDict[varKey]
         
         counter+=1
         if counter%10000==0:
             Helper.status('%s mm parsed ' % counter ,self.logFile, self.textField,"grey")
         
         keepSNP=False
         varPos=variant.position-1
         iter = bamFile.pileup(variant.chromosome, variant.position-1, variant.position)
         #walks up the region wich overlap this position
         for x in iter:
             if x.pos == varPos:
                 for pileupread in x.pileups: #walk through the single reads
                     if not pileupread.is_del and not pileupread.is_refskip:
                         distance=abs(pileupread.alignment.alen-pileupread.query_position) if pileupread.alignment.is_reverse else pileupread.query_position
                         if distance >= minDistance:
                             #check readBase and Base Quality
                             if pileupread.alignment.query_sequence[pileupread.query_position] == variant.alt and pileupread.alignment.query_qualities[pileupread.query_position]>=minBaseQual:
                             #if pileupread.alignment.query_sequence[pileupread.query_position] == variant.alt:
                                 keepSNP=True
                                 
         if keepSNP==False:
             j+=1
             del self.variantDict[varKey]
     
     Helper.status('%s of %svariants were deleted' % (j,num_lines), self.logFile, self.textField,"black") 
     Helper.printTimeDiff(startTime, self.logFile, self.textField)
     bamFile.close()

# 需要导入模块: from pysam import Samfile [as 别名]
# 或者: from pysam.Samfile import pileup [as 别名]
def BamFile(bam_path):
  """Return enclosed function to read, read depths from the "bam_path".

  .. code-block:: python

    >>> from chanjo.depth_reader import BamFile
    >>> read_depths = BamFile('./alignment.bam')

  Args:
    bam_path (path): path to alignment BAM-file

  Returns:
    function: function to read from the BAM-file
  """
  # raise an error if the file doesn't exist
  if not os.path.exists(bam_path):
    raise OSError(errno.ENOENT, bam_path)

  bam = Samfile(bam_path)

  try:
    bam.pileup()
  except ValueError:
    # catch error when BAM-file isn't indexed (+ ".bai" file)
    raise OSError(
      errno.ENOENT,
      "BAM-file (%s) must be indexed." % os.path.basename(bam_path)
    )

  def reader(contig, start, end):
    """Generate a list of read depths for each position (start, end).

    The `numpy` array is used to optimize performance when building and
    slicing the list.

    This function depends on `Pysam` >=0.7.5 since the ``truncate``
    option wasn't available in previous versions.

    .. code-block:: python

      >>> read_depths = BamFile('./alignment.bam')
      >>> read_depths('17', 1, 5)
      array([3., 4., 4., 5., 4.])

    .. note::

      Positions are expected to be 1:1-based. In other words; if
      start=1, end=9 you should expect read depths for base pair
      positions 1-9 to be returned.

    Args:
      contig (str): contig/chromosome id (str) of interest
      start (int): first position of the interval (1-based)
      end (int): last position of the interval (1-based)

    Returns:
      list or numpy.array: array of read depths for *each* position in
        the interval
    """
    # convert start to 0-based since this is what pysam expects!
    pysam_start = start - 1

    # pysam expects contig as bytes in Python 2
    pysam_contig = str(contig)

    # check that we don't have a negative start position
    if pysam_start < 0:
      raise ValueError("Start position must be > 0, not %d" % start)

    # preallocate an array of 0 read depth for each position
    # pysam excludes positions with 0 read depth
    read_depths = prealloc_func(end - pysam_start)

    try:
      # overwrite read-covered positions (>0 read depth)
      # ``truncate`` ensures it starts and ends on the gives positions
      # note: ``col.pos`` is 0-based, as is ``pysam_start``
      for col in bam.pileup(pysam_contig, pysam_start, end, truncate=True):
        read_depths[col.pos - pysam_start] = col.n

    except ValueError as ve:
      # catch errors where the contig doesn't exist in the BAM-file
      raise ValueError(
        "Must use contig that exist in the Bam-file. Error: %s" % ve)

    return read_depths

  return reader