【直播】我的基因组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科技评论

开发 | 一文详解英伟达刚发布的 Tesla V100 究竟牛在哪?

AI科技评论按:很多读者在思考,“我和AI科技评论的距离在哪里?”答案就是:一封求职信。 5 月 11 日,在加州圣何塞举办的的 2017 年度 GPU 技术大...

333130
来自专栏数据小魔方

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

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

410130
来自专栏生信宝典

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

前面两篇文章(高通量数据分析必备|基因组浏览器使用介绍 - 1和高通量数据分析必备|基因组浏览器使用介绍 - 2)介绍了EPGG的基本使用、各部分特征、Trac...

13750
来自专栏思影科技

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

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

35480
来自专栏生信技能树

lncRNA数据分析传送门

step1: 计算资源的准备 如果有差不多配置的服务器,就可以从SRA/FASTQ格式数据开始走全套流程。不懂配置,请看前面转录组和表观组的传送门。 如果只有个...

65990
来自专栏量子位

老司机养成:教神经网络变身《马里奥赛车》高手 | 论文+代码

问耕 编译整理 量子位 出品 | 公众号 QbitAI ? 神经网络持续在游戏界立功,这次拿下的是经典游戏:《马里奥赛车64》,而且只需要很小的计算力就能完成。...

36360
来自专栏生信技能树

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

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

88890
来自专栏生信宝典

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

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

52650
来自专栏机器人网

线性执行元件的工作方式及分类

线性执行元件是一种以直线为基础进行能量转换的一种元件。线性执行元件可以根据应用者的要求而改变控制对象的状态,这种独特性能吸引着越来越多的人发现和应用它。线性执行...

30150
来自专栏数据小魔方

数据地图入门篇——素材获取!

从今天开始要跟大家分享新的专题——数据地图! 这一篇先讲一些准备性的操作,教大家怎么获取矢量地图素材,以及素材的编辑、加工和整理! 也许你曾见到过一些高大上的p...

67560

扫码关注云+社区

领取腾讯云代金券