由于实验人员、技术、环境、时间点、芯片处理等各种因素原因，得到的芯片表达矩阵是含有非生物学因素差异的。
有时候想分析两组间差异，但是一个数据集下面的样本数量略少，因此想整合多个样本的数据进行大样本分析。但是在一个数据集里面，即使是统一平台的芯片数据之间都是有批次效应和偏差的，整合在一起就更加放大了。（对于seq数据来说好像影响没这么大，TPM、FPKM等归一方法效果拔群，因此经常有TCGA+GTEx的联合分析）
这里记录芯片数据的批次效应及多数据集整合。参考资料包括：https://www.jianshu.com/p/454fbb13af28、https://www.jianshu.com/p/d4ff373d3d3f、https://www.jianshu.com/p/dce79ecf52ea和https://blog.csdn.net/yayiling/article/details/113664434
有一篇18年的sci report，用的sva来batch normalize，一言难尽，这真的能筛出来吗？anyway，先记录，再比较，看结果。（也许包的功能很好，只是文章让人怀疑）

什么是批次效应？

Leek et. al (2010) define batch effects as follows: Batch effects are sub-groups of measurements that have qualitatively different behaviour across conditions and are unrelated to the biological or scientific variables in a study. For example, batch effects may occur if a subset of experiments was run on Monday and another set on Tuesday, if two technicians were responsible for different subsets of the experiments, or if two different lots of reagents, chips or instruments were used.
批次效应是测量时在条件间有质量差别行为的，在实验中与生物学或科学变量不相关的一小群。举个例子，如果一部分实验在周一进行，另一部分在周二进行，批次效应可能发生。或者两个技术员分别做同一实验的不同部分，或者使用了不同的试剂、芯片或者使用手册都会参数batch effect。

如何判断批次效应？

1、数据集种明确说明
2、Boxplot图
3、聚类树
4、PCA分析
下面代码节选自JIMMY的教程，有兴趣的可以自行学习

## hclust
  colnames(exprSet)=paste(group_list,1:ncol(exprSet),sep='_')
  # Define nodePar
  nodePar <- list(lab.cex = 0.6, pch = c(NA, 19),
                  cex = 0.7, col = "blue")
  hc=hclust(dist(t(exprSet)))
  par(mar=c(5,5,5,10))
  png('hclust.png',res=120)
  plot(as.dendrogram(hc), nodePar = nodePar, horiz = TRUE)
  dev.off()
   ## PCA
  library(ggfortify)
  df=as.data.frame(t(exprSet))
  df$group=group_list
  png('pca.png',res=120)
  autoplot(prcomp( df[,1:(ncol(df)-1)] ), data=df,colour = 'group',frame.type = 'norm')+theme_bw()
  dev.off()

芯片数据的处理

1.批次效应不能被消除，只有尽可能的降低。
2.是否会矫枉过正？
① 联合分析正常组为一个数据集(一种芯片)；
② 实验组为另一个数据集（另一种芯片）；
批次因素和分组因素可能重叠，所以直接对原数据矫正批次可能会抵消一部分真实生物学因素
3.使用removeBatchEffect (limma) 或者ComBat (sva) 函数后得到的表达数据，仅可用于衔接可视化（如聚类、PCA等），可视化展示，不能将去批次后的数据用于差异分析！
问：那我去批次效应整合数据集有什么意义呢？
4.如果想要在鉴定差异基因的过程中降低批次效应，将批次加入到design中

情况1 【数据集(1个)明确给出批次效应】

使用sva包的ComBat函数
有些GEO数据里面包含了实验的批次，就可以用这个包来去除。

情况2 【数据集(1个)压根没提批次效应】

仅判断是否需要log2和normalization
使用limma包中normalizeBetweenArrays函数 （仅在同一个数据集里面使用）

boxplot(rt)
rt <- normalizeBetweenArrays(rt)
boxplot(rt)

问：这些所谓的去除批次效应的方法，和quantile normalization的区别在哪里呢？是否目的是一致的？

情况3 【多数据集合并】

##########################合并数据#########################
#了解到芯片数据每个都平台保留的基因不一样，要全都有的才留下
#setwd('path_to_data_you_download')
files <- unlist(lapply(strsplit(list.files(pattern = 'matrix'),'_'),function(x) {x[[1]]}))
exprslist <- list()
pdatalist <- list()
count <- c()
for (i in files){
  #批量生成数据框
  #assign(paste0(i,'_exprs'),read.table(paste0(i,'_matrix.txt'),sep='\t',row.names = 1))
 #assign(paste0(i,'_exprs'),mutate(get(paste0(i,"_exprs")),ID=rownames(get(paste0(i,"_exprs")))))
 #生成一个列表
  rt <- read.table(paste0(i,'_matrix.txt'),sep='\t',row.names = 1)
  count <- c(count,ncol(rt))
  rt <- as.data.frame(normalizeBetweenArrays(rt))
  rt$ID <- rownames(rt)
  exprslist[[i]] <- rt
  #assign(paste0(i,'_pdata'),read.delim(paste0(i,'_pdata.txt'),sep='\t',row.names = NULL))
  rt <- read.delim(paste0(i,'_pdata.txt'),sep='\t',row.names = NULL)
  pdatalist[[i]] <- rt
  }
