search

ixxmu / mp_duty

抓取网络文章到github issues保存
https://archives.duty-machine.now.sh/
110 stars 30 forks source link

批量读取转录组数据，绘制感兴趣的GSEA通路 #3104

Closed ixxmu closed 1 year ago

ixxmu commented 1 year ago

https://mp.weixin.qq.com/s/OOXvUVbQsGhYlrfma0i0Eg

ixxmu commented 1 year ago

批量读取转录组数据，绘制感兴趣的GSEA通路 by 生信菜鸟团

缘起

前两天，曾老师给了我一个3组每组4个样本的转录组数据，原文主要对于该数据进行了GSEA通路分析。老师想让我利用原文提供的FPKM矩阵，验证一下原文的通路变化。本打算对3组之间进行相互之间差异分析的，但是由于作者绘制的GSEA通路分析仅针对WT与KO两组，所以此处只针对这两组进行差异分析，绘制相应的GSEA通路。存在一个小插曲的是，文中绘制的5条GSEA通路中，小编只看见两条，这两条一条是通过gseKEGG获得，另一条是通过自定义通路进行GSEA分析，其余三条在文中没有找到，感兴趣的小伙伴可以尝试下。此外，本文挺有意思的点是，作者提供的基因表达矩阵是按照每个样本提供的，每组样本间的基因数量是不同的，所以需要分别批量读取三个组别，最终将三个组别的基因名取交集，获得相应的基因表达矩阵。

探究

今天，我们使用的转录组数据集是2020年发表在CANCER RESEARCH|GENOME AND EPIGENOME杂志上，文献名称为Loss of Apc Rapidly Impairs DNA Methylation Programs  and Cell Fate Decisions in Lgr5þ Intestinal Stem Cells。该文主要对三个组别进行了转录组数据分析，具体介绍如下。

注意哈，咱们今天要绘制的通路主要是如下原文结果中框起来的两条通路：

转录组数据集介绍

该数据集提交在GEO官网，其GSE号是GSE123006。该数据集转录组数据有3个组别，分别是WT组、HT组与KO组，每个组别含有4个样本。感兴趣的小伙伴可以尝试下，具体数据下载链接见 https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE123006&format=file。注意该下载链接中除了转录组数据，还含有RRBS数据，分析转录组数据时，需要将RRBS数据剔除。

正式分析

1.清空环境，加载包，加载数据

rm(list=ls())
library(data.table)
samples=list.files("GSE123005_RAW/")
dir="GSE123005_RAW/"
#获取三组样本
wt=samples[grep(samples,pattern = "wt")]
ht=samples[grep(samples,pattern = "ht")]
ko=samples[grep(samples,pattern = "ko")]

## 定义一个函数，直接获得每组的基因表达矩阵
group_exp=function(x){
  samples=x #这里锁定组别，应该就可以了
  df=list()
  for(i in 1:4){
    df[[i]]=fread(paste0(dir,samples[i]))[,c(1,5)]
  }
  
  exp=do.call(cbind,df) #；将list转换为数据框，do.call只能取一列
  exp=as.data.frame(exp)
  exp=exp[,c(1,2,4,6,8)]
  table(duplicated(exp$gene_short_name))
  ## 重复数目有些多，按照最大值去重
  exp$meadian=apply(exp[,c(2,5)],1,median) #1.先求中位数
  exp=exp[order(exp$gene_short_name,abs(exp$meadian),decreasing = T),] #2.按order排降序
  exp=exp[!duplicated(exp$gene_short_name),] #3.去重
  rownames(exp)=exp$gene_short_name #4.命基因名
  exp=exp[,2:5]
  return(exp) #最终获取exp的结果
}

wt_exp=group_exp(wt)
ht_exp=group_exp(ht)
ko_exp=group_exp(ko)

#将每组的基因表达矩阵整合成为一个大的基因表达矩阵
intersect_gene=intersect(intersect(rownames(wt_exp),rownames(ht_exp)),rownames(ko_exp))
wt_exp=wt_exp[intersect_gene,]
ht_exp=ht_exp[intersect_gene,]
ko_exp=ko_exp[intersect_gene,]
all_exp=cbind(wt_exp,ht_exp,ko_exp)

2.对WT组与KO组进行差异分析，获得GSEA排序所需要的logFC

