comprehensive/intergrated/meta express data analysis

通过整合多个数据集来构建大数据基础，进行深入的数据挖掘

推荐文献：Transcriptomic profiling of skeletal muscle adaptations to exercise and inactivity

当使用多数据集的时候，注意经过两个处理步骤：

前处理：归一化；背景校正等（具体见上一篇差异基因分析，使用affy包的justMA函数进行的CEL格式原始数据的预处理）前处理：去除批次效应（批间差，batch effects）之后再进行差异表达分析（limma包）

批次效应：不同人员，不同环境，不同时间，不同机器以及不同批次都会导致batch effects

去除批次效应常用的R包"sva" ；常用函数combat()

sva将整个实验的变量分为adjustment variables和variables of interest,我简称为待修正变量（非关注变量）和关注变量。其中关注变量是实验设计中关注的差异效应（CK-VS-T），而非关注变量就是其他无关的，但是在实验过程引入的变量，例如实验的批次，试剂的批次，人员的批次，测序的批次，甚至包括材料本身的批次（背景噪音）。

文档来源：

处理批次效应连续剧第三集（R包基本OK）​mp.weixin.qq.com/s/ssHlEUNOGxmiiWL1zCTxWQComBat function - RDocumentation​www.rdocumentation.org/packages/sva/versions/3.20.0/topics/ComBat

推荐文献：

Removing Batch Effects in Analysis of Expression Microarray Data: An Evaluation of Six Batch Adjustment Methods

实操演示：

开始前准备三份文件：

数据原始文件（CEL格式，把所有需要的数据集样本文件放在一个文件夹下）Annotation文件（芯片的注释信息文件）Target文件（样本分组信息文件，包括三列：样本名（GSMxxx）、分组信息（control / Test）、批次信息（数据集名GSExxx，几个数据集就有几种批次信息）

实操代码：

##加载affy包和limma包

library(affy)

library(limma)

##设置工作路径

getwd()

setwd()

##import phenotype data 输入样本的分组信息（文件3）

phenoData = read.AnnotatedDataFrame(data/target.txt)

pheno = pData(phenoData)

View(pheno)

#import Annotaion 输入数据集的注释文件（文件2）

anno = read.csv("./data/Annotation.csv",head=T)

View(head(anno))

##第一步：RMA normalization 使用RMA函数进行标准化

#读取CEL格式原始文件，#eset.rma = RMA(Data)

eset.rma <- justRMA(filenames=paste(rownames(pheno),.CEL,sep=), celfile.path=./data/CEL/)

#标准化以后的datExpr

datExpr = exprs(eset.rma)

#缺失值处理（利用impute函数）

library(impute)

#KNN法计算缺失值（不同的数据可能需要不同的算法计算）

imputed_gene_exp = impute.knn(datExpr,k=10,rowmax = 0.5,

colmax=0.8,maxp =3000, rng.seed=362436069)

#补充缺失值以后的datExpr2

datExpr2 = imputed_gene_exp$data

##第二步：去除批间差（batch effects）

#sva--combat

if (!requireNamespace("BiocManager", quietly = TRUE))

install.packages("BiocManager")

BiocManager::install("sva")

library(sva)

library(pamr)

batch = pheno[,c(batch)]

View(batch)

modcombat = model.matrix(~1,data=pheno)

##查询combat参数help("ComBat")

combat_edata = ComBat(dat=datExpr2,

batch=batch,

mod=modcombat,

par.prior=T,

prior.plot=T) ##datExpr2为补充缺失值以后的数据；batch为定义的批次

write.table(combat_edata,file="Expdata.txt",sep="\t")

##第三步：limma函数

#样本分组

Group = factor(pheno$group,levels=c(Responder,Nonresponder))

design = model.matrix(~0+Group)

colnames(design) <- c(Responder,Nonresponder)

design

#线性模型拟合

fit <- lmFit(combat_edata, design)

#构建比对模型，比较两个条件下的表达数据

contrast.matrix <- makeContrasts(Responder-Nonresponder,

levels=design)

#####################

####################

install.packages("openxlsx") # Install just once

install.packages("dplyr")

library(openxlsx) # Load on each new session

library(futile.logger)

#比对模型进行差值计算

fit2 <- contrasts.fit(fit, contrast.matrix)

#贝叶斯检验

fit2 <- eBayes(fit2)

#找出差异基因并输出符合阈值的结果

diff = topTable(fit2,adjust.method="fdr",coef=1,p.value=0.01,

lfc=log(1.5,2),number=5000,Future home of sort.by = logFC)

diff

#ID转换

if(nrow(diff) == 0){next}

diff$Gene = anno$Gene.Symbol[match(rownames(diff),anno$ID_REF)]

diff$ID_REF = rownames(diff)

diff = diff[,c(8,7,1:6)]

diff = diff[diff$Gene != ---,]

View(diff)

##输出表格

write.xlsx(diff,

DEG.xlsx,

sheetName=colnames(contrast.matrix)[1],

col.names=T,

row.names=F,

append=T)

日期：2021/4/30