2021-07-16

2021-07-16 本文已影响0人学习生信的小兔子

## 加载R包
library(GEOquery)
## 下载数据，如果文件夹中有会直接读入
gset = getGEO('GSE32575', destdir=".",getGPL = F)
## 获取ExpressionSet对象，包括的表达矩阵和分组信息
gset=gset[[1]]
## 获取分组信息,点开查阅信息
pdata=pData(gset)
## 只要后36个,本次选择的这36个是配对的。
## 所以，别人的芯片我们也不是全要，我们只要适合自己的数据
group_list=c(rep('before',18),rep('after',18))
group_list=factor(group_list)
## 强制限定顺序
group_list <- relevel(group_list, ref="before")
exprSet=exprs(gset)
boxplot(exprSet[,-c(1:12)],outline=FALSE, notch=T,col=group_list, las=2)
library(limma) 
exprSet=normalizeBetweenArrays(exprSet)
boxplot(exprSet,outline=FALSE, notch=T,col=group_list, las=2)

未校正前

library(limma) 
exprSet=normalizeBetweenArrays(exprSet)
boxplot(exprSet,outline=FALSE, notch=T,col=group_list, las=2)

校正后

我们通过分组信息已经知道，前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")}
[1] "log2 transform finished"

获取注释信息

platformMap <- data.table::fread("platformMap.txt")
##可知注释包为illuminaHumanv2.db
library(illuminaHumanv2.db)

获取探针对应的symbol信息

## 获取表达矩阵的行名，就是探针名称
probeset <- rownames(exprSet)
## 使用lookup函数，找到探针在illuminaHumanv2.db中的对应基因名称
## 如果分析别的芯片数据，把illuminaHumanv2.db更换即可
SYMBOL <-  annotate::lookUp(probeset,"illuminaHumanv2.db", "SYMBOL")
## 转换为向量
symbol = as.vector(unlist(SYMBOL))
probe2symbol = data.frame(probeset,symbol,stringsAsFactors = F)

一个基因会对应对个探针，有些基因名称会是重复的，这些都需要处理。
对于，多个探针，我们选取在样本中平均表达量最高的探针作为对应基因的表达量。一下代码完成所有事情，而且可以复用。

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")

现在数据变成这个样子

探针去重与转换

而且，行数从原来的48701变成18948行。

差异分析

如果没有配对信息，差异分析这样做：

##差异分析没有配对
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) 
dim(allDiff)
[1] 6798    6

其中allDiff得到6798行。
因为我们这里实际上是有配对信息的，差异分析应该这样做：

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)
[1] 7648    6

得到的差异分析结果allDiff_pair有7648个，比直接做差异分析要多接近1000个。

这里配对和非配对的差异分析，只相差一步，也就是是否有配对信息

作图验证

首先基因名称需要在列，所以需要用t函数，行列转置。

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

data_plot局部

作图展示:
挑选了一个基因CAMKK2

library(ggplot2)
ggplot(data_plot, aes(group,CAMKK2,fill=group)) +
  geom_jitter(aes(colour=group), size=2, alpha=0.5)+
  xlab("") +
  ylab(paste("Expression of ","CAMKK2"))+
  theme_classic()+
  theme(legend.position = "none")

无配对信息

加入配对信息

library(ggplot2)
ggplot(data_plot, aes(group,CAMKK2,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 ","CAMKK2"))+
  theme_classic()+
  theme(legend.position = "none")

加入配对信息

再来看看真正有差异的基因

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")

CH25H

批量画出差异最大的8个基因：

![差异最大的前8个基因](https://img.haomeiwen.com/i25205178/416bf599ce4c5d00.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)

热图

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)

火山图

火山图画起来简单，难的是如何定制化展示。比如在图上有不同的颜色，不同的点来表示数据。有什么意义呢？我其实理解的也不是很透彻，但是考虑到他的横坐标是变化倍数，纵坐标是p值的负数，那么p值越小，倍数越大的基因就会出现在左右上方。
数据是死的，而图形化展示是活的，火山图可能是大家觉得比较好的展示差异基因的方法。从这个图中我们还可以看出，两边是否对称，比如有的药物就是抑制基因转录的，那么加了药物的两组，可能会出现，上调和下调基因不对称的情形，这个从数据中可以分析出来，但是用图就一目了然了。
此外还有耳熟能详的GO分析，KEGG分析，我觉得都是名气很大用的很少的聚类方法。
做这个的意义在于，我们获得了几千个差异基因，但是我们确定最重要的基因呢？因人而异，不同的课题组，想的不一样。但是大体的思路是聚类，假如P53通路有100个基因，对于人类20000个基因而言，随便抓一把基因，出现P53通路基因的概率是100/20000,是0.005，现在3000个差异基因中就出现了50个P53通路的基因，这个概率是0.016，明显超过了他随机出现的概率，我们就可以说，p53通路被富集了，这个通路可能在你研究的课题中有重要作用。
Y叔的神包clusterprofiler可以做出特别好看的图。
比如GO分析，有三个组成部分，分别是CC，cellular compartment 细胞组分，BP，biological process，生物进程，MF，molecular function，分子功能。

GO分析

## 加载R包
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
....电脑卡住了 hhhhh 卑微.jpg

GO分析

如何看这个图，如何解释？看图容易解释难。简单说来，一看颜色，二看大小。颜色越红p值越小越可信，点越打表示富集的基因比例越多，越有可能被验证。

举个例子，BP那一项中，有个叫neutrophil mediate immunity的GO通路被富集了，看字面意思应该是中性粒细胞介导的免疫反应，联系到这到实验设计，就是一开始的summary和文章部分，我觉得这个是合理的，因为他研究的就是长期慢性炎症对肥胖的影响，这里我们发现，炎症导致的免疫反应可能可以解释这个事情。还有很多其他的GO通路，需要自己根据背景去解释。去发现线索。

KEGG分析

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

假如我们对凋亡感兴趣，我们还可以看看我们有多少基因被富集在了凋亡通路，Y叔的神包clusterprofiler也是可以做的。而且做的很云淡风轻。
首先我们把富集的结果变成数据框，便于自己查看。

test <- data.frame(EGG)

我们发现想看的凋亡通路牌号是，hsa04210，那我们就用下面的代码浏览一下，会跳出一个网页。

browseKEGG(EGG, 'hsa04210')

2021-07-16

获取注释信息

差异分析

作图验证

热图

火山图

GO分析

KEGG分析

猜你喜欢

热点阅读