从头 使用核糖体分析数据识别主动翻译的开放阅读框

摘要

翻译核糖体将每个密码子的三个核苷酸解码为肽。它们沿着mRNA的运动，通过核糖体分析捕获，产生表现出特征性三重态周期性的足迹。该协议描述了如何使用RiboCode从核糖体分析数据中破译这一突出特征，以识别全转录组水平上主动翻译的开放阅读框。

摘要

识别开放阅读框（ORF），特别是那些编码小肽并在特定生理环境下被积极翻译的阅读框，对于上下文依赖性翻译组的全面注释至关重要。核糖体分析是一种检测RNA上翻译核糖体的结合位置和密度的技术，为快速发现全基因组范围内翻译发生的位置提供了一条途径。然而，在生物信息学中，高效、全面地鉴定用于核糖体分析的转化ORFs并非易事。这里描述的是一个易于使用的包，名为RiboCode，旨在从核糖体分析数据中的失真和模糊信号中搜索任何大小的ORF。本文以我们之前发布的数据集为例，提供了整个RiboCode管道的分步说明，从原始数据的预处理到最终输出结果文件的解释。此外，为了评估注释ORF的平移率，还详细描述了每个ORF上核糖体密度的可视化和定量程序。综上所述，本文是对翻译、小ORF和肽相关研究领域的有用和及时的指导。

引言

最近，越来越多的研究表明，从编码基因的ORF和先前注释的基因翻译的肽被广泛生产为非编码基因，例如长非编码RNA（lncRNA）^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（即蛋白质或肽）的产物具有挑战性^，因为这两种方法的检测效率都取决于所产生的蛋白质或肽的长度，丰度和氨基酸组成。核糖体分析是一种在核苷酸分辨率下鉴定mRNA上核糖体占用的技术，它的出现为评估不同转录本的编码潜力提供了一种精确的方法^3，20，21，无论它们的长度和组成如何。使用核糖体分析鉴定主动翻译ORF的一个重要且常用的特征是核糖体从起始密码子到停止密码子在mRNA上的足迹的三核苷酸（3-nt）周期性。然而，核糖体分析数据通常存在几个问题，包括沿ORF的低和稀疏测序读数，高测序噪声和核糖体RNA（rRNA）污染。因此，这些数据产生的扭曲和模糊信号削弱了核糖体在mRNA上足迹的3-nt周期模式，最终使得高置信翻译ORFs的鉴定变得困难。

一个名为"RiboCode"的软件包采用了改进的Wilcoxon签名秩测试和P值积分策略，以检查ORF是否比帧外RPF具有更多的帧内核糖体保护片段（RPM）²²。它被证明对于模拟和真实核糖体分析数据中翻译组的 从头 注释是高效，灵敏和准确的。在这里，我们描述了如何使用该工具从先前研究生成的原始核糖体分析测序数据集中检测潜在的转化^ORF23。这些数据集用于通过比较MCF-10A细胞的核糖体占用谱来探索EIF3亚基"E"（EIF3E）在翻译中的功能，这些细胞转染对照（si-Ctrl）和 EIF3E （si-eIF3e）小干扰RNA（siRNA）。通过将RiboCode应用于这些示例数据集，我们检测到5，633个可能编码小肽或蛋白质的新型ORF。这些ORF根据其相对于编码区域的位置分为各种类型，包括上游ORF（uORFs），下游ORF（dORFs），重叠ORF，来自新型蛋白质编码基因（新型PCG）的ORF以及来自新型非蛋白编码基因（新型NonPCGs）的ORF。与对照细胞相比，EIF3E缺陷细胞中uORFs上的RPF读数密度显着增加，这可能至少部分是由主动翻译核糖体的富集引起的。EIF3E缺陷细胞第25^~ 75^个 密码子区域的局部核糖体积累表明早期翻译伸长受阻。该协议还展示了如何可视化所需区域的RPF密度，以检查已识别ORF上核糖体足迹的3-nt周期模式。这些分析证明了RiboCode在识别翻译ORF和研究翻译监管方面的强大作用。

研究方案

1. 环境设置和 RiboCode 安装

1. 打开一个 Linux 终端窗口并创建一个 conda 环境:
   conda create -n RiboCode python=3.8
2. 切换到创建的环境并安装 RiboCode 和依赖项:
   康达激活核糖代码
   conda install -c bioconda ribocode ribominer sra-tools fastx_toolkit cutadapt bowtie star samtools

