limma分析差异基因
数据载入
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
## [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)
image
## 可视化结果
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
)
image
层次聚类图-聚类结果也与PCA相似，结果并不好
聚类分析的结果也同样可以进一步美化，但这里不做
计算距离时同样需进行转置，但在前一步PCA分析中的data已经经过转置，故未重复
dd <- dist(data, method = "euclidean")##data是经过行列转换的
hc <- hclust(dd, method = "ward.D2")
plot(hc)
image
##对结果进行美化
# 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")
image
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))
image
高级的画法
借助于有人开发的更高级的包，用于完成某些特殊的功能，或者更美观
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)
image
绘制热图
热图的使用比较频繁，得到差异基因后可以直接绘制热图
相对简单好用的要属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)
image
稍微调整细节
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)
image
R语言GEO数据挖掘03-limma分析差异基因
