R语言:GSEA分析

2024-05-13 20:04
文章标签 语言 分析 gsea

本文主要是介绍R语言:GSEA分析,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

#安装软件包

> if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
> BiocManager::install("limma")
> BiocManager::install("org.Hs.eg.db")
> BiocManager::install("DOSE")
> BiocManager::install("clusterProfiler")
> BiocManager::install("enrichplot")

#加载软件包

> library(limma)
> library(org.Hs.eg.db)
> library(clusterProfiler)
> library(enrichplot)

#设置变量

> gene="CRTAC1"

> expFile="combined_RNAseq_counts.txt"

> gmtFile="c5.go.v7.4.symbols.gmt"

> rt=read.table(expFile, header=T, sep="\t", check.names=F)
> head(rt)

head(rt)id TCGA-E2-A1L7-11A TCGA-E2-A1L7-01A TCGA-AR-A0U0-01A TCGA-BH-A28O-01A TCGA-A2-A0D4-01A TCGA-E9-A1R4-01A TCGA-AO-A1KQ-01ATCGA-AC-A62V-01A TCGA-D8-A143-01A TCGA-A2-A0SV-01A TCGA-AN-A0XW-01A TCGA-D8-A1XV-01A TCGA-A2-A4RW-01A TCGA-A7-A0CD-01A TCGA-E2-A1IG-11ATCGA-D8-A1XB-01A TCGA-C8-A134-01A TCGA-BH-A0BS-11A TCGA-AR-A2LE-01A TCGA-A2-A0CO-01A TCGA-E9-A1NA-11A TCGA-AN-A0AK-01A TCGA-E9-A1NA-01ATCGA-A7-A0DA-01A TCGA-E2-A572-01A TCGA-A2-A259-01A TCGA-BH-A28Q-01A TCGA-E2-A1IO-01A TCGA-AQ-A7U7-01A TCGA-AN-A0FD-01A TCGA-A8-A07G-01ATCGA-AO-A0JL-01A TCGA-B6-A0IM-01A TCGA-B6-A0IP-01A TCGA-GM-A2DF-01A TCGA-A2-A25B-01A TCGA-BH-A0B0-01A TCGA-AO-A0JD-01A TCGA-AN-A0FL-01ATCGA-E2-A14V-01A TCGA-AN-A0FF-01A TCGA-C8-A138-01A TCGA-E2-A14R-01A TCGA-AC-A2BM-01A TCGA-A1-A0SP-01A TCGA-A2-A0CQ-01A TCGA-A8-A08J-01ATCGA-BH-A6R8-01A TCGA-E9-A1QZ-01A TCGA-A8-A0AB-01A TCGA-BH-A0H9-11A TCGA-AC-A3W7-01A TCGA-B6-A0IE-01A TCGA-A8-A07I-01A TCGA-BH-A0BQ-11A

> rt=as.matrix(rt)
> rownames(rt)=rt[,1]
> exp=rt[,2:ncol(rt)]
> dimnames=list(rownames(exp),colnames(exp))
> data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
> data=avereps(data)
> data=data[rowMeans(data)>0,]

> group=sapply(strsplit(colnames(data),"\\-"), "[", 4)
> group=sapply(strsplit(group,""), "[", 1)
> group=gsub("2", "1", group)
> data=data[,group==0]
> data=t(data)
> rownames(data)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-.*", "\\1\\-\\2\\-\\3", rownames(data))
> data=t(avereps(data))

> dataL=data[,data[gene,]<=median(data[gene,]),drop=F]
> dataH=data[,data[gene,]>median(data[gene,]),drop=F]
> meanL=rowMeans(dataL)
> meanH=rowMeans(dataH)
> meanL[meanL<0.00001]=0.00001
> meanH[meanH<0.00001]=0.00001
> logFC=log2(meanH)-log2(meanL)

#排序
> logFC=sort(logFC,decreasing=T)
> genes=names(logFC)

> gmt=read.gmt(gmtFile)

#GESA分析
> kk=GSEA(logFC, TERM2GENE=gmt, pvalueCutoff = 1)


> kkTab=as.data.frame(kk)
> kkTab=kkTab[kkTab$pvalue<0.05,]
> write.table(kkTab,file="GSEA.result-GO.txt",sep="\t",quote=F,row.names = F)
    

> termNum=5   
> if(nrow(kkTab)>=termNum){
    showTerm=row.names(kkTab)[1:termNum]
    gseaplot=gseaplot2(kk, showTerm, base_size=8, title="")
    pdf(file="GSEA-GO.pdf", width=10, height=8)
    print(gseaplot)
    dev.off()
}

> my=read.table("my.txt", header=T, sep="\t", check.names=F)
> my=as.matrix(my)
> rownames(my)=my[,1]
> mys=my[,2:ncol(my)]
> showmy=row.names(mys)
> myplot=gseaplot2(kk, showmy, base_size=8, title="")
> pdf(file="GSEA-GO-myself.pdf", width=10, height=8)
> print(myplot)
> dev.off()

