查看: 19618|回复: 21

RNA-seq差异表达基因分析方法探索[sas]

  [复制链接]
新浪微博达人勋 pypy 未实名认证
论坛徽章:
2
R研习者中级
日期:2012-04-29 00:26:34R研习者中级
日期:2013-06-13 18:59:20
发表于 2012-7-20 08:51 | 显示全部楼层 |阅读模式
本帖最后由 pypy 于 2012-7-20 09:04 编辑


RNA-seq差异表达基因分析方法探索
By: hitchpy
前言:
    RNA-seq技术是通过对细胞中的RNA进行序列测定并且通过比对到已知的注释了功能的基因组上,可以测得不同的基因的表达数目大小。因此可以通过这项技术研究生物体在不同的环境下,或者不同生理状态下的基因表达量之间的差异,从而能够了解身体的反应机制以及构建细胞内的调节网络。
但是并不是所有的基因都会在不同状态下有差异的表达,而且同一个基因在同一个状态下,由于实验操作或者样本间的差异,测得的表达量读数也不一定完全相等。因此,如何能够鉴定出真正为差异表达的基因,能够更加较精确的构建调控网络,或者在医学诊断上,通过对少数几个基因的表达进行测定即能够达到诊断效果,将是一个重要的目标。
而在公式当中的表示,每一个基因相当于一个变量,测定的表达量即为其一次取值,而并不是差异表达的基因为辅助的变量,而我们的目标就是降维作用,找出主要的贡献基因。
现在存在的R程序中即有分析基因差异表达的包,如bioconductor中的edgeR,bayseq等,他们都是通过拟定基因在随机的抽样测序中是按照负二项分布的假设而通过数据进行分布参数的估计而得到具体分布,并且计算各个基因间是否有显著的差异,依据FDR(false discovery rate)等为标准。具体各个包间差异可能包括模型的选择差别,参数估计方式等的不同。
这里,我试图利用SAS中判别分析中的stepwise过程选择变量值的方式,主成分分析的方式选择贡献较高的基因,并与通过几种软件判定为差异表达的基因做比较。并且能够利用判定为重要的基因,利用DISCRIM过程,构建判别函数,能够辅助以后医学方面疾病的鉴定等的功能.


数据:
    数据为中山大学生命科学学院的某实验室内对果蝇的一个基因片段进行敲除,与对照组一同进行RNA-seq测序,并通过htseq-count软件得到的每个基因片段的表达量读数值。其中敲除组有两个样本,对照组也有两个样本。即包含14283个基因,每个基因有四个情况的表达读数的文件。

jj.png

文图显示各个不同软件所鉴定为差异表达的基因间的关系


1,stepwise挑选出有重要贡献基因
    为了构建判别一个未知个体为两种不同状态中的哪一种,可以根据已知的变量进行推测,而原始的数据包含一万多个的变量,一开始我尝试了将所有变量均带入进行计算,但是系统很快就崩溃了,内存和存储空间,以及SAS无法自动利用到多线程处理,使得在PC上无法运算。于是,我将数据带入到不同的软件当中,并且分别得到判定为差异表达的基因,对其求得交集,即有53个基因被认为是差异表达。将其带入stepwise过程,结果只能留下三个基因,而再随机添加其它的100个基因组成153个基因,得到的结果如下图,选择的基因并不是其他软件共同判定为差异的基因。因此这个方法我无法进行下去,宣布失败(希望大家能给些建议为什么会这样,有什么方法能改进,感激)
n1.png

具体实现SAS语句:
DATAhome.new;
infile"C:\Users\PY\Desktop\result.txt";
inputname $ a1 a2 b1 b2;
run;
datahome.ran;
infile"C:\Users\PY\Desktop\random.txt";
inputname $ a1 a2 b1 b2;
run;

datahome.merge;
sethome.new home.ran;
run;
proctranspose data=home.newout=home.new2 name=classlabel=class;

run;

分别读取被其他软件鉴定为差异表达的基因,以及选取的100个随机基因,进行合并,然后进行转置。
procstepdisc data=home.test;
classclass;
run;

利用stepdisc过程进行变量选择。



2,利用主成分分析确定重要的变量
    主成分分析也是一种常用的降维的方式,通过将各个变量进行正交分解,取定其中权重较大的部分而实现降维。而通过计算原始变量与各个主成分间的相关系数,能够确定各个原始变量是否对新变量有显著的贡献,也即是有显著差异的基因。同样,这里对153个基因进行分析,得到主成分分析的结果以及计算原始变量与主成分间的相关系数,同时SAS中提供了各个相关系数的显著性检验。
n8.png

从主成分分析可以得出前三个主成分能够解释全部的变量。

n4.png
而其相应各个原始的变量的系数如下图所示
n5.png

