用R語(yǔ)言對(duì)vcf文件進(jìn)行數(shù)據(jù)挖掘.11 CNV分析

目錄

  1. 前言
  2. 方法簡(jiǎn)介
  3. 從vcf文件里提取有用信息
  4. tidy vcfR
  5. vcf可視化1
  6. vcf可視化2
  7. 測(cè)序深度覆蓋度
  8. 窗口縮放
  9. 如何單獨(dú)分離染色體
  10. 利用vcf信息判斷物種染色體倍數(shù)
  11. CNV分析

在之前的文章里介紹了如何通過直方圖來可視化等位雜合堿基的比例來判斷物種的染色體倍數(shù)性。在本文里會(huì)繼續(xù)向下挖掘，介紹如何可視化染色體上的拷貝數(shù)變化(CNVs)。

數(shù)據(jù)導(dǎo)入

和前文一樣的操作，使用包自帶的數(shù)據(jù)。

# 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
## *****        *****         *****

根據(jù)深度進(jìn)行過濾

我們需要去除過高和過低深度的數(shù)據(jù)。和前文的操作一樣，提取vcf文件里的深度數(shù)據(jù)"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%的數(shù)據(jù)，當(dāng)然此處可以根據(jù)實(shí)際情況進(jìn)行微調(diào)。然后對(duì)第一種出現(xiàn)頻率最高的堿基進(jìn)行可視化。(一般情況下一個(gè)位點(diǎn)上會(huì)有兩種堿基，具體參考前文。)

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)

同樣也可以對(duì)出現(xiàn)頻率第二高的堿基進(jìn)行同樣的操作，這里節(jié)約篇幅就省略了。

尋找峰值

為了避免復(fù)雜的基于AD比例的模型假設(shè)，程序里設(shè)計(jì)了非參數(shù)估計(jì)法來計(jì)算峰值。計(jì)算完了以后可以直接對(duì)染色體進(jìn)行拆分以后可視化進(jìn)行校驗(yàn)。

# 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)

根據(jù)尺寸把染色體分割成合適的大小

# 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函數(shù)計(jì)算峰值。并對(duì)數(shù)據(jù)進(jìn)行處理，去掉負(fù)數(shù)和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

計(jì)算到此為止，可以可視化實(shí)際數(shù)據(jù)來驗(yàn)證計(jì)算的正確性。

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))

峰值計(jì)算結(jié)果可視化

仔細(xì)想一下，峰值計(jì)算的結(jié)果其實(shí)就是CNV的結(jié)果。這里根據(jù)窗口大小把染色體分成了若干段。(那么是不是可以給每一段 CDS進(jìn)行細(xì)分然后計(jì)算出每一個(gè)CDS的具體數(shù)字呢????)

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")

當(dāng)然也可以把所有樣本組合到一起。

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)
}

全文終

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容