【直播】我的基因组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
图文编辑:吃瓜群众
- HDU 1010 Tempter of the Bone【DFS经典题+奇偶剪枝详解】
- Ethereum - 以太坊项目
- COGS 144. [USACO Dec07] 魅力手镯【01背包复习】
- SQL Server 使用全文索引进行页面搜索
- HDU 1003 Max Sum【动态规划求最大子序列和详解 】
- HDU 1005 Number Sequence【多解,暴力打表,鸽巢原理】
- HDU 1019 Least Common Multiple【gcd+lcm+水+多个数的lcm】
- HDU 1017 A Mathematical Curiosity【水,坑】
- 比特币项目
- HDU 1014 Uniform Generator【GCD,水】
- 【AlphaGo Zero 核心技术-深度强化学习教程代码实战05】SARSA(λ)算法实现
- 区块链应用场景:物联网和物流供应链
- HDU 1012 u Calculate e【暴力打表,水】
- Gym 100952C&&2015 HIAST Collegiate Programming Contest C. Palindrome Again !!【字符串,模拟】
- JavaScript 教程
- JavaScript 编辑工具
- JavaScript 与HTML
- JavaScript 与Java
- JavaScript 数据结构
- JavaScript 基本数据类型
- JavaScript 特殊数据类型
- JavaScript 运算符
- JavaScript typeof 运算符
- JavaScript 表达式
- JavaScript 类型转换
- JavaScript 基本语法
- JavaScript 注释
- Javascript 基本处理流程
- Javascript 选择结构
- Javascript if 语句
- Javascript if 语句的嵌套
- Javascript switch 语句
- Javascript 循环结构
- Javascript 循环结构实例
- Javascript 跳转语句
- Javascript 控制语句总结
- Javascript 函数介绍
- Javascript 函数的定义
- Javascript 函数调用
- Javascript 几种特殊的函数
- JavaScript 内置函数简介
- Javascript eval() 函数
- Javascript isFinite() 函数
- Javascript isNaN() 函数
- parseInt() 与 parseFloat()
- escape() 与 unescape()
- Javascript 字符串介绍
- Javascript length属性
- javascript 字符串函数
- Javascript 日期对象简介
- Javascript 日期对象用途
- Date 对象属性和方法
- Javascript 数组是什么
- Javascript 创建数组
- Javascript 数组赋值与取值
- Javascript 数组属性和方法
- docker安装伏羲扫描器fuxi-scanner
- 基于深度学习的文本分类应用!
- 表驱动法
- mysql将表结构导出excel
- 为什么会是Docker?
- 浅析http报文
- MySQL explain 中的 rows 究竟是如何计算的?
- SwiftUI: 使用 Touch ID 和 Face I
- Linux 系统中查找正在运行的进程的完整命令、当前工作目录等信息的方法
- Go by Example 中文:通道方向
- mycat数据库集群系列之mysql主从同步设置
- Tun/Tap接口使用指导
- Swift中? 、! 和 ??
- 故障分析 | 记一次 MySQL 主从双写导致的数据丢失问题
- 集成 SpringBoot 2.3.2 + Shiro 1.5.3 + jwt (无状态)