1

0.安裝R包

options("repos"="https://mirrors.ustc.edu.cn/CRAN/")
if(!require("BiocManager")) install.packages("BiocManager",update = F,ask = F)
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")

cran_packages <- c('tidyr',
                   'tibble',
                   'dplyr',
                   'stringr',
                   'ggplot2',
                   'ggpubr',
                   'factoextra',
                   'FactoMineR',
                   'devtools',
                   'patchwork') 
Biocductor_packages <- c('GEOquery',
                         'hgu133plus2.db',
                         'ggnewscale',
                         "KEGG.db",
                         "limma",
                         "impute",
                         "GSEABase",
                         "GSVA",
                         "clusterProfiler",
                         "org.Hs.eg.db",
                         "preprocessCore",
                         "enrichplot",
                         "ggplotify")

for (pkg in cran_packages){
  if (! require(pkg,character.only=T) ) {
    install.packages(pkg,ask = F,update = F)
    require(pkg,character.only=T) 
  }
}


for (pkg in Biocductor_packages){
  if (! require(pkg,character.only=T) ) {
    BiocManager::install(pkg,ask = F,update = F)
    require(pkg,character.only=T) 
  }
}

#前面的所有提示和報錯都先不要管。主要看這里
for (pkg in c(Biocductor_packages,cran_packages)){
  require(pkg,character.only=T) 
}
#沒有error就是成功！

#哪個報錯，就回去安裝哪個。如果你沒有安裝xx包，卻提示你xx包不存在，這也正常，是因為復(fù)雜的依賴關(guān)系，缺啥補啥。
if(!require(AnnoProbe))devtools::install_local("./AnnoProbe-master.zip",upgrade = F)
library(AnnoProbe)

#上次遇到的報錯解決：
BiocManager::install("graphlayouts")
library(clusterProfiler)
library(enrichplot)

1.下載數(shù)據(jù)，提取表達矩陣和臨床信息

#數(shù)據(jù)下載
rm(list = ls())#清空環(huán)境
options(stringsAsFactors = F)#避免因子影響
library(GEOquery)#加載要用的包
gse_number = "GSE56649"#給gse編號一個變量，后面只需修改這里就好
eSet <- getGEO(gse_number, 
               destdir = '.', 
               getGPL = F)#下載且讀取，destdir從哪讀取文件"."代表當(dāng)前目錄;getGPL意思是下載matrix文件時要不要一起下載GPL注釋
class(eSet)#看是什么數(shù)據(jù)類型
length(eSet)#看長度
eSet = eSet[[1]]#簡化并重新賦值為簡化版
#eSet@phenoData#提取臨床信息
eSet@annotation#返回這個表達矩陣用哪個平臺測序的

#(1)提取表達矩陣exp
exp <- exprs(eSet)
exp[1:4,1:4]
#range(exp)#可以查看矩陣取值范圍，看有沒有l(wèi)og過
exp = log2(exp+1)#確認好有沒有取過log，沒取過再取，避免矩陣?yán)镉辛慵拥?（不改）甲基化要加更小的值
boxplot(exp)#這里也用來看范圍
#(2)提取臨床信息
pd <- pData(eSet)
#(3)調(diào)整pd的行名順序與exp列名完全一致
p = identical(rownames(pd),colnames(exp));p#判斷pd的行名和exp是否一致
if(!p) exp = exp[,match(rownames(pd),colnames(exp))]
#(4)提取芯片平臺編號
gpl_number <- eSet@annotation
save(gse_number,pd,exp,gpl_number,file = "step1output.Rdata")#輸出結(jié)果

• 檢查數(shù)據(jù)下載的完整性

  
  
  看兩個KB或MB大小是否一樣

2.分組信息&探針注釋

# Group(實驗分組)和ids(探針注釋)
rm(list = ls())  
load(file = "step1output.Rdata")
library(stringr)
# 1.Group----
# 第一類，有現(xiàn)成的可以用來分組的列
if(F) Group = pd$`disease state:ch1` 

