Day1—差异分析，我该如何选择？limma DESeq edge ?

目前，转录组学差异分析的三种常用流程：DEseq2、edgeR、limma。

在具体选择哪一个的问题上，有几个点需要明确：

1. 在进行差异分析之前，首先明确你处理的是哪种测序数据。芯片测序只能用limma包分析，高通量测序（RNAseq）三者通吃。

2. 此外，imma包在生物信息学分析中，还有一个很重要的功能是处理批次效应,特别是进行TCGA与GTExs数据合并

3. edgeR包有个优势，就是可以处理无生物学重复的样本，而DEseq2、limma不可。

一.limma

1.  GSE42872数据获取

 ### 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)

二.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标准化后的数据