您所在位置： &
&nbsp&&nbsp&nbsp&&nbsp
基因表达数据实验指导解决方案.doc 40页
本文档一共被下载：
次 ,您可全文免费在线阅读后下载本文档。
下载提示
1.本站不保证该用户上传的文档完整性，不预览、不比对内容而直接下载产生的反悔问题本站不予受理。
2.该文档所得收入(下载+内容+预览三)归上传者、原创者。
3.登录后可充值，立即自动返金币，充值渠道很便利
需要金币：350 &&
你可能关注的文档：
··········
··········
基因表达数据分析实验指导
1.?实验基本情况
2.?实验方法：
2.1?表达谱数据的下载
2.2?将表达谱数据导入matlab软件
2.3?补缺失值
2.4?数据标准化
2.5?差异表达基因筛选
2.6?选择差异表达的基因
2.7对差异表达基因送入功能注释
附?-- Matlab的Microarray Data Analysis
1.?实验基本情况
实验目的：
掌握和了解常用的基因表达分析过程，包括数据下载、数据预处理、差异表达分析和基因功能注释。了解GEO、SMD、Matlab软件和WebGestalt数据库的使用。
实验方法：
详见下面的描述。
实验作业：
每位同学从GEO或SMD数据库上下载一套表达谱数据，进行数据预处理，差异表达基因分析或聚类分析等数据分析过程（依据具体问题操作，arraytool或matlab或其他软件均可），基因功能注释（WebGestalt、GO、KEGG等数据库）。
实验实例分析
=====================================================================
2.?实验方法：
2.1?表达谱数据的下载
2.1.1?从GEO数据库上下载表达谱数据
1)?网址及数据库概述
GEO主页：/geo/
GEO数据库中包含四种类型的条目，分别以GPLXXXX（检测平台），GSMXXXX（生物样本），GSEXXXX（基因表达系列），GDSXXXX（基因表达数据集）表示。其中GPLXXXX有SAGE、MPSS、单色芯片（Affymetrix）、双色芯片（spotcDNA/DNA）几种；GSEXXXX与GDSXXXX的区别在于：GSE是实验者一次一起提交的数据集，包含原始的数据文件，而GDS是GEO数据库的维护者根据样本和实验平台的特性进行整理的，与原有的GSE数据可能有样本量上的差异；一般GDS都有对应的GSE数据；GDS不包含单独的原始数据，如果想获得其原始数据，需要链接到他的GSE网页上下载；GDS样本间的可比性更强，如果有GDS就先分析GDS。
2）数据下载
GEO可提供两种数据的下载，一种是整理好的soft格式数据，是一个数据矩阵，包含多个基因在多个条件下的表达值，如GDS2220.soft；另一种是单独的数据文件，每张芯片一个数据表格，如GSE3519_family.xml文件夹里的文件，就是对应GDS2220这次实验的原始数据。另外还有一个GDS2220.annot数据是提供基因描述的。具体的下载方式如下：
在GEO主页上（图1），可以通过浏览（browse）或query中输入疾病名字，如风湿性关节炎（rheumatoid?arthritis）在Datasets后，点击go进行搜索，结果如图2。
图1. GEO的主页
图2.GEO的搜索结果
之后点击你感兴趣的GDS集合，如GDS2220，就进入每套数据的页面了（图3）。
图3.GDS2220数据的浏览界面
在图3中，点击下拉菜单中的DataSet?SOFT file，就能下载GDS2220.soft文件；点击Annotation SOFT file就可以下载GDS2220.annot文件；点击seriers?family?miniml?file就可以下载GSE3519_family.xml文件夹,但这个速度较慢，这里有个小窍门，大家可以在迅雷中新建一个下载任务，粘贴地址：?/pub/geo/DATA/MINiML/by_series/GSE139/GSE139_family.xml.tgz?，这里GSE139是可以替换的，比如要下载GDS2220配套的数据，就直接把两个GSE139都替换成GSE3519就可以直接下载了；点击series family soft file下载的文件与GDS2220.soft类似，只是样本是GSE3519的数据，可能和GDS2220的样本不同，这里是相同的。
也可以通过以下方式寻找特殊平台的数据。
3)?文件描述
（a）GDS22.soft
该文件从上到下分为三个部分：第一部分，数据集合基本描述，文字形式，以！或#开头；第二部分，表格的表头，如“ID_REF??????????????????IDENTIFIER??????GSM80309?????????GSM80310?????????GSM80311?????????GSM80312?????????GSM80313?????????GSM80314?????????GSM80315????????GSM80316?????????GSM80317”，以tab键分割，表示下面的数据部分每一列的含义；第三部分，数据，如GDS2220.soft中第一列为每一个基因的编号，第二列是基因
正在加载中，请稍后... 上传我的文档
 下载
 收藏