而利用corr过程计算原始变量与主成分间的相关系数,结果如下图所示,其中的前53列为其他软件标识为差异表达的基因,而后100个为随机挑选的基因。从图中可以看到,前53个基因,每一个至少与一个主成分的相关系数的p值小于0.05,因此有显著的贡献,而后面的随机基因则可能对三个主成分并没有显著的贡献,也即其对于组的区分无显著贡献,可能为非差异表达基因。
n2.png

n3.png



实现的SAS代码:
procprincomp data=merge2out=final1 outstat=final2standard;
run;
proccorr data=home.final1;
varcol1-col153;
withprin1-prin3;
run;

3,判别函数的构建
    根据主成分分析以及其他的软件的鉴定,可以以一定的概率说明共同的53个基因能够作为分类的标准,因此以此作为分类器的分类标准,利用discrim过程得到分类函数。由图可以看出对于原数据都能够正确分类,虽然样本非常小,并且没有测试数据集。结果如下:
n6.png

n7.png


实现的SAS代码如下:
procdiscrim distancemanova data=home.new2;
classclass;
varcol1-col53;
run;


总结:
本文章的分析是根据之前利用bioconductor上的bayseq,deseq,edgeR以及cuffdiff四个软件对数据的分析,得到的53个显著性基因作为基础而进行的进一步的探索分析。由于本身的知识水平的限制,对于许多现有的更为较精确的模型无法更好的采用。因此在这里,希望通过测试比较通用基础的主成分分析方法和变量筛选方法来得到一定的分析结果。从结果来看,主成分分析也能较好的辨别出较为显著差异表达的基因。虽然这个地方并没有对所有的14000多个的基因进行分析。基因的分析是现在生物信息学中的最主要的目标,而差异基因的辨别也是其中重要的一个分支,已经有很多工作在这个方面展开。并且可能利用如SVM等的工具,能够从新的角度利用所有基因表达数据的信息,避免所谓的维数灾难,而采用的降维的方式。
除此之外,之所以现在差异基因的判别出现难度,也是因为现有的技术无法支持高昂的多样本实验,所以以后技术的进步也将能够显著的提高统计的准确性。
P.S建议阅读leo breiman的Statistical Modeling: The Two Cultures。网上可以搜索到PDF,很有启发意义,与大家共同分享学习!

参考文献:
[1]Robinson, M. D. and G. K. Smyth Moderated statistical tests for assessingdifferences in tag abundance. Bioinformatics  2007 23(21): 2881-2887.
[2]Anders, S. and W. Huber Differential expression analysis for sequence countdata. Genome Biology  2010 11(10).
[3]Hardcastle, T. J. and K. A. Kelly  baySeq: Empirical Bayesian methods foridentifying differential expression in sequence count data. BmcBioinformatics 2010 11.
[4]董大均.SAS统计分析应用.[M]电子工业出版社,2010297-302.
py_data.zip (129.48 KB)
回复

使用道具 举报

新浪微博达人勋 pypy 未实名认证
论坛徽章:
2
R研习者中级
日期:2012-04-29 00:26:34R研习者中级
日期:2013-06-13 18:59:20
 楼主| 发表于 2012-7-23 18:46 | 显示全部楼层
再后续贴上一个自己制作的R包,其实这些都是在我的毕业论文的基础上做的工作……上面的是对数据再进一步的分析,而这个包主要涉及之前数据在三个软件跑的时候的内部操作。
   把这个过程做成一个包,主要是因为之前要求做R包的时候,我没有完成……所以这段时间来都一直受到良心的煎熬,能让自己安心的就是做一个包出来。
   这个包非常简单,基本上就是把其他程序的语句贴在里边,进行了细微的整合。而且参数都很简单,如果刚好遇见需要的人,会很方便,不过……99.9%的人都不会用到这个包。
   之前的构思是做一个对各个数据框进行操作的包,比如根据某列的指标选出符合要求的行,合并等操作,毕竟一条一条R指令还是比较烦,过段时间可以完成这个顺手小工具,这样excel就大部分被淘汰了!没有做那个包是因为要求要是关于某行业的统计问题……我是做了统计问题,只是涉及统计的都是包装好了的,算是钻个空子……
    总之,亲手做过一个包才有感觉,不尝试永远心虚。(PS……还有很大一部分原因是怕缺少一个大作业期末总成绩会悲剧!所以,希望老师高抬贵手,那段时间做毕业论文,发烧,准备期末……还有我懒……)
    好了,下边是附件R包,还有check和例子运行结果截图。 n1.png n2.png

cdgene_1.0.zip (12.22 KB)
回复 支持 反对

使用道具 举报

