文件读写

文件读写句柄

bam_file : bam文件路径

读取文件

1	samfile = pysam.AlignmentFile(bam_file) # 读取文件句柄

写入文件

WriteBam = pysam.AlignmentFile(bam_file , "wb", header=pysamFile.header)
for read in samfile.fetch():
    if read.is_paired:
        WriteBam.write(read)

pairedreads.close()
samfile.close()

Bam操作

提取特定区域的Reads

# 获取染色体chr1 上 100~120bp的reads
for read in samfile.fetch('chr1', 100, 120):
    print(read)

samfile.close()

pileup(self, contig=None, start=None, stop=None, region=None, reference=None, end=None, **kwargs)

统计特定区域的碱基分布

1
2
3

>>> base = bam.count_coverage('chr15', 43997074,43997076,quality_threshold=0)

(array('L', [1, 0]), array('L', [0, 143]), array('L', [0, 0]), array('L', [141, 0]))

返回一个4维的属组，分别对应 A(base[0][site_index]) C(base[1][site_index]) G(base[2][site_index]) T(base[3][site_index]) ,
注意：pysam 提供的参数 quality_threshold 只针对 base_quality 进行过滤，不支持 map_quality 的过滤

文件函数

# region : 提取的区域
# mapq  : 比对质量
# baseq ：碱基质量
# ref   :参考基因组



for pileup_column in bam_file.pileup(region="chr1:1-10",mapq=mapq , baseq = baseq,
                                     stepper=stepper, fastaFile=ref, max_depth=200000, **{"truncate": True}):
#
        for pileup_read in pileup_column.pileups:
            aln = pileup_read.alignment
            read_name = aln.query_name
            pair = 'pe1' if aln.is_read1 else 'pe2'
            strand = '-' if aln.is_reverse else '+'
            read = Read(read_name, pair, strand)
            if pileup_read.is_del or pileup_read.is_refskip or (aln.flag > 1024) or (aln.mapping_quality < mapq) or \
                    aln.query_qualities[pileup_read.query_position] < baseq:
                continue
            start_reads[read] = [pileup_read.query_position, aln]

Reference

readthedoc