「生信技能樹」三陰性乳腺癌表達(dá)矩陣探索 系列筆記
GEO挖掘?qū)崙?zhàn)一恩沛、初步探索數(shù)據(jù) - 簡(jiǎn)書
GEO挖掘?qū)崙?zhàn)二、差異分析及富集分析 - 簡(jiǎn)書
GEO挖掘?qū)崙?zhàn)三已卷、GSVA - 簡(jiǎn)書
GEO挖掘?qū)崙?zhàn)四、TNBC相關(guān)探索 - 簡(jiǎn)書

在之前的步驟中，已經(jīng)準(zhǔn)備好表達(dá)矩陣與分組信息腐螟；并且在兩組分組存在一定區(qū)分度。接下來就可進(jìn)行差異分析困后。得到差異基因結(jié)果后可繪制火山圖乐纸、熱圖；進(jìn)行富集分析(ORF摇予，GSEA)等汽绢。

4、limma差異分析

芯片數(shù)據(jù)的差異分析一般使用limma包

4.1

#可選
#先隨便選擇一個(gè)探針在兩組樣本間的分布進(jìn)行可視化
rm(list = ls())
load('exp_group.Rdata')
library(ggpubr)
bp=function(g){
  df <- data.frame(gene=g,group=group_list)
  ggboxplot(df, x='group', y='gene',
            color = "group", palette = 'jco',
            add = "jitter") +  #散點(diǎn)圖抖動(dòng)
    stat_compare_means()  #添加p值
}
bp(exp[2,])

4-1

4.2侧戴、limma差異分析

library(limma)
design <- model.matrix(~factor(group_list))
fit=lmFit(exp, design)
fit=eBayes(fit)
options(digits = 4)
deg <- topTable(fit,coef=2,adjust='BH',number = Inf)
head(deg,3)
t.test(exp[rownames(deg)[1],]~group_list)
bp(exp[rownames(deg)[1],])

4-2

如上圖宁昭，logFC即logFoldchange，為raw counts變化倍數(shù)的log值酗宋。由于數(shù)據(jù)已經(jīng)是log后的所有是兩組對(duì)應(yīng)的均值的差值积仗。
此外根據(jù)logFC的正負(fù)性，我們能夠判斷這里是TNBC與nonTNBC相比蜕猫。若想修改比較順序斥扛，簡(jiǎn)單的方法是直接取logFC的相反數(shù)即可

4-3

TCGA學(xué)習(xí)02：差異分析 - 簡(jiǎn)書的筆記中是limma配合edgeR包差異分析的pipeline，二者區(qū)別暫時(shí)還未弄清。

4.3稀颁、探針名轉(zhuǎn)化

• 差異分析后芬失，圍繞差異基因，探究其背后的生物意義匾灶。這里就要即將探針名轉(zhuǎn)換成有意義的基因ID
• 芯片探針I(yè)D轉(zhuǎn)換可參考用R獲取芯片探針與基因的對(duì)應(yīng)關(guān)系三部曲-bioconductor | 生信菜鳥團(tuán)

• 首先了解芯片平臺(tái)信息棱烂，再下載對(duì)應(yīng)R包，最后轉(zhuǎn)化ID即可阶女。

  
  
  4.3-1
  
  
  4.3-2

if (!require("hgu133plus2.db"))
  BiocManager::install("hgu133plus2.db")
ls('package:hgu133plus2.db')
ids <- toTable(hgu133plus2SYMBOL)
#得到這個(gè)平臺(tái)里颊糜，所有探針與基因名（symbol格式）的對(duì)應(yīng)關(guān)系
ide <- toTable(hgu133plus2ENTREZID)
head(ids);dim(ids)
head(ide);dim(ide)
deg$probe_id=rownames(deg)
deg <- merge(deg, ids, by='probe_id')
deg <- merge(deg, ide, by='probe_id')
head(deg)

save(deg, file = 'deg.Rdata')
rm(list = ls())

