从头合成 利用核糖体谱分析数据鉴定活跃翻译的开放阅读框

摘要

核糖体在翻译过程中以每密码子三个核苷酸的规律解码生成肽链。核糖体沿mRNA移动的过程通过核糖体谱型分析（ribosome profiling）得以捕获，产生具有典型三联体周期性的足迹信号。本实验方案介绍如何利用RiboCode从核糖体谱型数据中解析这一显著特征，从而在全转录组水平上鉴定处于活跃翻译状态的开放阅读框。

摘要

鉴定开放阅读框（ORF），特别是那些编码小肽且在特定生理条件下处于活跃翻译状态的ORF，对于全面注释依赖于特定条件的翻译组至关重要。 核糖体谱分析（ribosome profiling）是一种用于检测核糖体在RNA上的结合位置和密度的技术，可实现全基因组范围内翻译事件发生位置的快速发现。然而，在生物信息学分析中，如何高效且全面地从核糖体谱数据中识别正在翻译的ORF仍非易事。本文介绍了一个易于使用的分析工具包RiboCode，该工具包旨在从核糖体谱数据中扭曲且模糊的信号里搜索任意大小的活跃翻译ORF。以我们先前发表的数据集为例，本文逐步演示了RiboCode完整分析流程的操作步骤，涵盖从原始数据预处理到最终输出结果文件的解读。此外，为评估已注释ORF的翻译水平，本文还详细描述了对每个ORF上核糖体密度进行可视化和定量分析的具体方法。总之，本文为翻译、小ORF及小肽相关研究领域提供了实用且及时的技术指南。

引言

近年来，越来越多的研究揭示了由编码基因的开放阅读框（ORFs）以及先前被注释为非编码的基因（如长链非编码RNA（lncRNAs））广泛翻译产生多肽的现象^{1,2,3,4,5,6,7,8}。这些被翻译的ORF受到细胞调控或诱导，以响应环境变化、应激以及细胞分化过程^{1,8,9,10,11,12,13}。一些ORF的翻译产物已被证实可在发育和生理过程中参与多种生物学过程的调控。例如，Chng等人¹⁴发现了一种名为Elabela（Ela，亦称为Apela/Ende/Toddler）的肽类激素，其对心血管系统的发育至关重要。Pauli等人提出，Ela还具有促有丝分裂原功能，可促进早期鱼类胚胎中的细胞迁移¹⁵。Magny等人报道了两种长度小于30个氨基酸的微肽，它们可调控果蝇心脏中的钙离子转运并影响正常的肌肉收缩¹⁰。

目前尚不清楚基因组编码了多少此类肽段，以及它们是否具有生物学意义。因此，系统性地鉴定这些潜在的编码开放阅读框（ORF）显得尤为重要。然而，使用传统的研究方法（如进化保守性^16,17和质谱分析^18,19）直接确定这些ORF的产物（即蛋白质或肽段）仍面临挑战，因为这两种方法的检测效率均依赖于所产生蛋白质或肽段的长度、丰度及氨基酸组成。核糖体谱分析（ribosome profiling）技术的出现为在核苷酸分辨率上识别mRNA上的核糖体占据情况提供了精确手段，从而能够不受转录本长度和组成限制地评估其编码潜力^3,20,21。利用核糖体谱分析鉴定活跃翻译ORF的一个重要且常用特征，是核糖体在mRNA上从起始密码子到终止密码子之间所形成的三核苷酸（3-nt）周期性足迹模式。然而，核糖体谱分析数据常存在若干问题，包括沿ORF分布的测序读段稀少且覆盖度低、测序噪声较高以及核糖体RNA（rRNA）污染等。因此，此类数据产生的畸变和模糊信号会削弱核糖体在mRNA上足迹的3-nt周期性模式，最终导致高置信度翻译ORF的鉴定变得困难。

