目录
在之前的文章里介绍了如何通过直方图来可视化等位杂合碱基的比例来判断物种的染色体倍数性。在本文里会继续向下挖掘，介绍如何可视化染色体上的拷贝数变化(CNVs)。
数据导入
和前文一样的操作，使用包自带的数据。
# Load libraries
library(vcfR)
library(pinfsc50)
# Determine file locations
vcf_file <- system.file("extdata", "pinf_sc50.vcf.gz",
package = "pinfsc50")
# Read data into memory
vcf <- read.vcfR(vcf_file, verbose = FALSE)
vcf
## ***** Object of Class vcfR *****
## 18 samples
## 1 CHROMs
## 22,031 variants
## Object size: 20.9 Mb
## 7.929 percent missing data
## ***** ***** *****
根据深度进行过滤
我们需要去除过高和过低深度的数据。和前文的操作一样，提取vcf文件里的深度数据"AD"。
ad <- extract.gt(vcf, element = 'AD')
allele1 <- masplit(ad, record = 1)
allele2 <- masplit(ad, record = 2)
ad1 <- allele1 / (allele1 + allele2)
ad2 <- allele2 / (allele1 + allele2)
然后过滤出10%~90%的数据，当然此处可以根据实际情况进行微调。然后对第一种出现频率最高的碱基进行可视化。(一般情况下一个位点上会有两种碱基，具体参考前文。)
dp <- allele1
#sums <- apply(dp, MARGIN=2, quantile, probs=c(0.15, 0.95), na.rm=TRUE)
sums <- apply(dp, MARGIN=2, quantile, probs=c(0.1, 0.9), na.rm=TRUE)
par(mfrow=c(4,3))
par(mar=c(2,2,1,1))
par(oma=c(1,1,0,0))
for(i in 1:12){
hist(allele1[,i], breaks = seq(0,1e3,by=1), xlim=c(0,100), col=8, main="", xlab="", ylab="")
title(main = colnames(allele1)[I])
abline(v=sums[,i], col=2)
}
title(xlab = "Depth", line=0, outer = TRUE, font=2)
title(ylab = "Count", line=0, outer = TRUE, font=2)
同样也可以对出现频率第二高的碱基进行同样的操作，这里节约篇幅就省略了。
寻找峰值
为了避免复杂的基于AD比例的模型假设，程序里设计了非参数估计法来计算峰值。计算完了以后可以直接对染色体进行拆分以后可视化进行校验。
# Filter on depth quantiles.
sums <- apply(allele1, MARGIN=2, quantile, probs=c(0.1, 0.9), na.rm=TRUE)
# Allele 1
dp2 <- sweep(allele1, MARGIN=2, FUN = "-", sums[1,])
#allele1[dp2 < 0] <- NA
vcf@gt[,-1][ dp2 < 0 & !is.na(vcf@gt[,-1]) ] <- NA
dp2 <- sweep(allele1, MARGIN=2, FUN = "-", sums[2,])
#allele1[dp2 > 0] <- NA
vcf@gt[,-1][dp2 > 0] <- NA
# Allele 2
dp2 <- sweep(allele2, MARGIN=2, FUN = "-", sums[1,])
vcf@gt[,-1][ dp2 < 0 & !is.na(vcf@gt[,-1]) ] <- NA
dp2 <- sweep(allele2, MARGIN=2, FUN = "-", sums[2,])
vcf@gt[,-1][dp2 > 0] <- NA
# Censor homozygotes.
gt <- extract.gt(vcf, element = 'GT')
hets <- is_het(gt)
is.na( vcf@gt[,-1][ !hets ] ) <- TRUE
# Extract allele depths
ad <- extract.gt(vcf, element = 'AD')
allele1 <- masplit(ad, record = 1)
allele2 <- masplit(ad, record = 2)
ad1 <- allele1 / (allele1 + allele2)
ad2 <- allele2 / (allele1 + allele2)
根据尺寸把染色体分割成合适的大小
# Parameters
#winsize <- 1e5
#
winsize <- 2e5
#bin_width <- 0.1
#bin_width <- 0.05
#bin_width <- 0.025
#
bin_width <- 0.02
#bin_width <- 0.01
然后用freq_peak函数计算峰值。并对数据进行处理，去掉负数和Na值。
# Find peaks
freq1 <- ad1/(ad1+ad2)
freq2 <- ad2/(ad1+ad2)
myPeaks1 <- freq_peak(freq1, getPOS(vcf), winsize = winsize, bin_width = bin_width)
#myCounts1 <- freq_peak(freq1, getPOS(vcf), winsize = winsize, bin_width = bin_width, count = TRUE)
is.na(myPeaks1$peaks[myPeaks1$counts < 20]) <- TRUE
myPeaks2 <- freq_peak(freq2, getPOS(vcf), winsize = winsize, bin_width = bin_width, lhs = FALSE)
#myCounts2 <- freq_peak(freq2, getPOS(vcf), winsize = winsize, bin_width = bin_width, count = TRUE)
is.na(myPeaks2$peaks[myPeaks2$counts < 20]) <- TRUE
计算到此为止，可以可视化实际数据来验证计算的正确性。
par(mfrow=c(4,4))
par(mar=c(2,2,1,1))
par(oma=c(1,1,0,0))
mySample <- "BL2009P4_us23"
for(i in 1:4){
hist(freq1[ myPeaks1$wins[i,'START_row']:myPeaks1$wins[i,'END_row'], mySample ],
breaks = seq(0,1,by=bin_width), xlim=c(0,1), col=8, main = "", xaxt='n')
axis(side=1, at=c(0,0.25,0.333,0.5,0.666,0.75,1),
labels=c(0,'1/4','1/3','1/2','2/3','3/4',1), las=1)
abline(v=myPeaks1$peaks[i,mySample], col=2)
if(i==2){ title(main=mySample) }
}
mySample <- "DDR7602"
for(i in 1:4){
hist(freq1[ myPeaks1$wins[i,'START_row']:myPeaks1$wins[i,'END_row'], mySample ],
breaks = seq(0,1,by=bin_width), xlim=c(0,1), col=8, main = "", xaxt='n')
axis(side=1, at=c(0,0.25,0.333,0.5,0.666,0.75,1),
labels=c(0,'1/4','1/3','1/2','2/3','3/4',1), las=1)
abline(v=myPeaks1$peaks[i,mySample], col=2)
if(i==2){ title(main=mySample) }
}
mySample <- "IN2009T1_us22"
for(i in 1:4){
hist(freq1[ myPeaks1$wins[i,'START_row']:myPeaks1$wins[i,'END_row'], mySample ],
breaks = seq(0,1,by=bin_width), xlim=c(0,1), col=8, main = "", xaxt='n')
axis(side=1, at=c(0,0.25,0.333,0.5,0.666,0.75,1),
labels=c(0,'1/4','1/3','1/2','2/3','3/4',1), las=1)
abline(v=myPeaks1$peaks[i,mySample], col=2)
if(i==2){ title(main=mySample) }
}
mySample <- "P17777us22"
for(i in 1:4){
hist(freq1[ myPeaks1$wins[i,'START_row']:myPeaks1$wins[i,'END_row'], mySample ],
breaks = seq(0,1,by=bin_width), xlim=c(0,1), col=8, main = "", xaxt='n')
axis(side=1, at=c(0,0.25,0.333,0.5,0.666,0.75,1),
labels=c(0,'1/4','1/3','1/2','2/3','3/4',1), las=1)
abline(v=myPeaks1$peaks[i,mySample], col=2)
if(i==2){ title(main=mySample) }
}
par(mfrow=c(1,1))
par(mar=c(5,4,4,2))
par(oma=c(0,0,0,0))
峰值计算结果可视化
仔细想一下，峰值计算的结果其实就是CNV的结果。这里根据窗口大小把染色体分成了若干段。(那么是不是可以给每一段 CDS进行细分然后计算出每一个CDS的具体数字呢????)
i <- 2
layout(matrix(1:2, nrow=1), widths = c(4,1))
par(mar=c(5,4,4,0))
mySample <- colnames(freq1)[I]
plot(getPOS(vcf), freq1[,mySample], ylim=c(0,1), type="n", yaxt='n',
main = mySample, xlab = "POS", ylab = "Allele balance")
axis(side=2, at=c(0,0.25,0.333,0.5,0.666,0.75,1),
labels=c(0,'1/4','1/3','1/2','2/3','3/4',1), las=1)
abline(h=c(0.25,0.333,0.5,0.666,0.75), col=8)
points(getPOS(vcf), freq1[,mySample], pch = 20, col= "#A6CEE344")
points(getPOS(vcf), freq2[,mySample], pch = 20, col= "#1F78B444")
segments(x0=myPeaks1$wins[,'START_pos'], y0=myPeaks1$peaks[,mySample],
x1=myPeaks1$wins[,'END_pos'], lwd=3)
segments(x0=myPeaks1$wins[,'START_pos'], y0=myPeaks2$peaks[,mySample],
x1=myPeaks1$wins[,'END_pos'], lwd=3)
bp1 <- hist(freq1[,mySample], breaks = seq(0,1,by=bin_width), plot = FALSE)
bp2 <- hist(freq2[,mySample], breaks = seq(0,1,by=bin_width), plot = FALSE)
par(mar=c(5,1,4,2))
barplot(height=bp1$counts, width=0.02, space = 0, horiz = T, add = FALSE, col="#A6CEE3")
barplot(height=bp2$counts, width=0.02, space = 0, horiz = T, add = TRUE, col="#1F78B4")
当然也可以把所有样本组合到一起。
par(mfrow=c(4,3))
par(mar=c(2,2,1,1))
par(oma=c(1,1,0,0))
for(i in 1:12){
mySample <- colnames(freq1)[I]
plot(getPOS(vcf), freq1[,mySample], ylim=c(0,1), type="n", yaxt='n',
main = mySample, xlab = "POS", ylab = "Allele balance")
axis(side=2, at=c(0,0.25,0.333,0.5,0.666,0.75,1),
labels=c(0,'1/4','1/3','1/2','2/3','3/4',1), las=1)
abline(h=c(0.25,0.333,0.5,0.666,0.75), col=8)
points(getPOS(vcf), freq1[,mySample], pch = 20, col= "#A6CEE344")
points(getPOS(vcf), freq2[,mySample], pch = 20, col= "#1F78B444")
segments(x0=myPeaks1$wins[,'START_pos'], y0=myPeaks1$peaks[,mySample],
x1=myPeaks1$wins[,'END_pos'], lwd=3)
segments(x0=myPeaks1$wins[,'START_pos'], y0=myPeaks2$peaks[,mySample],
x1=myPeaks1$wins[,'END_pos'], lwd=3)
}
全文终
