R语言实现拷贝数变异评估预测
大家对拷贝数变异很熟悉，为了对样本进行更有意义的拷贝数变异评估，有很多学者建立了很多算法去评估拷贝数。我们今天介绍一个和拷贝数评估相关的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')
导出的数据格式如下：