expMatrix=all_exp[,c(1:4,9:12)] #分析WT组与KO组就可以了
# 2.将FPKM值转换为TPM值
fpkmToTpm <- function(fpkm)
{
  exp(log(fpkm) - log(sum(fpkm)) + log(1e6))
}
class(expMatrix)
## [1] "data.frame"
## [1] "data.frame"
expMatrix=na.omit(expMatrix)
numer=function(x){as.numeric(x)}
expMatrix=apply(expMatrix, 1, numer) 
##必须将矩阵中的字符类型由向量型改为字符型
expMatrix=t(expMatrix)
tpms <- apply(expMatrix,2,fpkmToTpm) 
tpms[1:3,]
##            [,1]     [,2]     [,3]      [,4]      [,5]      [,6]      [,7]
## Zzz3   42.12465 23.35684 34.31037  28.13509  54.89502  41.10471  47.76531
## Zzef1  40.00925 28.08166 43.80476  42.77359  18.23898  12.34853  15.04639
## Zyx   103.04799 83.06527 73.74318 103.07322 320.04760 190.79061 225.77413
##             [,8]
## Zzz3   25.362704
## Zzef1   9.756136
## Zyx   192.170262
colSums(tpms)
## [1] 1e+06 1e+06 1e+06 1e+06 1e+06 1e+06 1e+06 1e+06
colnames(tpms)=c("WT1","WT2","WT3","WT4","KO1","KO2","KO3","KO4")
## 组别设置
group_list=c(rep('WT',4),rep('KO',4))
group_list <- factor(group_list,levels = c("WT","KO"),ordered = F)
## 表达矩阵矫正
exprSet <- tpms
boxplot(exprSet,outline=FALSE, notch=T,col=group_list, las=2)
library(limma) 
head(exprSet)
##               WT1       WT2       WT3        WT4        KO1        KO2
## Zzz3    42.124648 23.356837 34.310374  28.135094  54.895016  41.104711
## Zzef1   40.009245 28.081663 43.804760  42.773594  18.238978  12.348534
## Zyx    103.047985 83.065266 73.743177 103.073219 320.047602 190.790613
## Zyg11b  19.859384 12.025697 17.759368  14.845641  22.107926  15.621464
## Zxdc    20.371582 14.089519 17.435818  14.917163  19.211321  12.538796
## Zxdb     5.170714  5.482996  4.413419   3.496625   7.223751   3.526581
##               KO3        KO4
## Zzz3    47.765309  25.362704
## Zzef1   15.046395   9.756136
## Zyx    225.774133 192.170262
## Zyg11b  20.294522   8.836657
## Zxdc    16.861710   9.471667
## Zxdb     7.879218   2.424491
exprSet=normalizeBetweenArrays(exprSet)
boxplot(exprSet,outline=FALSE, notch=T,col=group_list, las=2)
range(exprSet)
## [1]     0.0 16468.1
## 判断数据是否需要转换；处理同芯片，20以内与否
exprSet <- log2(exprSet+1)
range(exprSet)
## [1]  0.00000 14.00747
## 差异分析：
dat <- exprSet
design=model.matrix(~factor( group_list ))
fit=lmFit(dat,design)
fit=eBayes(fit)
options(digits = 4)
topTable(fit,coef=2,adjust='BH')
##           logFC AveExpr      t   P.Value adj.P.Val     B
## Slc1a1  -11.170   5.682 -52.12 8.303e-12 1.198e-07 16.36
## Sox17     7.471   4.029  44.73 2.971e-11 2.144e-07 15.60
## Msx1      6.912   4.405  40.46 6.846e-11 3.293e-07 15.04
## Notum     8.651   5.182  35.94 1.836e-10 5.677e-07 14.32
## Asb4      6.560   3.802  35.14 2.212e-10 5.677e-07 14.18
## Bpifb5    5.543   2.817  34.87 2.360e-10 5.677e-07 14.13
## Sbspon    5.064   3.056  33.05 3.686e-10 7.598e-07 13.78
## Wif1      5.249   3.437  31.96 4.868e-10 8.151e-07 13.56
## Slc38a4   5.391   4.767  31.58 5.380e-10 8.151e-07 13.47
## Prr18     4.653   3.414  31.39 5.648e-10 8.151e-07 13.43
deg=topTable(fit,coef=2,adjust='BH',number = Inf)
head(deg)
##          logFC AveExpr      t   P.Value adj.P.Val     B
## Slc1a1 -11.170   5.682 -52.12 8.303e-12 1.198e-07 16.36
## Sox17    7.471   4.029  44.73 2.971e-11 2.144e-07 15.60
## Msx1     6.912   4.405  40.46 6.846e-11 3.293e-07 15.04
## Notum    8.651   5.182  35.94 1.836e-10 5.677e-07 14.32
## Asb4     6.560   3.802  35.14 2.212e-10 5.677e-07 14.18
## Bpifb5   5.543   2.817  34.87 2.360e-10 5.677e-07 14.13
fpkm_deg=deg

3.进行gseKEGG分析，绘制Wnt信号通路

