搜索
查看: 2360|回复: 0

【直播】我的基因组60:CNV初步探索

[复制链接]

103

主题

133

帖子

860

积分

版主

Rank: 7Rank: 7Rank: 7

积分
860
发表于 2017-3-10 21:53:06 | 显示全部楼层 |阅读模式
本帖最后由 zckoo007 于 2017-3-10 21:54 编辑

【直播】我的基因组60:CNV初步探索

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

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


搜索一些CNV的简单资料放在这里吧
f5db8836684d2985beedba66a6d185bc.jpg
159e93f5faa89a52ac2af41d06162f95.jpg
参考文献:

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


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

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

[url=]http://bioconductor.org/packages/release/bioc/html/DNAcopy.html[/url]
说明书非常通俗易懂,就是接收每个探针对应区域的染色体号,探针坐标,以及该探针检测到的信号值。
03acef96cb55b516570afa7e8c311863.png

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


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

  • [AppleScript] 纯文本查看 复制代码
    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这个矩阵如下:
41dcd1e12b13a5928ada4aa544247576.png

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


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

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


  • [AppleScript] 纯文本查看 复制代码
    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文本文件结果给大家看到。

9ee8b19ce5f23aa1d062f2a8fc74b596.jpg

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

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


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

测试的代码如下:
[url=]http://bioconductor.org/packages/release/bioc/vignettes/DNAcopy/inst/doc/DNAcopy.R[/url]


文:Jimmy

图文编辑:吃瓜群众





上一篇:华大基因的BGI-SEQ500 公开GIAB标准品数据 开放下载!
下一篇:【直播】我的基因组61:scalpel软件找indel
基因组,转绿组,肿瘤信息,生物统计,Python, Linux.
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-12-16 11:36 , Processed in 0.027979 second(s), 30 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.