【直播】我的基因组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 条评论
登录 后参与评论

相关文章

来自专栏CSDN技术头条

Apache Kylin 深入Cube和查询优化

近几年,Apache Kylin作为一个高速的开源分布式大数据查询引擎正在迅速崛起。它充分发挥Hadoop、Spark、HBase等技术的优势,通过对超大规模数...

3038
来自专栏腾讯大数据的专栏

TDW千台Spark千亿节点对相似度计算

相似度计算在信息检索、数据挖掘等领域有着广泛的应用,是目前推荐引擎中的重要组成部分。随着互联网用户数目和内容的爆炸性增长,对大规模数据进行相似度计算的需求变得...

30210
来自专栏PPV课数据科学社区

【学习】七天搞定SAS(三):基本模块调用(格式、计数、概要统计、排序等)(下)

SAS里面总结数据:MEANS SAS当然还有类似于excel的数据透视表和R的data.table的模块,就是MEANS。可以输出的summary stat...

3056
来自专栏生信技能树

单细胞转录组3大R包之Seurat

牛津大学的Rahul Satija等开发的Seurat,最早公布在Nature biotechnology, 2015,文章是; Spatial reconst...

2.9K19
来自专栏数据小魔方

带实际执行进度的甘特图

今天要跟大家分享的图标是带实际执行进度的甘特图! ▽▼▽ 由于本图所用到的技巧和思路特别复杂,过程相对繁琐,所以本案例的介绍会省略掉很多细节性的步骤,否则图文会...

3235
来自专栏蓝鸟资源分享网

关于服务器性能的一些思考

平常的工作中,在衡量服务器的性能时,经常会涉及到几个指标,load、cpu、mem、qps、rt,其中load、cpu、mem来衡量机器性能,qps、rt来衡量...

1134
来自专栏大数据和云计算技术

Automatic Management of Data and Computation in Datacenters

image.png 最近在研究数据中心的数据管理和性能优化,看了一篇2010的论文Nectar:Automatic Management of Data and...

30810
来自专栏逍遥剑客的游戏开发

挑战蓝龙Kalecgos

922
来自专栏CSDN技术头条

Apache Spark 1.5新特性介绍

Apache Spark社区刚刚发布了1.5版本,大家一定想知道这个版本的主要变化,这篇文章告诉你答案。 DataFrame执行后端优化(Tungsten第一阶...

1869
来自专栏祝威廉

谷歌BigQuery ML VS StreamingPro MLSQL

今天看到了一篇 AI前线的文章谷歌BigQuery ML正式上岗,只会用SQL也能玩转机器学习!。正好自己也在力推 StreamingPro的MLSQL。 今...

863

扫码关注云+社区