使用pysam操作BAM文件

最新推荐文章于 2024-08-13 08:47:56 发布

生信修炼手册

最新推荐文章于 2024-08-13 08:47:56 发布

阅读量2.7k

点赞数 1

文章标签：编程语言 python 数据分析封装 eclipse

欢迎关注”生信修炼手册”!

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·

—如果喜欢，快分享给你的朋友们吧—

原创不易，欢迎收藏，点赞，转发！生信知识浩瀚如海，在生信学习的道路上，让我们一起并肩作战！

本公众号深耕耘生信领域多年，具有丰富的数据分析经验，致力于提供真正有价值的数据分析服务，擅长个性化分析，欢迎有需要的老师和同学前来咨询。