名为 "RiboCode" 采用了一种改进的 Wilcoxon 符号秩检验与 P 值整合策略，以检验开放阅读框（ORF）中的同相核糖体保护片段（RPFs）是否显著多于异相 RPFs²²已证明该方法高效、灵敏且准确。 从头合成 在模拟和真实核糖体谱型数据中对翻译组进行注释。本文介绍了如何利用该工具从先前研究产生的原始核糖体谱型测序数据集中检测潜在的翻译开放阅读框²³这些数据集已被用于探索EIF3亚基的功能 "E" 通过比较转染了对照（si-Ctrl）和 EIF3E 的 MCF-10A 细胞的核糖体占据谱，研究 EIF3E 在翻译过程中的作用 EIF3E (si-eIF3e) 小干扰RNA（siRNA）。通过将RiboCode应用于这些示例数据集，我们检测到5,633个可能编码小肽或蛋白质的新型开放阅读框（ORF）。这些ORF根据其相对于编码区的位置被分为多种类型，包括上游开放阅读框（uORF）、下游开放阅读框（dORF）、重叠开放阅读框、来自新型蛋白编码基因（novel PCGs）的开放阅读框以及来自新型非蛋白编码基因（novel NonPCGs）的开放阅读框。在EIF3E缺陷细胞中，uORF上的核糖体保护片段（RPF）读数密度相较于对照细胞显著增加，这可能至少部分归因于活跃翻译核糖体的富集。核糖体在25^th 至 75^th EIF3E 缺陷细胞中的密码子表明翻译延伸在早期阶段受阻。本方案还展示了如何可视化目标区域的 RPF 密度，以检测核糖体足迹在已鉴定 ORF 上的 3-nt 周期性模式。这些分析展示了 RiboCode 在鉴定翻译 ORF 及研究翻译调控方面的强大作用。

方案

1. 环境配置与 RiboCode 安装

1. 打开一个 Linux 终端窗口，并创建一个 conda 环境：
   conda create -n RiboCode python=3.8
2. 切换到所创建的环境，并安装 RiboCode 及其依赖项：
   conda activate RiboCode
   conda install -c bioconda ribocode ribominer sra-tools fastx_toolkit cutadapt bowtie star samtools

2. 数据准备

1. 获取基因组参考文件。
   1. 对于参考序列，请访问 Ensemble 网站 https://www.ensembl.org/index.html, 点击顶部菜单 下载 和左侧菜单 FTP 下载在所示表格中，单击 FASTA 在柱中 DNA（FASTA） 以及所在的行 物种 是 人类在打开的页面中，复制该链接 Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz，然后在终端中下载并解压：
      wget -c \
      http://ftp.ensembl.org/pub/release-104/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
      gzip -d Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
   2. 用于参考注释，请右键单击 GTF 在柱中 基因集 在最后打开的网页中。复制该链接 Homo_sapiens.GRCh38.104.gtf.gz 并使用以下方法下载：
      wget -c \
      http://ftp.ensembl.org/pub/release-104/gtf/homo_sapiens/Homo_sapiens.GRCh38.104.gtf.gz
      gzip -d Homo_sapiens.GRCh38.104.gtf.gz
      注意：建议从Ensemble网站获取GTF文件，因为它包含以三级层次结构组织的基因组注释，即每个基因包含转录本，而转录本又包含外显子以及可选的翻译序列（例如., 编码序列 [CDS]、翻译起始位点、翻译终止位点）。当基因或转录本的注释信息缺失时，例如从 UCSC 或 NCBI 获取的 GTF 文件，应使用 GTFupdate 生成一个具有完整父子层级注释的更新版GTF： GTFupdate original.gtf > 更新的.gtf文件对于 .gff 格式的注释文件，使用 AGAT 工具包²⁴ 或使用其他工具转换为 .gtf 格式。
