轉(zhuǎn)自果子學(xué)生信
1.加載R包獲取數(shù)據(jù)

library(GEOquery)
gset = getGEO('GSE32575', destdir=".",getGPL = F)
gset=gset[[1]]

2.通過(guò)pData函數(shù)獲取分組信息

pdata=pData(gset)
group_list=c(rep('before',18),rep('after',18))
group_list=factor(group_list)
## 強(qiáng)制限定順序
group_list <- relevel(group_list, ref="before")

3.通過(guò)exprs函數(shù)獲取表達(dá)矩陣并校正

exprSet=exprs(gset)
##查看整體樣本的表達(dá)情況
boxplot(exprSet[,-c(1:12)],outline=FALSE, notch=T,col=group_list, las=2)
##整體表達(dá)不整齊项钮，使用limma包內(nèi)置函數(shù)人工校正
library(limma) 
exprSet=normalizeBetweenArrays(exprSet)
boxplot(exprSet,outline=FALSE, notch=T,col=group_list, las=2)

4.判斷是否需要進(jìn)行數(shù)據(jù)轉(zhuǎn)換

##根據(jù)分組信息绒极，去除前12個(gè)樣本
exprSet = as.data.frame(exprSet)[,-seq(1,12)]
##表達(dá)量很大燕酷，log轉(zhuǎn)換(選log2)
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
exprSet <- log2(ex)
print("log2 transform finished")}else{print("log2 transform not needed")}

5.注釋基因

##導(dǎo)入R包和平臺(tái)的注釋信息對(duì)應(yīng)關(guān)系表格 platformMap
platformMap <- data.table::fread("platformMap.txt")
##獲取注釋平臺(tái)
index = gset@annotation
##使用代碼自動(dòng)獲取對(duì)應(yīng)注釋包
platformDB = paste0(platformMap$bioc_package[grep(index,platformMap$gpl)],".db")
##安裝、加載包
if(length(getOption("BioC_mirror"))==0) options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
if(!require("illuminaHumanv2.db")) BiocManager::install("illuminaHumanv2.db",update = F,ask = F)
library(illuminaHumanv2.db)
##獲取探針對(duì)應(yīng)的symbol信息
probeset <- rownames(exprSet)
## 使用lookup函數(shù)，找到探針在illuminaHumanv2.db中的對(duì)應(yīng)基因名稱
SYMBOL <-  annotate::lookUp(probeset,"illuminaHumanv2.db", "SYMBOL")
## 轉(zhuǎn)換為向量
symbol = as.vector(unlist(SYMBOL))
##制作probe2symbol轉(zhuǎn)換文件
probe2symbol = data.frame(probeset,symbol,stringsAsFactors = F)

6.探針轉(zhuǎn)換與基因去重

library(dplyr)
library(tibble)
exprSet <- exprSet %>% 
  rownames_to_column(var="probeset") %>% 
  #合并探針的信息
  inner_join(probe2symbol,by="probeset") %>% 
  #去掉多余信息
  select(-probeset) %>% 
  #重新排列
  select(symbol,everything()) %>% 
  #求出平均數(shù)(這邊的點(diǎn)號(hào)代表上一步產(chǎn)出的數(shù)據(jù))
  mutate(rowMean =rowMeans(.[grep("GSM", names(.))])) %>% 
  #去除symbol中的NA
  filter(symbol != "NA") %>% 
  #把表達(dá)量的平均值按從大到小排序
  arrange(desc(rowMean)) %>% 
  # symbol留下第一個(gè)
  distinct(symbol,.keep_all = T) %>% 
  #反向選擇去除rowMean這一列
  select(-rowMean) %>% 
  # 列名變成行名
  column_to_rownames(var = "symbol")
現(xiàn)在數(shù)據(jù)變成這個(gè)樣子

7.差異分析

##如果沒(méi)有配對(duì)信息
design=model.matrix(~ group_list)
fit=lmFit(exprSet,design)
fit=eBayes(fit) 
allDiff=topTable(fit,adjust='fdr',coef="group_listafter",number=Inf,p.value=0.05) 
##如果有配對(duì)信息
pairinfo = factor(rep(1:18,2))
design=model.matrix(~ pairinfo+group_list)
fit=lmFit(exprSet,design)
fit=eBayes(fit) 
allDiff_pair=topTable(fit,adjust='BH',coef="group_listafter",number=Inf,p.value=0.05)

分析結(jié)果的各列數(shù)據(jù)含義：
　　“l(fā)ogFC”是兩組表達(dá)值之間以2為底對(duì)數(shù)化的的變化倍數(shù)锥债，一般表達(dá)相差2倍以上有意義岸售；
　　“AveExpr”是該探針組所在所有樣品中的平均表達(dá)值；
　　“t”是貝葉斯調(diào)整后T 檢驗(yàn)的 t 值翰灾；
　　“P.Value”是貝葉斯檢驗(yàn)的 P 值缕粹；
　　“adj.P.Val”是調(diào)整后的 P 值，更有參考價(jià)值纸淮；
　　“B”是經(jīng)驗(yàn)貝葉斯得到的標(biāo)準(zhǔn)差的對(duì)數(shù)化值平斩。

8.作圖驗(yàn)證(非必要)
轉(zhuǎn)換為ggplot2喜歡的數(shù)據(jù)格式，行是觀測(cè)咽块，列是變量绘面，即清潔數(shù)據(jù)