论坛徽章:
4
SAS研习者初级
日期:2012-08-20 22:27:43R研习者初级
日期:2012-09-17 19:50:38Hadoop研习者初级
日期:2012-09-28 23:48:40Hadoop研习者初级
日期:2013-10-21 22:39:48
发表于 2012-7-29 02:00 | 显示全部楼层
读了楼主的帖子很有收获,最近也在考虑用sas做类似的分析,我想问一个关于测序的技术问题,请问两组的测序深度是怎样了,覆盖率怎么样,是否可以把生成的原始文件fastq,或者 eland格式的文件切成同样大小,在进行后续分析? 有人说这样做可以因为read都是随机的,所以可以随便切文件。你觉得哪?
回复 支持 反对

使用道具 举报

新浪微博达人勋 pypy 未实名认证
论坛徽章:
2
R研习者中级
日期:2012-04-29 00:26:34R研习者中级
日期:2013-06-13 18:59:20
 楼主| 发表于 2012-7-29 13:46 | 显示全部楼层
sas-67 发表于 2012-7-29 02:00
读了楼主的帖子很有收获,最近也在考虑用sas做类似的分析,我想问一个关于测序的技术问题,请问两组的测序深 ...

well……当时数据是老师直接给我的,因此覆盖率我也不太清楚,不过是不是能够通过所有基因的总读数大概估计出为600万。应该是的吧,读数机器操作是random的,因此可以随意切吧,可是能够利用全部的读取结果不会是更好的么?毕竟深度越大会越较精确?
回复 支持 反对

使用道具 举报

论坛徽章:
4
SAS研习者初级
日期:2012-08-20 22:27:43R研习者初级
日期:2012-09-17 19:50:38Hadoop研习者初级
日期:2012-09-28 23:48:40Hadoop研习者初级
日期:2013-10-21 22:39:48
发表于 2012-8-1 12:33 | 显示全部楼层
我也是这么认为的,不过周围有的同事就是这么做的,不管来的数据多大都一律切成一边大,说深度不会对DE分析造成太大影响,理由是搞表达的基因只会因为深度的增加更加高,低表达的变化不大,所以高的只会更高,低的还是很低
回复 支持 反对

使用道具 举报

论坛徽章:
4
SAS研习者初级
日期:2012-08-20 22:27:43R研习者初级
日期:2012-09-17 19:50:38Hadoop研习者初级
日期:2012-09-28 23:48:40Hadoop研习者初级
日期:2013-10-21 22:39:48
发表于 2012-8-1 12:40 | 显示全部楼层
还有一个菜鸟问题,请问你第一个那个交集的图是怎么弄得?很好看,我也想做个类似的
谢谢了
回复 支持 反对

使用道具 举报

新浪微博达人勋 pypy 未实名认证
论坛徽章:
2
R研习者中级
日期:2012-04-29 00:26:34R研习者中级
日期:2013-06-13 18:59:20
 楼主| 发表于 2012-8-1 15:36 | 显示全部楼层
sas-67 发表于 2012-8-1 12:40
还有一个菜鸟问题,请问你第一个那个交集的图是怎么弄得?很好看,我也想做个类似的
谢谢了

http://bioinfogp.cnb.csic.es/tools/venny/index.html
这个网站直接所见即所得,不过R里边也有包可以做,我没仔细去找
回复 支持 反对

使用道具 举报

论坛徽章:
4
SAS研习者初级
日期:2012-08-20 22:27:43R研习者初级
日期:2012-09-17 19:50:38Hadoop研习者初级
日期:2012-09-28 23:48:40Hadoop研习者初级
日期:2013-10-21 22:39:48
发表于 2012-8-3 03:36 | 显示全部楼层
非常感谢,对我很有帮助
回复 支持 反对

使用道具 举报

论坛徽章:
0
发表于 2012-9-25 13:21 | 显示全部楼层
RNA-seq是热门题目, 谢谢楼主详细的案例分享;

下面是生物通RNA-seq找癌症融合基因的简介
http://124.172.245.86/newsf/read.asp?page=2012424134839300
回复 支持 反对

使用道具 举报

新浪微博达人勋 pypy 未实名认证
论坛徽章:
2
R研习者中级
日期:2012-04-29 00:26:34R研习者中级
日期:2013-06-13 18:59:20
 楼主| 发表于 2012-10-13 20:00 | 显示全部楼层
genesquared 发表于 2012-9-25 13:21
RNA-seq是热门题目, 谢谢楼主详细的案例分享;

下面是生物通RNA-seq找癌症融合基因的简介

惭愧,一直都没有关注了,现在还是在打数学基础,那个是本科水出来的论文。不过以后还是希望能够继续深造生物信息方向的。这几年来闲时看的论文也挺多的,可惜自己动手能力还是要跟进
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

 

GMT+8, 2019-8-21 21:41 , Processed in 0.154233 second(s), 52 queries .