寻找差异表达的基因

合集下载
  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

➢ 基因表达谱数据

基因表达谱可以用一个矩阵来表示,每一行代表一个基因,每一列代表一个样本(如图1)。所有基因的表达谱数据在“gene_exp.txt ”文件中存储,第一列为基因的entrez geneid ,第2~61列是疾病样本的表达,第62~76列是正常样本的表达。

图1 基因表达谱的矩阵表示

➢ 寻找差异表达的基因:

原理介绍:

差异表达分析是目前比较常用的识别疾病相关miRNA 以及基因的方法,目前也有很多差异表达分析的方法,但比较简单也比较常用的是Fold change 方法。它的优点是计算简单直观,缺点是没有考虑到差异表达的统计显著性;通常以2倍差异为阈值,判断基因是否差异表达。Fold change 的计算公式如下:

normal

Disease

x x c Fold =

_

即用疾病样本的表达均值除以正常样本的表达均值。

差异表达分析的目的:识别两个条件下表达差异显著的基因,即一个基因在两个条件中的表达水平,在排除各种偏差后,其差异具有统计学意义。我们利用一种比较常见的T 检验(T-test )方法来寻找差异表达的miRNA 。T 检验的主要原理为:对每一个miRNA 计算一个T 统计量来衡量疾病与正常情况下miRNA 表达的差异,然后根据t 分布计算显著性p 值来衡量这种差异的显著性,T 统计量计算公式如下:

n

s n s x x t normal Disease normal

Disease miRNA //22+-=

对于得到的显著性p 值,我们需要进行多重检验校正(FDR ),比较常用的是BH 方法(Benjamini and Hochberg, 1995)。

1+

=N

v

t分布

程序实现:

●基因表达谱数据--- gene_exp.txt

●Matlab软件实现mRNA差异表达分析:

MATLAB软件安装好之后,双击系统桌面的MATLAB图标,或在开始菜单的程序选项中选择MATLAB快捷方式,即开始启动MATLAB。初次启动MATLAB后,将进入MATLAB默认设置的桌面平台。桌面平台包括命令窗口、历史窗口、当前目录窗口和工作间管理窗口等窗口(如图2)。

图2 matlab窗口简介

工作空间主要包含了目前用户定义的一些变量,用户可以在命令窗口执行一些特定的命令操作来完成特定的功能。我们首先将工作目录选择到我们数据存放的硬盘目录下,然后导入要分析的基因表达谱数据,进行差异表达分析。

在命令窗口输入main_MTDN_end.m程序中的1-21行命令(注意要将程序中的目录改变到自己数据的存储目录下),即可得到差异表达的基因。这段程序主要包含两个函数:mattest和mafdr。

mattest函数是进行t检验的,输入的数据为疾病和正常的表达谱数据,返回每个miRNA的T统计量和对应的p值。这个参数还可以利用‘Permute’参数进行随机扰动,'Showhist'参数用来显示T统计量和p值的分布。

mafdr函数是用来计算FDR的函数,可以利用参数来选择计算FDR的方法,这里我们利用“BHFDR”参数来选择BH方法对p值进行校正,利用'showplot'参数来显示FDR的图示结果。

结果可以在工作空间窗口中通过双击变量进行查看。

结果展示:

T-统计量和p值的分布图以及FDR:

图3 T-score,P-values以及FDR的分布

●差异表达mRNA:我们卡的阈值为FDR<0.1;2倍fold change

(Fold_c>2 or <1/2 ),我们识别了11个下调的mRNA和6个上调的mRNA。

差异表达基因的层次聚类分析

➢mRNA表达谱数据:差异表达17个mRNA的表达数据

➢程序实现:

我们接下来利用差异表达mRNA的表达谱进行聚类分析,在命令窗口输入main_MTDN_end.m程序中的23-30行命令,结果会输出利用差异表达mRNA聚类分析的结果。这部分主要是利用一个现有的函数clustergram进行聚类分析,函数的输入数据是差异表达mRNA的表达谱。之后可以利用set 函数对行的符号和列的符号进行设定。

➢聚类分析结果展示:

➢聚类做heatmap,我比较喜欢用pheamap,简单又好看,但是很多做heatmap

的函数都不带输出聚类后基因名字的功能。heatmap旁标注基因是很有用的信息,

论文中经常会用到,所以我们可以更改pheatmap的源代码,让它输出基因列表,

其实如果能够给出基因list,在heatmap旁边标注出list中的基因就好了,但有了基因列表也可以做这个事情。

➢从cran上下载pheatmap的源代码,打开pheatmap的R文件夹中pheatmap.R文件,在一大串#上面添加write_matrix = function( mat,

out_file ){

➢write.table(as.data.frame(mat),sep="\t",quote=FALSE, file=out_file)

➢}

➢在一大串#下面的pheatmap中添加out_file = NA,此乃默认参数设定。➢

➢在hclust之后,就是当cluster_mat函数处理了mat矩阵后,添加

➢if( !is.na(out_file) ){

➢write_matrix( mat, out_file )

➢}

➢**************************************************昏割线

********************************************************

➢打开Rstudio,tools--install packages--选择那个压缩包ok啦

➢用法:

➢>setwd("F:/project/PTEN/01.RPKM/correlation")

相关文档
最新文档