目前，转录组学差异分析的三种常用流程：DEseq2、edgeR、limma。
在具体选择哪一个的问题上，有几个点需要明确：
在进行差异分析之前，首先明确你处理的是哪种测序数据。芯片测序只能用
limma包分析，高通量测序（RNAseq）三者通吃。此外，
imma包在生物信息学分析中，还有一个很重要的功能是处理批次效应,特别是进行TCGA与GTExs数据合并edgeR包有个优势，就是可以处理无生物学重复的样本，而DEseq2、limma不可。
一.limma
1. GSE42872数据获取
<pre class="md-fences md-end-block md-fences-with-lineno ty-contain-cm modeLoaded" spellcheck="false" lang="R" cid="n17" mdtype="fences" style="box-sizing: border-box; overflow: visible; font-family: Monaco, Consolas, "Andale Mono", "DejaVu Sans Mono", monospace; margin-top: 0px; margin-bottom: 20px; font-size: 0.9rem; display: block; break-inside: avoid; text-align: left; white-space: normal; background: rgb(51, 51, 51); position: relative !important; padding: 10px 10px 10px 0px; width: inherit; color: rgb(184, 191, 198); font-style: normal; font-variant-ligatures: normal; font-variant-caps: normal; font-weight: 400; letter-spacing: normal; orphans: 2; text-indent: 0px; text-transform: none; widows: 2; word-spacing: 0px; -webkit-text-stroke-width: 0px; text-decoration-thickness: initial; text-decoration-style: initial; text-decoration-color: initial;"> ### 1.读入数据(以数据GSE42872为例)
#方法1：直接借助getGEO函数读取
library(GEOquery)
library(limma)
exprSet<-getGEO('GSE42872',destdir='.',
AnnotGPL=F,
getGPL=F)
#方法2：在GEO官网界面，下载数据Series Matrix File(s)，然后本地读取
exprSet <- read.table("data/GSE42872_series_matrix.txt",
comment.char="!",
stringsAsFactors=F,
header=T)
rownames(exprSet) <- exprSet[,1]
exprSet <- exprSet[,-1]
### 2.数据预处理，探针ID转换，探针去重
#注意此时数据的格式是矩阵
ex <- exprSet
qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
LogC <- (qx[5] > 100) ||
(qx[6]-qx[1] > 50 && qx[2] > 0) ||
(qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2)
# 开始判断
if (LogC) {
ex[which(ex <= 0)] <- NaN
## 取log2
exprSet <- log2(ex)
print("log2 transform finished")
}else{
print("log2 transform not needed")
}
boxplot(exprSet,outline=FALSE, notch=T, las=2)
# 该函数默认使用quntile 矫正差异
exprSet=normalizeBetweenArrays(exprSet)
boxplot(exprSet,outline=FALSE, notch=T, las=2)
exprSet <- as.data.frame(exprSet)
### 3.探针基因名转换
### platformMap 中有常见的平台个R注释包的对应关系，这是果子老师整理的。当然，平台对应的R注释包自行检索，也不麻烦。
platformMap <- data.table::fread("resource/platformMap.txt",data.table = F)
index <- "GPL6244" #平台的名称在GSE42872中查看
paste0(platformMap$bioc_package[grep(index,platformMap$gpl)],".db") #输出结果即为R注释包
library(hugene10sttranscriptcluster.db) #hugene10sttranscriptcluster.db为上一步所得的结果
probe2symbol_df <- toTable(get("hugene10sttranscriptclusterSYMBOL")) #probe2symbol_df即为包含探针和基因名的数据
### 4.探针转换以及去重
library(dplyr)
library(tibble)
exprSet <- exprSet %>%
rownames_to_column("probe_id") %>%
inner_join(probe2symbol_df,by="probe_id") %>%
select(-probe_id) %>%
select(symbol,everything()) %>%
mutate(rowMean =rowMeans(.[,-1])) %>%
arrange(desc(rowMean)) %>%
distinct(symbol,.keep_all = T) %>%
select(-rowMean) %>%
column_to_rownames("symbol")
save(exprSet,file = "output/exprSet_rmdup.Rdata")
2. limma包进行差异分析流程
### 1.创建分组
group <- c(rep("con",3),rep("treat",3))
group <- factor(group,levels = c("con","treat"))
### 2.构建比较矩阵
design <- model.matrix(~0+group)
colnames(design) <- levels(group)
design
contrast.matrix <- makeContrasts(treat - con,levels=design)
contrast.matrix
### 3.接下来是一顿常规操作
fit <- lmFit(exprSet,design)#线性模型拟合
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit)#贝叶斯检验
### 4.最终得到差异分析结果
allDiff=topTable(fit2,adjust='fdr',coef=1,number=Inf)</pre>
二.DESeq
### 1.读入表达矩阵文件 （数据链接：https://pan.baidu.com/s/1MmO2_oPMkqYZjoulna527Q 提取码：djft）
exprSet <- data.table::fread("data/exprSet_counts.csv",data.table = F)
### 2.读入样本信息文件（数据链接：https://pan.baidu.com/s/1hGfBDTnD11_MJrNU4WxUTw 提取码：ci1h）
metadata <- data.table::fread("data/metadata.txt",data.table = F,header = F)
rownames(metadata) <- metadata[,1]
metadata <- metadata[colnames(exprSet)[-1],]
# 添加分组信息
metadata$group <- ifelse(grepl("ESR1",metadata$V2),"treat","control")
metadata = metadata[,c(1,3)]
colnames(metadata) <- c("sample","group")
### 3.核心环节，构建dds对象，前面的操作都是铺垫
library(DESeq2)
dds <-DESeqDataSetFromMatrix(countData=exprSet,
colData=metadata,
design=~group,
tidy=TRUE)# 第一列如果是基因名称，需要自动处理，设置参数tidy=TRUE
nrow(dds)
rownames(dds)
# 筛选样本，counts函数提取表达量数据
dds <- dds[rowSums(counts(dds))>1,]
nrow(dds)
### 4.数据质量判断
vsd <- vst(dds, blind = FALSE) # vst标准化处理
plotPCA(vsd, "group")# 内置函数plotPCA进行主成分分析画图
### 5.导出标准化后的表达数据
exprSet_vst <- as.data.frame(assay(vsd))# assay函数提取vst标准化后的数据