2. 获取rRNA序列。
   1. 打开 UCSC 基因组浏览器 https://genome.ucsc.edu 并单击 工具 | 表格浏览器 在下拉列表中。
   2. 在打开的页面中，指定 哺乳动物 对于分支群， 人类 基因组 所有表格 对于组， rmask 用于表格，以及 基因组 用于区域。如需筛选，请点击 创建 跳转到新页面并将 repClass 设置为 与rRNA匹配.
   3. 点击 提交 然后将输出格式设置为 序列 以及输出文件名作为 hg38_rRNA.fa最后，单击 获取输出 | 获取序列 以获取rRNA序列。
3. 从序列读长存档（Sequence Read Archive，SRA）获取核糖体谱分析数据集。
   1. 下载si-eIF3e处理组的重复样本，并将其重命名：
      fastq-dump SRR9047190 SRR9047191 SRR9047192
      mv SRR9047190.fastq si-eIF3e-1.fastq
      mv SRR9047191.fastq si-eIF3e-2.fastq
      mv SRR9047192.fastq si-eIF3e-3.fastq
   2. 下载对照组的重复样本并重命名：
      fastq-dump SRR9047193 SRR9047194 SRR9047195
      mv SRR9047193.fastq si-Ctrl-1.fastq
      mv SRR9047194.fastq si-Ctrl-2.fastq
      mv SRR9047195.fastq si-Ctrl-3.fastq
      注意：这些示例数据集的 SRA 登录号来自基因表达综合数据库（Gene Expression Omnibus，GEO）网站²⁵ 通过搜索GSE131074获取。

3. 切除接头并去除rRNA污染

1. （可选）从测序数据中去除接头序列。如果接头序列已被修剪（如此处所示），则跳过此步骤。否则，使用 cutadapt 从读段中修剪接头。
   for i in si-Ctrl-1 si-Ctrl-2 si-Ctrl-3 si-eIF3e-1 si-eIF3e-2 si-eIF3e-3
   do
   cutadapt -m 15 --match-read-wildcards -a CTGTAGGCACCATCAAT \
   -o ${i}_trimmed.fastq ${i}.fastq
   done
   注意：参数 -a 后的接头序列会因 cDNA 文库制备方法的不同而有所变化。长度小于 15（由 -m 指定）的读段将被丢弃，因为核糖体保护片段通常长于该长度。
2. 使用以下步骤去除 rRNA 污染：
   1. 构建 rRNA 参考序列的索引：
      bowtie-build -f hg38_rRNA.fa hg38_rRNA
   2. 将读段比对至 rRNA 参考序列，以排除来源于 rRNA 的读段：
      for i in si-Ctrl-1 si-Ctrl-2 si-Ctrl-3 si-eIF3e-1 si-eIF3e-2 si-eIF3e-3
      do
      bowtie -n 0 -y -a --norc --best --strata -S -p 4 -l 15 \
      --un=./${i}_noncontam.fastq hg38_rRNA -q ${i}.fastq ${i}.aln
      done
      -p 指定并行运行任务所使用的线程数。考虑到 RPF 读段相对较小，应设置其他参数（如 -n、-y、-a、--norc、--best、--strata 和 -l），以确保报告的比对结果为最优比对。更多详细信息，请参阅 Bowtie 官方网站²⁶。

4. 将干净的测序读段比对到基因组

1. 创建基因组索引。
   mkdir STAR_hg38_genome
   STAR --runThreadN 8 --runMode genomeGenerate --genomeDir ./STAR_hg38_genome --genomeFastaFiles Homo_sapiens.GRCh38.dna.primary_assembly.fa --sjdbGTFfile Homo_sapiens.GRCh38.104.gtf
