对SNP的VCF文件质控筛选
1.vcfR简介
vcfR is a package intended to help visualize, manipulate and quality filter data in VCF file.
参考资料链接
2.安装vcfR
直接通过install.packages()函数就可以安装
install.packages("vcfR")
3.vcf格式文件
VCF文件分为两部分内容：前面以“#”开头的注释部分（非重点），和后面整齐规则的没有“#”开头的突变信息部分。
3.1 注释信息
注释主要包括一些软件参数和文件信息，例如vcf注释信息1vcf注释信息2 可以看出，注释中标出了样本名，筛选参数（Hard—Fliter），筛选用的个软件以及对应的genome信息文件，分析路径，snpeff注释文件的一些信息。
vcf突变信息1
vcf突变信息2
突变信息一般是10列，包括
- CHROM ： 染色体信息
- POS ：发生变异的位置的第一个碱基所在的位置
- ID ：variant位点对应dbSNP数据库中的ID，若没有，则默认使用‘.’
- REF ： genome.fa该位置的碱基类型及碱基数量
- ALT ：测序得到结果中，所支持的碱基类型及碱基数量；对于SNP来说是单个碱基类型的编号，而对于Indel来说是指碱基个数的添加或缺失，以及碱基类型的变化，即variant碱基类型
- QUAL ： variant的质量。Phred格式的数值，代表着此位点是纯合的概率，此值越大，则概率越低，代表着次位点是variants的可能性越大（表示变异碱基的可能性）
- FILTER ： 次位点是否要被过滤掉。如果是PASS，则表示此位点可以考虑为variant。
- INFO ： variant的相关信息
- FORMAT ： variants的格式，例如GT:AD:DP:GQ:PL
- Sample：由BAM文件中的@RG下的SM标签所决定，这些值对应着第9列的各个格式，不同格式的值用冒号分开，每一个sample对应着1列；多个samples则对应着多列，这种情况下列的数大于10列。
4.vcf文件第8列信息
第八列信息是通过注释软件对variant进行注释，增加的基因信息和位置信息，主要关注三部分信息
DP ：样本在这个位置的reads覆盖度。是一些reads被过滤掉后的覆盖度，DP4:高质量测序碱基，位于REF或者ALT前后
QD：通过深度来评估一个变异的可信度
ANN：基因的变化情况，突变碱基的位置，基因，蛋白变化情况，可信度等级
5.vcf文件的基因型信息
VCF文件第9列是基因型信息的多个标签，各个标签的意义展示如下：
- 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一致。- GQ：即第二可能的基因型的PL值，相对于最可能基因型的PL值（其PL=0）而言，大于99时，其信息量已不大，因此大于99的全部赋值99。当GQ值很小时，意味着第二可能基因型与最可能基因型差别不大。
- GL：三种基因型（RR RA AA）出现的可能性，R表示参考碱基，A表示变异碱基
- DV：高质量的非参考碱基
- SP：phred的p值误差线
- PL：指定的三种基因型的可能性(provieds the likelihoods of the given genotypes)。这三种指定的基因型为(0/0,0/1,1/1)，这三种基因型的概率总和为1。和之前不一致，该值越大，表明为该种基因型的可能性越小。
Phred值 = -10 * log (p) p为基因型存在的概率。
6.数据处理分析可视化
6.1 质控
library(vcfR)
data(vcfR_example)
chrom <- create.chromR('sc50', seq=dna, vcf=vcf, ann=gff)
head(chrom)
chrom
plot(chrom)
vcf质控1
chrom <- masker(chrom, min_QUAL = 1, min_DP = 300, max_DP = 700, min_MQ = 59, max_MQ = 61)
chrom <- proc.chromR(chrom, win.size=1000)
chrom <- proc.chromR(chrom, verbose=TRUE)
plot(chrom)
vcf质控2
vcf质控3
vcf质控4
6.2 Allele分析
data(vcfR_example)
gt <- extract.gt(vcf)
hets <- is_het(gt)
# Censor non-heterozygous positions.
is.na(vcf@gt[,-1][!hets]) <- TRUE
# Extract allele depths.
ad <- extract.gt(vcf, element = "AD")
ad1 <- masplit(ad, record = 1)
ad2 <- masplit(ad, record = 2)
freq1 <- ad1/(ad1+ad2)
freq2 <- ad2/(ad1+ad2)
myPeaks1 <- freq_peak(freq1, getPOS(vcf))
is.na(myPeaks1$peaks[myPeaks1$counts < 20]) <- TRUE
myPeaks2 <- freq_peak(freq2, getPOS(vcf), lhs = FALSE)
is.na(myPeaks2$peaks[myPeaks2$counts < 20]) <- TRUE
freq_peak_plot(pos = getPOS(vcf), ab1 = freq1, ab2 = freq2, fp1 = myPeaks1, fp2=myPeaks2)
等位基因分布
6.3 各样本序列情况
library(ape)
data(vcfR_test)
# Create an example reference sequence.
nucs <- c('a','c','g','t')
set.seed(9)
myRef <- as.DNAbin(matrix(nucs[round(runif(n=50, min=0.5, max=4.5))], nrow=1))
# Recode the POS data for a smaller example.
set.seed(99)
vcfR_test@fix[,'POS'] <- sort(sample(10:20, size=length(getPOS(vcfR_test))))
# Just vcfR
myDNA <- vcfR2DNAbin(vcfR_test)
seg.sites(myDNA)
image(myDNA)
# ref.seq, no start.pos
myDNA <- vcfR2DNAbin(vcfR_test, ref.seq = myRef)
seg.sites(myDNA)
image(myDNA)
# ref.seq, start.pos = 4.
# Note that this is the same as the previous example but the variants are shifted.
myDNA <- vcfR2DNAbin(vcfR_test, ref.seq = myRef, start.pos = 4)
seg.sites(myDNA)
image(myDNA)
myDNA <- vcfR2DNAbin(vcfR_test, unphased_as_NA = FALSE, ref.seq = myRef)
seg.sites(myDNA)
image(myDNA)
样本序列情况
还有很多其他功能，SNP的热图，位点pop富集等。。。