> gmtFile="c2.cp.kegg.v7.4.symbols.gmt"     
> gmt=read.gmt(gmtFile)
> kk=GSEA(logFC, TERM2GENE=gmt, pvalueCutoff = 1)
> kkTab=as.data.frame(kk)
> kkTab=kkTab[kkTab$pvalue<0.05,]
> write.table(kkTab,file="GSEA.result-KEGG.txt",sep="\t",quote=F,row.names = F)


> termNum=5    
> if(nrow(kkTab)>=termNum){
  showTerm=row.names(kkTab)[1:termNum]
  gseaplot=gseaplot2(kk, showTerm, base_size=8, title="")
  pdf(file="GSEA-KEGG.pdf", width=10, height=8)
  print(gseaplot)
  dev.off()
}

一起学习交流。

这篇关于R语言:GSEA分析的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/986685

相关文章

C 语言中enum枚举的定义和使用小结

《C语言中enum枚举的定义和使用小结》在C语言里,enum(枚举)是一种用户自定义的数据类型,它能够让你创建一组具名的整数常量,下面我会从定义、使用、特性等方面详细介绍enum,感兴趣的朋友一起看... 目录1、引言2、基本定义3、定义枚举变量4、自定义枚举常量的值5、枚举与switch语句结合使用6、枚

Java程序进程起来了但是不打印日志的原因分析

《Java程序进程起来了但是不打印日志的原因分析》:本文主要介绍Java程序进程起来了但是不打印日志的原因分析,具有很好的参考价值,希望对大家有所帮助,如有错误或未考虑完全的地方,望不吝赐教... 目录Java程序进程起来了但是不打印日志的原因1、日志配置问题2、日志文件权限问题3、日志文件路径问题4、程序

Java字符串操作技巧之语法、示例与应用场景分析

《Java字符串操作技巧之语法、示例与应用场景分析》在Java算法题和日常开发中,字符串处理是必备的核心技能,本文全面梳理Java中字符串的常用操作语法,结合代码示例、应用场景和避坑指南,可快速掌握字... 目录引言1. 基础操作1.1 创建字符串1.2 获取长度1.3 访问字符2. 字符串处理2.1 子字

Python 迭代器和生成器概念及场景分析

《Python迭代器和生成器概念及场景分析》yield是Python中实现惰性计算和协程的核心工具,结合send()、throw()、close()等方法,能够构建高效、灵活的数据流和控制流模型,这... 目录迭代器的介绍自定义迭代器省略的迭代器生产器的介绍yield的普通用法yield的高级用法yidle

Go 语言中的select语句详解及工作原理

《Go语言中的select语句详解及工作原理》在Go语言中,select语句是用于处理多个通道(channel)操作的一种控制结构,它类似于switch语句,本文给大家介绍Go语言中的select语... 目录Go 语言中的 select 是做什么的基本功能语法工作原理示例示例 1:监听多个通道示例 2:带

C++ Sort函数使用场景分析

《C++Sort函数使用场景分析》sort函数是algorithm库下的一个函数,sort函数是不稳定的,即大小相同的元素在排序后相对顺序可能发生改变,如果某些场景需要保持相同元素间的相对顺序,可使... 目录C++ Sort函数详解一、sort函数调用的两种方式二、sort函数使用场景三、sort函数排序

C语言函数递归实际应用举例详解

《C语言函数递归实际应用举例详解》程序调用自身的编程技巧称为递归,递归做为一种算法在程序设计语言中广泛应用,:本文主要介绍C语言函数递归实际应用举例的相关资料,文中通过代码介绍的非常详细,需要的朋... 目录前言一、递归的概念与思想二、递归的限制条件 三、递归的实际应用举例(一)求 n 的阶乘(二)顺序打印

kotlin中const 和val的区别及使用场景分析

《kotlin中const和val的区别及使用场景分析》在Kotlin中,const和val都是用来声明常量的,但它们的使用场景和功能有所不同,下面给大家介绍kotlin中const和val的区别,... 目录kotlin中const 和val的区别1. val:2. const:二 代码示例1 Java

Go标准库常见错误分析和解决办法

《Go标准库常见错误分析和解决办法》Go语言的标准库为开发者提供了丰富且高效的工具,涵盖了从网络编程到文件操作等各个方面,然而,标准库虽好,使用不当却可能适得其反,正所谓工欲善其事,必先利其器,本文将... 目录1. 使用了错误的time.Duration2. time.After导致的内存泄漏3. jsO

Spring事务中@Transactional注解不生效的原因分析与解决

《Spring事务中@Transactional注解不生效的原因分析与解决》在Spring框架中,@Transactional注解是管理数据库事务的核心方式,本文将深入分析事务自调用的底层原理,解释为... 目录1. 引言2. 事务自调用问题重现2.1 示例代码2.2 问题现象3. 为什么事务自调用会失效3