#第二類，自己生成
if(F){
  Group=c(rep("RA",times=13),
          rep("control",times=9))
}#（1）
rep(c("RA","control"),times = c(13,9))#與(1)相同，簡單寫法

#第三類，**匹配關(guān)鍵詞，自行分類
Group=ifelse(str_detect(pd$source_name_ch1,
            "control"),"control","RA")

#設(shè)置參考水平，指定levels，對照組在前，處理組在后
Group = factor(Group,
               levels = c("control","RA"))
Group

# 注意levels與因子內(nèi)容必須對應(yīng)一致
# Group = pd$`disease state:ch1`
# Group = factor(Group,
#                 levels = c("healthy control","rheumatoid arthritis"))

#2.ids-----------------  

#方法1 BioconductorR包(最常用)
gpl_number #用這個GPL編號在以下網(wǎng)址搜索，第三列包名字的前綴，后加.db
#http://www.bio-info-trainee.com/1399.html

#安裝并加載這個包
if(!require(hgu133plus2.db))BiocManager::install("hgu133plus2.db")
library(hgu133plus2.db)
ls("package:hgu133plus2.db")#看包里有什么
ids <- toTable(hgu133plus2SYMBOL)#toTable是變成數(shù)據(jù)框
#SYMBOL代表探針和symbol之間的關(guān)系
head(ids)

