yhlhhhhh yhlhhhhh - anthrogenica主页:Sansui Yang CSDN:yhlhhhhh blogger:yhlgeneticslab 全基因组测序

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

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

WechatIMG24.jpeg

 
步骤:

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表示突变
代码:

截屏2021-07-21_下午9.30_.34_.png

 
结果展示:

截屏2021-07-19_上午9.41_.37_.png


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

截屏2021-07-21_下午9.32_.59_.png

 

截屏2021-07-21_下午9.33_.10_.png

 
2021-07-21 • 北京, 北京市, 中国
按热门排序    按默认排序

3 个回复

大灰狼 - Don't worry, Be happy~
用plink --vcf in.vcf --recode 23 就行了
啊,我想搞魔方的
鬼王 - 팬지가 피다
越来越深奥了+++

要回复问题请先登录注册