limma分析差异基因

在经过了前两期中的数据下载，数据基本处理之后，解决了一个探针对应多个基因数的
以及多个探针对应一个基因求平均值，在此基础上运用limma包分析差异基因
除此以外，包括绘制火山图，热图，PCA等，都在本文中解决

数据载入

if(T){
Sys.setlocale('LC_ALL','C')
library(dplyr)
##
if(T){
load("expma.Rdata")
load("probe.Rdata")
}
expma[1:5,1:5]
boxplot(expma)##看下表达情况
metdata[1:5,1:5]
head(probe)

## 查看Gene Symbol是否有重复
table(duplicated(probe$`Gene Symbol`))##12549 FALSE

## 整合注释信息到表达矩阵
ID<-rownames(expma)
expma<-apply(expma,1,function(x){log2(x+1)})
expma<-t(expma)
eset<-expma[ID %in% probe$ID,] %>% cbind(probe)
eset[1:5,1:5]
colnames(eset)
}

image.png

## [1] "GSM188013"      "GSM188014"      "GSM188016"      "GSM188018"     
## [5] "GSM188020"      "GSM188022"      "ID"             "Gene Symbol"   
## [9] "ENTREZ_GENE_ID"

多个探针求平均值

test1<-aggregate(x=eset[,1:6],by=list(eset$`Gene Symbol`),FUN=mean,na.rm=T) 
test1[1:5,1:5]##与去重结果相吻合
##   Group.1 GSM188013 GSM188014 GSM188016 GSM188018
## 1          8.438846  8.368513  7.322442  7.813573
## 2    A1CF 10.979025 10.616926  9.940773 10.413311
## 3     A2M  6.565276  6.422112  8.142194  5.652593
## 4  A4GALT  7.728628  7.818966  8.679885  7.048563
## 5   A4GNT 10.243388 10.182382  9.391991  8.779887
dim(test1)##
## [1] 12549     7
colnames(test1)
## [1] "Group.1"   "GSM188013" "GSM188014" "GSM188016" "GSM188018" "GSM188020"
## [7] "GSM188022"
rownames(test1)<-test1$Group.1
test1<-test1[,-1]
eset_dat<-test1

PCA plot

Principal Component Analysis (PCA)分析使用的是基于R语言的 prcomp() and princomp()函数.
完成PCA分析一般有两种方法： princomp()使用的是一种的spectral decomposition方法
prcomp() and PCA()[FactoMineR]使用的是SVD法

data<-eset_dat
data[1:5,1:5]##表达矩阵
##        GSM188013 GSM188014 GSM188016 GSM188018 GSM188020
##         8.438846  8.368513  7.322442  7.813573  7.244615
## A1CF   10.979025 10.616926  9.940773 10.413311  9.743305
## A2M     6.565276  6.422112  8.142194  5.652593  5.550033
## A4GALT  7.728628  7.818966  8.679885  7.048563  5.929258
## A4GNT  10.243388 10.182382  9.391991  8.779887  9.431585
metdata[1:5,1:5]
##                                                              title
## GSM188013   DMSO treated MCF7 breast cancer cells [HG-U133A] Exp 1
## GSM188014 Dioxin treated MCF7 breast cancer cells [HG-U133A] Exp 1
## GSM188016   DMSO treated MCF7 breast cancer cells [HG-U133A] Exp 2
## GSM188018 Dioxin treated MCF7 breast cancer cells [HG-U133A] Exp 2
## GSM188020   DMSO treated MCF7 breast cancer cells [HG-U133A] Exp 3
##           geo_accession                status submission_date
## GSM188013     GSM188013 Public on May 12 2007     May 08 2007
## GSM188014     GSM188014 Public on May 12 2007     May 08 2007
## GSM188016     GSM188016 Public on May 12 2007     May 08 2007
## GSM188018     GSM188018 Public on May 12 2007     May 08 2007
## GSM188020     GSM188020 Public on May 12 2007     May 08 2007
##           last_update_date
## GSM188013      May 11 2007
## GSM188014      May 11 2007
## GSM188016      May 11 2007
## GSM188018      May 11 2007
## GSM188020      May 11 2007
##构建group_list
group_list<-rep(c("Treat","Control"),3)
colnames(data)<-group_list
library(factoextra)
## Warning: package 'factoextra' was built under R version 3.5.3
## Loading required package: ggplot2
## Welcome! Related Books: `Practical Guide To Cluster Analysis in R` at https://goo.gl/13EFCZ
## 计算PCA
data<-t(data)##转换数据至行为sample,列为gene
data<-as.data.frame(data)##注意数据要转换为数据框
res.pca <- prcomp(data, scale = TRUE)
##展示主成分对差异的贡献
fviz_eig(res.pca)

Fig2

## 可视化结果
fviz_pca_ind(res.pca,
             col.ind = group_list, # 颜色对应group信息
             palette = c("#00AFBB",  "#FC4E07"),
             addEllipses = TRUE, # Concentration ellipses
             ellipse.type = "confidence",
             legend.title = "Group",## Legend名称
             repel = TRUE
             )

Fig3

层次聚类图-聚类结果也与PCA相似，结果并不好

聚类分析的结果也同样可以进一步美化，但这里不做
计算距离时同样需进行转置，但在前一步PCA分析中的data已经经过转置，故未重复

dd <- dist(data, method = "euclidean")##data是经过行列转换的
hc <- hclust(dd, method = "ward.D2")
plot(hc)

Fig4

##对结果进行美化
# Convert hclust into a dendrogram and plot
hcd <- as.dendrogram(hc)
# Define nodePar
nodePar <- list(lab.cex = 0.6, pch = c(NA, 19), 
                cex = 0.7, col = "blue")
