Differential Expression of RNA-seq data
本文将要介绍的是在R中进行RNA-seq 数据基因表达差异分析的实战代码.
原文地址：https://bioinformatics-core-shared-training.github.io/RNAseq-R/rna-seq-de.nb.html
书接上文（用R也可以完成的RNA-Seq下游分析-2）
在预处理生成了表达矩阵且标准化后，接下来我们要做的就是差异分析了。
本次流程需要的包
同时载入上次预处理完的数据preprocessing.Rdata
library(edgeR)
library(limma)
library(Glimma)
library(gplots)
library(org.Mm.eg.db)
load("preprocessing.Rdata")
在edgeR包的分析流程中，构建了表达矩阵后，便需要使用model.matrix函数创建保存有分组信息的矩阵design matrix。该矩阵以的方式存储了分组信息。本次分析中的分组信息包括了：小鼠状态与细胞类型，现在我们先在两种因素并不存在相互作用的假设下拟合线性模型。
> group <- as.character(group)
> group
[1] "basal.virgin" "basal.virgin" "basal.pregnant" "basal.pregnant" "basal.lactate"
[6] "basal.lactate" "luminal.virgin" "luminal.virgin" "luminal.pregnant" "luminal.pregnant"
[11] "luminal.lactate" "luminal.lactate"
group的分组信息包括了状态和细胞类型，因此我们使用strsplit函数，以‘.’为分隔把细胞类型信息和状态信息分别提取出来。
> type <- sapply(strsplit(group, ".", fixed=T), function(x) x[1])
> status <- sapply(strsplit(group, ".", fixed=T), function(x) x[2])
> type
[1] "basal" "basal" "basal" "basal" "basal" "basal" "luminal" "luminal" "luminal"
[10] "luminal" "luminal" "luminal"
> status
[1] "virgin" "virgin" "pregnant" "pregnant" "lactate" "lactate" "virgin" "virgin"
[9] "pregnant" "pregnant" "lactate" "lactate"
成功提取后，构建分组信息矩阵
> design <- model.matrix(~ type + status)
> design
(Intercept) typeluminal statuspregnant statusvirgin
1 1 0 0 1
2 1 0 0 1
3 1 0 1 0
4 1 0 1 0
5 1 0 0 0
6 1 0 0 0
7 1 1 0 1
8 1 1 0 1
9 1 1 1 0
10 1 1 1 0
11 1 1 0 0
12 1 1 0 0
attr(,"assign")
[1] 0 1 2 2
attr(,"contrasts")
attr(,"contrasts")$type
[1] "contr.treatment"
attr(,"contrasts")$status
[1] "contr.treatment"
Estimating the dispersion
> dgeObj <- estimateCommonDisp(dgeObj)
#Then we estimate gene-wise dispersion estimates, allowing a possible trend with averge count size:
> dgeObj <- estimateGLMTrendedDisp(dgeObj)
> dgeObj <- estimateTagwiseDisp(dgeObj)
#Plot the estimated dispersions
> plotBCV(dgeObj)
Testing for differential expression
> fit <- glmFit(dgeObj, design)
# 看一眼系数
> head(coef(fit))
(Intercept) typeluminal statuspregnant statusvirgin
497097 -11.187922 -7.58804851 -0.7085514 -0.09305118
20671 -12.715063 -1.85287334 0.2269001 0.49554506
27395 -11.221391 0.56368066 -0.1415910 -0.29221577
18777 -10.146793 0.08280255 -0.1845489 -0.48795441
21399 -9.909825 -0.24195503 0.1753606 0.13494615
58175 -16.310131 3.09936215 1.1975518 0.84742701
# Conduct likelihood ratio tests for luminal vs basal and show the top genes:
> lrt.BvsL <- glmLRT(fit, coef=2)
# 查看前10个最显著的差异表达基因
> topTags(lrt.BvsL)
Coefficient: typeluminal
logFC logCPM LR PValue FDR
110308 -8.940579 10.264297 24.89789 6.044844e-07 0.0004377961
50916 -8.636503 5.749781 24.80037 6.358512e-07 0.0004377961
12293 -8.362247 6.794788 24.68526 6.749827e-07 0.0004377961
56069 -8.419433 6.124377 24.41532 7.764861e-07 0.0004377961
24117 -9.290691 6.757163 24.32506 8.137331e-07 0.0004377961
12818 -8.216790 8.172247 24.24233 8.494462e-07 0.0004377961
22061 -8.034712 7.255370 24.16987 8.820135e-07 0.0004377961
12797 -9.001419 9.910795 24.12854 9.011487e-07 0.0004377961
50706 -7.697022 10.809629 24.06926 9.293193e-07 0.0004377961
237979 -8.167451 5.215921 24.03528 9.458678e-07 0.0004377961
最后，保存差异分析的结果（实际分析时可不保存，可以直接代码分析）
save(lrt.BvsL,dgeObj,group,file="DE.Rdata")
