使用pysam操作VCF/BCF文件
读取和写出
from pysam import VariantFile
bcf_in = VariantFile("test_in.vcf", "r")
bcf_out = VariantFile("test_out.vcf", "w", header=bcf_in.header)
for rec in bcf_in.fecth():
bcf_out.write(rec)
VariantFile函数得到的是 pysam.libcbcf.VariantFile 对象, 这是一个可遍历对象, 通过dir()可以发现它有__iter__和__next__方法。因此如果仅仅是遍历全部记录，那么__iter__等价于fecth.
type(bcf_in) # 对象类型
dir(bcf_out) # 方法
VCF格式分为Header和Record两个部分. record记录每个变异位点的具体信息，为了从中提取所需数据，需要理解Pysam的解析策略。
rec1 = bcf_in.__next__()
dir(rec1)
vcf的record每一行都是9列+N列样本(CHROM, POS, ID, REF, ALT, QUAL, FILTER, INFO, FORMAT, sample1, sample2,..), 解析之后就是如下方法
- .chrom: 返回字符串
- .pos: 返回数值。 这个是以0为基, 可以用.start和.stop
- .id: 如果无记录, 就是NoneType
- .ref: 返回字符串
- .alt: 返回元祖(tuple), 因为一个位点上可以有多个变异类型
- .qual: 返回数值
- .filter: 返回pysam.libcbcf.VariantRecordFilter对象, 类似于字典
- .info: 返回pysam.libcbcf.VariantRecordInfo对象，类似于字典, 存放所有样本的统计信息
- .format: 返回pysam.libcbcf.VariantRecordFormat，类似于字典, 存放后续每个样本数据存放顺序和数据类型
- .samples: 返回pysam.libcbcf.VariantRecordSamples, 类似于字典, 存放每一个样本的具体信息
.filter, .info, .format, .samples虽然都能返回类字典(或者说哈希表)数据结果，但是在方法上存在差别。
VariantRecordFilter对象可以通过.filter.add增加过滤类型, 当然需要事先在header中添加元信息，如下:
bcf_in.header.filters.add(id="ugly",number=None, type=None,description="i don't likt it") #增加员信息
rec = bcf_in.__next__()
rec.filter.add("ugly") # 增加过滤条件
rec.filter.keys() # 查看
VariantRecordInfo对象可以删除一个键值对(pop),可以更新已有的键值对。
rec.info.pop('TYPE') # 删除TYPE
rec.info['ODDS'] # 变更前
rec.info.update({'ODDS':12}) #变更
rec.info['ODDS'] # 变更后
VariantRecordFormat和VariantRecordSamples关系比较紧密，但前者只能查看不提供方法进行修改, 而VariantRecordSamples和VariantRecordInfo一致。由于可以有多个样本，提取数据的时候就需要多层迭代，例如提取所有样本的GT
for key,value in rec.samples.iteritems():
print(key, value['GT'])
例如只有两个样本，我想比较这两个样本的GT是否相同
GT = [value['GT'] for value in rec.samples.values()]
GT[0].__eq__(GT[-1])
综上，就可以在Python中写出一个过滤器剔除缺失基因组记录，保留其中样本基因组纯合但不同的记录
import sys
from pysam import VariantFile as vcf
if len(sys.argv) < 3:
sys.exit(1)
else:
in_name = sys.argv[1]
out_name = sys.argv[2]
bcf_in = vcf(in_name)
# add metadata
command = "##pysamCommand=GT[0].__ne__((None,)) and GT[-1].__ne__((None,)) and GT[0].__ne__(GT[-1]) and GT[0].__ne__((0,1)) and GT[-1].__ne__((0,1))"
bcf_in.header.add_line(command)
bcf_out = vcf(out_name, "w", header=bcf_in.header)
for rec in bcf_in.__iter__():
GT = [value['GT'] for value in rec.samples.values()]
if GT[0].__ne__((None,)) and GT[-1].__ne__((None,)) and \
GT[0].__ne__((0,1)) and GT[-1].__ne__((0,1)) and \
GT[0].__ne__(GT[-1]):
bcf_out.write(rec)
