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语言函数递归实际应用举例详解

《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

找不到Anaconda prompt终端的原因分析及解决方案

《找不到Anacondaprompt终端的原因分析及解决方案》因为anaconda还没有初始化,在安装anaconda的过程中,有一行是否要添加anaconda到菜单目录中,由于没有勾选,导致没有菜... 目录问题原因问http://www.chinasem.cn题解决安装了 Anaconda 却找不到 An

Spring定时任务只执行一次的原因分析与解决方案

《Spring定时任务只执行一次的原因分析与解决方案》在使用Spring的@Scheduled定时任务时,你是否遇到过任务只执行一次,后续不再触发的情况?这种情况可能由多种原因导致,如未启用调度、线程... 目录1. 问题背景2. Spring定时任务的基本用法3. 为什么定时任务只执行一次?3.1 未启用

C语言中的数据类型强制转换

《C语言中的数据类型强制转换》:本文主要介绍C语言中的数据类型强制转换方式,具有很好的参考价值,希望对大家有所帮助,如有错误或未考虑完全的地方,望不吝赐教... 目录C语言数据类型强制转换自动转换强制转换类型总结C语言数据类型强制转换强制类型转换:是通过类型转换运算来实现的,主要的数据类型转换分为自动转换

利用Go语言开发文件操作工具轻松处理所有文件

《利用Go语言开发文件操作工具轻松处理所有文件》在后端开发中,文件操作是一个非常常见但又容易出错的场景,本文小编要向大家介绍一个强大的Go语言文件操作工具库,它能帮你轻松处理各种文件操作场景... 目录为什么需要这个工具?核心功能详解1. 文件/目录存javascript在性检查2. 批量创建目录3. 文件

C语言实现两个变量值交换的三种方式

《C语言实现两个变量值交换的三种方式》两个变量值的交换是编程中最常见的问题之一,以下将介绍三种变量的交换方式,其中第一种方式是最常用也是最实用的,后两种方式一般只在特殊限制下使用,需要的朋友可以参考下... 目录1.使用临时变量(推荐)2.相加和相减的方式(值较大时可能丢失数据)3.按位异或运算1.使用临时

使用C语言实现交换整数的奇数位和偶数位

《使用C语言实现交换整数的奇数位和偶数位》在C语言中,要交换一个整数的二进制位中的奇数位和偶数位,重点需要理解位操作,当我们谈论二进制位的奇数位和偶数位时,我们是指从右到左数的位置,本文给大家介绍了使... 目录一、问题描述二、解决思路三、函数实现四、宏实现五、总结一、问题描述使用C语言代码实现:将一个整