2. 将干净的测序数据（无rRNA污染）比对至构建的参考序列。
   对于 si-Ctrl-1、si-Ctrl-2、si-Ctrl-3、si-eIF3e-1、si-eIF3e-2、si-eIF3e-3
   执行
   STAR --runThreadN 8 --outFilterType 正常 --outWigType wiggle --outWigStrand 链特异性 --outWigNorm RPM --outFilterMismatchNmax 1 --outFilterMultimapNmax 1 --genomeDir STAR_hg38_genome --readFilesIn ${i}_非污染.fastq --outFileNamePrefix ${i}. --outSAMtype BAM 坐标排序 --quantMode 转录组SAM 基因计数 --outSAMattributes 全部
   完成
   注意：逆转录酶常在每条读段的5'端添加一个非模板依赖的核苷酸²⁷，STAR 默认会通过软剪切（soft-clipping）高效地将其切除。STAR 的参数说明详见 STAR 用户手册²⁸.
3. 排序并索引比对文件。
   对于 si-Ctrl-1、si-Ctrl-2、si-Ctrl-3、si-eIF3e-1、si-eIF3e-2、si-eIF3e-3 中的每一个
   做
   samtools sort -T ${i}.Aligned.toTranscriptome.out.sorted
   -o ${i}.Aligned.toTranscriptome.out.sorted.bam \
   ${i}.与转录组比对后的输出BAM文件
   samtools index ${i}.Aligned.toTranscriptome.out.sorted.bam
   samtools index ${i}.Aligned.sortedByCoord.out.bam
   完成

5. RPF 片段的大小筛选及其 P 位点的鉴定

1. 准备转录本注释。
   prepare_transcripts -g Homo_sapiens.GRCh38.104.gtf \
   -f Homo_sapiens.GRCh38.dna.primary_assembly.fa -o RiboCode_annot
   注意：该命令从 GTF 文件中收集 mRNA 转录本所需的注释信息，并从 FASTA 文件中提取所有 mRNA 转录本的序列（每个转录本根据 GTF 文件中定义的结构通过合并外显子组装而成）。
2. 选择特定长度的 RPFs 并确定其 P 位点位置。
   for i in si-Ctrl-1 si-Ctrl-2 si-Ctrl-3 si-eIF3e-1 si-eIF3e-2 si-eIF3e-3
   do
   metaplots -a RiboCode_annot -r ${i}.Aligned.toTranscriptome.out.bam \
   -o ${i} -f0_percent 0.35 -pv1 0.001 -pv2 0.001
   done
   注意：该命令绘制每个长度的测序读段在已注释翻译起始（或终止）密码子周围的 5' 端的总体分布图。通过检查主要读段 5' 端与起始密码子之间偏移距离的分布图（例如 图 1B），可手动确定依赖读段长度的 P 位点。RiboCode 还为每个样本生成一个配置文件，其中自动确定了具有显著三核苷酸周期性模式的读段的 P 位点位置。参数 -f0_percent、-pv1 和 -pv2 定义了选择在阅读框中富集的 RPF 读段的比例阈值和 p 值截断值。在本示例中，对于 29、30 和 31 nt 的读段，其 5' 端的 +12、+13 和 +13 核苷酸位置在每个配置文件中被手动定义。
3. 编辑每个样本的配置文件并将其合并
   注意：为生成一组一致的唯一开放阅读框（ORF），并确保后续分析有足够的读段覆盖度，需将上一步中所有样本的选定读段进行合并。在下一步评估 ORF 翻译潜能时，将使用 merged_config.txt 文件（补充文件 1）中定义的特定长度读段及其 P 位点信息。

6. 从头合成 注释翻译中的开放阅读框

1. 运行 RiboCode。
   RiboCode -a RiboCode_annot -c merged_config.txt -l yes -g \
   -o RiboCode_ORFs_result -s ATG -m 5 -A CTG,GTG,TTG
   该命令的重要参数如下：
   -c，包含输入文件路径以及所选读段及其P位点信息的配置文件。
   -l，对于在终止密码子上游具有多个起始密码子的转录本，是否使用最长的开放阅读框（从最远端的起始密码子到终止密码子的区域）来评估其翻译潜能。若设置为 no，起始密码子将被自动确定。
   -s，用于开放阅读框识别的标准起始密码子。
   -A，（可选）用于开放阅读框识别的非标准起始密码子（例如，人类中的CTG、GTG和TTG），在其他物种的线粒体或细胞核中可能有所不同²⁹。
   -m，开放阅读框的最小长度（即氨基酸数目）。
   -o，输出文件名的前缀，包含预测的开放阅读框的详细信息（补充文件 2）。
   -g 和 -b，分别将预测的开放阅读框输出为 gtf 或 bed 格式。