# 方法2 讀取GPL平臺的soft文件，按列取子集
##https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL570
if(F){
  #注：soft文件列名不統(tǒng)一，活學(xué)活用，有的GPL平臺沒有提供注釋，如GPL16956
  a = getGEO(gpl_number,destdir = ".")
  b = a@dataTable@table
  colnames(b)
  ids2 = b[,c("ID","Gene Symbol")]
  colnames(ids2) = c("probe_id","symbol")
  ids2 = ids2[ids2$symbol!="" & #空字符串的去掉
      !str_detect(ids2$symbol,"http:///"),]#一個探針對應(yīng)多個基因?qū)儆诜翘禺愋蕴结槪サ?}
# 方法3 官網(wǎng)下載,文件讀取
##http://www.affymetrix.com/support/technical/byproduct.affx?product=hg-u133-plus
# 方法4 自主注釋 
#https://mp.weixin.qq.com/s/mrtjpN8yDKUdCSvSUuUwcA
save(exp,Group,ids,gse_number,file = "step2output.Rdata")

• 當(dāng)一列中一部分相同，一部分不相同，取不同部分時

pd$characteristics_ch1
#兩種方法
str_remove(pd$characteristics_ch1,"disease state:")
str_split(pd$characteristics_ch1,": ",simplify = T)[,2]

3.數(shù)據(jù)探索(pca)和熱圖

rm(list = ls())  
load(file = "step1output.Rdata")
load(file = "step2output.Rdata")
#輸入數(shù)據(jù)：exp和Group
#Principal Component Analysis
#http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/112-pca-principal-component-analysis-essentials

# 1.PCA 圖----展示兩組之間的差異
dat=as.data.frame(t(exp))#把數(shù)據(jù)轉(zhuǎn)置并變成數(shù)據(jù)框
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"))#保存圖片
save(pca_plot,file = "pca_plot.Rdata")#保存Rdata，為了后面合圖

# 2.top 1000 sd 熱圖---- 從總基因中挑變化比較大的一部分基因
cg=names(tail(sort(apply(exp,1,sd)),1000))
n=exp[cg,]

# 直接畫熱圖，對比不鮮明
library(pheatmap)
annotation_col=data.frame(group=Group)
rownames(annotation_col)=colnames(n) 
pheatmap(n,
         show_colnames =F,
         show_rownames = F,
         annotation_col=annotation_col
)

# 用標(biāo)準(zhǔn)化的數(shù)據(jù)畫熱圖，兩種方法的比較：https://mp.weixin.qq.com/s/jW59ujbmsKcZ2_CM5qRuAg

## 1.使用熱圖參數(shù)
pheatmap(n,
         show_colnames =F,
         show_rownames = F,
         annotation_col=annotation_col,
         scale = "row",#按行標(biāo)準(zhǔn)化，只比較列與列的區(qū)別
         breaks = seq(-3,3,length.out = 100)
         ) #breaks 參數(shù)解讀在上面鏈接，設(shè)置顏色分配范圍
#breaks：小于-3的值顯示和-3一個顏色，大于3顯示和3一個顏色
#        常用-3到3 或 -2到2
dev.off()

## 2.自行標(biāo)準(zhǔn)化再畫熱圖，與上面1畫出的圖沒有區(qū)別
n2 = t(scale(t(n)))#scale只能按列標(biāo)準(zhǔn)化，所以先轉(zhuǎn)置
pheatmap(n2,
         show_colnames =F,
         show_rownames = F,
         annotation_col=annotation_col,
         breaks = seq(-3,3,length.out = 100)
)
dev.off()

# 關(guān)于scale的進一步探索：zz.scale.R

# 3.相關(guān)性熱圖---- 

pheatmap::pheatmap(cor(exp),#看表達矩陣組成新矩陣列與列之間的相關(guān)性
                   annotation_col = annotation_col)
#一般選第二個，或者做完差異分析后拿差異基因做相關(guān)性熱圖比較明顯
pheatmap::pheatmap(cor(n),#top1000差異基因的相關(guān)性，比表達矩陣有意義一些
                   annotation_col = annotation_col)

pheatmap::pheatmap(cor(n2),#標(biāo)準(zhǔn)化后的top1000的相關(guān)性
                   annotation_col = annotation_col
                   )

dev.off()

# 關(guān)于相關(guān)性背后的故事：https://mp.weixin.qq.com/s/IqMW6Qjf64dn30F4RQg5kQ

2

2

數(shù)據(jù)框、矩陣轉(zhuǎn)置后都會變成矩陣，as.data.frame把矩陣變成數(shù)據(jù)框

4.差異分析

rm(list = ls()) 
load(file = "step2output.Rdata")

#差異分析，用limma包來做
#需要表達矩陣和Group，不需要改
library(limma)
design=model.matrix(~Group)#group要求對照組在前，實驗組在后
fit=lmFit(exp,design)#線性擬合，數(shù)據(jù)：表達矩陣和根據(jù)分組信息生成的比較矩陣design
fit=eBayes(fit)#貝葉斯檢驗
deg=topTable(fit,coef=2,number = Inf)

#為deg數(shù)據(jù)框添加幾列
#1.加probe_id列（探針id），把行名變成一列
library(dplyr)
deg <- mutate(deg,probe_id=rownames(deg))
head(deg)

#2.加上探針注釋，每個探針對應(yīng)的基因
table(!duplicated(ids$probe_id))
table(!duplicated(ids$symbol))
#按symbol列去重，常見標(biāo)準(zhǔn)有3個：最大值/平均值/隨機去重
#隨機去重，另兩個見zz.去重方式.R
ids = ids[!duplicated(ids$symbol),]

deg <- inner_join(deg,ids,by="probe_id")#對deg和ids按照基因取交集
head(deg)
nrow(deg)#基因的數(shù)量

#3.加change列,標(biāo)記上下調(diào)基因
logFC_t=1#設(shè)置logFC的閾值
P.Value_t = 0.01#設(shè)置P.Value的閾值
k1 = (deg$P.Value < P.Value_t)&(deg$logFC < -logFC_t)#下調(diào)基因
k2 = (deg$P.Value < P.Value_t)&(deg$logFC > logFC_t)#上調(diào)基因
change = ifelse(k1,
                "down",
                ifelse(k2,"up","stable"))#標(biāo)記上下調(diào)基因
deg <- mutate(deg,change)#添加到deg表格中

#4.加ENTREZID列，用于富集分析（symbol轉(zhuǎn)entrezid，然后inner_join）
library(clusterProfiler)
library(org.Hs.eg.db)
s2e <- bitr(deg$symbol, 
            fromType = "SYMBOL",
            toType = "ENTREZID",
            OrgDb = org.Hs.eg.db)#人類
#其他物種http://bioconductor.org/packages/release/BiocViews.html#___OrgDb
dim(deg)
deg <- inner_join(deg,s2e,by=c("symbol"="SYMBOL"))
dim(deg)
length(unique(deg$symbol))
save(Group,deg,logFC_t,P.Value_t,gse_number,file = "step4output.Rdata")

5.可視化：火山圖和熱圖

rm(list = ls()) 
load(file = "step1output.Rdata")
load(file = "step4output.Rdata")
#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

#在圖上標(biāo)記關(guān)心的基因
if(T){
  #自選基因
  for_label <- dat%>% 
    filter(symbol %in% c("HADHA","LRRFIP1")) 
}
if(F){
  #p值最小的10個
  for_label <- dat %>% head(10)
}
if(F) {
  #p值最小的前3下調(diào)和前3上調(diào)
  x1 = dat %>% 
    filter(change == "up") %>% 
    head(3)
  x2 = dat %>% 
    filter(change == "down") %>% 
    head(3)
  for_label = rbind(x1,x2)
}

volcano_plot <- p +
  geom_point(size = 3, shape = 1, data = for_label) +
  ggrepel::geom_label_repel(#添加標(biāo)記基因的圖層
    aes(label = symbol),
    data = for_label,
    color="black"
  )
volcano_plot
ggsave(plot = volcano_plot,filename = paste0(gse_number,"_volcano.png"))

#2.差異基因熱圖----

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)
}
n=exp[cg,]
dim(n)

