一、背景了解 芯片數(shù)據(jù)：首選limma 。
二代測序（RNA_seq）：如果是counts 可選擇limma的voom算法或者edgeR或者DESeq2。 如果是FPKM或TPM可選擇limma，注意： 注意：
if (!require ("BiocManager" )) install.packages("BiocManager" ,update = F ,ask = F )"https://mirrors.ustc.edu.cn/bioc/" )'tidyr' ,'tibble' ,'dplyr' ,'stringr' ,'ggplot2' ,'ggpubr' ,'factoextra' ,'FactoMineR' ,'devtools' ,'cowplot' ,'patchwork' ,'basetheme' ,'paletteer' ,'AnnoProbe' ,'ggthemes' ,'VennDiagram' ,'tinyarray' ) 'GEOquery' ,'hgu133plus2.db' ,'ggnewscale' ,"limma" ,"impute" ,"GSEABase" ,"GSVA" ,"clusterProfiler" ,"org.Hs.eg.db" ,"preprocessCore" ,"enrichplot" ,"ggplotify" )for (pkg in cran_packages){if (! require (pkg,character.only=T ) ) {F ,update = F )require (pkg,character.only=T ) for (pkg in Biocductor_packages){if (! require (pkg,character.only=T ) ) {F ,update = F )require (pkg,character.only=T ) #前面的所有提示和報(bào)錯(cuò)都先不要管。主要看這里 for (pkg in c(Biocductor_packages,cran_packages)){require (pkg,character.only=T ) #沒有任何提示就是成功了，如果有warningxx包不存在，用library檢查一下。 #library報(bào)錯(cuò)，就單獨(dú)安裝。 #數(shù)據(jù)下載 library (GEOquery)#先去網(wǎng)頁確定是否是表達(dá)芯片數(shù)據(jù)，不是的話不能用本流程。 "GSE56649" '.' , getGPL = F )1 ]]#(1)提取表達(dá)矩陣exp 1 :4 ,1 :4 ]#檢查矩陣是否正常，如果是空的就會(huì)報(bào)錯(cuò)，空的和有負(fù)值的、有異常值的矩陣需要處理原始數(shù)據(jù)。 #如果表達(dá)矩陣為空，大多數(shù)是轉(zhuǎn)錄組數(shù)據(jù)，不能用這個(gè)流程(后面另講)。 #自行判斷是否需要log 1 )# 將原始數(shù)據(jù)標(biāo)準(zhǔn)化（針對畫出的圖水平不一致） #(2)提取臨床信息 #(3)讓exp列名與pd的行名順序完全一致 if (!p) exp = exp[,match(rownames(pd),colnames(exp))]#(4)提取芯片平臺編號 "step1output.Rdata" )Q：如何進(jìn)行分組？
A：三種方法：芯片中最常用的是str_detect（）函數(shù)；轉(zhuǎn)錄組數(shù)據(jù)中最常用的是Group = c(rep(“RA”,times=13),rep(“control”,times=9))注意 ：需要把Group轉(zhuǎn)換成因子，并設(shè)置參考水平，指定levels，對照組在前，處理組在后
# Group(實(shí)驗(yàn)分組)和ids(探針注釋) "step1output.Rdata" )library (stringr)# 標(biāo)準(zhǔn)流程代碼是二分組，多分組數(shù)據(jù)的分析后面另講 # 生成Group向量的三種常規(guī)方法，三選一，選誰就把第幾個(gè)邏輯值寫成T，另外兩個(gè)為F。如果三種辦法都不適用，可以繼續(xù)往后寫else if if (F ){# 1.Group---- # 第一種方法，有現(xiàn)成的可以用來分組的列 else if (F ){# 第二種方法，自己生成 "RA" ,times=13 ),"control" ,times=9 ))"RA" ,"control" ),times = c(13 ,9 ))else if (T ){# 第三種方法，使用字符串出理的函數(shù)獲取分組 "control" ),"control" ,"RA" )# 需要把Group轉(zhuǎn)換成因子，并設(shè)置參考水平，指定levels，對照組在前，處理組在后 "control" ,"RA" ))#2.探針注釋的獲取----------------- #捷徑 library (tinyarray) # 這個(gè)包為生信技能樹老師所寫 'GPL570' )#四種方法，方法1里找不到就從方法2找，以此類推。 #方法1 BioconductorR包(最常用) #http://www./1399.html if (!require (hgu133plus2.db))BiocManager::install("hgu133plus2.db" )library (hgu133plus2.db)"package:hgu133plus2.db" )# 方法2 讀取GPL網(wǎng)頁的表格文件，按列取子集 ##https://www.ncbi.nlm./geo/query/acc.cgi?acc=GPL570 if (F ){#注：表格讀取參數(shù)、文件列名不統(tǒng)一，活學(xué)活用，有的表格里沒有symbol列，也有的GPL平臺沒有提供注釋表格 "GPL570-55999.txt" ,header = T ,"\"" ,sep = "\t" ,check.names = F )"ID" ,"Gene Symbol" )]"probe_id" ,"symbol" )"" & !str_detect(ids2$symbol,"http:///" ),]# 方法3 官網(wǎng)下載注釋文件并讀取 ##http://www./support/technical/byproduct.affx?product=hg-u133-plus # 方法4 自主注釋 #https://mp.weixin.qq.com/s/mrtjpN8yDKUdCSvSUuUwcA "step2output.Rdata" )
Q：畫PCA和熱圖需要使用什么樣的數(shù)據(jù)，使用什么函數(shù)呢？
A：（1）PCA：加載FactoMineR和factoextra包，使用PCA（）和 fviz_pca_ind（）函數(shù)；數(shù)據(jù)：需要對exp矩陣進(jìn)行t轉(zhuǎn)換，將行名設(shè)置為樣本名，列名設(shè)置為基因名，并轉(zhuǎn)換成數(shù)據(jù)框的形式，此外，這里還需要輸入分組信息。
（2）熱圖：加載pheatmap()包，數(shù)據(jù)一般使用log(count+1),這樣畫出來的圖較顯著。
rm(list = ls()) "step1output.Rdata" )"step2output.Rdata" )#輸入數(shù)據(jù)：exp和Group #Principal Component Analysis #http://www./english/articles/31-principal-component-methods-in-r-practical-guide/112-pca-principal-component-analysis-essentials # 1.PCA 圖---- library (FactoMineR)library (factoextra) FALSE )"point" , # show points only (nbut not "text") # color by groups "#00AFBB" , "#E7B800" ),TRUE , # Concentration ellipses "Groups" "pca_plot.Rdata" )# 2.top 1000 sd 熱圖---- 1 ,sd)),1000 ))# 直接畫熱圖，對比不鮮明 library (pheatmap)F ,F ,# 按行標(biāo)準(zhǔn)化 F ,F ,"row" ,3 ,3 ,length.out = 100 )# 關(guān)于scale的進(jìn)一步學(xué)習(xí)：zz.scale.R 上面的因?yàn)樾忻腔?，所以對行進(jìn)行標(biāo)準(zhǔn)化，是為了讓基因在不同的樣本中進(jìn)行標(biāo)準(zhǔn)化。
require (stats)# 1.示例數(shù)據(jù) 1 :30 ,30 ), ncol = 6 )"gene" ,1 :5 )"sample" ,1 :6 )# 2.標(biāo)準(zhǔn)化 #函數(shù)只能按列標(biāo)準(zhǔn)化，但是我們需要按行 # 3.標(biāo)準(zhǔn)化前后，某gene的表達(dá)量點(diǎn)圖比較，大小趨勢不變。 2 ,2 ))1 ,])1 ,])2 ,])2 ,])Q：為什么要去重，各種去重方法對結(jié)果有什么影響 A：去重是因?yàn)樾忻胁荒苡兄貜?fù)值，所以需要先將基因名去重，然后再將這個(gè)基因名作為行名。各種去重方法沒有好壞的定論，一般都可以使用
rm(list = ls())"step2output.Rdata" )#保留最大值 T ),]# 拿這個(gè)ids去inner_join #求平均值 "step2output.Rdata" )1 :4 ,1 :4 ]# 此時(shí)拿到的exp4已經(jīng)是一個(gè)基因?yàn)樾忻谋磉_(dá)矩陣，直接差異分析，不再需要inner_join 在這個(gè)部分才進(jìn)行id轉(zhuǎn)換，不過也可以提到熱圖之前，不過在求差異基因后，再進(jìn)行ID轉(zhuǎn)換，熱圖可以直接畫上調(diào)和下調(diào)基因了
rm(list = ls()) "step2output.Rdata" )#差異分析，用limma包來做 #需要表達(dá)矩陣和Group，不需要改 library (limma)2 ,number = Inf )#為deg數(shù)據(jù)框添加幾列 #1.加probe_id列，把行名變成一列 library (dplyr)#2.加上探針注釋 #其他去重方式在zz.去重方式.R "probe_id" )#3.加change列,標(biāo)記上下調(diào)基因 1 0.05 "down" ,ifelse(k2,"up" ,"stable" )))#4.加ENTREZID列，用于富集分析（symbol轉(zhuǎn)entrezid，然后inner_join） library (clusterProfiler)library (org.Hs.eg.db)"SYMBOL" ,"ENTREZID" ,#人類 #其他物種http:///packages/release/BiocViews.html#___OrgDb "symbol" ="SYMBOL" ))"step4output.Rdata" )Q1：畫火山圖需要什么數(shù)據(jù) A1：需要差異分析后的數(shù)據(jù)，即DESeq2、edgeR、limma分析后的數(shù)據(jù)，需要使用logFC、P.Value。需要加載ggplot2包
Q2：如何畫基因的相關(guān)性圖？ A2：需要加載corrplot包，然后篩選自己想要的基因和它在各組的表達(dá)量，M = cor(t(exp[g,]))，具體看代碼
Q3：如何拼圖？ A3：如果使用ggplot2畫出來的圖，可以加載patchwork包，如果是其他，可以使用plot_grid（）函數(shù)，具體如下
rm(list = ls()) "step1output.Rdata" )"step4output.Rdata" )#1.火山圖---- library (dplyr)library (ggplot2)0.4 , size=3.5 , "-log10(Pvalue)" )+"blue" , "grey" ,"red" ))+4 ,col="black" ,lwd=0.8 ) +4 ,col="black" ,lwd=0.8 ) +in % c("HADHA" ,"LRRFIP1" ))3 , shape = 1 , data = for_label) +"black" #2.差異基因熱圖---- 'step2output.Rdata' )# 表達(dá)矩陣行名替換 if (F ){#全部差異基因 "stable" ]else {#取前10上調(diào)和前10下調(diào) library (dplyr)"stable" ) %>%10 ),10 ))#差異基因熱圖 library (pheatmap)F ,"row" ,#cluster_cols = F, 3 ,3 ,length.out = 100 )#拼圖 library (patchwork)library (ggplotify)# f <- volcano_plot +as.ggplot(heatmap_plot) # ggsave(f,filename = "./pin.pdf",width = 15,height = 15) # 3.感興趣基因的相關(guān)性---- library (corrplot)1 :10 ] # 換成自己感興趣的基因 library (paletteer)"RColorBrewer::RdYlBu" ))10 )"upper" ,"pie" ,"hclust" , "black" , 45 )library (cowplot)# 拼圖 "pca_plot.Rdata" )#保存 "deg.pdf" )#富集分析所有圖表默認(rèn)都是用p.adjust,富集不到可以退而求其次用p值，在文中說明即可
Ｑ：富集分析需要使用到什么數(shù)據(jù)
Ａ：（１）如果僅僅是畫條帶圖和氣泡圖，那么只使用基因的ENTREZID值即可；（２）如果需要畫弦圖（Goplot包），需要Term，logFC，symbol
（3）畫Heatmap-like functional classification，Gene-Concept Network，需要geneList 和ego。
（4）畫。
rm(list = ls()) 'step4output.Rdata' )library (clusterProfiler)library (ggthemes)library (org.Hs.eg.db)library (dplyr)library (ggplot2)library (stringr)library (enrichplot)# 1.GO 富集分析---- #(1)輸入數(shù)據(jù) 'up' ] 'down' ] #(2)富集 #以下步驟耗時(shí)很長，設(shè)置了存在即跳過 if (!file.exists(paste0(gse_number,"_GO.Rdata" ))){"ALL" ,TRUE )"BP" ,TRUE )#ont參數(shù)：One of "BP", "MF", and "CC" subontologies, or "ALL" for all three. "_GO.Rdata" ))"_GO.Rdata" )) ## 這個(gè)很有用 #(3)可視化 #條帶圖 #氣泡圖 "ONTOLOGY" , font.size = 10 , 5 ) + "free_y" ,scales = "free_y" ) + function (x) str_wrap(x, width = 45 ))#geneList 用于設(shè)置下面圖的顏色 #(3)展示top通路的共同基因，要放大看。 #Gene-Concept Network "pvalue" , foldChange=geneList,colorEdge = TRUE )3 ,foldChange=geneList, circular = TRUE , colorEdge = TRUE )#Enrichment Map,這個(gè)函數(shù)最近更新過，版本不同代碼會(huì)不同 "enrichplot" )if (T ){#新版本 else {#舊版本 #(4)展示通路關(guān)系 https://zhuanlan.zhihu.com/p/99789859 #goplot(ego) #(5)Heatmap-like functional classification 8 )# 2.KEGG pathway analysis---- #上調(diào)、下調(diào)、差異、所有基因 #（1）輸入數(shù)據(jù) 'up' ,'ENTREZID' ] 'down' ,'ENTREZID' ] #（2）對上調(diào)/下調(diào)/所有差異基因進(jìn)行富集分析 if (!file.exists(paste0(gse_number,"_KEGG.Rdata" ))){'hsa' )'hsa' )'hsa' )"_KEGG.Rdata" ))"_KEGG.Rdata" ))#(3)看看富集到了嗎？https://mp.weixin.qq.com/s/NglawJgVgrMJ0QfD-YRBQg 0.05 )0.05 )0.05 )#(4)雙向圖 # 富集分析所有圖表默認(rèn)都是用p.adjust,富集不到可以退而求其次用p值，在文中說明即可 0.05 ) %>% #篩選行 1 ) #新增列 0.05 ) %>%1 )source ("kegg_plot_function.R" )#g_kegg +scale_y_continuous(labels = c(2,0,2,4,6)) # 3.GSEA作kegg和GO富集分析---- #https://www./xiaojiewanglezenmofenshen/dbwkg1/ytawgg #(1)查看示例數(shù)據(jù) "DOSE" )#(2)將我們的數(shù)據(jù)轉(zhuǎn)換成示例數(shù)據(jù)的格式 T )#(3)GSEA富集 'hsa' ,FALSE )0.05 & kk_gse$enrichmentScore < 0 ,];down_kegg$group=-1 0.05 & kk_gse$enrichmentScore > 0 ,];up_kegg$group=1 #(4)可視化 # 4.能看懂的資料越來越多---- # GSEA學(xué)習(xí)更多：https://www./docs/share/a67a180f-dd2b-4f6f-96c2-68a4b86fe862?# # 富集分析學(xué)習(xí)更多：http:///clusterProfiler-book/index.html # 弦圖：https://www./xiaojiewanglezenmofenshen/dbwkg1/dgs65p # GOplot：https://mp.weixin.qq.com/s/LonwdDhDn8iFUfxqSJ2Wew # 網(wǎng)上的資料和寶藏?zé)o窮無盡，學(xué)好R語言慢慢發(fā)掘~ # 弦圖 1 :10 ,c(1 ,2 ,8 ,6 )] "/" ,"," ) "ID" ,"Term" ,"Genes" ,"adj_pval" )"BP" # 開始畫弦圖 library (GOplot)# 0.01 , gene.order = 'logFC' , gene.space = 0.15 , gene.size = 3 )"./xiantu.pdf" ,width = 17 ,height = 18 )kegg_plot_function.R里面的代碼
### --------------- ### ### Create: Jianming Zeng ### Date: 2018-07-09 20:11:07 ### Email: jmzeng1314@163.com ### Blog: http://www./ ### Forum: http://www./thread-1376-1-1.html ### CAFS/SUSTC/Eli Lilly/University of Macau ### Update Log: 2018-07-09 First version ### ### Update: Xiaojie Liang ### Date: 2021-07-12 23:38:07 ### Email: liangxiaojiecom@163.com ### --------------- library (ggthemes)library (ggplot2)function (up.data,down.data){F ),]0 )+1 )),stat="identity" , width=0.7 , position=position_dodge(0.7 )) +"#0072B2" , "#B20072" ), guide="none" ) +"" , y="" ) +#axis.ticks.x = element_blank(), 0.3 , colour = "black" ),#x軸連線 0.20 , "cm" ),#修改x軸刻度的高度，負(fù)號表示向上 0.3 , unit = "cm" )),##線與數(shù)字不要重疊 "black" ,size = 0.3 ) ,#修改x軸刻度的線 0 ),NULL , colour = 'white' )
Heatmap-like functional classification
KEGG富集分析（針對下調(diào)和上調(diào)基因）
rm(list = ls())"step4output.Rdata" )'up' ,'symbol' ] 'down' ,'symbol' ]# 1.制作string的輸入數(shù)據(jù) "diffgene.txt" ,F ,F ,F )# 從string網(wǎng)頁獲得string_interactions.tsv # 2.準(zhǔn)備cytoscape的輸入文件 "stable" ,"symbol" ,"logFC" )]"deg.txt" ,"\t" ,F ,F )# string_interactions.tsv是網(wǎng)絡(luò)文件 # deg.txt是屬性表格 