7. （可选）ORF 定量与统计分析

1. 统计每个 ORF 中的 RPF 读段数量。
   for i in si-Ctrl-1 si-Ctrl-2 si-Ctrl-3 si-eIF3e-1 si-eIF3e-2 si-eIF3e-3
   do
   ORFcount -g RiboCode_ORFs_result_collapsed.gtf \
   -r ${i}.Aligned.sortedByCoord.out.bam -f 15 -l 5 -m 25 -M 35 \
   -o ${i}_ORF.counts -s yes -c intersection-strict
   done
   注意：为排除核糖体在 ORF 起始端和终止端附近堆积的潜在影响，不计入前 15 个（由 -f 指定）和最后 5 个密码子（由 -l 指定）内的读段数量。可选地，可将用于计数的 RPF 长度限制在 25 至 35 nt 范围内（为 RPF 的常见大小）。
2. 使用 RiboCode 计算所检测 ORF 的基本统计信息：
   Rscript RiboCode_utils.R
   注意：RiboCode_utils.R（补充文件 3）提供了一系列针对 RiboCode 输出结果的统计功能，例如统计已鉴定 ORF 的数量、查看 ORF 长度分布，以及计算标准化的 RPF 密度（即 RPKM，每千碱基每百万比对读段的读段数）。

8. （可选）预测 ORF 的可视化

1. 从 RiboCode_ORFs_result_collapsed.txt（补充文件 3）中获取目标 ORF（例如，ENSG00000100902_35292349_35292552_67）在其转录本上的起始密码子与终止密码子的相对位置。然后，绘制该 ORF 的 RPF 读段密度图：
   plot_orf_density -a RiboCode_annot -c merged_config.txt -t ENST00000622405 \
   -s 33 -e 236 --start-codon ATG -o ENSG00000100902_35292349_35292552_67
   其中，-s 和 -e 指定待绘制 ORF 的翻译起始与终止位置。--start-codon 定义 ORF 的起始密码子，该密码子将显示在图的标题中。-o 定义输出文件名的前缀。

9. （可选）使用 RiboMiner 进行基因组平均分析

注意：按照以下步骤进行基因集分析，以评估EIF3E基因敲低对已鉴定注释ORF翻译的影响：

1. 为 RiboMiner 生成转录本注释，该步骤基于 RiboCode 生成的注释文件（步骤 5.1）提取每个基因的最长转录本。
   OutputTranscriptInfo -c RiboCode_annot/transcripts_cds.txt \
   -g Homo_sapiens.GRCh38.104.gtf -f RiboCode_annot/transcripts_sequence.fa \
   -o 最长转录本信息.txt -O 所有转录本信息.txt
2. 准备 RiboMiner 的配置文件。复制由程序生成的配置文件 元图 RiboCode 命令（步骤 5.4）并重命名它 "RiboMiner_config.txt" 然后，根据所示格式对其进行修改 补充文件 4.
3. 使用 RiboMiner 进行元基因分析
   1. 使用 MetageneAnalysis 用于 生成转录本上RPFs密度的聚合与平均谱型。
      MetageneAnalysis -f RiboMiner_config.txt -c longest.transcripts.info.txt
      -o MA_normed -U 密码子 -M RPKM -u 100 -d 400 -l 100 -n 10 -m 1 -e 5 --norm yes \
      -y 100 --type UTR
      其中重要参数为： --类型，分析任一 CDS 或 非翻译区 区域； --norm，是否对读段密度进行了归一化； -y每条转录本所使用的密码子数量； - U，绘制RPF密度于 密码子 水平或 nt 水平； -u 和 - d，定义相对于起始密码子或终止密码子的分析区域范围；-l，CDS的最小长度（即密码子数量）； - M，转录本过滤的模式，二者之一 计数 或 RPKM; - n 分析中CDS区域的最小计数或RPKM值。 - m 标准化区域中CDS的最小计数或RPKM； - e，标准化区域中被排除的密码子数量。
   2. 生成一组PDF文件，用于比较对照细胞和eIF3缺陷细胞中mRNA上的核糖体占据情况。
      PlotMetageneAnalysis -i MA_normed_dataframe.txt -o MA_normed
      -g si-Ctrl,si-eIF3e -r si-Ctrl-1,si-Ctrl-2,si-Ctrl-3__si-eIF3e-1,si-eIF3e-2,si-eIF3e-3 -u 100 -d 400 --mode mean
      注意： 绘制Metagene分析图 生成一组 PDF 文件。有关使用方法的详细信息 元基因分析 和 绘图元基因分析 可在 RiboMiner 网站获取³⁰.

结果

示例核糖体谱分析数据集已存入 GEO 数据库，登录号为 GSE131074。本方案中使用的所有文件和代码均可从以下获取： 补充文件 1-4通过将 RiboCode 应用于一组已发表的核糖体谱型分析数据集²³，我们在经对照处理的MCF-10A细胞中鉴定了新开放阅读框的活跃翻译 EIF3E siRNA。为了筛选出最可能由翻译核糖体结合的RPF测序读段，对测序读段的长度进行了分析，并利用比对到已知翻译基因上的RPF进行全基因平均分析（metagene分析）。读段长度的频率分布显示，大多数RPF的长度为25–35 nt图1A），对应于核糖体所覆盖的核苷酸序列，符合预期。通过分别检查不同长度RPFs的5'端到注释的起始密码子和终止密码子之间的距离，确定其P位点位置（图1B）。28-32nt范围内的RPF读段显示出强烈的3-nt周期性，其P位点位于+12位置

讨论

核糖体谱分析为在基因组尺度上研究细胞中核糖体的活动提供了前所未有的机会。精确解读核糖体谱数据所携带的信息，有助于揭示基因或转录本的哪些区域正处于活跃翻译状态。本逐步操作方案详细指导如何使用 RiboCode 分析核糖体谱数据，包括软件包安装、数据准备、命令执行、结果解释以及数据可视化。RiboCode 的分析结果表明，翻译现象广泛存在，不仅发生在已知编码基因的未注释开放阅读框（ORF）上，也发生在许多此前被认为是非编码的转录本上。下游分析进一步提供了证据，显示核糖体在预测的 ORF 上以每三个核苷酸为一步的方式移动，符合翻译过程的特征；然而，目前尚不清楚这一翻译过程本身或所产生的多肽是否具有生物学功能。尽管如此，对基因组上正在翻译的 ORF 进行准确注释，将为探索此前未被表征的转录本功能带来令人振奋的新机遇³¹。

利用核糖体谱分析数据预测每个开放阅读框（ORF）的编码潜力，高度依赖于从ORF起始密码子到终止密码子范围内每个密码子上P位点密度的3-核苷酸周期性。因此，该方法需要精确确定不同长度测序读段的P位点位置。此类信息并非核糖体谱分...

材料

本文使用的材料清单

姓名	公司	目录编号	评论
运行 Linux 的计算机/服务器	任意	-	-
Anaconda 或 Miniconda	Anaconda	-	Anaconda: https://www.anaconda.com; Miniconda:https://docs.conda.io/en/latest/miniconda.html
R	R 基金会	-	https://www.r-project.org/
Rstudio	Rstudio	-	https://www.rstudio.com/