利用WeGene WGS给出的VCF文件输出类似WeGene芯片数据txt

概述：

---

编程语言：Python3.8
模块：pyvcf csv
可选：jupyter
整体思路：识别WeGene芯片数据txt的文件特征，读取vcf文件并根据其中内容获取所需数据并写入到txt中
前排提示：强烈建议买一个读写速度快一点而且至少是128GB或以上的U盘，当然我是直接买了个1T的移动硬盘

步骤：

---

1. 通过观察微基因芯片测试txt结果，我们可以得知重要信息分别为：RSID chromosome position genotype，并且我们可以发现每列之间的分隔符是一个tab缩进，所以接下来的问题就转化成如何获取这四种信息并且以分隔符为tab缩进形式输出
2. 众所周知vcf文件里的数据包括ID CHROM POS 参考序列和突变方向，而vcf中又有GT标签来表示是杂合还是纯合，是突变还是未突变。所以得出结论：我们可以直接调用pyvcf模块，并读取vcf，遍历每行的内容，同时直接用ID POS CHROM标签来获取所需的RSID chromosome position三种数据，再通过读取GT标签，原始型REF和突变型ALT来确定这个位点的genotype。注意：这几种标签的数据类型如下表

标签	类型
ID	str
CHROM	int
POS	int
ID	str
GT	str
REF	str
ALT	list

GT标签说明：示例结构'0/1'，其中0表示原始型，1表示突变

代码：

import vcf
# 读取vcf文件
vcf_data = vcf.Reader(open('input.cram', 'r'))
#设置待输出txt的表头
data = 'rsid\tchromsome\tposition\tgenotype\n'
for i in vcf_data:
	# 首先先排除没有rsid的情况
	if str(i.ID) != 'None':
		# 定义一个能够通过判断GT标签而输出基因型的函数
		def genotype(num_str):
			if num_str == '1':
				return str(i.REF)
			elif num_str == '0':
				return str(i.ALT)
		data += f'{i.ID}\t{i.CHROM}\t{i.POS}\t'
		# 获取GT标签并将DT标签的字符串以'/'为分割切割为list
		data_list = i.samples[0]['GT'].split('/')
		data += ''.join(list(map(genotype, data_list))) + '\n'

# 将数据写入txt
with open('output.txt', 'w') as f:
	f.write(data)

结果展示：

3. 从我们输出的结果可以看出还有些问题，由于ALT标签输出的是list所以文件中会存在[]字样，并且由于读取的位点中并不是所有位点的genotype长度都是2，而微基因的数据中的genotype长度都是2，所以去除[]的同时还要去除那些genotype长度非2的位点

代码：

# 去除[]
import csv
# 以\t为分割符读取刚才输出的txt，使读取数据直接变为list
with open('input.txt', 'r') as f:
	data = csv.reader(f, delimiter = '\t')
	new = []
	for row in data:
		n_row = []
		if len(row) == 4 and row[-1][-1] != 'e' and len(row[-1]) == 2 or len(row[-1]) == 4:
			if '[' in row[-1]:
				genotype = row[-1].replace('[', '').replace(']', '')
			else:
				genotype = row[-1]
			n_row.append(row[0])
			n_row.append(row[1])
			n_row.append(row[2])
			n_row.append(genotype)
			new.append('\t'.join(n_row))
		data = '\n'.join(new)
# 写
with open('output.txt', 'w') as f:
	f.write(data)

# 将genotype长度非2的位点去除
import csv
# 以\t为分割符读取刚才输出的txt，使读取数据直接变为list
with open('input.txt', 'r') as f:
	data = csv.reader(f, delimiter = '\t')
	new = []
	for row in data:
		n_row = []
		if len(row[-1]) == 2:
			n_row.append(row[0])
			n_row.append(row[1])
			n_row.append(row[2])
			n_row.append(row[-1])
			new.append('\t'.join(n_row))
		data = '\n'.join(new)
# 写
with open('output.txt', 'w') as f:
	f.write(data)