#差異基因熱圖
library(pheatmap)
annotation_col=data.frame(group=Group)#變成輸入數(shù)據(jù)要求的樣子
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"))
load("pca_plot.Rdata")
library(patchwork)
library(ggplotify)
(pca_plot + volcano_plot +as.ggplot(heatmap_plot))

6.富集分析

rm(list = ls())  
load(file = 'step4output.Rdata')
library(clusterProfiler)
library(dplyr)
library(ggplot2)
library(stringr)
library(enrichplot)

# 1.GO 富集分析----

#(1)輸入數(shù)據(jù)
gene_up = deg[deg$change == 'up',
              'ENTREZID'] #上調(diào)基因?qū)?yīng)的ENTREZID
gene_down = deg[deg$change == 'down','ENTREZID'] 
gene_diff = c(gene_up,gene_down)
gene_all = deg[,'ENTREZID']

#(2)富集
#以下步驟耗時很長，設(shè)置了存在即跳過
if(!file.exists(paste0(gse_number,"_GO.Rdata"))){
  ego <- enrichGO(gene = gene_diff,#輸入數(shù)據(jù)
                  OrgDb= org.Hs.eg.db,
                  ont = "ALL",
                  readable = TRUE)
  #ont參數(shù)：One of "BP", "MF", and "CC" subontologies, or "ALL" for all three.
  save(ego,file = paste0(gse_number,"_GO.Rdata"))
}
load(paste0(gse_number,"_GO.Rdata"))

#(3)可視化
#條帶圖
barplot(ego)
#氣泡圖
dotplot(ego)

dotplot(ego, split = "ONTOLOGY", font.size = 10, 
        showCategory = 5) + facet_grid(ONTOLOGY ~ ., scale = "free") + 
  scale_y_discrete(labels = function(x) str_wrap(x, width = 45))

#geneList 用于設(shè)置下面圖的顏色
geneList = deg$logFC
names(geneList)=deg$ENTREZID
geneList = sort(geneList,decreasing = T)

#(3)展示top通路的共同基因，要放大看。
#Gene-Concept Network
cnetplot(ego,categorySize="pvalue", foldChange=geneList,colorEdge = TRUE)
cnetplot(ego, showCategory = 3,foldChange=geneList, circular = TRUE, colorEdge = TRUE)
#Enrichment Map,這個函數(shù)最近更新過，版本不同代碼會不同
#showCategory = 3，展示條目的數(shù)量，默認是五
Biobase::package.version("enrichplot")