2. 数据准备

1. 获取基因组参考文件。
   1. 有关参考序列，请转到位于 https://www.ensembl.org/index.html 的 Ensemble 网站，单击顶部菜单"下载"和左侧菜单"FTP 下载"。在显示的表中，单击"DNA（FASTA）"列中的"FASTA"和"物种为人类"行中的"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法.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> updated.gtf。对于 .gff 格式的注释文件，请使用 AGAT 工具包²⁴ 或任何其他工具转换为 .gtf 格式。
2. 获取 rRNA 序列。
   1. 在 https://genome.ucsc.edu 打开 UCSC 基因组浏览器，然后单击工具|下拉列表中的表浏览器 。
   2. 在打开的页面上，指定 哺乳动物 为分支， 人类 为基因组， 所有表 为组， rmask 为表， 基因组 为区域。对于过滤器，单击" 创建 "转到新页面，并将 repClass 设置为 匹配 rRNA。
   3. 单击" 提交 "，然后将输出格式设置为 序列 ，并将文件名输出为 hg38_rRNA.fa。最后，单击" 获取输出|获取序列 以检索rRNA序列。
3. 从序列读取存档 （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加入ID是通过搜索GSE131074从Gene Expression Omnibus（GEO）网站²⁵ 获得的。

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
   做
   cutadapt -m 15 --match-read-wildcards -a CTGTAGGCACCATCAAT \
   -o ${i}_trimmed.fastq ${i}.fastq
   做
   注意: -a 参数之后的适配器序列将因 cDNA 文库制备情况而异。短于15的读数（由 -m给出）将被丢弃，因为受核糖体保护的片段通常长于此大小。
2. 使用以下步骤去除 rRNA 污染:
   1. rRNA参考序列索引:
      领结-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
      做
      领结 -n 0 -y -a --norc --最佳 --地层 -S -p 4 -l 15 \
      --un=./${i}_noncontam.fastq hg38_rRNA -q ${i}.fastq ${i}.aln
      做
      -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污染）与创建的参考对齐。
   for i in si-Ctrl-1 si-Ctrl-2 si-Ctrl-3 si-eIF3e-1 si-eIF3e-2 si-eIF3e-3
   做
   STAR --runThreadN 8 --outFilterType Normal --outWigType wiggle --outWigStrand Stranded --outWigNorm RPM --outFilterMismatchNmax 1 --outFilterMultimapNmax 1 --genomeDir STAR_hg38_genome --readFilesIn ${i}_noncontam.fastq --outFileNamePrefix ${i}.--outSAMtype BAM SortbyCoordinate --quantMode TranscriptomeSAM GeneCounts --outSAMattributes All
   做
   注意:逆转录酶²⁷经常将未模板化的核苷酸添加到每次读取的5'末端，STAR将在默认情况下执行软剪辑时有效地修剪掉。STAR 的参数在 STAR 手册²⁸ 中进行了描述。
3. 对齐文件排序和索引。
   for i in 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}.Aligned.toTranscriptome.out.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法 -o RiboCode_annot
   注意:此命令从GTF文件收集mRNA转录本的所需信息，并从FASTA文件中提取所有mRNA转录本的序列（每个转录本都是根据GTF文件中定义的结构通过合并外显子来组装的）。
2. 选择特定长度的 RPM 并标识其 P 站点位置。
   for i in si-Ctrl-1 si-Ctrl-2 si-Ctrl-3 si-eIF3e-1 si-eIF3e-2 si-eIF3e-3
   做
   metaplots -a RiboCode_annot -r ${i}.Aligned.toTranscriptome.out.bam \
   -o ${i} -f0_percent 0.35 -pv1 0.001 -pv2 0.001
   做
   注意:此命令绘制每个长度的对齐读取的 5' 端的聚合配置文件，这些读数围绕带注释的平移开始（或停止）密码子。读取长度相关的P位点可以通过检查主读数5'端和起始密码子之间的偏移距离的分布图（例如， 图1B）手动确定。RiboCode还为每个样本生成一个配置文件，其中自动确定显示显着3-nt周期模式的读数的P位点位置。参数 -f0_percent、 -pv1 和 -pv2 定义了比例阈值和 p 值截止值，用于选择读取帧中丰富的 RPF 读数。在此示例中，每个配置文件中手动定义了来自 29、 30 和 31 nt 读取的 5' 端的 + 12、+ 13 和 + 13 核苷酸。
