干扰MYC-WWP1通路重新激活PTEN的抑癌活性——3步搞定GSEA分析
你好啊，欢迎来到周一数据挖掘专栏，我是新学徒Christine，很高兴和你一起学习！
以后将由我负责带领大家进行数据挖掘实战训练。
菜鸟团（周一数据挖掘专栏）成果展 （如果大家感兴趣前面学徒的教程，可以收藏学习）
继上周的免疫治疗明星分子PD-L1，这周我们来看癌症中另一个“闪耀的明星”PTEN的新研究，文章标题是：Reactivation of PTEN tumor suppressor for cancer treatment through inhibition of a MYC-WWP1 inhibitory pathway，上个月发表于Science。。
文章简介
PTEN（Phosphatase And Tensin Homolog）是一个重要的抑癌基因，编码的蛋白具有蛋白磷酸酶和脂质磷酸酶活性，能拮抗PIK3-AKT信号通路，调控细胞增殖、生长和代谢。PTEN在多种肿瘤中发生高频突变，但通常没有发生完全失去活性，多表现为单等位基因上的功能缺失，亚细胞定位异常，或者发生特殊的翻译后修饰，因为完全的PTEN缺失会给癌细胞带来衰老问题，当然这也给靶向PTEN的治疗带来了机会。
研究目的：已知PTEN发挥功能需要细胞质膜上形成二聚体，但原因未知。作者希望找到调控PTEN二聚化及膜定位的机制，进而研究出恢复PTEN活性的方式。
文章大致思路：
作者对PTEN进行共免疫沉淀，找到一个与PTEN直接物理互作的蛋白WWP1（E3泛素连接酶）
证明WWP1引起PTEN聚泛素化，从而抑制PTEN的二聚化和膜定位，使其不能发挥抑癌作用
发现并验证MYC-WWP1-PTEN调控信号链
通过结构模拟和生化分析得到一个WWP1的潜在抑制剂IC3（来源于十字花科蔬菜）
本次我们复现的图是Fig4F：野生型和Wwp1敲除小鼠RNA-seq的GSEA分析图，用于说明Wwp1敲除影响PI3K-AKT信号通路。
Step1. 下载数据
野生型和Wwp1敲除小鼠的RNA-seq数据存放在GSE126789，但是作者没有给整理好的Matrix文件，需要下载GSE126789_RAW.tar自己整理一下。
rm(list = ls())
options(stringsAsFactors = F)
# 下载GSE126789_RAW.tar
# 批量读取文件处理为表达矩阵
dir <- "./raw_data/GSE126789_RAW/"
files <- list.files(path = dir, pattern = "*wwp1.*.txt.gz", recursive = T) #找到wwp1的WT和敲除样本
expr <- lapply(files,
function(x){
# 只读取基因名和count
expr <- read.table(file = file.path(dir,x), sep = "\t", header = T, stringsAsFactors = F)[,c(1,7)]
return(expr)
})
df <- do.call(cbind, expr)
rownames(df) <- df[,1]
df <- df[,seq(2,ncol(df),by=2)] #去掉重复读取的基因名
colnames(df) <- substr(files,1,10) #列名为样本名
df <- df[apply(df, 1, sum)!=0,] #去掉在所有样本中count都为0的基因
grouplist <- c(rep("KO",3),rep("WT",3)) #记录下分组信息
expr <- df
save(expr,grouplist,file = "./expr_group.Rdata")
Step2. 差异表达分析
上一步中得到的是原始的count，可以直接用Deseq2进行差异分析，用edgeR或limma中的voom方法也是可以的，但是要记得标准化表达数据。
## 差异分析
#原始的count可以用DESeq2做差异分析，表达矩阵需要为integer
exprSet <- apply(expr, 2, function(x){as.integer(x)})
rownames(exprSet) <- rownames(expr)
library(DESeq2)
colData <- data.frame(row.names=colnames(exprSet),
dds <- DESeqDataSetFromMatrix(countData = exprSet,
colData = colData,
design = ~ group_list)
#使用DESeq函数进行估计离散度，然后进行标准的差异表达分析，得到res对象结果
dds <- DESeq(dds)
res <- results(dds,
contrast=c("group_list","KO", "WT")) #提取差异分析结果
DEG <- as.data.frame(res)
DEG <- na.omit(DEG)
nrDEG <- DEG[,c(2,6)] #第6列是padj
Step3. GSEA分析
GSEA：Gene Set Enrichment Analysis，一种基于基因集的富集方法，简单来说就是使用预定义的基因集（一般来自MSigDB数据库），将基因按照某种指标排序，查看这些基因在关注的基因集中是否出现显著统计学富集。
关于GSEA需要了解的内容很多，我把一些相关的资料放在了文末，大家一起学习！
这里我直接用了clusterProfiler包，按照log2FoldChange对基因排序，结果看起来还可以。
library(clusterProfiler)
BiocManager::install("org.Mm.eg.db") # 下载小鼠的OrgDB包
library("org.Mm.eg.db")
#将ENSEMBL ID转换为ENTREZID，用于GSEA分析
eg <- bitr(geneID = rownames(nrDEG), fromType = "ENSEMBL", toType = "ENTREZID",
OrgDb = org.Mm.eg.db)
nrDEG$EZTREZID <- eg[match(rownames(nrDEG),eg$ENSEMBL),2]
nrDEG <- na.omit(nrDEG[order(nrDEG$log2FoldChange,decreasing = T),])
# 得到一个按log2FoldChange排序的genelist
gene_list <- nrDEG$log2FoldChange
names(gene_list) <- nrDEG$EZTREZID
# GSEA分析
kk <- gseKEGG(gene_list, organism = "mmu")
gesa_res <- kk@result
# 画PI3K-AKT的GSEA图，编号为“mmu04151”
gseaplot(kk, "mmu04151", by = "runningScore",
title = "KO vs WT \nKEGG PI3K/AKT Signaling Pathway")
和文章中的图基本一致，文中Fig5I还有2张类似的GSEA图，只是选用的样本不同，有兴趣的小伙伴可以做一下试试~
最后，附上GSEA分析的资料，和你一起学习 ~
GSEA的统计原理比较复杂：GSEA的统计学原理试讲
运行GSEA的方法很多：GSEA分析一文就够（单机版+R语言版）
GSEA中的基因排序也是很有讲究的：
到底该选什么样的基因做GSEA：GSEA分析合理性讨论
生信爆款入门培训
正是因为简单，所以需要指导，毕竟万事开头难，那么你肯定会需要最正宗 生信工程师级别培训！
全国巡讲武汉和成都站报名继续
请联系技能树小助手，扫描下方二维码添加微信~