使用pysam操作BAM文件
pysam模块对samtools和tabix进行了封装，可以在python程序内部来操作和访问相关的文件，具体地，支持以下4种文件
1. Fasta/Fastq
2. VCF
3. Tabix file
4. BAM/CRAM/SAM
对于samtools的封装，提现在操作bam文件上，既可以通过编程来读取bam文件中的内容，也可以实现samtools的调用；对tabix的封装，体现在利用索引来提取对应区域的record。
1. Fasta和Fastq
Fasta和Fastq，也常称为fastx格式，对于读取而言，pysam提供了以下接口
>>> with pysam.FastxFile(fastq) as f:
... for entry in f:
... print(entry.name)
... print(entry.sequence)
... print(entry.comment)
... print(entry.quality)该接口同时适用于fastq和fasta,只不过对于fasta而言，没有quality和comment属性。
对于有fai索引的fasta文件，还可以通过fetch函数来提取对应region的碱基，此时的读取方式如下
>>> import pysam
>>> fasta = pysam.FastaFile('input.fasta')
>>> fasta.references
['chr1', 'chr2', 'chr3', 'chr4', 'chr5']
>>> fasta.filename
'input.fasta'
>>> fasta.nreferences
5
>>> fasta.lengths
[248956422, 242193529, 198295559, 190214555, 181538259]
>>> fasta.get_reference_length('chr1')
248956422通过对应的属性，可以方便获取染色体名称，长度，个数等属性，fetch通过指定染色体，起始和终止位置来定义region,用法如下
>>> region = fasta.fetch('chr1', 20000, 20100)
>>> region
'CCTGGTGCTCCCACAAAGGAGAAGGGCTGATCACTCAAAGTTGCGAACACCAAGCTCAACAATGAGCCCTGGAAAATTTCTGGAATGGATTATTAAACAG'2. VCF
对于VCF文件，可以通过如下方式遍历所有的记录
>>> vcf = pysam.VariantFile('dbsnp.hg19.vcf')
>>> for i in vcf:
... print(i.chrom)
... print(i.id)
... print(i.ref)对于vcf的头文件，可以通过如下方式访问
>>> vcf.header
<pysam.libcbcf.VariantHeader object at 0x7f4cf923f190>
>>> print(vcf.header)
##fileformat=VCFv4.1
##FILTER=<ID=PASS,Description="All filters passed">
>>> list(vcf.header.contigs)
['chrM', 'chr1', 'chr2', 'chr3', 'chr4', 'chr5']
>>> list(vcf.header.filters)
['PASS', 'NC']
>>> list(vcf.header.info)
['ASP', 'ASS', 'CAF', 'CDA', 'CFL', 'CLNACC', 'CLNALLE', 'CLNDBN', 'CLNDSDB', 'CLNDSDBID', 'CLNHGVS', 'CLNORIGIN', 'CLNSIG', 'CLNSRC', 'CLNSRCID', 'COMMON', 'DSS', 'G5', 'G5A', 'GENEINFO', 'GNO', 'HD', 'INT', 'KGPROD', 'KGPhase1', 'KGPilot123', 'KGValidated', 'LSD', 'MTP', 'MUT', 'NOC', 'NOV', 'NSF', 'NSM', 'NSN', 'OM', 'OTH', 'OTHERKG', 'PH3', 'PM', 'PMC', 'R3', 'R5', 'REF', 'RS', 'RSPOS', 'RV', 'S3D', 'SAO', 'SLO', 'SSR', 'SYN', 'TPA', 'U3', 'U5', 'VC', 'VLD', 'VP', 'WGT', 'WTD', 'dbSNPBuildID']对于有tbi索引的vcf文件，还可以通过fetch函数来访问特定region的记录，用法如下
>>> for i in vcf.fetch('chr1', 1000, 2000):
... print(i.id)3. Tabix
tabix支持对bed, gff, bam, vcf等多种文件建立索引，这里的Tabix的意思是专指对于bed, gff这两种纯文本格式的文件的处理，主要功能是使用fetch来提取对应region的记录，用法如下
>>> bed = pysam.TabixFile("hg19.bed")
>>> for row in bed.fetch('chr1', 1000, 2000):
... print(str(row))4. BAM
对于Bam文件，遍历行的操作如下
>>> bam = pysam.AlignmentFile('input.bam')
>>> for i in bam:
... print(i.qname)
... print(i.flag)同时，还可以通过fetch和pileup两种方式来访问，fetch访问区域的alignment，用法如下
>>> for i in bam.fetch('chr1', 10000, 20000):
... print(i.qname)
... print(i.flag)pileup访问基因组每个碱基的比对情况，用法如下
>>> for i in bam.pileup('chr1', 10000, 20000):
... print(i.reference_name)
... print(i.pos)
... print(i.nsegments)从形式上看，两种方法差不多，但是其返回值为不同Class的对象，可以根据API来访问具体的属性和方法。
除了访问操作，也可以调用samtools的功能，因为pysam是对samtools的封装，所以samtools的子命令在该模块中，可以通过函数形式来调用，用法如下
>>> print(pysam.view.usage())
Usage: samtools view [options] <in.bam>|<in.sam>|<in.cram> [region ...]
Options:
-b output BAM
-C output CRAM (requires -T)
-1 use fast BAM compression (implies -b)
-u uncompressed BAM output (implies -b)
-h include header in SAM output
-H print SAM header only (no alignments)
------
>>> pysam.view('-o', 'out.bam', 'accepted_hits.bam')如果需要对上述几种文件根据指定区域提取子集，或者针对bam文件进行更加个性化的统计处理，可以使用pysam来实现，集成到python开发环境中，实现更加复杂的逻辑处理，会更加的高效。
·end·
—如果喜欢，快分享给你的朋友们吧—
腾讯云开发者
