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

相关文章

来自专栏点滴积累

geotrellis使用(二十七)栅格数据色彩渲染

目录 前言 复杂原因及思路分析 实现过程 总结 一、前言        今天我们来研究一下这个看似简单的问题,在地理信息系统中颜色渲染应当是最基本的操作和功能,...

3645
来自专栏葡萄城控件技术团队

位图和SVG用法比较

位图,亦称为点阵图像或绘制图像,是由称作像素(图片元素)的单个点组成的。这些点可以进行不同的排列和染色以构成图样。当放大位图时,可以看见赖以构成整个图像的无数单...

1866
来自专栏互联网数据官iCDO

A/B测试时间计算器&显著性计算器(Excel版)

译者:吕东昊 译者前言: 进行AB测试,不是直接做了比较谁好谁差那么简单。如果希望对比分析有意义,必须符合统计学意义上的真正的具有显著性的差异。不过,这个差异让...

3575
来自专栏大数据挖掘DT机器学习

PageRank算法(1):PageRank算法原理入门

一、PageRank简介 大名鼎鼎的PageRank算法是Google排名运算法则(排名公式)的一个非常重要的组成部分,其用于衡量一个网站好坏的标准。在揉合了...

2926
来自专栏机器学习养成记

《无问西东》豆瓣短评分析

《无问西东》讲述了四代人清华人,在矛盾与期待中不断找寻自我、砥砺前行的故事。上映后得到了广泛的关注,也引发了强烈的讨论。本文首先通过爬虫技术,从豆瓣上获得相关短...

3349
来自专栏AI研习社

如何用TensorFlow快速搭建神经网络?来看看DeepMind新开源工具Sonnet!

去年 DeepMind 作出决定,将全部研究搬到 TensorFlow 框架上进行。 近一年时间过去,回头来看,DeepMind 认为这项选择十分正确——许多模...

3379
来自专栏MixLab科技+设计实验室

从线条艺术到DIY实现一个网状体Net的js库

今天无意中看到一个可视化作品: WHAT MADE ME? INTERACTIVE PUBLIC INSTALLATION Most Original Exh...

2596
来自专栏深度学习之tensorflow实战篇

python 实现数据降维推荐系统(附Python源码)

主成分分析原理:请点击PCA查看 #!usr/bin/env python #_*_ coding:utf-8 _*_ import pandas as pd ...

3514
来自专栏程序生活

斯坦福tensorflow教程(四) 贪婪执行Eager Execution

1425
来自专栏数据小魔方

竟然是一个升级版的数据透视表,Tableau真的没有那么神秘~

前一篇跟大家详述了关于Excel在作图理念上的诸多细节,今天让我们把聚焦于多分类维度的数据呈现问题。 当然今天这篇我会把视角从Excel切换到Tableau,没...

5757

扫描关注云+社区