拷贝数的概念大家并不陌生,今天我们给大家介绍一个可以计算array-CGH平台测拷贝数的数据分析,拷贝数评估R包DNAcopy。
首先是包的安装,其还是存在于bioconductor数据库中。具体的安装请参见前面的教程。
接下来我们看下具体的评估流程:
1. 数据的获取。我们在这里实例是从包所带的例子展开,数据为脑皮质细胞的测序数据。数据结构如下:
2. 构建CAN对象方便后期的数据分析,函数为CNA。
其中的参数maploc指的探针在染色体上的位置。CGH主要的数据类型是log ratio;如果ROMA就是binary。如果染色体数据没有排序那么需要presorted设为FALSE。
CNA.object
coriell$Chromosome,coriell$Position,
data.type="logratio",sampleid=c("c05296","c13330"))
运行后如果数据中有重复的探针出现,会有个warnings,但是不影响后续的结果。
3. smooth.CNA 离散值的检测,同时对离散数据进行平滑处理。
具体的参数我们就直接取默认值了,毕竟包开发者不会把不合理的值作为默认值。
smoothed.CNA.object
4. segment 利用CBS算法进行数据的分割。循环二元分割算法(circular binary segmentation, CBS)是目前常用的SNParray 数据分段算法, 其优势在于利用相邻待测区间的数据均值差构建 t 统计量, 进而精确检测不同变异区域间的分段点。
其中verbose参数比较有趣,是为了打印计算结果的,0不产生打印;1打印样本名称;2打印染色体信息;3打印当前的分割情况。默认状态是1。Undo.splits主要是为了确定是否设置分割的阈值。后面接着给出所设阈值的大小。
segment.smoothed.CNA.object
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/
欢迎大家学习交流!
领取专属 10元无门槛券
私享最新 技术干货