Skip to content

This issue was moved to a discussion.

You can continue the conversation there. Go to discussion →

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

基因芯片基因差异表达分析流程示例与讨论 #10

Closed
ShixiangWang opened this issue Aug 15, 2018 · 1 comment
Closed

基因芯片基因差异表达分析流程示例与讨论 #10

ShixiangWang opened this issue Aug 15, 2018 · 1 comment
Assignees
@ShixiangWang

Comments

@ShixiangWang
Copy link
Owner

基因芯片的差异表达分析主要有 构建基因表达矩阵、构建实验设计矩阵、构建对比模型（对比矩阵）、线性模型拟合、贝叶斯检验和生成结果报表 六个关键步骤。

下面是模拟的一个示例：

# Simulate gene expression data for 100 probes and 6 microarraexprSets
# MicroarraexprSet are in two groups
# First two probes are differentiallexprSet expressed in second group
# Std deviations varexprSet between genes with prior df=4

# 构建模拟的表达矩阵，实际处理时换成自己的表达矩阵即可
sd <- 0.3*sqrt(4/rchisq(100,df=4))
exprSet <- matrix(rnorm(100*6,sd=sd),100,6)
rownames(exprSet) <- paste("Gene",1:100)
colnames(exprSet) <- c(paste0("con-",1:3), paste0("G3-",1:3))
exprSet[1:2,4:6] <- exprSet[1:2,4:6] + 2

library(limma)
# 构建实验设计矩阵
group_list = c(rep("con",3), rep("G3",3))
# 这里根据实际的情况设置（表型）分组，对应表达矩阵的列：样本

design <- model.matrix(~0+factor(group_list))
design

colnames(design)=levels(factor(group_list))
rownames(design)=colnames(exprSet)
design
# 实验设计矩阵的每一行对应一个样品的编码，
# 每一列对应样品的一个特征。这里只考虑了一个因素两个水平，
# 如果是多因素和多水平的实验设计，会产生更多的特征，需要参考文档设计。

# 构建对比模型，比较两个实验条件下表达数据
contrast.matrix<-makeContrasts(G3-con,levels = design)
#contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)
contrast.matrix ##这个矩阵声明，我们要把G3组跟con组进行差异分析比较


##### 差异分析
##step1 线性模型拟合
fit <- lmFit(exprSet,design)
##step2 根据对比模型进行差值计算 
fit2 <- contrasts.fit(fit, contrast.matrix) 
##step3 贝叶斯检验
fit2 <- eBayes(fit2) 
##step4 生成所有基因的检验结果报告
tempOutput = topTable(fit2, coef=1, n=Inf)
##step5 用P.Value进行筛选，得到全部差异表达基因
dif <- tempOutput[tempOutput[, "P.Value"]<0.01,]
# 显示一部分报告结果
head(dif)

参考：

更新资料：

相关问题请下面留言讨论。
@ShixiangWang ShixiangWang self-assigned this Aug 15, 2018
@ShixiangWang
Copy link
Owner Author

下面的知识点有点重要，我copy一下jimmy的 差异分析是否需要比较矩阵

差异分析是否需要比较矩阵

最流行的差异分析软件就是limma了，它现在更新了一个voom的算法，所以既可以对芯片数据，也可以对转录组高通量测序数据进行分析，其它所有的差异分析软件其实都是模仿这个的。

我以前讲到过做差异分析，需要三个数据：

  • 表达矩阵

  • 分组矩阵

  • 差异比较矩阵

前面两个肯定是必须的，有表达矩阵，样本必须进行分组，才能分析，但是我看到过好几种例子，有的有差异比较矩阵，有的没有。

后来我仔细研究了一下limma包的说明书，发现这其实是一个很简单的问题。

大家仔细观察下面的两个代码

首先是不需要差异比较矩阵的

    library(CLL)
    data(sCLLex)
    library(limma)
    design=model.matrix(~factor(sCLLex$Disease))
    fit=lmFit(sCLLex,design)
    fit=eBayes(fit)
    options(digits = 4)
    #topTable(fit,coef=2,adjust='BH') 
    > topTable(fit,coef=2,adjust='BH')
               logFC AveExpr      t   P.Value adj.P.Val     B
    39400_at  1.0285   5.621  5.836 8.341e-06   0.03344 3.234
    36131_at -0.9888   9.954 -5.772 9.668e-06   0.03344 3.117
    33791_at -1.8302   6.951 -5.736 1.049e-05   0.03344 3.052
    1303_at   1.3836   4.463  5.732 1.060e-05   0.03344 3.044
    36122_at -0.7801   7.260 -5.141 4.206e-05   0.10619 1.935
    36939_at -2.5472   6.915 -5.038 5.362e-05   0.11283 1.737
    41398_at  0.5187   7.602  4.879 7.824e-05   0.11520 1.428
    32599_at  0.8544   5.746  4.859 8.207e-05   0.11520 1.389
    36129_at  0.9161   8.209  4.859 8.212e-05   0.11520 1.389
    37636_at -1.6868   5.697 -4.804 9.355e-05   0.11811 1.282

然后是需要差异比较矩阵的

    library(CLL)
    data(sCLLex)
    library(limma)
    design=model.matrix(~0+factor(sCLLex$Disease))
    colnames(design)=c('progres','stable')
    fit=lmFit(sCLLex,design)
    cont.matrix=makeContrasts('progres-stable',levels = design)
    fit2=contrasts.fit(fit,cont.matrix)
    fit2=eBayes(fit2)
    options(digits = 4)
    topTable(fit2,adjust='BH')
     
               logFC AveExpr      t   P.Value adj.P.Val     B
    39400_at -1.0285   5.621 -5.836 8.341e-06   0.03344 3.234
    36131_at  0.9888   9.954  5.772 9.668e-06   0.03344 3.117
    33791_at  1.8302   6.951  5.736 1.049e-05   0.03344 3.052
    1303_at  -1.3836   4.463 -5.732 1.060e-05   0.03344 3.044
    36122_at  0.7801   7.260  5.141 4.206e-05   0.10619 1.935
    36939_at  2.5472   6.915  5.038 5.362e-05   0.11283 1.737
    41398_at -0.5187   7.602 -4.879 7.824e-05   0.11520 1.428
    32599_at -0.8544   5.746 -4.859 8.207e-05   0.11520 1.389
    36129_at -0.9161   8.209 -4.859 8.212e-05   0.11520 1.389
    37636_at  1.6868   5.697  4.804 9.355e-05   0.11811 1.282

大家运行一下这些代码就知道，两者结果是一模一样的。

而差异比较矩阵的需要与否，主要看分组矩阵如何制作的！

design=model.matrix(~factor(sCLLex$Disease))

design=model.matrix(~0+factor(sCLLex$Disease))

有本质的区别！！！

前面那种方法已经把需要比较的组做出到了一列，需要比较多次，就有多少列，第一列是截距不需要考虑，第二列开始往后用coef这个参数可以把差异分析结果一个个提取出来。

而后面那种方法，仅仅是分组而已，组之间需要如何比较，需要自己再制作差异比较矩阵，通过makeContrasts函数来控制如何比较！

@ShixiangWang ShixiangWang pinned this issue Jan 11, 2019
@ShixiangWang ShixiangWang transferred this issue from ShixiangWang/MessageBoard Jan 27, 2021
Repository owner locked and limited conversation to collaborators Jan 27, 2021

This issue was moved to a discussion.

You can continue the conversation there. Go to discussion →

Assignees

@ShixiangWang ShixiangWang

Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
1 participant
@ShixiangWang