用R语言对vcf文件进行数据挖掘.6 vcf可视化2

目录

1. 前言
2. 方法简介
3. 从vcf文件里提取有用信息
4. tidy vcfR
5. vcf可视化1
6. vcf可视化2
7. 测序深度覆盖度
8. 窗口缩放
9. 如何单独分离染色体
10. 利用vcf信息判断物种染色体倍数
11. CNV分析

在上一篇文章vcf可视化1里，主要介绍了针对数据全体summary的可视化方法。这篇文章里会介绍针对个别样本的可视化。

使用数据

和之前一样，我们使用自带数据包

library(vcfR)

vcf_file <- system.file("extdata", "pinf_sc50.vcf.gz", package = "pinfsc50")
dna_file <- system.file("extdata", "pinf_sc50.fasta", package = "pinfsc50")
gff_file <- system.file("extdata", "pinf_sc50.gff", package = "pinfsc50")

vcf <- read.vcfR(vcf_file, verbose = FALSE)
dna <- ape::read.dna(dna_file, format = "fasta")
gff <- read.table(gff_file, sep="\t", quote="")

chrom <- create.chromR(name="Supercontig", vcf=vcf, seq=dna, ann=gff, verbose=FALSE)
chrom <- masker(chrom, min_DP = 300, max_DP = 700)
chrom <- proc.chromR(chrom, verbose = FALSE)

Genotype数据

同样在之前的文章里介绍了如何用extract.gt()提取Genotype数据。可以通过同样的方法提取genotype数据的深度(DP)，然后对其进行可视化。

dp <- extract.gt(chrom, element="DP", as.numeric=TRUE)
rownames(dp) <- 1:nrow(dp)
head(dp)

##   BL2009P4_us23 DDR7602 IN2009T1_us22 LBUS5 NL07434 P10127 P10650 P11633
## 1             7       6             8     6      12      6      4      6
## 2            12      20            16    20      28      9      8     11
## 3            27      26            23    26      39     22      8     11
## 4            29      27            32    27      38     22      7     10
## 5            26      30            41    30      44     18     11     11
## 6            23      36            35    36      54     18     20     18
##   P12204 P13527 P1362 P13626 P17777us22 P6096 P7722 RS2009P1_us8 blue13
## 1      1      7    NA     13         14     6    NA            6     16
## 2      6     29     1     41         33    11    NA           21     31
## 3      6     47     3     58         58    11    NA           28     71
## 4      4     44     5     70         62    10    NA           35     63
## 5     NA     29     3     62         49    11    NA           38     60
## 6      1     37     4     48         52    18    NA           27     68
##   t30-4
## 1     2
## 2     3
## 3     2
## 4     2
## 5     4
## 6     7

我们可以清楚的看到每一个位点的测序深度，可以对这些数据进行可视化，看一下大概的深度变化。从而估计拷贝数的变化。
用heatmap.bp()来上一个热图来可视化测序深度。

heatmap.bp(dp[1001:1500,])

棒状图分别代表行或者列的和。右边的渐变色图例里，黄色代表高数量的高质量碱基，蓝色代表低数量的高质量碱基。如果指定的碱基长度过长，就会覆盖和忽略掉很多细节，所以建议热图的碱基长度不要超过1000。
当然，虽然热图的显示有限制只能显示1000bp以内的部分区域，但是可以单独显示全局棒状图。

par(mar=c(8,4,4,2))
barplot(apply(dp, MARGIN=2, mean, na.rm=TRUE), las=3)

同理，是不是可以按照这个办法按照染色体来显示拷贝数(CNVs)变化呢，之后需要来尝试一下。