R语言实现VCF文件的处理可视化
基因突变数据大家应该很熟悉，作为突变信息的存储文件VCF文件，记录了突变的位点以及对应的突变信息。文件分为三个部分 ‘#’号开头行——meta， 非#号开头行分为fix和gt两个部分。fix部分存储vcf文件中非#号开头行的前7列，分别是染色体编号、碱基位置、ID、参考碱基、变异碱基、质量值、是否过滤；gt 部分存储两部分内容format、样本基因型。今天给大家介绍下在R语言中处理vcf文件的包vcfR。首先看下包的安装：
install.packages("vcfR")
install.packages('adegenet')
install.packages('poppr')接下来通过实例来看下具体的操作：
###实例数据位置
pkg <-"pinfsc50"
vcf_file <-system.file("extdata", "pinf_sc50.vcf.gz", package = pkg)
dna_file <-system.file("extdata", "pinf_sc50.fasta", package = pkg)
gff_file <-system.file("extdata", "pinf_sc50.gff", package = pkg)##数据读入
library(vcfR)
library(poppr)
library('adegenet')
library(ape)
vcf <-read.vcfR( vcf_file, verbose = FALSE )
dna <-ape::read.dna(dna_file, format = "fasta")
gff <-read.table(gff_file, sep="\t", quote="")###创建数据对象，其中dna和ann主要是注释vcf文件的信息
chrom <-create.chromR(name='Supercontig', vcf=vcf, seq=dna, ann=gff)###可视化对象
plot(chrom)
图中，Read depth（DP）测序深度（reads数）指的是不同位置频率的密度分布，从图中来看每个基因组的大部分都是在某个倍体水平进行测序的。在这里我们看到了一个峰值，这可能代表了那个基倍体区域，但我们也看到了一个长尾，这可能代表了拷贝数变异。从MQ图可以看出映射质量(MQ)都在60左右达到峰值。这个值会根据不同的方法有所差异。
##低质量数据的筛选
chrom <-masker(chrom, min_QUAL = 1, min_DP = 300, max_DP = 700, min_MQ = 59.9, max_MQ = 60.1)
plot(chrom)
##获取突变信息，右下角则为对应每个位点突变的频数
chrom <-proc.chromR(chrom, verbose=TRUE)
plot(chrom)
##多信息合并展示
chromoqc(chrom,dp.alpha=20)
##放大局部区域
chromoqc(chrom,xlim=c(5e+05, 6e+05))
VCF文件中基因型数据包括：
GT：样品的基因型（genotype）。两个数字中间用’/'分开，这两个数字表示双倍体的sample的基因型。0 表示样品中有ref的allele；1 表示样品中variant的allele；2表示有第二个variant的allele。因此：0/0 表示sample中该位点为纯合的，和ref一致；0/1 表示sample中该位点为杂合的，有ref和variant两个基因型；1/1 表示sample中该位点为纯合的，和variant一致。
AD 和 DP：AD(Allele Depth)为sample中每一种allele的reads覆盖度,在diploid（二倍体）中则是用逗号分割的两个值，前者对应ref基因型，后者对应variant基因型；DP（Depth）为sample中该位点的覆盖度。
GQ：基因型的质量值(Genotype Quality)。Phred格式(Phred_scaled)的质量值，表示在该位点该基因型存在的可能性；该值越高，则Genotype的可能性越大；计算方法：Phred值 = -10 * log (1-p) p为基因型存在的概率。
PL：指定的三种基因型的质量值(provieds the likelihoods of the given genotypes)。这三种指定的基因型为(0/0,0/1,1/1)，这三种基因型的概率总和为1。和之前不一致，该值越大，表明为该种基因型的可能性越小。Phred值 = -10 * log (p) p为基因型存在的概率。
##获取基因型数据中的测序深度
dp <-extract.gt(chrom, element = "DP", as.numeric=TRUE)
boxplot(dp,col=2:8, las=3)
##热图绘制
heatmap.bp(dp[1001:1500,])
##缺失信息筛选
myMiss <-apply(dp, MARGIN = 1, function(x){ sum( is.na(x) ) } )
myMiss <-myMiss / ncol(dp)
vcf <-vcf[myMiss < 0.2, ]##导出vcf文件
write.vcf( chrom,file = "vcfR_test.vcf.gz" )##基因型数据转化
gt <-extract.gt(chrom)
gt[gt=="0|0"]<- 0
gt[gt=="0/0"]<- 0
gt[gt=="1/1"]<- 1
gt[gt=="2/2"]<- 2
head(gt)
##数值化基因型
myHQ1 <-masplit(gt[,1:2], sort = 0)
##各样本序列信息
myDNA <-vcfR2DNAbin(chrom, unphased_as_NA = FALSE, ref.seq = dna)
seg.sites(myDNA)
image(myDNA)
##导出序列数据
write.dna( myDNA,file = 'my_gene.fasta', format = 'fasta' )欢迎大家互相学习！