此外知道symbol等標(biāo)準(zhǔn)基因ID后，也可通過clusterProfiler包進(jìn)行基因名間的轉(zhuǎn)換秃踩。例如

4.3-3

5衬鱼、差異分析結(jié)果可視化

5.1、火山圖

load('deg.Rdata')
head(deg)
deg$change=as.factor(
  ifelse(deg$P.Value>0.01, 'stable', 
          ifelse(deg$logFC>1.5,'UP',
                 ifelse(deg$logFC< -1.5,'DOWN','stable'))
         )
  )
table(deg$change)
deg$'-log10(pvalue)' <- -log10(deg$P.Value)
head(deg)
ggscatter(deg, x="logFC", y="-log10(pvalue)",
          color="change",size = 0.7,
          label = "symbol", repel = T, #自動(dòng)調(diào)整標(biāo)簽位置
          label.select = c("AGR3","CYP4Z1","PROM1","SHC4"))

save(deg,file = "deg.Rdata")
rm(list = ls())

5.1

此外也可利用ggplot繪制憔杨，可參考TCGA學(xué)習(xí)02：差異分析 - 簡(jiǎn)書

5.2鸟赫、熱圖

• 繪制熱圖根據(jù)差異探針，選取表達(dá)矩陣的子集消别；
• 這里我們分別選取上調(diào)抛蚤、下調(diào)的top100基因繪圖。

load("deg.Rdata")
head(deg)
tmp <- deg$logFC; names(tmp) <- deg$probe_id
head(tmp)
cg_up <- names(head(sort(tmp,decreasing = T),100))
cg_down <- names(head(sort(tmp),100))

load("exp_group.Rdata")
n=t(scale(t(exp[c(cg_up,cg_down),])))
n[n>2]=2; n[n< -2]= -2
n[1:4,1:4]
group_dat <- data.frame(group=group_list, row.names = colnames(exp))
library(pheatmap)
pheatmap(n,
         show_rownames = F,
         show_colnames = F,
         #cluster_cols = F,
         annotation_col = group_dat)

5.2

6寻狂、富集分析

之前學(xué)習(xí)RNA-seq轉(zhuǎn)錄組學(xué)習(xí)時(shí)岁经，對(duì)富集分析的概念與流程有過一定的了解。主要分為ORF與GESA兩類蛇券，都可用clusterProfiler包完成缀壤。在曾老師的視頻中后者是使用了MsigDB的數(shù)據(jù)集進(jìn)行分析的。
-RNA-seq學(xué)習(xí)：No.5富集分析--ORF過表達(dá) - 簡(jiǎn)書
-RNA-seq學(xué)習(xí)：No.6富集分析--GESA - 簡(jiǎn)書

6.1 ORF過表達(dá)富集分析

主要需要上下調(diào)基因的ENTREZID

rm(list=ls())
library(clusterProfiler)
load('deg.Rdata')
gene_up <- deg[deg$change=="UP","gene_id"]
gene_down <- deg[deg$change=="DOWN","gene_id"]

##kegg.up
kk.up <- enrichKEGG(gene = gene_up,
                    organism = "hsa",
                    pvalueCutoff = 0.9,
                    qvalueCutoff = 0.9)
head(kk.up)
kegg_up_dt <- as.data.frame(kk.up)
##kegg.down
kk.down <- enrichKEGG(gene = gene_down,
                    organism = "hsa",
                    pvalueCutoff = 0.9,
                    qvalueCutoff = 0.9)
head(kk.down)
kegg_down_dt <- as.data.frame(kk.down)

