R语言实现拷贝数评估计算

拷贝数的概念大家并不陌生，今天我们给大家介绍一个可以计算array-CGH平台测拷贝数的数据分析，拷贝数评估R包DNAcopy。

首先是包的安装，其还是存在于bioconductor数据库中。具体的安装请参见前面的教程。

接下来我们看下具体的评估流程：

1. 数据的获取。我们在这里实例是从包所带的例子展开，数据为脑皮质细胞的测序数据。数据结构如下：

2. 构建CAN对象方便后期的数据分析，函数为CNA。

其中的参数maploc指的探针在染色体上的位置。CGH主要的数据类型是log ratio；如果ROMA就是binary。如果染色体数据没有排序那么需要presorted设为FALSE。

CNA.object <-CNA(cbind(coriell$Coriell.05296,coriell$Coriell.13330),

coriell$Chromosome,coriell$Position,

data.type="logratio",sampleid=c("c05296","c13330"))

运行后如果数据中有重复的探针出现，会有个warnings，但是不影响后续的结果。

3. smooth.CNA 离散值的检测，同时对离散数据进行平滑处理。

具体的参数我们就直接取默认值了，毕竟包开发者不会把不合理的值作为默认值。

smoothed.CNA.object <-smooth.CNA(CNA.object)

4. segment 利用CBS算法进行数据的分割。循环二元分割算法(circular binary segmentation, CBS)是目前常用的SNParray 数据分段算法, 其优势在于利用相邻待测区间的数据均值差构建 t 统计量, 进而精确检测不同变异区域间的分段点。

其中verbose参数比较有趣，是为了打印计算结果的，0不产生打印；1打印样本名称；2打印染色体信息；3打印当前的分割情况。默认状态是1。Undo.splits主要是为了确定是否设置分割的阈值。后面接着给出所设阈值的大小。

segment.smoothed.CNA.object <-segment(smoothed.CNA.object, verbose=1)

5. plot 可视化分割的结果。

其中主要的是参数plot.type。whole(w)指的绘制完整基因组；plateau(p)对全基因组进行排序；samplebychrom(s) 绘制每一个染色体图像;chrombysample(c)将每个染色体的所有样本展示在一起。

plot(segment.smoothed.CNA.object,plot.type="w")

plot(segment.smoothed.CNA.object, plot.type="p")

plot(segment.smoothed.CNA.object, plot.type="s")

plot(segment.smoothed.CNA.object,plot.type="c")

最后我们可以看下我们输出的结果：

head(segment.smoothed.CNA.object$output)

其中seg.mean就是我们所要的评估值。具体的分割阈值我们可以参考TCGA给出的官方文档。网址：https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/CNV_Pipeline/