GEO下载与分析

转自果子学生信
1.加载R包获取数据

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

2.通过pData函数获取分组信息

pdata=pData(gset)
group_list=c(rep('before',18),rep('after',18))
group_list=factor(group_list)
## 强制限定顺序
group_list <- relevel(group_list, ref="before")

3.通过exprs函数获取表达矩阵并校正

exprSet=exprs(gset)
##查看整体样本的表达情况
boxplot(exprSet[,-c(1:12)],outline=FALSE, notch=T,col=group_list, las=2)
##整体表达不整齐，使用limma包内置函数人工校正
library(limma) 
exprSet=normalizeBetweenArrays(exprSet)
boxplot(exprSet,outline=FALSE, notch=T,col=group_list, las=2)

4.判断是否需要进行数据转换

##根据分组信息，去除前12个样本
exprSet = as.data.frame(exprSet)[,-seq(1,12)]
##表达量很大，log转换(选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.注释基因

##导入R包和平台的注释信息对应关系表格 platformMap
platformMap <- data.table::fread("platformMap.txt")
##获取注释平台
index = gset@annotation
##使用代码自动获取对应注释包
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)
##获取探针对应的symbol信息
probeset <- rownames(exprSet)
## 使用lookup函数，找到探针在illuminaHumanv2.db中的对应基因名称
SYMBOL <-  annotate::lookUp(probeset,"illuminaHumanv2.db", "SYMBOL")
## 转换为向量
symbol = as.vector(unlist(SYMBOL))
##制作probe2symbol转换文件
probe2symbol = data.frame(probeset,symbol,stringsAsFactors = F)

6.探针转换与基因去重

library(dplyr)
library(tibble)
exprSet <- exprSet %>% 
  rownames_to_column(var="probeset") %>% 
  #合并探针的信息
  inner_join(probe2symbol,by="probeset") %>% 
  #去掉多余信息
  select(-probeset) %>% 
  #重新排列
  select(symbol,everything()) %>% 
  #求出平均数(这边的点号代表上一步产出的数据)
  mutate(rowMean =rowMeans(.[grep("GSM", names(.))])) %>% 
  #去除symbol中的NA
  filter(symbol != "NA") %>% 
  #把表达量的平均值按从大到小排序
  arrange(desc(rowMean)) %>% 
  # symbol留下第一个
  distinct(symbol,.keep_all = T) %>% 
  #反向选择去除rowMean这一列
  select(-rowMean) %>% 
  # 列名变成行名
  column_to_rownames(var = "symbol")
现在数据变成这个样子

7.差异分析

##如果没有配对信息
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) 
##如果有配对信息
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)

分析结果的各列数据含义：
　　“logFC”是两组表达值之间以2为底对数化的的变化倍数，一般表达相差2倍以上有意义；
　　“AveExpr”是该探针组所在所有样品中的平均表达值；
　　“t”是贝叶斯调整后T 检验的 t 值；
　　“P.Value”是贝叶斯检验的 P 值；
　　“adj.P.Val”是调整后的 P 值，更有参考价值；
　　“B”是经验贝叶斯得到的标准差的对数化值。

8.作图验证(非必要)
转换为ggplot2喜欢的数据格式，行是观测，列是变量，即清洁数据

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

以CAMKK2为例做配对图

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个基因

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.后续分析
①热图:

library(pheatmap)
## 设定差异基因阈值，减少差异基因用于提取表达矩阵
allDiff_pair=topTable(fit,adjust='BH',coef="group_listafter",number=Inf,p.value=0.05,lfc =0.5) 
##提前部分数据用作热图绘制
heatdata <- exprSet[rownames(allDiff_pair),]
##制作一个分组信息用于注释
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))

画热图的意义:
第一看样本质量：本来before和after两组应该完全分开的，但是热图里面after有两个样本跟bfefore分不开，要考虑是不是测量失误，还是本身样本就特殊；
第二看差异基因：差异基因提取出来的热图，就应当呈现横竖两条线，把表格分成四个象限，也就是差异基因有高有低，这才符合常识。

②火山图

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)
#基因名称转换，返回的是数据框
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)

把富集的结果变成数据框，查看凋亡通路hsa04210：

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