kegg_plot <- function(up_kegg,down_kegg){
dat=rbind(up_kegg,down_kegg)
dat$pvalue <- -log10(dat$pvalue)
dat$pvalue <- dat$pvalue*dat$group
dat=dat[order(dat$pvalue,decreasing = F),]
ggplot(dat, aes(x=reorder(Description,order(pvalue,decreasing= F)),y=pvalue, fill=group)) + 
  #x軸按對(duì)應(yīng)的pvalue值從大到小排列pathway的Description
  geom_bar(stat='identity') +
  #設(shè)置柱狀圖高低直接為數(shù)值大小纠亚，而不是counts
  scale_fill_gradient(low="blue",high = "red", guide = F) +
  scale_x_discrete(name="pathway names") +
  #針對(duì)字符型軸標(biāo)簽注釋
  scale_y_continuous(name="log10P-value") +
  #針對(duì)連續(xù)型軸標(biāo)題注釋
  coord_flip() + theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
  ggtitle("pathway enrichment")
}
up_kegg <- kegg_up_dt[kegg_up_dt$pvalue<0.05,];up_kegg$group <- -1
down_kegg <- kegg_down_dt[kegg_down_dt$pvalue<0.05,];down_kegg$group <- 1
kegg_plot(up_kegg,down_kegg)

• 分別查看上下調(diào)基富集到哪些通路里诉位。

  
  
  6.1

此外還有 go的富集分析大同小異，這里不做介紹菜枷，可參考上面的鏈接。

6.2 GSEA的FCS富集分析（基于MsigDB數(shù)據(jù)庫）

（1）genelist

需要準(zhǔn)備genelist數(shù)值型字符串叁丧，即為logFC值啤誊，從大到小排列；并以ENTREZID/SYMBOL命名拥娄。

rm(list = ls())
load("deg.Rdata")
head(deg);table(deg$change)
dat <- deg[deg$change %in% c("DOWN","UP"), c('logFC','gene_id')]
dim(dat);head(dat)
genelist <- dat$logFC
names(genelist) <- dat$gene_id
genelist <- sort(genelist, decreasing = T)

（2）下載MsigDB

• GSEA | MSigDB https://www.gsea-msigdb.org/gsea/msigdb/index.jsp
  需要注冊(cè)一個(gè)郵箱蚊锹，才能下載。這里我們用的ENTREZID稚瘾，故下載對(duì)應(yīng)的版本即可牡昆，保存到自己工作目錄的一個(gè)子目錄即可。

d <- "./msigdb/"
gmts <- list.files(d, pattern = 'all')
gmts
#[1] "c1.all.v7.1.entrez.gmt" "c2.all.v7.1.entrez.gmt" "c3.all.v7.1.entrez.gmt"
#[4] "c4.all.v7.1.entrez.gmt" "c5.all.v7.1.entrez.gmt" "c6.all.v7.1.entrez.gmt"
#[7] "h.all.v7.1.entrez.gmt" 

（3）富集分析

if (!require("BiocManager"))
  install.packages("BiocManager")
if (!require("GSEABase"))
  BiocManager::install("GSEABase")
gsea_results <- lapply(gmts, function(gmtfile){
  geneset <- read.gmt(file.path(d,gmtfile))
  print(paste0("Now process the ", gmtfile))
  egmt <- GSEA(genelist, 
               TERM2GENE = geneset, 
               pvalueCutoff = 0.05) #默認(rèn)值，導(dǎo)致有的數(shù)據(jù)庫就沒有富集到丢烘，可適當(dāng)調(diào)大)
  head(egmt)
  return(egmt)
})
gsea_results_df <- lapply(gsea_results, function(x){
  cat(paste(dim(x@result)), '\n')
  x@result
}) #得到的list包含各個(gè)result

gsea_results_df <- do.call(rbind,gsea_results_df) #合并list里的所有result

gseaplot(gsea_results[[2]],gsea_results[[2]]@result$ID[3])
gsea_results[[2]]@result$ID[3]

如下圖柱宦，可看出差異基因在這個(gè)geneset里主要富集下調(diào)了。更多解讀可見上文的筆記鏈接播瞳。此外clusterfile包也可做GSEA分析掸刊，主要是go、kegg數(shù)據(jù)庫赢乓。

6.2