## 按logFC排序，获得genelist,进行gseKEGG分析
DEG_DESeq2=fpkm_deg
head(DEG_DESeq2)
##          logFC AveExpr      t   P.Value adj.P.Val     B
## Slc1a1 -11.170   5.682 -52.12 8.303e-12 1.198e-07 16.36
## Sox17    7.471   4.029  44.73 2.971e-11 2.144e-07 15.60
## Msx1     6.912   4.405  40.46 6.846e-11 3.293e-07 15.04
## Notum    8.651   5.182  35.94 1.836e-10 5.677e-07 14.32
## Asb4     6.560   3.802  35.14 2.212e-10 5.677e-07 14.18
## Bpifb5   5.543   2.817  34.87 2.360e-10 5.677e-07 14.13
nrDEG=DEG_DESeq2[,c(1,4)]
colnames(nrDEG)=c('log2FoldChange','pvalue')
head(nrDEG) 
##        log2FoldChange    pvalue
## Slc1a1        -11.170 8.303e-12
## Sox17           7.471 2.971e-11
## Msx1            6.912 6.846e-11
## Notum           8.651 1.836e-10
## Asb4            6.560 2.212e-10
## Bpifb5          5.543 2.360e-10
library(org.Hs.eg.db)
library(org.Mm.eg.db)
library(clusterProfiler)
R.utils::setOption("clusterProfiler.download.method",'auto')
gene <- bitr(rownames(nrDEG), fromType = "SYMBOL",
             toType =  "ENTREZID",
             OrgDb = org.Mm.eg.db)
gene$logfc <- nrDEG$log2FoldChange[match(gene$SYMBOL,rownames(nrDEG))]
geneList=gene$logfc
names(geneList)=gene$ENTREZID 
geneList=sort(geneList,decreasing = T)
## 进行gseKEGG分析
kk_gse <- gseKEGG(geneList     = geneList,
                  organism     = 'mmu',
                  nPerm        = 1000,
                  minGSSize    = 10,
                  pvalueCutoff = 0.9,
                  verbose      = FALSE)
library(enrichplot)

## 5个通路中只找到Wnt信号通路---对应的通路ID号为mmu04310
paths="mmu04310"
gseaplot2(kk_gse, paths, color = colorspace::rainbow_hcl(4),
          pvalue_table = TRUE,title=("Wnt signaling pathway"))

4. 自定义通路，进行GSEA分析

文中作者根据以前发表的文献，自定义了一条肠干细胞通路；为文章的附表4，感兴趣的同学可以下载,尝试绘制下。

custome=readxl::read_xls("custome_gene_signature.xls")
head(custome)
## # A tibble: 6 x 4
##   `Table S4 : custom gene set signatures`                     ...2   ...3  ...4 
##   <chr>                                                       <chr>  <chr> <chr>
## 1 GSEA INTESTINAL STEM CELL SIGNATURE (Mu<U+00F1>oz, J. et al. 2012) <NA>   <NA>  <NA> 
## 2 PROBE                                                       DESCR~ GENE~ GENE~
## 3 Sycn                                                        NM_02~ SYCN  sync~
## 4 Fras1                                                       NM_17~ FRAS1 Fras~
## 5 Vdr                                                         NM_00~ VDR   vita~
## 6 Sfrp5                                                       NM_01~ SFRP5 secr~
geneset=custome$`Table S4 : custom gene set signatures`
geneset=geneset[2:length(geneset)]
# 建立一个符合gmt文件转换的list
ct_lis=list()
ct_lis[[1]]=geneset
names(ct_lis)="Interstinal stem cells"
# 获取gmt文件
file="gmt.txt"
gs=ct_lis
write.gmt <- function(gs,file){
  sink(file)
  lapply(names(ct_lis), function(i){
    cat( paste(c(i,'tmp',gs[[i]]),collapse='\t') )
    cat('\n')
  })
  sink()
}
write.gmt(gs,file)

names(geneList)=gene$SYMBOL 
geneset <- read.gmt('gmt.txt') #本质只有两列
egmt <- GSEA(geneList, TERM2GENE=geneset, 
             verbose=FALSE)
## 绘制gsea分析的结果
gseaplot2(egmt, color = colorspace::rainbow_hcl(4),
          geneSetID = rownames(egmt[1,]),title=("Interstinal stem cells"), pvalue_table = TRUE)

以上演示了批量读取单个样本的转录组数据，将其整合，进行差异分析以及绘制GSEA通路图。文中GSEA通路图有5副，但是由于小编在原文中只看见了两幅，并且两幅一定程度上代表了GSEA分析，gseKEGG与自定义通路，所以小编在此展现两幅。当然，可以导入其它的通路以及其它包进行分析，如Msigdb等。感兴趣的小伙伴可以阅读原文尝试下，若看见了其余三条通路，也可以一起沟通交流下。 

小编绘制的两幅通路图的结果与作者的结果是一致的，相较于WT组，两个通路在KO组均显著富集。