data_plot = as.data.frame(t(exprSet))
data_plot = data.frame(pairinfo=rep(1:18,2),
                       group=group_list,
                       data_plot,stringsAsFactors = F)

以CAMKK2為例做配對(duì)圖

library(ggplot2)
ggplot(data_plot, aes(group,CH25H,fill=group)) +
  geom_boxplot() +
  geom_point(size=2, alpha=0.5) +
  geom_line(aes(group=pairinfo), colour="black", linetype="11") +
  xlab("") +
  ylab(paste("Expression of ","CH25H"))+
  theme_classic()+
  theme(legend.position = "none")

批量畫出差異最大的8個(gè)基因

library(dplyr)
library(tibble)
allDiff_arrange <- allDiff_pair %>% 
  rownames_to_column(var="genesymbol") %>% 
  arrange(desc(abs(logFC)))
genes <- allDiff_arrange$genesymbol[1:8]

plotlist <- lapply(genes, function(x){
  data =data.frame(data_plot[,c("pairinfo","group")],gene=data_plot[,x])
  ggplot(data, aes(group,gene,fill=group)) +
    geom_boxplot() +
    geom_point(size=2, alpha=0.5) +
    geom_line(aes(group=pairinfo), colour="black", linetype="11") +
    xlab("") +
    ylab(paste("Expression of ",x))+
    theme_classic()+
    theme(legend.position = "none")
})
library(cowplot)
plot_grid(plotlist=plotlist, ncol=4,labels = LETTERS[1:8])

9.后續(xù)分析
①熱圖:

library(pheatmap)
## 設(shè)定差異基因閾值，減少差異基因用于提取表達(dá)矩陣
allDiff_pair=topTable(fit,adjust='BH',coef="group_listafter",number=Inf,p.value=0.05,lfc =0.5) 
##提前部分?jǐn)?shù)據(jù)用作熱圖繪制
heatdata <- exprSet[rownames(allDiff_pair),]
##制作一個(gè)分組信息用于注釋
annotation_col <- data.frame(group_list)
rownames(annotation_col) <- colnames(heatdata)
pheatmap(heatdata, 
         cluster_rows = TRUE,
         cluster_cols = TRUE,
         annotation_col =annotation_col, 
         annotation_legend=TRUE, 
         show_rownames = F,
         show_colnames = F,
         scale = "row",
         color =colorRampPalette(c("blue", "white","red"))(100))

畫熱圖的意義:
第一看樣本質(zhì)量：本來(lái)before和after兩組應(yīng)該完全分開的侈沪，但是熱圖里面after有兩個(gè)樣本跟bfefore分不開揭璃，要考慮是不是測(cè)量失誤，還是本身樣本就特殊亭罪；
第二看差異基因：差異基因提取出來(lái)的熱圖瘦馍，就應(yīng)當(dāng)呈現(xiàn)橫豎兩條線，把表格分成四個(gè)象限应役，也就是差異基因有高有低情组，這才符合常識(shí)。

②火山圖

library(ggplot2)
library(ggrepel)
library(dplyr)
data <- topTable(fit,adjust='BH',coef="group_listafter",number=Inf) 
data$significant <- as.factor(data$adj.P.Val<0.05 & abs(data$logFC) > 0.5)
data$gene <- rownames(data)
ggplot(data=data, aes(x=logFC, y =-log10(adj.P.Val),color=significant)) +
  geom_point(alpha=0.8, size=1.2,col="black")+
  geom_point(data=subset(data, logFC > 0.5),alpha=0.8, size=1.2,col="red")+
  geom_point(data=subset(data, logFC < -0.5),alpha=0.6, size=1.2,col="blue")+
  labs(x="log2 (fold change)",y="-log10 (adj.P.Val)")+
  theme(plot.title = element_text(hjust = 0.4))+
  geom_hline(yintercept = -log10(0.05),lty=4,lwd=0.6,alpha=0.8)+
  geom_vline(xintercept = c(0.5,-0.5),lty=4,lwd=0.6,alpha=0.8)+
  theme_bw()+
  theme(panel.border = element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),   
        axis.line = element_line(colour = "black")) +
  geom_point(data=subset(data, abs(logFC) > 1),alpha=0.8, size=3,col="green")+
  geom_text_repel(data=subset(data, abs(logFC) > 1), 
                  aes(label=gene),col="black",alpha = 0.8)

③clusterprofiler作圖
GO分析：

suppressMessages(library(clusterProfiler))
#獲得基因列表
gene <- rownames(allDiff)
#基因名稱轉(zhuǎn)換箩祥，返回的是數(shù)據(jù)框
gene = bitr(gene, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")
de <- gene$ENTREZID
## GO分析
go <- enrichGO(gene = de, OrgDb = "org.Hs.eg.db", ont="all")
library(ggplot2)
p <- dotplot(go, split="ONTOLOGY") +facet_grid(ONTOLOGY~., scale="free")
p

KEGG通路富集分析:

EGG <- enrichKEGG(gene= gene$ENTREZID,
                  organism     = 'hsa',
                  pvalueCutoff = 0.05)
dotplot(EGG)

把富集的結(jié)果變成數(shù)據(jù)框院崇，查看凋亡通路hsa04210：

test <- data.frame(EGG)
browseKEGG(EGG, 'hsa04210')