安装

BiocManager::install("vcfR")

测试数据

该软件自带测试数据，可进行操练

# 读取文件
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)
vcf <- read.vcfR( vcf_file, verbose = FALSE )
dna <- ape::read.dna(dna_file, format = "fasta")
gff <- read.table(gff_file, sep="\t", quote="")

读取的vcf为一个 $\color{red}{S4}$ 对象

> vcf
***** Object of Class vcfR *****
18 samples
1 CHROMs
22,031 variants
Object size: 22.4 Mb
7.929 percent missing data
*****        *****         *****

## 包含如下3个元素
meta
character vector for the meta information
fix
matrix for the fixed information
gt
matrix for the genotype information

3个元素包含的信息可以简单查阅

# meta 头信息
> head(vcf@meta)
[1] "##fileformat=VCFv4.1"                                                                                                            
[2] "##source=\"GATK haplotype Caller, phased with beagle4\""                                                                         
[3] "##FILTER=<ID=LowQual,Description=\"Low quality\">"                                                                               
[4] "##FORMAT=<ID=AD,Number=.,Type=Integer,Description=\"Allelic depths for the ref and alt alleles in the order listed\">"           

# fix 1-9列的信息
> head(vcf@fix[,1:7])
     CHROM              POS   ID REF  ALT QUAL      FILTER
[1,] "Supercontig_1.50" "41"  NA "AT" "A" "4784.43" NA    
[2,] "Supercontig_1.50" "136" NA "A"  "C" "550.27"  NA    
[3,] "Supercontig_1.50" "254" NA "T"  "G" "774.44"  NA    
[4,] "Supercontig_1.50" "275" NA "A"  "G" "714.53"  NA    

# gt 各个样本的基因型信息
 > head(vcf@gt[,1:3])
     FORMAT           BL2009P4_us23              DDR7602                    
[1,] "GT:AD:DP:GQ:PL" "1|1:0,7:7:21:283,21,0"    "1|1:0,6:6:18:243,18,0"    
[2,] "GT:AD:DP:GQ:PL" "0|0:12,0:12:36:0,36,427"  "0|0:20,0:20:60:0,60,819"  
[3,] "GT:AD:DP:GQ:PL" "0|0:27,0:27:81:0,81,1117" "0|0:26,0:26:78:0,78,1077" 
[4,] "GT:AD:DP:GQ:PL" "0|0:29,0:29:87:0,87,1243" "0|0:27,0:27:81:0,81,1158"

既然已经知道怎么调取相应信息，则可将其存为数据框进行任意操作

# 将fix内容存为 info2
info2 <- vcf@fix
info2 <- as.data.frame(info2)

# 调取QUAL值存为qa
qa <- info2$QUAL
# 查看qual >30 的个数
qa <- as.numeric(as.character(qa))
table(qa > 30)
 TRUE 
22031

$\color{red}{将因子转变为数值时现将其转变为字符，然后在转变为数值}$

建立chromR

chrom <- create.chromR(name='Supercontig', vcf=vcf, seq=dna, ann=gff)
# 画图
plot(chrom)

该图可以看出，深度分布非常突出，峰值可以表示测序深度的一个区域，也能看出有很长的尾巴，表示存在变异。同时MQ值都在60左右达到了峰值，所以可以根据改值进行过滤。QUAL不太好看出时从那里开始逐渐到峰值

还可以根据 $\color{red}{masker()}$ 进行过滤

chrom <- masker(chrom, min_QUAL = 1, min_DP = 300, max_DP = 700, min_MQ = 59.9,  max_MQ = 60.1)
plot(chrom)

一旦我们对要考虑的高质量变体感到满意，就可以使用 $\color{red}{proc.chromR（）}$ 处理chromR对象

chrom <- proc.chromR(chrom, verbose=TRUE)
plot(chrom)

刚才我们已经得到了变异文件，fasta文件，gff文件，并且过滤掉我们不满意的数据，现在我们来对其进行可视化
使用 $\color{red}{chromoqc()}$

chromoqc(chrom, dp.alpha=20)

也可以进行调整X坐标，筛选出我们想看的区域

chromoqc(chrom, xlim=c(5e+05, 6e+05))

最后导出数据可以使用 $\color{red}{write.vcf()}$

交流欢迎扫码

最后编辑于：2020.06.12 15:01:29