【直播】我的基因组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代码如下:

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文本文件结果给大家看到。

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

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

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

测试的代码如下:

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

文:Jimmy

图文编辑:吃瓜群众

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2017-03-09

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏AI研习社

P图现在都这么狠了吗?英伟达新算法,轻松帮你P掉前男友

阳光好,阳光不好;心情好,心情糟;今天的妆一级棒,今天的黑眼圈快要掉下来了!时时刻刻都想用照片记录,我的天,前男友和发黄的墙纸怎么还不消失?

1323
来自专栏生信宝典

高通量数据分析必备|基因组浏览器使用介绍 - 1

基因组浏览器是高通量测序分析的一个重要的可视化工具。相比于最终提供的表格,基因组浏览器可以提供更多的信息,如直观展示突变位点、查看有无新转录本或新的可变剪接形式...

1222
来自专栏生信技能树

【直播】我的基因组47:测序深度和GC含量的关系

在前面我们用 ChIP-seq 的分析方法可视化了一下我的 WGS数据,结果我们的测序深度分布居然是跟基因组的genomic feature相关。 比如在TSS...

6079
来自专栏数据小魔方

R语言可视化——ggplot图表系统中的辅助线

在之前的推送中,曾经有过一篇介绍excel图表辅助线的制作方法,其中用到的技巧五花八门、令人眼花缭乱。 而ggplot图表系统中的辅助线添加起来却异常简单,非常...

37713
来自专栏深度学习入门与实践

【原】Spark之机器学习(Python版)(二)——分类

  写这个系列是因为最近公司在搞技术分享,学习Spark,我的任务是讲PySpark的应用,因为我主要用Python,结合Spark,就讲PySpark了。然而...

2316
来自专栏思影科技

异质脑:自闭症谱系障碍患者自发连接模式畸变

来自以色列魏茨曼科学研究学院的Avital Hahamy等人在Nature neuroscience上发表文章,发现自闭症谱系障碍(Autism spectru...

3398
来自专栏数据小魔方

R语言可视化——图表美化与套用主题(下)

昨天的分享跟大家简单介绍了关于柱形图图表元素美化的思路,今天接着分享关于套用主题。 因为单独使用代码来调整单个图表元素,实在是太费劲了,更何况图表的细节元素有那...

2766
来自专栏北京马哥教育

干货 | 10 行 Python 代码创建可视化地图

当我开始建造Vincent时, 我的一个目的就是使得地图的建造尽可能合理化. 有一些很棒的python地图库-参见Basemap 和 Kartograph能让地...

4107
来自专栏生信宝典

生信宝典之傻瓜式(四)蛋白蛋白互作网络在线搜索

傻瓜系列重启了,今天要介绍的是一款在线查询蛋白-蛋白互作网络的工具 STRING (https://string-db.org/)。 STRING数据库收录了2...

3725
来自专栏AI研习社

官方博客:英伟达的新卡如何从硬件上支持了深度学习

编者按:5 月 11 日,在加州圣何塞举办的 2017 年度 GPU 技术大会上,英伟达发布了 Tesla V100,号称史上最强的 GPU 加速器。发布之后,...

3185

扫码关注云+社区