该文档贡献者很忙，什么也没留下。
 下载此文档
MATLAB7_X生物信息工具箱的应用_基因表达图谱分析_4_杨英杰
下载积分：100
内容提示：MATLAB7_X生物信息工具箱的应用_基因表达图谱分析_4_杨英杰
文档格式：PDF|
浏览次数：24|
上传日期： 11:21:28|
文档星级：
全文阅读已结束，如果下载本文需要使用
 100 积分
下载此文档
该用户还上传了这些文档
MATLAB7_X生物信息工具箱的应用_基因表达图谱分析_4_杨英杰
关注微信公众号生物秀人才网[]
打造生物医药领域人才求职与企业招聘的专业平台
欢迎加入最专业的生命科学交流社区，结交业内好友，体验更多功能。
才可以下载或查看，没有帐号？
比如我用的一个GSE50081这个数据，我用R里面的getGEO下载下来的数据，exprs(gse)之后得到的基因表达矩阵的行名是探针的ID，如果要找对应的基因symbol的话我就去featureData(gse)@data$`Gene Symbol`，结果发现只有500多个，其他的全部都是NA。但是如果我去soft文件里面就可以找到探针ID对应的基因symbol，这到底是怎么回事呢？
美开发细菌制造生物燃料技术电泳单引物跑出来为什么是这个样子电泳单引实时荧光定量PCR出现很多乱峰，求大神们帮晒动图——活体细胞4D成像新手段，细胞筛选《自然遗传学》超25万人基因组数据，揭示糖
12345678910
wx_非凡人_lDsv9大叔级摄影爱好者，喜欢分享murongxiaoxiao积极有责任心，热心公益事业绝对零度红米达人，爱拍照的北京女孩
逛了这许久，何不进去瞧瞧？基因表达谱的生物信息学_图文_百度文库
两大类热门资源免费畅读
续费一年阅读会员，立省24元！
基因表达谱的生物信息学
&&数据挖掘 生物信息 模式识别 DNA微阵列
阅读已结束，下载文档到电脑
想免费下载更多文档？
定制HR最喜欢的简历
下载文档到电脑，方便使用
还剩4页未读，继续阅读
定制HR最喜欢的简历
你可能喜欢library(affy)
library(tcltk)
filters&-matrix(c(&CEL file&, &.[Cc][Ee][Ll]&, &All&,&.*&), ncol = 2,byrow = T)
cel.files &-tk_choose.files(caption = &Select CELs&, multi =TRUE,filters= filters, index = 1)
data.raw&-ReadAffy(filenames = cel.files)
n.cel &-length(cel.files)
sampleNames(data.raw)
==============================================
#生成0h，1h，24h，7d四个值依次重复两次所组成的数列，数列命名为treatment
pData(data.raw)
#指针pData函数读取文件
==============================================
1、计算基因表达量
eset.rma &-rma(data.raw)
eset.mas5&-mas5(data.raw)
%%注意rma的eset结果是经过对数变换的，而mas5的eset结果是原始信号强度。虽然表达量是用对数变换的信号值表示的，但是有些计算过程要用到未经变换的原始值，应该把它们都计算出来：
emat.rma.nologs&-2^emat.rma.log2
emat.mas5.log2 &-log2(emat.mas5.nologs)
results.rma&-data.frame((emat.rma.log2[,c(1,3,5,7)] + emat.rma.log2[,c(2,4,6,8)])/2)
#计算平均值，并转换为数据框格式#计算表达量差异倍数
subset.logic &-results.rma$fc.7d&0
subset.data &- results.rma[subset.logic,]
#要注意的是逻辑向量的长度要和相应维度的数据长度一致，逻辑向量中为TRUE的就保留，FALSE的就丢弃。
2、选取表达基因
%选取“表达”基因的方法常见的有两种，一是使用genefilter软件包，另外一种是调用affy包的mas5calls()函数。使用 genefilter需要设定筛选阈值，不同的人可能有不同的标准。mas5calls方法使用探针水平数据（AffyBatch类型数据）进行处理，一般使用没经过预处理的芯片数据通用性强些，其他参数用默认就可以。
data.mas5calls &-mas5calls(data.raw)
eset.mas5calls&-exprs(data.mas5calls)
#继续用exprs计算“表达”量，得到的数据只有三个值P/M/A。对于这三个值的具体解释可以用?mas5calls查看帮助。P为present，A为absent，M为marginal（临界值）。
#把至少一个芯片中有表达的基因选出来！！！！
present.probes &-names(AP[AP])
results.present &- results.rma[present.probes,]
#present.probes是名称向量，用它进行数据子集提取。
3、获取差异表达基因
%生物学数据分析时的&差异&应该有两个意思，一是统计学上的差异，另外一个是生物学上的差异。一个基因在两个条件下的表达量分别有3个测量值：99,100,101 和 102,103,104。统计上两种条件下的基因表达数值是有差异的，后者比前者表达量要大。但生物学上有意义吗？未必。按平均值计算表达变化上升了3%，能产生什么样的生物学效应？这得看是什么基因了。所以差异表达基因的选取一般设置至少两个阈值：基因表达变化量和统计显著性量度（p值、q值等）。
3.1、t-测验
%经常使用的筛选阈值是表达量变化超过2倍，即|log2(fc)|&=log(2)。
==================================================
apply(abs(results.present[,5:7]),2, max)
%apply是一个很有用的函数，它对数据按某个维度批量应用一个函数进行计算。第一个参数为向量或矩阵（或者是能转成向量或矩阵的数据，如数据框），第三个参数表示要使用的函数，第二个参数为应用的维度。上面语句的意思是对数据 abs(results.present[,5:7]) 按列（第二维）使用统计函数max（计算最大值）。表达变化超过2倍的基因共有842个：
===========================================
## [1] 842, 选出这842个基因：
results.st &- results.present[abs(results.present$fc.7d)&=log2(2),]
sel.genes &- row.names(results.st)
, 1,function(x){t.test(x[1:2], x[7:8])$p.value})
results.st$p.value &- p.value
results.st &- results.st[, c(1,4,7,8)]
results.st &- results.st[p.value&0.05,]
head(results.st, 2)
## & & & & & & &S1 & &S7&fc.7d p.value
## 245042_at 8.153 7.021 -1.133 0.01004
## 245088_at 7.041 5.419 -1.622 0.03381
## t测验，并选出p&0.05的差异表达基因：
3.2 SAM（Significance Analysis of Microarrays）
library(BiocInstaller)
biocLite(&samr&)
library(samr)
samfit &- SAM(emat.rma.nologs[present.probes,c(1,2,7,8)],c(1,1,2,2),&resp.type=&Twoclass unpaired&, genenames=present.probes)
SAM函数返回值一个列表结构，可以自己用?SAM看看。差异表达基因的数据在siggenes.table中，也是一个列表结构：
str(samfit$siggenes.table)
%上调基因在siggenes.table的genes.up中，下调基因在genes.lo中。从上面的数据结构显示还可以看到差异表达基因的数量： ngenes.up和ngenes.lo。提取差异表达基因数据：
results.sam&-data.frame(rbind(samfit$siggenes.table$genes.up,samfit$siggenes.table$genes.lo),
&row.names=1, stringsAsFactors=FALSE)
for(i in 1:ncol(results.sam)) results.sam[,i] &- as.numeric(results.sam[,i])
head(results.sam, 2)
3.3、Wilcoxon's signed-rank test
%这个方法发表在 Liu, W.-m. etal, Analysis of high density expression microarrays with signed-rank callalgorithms. Bioinformatics, 93-1599。R软件包simpleaffy的detection.p.val函数有实现，可以通过parison函数调用：
library(simpleaffy)
#注意下面语句中的数据顺序
sa.fit &- parison(eset.rma, &treatment&,c(&7d&, &0h&))
%parison返回的数据为simpleaffy自定义的&PairComp&类型，提取数据要用它专门的函数：平均值用means函数获得，变化倍数（log2）用fc函数获得，t测验的p值用tt函数获得：
results.sa &- data.frame(means(sa.fit), fc(sa.fit), tt(sa.fit))
#选择有表达的基因
results.sa &- results.sa[present.probes,]
head(results.sa, 2)
## & & & & & & X7d & X0h fc.sa.fit.tt.sa.fit.
## 244901_at 4.047 4.203 & &-0.1562 & &0.43982
## 244902_at 3.938 4.295 & &-0.3570 & &0.05824
应用表达倍数筛选得到表达倍数超过2倍的基因数量有862个，应用p值筛选后得到562个差异表达基因：
results.sa &- results.sa[abs(results.sa$fold.change)&=log2(2), ];nrow(results.sa)
## [1] 862
results.sa &- results.sa[results.sa$p.val&0.05,]; nrow(results.sa)
## [1] 562
3.4、Moderated Tstatistic——经验贝叶斯
%这种方法在R软件包limma里面实现得最好。limma最初主要用于双色（双通道）芯片的处理，现在不仅支持单色芯片处理，新版还添加了对RNAseq数据的支持，很值得学习使用。安装方法同前面其他Bioconductor软件包的安装。载入limm软件包后可以用limmaUsersGuide()函数获取pdf格式的帮助文档。limma需要先产生一个design矩阵，用于描述RNA样品：
library(limma)
design &- model.matrix(~ 0 + treatment)
colnames(design) &- c(&C0h&, &T1h&, &T24h&,&T7d&)
## & C0h T1h T24h T7d
## 1 & 1 & 0 & &0 & 0
## 2 & 1 & 0 & &0 & 0
## 3 & 0 & 1 & &0 & 0
## 4 & 0 & 1 & &0 & 0
## 5 & 0 & 0 & &1 & 0
## 6 & 0 & 0 & &1 & 0
## 7 & 0 & 0 & &0 & 1
## 8 & 0 & 0 & &0 & 1
##可以看到：矩阵的每一行代表一张芯片，每一列代表一种RNA来源（或处理）。此外，你可能还需要另外一个矩阵，用来说明你要进行哪些样品间的对比分析：
contrast.matrix &- makeContrasts(T1h-C0h, T24h-C0h, T7d-C0h, levels=design)
contrast.matrix
## & & & Contrasts
## Levels T1h - C0h T24h - C0h T7d - C0h
## & C0h & & & & -1 & & & & -1& & & &-1
## & T1h & & & & &1 & & & &&0 & & & & 0
## & T24h & & & & 0 & & & &&1 & & & & 0
## & T7d & & & & &0 & & & &&0 & & & & 1
fit &- lmFit(eset.rma[present.probes,], design)
fit1 &- contrasts.fit(fit, contrast.matrix)
fit2 &- eBayes(fit1)
#首先进行线性模型拟合，根据对比模型进行差值计算，最后是贝叶斯检验
先统计一下数量。可以看到：对于T7d-C0h比较组（coef=3），表达变化超过2倍（lfc参数）的基因数为842个，而变化超过2倍、p&0.05的基因总数为740个：
nrow(topTable(fit2, coef=3, adjust.method=&fdr&, lfc=1,number=30000))
## [1] 842
nrow(topTable(fit2, coef=3, adjust.method=&fdr&, p.value=0.05, lfc=1,number=30000))
## [1] 740
##把toTable函数的返回结果存到其他变量就可以了，它是数据框类型数据，可以用write或write.csv函数保存到文件：
results.lim &- topTable(fit2, coef=3, adjust.method=&fdr&,p.value=0.05, lfc=1, number=30000)
%为什么以上几种方法仅用表达倍数（2倍）筛选得到的数字不大一样？limma和直接计算的结果都是842个，而simpleaffy和SAM为862/861个。这是对eset信号值取对数和求平均值的先后导致的，limma先取对数再求平均值，而simpleaffy和SAM是先求平均值再取对数。
3.5 其他方法：
如Rank products方法，在R软件包RankProd里实现，方法文献为：Breitling R, Armengaud P, Amtmann A, etal. Rank products: a simple, yet powerful, new method to detect differentiallyregulated genes in replicated microarray experiments[J]. FEBS letters, ): 83-92.
最后我们保存部分数据备下次使用：
#所有表达基因的名称
write(present.probes, &genes.expressed.txt&)
#处理7天的差异表达基因
write.csv(results.lim, &results.lim.7d.csv&)
#emat.rma.log2
write.csv(emat.rma.log2[present.probes,], &emat.rma.log2.csv&)
如果要全部结果：
results.lim.all &- topTable(fit2, coef=1:3, adjust.method=&fdr&,p.value=1, lfc=0, number=30000)
## 仅保存logFC供后面的分析使用
results.lim.all &- results.lim.all[, 1:3]
colnames(results.lim.all) &- c('T1h', 'T24h', 'T7d')
head(results.lim.all, 3)
## & & & & & & & &T1h &T24h &T7d
## 254818_at &0. 6.215
## 245998_at -0. 2.778
## 265119_at -0. 4.380
write.csv(results.lim.all, 'results.lim.all.csv')
4、Session Info
sessionInfo()
## R version 3.1.0 ()
## Platform: x86_64-pc-linux-gnu (64-bit)
## locale:
## &[1] LC_CTYPE=zh_CN.UTF-8 & & & LC_NUMERIC=C && & & & & &
## &[3] LC_TIME=zh_CN.UTF-8 & & &&LC_COLLATE=zh_CN.UTF-8 & &
## &[5] LC_MONETARY=zh_CN.UTF-8 & &LC_MESSAGES=zh_CN.UTF-8&&
## &[7] LC_PAPER=zh_CN.UTF-8 & & & LC_NAME=C & && & & & & &&
## &[9] LC_ADDRESS=C & & & & & & &LC_TELEPHONE=C & & & & & &
## [11] LC_MEASUREMENT=zh_CN.UTF-8 LC_IDENTIFICATION=C & &&&
## attached base packages:
## [1] parallel &tcltk & & stats & & graphics &grDevicesutils & & datasets&
## [8] methods & base & &&
## other attached packages:
## &[1] limma_3.21.1 & & & & simpleaffy_2.41.0 &&gcrma_2.37.0 & & & &
## &[4] genefilter_1.47.0 & &samr_2.0 & & && & & matrixStats_0.8.14 &
## &[7] impute_1.39.0 & & & &ath1121501cdf_2.14.0affy_1.43.0 & & & &&
## [10] Biobase_2.25.0 & & & BiocGenerics_0.11.0&zblog_0.1.0 & & & &&
## [13] knitr_1.5 & & & & &&
## loaded via a namespace (and not attached):
## &[1] affyio_1.33.0 & & & & annotate_1.43.2 && & AnnotationDbi_1.27.3&
## &[4] BiocInstaller_1.15.2 &Biostrings_2.33.3 & &DBI_0.2-7 & & & & & &
## &[7] evaluate_0.5.3 & & & &formatR_0.10 && & & &GenomeInfoDb_1.1.2 &&
## [10] highr_0.3 & & & & & & IRanges_1.99.2& & & &preprocessCore_1.27.0
## [13] R.methodsS3_1.6.1 & & RSQLite_0.11.4 & & &&S4Vectors_0.0.2 & & &
## [16] splines_3.1.0 & & & & stats4_3.1.0 & && & &stringr_0.6.2 & & & &
## [19] survival_2.37-7 & & & tools_3.1.0 & & && & XML_3.98-1.1 & & & &&
## [22] xtable_1.7-3 & & & & &XVector_0.5.3 && & & zlibbioc_1.11.1实验思路：
本文已收录于以下专栏：
相关文章推荐
1、数据读取
1）读取excel文件——GEO中Series Matrix File(s)是预处理过的基因表达矩阵，用excel打开删掉注释信息，获得行为探针，列为样本的基因表达矩阵。
芯片质量分析
芯片数据预处理
获取差异表达基因
GO和KEGG分析
(本文于更新）
“差异”是个统计学概念，获取差异表达基因就要用统计方法，R的统计功能...
source(&http://bioconductor.org/biocLite.R&)
biocLite(&ArrayExpress&)
library(ArrayExpress)
library(affy)%把affy包载入R中
library(tcltk)%把tcltk包载入R中
source(&http://bioconductor.org/biocLite.R&)
biocLite(&ArrayExpress&)
library(ArrayExpress)
芯片质量分析芯片数据预处理获取差异表达基因GO和KEGG分析聚类分析
(本文于更新）
“差异”是个统计学概念，获取差异表达基因就要用统计方法，R的统计功能很强大，适合做这样的事...
芯片质量分析芯片数据预处理获取差异表达基因GO和KEGG分析聚类分析
(本文于更新）
基因芯片技术的特点是使用寡聚核苷酸探针检测基因。前一节使用ReadAffy函数读取...
芯片质量分析芯片数据预处理获取差异表达基因GO和KEGG分析聚类分析
(本文于更新）
TAIR，NASCarray 和 EBI 都有一些公开的免费芯片数据可以下载。本专...
从CEL探针灰度信息文件开始处理，调用R包，进行数据输入、质量控制、背景校正、标准化、汇总、选取差异表达基因、注释、统计分析及可视化，进而将样本分为测试集和验证集，构建Knn、svm分类器，并用留一法...
他的最新文章
讲师：刘文志
讲师：陈伟
您举报文章：
举报原因：
原文地址：
原因补充：
(最多只允许输入30个字)