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

用户7010445

发布于 2020-03-03 14:53:44

1.5K0

发布于 2020-03-03 14:53:44

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

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

~/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脚本

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()

使用方法

python3 convert_vcf_to_vennplot.py 5_sample.vcf

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

WS-12.txt  WS-17.txt  WS-2.txt  WS-4.txt  WS-5.txt

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

韦恩图R代码

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

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

简单用法

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")