3. 编辑每个示例的配置文件并合并它们
   注意:为了生成一组唯一 ORF 的共识，并确保足够的读取覆盖率以执行后续分析，将合并上一步中所有样本的选定读取。 merged_config.txt文件 （补充文件 1）中定义的特定长度的读取及其 P 站点信息用于在下一步中评估 ORF 的转换潜力。

6 . 从头开始 注释翻译ORF

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，对于在终止密码子上游具有多个起始密码子的转录本，是否使用最长的ORF（从最远端的起始密码子到停止密码子的区域）来评估其翻译潜力。如果设置为 no，将自动确定起始密码子。
   -s，用于ORF识别的规范起始密码子。
   -A，（可选）用于ORF鉴定的非规范起始密码子（例如，CTG，GTG和人类TTG），其线粒体或其他物种的细胞核可能不同²⁹。
   -m，ORF的最小长度（即氨基酸）。
   -o，包含预测 ORF 详细信息的输出文件名的前缀（补充文件 2）。
   -g 和 -b 分别将预测的 ORF 输出为 gtf 或 床 格式。

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
   做
   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
   做
   注意:为了排除ORF开始和结束周围潜在的累积核糖体，在前 15 个（由 -f指定）和最后 5 个密码子（由 -l特异性）中分配的读取次数不计算在内。（可选）计数的 RPM 的长度限制为 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 --起始密码子 ATG -o ENSG00000100902_35292349_35292552_67
   其中 -s 和 -e 指定绘制 ORF 的平移开始和停止位置。 --start-codon 定义了 ORF 的起始密码子，它将出现在图标题中。 -o 定义输出文件名的前缀。

9. （可选）使用核糖矿工进行元基因分析

注意:执行元基因分析，以评估 EIF3E 敲低对已识别的注释ORF的翻译的影响，请按照以下步骤进行:

1. 为RiboMiner生成转录本注释，该注释根据RiboCode生成的注释文件提取每个基因的最长转录本（步骤5.1）。
   输出转录信息 -c RiboCode_annot/transcripts_cds.txt \
   -g Homo_sapiens。GRCh38.104.gtf -f RiboCode_annot/transcripts_sequence.fa \
   -o longest.transcripts.info.txt -O all.transcripts.info.txt
2. 为 RiboMiner 准备配置文件。复制由 RiboCode 的 metaplots 命令生成的配置文件（步骤 5.4），并将其重命名为"RiboMiner_config.txt"。然后，根据 补充文件4中显示的格式对其进行修改。
3. 使用核糖矿工进行元基因分析
   1. 使用 MetageneAnalysis 生成 转录本中 RPF 密度的聚合和平均配置文件。
      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 --范数 是 \
      -y 100 --型 UTR
      其中重要参数为: --type，分析 CDS 或 UTR 区域; --范数，是否归一化读取密度; -y，用于每个转录本的密码子数; -U，在 密码子 水平或 nt 水平上绘制RPF密度; -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
      注意: PlotMetageneAnalysis 生成一组 pdf 文件。有关 MetageneAnalysis 和 PlotMetageneAnalysis 使用的详细信息，请访问RiboMiner网站³⁰。

结果

将示例核糖体分析数据集存入GEO数据库，加入号为GSE131074。此协议中使用的所有文件和代码均可从补充文件 1-4 获得。通过将RiboCode应用于一组已发表的核糖体分析数据集²³，我们确定了在用对照和EIF3E siRNA处理的MCF-10A细胞中主动翻译的新型ORF。为了选择最有可能被翻译核糖体结合的RPF读数，检查了测序读数的长度，并使用映射在已知?...

讨论

核糖体分析为在基因组尺度上研究核糖体在细胞中的作用提供了前所未有的机会。精确破译核糖体分析数据携带的信息可以深入了解基因或转录本的哪些区域正在积极翻译。此分步协议提供了有关如何使用 RiboCode 详细分析核糖体分析数据的指导，包括软件包安装、数据准备、命令执行、结果说明和数据可视化。RiboCode的分析结果表明，翻译是普遍存在的，并且发生在编码基因的未注释ORF和许多先?...

材料

Name	Company	Catalog Number	Comments
A computer/server running Linux	Any	-	-
Anaconda or Miniconda	Anaconda	-	Anaconda: https://www.anaconda.com; Miniconda:https://docs.conda.io/en/latest/miniconda.html
R	R Foundation	-	https://www.r-project.org/
Rstudio	Rstudio	-	https://www.rstudio.com/