PCA和差異基因圖是生信技能樹[生信爆款入門課程]GEO數(shù)據(jù)挖掘的重點。為拓展課堂所學知識，現(xiàn)在找一個數(shù)據(jù)集對他們做下練習總結(jié)。
1.主成分PCA 圖----
> dat=as.data.frame(t(exp))
> library(FactoMineR)
> library(factoextra)
> dat.pca <- PCA(dat, graph = FALSE)
> pca_plot <- fviz_pca_ind(dat.pca,
+ geom.ind = "point", # show points only (nbut not "text")
+ col.ind = Group, # color by groups
+ palette = c("#00AFBB", "#E7B800"),
+ addEllipses = TRUE, # Concentration ellipses
+ legend.title = "Groups"
+ )
> pca_plot
> ggsave(plot = pca_plot,filename = paste0(gse_number,"_PCA.png"))
Saving 6.4 x 3.77 in image
> save(pca_plot,file = "pca_plot.Rdata")
image.png
2.top 1000 sd 熱圖----
> cg=names(tail(sort(apply(exp,1,sd)),1000))
> n=exp[cg,]
> library(pheatmap)
Warning message:
程輯包‘pheatmap’是用R版本4.0.3 來建造的
> annotation_col=data.frame(group=Group)
> rownames(annotation_col)=colnames(n)
> ## 2.自行標準化再畫熱圖
> n2 = t(scale(t(n)))
> pheatmap(n2,
+ show_colnames =F,
+ show_rownames = F,
+ cluster_cols = F,
+ annotation_col=annotation_col,
+ breaks = seq(-3,3,length.out = 100)
+ )
image.png
3.相關性熱圖
> pheatmap::pheatmap(cor(exp),
+ annotation_col = annotation_col)
>
image.png
- 差異基因火山圖
1.火山圖----
> library(dplyr)
> library(ggplot2)
> dat = deg
> p <- ggplot(data = dat,
+ aes(x = logFC,
+ y = -log10(P.Value))) +
+ geom_point(alpha=0.4, size=3.5,
+ aes(color=change)) +
+ ylab("-log10(Pvalue)")+
+ scale_color_manual(values=c("blue", "grey","red"))+
+ geom_vline(xintercept=c(-logFC_t,logFC_t),lty=4,col="black",lwd=0.8) +
+ geom_hline(yintercept = -log10(P.Value_t),lty=4,col="black",lwd=0.8) +
+ theme_bw()
> p
> load(file = 'step2output.Rdata')
> if(T){
+ #全部差異基因
+ cg = deg$probe_id[deg$change !="stable"]
+ length(cg)
+ }else{
+ #取前30上調(diào)和前30下調(diào)
+ x=deg$logFC[deg$change !="stable"]
+ names(x)=deg$probe_id[deg$change !="stable"]
+ cg=names(c(head(sort(x),30),tail(sort(x),30)))
+ length(cg)
+ }
[1] 1616
> n=exp[cg,]
> dim(n)
[1] 1616 22
image.png
差異基因熱圖
> library(pheatmap)
> annotation_col=data.frame(group=Group)
> rownames(annotation_col)=colnames(n)
> heatmap_plot <- pheatmap(n,show_colnames =F,
+ show_rownames = F,
+ scale = "row",
+ cluster_cols = F,
+ annotation_col=annotation_col,
+ breaks = seq(-3,3,length.out = 100)
+ )
> heatmap_plot
> ggsave(heatmap_plot,filename = paste0(gse_number,"_heatmap.png"))
Saving 6.4 x 3.77 in image
> load("pca_plot.Rdata")
image.png
3.拼圖
> library(patchwork)
Warning message:
程輯包‘patchwork’是用R版本4.0.3 來建造的
> library(ggplotify)
Warning message:
程輯包‘ggplotify’是用R版本4.0.3 來建造的
> (pca_plot + p +as.ggplot(heatmap_plot))
>
image.png
