本帖最后由 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个基因，每个基因有四个情况的表达读数的文件。
文图显示各个不同软件所鉴定为差异表达的基因间的关系
1,stepwise挑选出有重要贡献基因 为了构建判别一个未知个体为两种不同状态中的哪一种，可以根据已知的变量进行推测，而原始的数据包含一万多个的变量，一开始我尝试了将所有变量均带入进行计算，但是系统很快就崩溃了，内存和存储空间，以及SAS无法自动利用到多线程处理，使得在PC上无法运算。于是，我将数据带入到不同的软件当中，并且分别得到判定为差异表达的基因，对其求得交集，即有53个基因被认为是差异表达。将其带入stepwise过程，结果只能留下三个基因，而再随机添加其它的100个基因组成153个基因，得到的结果如下图，选择的基因并不是其他软件共同判定为差异的基因。因此这个方法我无法进行下去，宣布失败（希望大家能给些建议为什么会这样，有什么方法能改进，感激）
具体实现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中提供了各个相关系数的显著性检验。
从主成分分析可以得出前三个主成分能够解释全部的变量。
而其相应各个原始的变量的系数如下图所示
而利用corr过程计算原始变量与主成分间的相关系数，结果如下图所示，其中的前53列为其他软件标识为差异表达的基因，而后100个为随机挑选的基因。从图中可以看到，前53个基因，每一个至少与一个主成分的相关系数的p值小于0.05，因此有显著的贡献，而后面的随机基因则可能对三个主成分并没有显著的贡献，也即其对于组的区分无显著贡献，可能为非差异表达基因。
实现的SAS代码： procprincomp data=merge2out=final1 outstat=final2standard; run; proccorr data=home.final1; varcol1-col153; withprin1-prin3; run;
3,判别函数的构建 根据主成分分析以及其他的软件的鉴定，可以以一定的概率说明共同的53个基因能够作为分类的标准，因此以此作为分类器的分类标准，利用discrim过程得到分类函数。由图可以看出对于原数据都能够正确分类，虽然样本非常小，并且没有测试数据集。结果如下：
实现的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]电子工业出版社，2010，297-302. | 
[复制链接]
发表于 2012-7-20 08:51
楼主
发表于 2012-9-25 13:21