rm(rt)
in_jo <- function(x,y){
  inner_join(x,y,by='ID')
}
merged_exprs <- Reduce(in_jo,exprslist)
rownames(merged_exprs) <- merged_exprs$ID
merged_exprs <- select(merged_exprs,-ID)
#有些行里面包含了na，需要去掉
merged_exprs <- merged_exprs[complete.cases(merged_exprs),]
Listrb <- function(x,y){
  rbind(x,y)
}
merged_pdata<- Reduce(Listrb,pdatalist)
merged_pdata$group <- ifelse(grepl('Primary',merged_pdata$title),'Primary','Metastasis')
rm(exprslist,pdatalist,i,in_jo,Listrb)
#最终我们得到了合并后的exprs和pdata，可以用来batch了,count包含每个batch的数量   

（1）可去除：使用sva包/limma包；

①limma包removeBatchEffect函数

library(limma)
boxplot(merged_exprs)
batch <- rep(files,times=count)
batch <- as.factor(batch)
rm(files,count)
design <- model.matrix(~0 + batch)
batchremove_limma <- removeBatchEffect(merged_exprs,
                                batch = batch)
boxplot(batchremove_limma)

library(ggfortify)
limma <- as.data.frame(t(batchremove_limma))
limma$group <- merged_pdata$group
limma$batch <- batch
autoplot(prcomp(limma[,1:(ncol(originaldata)-2)] ), 
         data=limma,colour = 'group',
         frame.type = 'norm')+
  theme_bw()
autoplot(prcomp(limma[,1:(ncol(originaldata)-2)] ), 
         data=limma,colour = 'batch',
         frame.type = 'norm')+
  theme_bw()

②sva包ComBat函数

library(sva)
batchremove_combat <- ComBat(dat = as.matrix(merged_exprs), batch = batch)
boxplot(batchremove_combat)

③MatchMixeR

坑

（2）不可去除：RobustRankAggreg包整合分析

RobustRankAggreg不再需要合并原始数据，只需要按照流程得到每个数据集的差异表达基因，按照logFC从大到小的顺序排列好即可。

padj=0.05
logFC=1

if(TRUE){
files=list.files()
upList=list()
downList=list()
allFCList=list()

for(i in 1:length(files)){
  inputFile=files[i]
  rt=read.table(inputFile,header=T,sep = '\t') # 注意文件读取
  rt <- rt[order(rt$logFC),]
  header=unlist(strsplit(unlist(strsplit(inputFile,"_"))[[2]],".txt"))# 新数据需要修改
  downList[[header[1]]]=as.vector(rt[,1])
  upList[[header[1]]]=rev(as.vector(rt[,1]))
  fcCol=rt[,1:2]
  colnames(fcCol)=c("Gene",header[[1]])
  allFCList[[header[1]]]=fcCol
}

rm(fcCol,rt,files,header,i,inputFile)

mergeLe=function(x,y){
  merge(x,y,by="Gene",all=T)}
newTab=Reduce(mergeLe,allFCList)
#newTab = na.omit(newTab) 
newTab <- newTab[!duplicated(newTab$Gene),]
rownames(newTab)=newTab[,1]
newTab=newTab[,2:ncol(newTab)]
newTab[is.na(newTab)]=0
}

library(RobustRankAggreg)
if(TRUE){
upMatrix = rankMatrix(upList)
#upMatrix = rankMatrix(upList,full=T)

upAR = aggregateRanks(rmat=upMatrix)
colnames(upAR)=c("Name","Pvalue")
upAdj=p.adjust(upAR$Pvalue,method="bonferroni")#
upXls=cbind(upAR,adjPvalue=upAdj)
upFC=newTab[as.vector(upXls[,1]),]
upXls=cbind(upXls,logFC=rowMeans(upFC))

#write.table(upXls,file="up.xls",sep="\t",quote=F,row.names=F)
upSig=upXls[(upXls$adjPvalue<padj & upXls$logFC>logFC),]
#upSig=upXls[ upXls$logFC>logFC,]
#write.table(upSig,file="upSig.xls",sep="\t",quote=F,row.names=F)

downMatrix = rankMatrix(downList)
#downMatrix = rankMatrix(downList,full=T)

downAR = aggregateRanks(rmat=downMatrix)
colnames(downAR)=c("Name","Pvalue")
downAdj=p.adjust(downAR$Pvalue,method="bonferroni")
downXls=cbind(downAR,adjPvalue=downAdj)
downFC=newTab[as.vector(downXls[,1]),]
downXls=cbind(downXls,logFC=rowMeans(downFC))
#write.table(downXls,file="down.xls",sep="\t",quote=F,row.names=F)

downSig=downXls[(downXls$adjPvalue<padj & downXls$logFC< -logFC),]
#downSig=downXls[downXls$logFC< -logFC,]

#write.table(downSig,file="downSig.xls",sep="\t",quote=F,row.names=F)

allSig = rbind(upSig,downSig)
colnames(allSig)
allSig = allSig[,c("Name","logFC")]
#write.table(allSig,file = 'allSign.xls',sep = '\t',quote = F)
}

hminput=newTab[c(as.vector(upSig[1:20,1]),as.vector(downSig[1:20,1])),]

library(pheatmap)
pheatmap(hminput,display_numbers = TRUE,
         fontsize_row=10,
         fontsize_col=12,
         color = colorRampPalette(c("green", "white", "red"))(50),
         cluster_cols = FALSE,cluster_rows = FALSE, )

write.table(downSig,sep='\t',quote=F,row.names = F,file='downSig.txt')
write.table(upSig,sep='\t',quote=F,row.names = F,file='upSig.txt')

总结

本文以代码形式提供了整合芯片数据获得差异基因的几种方法，具体的算法流程需要各位亲自理解和摸索（因为我也不会）