# Customized plot; remove labels
plot(hcd, ylab = "Height", nodePar = nodePar, leaflab = "none")

Fig5

limma分析差异基因

limma包的j具具具体用法参考 limma Users Guide
构建分组信息，构建好比较矩阵是关键
注意这里的表达矩阵信息 eset_dat是经过处理后的，为转置，行为gene，列为sample

library(limma)
library(dplyr)
group_list
## [1] "Treat"   "Control" "Treat"   "Control" "Treat"   "Control"
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
design
##   Control Treat
## 1       0     1
## 2       1     0
## 3       0     1
## 4       1     0
## 5       0     1
## 6       1     0
## attr(,"assign")
## [1] 1 1
## attr(,"contrasts")
## attr(,"contrasts")$`factor(group_list)`
## [1] "contr.treatment"
## 比较信息
contrast.matrix<-makeContrasts("Treat-Control",
                               levels = design)
contrast.matrix##查看比较矩阵的信息，这里我们设置的是Treat Vs. control
##          Contrasts
## Levels    Treat-Control
##   Control            -1
##   Treat               1
## 拟合模型
fit <- lmFit(eset_dat,design)
fit2 <- contrasts.fit(fit, contrast.matrix) 
fit2 <- eBayes(fit2) 
DEG<-topTable(fit2, coef=1, n=Inf) %>% na.omit()  ## coef比较分组 n基因数
head(DEG)
##             logFC   AveExpr          t      P.Value adj.P.Val         B
## ALDH3A1 -3.227263 10.302323 -10.710306 4.482850e-05 0.3134585 -4.048355
## CYP1B1  -3.033684 13.287607 -10.505888 4.995753e-05 0.3134585 -4.049713
## CYP1A1  -9.003353 11.481268  -8.371476 1.762905e-04 0.7374232 -4.069681
## HHLA2   -1.550587  6.595658  -7.443431 3.337672e-04 0.9308066 -4.083411
## SLC7A5  -2.470333 13.628775  -7.298868 3.708688e-04 0.9308066 -4.085966
## TIPARP  -1.581274 12.764218  -7.024252 4.552834e-04 0.9522252 -4.091192
dim(DEG)
## [1] 12549     6
save(DEG,file = "DEG_all.Rdata")

绘制火山图

火山图其实仅仅是一种可视化的方式，能够从整体上让我们对整体的差异分析情况有个了解
筛选到差异基因后，可以直接绘制出火山图
火山图的横坐标为logFC, 纵坐标为-log10(pvalue)，因此其实理论上讲plot即可完成火山图绘制

最简单原始的画法

colnames(DEG)
## [1] "logFC"     "AveExpr"   "t"         "P.Value"   "adj.P.Val" "B"
plot(DEG$logFC,-log10(DEG$P.Value))

Fig6

高级的画法

借助于有人开发的更高级的包，用于完成某些特殊的功能，或者更美观

require(EnhancedVolcano)
EnhancedVolcano(DEG,
                
                lab = rownames(DEG),
                
                x = "logFC",
                
                y = "P.Value",
                
                selectLab = rownames(DEG)[1:5],
                
                xlab = bquote(~Log[2]~ "fold change"),
                
                ylab = bquote(~-Log[10]~italic(P)),
                
                pCutoff = 0.05,## pvalue阈值
                
                FCcutoff = 1,## FC cutoff
                
                xlim = c(-5,5),
                
                transcriptPointSize = 1.8,
                
                transcriptLabSize = 5.0,
                
                colAlpha = 1,
                
                legend=c("NS","Log2 FC"," p-value",
                         " p-value & Log2 FC"),
                
                legendPosition = "bottom",
                
                legendLabSize = 10,
                
                legendIconSize = 3.0)

Fig7

绘制热图

热图的使用比较频繁，得到差异基因后可以直接绘制热图
相对简单好用的要属pheatmap包了
管道中的常规提取需要加上特殊的占位符.

## 首先提取出想要画的数据
head(DEG)
## 提取FC前50
up_50<-DEG %>% as_tibble() %>% 
  mutate(genename=rownames(DEG)) %>% 
  dplyr::arrange(desc(logFC)) %>% 
  .$genename %>% .[1:50] ## 管道符中的提取
## FC低前50
down_50<-DEG %>% as_tibble() %>% 
  mutate(genename=rownames(DEG)) %>% 
  dplyr::arrange(logFC) %>% 
  .$genename %>% .[1:50] ## 管道符中的提取
index<-c(up_50,down_50)  

## 开始绘图-最简单的图
library(pheatmap)
pheatmap(eset_dat[index,],show_colnames =F,show_rownames = F)

Fig8

稍微调整细节

index_matrix<-t(scale(t(eset_dat[index,])))##归一化
index_matrix[index_matrix>1]=1
index_matrix[index_matrix<-1]=-1
head(index_matrix)

## 添加注释
anno=data.frame(group=group_list)
rownames(anno)=colnames(index_matrix)
anno##注释信息的数据框
pheatmap(index_matrix,
           show_colnames =F,
           show_rownames = F,
           cluster_cols = T, 
           annotation_col=anno)
   

Fig9

本期内容就到这里，我是白介素2，下期再见

广而告之

说一个事，鉴于简书平台在信息传播方面有不足之处，应粉丝要求，白介素2的个人微信平台已经开启，继续聊临床与科研的故事，R语言，数据挖掘，文献阅读等内容。当然也不要期望过高，微信平台目前的定位是作为自己的读书笔记，如果对大家有帮助最好。如果感兴趣， 可以扫码关注下。

qrcode_for_gh_9eaa04438675_258.jpg

如果没有时间精力学习代码，推荐了解：零代码数据挖掘课程

转载请注明出处