pysam - 多种格式基因组数据（sam/bam/vcf/bcf/cram/…）读写与处理模块（python）
2023年02月04日
现在又回头来处理fasta，读取可以用pysam，可以很省力，输出就用python，因为结构很简单。
import pysam
all_fasta = pysam.FastaFile("all.human.protein.fasta")
all_fasta.references[0]
all_fasta.fetch(all_fasta.references[0])
for tmp_protein_name in all_fasta.references:
tmp_protein_seq = all_fasta.fetch(tmp_protein_name)
# just use protein ID as file name, good for shell process
tmp_fasta_name = "individual/" + tmp_protein_name
# safe open and close file
with open(tmp_fasta_name, mode="w") as fout:
fout.write('>' + tmp_protein_name + '\n')
fout.write(tmp_protein_seq + '\n')
# break
# fout = open(tmp_fasta_name, mode='w')
# fout.close()
# https://www.cnblogs.com/leezx/p/5908767.html
# https://pysam.readthedocs.io/en/latest/api.html?highlight=fasta#fasta-files
# https://ucdavis-bioinformatics-training.github.io/2022-Feb-Introduction-To-Python-For-Bioinformatics/python/python5
在开发基因组相关流程或工具时，经常需要读取、处理和创建bam、vcf、bcf文件。目前已经有一些主流的处理此类格式文件的工具，如samtools、picard、vcftools、bcftools，但此类工具集成的大多是标准功能，在编程时如果直接调用的话往往显得不够灵活。
本文介绍的是一个处理基因组数据的python模块，它打包了htslib-1.3、samtools-1.3 和 bcftools-1.3的核心功能，能在编程时非常灵活的处理bam和bcf文件。
以下主要介绍pysam的安装和使用方法：
1. 安装
如果Linux上安装了pip，可以一键安装，在集群上的话，需要登录安装节点进行安装。
pip3 install pysam
检查是否安装成功
import pysam
2.读取bam文件(pysam.AlignmentFile)
bam是sam的二进制文件，因其占用空间少，所以都会使用bam进行存储和操作。
要读取bam文件,必须先创建一个AlignmentFile对象.
path_in = './test.bam' samfile = pysam.AlignmentFile(path_in, "rb")
之后就可以逐行读取和处理bam文件了(顺序读取),以下打印出了bam的一行.
for line in samfile: print(line) break
但顺序读取还不够灵活,我们有时需要随机读取(提示:sam不能随机读取),pysam的fetch方法提供了随机读取功能.
直接使用fetch会报错
ValueError: fetch called on bamfile without index
提示我们需要建立(.bai)索引
samtools index corrected.bam
fetch返回的是一个迭代器(iterator),可以迭代读取内容.
for read in samfile.fetch('chr6', 28478220, 28478222): ... print(read)
fetch方法的API如下,chr6为参考序列,后面数字分别为读取的起始和终止位置.
fetch(self, reference=None, start=None, end=None, region=None, tid=None, until_eof=False, multiple_iterators=False)
3.读取vcf/bcf文件(pysam.VariantFile)
读取方法同上,只是使用的是VariantFile方法:
gvcf = "./MHC.unified.g.vcf.gz" vcf_in = pysam.VariantFile(gvcf)
若想随机读取,仍然需要建立索引:
首先使用bgzip压缩vcf
bgzip -c MHC.g.vcf > MHC.g.vcf.gz
然后用bcftools建立索引
bcftools index -c MHC.g.vcf.gz
使用fetch读取
for rec in vcf_in.fetch('chr6', 28577796, 28577896): ... print(rec) ... break
3.随机读取fasta文件（faidx建立索引）
读取方法略有不同，fetch返回的本身就是一个字符串。
samtools faidx total_PacBio_reads.fasta
fasta_file = pysam.FastaFile(path) fasta_file.fetch("m160727_060737_42266_c101014182550000001823222610211695_s1_p0/110008/22268_22731")
4.创建并写入到新的bam或vcf文件
pysam的核心功能是可以随心所欲的读取数据,处理之后,写入到一个新建的bam或bcf文件里.
我们可以完全自定义一些内容,然后写入到一个新的bam文件里,如下:
header = { 'HD': {'VN': '1.0'},
'SQ': [{'LN': 1575, 'SN': 'chr1'},
{'LN': 1584, 'SN': 'chr2'}] }
with pysam.AlignmentFile(tmpfilename, "wb", header=header) as outf:
a = pysam.AlignedSegment()
a.query_name = "read_28833_29006_6945"
a.query_sequence="AGCTTAGCTAGCTACCTATATCTTGGTCTTGGCCG"
a.flag = 99
a.reference_id = 0
a.reference_start = 32
a.mapping_quality = 20
a.cigar = ((0,10), (2,1), (0,25))
a.next_reference_id = 0
a.next_reference_start=199
a.template_length=167
a.query_qualities = pysam.qualitystring_to_array("<<<<<<<<<<<<<<<<<<<<<:<9/,&,22;;<<<")
a.tags = (("NM", 1),
("RG", "L1"))
outf.write(a)
同理,我们也可以读取一个已有的bam文件,逐个修改以上的属性,然后存储到一个新的bam文件里.这里不再举例.
上面设置header可能有点麻烦,容易出错,但我们可以复制一个已有bam文件的header到一个新的bam文件里.
outf = pysam.AlignmentFile(path_out, "wb", template=samfile)
以上template参数指定了模板bam文件.
5. 关闭文件
outf.close()
总结:
pysam模块非常实用,有了pysam模块,我们就可以非常灵活的操纵bam/bcf文件,而不必依赖于samtools或bcftools. pysam可以随机读取bam/bcf文件,也可以将处理后的内容自定义输出到bam/bcf文件.
以上只介绍了pysam最常见的功能,更多pysam功能请参照:http://pysam.readthedocs.io/en/latest/index.html
