大家对拷贝数变异很熟悉,为了对样本进行更有意义的拷贝数变异评估,有很多学者建立了很多算法去评估拷贝数。我们今天介绍一个和拷贝数评估相关的R包CNAnorm。
首先是R包的安装:
source("https://bioconductor.org/biocLite.R")
biocLite("CNAnorm")
接下来我们看下具体的实例:
我们利用包自带的数据:
data(LS041)
LS041[1:5,]
数据的结构我们可以看出,有拷贝数的位置(Pos),测试的结果(Test),参照的结果(Norm),还有GC的含量。以上的数据可以通过处理sam/bam文件获得。使用的工具可以是bam2windows.pl脚本。
如果我们的数据是dataframe格式的数据我们需要使用R包自带的函数dataFrame2object进行处理。在此我们的数据需要进行处理:
CN <- dataFrame2object(LS041)
有些时候我们可能不需要所有染色体的GC值都计算,那么我们就可以将不需要计算的染色体过滤掉,利用下面的代码:
toSkip <- c("chrY","chrM")
CN <- gcNorm(CN, exclude = toSkip)#exclude参数就是我们要过滤的染色体名称
我们得到的上面的比率需要进一步的平滑处理,去除无用的噪声。代码如下:
CN <- addSmooth(CN, lambda = 7)#此处lambda默认是7,也不知道为啥,也许经验吧,我们就遵循默认值。
接下来就是评估拷贝数变异与否进行评估:
CN <- peakPloidy(CN, exclude = toSkip)
我们接下来对我们得到的评估结果进行可视化展示,主要展示其拷贝数的变化,以及reads的分布,代码为如下:
plotPeaks(CN, special1 = 'chrX', special2 ='chrY')
当然如果我们评估的拷贝数和实际拷贝数不相符,或者增加减少,我们也可以进行调整,代码如下:
CN <- validation(CN, ploidy =(sugg.ploidy(CN) + 1) )#假设拷贝数都多1
调整和评估的结果差异:
为了展示全基因组的拷贝数情况我们需要对以上的数据进行DNAcopy和discrete处理:
CN <- addDNACopy(CN)
CN <- discreteNorm(CN)
然后就是绘制图形:
plotGenome(CN, superimpose = 'DNACopy')
当然这个就是包的核心可视化图形了,我们可以对其进行更进一步的编辑,那就是只显示一部分染色体的分布情况:
toPlot <- c('chr10', 'chr11', 'chr12')
subSet <- chrs(CN) %in% toPlot
plotGenome(CN[subSet], superimpose ='DNACopy')
最后就是终极大招那就是,更细节的改变:
data(gPar)
gPar$genome$colors$gain.dot <-'darkorange'
gPar$genome$colors$grid <- NULL
gPar$genome$cex$gain.dot <- .4
gPar$genome$cex$loss.dot <- .4
plotGenome(CN[subSet], superimpose ='DNACopy', gPar = gPar, colorful = TRUE)
当然我们也是可以将评估的数据导出:
exportTable(CN, file ="CNAnorm_table.tab", show = 'ploidy')
导出的数据格式如下: