当前位置： 首页 > news >正文

可变剪接分析一步到位，这个 R 包够猛！

news 2024/11/15 9:06:22

生信碱移

ASpediaFI可变剪接

可变剪接（Alternative Splicing, AS）是基因表达过程中一种重要的调控机制，通过这种机制，单个基因可以产生多个不同的mRNA转录本，这些转录本通过不同的剪接方式（即选择性地包括或排除特定的外显子、内含子或外显子的一部分）生成。结果是，一个基因能够编码多种不同的蛋白质，从而极大地增加了蛋白质的多样性。

可变剪接的主要类型包括：

  • 外显子跳跃（Exon Skipping）：某些外显子被跳过，不被包含在成熟mRNA中。

  • 可变5'端剪接（Alternative 5' Splice Site）：在mRNA前端（5'）的剪接位点是不同的位置，导致外显子长度的变化。

  • 可变3'端剪接（Alternative 3' Splice Site）：在mRNA后端（3'）的剪接位点是不同的位置，导致外显子长度的变化。

  • 内含子保留（Intron Retention）：内含子没有被切除，而是保留在成熟mRNA中。

  • 互斥外显子（Mutually Exclusive Exons）：一对外显子在不同的转录本中互相排斥，只能选择其中一个。

图片

▲ 可变剪接的类型

在可变剪接的RNA-Seq分析中，rMATS和MAJIQ是最常用的两种工具，但是它们无法在windows系统上运行。此外，有数百种剪接因子调控AS事件，这对多种生物功能产生显著影响。然而，识别与剪接体相关的功能性AS事件以及探索与特定通路相对应的相互作用基因是一项挑战。为此，来自韩国国家癌症中心研究所的研究人员开发了一个R包ASpediaFI，于2022年1月25日发表于Genomics, Proteomics & Bioinformatics[IF: 11.5]，用于对可变剪接事件及其功能相互作用进行系统和综合的分析。

图片

▲ DOI: 10.1016/j.gpb.2021.10.004

ASpediaFI的工作流程基于R语言软件，其需要RNA-Seq生成的BAM文件和参考数据集，包括GTF文件、基因-基因交互网络和通路基因集。简单来讲，ASpediaFI通过获取比对结果定量可变剪接（AS）事件，随后构建一个异构网络，并使用带重启的判别随机游走（DRaWR）算法在该网络上对AS事件和通路进行排名，最终生成特定查询的子网络，并提供基因和特征的排名列表以便进一步分析。

图片

▲ ASpediaFI的分析工作流程。首先获取AS事件注释和定量。然后，使用AS事件和基因表达定量、基因-基因交互网络和通路基因集构建一个包含多种类型节点和边的异构网络。初始异构网络（如图B所示）由基因节点和两种类型的特征节点组成：AS事件和通路。下一步是在异构网络上运行DRaWR，以对AS事件和通路进行排名，评估它们与感兴趣的基因集（称为“查询”）的相关性。DRaWR算法在图B中展示：给定异构网络和一个查询基因集（用黄色标记），进行两阶段的带重启随机游走（RWR）。在第一阶段，RWR在初始网络上运行两次，一次使用查询基因集，另一次使用网络中的所有基因作为重启集。特征节点根据两次RWR中收敛的概率分布之间的差异进行排名。除了排名前k的节点外，所有特征节点都被移除，以重建一个特定于查询的子网络。在第二阶段，RWR在特定于查询的子网络上运行，以获得基因和特征的最终排名。对最终子网络的基因节点进行置换检验，以测试它们与查询基因集的关联。ASpediaFI为用户提供最终子网络以及基因和特征的排名列表，以便进行进一步分析，包括基因集富集分析和可视化。

本文简要展示ASpediaFI的使用，更多信息可以参见：

  • https://bioconductor.org/about/removed-packages/

0.安装

使用以下代码安装R包：

if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("ASpediaFI")

1.简要介绍

① ASpediaFI提供以下功能：

  1. AS事件检测与标注

  2. AS事件量化

  3. AS事件的功能交互分析

  4. AS事件和通路的可视化

② 需要注意的是，该软件包使用S4类ASpediaFI对象作为其功能（方法）的封装器，以及输入和输出（槽）的容器。

#Load the ASpediaFI package
library(ASpediaFI)
names(getSlots("ASpediaFI"))
#[1] "samples" "events" "psi" "gtf"
#[5] "network" "gene.table" "as.table" "pathway.table"

ASpediaFI类包含以下数据槽：

  • samples: 一个包含样本信息的数据框。前三列应为名称、BAM 文件路径和条件。

  • events: 从GTF文件中提取的AS事件列表。

  • psi: 一个包含AS事件定量的SummarizedExperiment对象。

  • gtf: 一个包含从GTF文件中提取的基因组特征的GRanges对象。

  • network: 一个包含特定查询子网络的igraph对象，作为DRaWR的结果。

  • gene.table, as.table, pathway.table: 包含基因节点、AS事件节点和通路节点的数据框。

③ ASpediaFI通过以下方法以逐步方式进行分析：

  1. annotateASevents: 从GTF文件中检测AS事件并将其保存在事件字段中。同时从GTF文件中提取特征并保存在gtf字段中。

  2. quantifyPSI: 使用样本字段中指定的BAM文件量化AS事件。

  3. analyzeFI: 构建基因、AS事件和通路的异构网络，并执行 DRaWR。

  4. visualize: 将事件或路径节点可视化。

  5. exportNetwork: 将与给定通路相关的子网络导出为GML格式，可以直接在Cytoscape中使用。

3.示例运行

① 首先，使用ASpediaFI函数创建ASpediaFI对象，该构造函数需要样本名称sample.names、RNA-Seq BAM文件的路径bam.files和样本分组conditions。然后，使用 annotateASevents方法从GRCh38 GTF文件中获取AS事件注释。由于文件大小限制，从包的extdata目录中提供的GRCh38 GTF文件的子集中提取AS事件注释。

#Create ASpediaFI object
bamWT <- system.file("extdata/GSM3167290.subset.bam", package = "ASpediaFI")
GSE114922.ASpediaFI <- ASpediaFI(sample.names = "GSM3167290",bam.files = bamWT, conditions = "WT")
#Detect and annotate AS events from a subset of the hg38 GTF file
gtf <- system.file("extdata/GRCh38.subset.gtf", package = "ASpediaFI")
GSE114922.ASpediaFI <- annotateASevents(GSE114922.ASpediaFI,gtf.file = gtf, num.cores = 1)
#[1] "-------------------Processing : chr11 -------------------"
sapply(events(GSE114922.ASpediaFI), length)
#A5SS A3SS SE MXE RI
#35 21 49 40 56head(events(GSE114922.ASpediaFI)$SE)
#EnsID Nchr Strand 1stEX DownEX
#1 "ENSG00000256269.10" "chr11" "+" "119089683-119089760" "119089217-119089272"
#2 "ENSG00000256269.10" "chr11" "+" "119089082-119089131" "119088635-119088848"
#3 "ENSG00000256269.10" "chr11" "+" "119089082-119089131" "119088635-119088707"
#4 "ENSG00000256269.10" "chr11" "+" "119089100-119089131" "119088635-119088707"
#5 "ENSG00000256269.10" "chr11" "+" "119092125-119092163" "119091413-119091526"
#6 "ENSG00000256269.10" "chr11" "+" "119092125-119092163" "119091861-119091874"
#UpEX
#1 "119089990-119090067"
#2 "119089217-119089272"
#3 "119089217-119089272"
#4 "119089217-119089272"
#5 "119092404-119092523"
#6 "119092404-119092523"
#EventID
#1 "HMBS:SE:chr11:119089217:119089272:119089683:119089760:119089990:119090067"
#2 "HMBS:SE:chr11:119088635:119088848:119089082:119089131:119089217:119089272"
#3 "HMBS:SE:chr11:119088635:119088707:119089082:119089131:119089217:119089272"
#4 "HMBS:SE:chr11:119088635:119088707:119089100:119089131:119089217:119089272"
#5 "HMBS:SE:chr11:119091413:119091526:119092125:119092163:119092404:119092523"
#6 "HMBS:SE:chr11:119091861:119091874:119092125:119092163:119092404:119092523"

annotateASevents方法识别五种AS事件：

  • A5SS (alternative 5’splice site)

  • A3SS (alternative 3’splice site)

  • SE (skipped exon)

  • MXE (mutually exclusive exon)

  • RI (retained intron)

AS事件注释列表包含Ensembl ID、染色体、链、外显子的基因组坐标和AS事件ID。AS事件ID的格式为[基因符号]:[事件类型]:[染色体]:[外显子边界的基因组坐标]，如ASpedia所定义（https://combio.snu.ac.kr/aspedia/help.html）。请注意，事件功能可用于访问AS事件注释。annotateASevents方法还从GTF文件中提取基因组特征，并将其保存为GRanges对象，以便可视化AS事件。

② 接下来，使用quantifyPSI方法从BAM文件中量化AS事件。quantifyPSI方法计算PSI（剪接包含百分比），即包含可变剪接外显子的mRNA的比例。用户需要指定RNA-Seq读取的类型（单端或双端）、读取长度、插入大小以及映射到给定外显子的最小读取数。在这一点上，从两个BAM文件的子集计算PSI值以进行演示。quantifyPSI方法将PSI值保存在psi槽中，作为带有样本信息的SummarizedExperiment对象。同样，可以使用psi函数访问PSI值。请注意，PSI值的行名称是AS事件ID（有了PSI，就可以对psi进行更多相关的分析来评估可变剪接的变化了）。

#Compute PSI values of AS events
GSE114922.ASpediaFI <- quantifyPSI(GSE114922.ASpediaFI, read.type = "paired",read.length = 100, insert.size = 300,min.reads = 3, num.cores = 1)
#[1] "Calculating PSI of SE events"
#[1] "Calculating PSI of MXE events"
#[1] "Calculating PSI of RI events"
#[1] "Calculating PSI of ALSS events"tail(assays(psi(GSE114922.ASpediaFI))[[1]])
#GSM3167290
#HMBS:RI:chr11:119088635:119088707:119089100:119089131 1.00
#HMBS:RI:chr11:119089683:119089760:119089990:119090067 0.49
#HMBS:RI:chr11:119092125:119092163:119092404:119092523 0.68
#HMBS:RI:chr11:119092125:119092163:119092758:119092811 1.00
#HMBS:RI:chr11:119092125:119092523:119092758:119092811 0.55
#HMBS:RI:chr11:119087987:119088078:119088255:119088308 0.96

由于接下来需要所有样本的PSI和基因表达谱来构建异质网络，因此我们在包中加载示例数据集。我们用示例数据集中存储的PSI值和样本信息更新psi和samples槽。

③ AS事件的功能交互分析：analyzeFI方法执行数据预处理、网络构建和DRaWR（带重启的判别随机游走）。由于DRaWR算法需要查询基因集作为输入，所以先使用limma包检测在SF3B1突变样本中差异表达的基因（也可以使用其他DEG分析工具，如edgeR和DESeq2）。假设DEG代表一个功能基因集，将其作为查询，以识别与SF3B1突变密切相关的AS事件和通路。如果查询作为字符向量给出，则查询中的所有基因具有相等的权重。当然，也可以通过提供一个包含第二列权重的数据框来赋予每个基因不同的权重(用logFC什么的)。

#Choose query genes based on differential expression
library(limma)design <- cbind(WT = 1, MvsW = samples(GSE114922.ASpediaFI)$condition == "MUT")
fit <- lmFit(log2(GSE114922.fpkm + 1), design = design)
fit <- eBayes(fit, trend = TRUE)
tt <- topTable(fit, number = Inf, coef = "MvsW")
query <- rownames(tt[tt$logFC > 1 &tt$P.Value < 0.1,])head(query)
#[1] "HBB" "FTL" "HBA2" "HBA1" "ATP6V0D2" "RPL37A"

接下来就运行analyzeFI方法。restart和num.feats定义了重启概率和最终子网络中要保留的特征数量。重启概率是跳回重启集（查询）的概率。如果重启概率较小，游走倾向于在查询节点的邻居之间移动。num.folds指定了DRaWR中交叉验证的折数。low.expr、low.var、prop.na和prop.extreme 是过滤 AS 事件的选项。cor.threshold定义了在异构网络中连接AS事件节点和基因节点的斯皮尔曼相关性的阈值。有关详细信息，可以参见help(analyzeFI)。

#Perform functional interaction analysis of AS events
GSE114922.ASpediaFI <- analyzeFI(GSE114922.ASpediaFI, query = query,expr = GSE114922.fpkm, restart = 0.9,num.folds = 5, num.feats = 200,low.expr = 1, low.var = 0, prop.na = 0.05,prop.extreme = 1, cor.threshold = 0.3)

图片

如下所示，analyzeFI方法将排名最高的AS事件和通路分别保存在as.table和pathway.table字段中。as.table包含有关AS事件的信息，包括AS事件ID、基因符号、AS事件类型、最终排名和静态概率。

#Table of AS nodes in the final subnetwork
as.table(GSE114922.ASpediaFI)[1:5,]
#EventID
#1 TIMM17B:SE:chrX:48896858:48896759:48895972:48895352:48895101:48895038
#2 NKTR:RI:chr3:42619664:42619708:42621429:42621516
#3 MCM7:RI:chr7:100099403:100099279:100099203:100099023
#4 COMMD3:RI:chr10:22318272:22318306:22318653:22318713
#5 MOV10:RI:chr1:112699894:112699982:112700219:112700340
#GeneSymbol EventType Rank StatP
#1 TIMM17B SE 1 0.00103
#2 NKTR RI 2 0.00099
#3 MCM7 RI 3 0.00099
#4 COMMD3 RI 4 0.00094
#5 MOV10 RI 5 0.00089#Table of GS nodes in the final subnetwork
pathway.table(GSE114922.ASpediaFI)[1:5,]
#Pathway Rank StatP Pvalue Adj.Pvalue EnrichmentScore 
#1 HALLMARK_HEME_METABOLISM 2 0.00089 1.4e-11 8.7e-09 0.55
#2 HALLMARK_COMPLEMENT 3 0.00066 1.0e-09 9.1e-08 0.53
#3 REACTOME_INNATE_IMMUNE_SYSTEM 4 0.00057 5.6e-03 3.5e-02 0.37
#4 HALLMARK_P53_PATHWAY 5 0.00054 1.6e-03 1.1e-02 0.41
#NormalizedEnrichmentScore Size Count AvgRank NumEvents
#1 1.4 933 618 1663 83
#2 2.0 200 133 1500 83
#3 1.9 200 119 1178 83
#4 1.4 279 153 1153 83

gs.table包含关于路径节点的以下信息:

  • Pathway: 通路名称

  • Rank: 最终排名

  • StatP: 平稳概率

  • Pvalue: GSEA P值

  • Adj.Pvalue: GSEA矫正后的P值

  • EnrichmentScore: 富集得分

  • NormalizedEnrichmentScore: 标准化的富集得分

  • Size: 通路基因集中的总基因数量

  • Count: 通路基因集中的基因数量，这些基因也存在于网络中

  • AvgRank: 通路基因集中的基因平均最终排名

  • NumEvents: 与通路基因集中的基因连接的最终子网络中的AS事件数量

  • Genes: 在通路基因集中存在的基因，同时也出现在网络中

analyzeFI方法将最终的查询特定子网络保存为igraph对象。用户可以通过以下方式探索AS事件与基因之间的相互作用：

#Extract AS-gene interactions from the final subnetwork
library(igraph)
edges <- as_data_frame(network(GSE114922.ASpediaFI))
AS.gene.interactions <- edges[edges$type == "AS", c("from", "to")]
head(AS.gene.interactions)
#from
#27032 A2M
#27033 A2M
#27034 A2M
#27035 A2M
#27036 A2M
#27037 A2M
#to
#27032 CLK4:SE:chr5:178623416:178623256:178620673:178619838:178618778:178618556
#27033 MRRF:A5SS:chr9:122285169:122285287:122285946:122286039:122286166
#27034 VPS29:SE:chr12:110499546:110499457:110497002:110496823:110496203:110496046
#27035 ERCC8:SE:chr5:60898400:60898276:60893643:60892037:60891086:60890889
#27036 CDK5RAP3:RI:chr17:47974400:47974448:47975159:47975337
#27037 RBM4B:SE:chr11:66668685:66668615:66666493:66666272:66665578:66665548

ASpediaFI还允许用户使用exportNetwork方法导出整个子网络或与特定通路相关的子网络。给定一个通路节点，提取特定于通路的网络并将其导出为GML格式，可以直接在Cytoscape中使用。如果未给定通路节点，则导出整个最终子网络。

#Export a pathway-specific subnetwork to GML format
exportNetwork(GSE114922.ASpediaFI, node = "HALLMARK_HEME_METABOLISM",
file = "heme_metabolism.gml")

④ 可视化，首先可以生成一个描述指定AS事件的基因组轨迹图和一个PSI值的箱线图

#Check if any event on the HMBS gene is included in the final subnetwork
as.nodes <- as.table(GSE114922.ASpediaFI)$EventID
HMBS.event <- as.nodes[grep("HMBS", as.nodes)]
#Visualize event
visualize(GSE114922.ASpediaFI, node = HMBS.event, zoom = FALSE)

图片

显示一个由高度排名的基因节点和与给定通路连接的AS事件节点组成的子网络，比如示例数据中血红素代谢的标志通路相关的子网络，该通路也与 MDS 中的SF3B1突变相关：

#Visualize nework pertaining to specific pathway
visualize(GSE114922.ASpediaFI, node = "HALLMARK_HEME_METABOLISM", n = 10)

图片

功能是不少的

可以整起

欢迎各位老哥老姐点击关注

就分享到这


http://www.mrgr.cn/news/35803.html

相关文章：

  • Tofu AI视频处理模块视频输入配置方法
  • 国标GB28181视频平台EasyCVR私有化部署视频平台对接监控录像机NVR时，录像机“资源不足”是什么原因？
  • WebGIS三维地图框架--Cesium
  • 丹摩征文活动｜FLUX.1 和 ComfyUI：从部署到上手，轻松驾驭！
  • kaggle 如何利用API下载数据集
  • C++组合复用中，委托的含义与作用
  • 表观项目文章速递，平均IF=9.7
  • 【tbNick专享】虚拟机域控、成员服务器、降级等管理
  • oracle 分表代码示例
  • Jetpack02-LiveData 数据驱动UI更新（类似EventBus）
  • 《AI设计类工具系列之三——Magic Design》
  • 7大网络安全机构真实测评：零基础转行网安这么选就对了~
  • 【记录】在返回值类型为BigDecimal情况下末尾小数位为0的会省略不显示
  • 邦芒干货：职场上快速适应工作新环境的十大举措
  • 表单多个输入框校验
  • 外网SpringBoot项目迁移到内网无法加载maven依赖
  • 逆概率加权（R和Python案例）
  • implementation ‘com.lxj:xpopup:2.2.1‘失败
  • 【Linux】驱动的基本架构和编译
  • TON生态系统开发指南：从零开始构建你的Web3应用
  • 手撕Transformer之Embedding Layer
  • 基于pdf.js实现对pdf预览、批注功能、 保存下载pdf,适配H5,平板 踩坑记录
  • 【鸿蒙】HarmonyOS NEXT开发快速入门教程之ArkTS语法装饰器（上）
  • ntpdate同步配置文件调整详解
  • 聊一聊软件系统性能测试的重要性
  • dspic33F Mplab IDE V8.92内存使用情况