使用pysam操作VCF/BCF文件

简介: 使用pysam操作VCF/BCF文件读取和写出from pysam import VariantFilebcf_in = VariantFile("test_in.

使用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)

目录
相关文章
|
3月前
|
Java
File常用的方法操作、在磁盘上创建File、获取指定目录下的所有文件、File文件的重命名、将数据写入File文件
这篇文章介绍了Java中`File`类的常用方法操作,包括文件和目录路径名的抽象表示、构造方法、常用方法以及代码实例。文章通过具体的代码示例,展示了如何在磁盘上创建文件和目录、处理目录存在的情况、文件重命名、获取指定目录下的所有文件以及将数据写入文件等操作。每个代码实例都配有测试结果截图,帮助读者更好地理解和掌握`File`类的应用。
File常用的方法操作、在磁盘上创建File、获取指定目录下的所有文件、File文件的重命名、将数据写入File文件
|
6月前
讲解:如何根据txt文本列出的文件名批量查找指定文件夹里的文件并复制到新的文件夹里 , 按照文件名批量复制文件 , 根据文件名批量复制 , 通过文件名批量复制文件
该文介绍了一款批量文件处理软件,用于解决三类问题:依据文件名清单批量复制图片、筛选PDF文件及删除指定文件。用户可从百度或腾讯下载链接获取软件。操作步骤包括加载文件夹、输入文件名清单、设置目标位置、选择操作类型(复制、剪切或删除)及匹配方式,然后开始查找。软件能快速处理大量文件,提高办公效率,避免手动操作。查找结果会显示在特定文件夹和日志中,记录每个操作详情。
657 5
|
6月前
|
Windows
推荐:如何批量根据PDF文件名批量查找PDF文件,复制到指定地方保存,通过文件名批量复制文件,按照文件名批量复制文件,根据文件名批量提取文件
该文介绍了一个批量查找PDF文件(不限于找PDF)的工具,用于在多级文件夹中快速查找并复制特定文件。通过下载提供的软件,用户可以加载PDF库,输入文件名列表,设置操作参数(如保存路径、复制或删除)及搜索模式。软件能高效执行,例如在1.1秒内完成对数千文件中的37个目标文件的复制,显著提升了工作效率,避免了手动逐个查找和复制的繁琐。
610 0
|
Web App开发 存储 数据可视化
VCF文件修改染色体格式的方法
VCF文件修改染色体格式的方法
ENVI_IDL:如何读取HDF4文件和HDF5文件的数据集和属性?
参考Documentation Center (l3harrisgeospatial.com) 由于时间和精力,这里并没有对其中的函数进行详细的说明,只是介绍了函数作用以及用法。
229 0
使用FILE结构操作文本文件
使用FILE结构操作文本文件
109 0
|
存储 编解码
ENVI_IDL:读取OMI数据(HDF5)并输出为Geotiff文件+详细解析
ENVI_IDL:读取OMI数据(HDF5)并输出为Geotiff文件+详细解析
412 1
ENVI_IDL: 文本文件的读取(主要是txt、csv文件)
ENVI_IDL: 文本文件的读取(主要是txt、csv文件)
451 0
|
存储
ENVI_IDL: 创建HDF5文件并写入数据(以将Geotiff文件写入HDF文件为例) + 详细解析
ENVI_IDL: 创建HDF5文件并写入数据(以将Geotiff文件写入HDF文件为例) + 详细解析
191 0
|
数据处理 索引 Python
Bioinfo|bedtools-操作VCF文件
Bioinfo|bedtools-操作VCF文件
261 0