if(F){
  emapplot(pairwise_termsim(ego)) #新版本
}else{
  emapplot(ego)#老版本
}
#(4)展示通路關(guān)系 https://zhuanlan.zhihu.com/p/99789859
#goplot(ego)

#(5)Heatmap-like functional classification
heatplot(ego,foldChange = geneList,showCategory = 8)

# 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)/所有差異基因進行富集分析

if(!file.exists(paste0(gse_number,"_KEGG.Rdata"))){
  kk.up <- enrichKEGG(gene         = gene_up,
                      organism     = 'hsa')#物種的縮寫
  kk.down <- enrichKEGG(gene         =  gene_down,
                        organism     = 'hsa')
  kk.diff <- enrichKEGG(gene         = gene_diff,
                        organism     = 'hsa')
  save(kk.diff,kk.down,kk.up,file = paste0(gse_number,"_KEGG.Rdata"))
}
load(paste0(gse_number,"_KEGG.Rdata"))

#(3)看看富集到了嗎？https://mp.weixin.qq.com/s/NglawJgVgrMJ0QfD-YRBQg
table(kk.diff@result$p.adjust<0.05)
table(kk.up@result$p.adjust<0.05)
table(kk.down@result$p.adjust<0.05)

#(4)按照pvalue篩選通路
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)可視化
source("kegg_plot_function.R")
g_kegg <- kegg_plot(up_kegg,down_kegg)
g_kegg
#g_kegg +scale_y_continuous(labels = c(4,2,0,2,4))#改橫坐標(biāo)用的，橫坐標(biāo)不應(yīng)該為負
ggsave(g_kegg,filename = 'kegg_up_down.png')

# 3.gsea作kegg和GO富集分析----
## http://www.itdecent.cn/p/c5b7b7dbf29b

#(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',
                  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)可視化
g2 = kegg_plot(up_kegg,down_kegg)
g2

# 4.能看懂的資料越來越多----
# GSEA學(xué)習(xí)更多：http://www.itdecent.cn/p/baf85b51752e
# 富集分析學(xué)習(xí)更多：http://yulab-smu.top/clusterProfiler-book/index.html
# 弦圖：http://www.itdecent.cn/p/e4bb41865b7f
# GOplot：https://mp.weixin.qq.com/s/LonwdDhDn8iFUfxqSJ2Wew
# 網(wǎng)上的資料和寶藏?zé)o窮無盡，學(xué)好R語言慢慢發(fā)掘~

* Annoprobe包三個函數(shù)的使用

#[1]
library(AnnoProbe)
geoChina("gse1009")#下載gse1009文件
load("GSE1009_eSet.Rdata")
eSet=gset[[1]]#即可對接標(biāo)準(zhǔn)流程

#[2]
ids=idmap("GPL570")#獲取ids的一種方法
head(ids)

#[3]
annoGene()#可以知道你提供的基因ID的具體信息是什么

*小潔老師自己寫的R包（簡化）

#devtools::install_github("xjsun1221/tinyarray")
library(tinyarray)
geo = geo_download("GSE56649")

View(geo$pd)
pd = geo$pd
library(stringr)
Group=ifelse(str_detect(pd$source_name_ch1,"control"),
             "control",
             "RA")

#設(shè)置參考水平，指定levels，對照組在前，處理組在后
Group = factor(Group,
               levels = c("control","RA"))
Group

find_anno(geo$gpl,install = T)

ids <- toTable(hgu133plus2SYMBOL)

geo$exp = log2(geo$exp+1)
deg = get_deg_all(geo$exp,Group,ids)
deg$plots

7.常見錯誤

3

8.自行探索部分

(1）配對數(shù)據(jù)分析

[1]

[2]

(2）多分組數(shù)據(jù)

[3]

[4]

(3)多數(shù)據(jù)聯(lián)合分析

[5]

[6]

[7]

[8]

(4)標(biāo)準(zhǔn)流程的后續(xù)

[9]

[10]

[11]

[12]

[13]

[14]

[15]

[16]