【直播】我的基因组59：CNV初步探索

好久不见，基因组直播又来了。这篇推送是对SNV进行一个初步探索。

单纯的一个样本来找CNV，总是不太准确的，但还是那句话，毕竟是自己的基因组，硬着头皮也要上。当然，分析的结果，我是不会拿来预测健康风险什么的，但是可以一步步的往前推，学习就是这样，慢慢来。

搜索一些CNV的简单资料放在这里吧

参考文献：

Statistical models for DNA copy number variation detection using
read-depth data from next generation sequencing experiments

好了，言归正传，我第一次分析CNV基于全基因组分窗口滑动的测序深度以及GC含量。

我在这里选择了一个bioconductor的包来做，叫做DNAcopy，

http://bioconductor.org/packages/release/bioc/html/DNAcopy.html

说明书非常通俗易懂，就是接收每个探针对应区域的染色体号，探针坐标，以及该探针检测到的信号值。

那么我的全基因组分窗口滑动的测序深度经过GC含量矫正之后与标准测序深度的偏差，就是信号值咯。

我处理数据的R代码如下：

  1. file <- 'raw-bam/GC_stat.10k.txt'

  2. dat <- read.table(file, sep = "\t", fill=TRUE,stringsAsFactors = F)

  3. a=dat

  4. a$GC = a[,4]/a[,3]

  5. a$depth = a[,5]/a[,3]

  6. #a = a[a$depth<100,]

  7. #a = a[a$depth>10,]

  8. #plot(a$GC,a$depth)

  9. chr=paste0('chr',1:22)

  10. a=a[a[,1] %in% 1:22,]

  11. #mean_depth = mean(a$depth,na.rm =T)

  12. a$seg= (a$depth-157*a$GC+32)/a$depth

  13. a$seg[a$seg<0.2 & a$seg>-0.2]=0

得到的a这个矩阵如下：

每一行是一个探针，第一列是染色体号，第二列是窗口的顺序编号，第3列是该窗口被测到的碱基数量，第4列是该窗口含有的GC碱基数量，第5列是该窗口所有碱基的测序深度总和。

因为我不是很明白GC含量跟测序深度的矫正关系，我把0.2以下的信号值全部归零。

这个数据就可以导入到DNAcopy这个R包了，它需要构建一个CNA.object对象，代码如下：

  1. CNA.object <- CNA(cbind(a$seg),

  2. a[,1],10000*(a[,2]),

  3. data.type="logratio",sampleid="jmzeng")

  4. CNA.object

  5. head(as.data.frame(CNA.object))

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

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

  8. pdf('tmp1.pdf');plot(segment.smoothed.CNA.object, plot.type="w");dev.off()

  9. pdf('tmp2.pdf');plot(segment.smoothed.CNA.object, plot.type="s") ;dev.off()

  10. pdf('tmp3.pdf');plot(segment.smoothed.CNA.object, plot.type="p");dev.off()

  11. sdundo.CNA.object <- segment(smoothed.CNA.object,

  12. undo.splits="sdundo",

  13. undo.SD=2,verbose=1)

  14. pdf('tmp4.pdf');plot(sdundo.CNA.object,plot.type="s");dev.off()    

因为隐私的问题，我只秀其中的一张图给大家看看，而且我不能把具体的CNV文本文件结果给大家看到。

可以看到我的X染色体有一个拷贝的完全缺失，因为我是男性，只有一条X染色体。

然后我大部分的染色体都是正常的2倍体，之所以中间的红线不是一条，而是在0.0附近，是因为我给信号值的时候简简单单的把0.2以内的归零，可能不够，我还需要在学习，今天就学到这里吧。

如果大家实在感兴趣这个CNV分析，可以直接去运行R包DNAcopy的测试数据即可：

测试的代码如下：

http://bioconductor.org/packages/release/bioc/vignettes/DNAcopy/inst/doc/DNAcopy.R

文：Jimmy

图文编辑：吃瓜群众

赞 (0)

相关推荐