使用DEseq2做转录组测序差异分析的时候顺便去除批次效应

# 1.构建表达矩阵 ----------------------------------------------------------------
dir.create("./data")

# 魔幻操作，一键清空
rm(list = ls()) 
options(stringsAsFactors = F)
# BiocManager::install('airway')
# 加载airway数据集并转换为表达矩阵
library(airway,quietly = T)
data(airway)
class(airway)

rawcount <- assay(airway)
colnames(rawcount)

# 查看表达谱
rawcount[1:4,1:4]

# 去除前的基因表达矩阵情况
dim(rawcount)

# 获取分组信息
group_list <- colData(airway)$dex
group_list

# 过滤在至少在75%的样本中都有表达的基因
keep <- rowSums(rawcount>0) >= floor(0.75*ncol(rawcount))
table(keep)

filter_count <- rawcount[keep,]
filter_count[1:4,1:4]
dim(filter_count)
 

# 保存表达矩阵和分组结果
save(filter_count, 
     group_list,
     file = "./data/Step01-airwayData.Rdata") 

rm(list = ls())
options(stringsAsFactors = F)

# 读取基因表达矩阵信息
lname <- load(file = "./data/Step01-airwayData.Rdata")
lname 

# 查看分组信息和表达矩阵数据
exprSet <- filter_count
dim(exprSet)
exprSet[1:6,1:6]
table(group_list)

# 加载包
library(DESeq2)
library(pheatmap)
# 每次都要检测数据 
exprSet[1:6,1:6]
table(group_list)
pheatmap(cor(exprSet))
batch=paste0('b',rep(1:2,each=4))
batch
dat=log2(edgeR::cpm(exprSet)+1)

ht_for_RNAcounts <- function(dat){
  dat[1:4,1:4] 
  colnames(dat)
  ac=data.frame(group=group_list,
                batch=batch)
  rownames(ac)=colnames(dat)
  pheatmap(cor(dat),annotation_col  =ac)
}
ht_for_RNAcounts(dat)

ct_with_batch = filter_count
x= sample(1:nrow(ct_with_batch))
y = which(batch=='b1')

ct_with_batch[x,y] = ct_with_batch[x,y]+100
dat=log2(edgeR::cpm(ct_with_batch)+1)
ht_for_RNAcounts(dat)

exprSet=ct_with_batch
suppressMessages(library(DESeq2))  
(colData <- data.frame(row.names=colnames(exprSet), 
                       group_list=group_list,
                       batch= batch ) )
dds <- DESeqDataSetFromMatrix(countData = exprSet,
                              colData = colData,
                              design = ~ group_list+batch )
dds <- DESeq(dds)
resultsNames(dds)
res <- results(dds, name=  "group_list_untrt_vs_trt")
summary(res)
save(res,file = 'group_list_untrt_vs_trt_rm_batch_deg.Rdata')

load(file = 'group_list_untrt_vs_trt_rm_batch_deg.Rdata')
resOrdered <- res[order(res$padj),]
DEG =as.data.frame(resOrdered)
rm_batch_deg = na.omit(DEG)  
head(rm_batch_deg)
# 筛选上下调，设定阈值
fc_cutoff <- 1
fdr <- 0.05 
rm_batch_deg$regulated <- "normal" 
loc_up <- intersect(which(rm_batch_deg$log2FoldChange>log2(fc_cutoff)),
                    which(rm_batch_deg$padj<fdr))
loc_down <- intersect(which(rm_batch_deg$log2FoldChange< (-log2(fc_cutoff))),
                      which(rm_batch_deg$padj<fdr))
rm_batch_deg$regulated[loc_up] <- "up"
rm_batch_deg$regulated[loc_down] <- "down" 
table(rm_batch_deg$regulated)
rm_batch_deg$ENSEMBL = rownames(rm_batch_deg)

# 这个文件Step03-DESeq2_nrDEG.Rdata 是 正常的表达量矩阵的常规差异分析 
load(file = "./data/Step03-DESeq2_nrDEG.Rdata")
head(deg_anno) 
table(deg_anno$regulated) 

tmp= merge(rm_batch_deg,deg_anno,by='ENSEMBL')
library(gplots)
table(tmp$regulated.x,tmp$regulated.y)
balloonplot(table(tmp$regulated.x,tmp$regulated.y)) 

> table(tmp$regulated.x,tmp$regulated.y)
        
          down normal    up
  down       0     36  1066
  normal   242  11899   419
  up       832    188     0

exprSet=ct_with_batch
# 加载包
library(DESeq2) 
# 第一步，构建DESeq2的DESeq对象
colData <- data.frame(row.names=colnames(exprSet),group_list=group_list)
dds <- DESeqDataSetFromMatrix(countData = exprSet,colData = colData,design = ~ group_list)

# 第二步，进行差异表达分析
dds2 <- DESeq(dds)

# 提取差异分析结果，trt组对untrt组的差异分析结果
tmp <- results(dds2,contrast=c("group_list","trt","untrt"))
DEG_DESeq2 <- as.data.frame(tmp[order(tmp$padj),])
head(DEG_DESeq2)

down normal     up 
 74  21718    138 

       down normal    up
  down      73      1     0
  normal  1070  14754  1388
  up         0      0   138