RNA-seq分析（四）pathway analysis

1. Gene_ID转换

手动把差异表格中基因那一列（如gene-Q0020，gene-替换掉 ），在gene_id那一列加上列名ENSEMBL，另存为csv格式再上传至服务器。

#进入R
load("diff.RData")
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("org.Sc.sgd.db")
BiocManager::install("clusterProfiler")
#首次使用需要install

这时候如果用系统的R，可能安装不了"org.Sc.sgd.db"，可以用conda下载R4.0到自己服务器，org.Sc.sgd.db就可以正常加载啦。如果把自己的R4.0加载到环境变量了，后续使用要注意R的环境，会不会和系统R时有冲突。

library(org.Sc.sgd.db)
library(clusterProfiler)
upTable <- read.csv("SY14_up_2.csv", header = TRUE)
head(upTable)
aTable <- bitr(upTable[,1],fromType = 'ENSEMBL', toType = c('ENTREZID','GENENAME'), OrgDb = 'org.Sc.sgd.db')
#转换
tTable <- merge(aTable,upTable, by= "ENSEMBL")
#合并
head(tTable)
write.csv(tTable,file ="up_symbol.csv")
#输出up_symbol，最后两列增加了ENTREZID，GENENAME
downTable <- read.csv("SY14_down_2.csv", header = TRUE)
head(downTable)
bTable <- bitr(downTable[,1],fromType = 'ENSEMBL', toType = c('ENTREZID','GENENAME'), OrgDb = 'org.Sc.sgd.db')
uTable <- merge(bTable,downTable, by="ENSEMBL")
head(uTable)
write.csv(uTable,file = "down_symbol.csv")
#同理输出down_symbol

2. GO富集分析

上调基因

vi go_up.R
#!bin/bash
library(org.Sc.sgd.db)
library(clusterProfiler)
UP <- read.csv("up_symbol.csv", header=TRUE)
head(UP)
a<- UP[,3]
head(a)
GO_UP <- enrichGO(UP[,3], keyType = "ENTREZID",  OrgDb = "org.Sc.sgd.db", ont = "BP", pAdjustMethod = "BH", pvalueCutoff  = 0.01, qvalueCutoff  = 0.05)
head(GO_UP)
#<0 rows> (or 0-length row.names)
#由于差异基因太少，未富集到

下调基因

padj由0.001调整为0.05，再做GO分析。

down <- read.csv("down_symbol.csv", header=TRUE)
head(down)
#a<- down[,3]
#head(a)
GO_down <- enrichGO(down[,3], keyType = "ENTREZID",  OrgDb = "org.Sc.sgd.db", ont = "BP", pAdjustMethod = "BH", pvalueCutoff  = 0.01, qvalueCutoff  = 0.05)
head(GO_down)
pdf("GO_DOWN.pdf",width = 25, height = 25)
dotplot(GO_down, showCategory=10,title="EnrichmentGO")
dev.off()
#showCategory指定展示的GO Terms的个数，默认展示显著富集的top10个，即p.adjust最小的10个。

3 GSEA (Gene Set Enrichment Analysis) 基因集富集分析

GSEA 是一种计算方法，使用预先定义的基因集gene set（比如想关注基因是否参与DNA REPAIR，可选择hallmark gene set），将基因按照在两类样本中的FoldChange排序得到gene list（按照差异表达倍数从大到小排序），观察预先定义的基因集是在这个gene list的顶端还是底端富集。若参与某通路的基因密集排列在gene list顶端（Leading edge subset），即显著上调，排列在底端即显著下调。

3.1 准备排序后的gene list

BiocManager::install("msigdbr")
#msigdbr 包提供多个物种的MSigDB (Explore the Molecular Signatures Database)数据，是注释基因集的集合
BiocManager::install("dplyr")
#dplyr是R语言的数据分析包，能对dataframe类型的数据做很方便的数据处理和分析操作。
library(clusterProfiler)
library(org.Sc.sgd.db)
library(msigdbr)
library(dplyr)
library(ggplot2)
library(stringr)
exp <- read.csv("SY14_VSBY4742_2.csv", header=TRUE)
head(exp)
gene_ID <- bitr(exp$ENSEMBL, fromType = 'ENSEMBL', toType = c('ENTREZID','GENENAME'), OrgDb = 'org.Sc.sgd.db')
#gene_ID转换 
head(gene_ID)
dim(gene_ID)
#5915    3
diff <-merge(exp,gene_ID, by = "ENSEMBL")
head(diff)
glist <- diff$log2FoldChange
head(glist)
names(glist) = diff$ENTREZID #定义好glist，再输入这一列，否则报错行数不一致
head(glist)
glist = sort(glist, decreasing = T)
head(glist) 
#准备好的genelist按ENTREZID和FoldChange排序

3.2 准备功能集 gene set

The MSigDB gene sets are divided into 9 major collections：H, C1 ~ C8.
H, hallmark gene sets, 聚合许多MSigDB基因集来表达明确的生物状态或过程而获得的一些特征。

Sc <- msigdbr(species = "Saccharomyces cerevisiae")
table(Sc$gs_cat)
#查看目录
Sc[str_detect(Sc$gs_name,"HALLMARK_DNA_REPAIR"),]
#查看Sc中HALLMARK_DNA_REPAIR基因集的gene
#A tibble: 97 × 18
gs_cat gs_subcat     gs_name            gene_symbol  entrez_gene ensembl_gene
<chr>    <chr>         <chr>             <chr>        <int>      <chr>
 1 H      ""        HALLMARK_DNA_REPAIR  AAH1        855581      YNL141W
 2 H      ""        HALLMARK_DNA_REPAIR  ADK2        856917      YER170W
 3 H      ""        HALLMARK_DNA_REPAIR  APT1        854986      YML022W
hallmark <- msigdbr(species = "Saccharomyces cerevisiae",category = "H") %>% dplyr::select(gs_name,entrez_gene)
#％>％来自dplyr包的管道函数，其作用是将前一步的结果直接传参给下一步的函数，
#select(gs_name,entrez_gene)，筛选gs_name和entrez_gene之间的所有列，
head(hallmark)
#gs_name               entrez_gene
  <chr>                       <int>
1 HALLMARK_ADIPOGENESIS      851013
2 HALLMARK_ADIPOGENESIS      852667
dim(hallmark)
#[1] 2759    2
gsea <- GSEA(glist, TERM2GENE = hallmark,verbose=FALSE, pvalueCutoff = 0.2)
head(gsea)
#查看富集到的geneset
library(enrichplot)
pdf("gsea_INTERFERON_ALPHA_RESPONSE.pdf",width = 20, height = 15)
#这里选INTERFERON_ALPHA_RESPONSE 基因集作图
gseaplot2(gsea, geneSetID="HALLMARK_INTERFERON_ALPHA_RESPONSE",color = "firebrick",rel_heights=c(1.5, 0.5, 1),subplots=1:3,pvalue_table = TRUE,ES_geom = "line" )
dev.off()

core_enrichment，即leading edge subset，定义其中对Enrichment score贡献最大的基因为核心基因。
我选的一个基因集结果~~

image.png