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

生信技能树

发布于 2018-03-08 11:49:11

1.5K1

发布于 2018-03-08 11:49:11

文章被收录于专栏：生信技能树生信技能树

好久不见，基因组直播又来了。这篇推送是对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代码如下：

file <- 'raw-bam/GC_stat.10k.txt'
dat <- read.table(file, sep = "\t", fill=TRUE,stringsAsFactors = F)
a=dat
a$GC = a[,4]/a[,3]
a$depth = a[,5]/a[,3]
#a = a[a$depth<100,]
#a = a[a$depth>10,]
#plot(a$GC,a$depth)
chr=paste0('chr',1:22)
a=a[a[,1] %in% 1:22,]
#mean_depth = mean(a$depth,na.rm =T)
a$seg= (a$depth-157*a$GC+32)/a$depth
a$seg[a$seg<0.2 & a$seg>-0.2]=0

得到的a这个矩阵如下：

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

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

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

CNA.object <- CNA(cbind(a$seg),
a[,1],10000*(a[,2]),
data.type="logratio",sampleid="jmzeng")
CNA.object
head(as.data.frame(CNA.object))
smoothed.CNA.object <- smooth.CNA(CNA.object)
segment.smoothed.CNA.object <- segment(smoothed.CNA.object, verbose=1)
pdf('tmp1.pdf');plot(segment.smoothed.CNA.object, plot.type="w");dev.off()
pdf('tmp2.pdf');plot(segment.smoothed.CNA.object, plot.type="s") ;dev.off()
pdf('tmp3.pdf');plot(segment.smoothed.CNA.object, plot.type="p");dev.off()
sdundo.CNA.object <- segment(smoothed.CNA.object,
undo.splits="sdundo",
undo.SD=2,verbose=1)
pdf('tmp4.pdf');plot(sdundo.CNA.object,plot.type="s");dev.off()

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