多样本vcf文件转换成R语言韦恩图输入格式

2020-03-03 14:53:44 浏览数 (1)

基因组重测序的论文中有些可能会用韦恩图来展示不同样本snp的交集和差异。那么如何将手头的vcf文件转换成R语言里做韦恩图要求的数据格式呢？想了几天有了一些想法，记录在这里。

从总vcf文件中提取出5个样本的信息重新组成一个vcf文件

代码语言：javascript复制
~/mvcf-subset --exclude-ref -c WS-2,WS-4,WS-5,WS-12,WS-17 412_all_cp.recode.eva.vcf > 5_sample.vcf

利用python脚本将数据转化为R语言里做韦恩图要求的格式

python脚本的基本原理就是判断样本的基因型，如果是0/0，则这个样本在这个位点不是变异，如果不是0/0，则在这个位点存在变异。又想到一种情况是如果这个位点在某个样本里是缺失数据如何处理，这个后续需要 考虑进去。可能还需要考虑的一个问题是把snp和indel分开。

  • python脚本
代码语言：javascript复制
import vcf
import sys


input_vcf = sys.argv[1]


records = vcf.Reader(filename=input_vcf)


vcf_dict = {}

for sample in records.samples:
    vcf_dict[sample] = []

i = 0
for record in records:
    i = i   1
    for sam in record.samples:
        if sam["GT"] == "0/0":
            vcf_dict[sam.sample].append(0)
        else:
            vcf_dict[sam.sample].append(i)

for sample,var_site in vcf_dict.items():
    fw = open(sample ".txt","w")
    fw.write(sample "n")
    for i in var_site:
        if i != 0:
            fw.write(str(i) "n")
    fw.close()
  • 使用方法
代码语言：javascript复制
python3 convert_vcf_to_vennplot.py 5_sample.vcf

脚本借助了pyvcf模块，如果没有这个模块需要使用pip命令来安装 在当前目录下就会多出来几个文件

代码语言：javascript复制
WS-12.txt  WS-17.txt  WS-2.txt  WS-4.txt  WS-5.txt

自己在win10上用这个脚本的时候会提示我连接不到某个动态库，但是也能获得结果，搜索了一下暂时还没有找到原因，在linux系统下没有遇到任何报错，如果自己的vcf文件比较大，运行的时间可能会长一点。

韦恩图R代码

参考 如何使用R来绘制韦恩图（Venn Diagram）

代码语言：javascript复制
setwd("../../vcf_handling/Rice_CP/")
A<-read.table("WS-2.txt",header=T)
B<-read.table("WS-4.txt",header=T)
C<-read.table("WS-5.txt",header=T)
D<-read.table("WS-12.txt",header=T)
E<-read.table("WS-17.txt",header=T)
library(VennDiagram)
help(package="VennDiagram")
venn.diagram(x=list(A=A$WS.2,
                    B=B$WS.4,
                    C=C$WS.5,
                    D=D$WS.12,
                    E=E$WS.17),
             filename="My5.png",
             fill =c("cornflowerblue","green","yellow","darkorchid1","red"),
             )

最近又发现一个新的R语言包用来做韦恩图 VennDetail

github 主页 https://github.com/guokai8/VennDetail

  • 简单用法
代码语言：javascript复制
install.packages("devtools")
devtools::install_github("guokai8/VennDetail")
help(package="VennDetail")
library(VennDetail)

res <- venndetail(list(A = A$WS.2, 
                       B = B$WS.4, 
                       C = C$WS.5,
                       D = D$WS.12,
                       E = E$WS.17))
plot(res, type = "venn")

其他功能有需要的时候再来探索！

本文中用到的vcf格式文件大家可以在论文中找到下载链接https://www.jianshu.com/p/f6b72450f589。

python r语言 https github

0 人点赞