本文主要参考：

1. //www.greatytc.com/p/ffb3b84c8a4c
2. //www.greatytc.com/p/e07a07a028e4

一. Agilent基因芯片数据：
原理可以参考：https://www.bilibili.com/video/BV1NJ411k7oJ
二. 数据的导入、背景噪音的消除、normalization
课题组自己送到公司测序后返回了rawdata，是荧光图片经软件Feature Extraction version 10.7.3及更新版本处理后得到的基因表达矩阵，每个样品都有单独的一个表达矩阵，仅包含了该样本的基因表达信息。下面就从单独的rawdata的读取开始

#读入原始文件路径
files = dir(path = ".",
               pattern = "txt",
               full.names = T,
               recursive = T)#来到存放rawdata的目录，此目录没有任何其他的文件，因此直接读取当前目录下格式为txt的文件的文件名，得到rawdata的路径

#读入文件
x <- read.maimages(files, source="agilent", green.only=TRUE) #本数据为单通道数据，双通道数据可见文首链接

#背景噪音矫正
data.bg <- backgroundCorrect(x,method = "normexp")#输出Large ElistRaw结果

#标准化：使用 quantile算法：https://www.medsci.cn/article/show_article.do?id=a77e693590b
data.norm <- normalizeBetweenArrays(data.bg,method = "quantile")
#在这里需要注意，data.norm为Large Elist结果，此时表达量为log2形式保存。
#因此此表达量可以认为近似符合正态分布，便于后续的分析。
#基因芯片标准化后的分布没有查到比较有说服力的论证，只看到别人自博客或文章里说，基因芯片表达数据为偏态分布，log2后符合正态分布：
#https://mp.weixin.qq.com/s?__biz=Mzg4MDc2MzUwMg==&mid=2247483703&idx=1&sn=61c9d3ec29d027fe17a5098b30611342&chksm=cf717636f806ff205425bd30fba444758c8b2a3e56bdcf733b9706007612420ea9a35b39e1c0#rd

三. ProbeName与Symbol的转换，参考https://zhuanlan.zhihu.com/p/121415080

#下载探针信息
library(biomaRt)
ensembl <- useMart("ensembl")
ensembl = useDataset("hsapiens_gene_ensembl",mart=ensembl)
searchFilters(mart = ensembl, pattern = "agilent.*")#在库中根据关键词agilent的正则表达式搜索
#我知道所用芯片为8x60k，因此我选择了agilent_sureprint_g3_ge_8x60k_v2（版本高可能包含更多的probe信息）
dfannot <- getBM(attributes=c('agilent_sureprint_g3_ge_8x60k_v2', 'entrezgene_id',
                                'hgnc_symbol'), 
                          values = data.norm$genes, 
                          mart = ensembl)

dfannot <- dfannot[!dfannot$agilent_sureprint_g3_ge_8x60k=="",]#去除空白行

#表达矩阵，并合并相同探针的表达量，使用平均值
data.exp <- as.data.frame(data.norm@.Data[[1]])
data.exp$ProbeName <- data.norm$genes$ProbeName
data.exp <- aggregate(data.exp,by = list(data.exp$ProbeName),FUN = mean)
data.exp <- data.exp[,-8]#去除掉已变为NA的ProbeName列，把名为Group.1的现在的ProbeName列修改为对应的名称

#匹配probename与symobol
probe2gene <- data.norm$genes
#在这里，测序公司返回了一个他们处理过的表达矩阵，其中含有probe与symbol的对应信息，
#因此在下面的代码中，我的probe2gene文件是从公司提供的表达矩阵中提取的，仅含有ProbeName和Symbol两列。
#在使用上面的方法制作probe2gene时最好把不需要的列删除掉方便后续合并。
data.exp <- merge(data.exp,probe2gene,by = "ProbeName",all = F)#这里的probe2gene仅含有ProbeName和Symbol两列

#根据ProbeName合并
data.exp <- merge(data.exp,probe2gene,by = "ProbeName",all = F)

