【发布时间】：2016-11-10 13:15:27
【问题描述】：
对于初学者，我是生物信息学的新手，尤其是编程新手，但我已经构建了一个脚本，该脚本将通过所谓的 VCF 文件（仅包括个人，一个 clumn = 一个个人），并使用搜索字符串找出每个变体（系）个体是纯合子还是杂合子。
这个脚本至少在小的子集上有效，但我知道它将所有内容都存储在内存中。我想在非常大的压缩文件（甚至整个基因组）上执行此操作，但我不知道如何将此脚本转换为逐行执行所有操作的脚本（因为我想计算整个列我只是不看看如何解决）。
所以输出是每个个体 5 个东西（总变体、纯合子数量、杂合子数量以及纯合子和杂合子的比例）。请看下面的代码：
#!usr/bin/env python
import re
import gzip
subset_cols = 'subset_cols_chr18.vcf.gz'
#nuc_div = 'nuc_div_chr18.txt'
gz_infile = gzip.GzipFile(subset_cols, "r")
#gz_outfile = gzip.GzipFile(nuc_div, "w")
# make a dictionary of the header line for easy retrieval of elements later on
headers = gz_infile.readline().rstrip().split('\t')
print headers
column_dict = {}
for header in headers:
column_dict[header] = []
for line in gz_infile:
columns = line.rstrip().split('\t')
for i in range(len(columns)):
c_header=headers[i]
column_dict[c_header].append(columns[i])
#print column_dict
for key in column_dict:
number_homozygotes = 0
number_heterozygotes = 0
for values in column_dict[key]:
SearchStr = '(\d)/(\d):\d+,\d+:\d+:\d+:\d+,\d+,\d+'
#this search string contains the regexp (this regexp was tested)
Result = re.search(SearchStr,values)
if Result:
#here, it will skip the missing genoytypes ./.
variant_one = int(Result.group(1))
variant_two = int(Result.group(2))
if variant_one == 0 and variant_two == 0:
continue
elif variant_one == variant_two:
#count +1 in case variant one and two are equal (so 0/0, 1/1, etc.)
number_homozygotes += 1
elif variant_one != variant_two:
#count +1 in case variant one is not equal to variant two (so 1/0, 0/1, etc.)
number_heterozygotes += 1
print "%s homozygotes %s" % (number_homozygotes, key)
print "%s heterozygotes %s" % (number_heterozygotes,key)
variants = number_homozygotes + number_heterozygotes
print "%s variants" % variants
prop_homozygotes = (1.0*number_homozygotes/variants)*100
prop_heterozygotes = (1.0*number_heterozygotes/variants)*100
print "%s %% homozygous %s" % (prop_homozygotes, key)
print "%s %% heterozygous %s" % (prop_heterozygotes, key)
任何帮助将不胜感激，因此我可以继续调查大型数据集， 谢谢你:)
顺便说一下，VCF 文件看起来像这样： 个人_1 个人_2 个人_3 0/0:9,0:9:24:0,24,221 1/0:5,4:9:25:25,0,26 1/1:0,13:13:33:347,33,0
然后是带有个人 ID 名称的标题行（我总共有 33 个具有更复杂 ID 标签的个人，我在这里简化了），然后我有很多具有相同特定模式的这些信息行。我只对带有斜线的第一部分感兴趣，因此是常规表达式。
【问题讨论】：
-
如果您可以编辑问题以包含 VCF 文件顶部的示例以及预期结果，这将有所帮助。
-
你用的是什么版本的 Python？
-
基本上你想要的是逐步解压缩一个csv文件并产生行
-
VCF 文件的示例会很有帮助。 TIA。
-
我在 VCF 文件中添加了简化的数据，但是通过查看下面的代码，我发现您已经自己弄清楚了 :) 通常每行也有常规信息（也在列中，在个人列之前） ，但我会在应用此脚本之前将它们过滤掉。
标签： python csv gzip bioinformatics vcf-variant-call-format