莎莎寫于2020.5.7?今天上午10點鐘準時就醒了，然后和室友磨嘰半天去實驗室打掃衛(wèi)生贷笛，然后簽了字就偷偷溜了??然后回到中校區(qū)買了超級甜的哈密瓜??還有火龍果应又，沃柑，一共是30￥ 關(guān)于今天要學習的內(nèi)容就是：昨天晚上做PPT的時候乏苦，我想知道一些母基因的Regulation情況丁频，但是現(xiàn)成的表格信息有很多NA，我就想著要不自己分析邑贴。然后我就找到一個示例數(shù)據(jù)集
GSE14520，它對應的兩個平臺是GPL571和GPL3921
#數(shù)據(jù)下載
rm(list = ls())
options(stringsAsFactors = F)
# 構(gòu)建結(jié)果儲存文件夾
if(!dir.exists("data")) dir.create("data")
if(!dir.exists("figures")) dir.create("figures")
suppressMessages(library(GEOquery))
gse = "GSE14520"
eSet <- getGEO(gse,
destdir = '.',
getGPL = F)
# class(eSet) #可以看一下是什么類型的數(shù)據(jù)：list
# length(eSet) #看這個list中有多少個個元素：2
# class(eSet[[1]]) #提取列表中的子集
image.png
image.png
可以看到上面有2個elements叔磷，有兩個元素
#(1)提取表達矩陣exp
exp <- exprs(eSet[[1]]) #GPL3921
exp1 <- exprs(eSet[[2]]) #GPL571
exp[1:4,1:4]
#dim(exp) #22268行 445列
#exp = log2(exp+1)
#(2)提取臨床信息
pd <- pData(eSet[[1]])
# class(pd) #是一個data.frame
# dim(pd) #6行 34列
#(3)調(diào)整pd的行名順序與exp列名完全一致：用identical函數(shù)
p = identical(rownames(pd),colnames(exp));p #輸出結(jié)果是TRUE就不需要調(diào)整
if(!p) exp = exp[,match(rownames(pd),colnames(exp))] #用match函數(shù)根據(jù)pd的行名來調(diào)整exp的列名
image.png
#(4)提取芯片平臺編號
gpl <- eSet[[1]]@annotation
gpl1 <- eSet[[2]]@annotation
image.png
可以看到通過提取列表中的子集就可以提取自己想要的平臺對應的數(shù)據(jù)拢驾，相當于兩個平臺當作兩個數(shù)據(jù)集來做就可以了:
exp <- exprs(eSet[[1]]) #GPL3921
exp1 <- exprs(eSet[[2]]) #GPL571
我需要的就是GPL3921,下面開始進行平臺注釋
平臺注釋有很多種方法：
#2.ids-----------------
#方法1 BioconductorR包
gpl
#http://www.bio-info-trainee.com/1399.html
if(!require(hugene10sttranscriptcluster.db))BiocManager::install("hugene10sttranscriptcluster.db")
library(hugene10sttranscriptcluster.db)
ls("package:hugene10sttranscriptcluster.db")
ids <- toTable(hugene10sttranscriptclusterSYMBOL)
head(ids)
# 方法2 讀取gpl頁面的soft文件，按列取子集
# 方法3 官網(wǎng)下載
# 方法4 自主注釋
我一開始嘗試了從GEO網(wǎng)站下載改基，但是最后壓縮文件是不完整的
image.png
然后就去了這個網(wǎng)站：http://www.bio-info-trainee.com/1399.html
正是我需要的
那么接下來就可以按照第一種bioconductor包平臺注釋方法來繼續(xù)往下走了
#2.ids-----------------
#如果鏡像一直出現(xiàn)問題導致下載包打不開繁疤，可以試試下面這個代碼,然后再下載
options(BioC_mirror="http://mirrors.cloud.tencent.com/bioconductor") #這個是Bioconductor的包下載之前設(shè)置的鏡像
options("repos" = c(CRAN="https://mirrors.aliyun.com/CRAN/")) #這個是CRAN上面的包下載之前要設(shè)置的鏡像
#方法1 BioconductorR包
gpl
#http://www.bio-info-trainee.com/1399.html
if(!require(hthgu133a.db))BiocManager::install("hthgu133a.db")
library(hthgu133a.db)
ls("package:hthgu133a.db")
ids <- toTable(hthgu133aSYMBOL)
head(ids)
內(nèi)心想法：這個 包下的可真不容易啊，鏡像也有問題秕狰，然后又設(shè)置了鏡像稠腊，最后才下載好!??
image.png
image.png
SYMBOL是我們需要的注釋信息ids
##下面的代碼就可以提取我們要的ids
ids <- toTable(hthgu133aSYMBOL)
head(ids)
image.png
然后保存這些變量，繼續(xù)往下走：開始畫PCA和熱圖
rm(list = ls())
load(file = "data/step1output.Rdata")
load(file = "data/step2output.Rdata")
#輸入數(shù)據(jù)：exp和group_list
#Principal Component Analysis:PCA詳細原理
#http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/112-pca-principal-component-analysis-essentials
dat=as.data.frame(t(exp)) #先轉(zhuǎn)置鸣哀，然后再as.data.frame
library(FactoMineR)#畫主成分分析圖需要加載這兩個包
library(factoextra)
# pca的統(tǒng)一操作走起
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_list, # color by groups
#palette = c("#00AFBB", "#E7B800"), #這里可以修改顏色架忌，網(wǎng)頁搜索RGB選擇子集喜歡的顏色替換即可
addEllipses = TRUE, # Concentration ellipses
legend.title = "Groups"
)
pca_plot
ggsave(plot = pca_plot,filename = paste0(gse,"PCA.png"))
save(pca_plot,file = "data/pca_plot.Rdata")
PCA圖
#熱圖
cg=names(tail(sort(apply(exp,1,sd)),1000)) #代碼從里向外看，先是apply三個參數(shù)設(shè)置好我衬，對每一行取標準差sd叹放，然后是sort排序，最后是取最大的1000個基因名
n=exp[cg,] #用cg作為行名挠羔，取子集,然后就得到了標準差最大的1000個基因井仰，接下來用于畫圖
#繪制熱圖
annotation_col=data.frame(group=group_list)
rownames(annotation_col)=colnames(n)
library(pheatmap)
pheatmap(n,
show_colnames =F,
show_rownames = F,
annotation_col=annotation_col,
scale = "row")
dev.off()
熱圖：可能是沒有標準化吧，所以看起來奇奇怪怪的
取差異基因：
rm(list = ls())
load(file = "data/step2output.Rdata")
#差異分析破加，用limma包來做
#需要表達矩陣和group_list俱恶，不需要改
suppressMessages(library(limma))
design=model.matrix(~group_list)
fit=lmFit(exp,design)
fit=eBayes(fit)
deg=topTable(fit,coef=2,number = Inf) #這個deg就是差異分析得到的結(jié)果
#class(deg) #是一個data.frame
#dim(deg) #22268 6
image.png
image.png
下面是做一些準備：
#為deg數(shù)據(jù)框添加幾列為后面的可視化和富集分析做準備
#1.加probe_id列，把行名變成一列范舀，與ids匹配做準備
#新增一列有兩種方法：deg$probe_id=rownames(deg)合是，第2種方法是下面的mutate
suppressMessages(library(dplyr))
deg <- mutate(deg,probe_id=rownames(deg))
head(deg) #取了交集之后，肯定會丟失一部分行（基因）尿背，就是沒有注釋的一些基因
#2.加symbol列端仰，火山圖要用
deg <- inner_join(deg,ids,by="probe_id") #probe_id是共有的列名
#上面這一步用inner_join和merge都可以，但是用inner_join的好處是可以保留deg的順序
head(deg)
#按照symbol列去重復：用duplicated根據(jù)邏輯值去重復田藐，注意如果用unique函數(shù)會出現(xiàn)大片NA值
deg <- deg[!duplicated(deg$symbol),]
#3.加change列,標記上下調(diào)基因
# 前兩行是定義閾值：這里的改變會直接影響差異基因的數(shù)量
logFC_t=1
P.Value_t = 0.01
k1 = (deg$P.Value < P.Value_t)&(deg$logFC < -logFC_t)
k2 = (deg$P.Value < P.Value_t)&(deg$logFC > logFC_t)
change = ifelse(k1,
"down",
ifelse(k2,"up","stable"))
#統(tǒng)計一下差異基因的個數(shù)：table(change)
deg <- mutate(deg,change) #把change列加到deg表中
之前是22268行荔烧，現(xiàn)在轉(zhuǎn)換ID加上去重復之后變成12403行
image.png
最后再加上ENTREZID:
#4.加ENTREZID列吱七，用于富集分析（symbol轉(zhuǎn)entrezid，然后inner_join）
suppressMessages(library(ggplot2))
suppressMessages(library(clusterProfiler))
suppressMessages(library(org.Hs.eg.db))
s2e <- bitr(deg$symbol, #bitr是轉(zhuǎn)換的函數(shù)
fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = org.Hs.eg.db)#人類
#class(s2e) #是一個data.frame
#其他物種http://bioconductor.org/packages/release/BiocViews.html#___OrgDb
deg <- inner_join(deg,s2e,by=c("symbol"="SYMBOL"))
image.png
下面就是可視化環(huán)節(jié)了，激不激動!!走到這一步真的很不容易啊~~~??
rm(list = ls())
load(file = "data/step1output.Rdata")
load(file = "data/step4output.Rdata")
#1.火山圖----
suppressMessages(library(dplyr))
suppressMessages(library(ggplot2))
dat = deg
for_label <- dat%>%
filter(abs(logFC)>logFC_t & P.Value < P.Value_t) %>% #這一步是取差異基因
head(10) #取差異最顯著的前10個基因
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
還是挺好看的呀哈哈哈~
#下面是一些進階性的需求：在火山圖上添加感興趣的基因名
volcano_plot <- p +
geom_point(size = 3, shape = 1, data = for_label) +
ggrepel::geom_label_repel(
aes(label = symbol),
data = for_label, #前面一步取的差異基因最顯著的10個基因
color="black"
)
volcano_plot
ggsave(plot = volcano_plot,filename = paste0(gse,"volcano.png"))
添加了差異最顯著的10個基因，也可以添加自己想添加的基因
#2.差異基因熱圖----
load(file = 'data/step2output.Rdata')
x=deg$logFC
names(x)=deg$probe_id
cg=c(names(head(sort(x),30)),names(tail(sort(x),30))) #取30個上調(diào)和30個下調(diào)
#cg = deg$probe_id[deg$change !="stable"] #這句是不想挑30亭敢，想把所有的都畫出來殉农，取不等于stable的所有行
n=exp[cg,]
dim(n)
#作熱圖
library(pheatmap)
annotation_col=data.frame(group=group_list)
rownames(annotation_col)=colnames(n)
library(ggplotify)
heatmap_plot <- as.ggplot(pheatmap(n,show_colnames =F,
show_rownames = F,
scale = "row",
#cluster_cols = F,
annotation_col=annotation_col))
heatmap_plot
ggsave(heatmap_plot,filename = paste0(gse,"heatmap.png"))
image.png
#下面是把兩個圖拼到一起
load("data/pca_plot.Rdata")
library(patchwork)
(pca_plot + volcano_plot +heatmap_plot)+ plot_annotation(tag_levels = "A")
#以3:1的比例來保存
ggsave(filename = "figures/test.png",width = 15,height = 5)
image.png
下面是富集分析：
rm(list = ls())
load(file = 'data/step4output.Rdata')
#富集分析考驗網(wǎng)速，因此給大家保存了Rdata
#上課運行示例數(shù)據(jù)無需修改宛裕，在做自己的數(shù)據(jù)時請注意把本行之后的load()去掉
library(clusterProfiler)
library(dplyr)
library(ggplot2)
source("kegg_plot_function.R")
#source的意思是在不打開某個腳本的情況下運行它，所以該腳本必須在工作目錄下才可以運行
#source表示運行整個kegg_plot_function.R腳本，里面是一個function
#以up_kegg和down_kegg為輸入數(shù)據(jù)做圖
#1.GO database analysis ----
#(1)輸入數(shù)據(jù):富集分析的輸入數(shù)據(jù)是EntrezID
gene_up = deg[deg$change == 'up','ENTREZID']
gene_down = deg[deg$change == 'down','ENTREZID']
gene_diff = c(gene_up,gene_down)
gene_all = deg[,'ENTREZID']
#(2)GO分析窜管，分三部分
#以下步驟耗時很長，實際運行時注意把if后面的括號里F改成T
if(T){
#細胞組分 CC
ego_CC <- enrichGO(gene = gene_diff, #輸入數(shù)據(jù)是差異基因
OrgDb= org.Hs.eg.db, #人類對應的數(shù)據(jù)庫
ont = "CC",
pAdjustMethod = "BH",
minGSSize = 1,
pvalueCutoff = 0.01,
qvalueCutoff = 0.01,
readable = TRUE)
#生物過程 BP
ego_BP <- enrichGO(gene = gene_diff,
OrgDb= org.Hs.eg.db,
ont = "BP",
pAdjustMethod = "BH",
minGSSize = 1,
pvalueCutoff = 0.01,
qvalueCutoff = 0.01,
readable = TRUE)
#分子功能： MM
ego_MF <- enrichGO(gene = gene_diff,
OrgDb= org.Hs.eg.db,
ont = "MF",
pAdjustMethod = "BH",
minGSSize = 1,
pvalueCutoff = 0.01,
qvalueCutoff = 0.01,
readable = TRUE)
save(ego_CC,ego_BP,ego_MF,file = "data/ego_GSE14520.Rdata")
}
load(file = "data/ego_GSE14520.Rdata")
#class(ego_BP) #數(shù)據(jù)類型是enrichResult,出自DOSE包
# tmp = ego_BP@result #看PPT中的表
# View(tmp)
然后進行可視化：
#(3)可視化
#條帶圖
barplot(ego_CC,showCategory=20)
#氣泡圖
dotplot(ego_CC) #不設(shè)置showCategory=幾的話就默認展示10個
條帶圖（barplot）
氣泡圖（dotplot）
#下面的圖需要映射顏色稚机，設(shè)置和示例數(shù)據(jù)一樣的geneList
# library(DOSE)
# data(geneList) #這個是DOSE包中的示例數(shù)據(jù)幕帆，不用運行
geneList = deg$logFC
names(geneList)=deg$ENTREZID
geneList = sort(geneList,decreasing = T) #排序是從大到小
#(3)展示top5通路的共同基因，要放大看赖条。
#Gene-Concept Network
cnetplot(ego_CC, categorySize="pvalue", foldChange=geneList,colorEdge = TRUE)
cnetplot(ego_CC, foldChange=geneList, circular = TRUE, colorEdge = TRUE)
#Enrichment Map
emapplot(ego_CC)
#(4)展示通路關(guān)系
goplot(ego_CC)
cnetplot
image.png
emapplot
展示通路關(guān)系
#(5)Heatmap-like functional classification
heatplot(ego_CC,foldChange = geneList)
#太多基因就會糊失乾。可通過調(diào)整比例或者減少基因來控制纬乍。
pdf("heatplot.pdf",width = 14,height = 5)
heatplot(ego_CC,foldChange = geneList)
dev.off()
因為基因太多了所以底下會糊碱茁，要適當調(diào)整比例
下面是KEGG:
# 2.KEGG pathway analysis----
#上調(diào)、下調(diào)仿贬、差異纽竣、所有基因
#（1）輸入數(shù)據(jù)
gene_up = deg[deg$change == 'up','ENTREZID']
gene_down = deg[deg$change == 'down','ENTREZID']
gene_diff = c(gene_up,gene_down)
gene_all = deg[,'ENTREZID']
#（2）對上調(diào)/下調(diào)/所有差異基因進行富集分析
#注意這里又有個F
if(T){
kk.up <- enrichKEGG(gene = gene_up,
organism = 'hsa',
universe = gene_all,
pvalueCutoff = 0.9,
qvalueCutoff = 0.9)
kk.down <- enrichKEGG(gene = gene_down,
organism = 'hsa',
universe = gene_all,
pvalueCutoff = 0.9,
qvalueCutoff =0.9)
kk.diff <- enrichKEGG(gene = gene_diff,
organism = 'hsa',
pvalueCutoff = 0.9)
save(kk.diff,kk.down,kk.up,file = "data/GSE14520kegg.Rdata")
}
load("data/GSE14520kegg.Rdata")
#(3)從富集結(jié)果中提取出結(jié)果數(shù)據(jù)框
kegg_diff_dt <- kk.diff@result
#(4)按照pvalue篩選通路
#在enrichkegg時沒有設(shè)置pvaluecutoff，在此處篩選
down_kegg <- kk.down@result %>%
filter(pvalue<0.05) %>% #篩選行
mutate(group=-1) #新增列
up_kegg <- kk.up@result %>%
filter(pvalue<0.05) %>%
mutate(group=1)
#(5)可視化
g_kegg <- kegg_plot(up_kegg,down_kegg)
#g_kegg +scale_y_continuous(labels = c(20,15,10,5,0,5)) #手動修改坐標軸的刻度
ggsave(g_kegg,filename = 'figures/kegg_up_down.png')
image.png
#gsea作kegg富集分析茧泪，可選----
#(1)查看示例數(shù)據(jù)
data(geneList, package="DOSE")
#(2)將我們的數(shù)據(jù)轉(zhuǎn)換成示例數(shù)據(jù)的格式
geneList=deg$logFC
names(geneList)=deg$ENTREZID
geneList=sort(geneList,decreasing = T)
#(3)富集分析
kk_gse <- gseKEGG(geneList = geneList,
organism = 'hsa',
nPerm = 1000,
minGSSize = 120,
pvalueCutoff = 0.9,
verbose = FALSE)
down_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore < 0,];down_kegg$group=-1
up_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore > 0,];up_kegg$group=1
#(4)可視化
kegg_plot(up_kegg,down_kegg)
ggsave('figures/kegg_up_down_gsea.png')
kegg_up_down_gsea
基本上一篇完整的生信分析退个，從下載數(shù)據(jù)，處理數(shù)據(jù)调炬，差異基因的分析语盈，以及到后面的可視化包括PCA/火山圖/熱圖，還有最后的富集分析的圖缰泡，就算是一整套流程走下來了刀荒。真的很不容易啊棘钞！??
想想真正開始自己寫簡書也就這段時間缠借，在賓館隔離開始，然后5.1號開始告訴自己要每天寫宜猜，然后現(xiàn)在回到宿舍也一直努力堅持著泼返，希望自己可以一直堅持，寫一些自己喜歡的東西姨拥，心情啊绅喉，想法啊渠鸽，學習筆記什么的。
然后真正學習生信是放寒假開始柴罐，1月10號凌晨的飛機到新橋機場徽缚，然后就開始了我的寒假，因為自己課題需要就好好學習生信革屠，學習R語言凿试，學習各種數(shù)據(jù)庫的使用，磕磕絆絆到現(xiàn)在似芝，4個月那婉，從一開始的代碼是什么意思，頻頻報錯党瓮，不知道怎么辦吧恃，瘋狂問別人，別人也不愛搭理麻诀，然后自己畫圖畫到半夜，到現(xiàn)在可以從頭到尾順下來整個流程傲醉，雖然代碼不是自己寫的蝇闭，但是代碼什么意思應該是可以說出一二的，我還記得之前我問一個朋友怎么學的生信硬毕，他就說總會有陣痛期的呻引，慢慢熬過去就好了，現(xiàn)在想來真的是這樣吐咳。很多東西逻悠，只有自己經(jīng)歷了，去認真學習了韭脊，體會過童谒，焦慮過，印象才最深刻沪羔。
剛開始逛CSDN帖子饥伊，看別人復現(xiàn)一篇又一篇文章，就覺得哇真的好厲害啊蔫饰，我要是也可以這樣就好了琅豆，這是2個月前的想法，由于自己一直不斷努力篓吁，到現(xiàn)在茫因，我也可以完整復現(xiàn)一篇文章，玩這個走完一套流程杖剪，自己學會了搜索冻押，學會解決報錯驰贷，最后出圖，然后在這里記錄我努力的過程翼雀，畫圖的結(jié)果饱苟，這樣的過程是快樂的，欣慰的狼渊。
我會一直努力下去的?加油??豬莎??
2020.6.22 后面要更一波通過芯片數(shù)據(jù)集提取生存分析的數(shù)據(jù)從而做生存分析箱熬，先在這里立一個Flag~
