Pysam pileup的使用_hustlxc

个人资料

微博

正文字体大小：大中小

Pysam pileup的使用

(2012-06-12 15:17:43)

标签：

分类： Pysam

1. 通过pysam包的pileup(...)方法来提取某个指定位置的pileup信息

# 导入pysam包，找不到pysam就退出

import sys

try:

import pysam

except Exception as e:

print >> sys.stderr, str(e)

sys.exit(-1)

hg="/human_hg19_Lilly/human_hg19_Lilly.fasta" # 已经通过BWA建立好index的人类基因组参考序列

samfile = pysam.Samfile("H318-T_80X/H318-T_80X-calibrated.bam", "rb"）# 打开bam文件，须建好index

fastafile = pysam.Fastafile(hg) # 创建Fastafile对象，方便随机检索基因组上某个位置的序列

cov = 0 # 覆盖度初始值

chr = 'chr1' # 染色体

pos = 631656 # 1-based

refbase = fastafile.fetch(chr, pos-1, pos) # 0-based index，pysam的index是从0开始，可以使用samtools region的表示方式，即region="chr1:631656-631656"（下标从1开始）

for pc in samfile.pileup(chr, pos-1, pos, stepper="all", fastafile=None): # 0-based index

if pc.pos != pos - 1: continue # 只提取位置631656（下标从1开始）的pileup信息，pysam的下标从0开始,须减1
depth = pc.n # Depth at this position

for pr in pc.pileups: # pr为PileupRead对象

cov += 1 # 该位置覆盖度加1
print pr.indel # indel length at this position

#print pr.alignment.seq[pr.qpos], pr.qpos, pr.is_tail, pr.level, pr.is_del, pr.indel, pr.alignment.qual[pr.qpos]

print pr.alignment.seq[pr.qpos], pr.qpos, pr.alignment.qual[pr.qpos] # 输出该位置reads的碱基，碱基位置和该碱基的质量值

print k, refbase # 输出质量值和位置631656的碱基序列

Further reading
http://thebird.nl/biolib/Adding_BioLib_BAM_SAM_Support.html

阅读┊ 收藏 ┊ 喜欢 ▼ ┊打印┊举报/Report