#合并不同探针对应的相同基因，同样使用平均值
summary(duplicated(data.exp$GeneSymbol))#可以看到由很多重复的symbol
exp <- as.data.frame(data.exp[,-1],row.names = data.exp[,1])
exp <- aggregate(exp,by=list(exp$GeneSymbol),FUN = mean)
exp <- as.data.frame(exp[,-c(1,8)],row.names = exp$Group.1)
colnames(exp) <- c("s2","s3","s4","s8","wt1","wt2")

hist(rowMeans(exp),breaks = seq(0,20,by=0.1))#查看表达值分布情况
write.csv(exp,"blood_normalized.csv")#保存

四、limma差异分析：参考jimmy大神

#加载包
library(limma)
library(ggplot2)
library(ggpubr)
library(pheatmap)
library(reshape2)
library(tidyverse)
library(ggrepel)

#导入数据
exp <- read.csv("blood_normalized.csv")
exp <- as.data.frame(exp[,-1],row.names = exp[,1])
exp <- as.matrix(exp)

#set filter
logFCfilter=1 
fdrFilter=0.05

#generate grouplist and design
group_list=c(rep("A",4), rep("C",2))
group_list <- factor(group_list,levels = c("A","C"),ordered = F)#默认排序为字母顺序，logFC中的FC为：后者的平均表达量/前者的平均表达量。
#如果要强制限定顺序的话，ordered=T
design=model.matrix(~group_list)

#deg
fit <- lmFit(exp,design)
fit <- eBayes(fit,trend = T,robust=T)
summary(decideTests(fit))
options(digits = 4)
topTable(fit,coef=2,adjust='BH')

deg=topTable(fit,coef=2,adjust='BH',number = Inf)
diffSig=deg[(deg$adj.P.Val < fdrFilter & (deg$logFC>logFCfilter | deg$logFC<(-logFCfilter))),]
write.csv(diffSig, file="diffsig.csv")
write.csv(deg, file="deg.csv")

#火山图
#为上调下调的top10基因做标记
top10up <- deg[deg$logFC>0,]
top10up <- top10up[order(top10up$P.Value),]
top10up <- top10up[1:10,]
top10down <- deg[deg$logFC<0,]
top10down <- top10down[order(top10down$adj.P.Val),]
top10down <- top10down[1:10,]
deg$label <- rep(NA,nrow(deg))

for(i in 1:nrow(deg)){
  if(rownames(deg)[i]%in% rownames(top10up)|rownames(deg)[i]%in% rownames(top10down)){
    deg[i,"label"] <- rownames(deg)[i]
  }else{deg[i,"label"] <- ""} 
}
#define significance
Significant=ifelse((deg$adj.P.Val<fdrFilter & abs(deg$logFC)>logFCfilter), ifelse(deg$logFC>logFCfilter,"Up","Down"), "Not")

#画图
p = ggplot(deg, aes(logFC, -log10(adj.P.Val)))+
  geom_point(aes(col=Significant))+
  theme_bw()+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+#blank background
  scale_color_manual(values=c("#2f5688", "#BBBBBB", "#CC0000"))+#设定颜色
  labs(title = " ")+
  theme(plot.title = element_text(size = 16, hjust = 0.5, face = "bold"))+
  geom_hline(aes(yintercept=-log10(fdrFilter)),colour="grey",linetype="dashed")+geom_vline(aes(xintercept=logFCfilter), colour="grey",linetype="dashed")+
  geom_vline(aes(xintercept=-logFCfilter), colour="grey",linetype="dashed")#添加阈值辅助线
#添加top10基因的标记
p=p+geom_text_repel(data = deg, 
                    aes(x = logFC,y = -log10(deg$adj.P.Val),label = label),
                    size = 3,box.padding = unit(0.5, "lines"),
                    point.padding = unit(0.8, "lines"), 
                    segment.color = "black",
                    max.overlaps = 20,
                    show.legend = FALSE)
ggsave(p,filename = "vol.pdf")