1 背景及分析流程

1.1 背景介绍

信使RNA（mRNA）转化为蛋白质以及将所得蛋白质折叠成活性形式是几乎每个细胞过程的先决条件，并不是所有mRNA都能与核糖体结合翻译出蛋白质，因此需要翻译组学方法鉴定正在翻译中的mRNA。mRNA的翻译通常从募集由核糖体小亚基、三元复合物和起始因子形成的起始前复合物（preinitiation complex）开始。在真核生物中，通过真核起始因子eIF4E识别5' m7G帽，将前起始复合物募集到mRNA的5'末端。然后，起始前复合体扫描5'非翻译区域 （5' UTR），直到识别起始位点。然而，起始位点可能绕过5' UTR 中的 AUG或者选择5' UTR 中的其他 AUG 样三联体作为起始密码子。一旦起始位点被识别，大多数起始因子被释放，核糖体大亚基加入起始复合物。在释放起始因子eIF5B后，起始核糖体进入延伸循环，真核延伸因子 eEF1 将带氨基酸的转移 RNA （tRNA） 递送至核糖体 A 位点，eEF2 催化核糖体易位。当核糖体到达终止密码子时翻译终止，与起始位点选择类似，终止密码子识别的保真度受顺式元件和反式调节因子的控制[1]。核糖体印记测序(Ribosome profile sequencing)，使用NGS测序鉴定核糖体保护的mRNA片段，可用于对体内翻译和偶联、共翻译事件进行全局分析[2]。

1.2 实验原理

核糖体分析依赖于核糖体足迹（ribosome footprints, RFPs；又称ribosome protected fragments, RPFs）的测序，核糖体足迹是 mRNA 的短片段（通常为约30个核苷酸 [nt]）。首先，用翻译抑制剂处理细胞或组织，使得细胞翻译活动停止。再用核酸酶处理核糖体新生肽链复合物，降解没有核糖体覆盖的mRNA片段，去除核糖体。RFPs片段被核糖体物理包裹并屏蔽，不被核酸酶消化，它们被转化为DNA片段库，检测合格后进行文库制备[3]。构建好的文库经过文库质控合格后，使用Illumina Platform进行测序，Ribo-Seq测序策略一般使用SE 150(单端150bp测序)。

图1.1 Ribo-Seq实验原理[4]

1.3 生信分析流程

信息分析主要分为7部分：1、数据预处理。去接头序列、污染序列、低质量碱基，获得 clean data序列，并进行相关数据统计；2、根据提供的rRNA参考序列，去除clean data的rRNA序列，之后定位到参考基因组上，得到 bam文件，并去除重复序列，保留唯一比对的序列；3、分析RPFs（核糖体保护片段）；4、鉴定样本中ORF特别是上游的uORF；5、RPFs表达量统计；6、差异表达基因分析及功能富集分析；7、与mRNA联合分析，计算各组翻译效率与差异以及各基因转录和翻译变化。

图1.2 Ribo-Seq分析流程

1.4 样本信息

进行分析的样本的信息，由客户或实验人员提供样本详细的处理操作信息。

表1.1 样本信息

2 数据处理

2.1 fastq文件说明

Illunima等高通量测序平台得到的原始图像数据文件，经碱基识别（Base Calling）转化为原始测序序列（Sequenced Reads），我们称之为 Raw Reads。Raw Reads 以 FASTQ 格式存储，包含序列以及对应的测序质量信息。Fastq 数据示例如下：

图2.1 FASTQ格式文件说明

每条序列由 4 行组成，第 1 行是序列名称；第 3 行也是序列名称，一般省略，只保留“+”；第 2 行是序列；第 4 行是序列的测序质量。第四行中每个字符对应的 ASCII 值减去 33，即为第二行对应碱基的测序质量值: Q = −10log10(E)。

测序错误率与对应字符换算示例如下所示：

表2.1 测序错误率与测序质量值对应关系

2.2 Raw Data 数据评估

对原始数据，使用软件 FastQC（version: 0.11.5） 质控处理，详细结果见 ./sup/01.qc/raw，结果如下所示。

2.3 clean Data 数据评估

使用fastp软件[5]进行数据过滤，截掉末尾的adapter序列，保留最短长度15bp。对过滤数据，使用软件 FastQC（version: 0.11.5）质控处理，详细结果见 ./sup/01.qc/clean，结果如下所示。

2.4 Clean 结果统计

对原始数据和过滤处理得到高质量的数据（clean reads）进行数据统计，得到数据基本信息。统计结果如下表：

表2.2 clean data数据统计表

2.5 质量评估Q&A

3 比对分析

3.1 去除rRNA

测序得到的reads中有很多rRNA序列，为了去除它们，我们使用软件bowtie（version :1.3.1）[7]将测序数据与参考rRNA序列进行比对，只保留没有比对成功的序列（即非rRNA序列）进行下一步分析。

3.2 比对参考基因组

过滤rRNA之后，我们将测序数据与参考基因组比对，使用软件为STAR（version:2.7.11a）[8]。我们将对比对得到的唯一比对reads用于后续分析。

表3.1 STAR比对情况

1. sample：样本名
2. after_remove_rRNA_reads：去除rRNA后reads总数
3. uniqMapReads：唯一比对无重复reads数量
4. mapReads：比对上参考基因组reads总数

3.3 reads在染色体上分布

为了直观的展示染色体长度和Reads总数以及Reads覆盖情况，将比对到基因组上各条染色体的Reads进行唯一比对和去重处理后进行统计，以200k长度为1个bin计算各个bin内部比对到参考基因组上的Reads数目的均值（作图时取对数）。通常情况下染色体越长定位到该染色体的Reads数目越多，DNA项目的分布较为均匀，RNA项目由于获得的多是exon区域，分布有一定间隔。如果染色体数目过多，则只会展示部分染色体。reads 在整个基因组上的分布详细结果见./sup/02.map/readsOnGenome，结果如图所示。

图3.1 reads on genome

3.4 reads 区域分布

统计reads在基因组上各类元件的分布，RNA项目reads主要分布在CDS和intron区，详细结果见./sup/02.map/mapRegion，结果如图所示。

表3.2 基因组各区域分布统计

图3.2 reads region

4 RPF分析

4.1 RPF长度统计

RPF长度普遍在28-30 nt左右，我们统计每个样本RPF长度并绘制直方图，如果RPF长度集中在过长或过短的区间，可能代表实验结果不理想。Frame0：从mRNA的第一个核苷酸开始翻译，每三个核苷酸为一组，形成一个阅读框；Frame1：从mRNA的第二个核苷酸开始翻译，同样每三个核苷酸为一组，形成另一个阅读框；Frame2：从mRNA的第三个核苷酸开始翻译，同样每三个核苷酸为一组，形成第三个阅读框。

图4.1 RPF长度直方图

表4.1 RPF长度统计

1. sample:样本名称
2. Length:RPF reads长度
3. proportion:该长度reads占比对上reads的百分比
4. f0_sum:frame0的数量
5. f1_sum:frame1的数量
6. f2_sum:frame2的数量
7. f0_percent：frame0所占百分比

4.2 P-site分析

使用软件RiboCode（version:1.2.13）[9]中的 metaplots 命令，选择最有可能来自翻译中核糖体的 RPF reads的长度范围，并通过已有注释中基因上比对到的RPF reads进行meta-gene analysis识别不同长度的 RPF 的 P-site。具体来说，对于具有特定长度的每组 RPF 读段，计算从它们的 5' 端到注释的起始和终止密码子的距离并总结为直方图。由于Ribo-seq亚密码子分辨率，核糖体分析揭示了 RPF reads中 80S 核糖体的肽基位点（P-site）的精确位置。因此，由翻译核糖体产生的 RPFs应沿ORF表现出 3-nt 周期性，且大多数reads应该在frame0处富集。详细结果见./sup/03.rpf/psite。

下图横坐标为以预测的P-site为中心左右50 bp RPF reads的分布，纵坐标为该位置比对到的reads数量。蓝色代表属于frame0的reads，红色代表frame1，绿色代表frame2。

图4.2 P-site区域RPFs数量

5 翻译停滞分析

5.1 元基因分析

我们使用软件RiboMiner（version:0.2）[10]对各样本进行元基因分析（Metagene analysis），可以检查核糖体占有率是否发生变化。如果编码区出现RPFs异常富集，说明此处翻译复合体没有顺利延伸，出现了翻译停滞。

各样本转录本元基因分析，横坐标为标准化后转录本位置，纵坐标为该位置reads密度。

图5.1 转录本元基因分析

各样本极性分数分布图，横坐标为核糖体极性分数，纵坐标为相关基因的数量。如果峰值不在0处，出现在负值说明核糖体向5`端移动，出现在正值说明核糖体向3`端移动。

图5.2 极性分数分布图

为了准确体现各样本在起始/终止密码子处核糖体密度分布，我们绘制了这两个区域核糖体足迹分布图，0坐标代表起始/终止密码子，纵坐标为该位置reads密度。

图5.3 起始密码子核糖体足迹分布

图5.4 终止密码子核糖体足迹分布

5.2 翻译停滞基因

上一步的分析显示了出现翻译停滞的区域，我们使用软件RiboMiner比较各个样本每个转录本翻译停滞区域密码子上的reads密度，选择reads比率大于1.5作为翻译停滞基因，详细结果见./sup/04.meta/stallgene。分析结果如下：

1. gene:基因ID
2. trans_id:转录本ID
3. codon：该转录本第几个密码子出现了差异富集
4. log2FoldChange：实验组与对照组之间log2化的差异倍数
5. pvalue：差异的显著性
6. status：实验组相较于对照组该位置上下调情况
7. 其他：DESeq2结果

6 ORF分析

我们使用软件RiboCode根据PRF鉴定ORF，ORF的结果可以分成7类：1.annotated，与已经注释到的ORF重合，且起始密码子相同；2.uORF，在已注释CDS上游且不与CDS重合；3.dORF，在已注释CDS下游且不与CDS重合；4.Overlap_uORF，在已注释CDS上游且与CDS重合；5.Overlap_dORF，在已注释CDS下游且与CDS重合；6.Internal，在已注释CDS的内部ORF，但在不同的ORF中；7.novel，来自非编码基因或编码基因的非编码转录物。详细结果请见./sup/05.orf/。sup/05.orf/gtf/内gtf文件为鉴定的ORF结果。

图6.1 ORF鉴定结果统计

6.1 uORF分析

在真核生物体内，uORF可能起着相似的翻译调节机制，它与下游的mORF竞争翻译起始复合物或影响mORF的复苏。而在uORF中可能存在并未以传统ATG为起始密码子而以其他碱基序列为密码子的情况，导致了翻译提前发生。

图6.2 uORF起始密码子统计

每个uORF起始密码子情况请见下表：

1. ORF_ID:本次分析时ORF编号
2. StartCodon：起始密码子
3. ORF_type：ORF类型
4. transcript_id：转录本ID
5. gene_id：基因ID
6. chrom：染色体
7. ORF_gstart：ORF在基因组起始坐标
8. ORF_gstop：ORF在基因组终止坐标

7 表达量分析

7.1 样本表达量分析

根据比对的reads，保留uniq比对（在基因组上仅比对到一个位置），去掉重复的reads(reads 片段的序列相同)，然后使用featureCounts(v2.0.1)[11]对样本进行表达量分析。详细结果见 ./sup/06.exp。

图7.1 样本表达量分布

7.2 PCA分析

根据基因在各个样本中的RPKM（Reads Per Kilobase per Million mapped reads,每千个碱基的转录每百万映射的reads）的概念对每一个基因进行定量，通过RPKM计算样本的PCA分析。详细结果见 ./sup/06.exp。

图7.2 样本PCA分析

7.3 样本相关性

根据基因在各个样本中的RPKM的概念对每一个基因进行定量，通过RPKM计算样本的相关性分析。详细结果见 ./sup/06.exp。

图7.3 样本相关性

7.4 表达量Q&A

8 差异表达分析

统计各个基因的reads数，根据reads，使用DESeq2[12]方法，分析各组的差异表达基因。参数一般选择FC(fold change) = 2，pValue = 0.05。差异基因列表见sup/07.deg/*_vs_*文件夹下各组的 *_ann_deg.txt 文件。

8.1 火山图

各组样本差异表达分析的火山图。

图8.1 火山图

8.2 差异基因统计

对各组差异基因进行汇总统计，详细结果见 ./sup/07.deg。

表8.1 差异基因统计表

图8.2 差异基因统计

8.3 差异基因样本相关性分析

根据差异表达的基因在各个样本中的RPKM（Reads Per Kilobase per Million mapped reads）的概念对每一个差异表达基因进行定量，通过RPKM计算样本的相关性分析。详细结果见 ./sup/07.deg。

图8.3 差异基因样本相关性

8.4 差异基因样本PCA分析

根据差异表达的基因在各个样本中的RPKM（Reads Per Kilobase per Million mapped reads）的概念对每一个差异表达基因进行定量，通过RPKM计算样本的PCA分析。详细结果见 ./sup/07.deg。

图8.4 差异基因样本PCA分析

8.5 差异基因样本heatMap分析

选取差异基因，对这些基因在各样本中的表达情况进行heatMap分析，详细结果见 ./sup/07.deg。

图8.5 差异基因样本heatMap分析

8.6 差异表达分析Q&A

9 与转录组联合分析

9.1 翻译效率分析

我们比较蛋白组数据和转录组数据时，往往会发现转录本丰度不等于蛋白质丰度，翻译效率（Translation Efficiency，TE）的改变是基因翻译水平调控的主要形式，直接影响蛋白质的产量，其计算公式为：TE＝(RPKM in RiboProfiling) / (RPKM in RNAseq)。完整结果请见 ./sup/08.te，各比较组翻译效率分析结果为*_TE.ann.txt，各比较组翻译效率差异基因如下表所示：

1. ID：基因ID
2. TEstatus：比较组中翻译效率上下调情况
3. log2FC_TE：翻译效率log2化倍数
4. pvalue：差异显著性p值
5. p.adjust：BH校正后差异显著性p值
6. Name：基因名称
7. GO：基因对应GO条目
8. KEGG：基因对应KEGG通路
9. Nr：NR数据库注释

各组样本差异翻译效率的火山图。

图9.1 TE差异基因火山图

9.2 转录与翻译差异基因交集

我们将mRNA表达差异基因和RPF表达差异基因取交集，绘制venn图。

为了进一步表现翻译组和转录组数据整合结果，我们根据mRNA和RPF差异倍数绘制转录翻译九象限图。下图横坐标为转录差异基因log2FC值，纵坐标为翻译差异基因log2FC值。从左往右第一条竖线左边为转录下调基因，第二条竖线右边为转录上调基因，两条竖线中间为转录没有差异的基因；从上往下第一条横线上面为翻译上调基因，第二条横线下面为翻译下调基因，两条横线中间为翻译没有差异基因。基因可以分为5类：

1.仅在mRNA水平上表达改变的基因（transcription only），基因仅在转录水平上存在差异，但它们的RPF水平没有显著差异；

2.仅在TE水平上表达变化的基因（translation only），基因的两个等位基因仅在RPF水平上存在差异，但在转录水平上没有显著差异；

3.反向变化基因（opposite change），转录和翻译水平上出现反向变化；

4.同向变化基因（homodirectional），转录和翻译水平上出现同向变化；

5.剩下的无差异基因。

图9.2 转录翻译九象限图

1. ID：基因ID
2. status：该基因转录翻译变化情况
3. mRNA_log2FC：mRNA中比较组log2差异倍数
4. RPF_log2FC:RPF中比较组log2差异倍数
5. Name：基因名称
6. GO：基因对应GO条目
7. KEGG：基因对应KEGG通路
8. Nr：NR数据库注释

9.3 与mRNA联合分析Q&A

10 差异基因富集分析

功能聚类分析介绍：

GO 功能分析

基因本体论（Gene Ontology）简称GO，是一个国际标准化的基因功能分类体系，包含生物学领域知识体系本质的表示形式，本体通常由一组类（或术语或概念）组成，分为3大类：分子功能（Molecular Function）、细胞组分（Cellular Component）、生物过程（Biological Process）。对研究的物种使用相关GO注释数据库或利用blast，得到每个基因对应的GO term，根据基因注释信息，统计每个基因所在的GO term ，根据每个GO term的基因数目，以及背景中此GO term的基因数目，用Fisher Exact Test分析每个GO term的显著性。

KEGG 功能分析

京都基因与基因组百科全书（Kyoto encyclopedia of genes and genomes）简称KEGG，是系统分析基因功能、基因组信息的数据库，它整合了基因组学、生物化学以及系统功能组学的信息，有助于研究者把基因及表达信息作为一个整体网络进行研究。对差异基因进行KEGG功能富集分析，研究的物种使用相关KEGG注释数据库或利用blast，得到每个基因对应的KEGG pathway。根据基因注释信息，统计每个基因所在的KEGG pathway，根据每个KEGG pathway的基因数目，以及背景中此KEGG pathway的基因数目，用Fisher Exact Test分析每个KEGG pathway的显著性。

10.1 RPF表达差异基因

我们对RPF表达差异基因进行GO与KEGG富集分析。基因注释结果请见./sup/10.gokegg/rpf/。

10.1.1 GO 功能分析

GO 功能聚类气泡图展示。纵坐标是GO Term 名称，横坐标是对应GO Term 中检出的基因占背景基因的个数，颜色代表显著性pvalue。

图10.1 GO功能聚类气泡图

GO 功能聚类bar图展示。纵坐标是GO Term 名称，横坐标是对应GO Term 中检出的基因个数，颜色代表显著性p.adjust。

图10.2 GO功能聚类bar图

GO 功能聚类wego图展示。横坐标坐标是GO Term 名称，纵坐标是对应GO Term 中检出的基因个数，取对数。根据GO分成了3大类，分别是分子功能（Molecular Function，MF）、细胞组分（Cellular Component，CC）、生物过程（Biological Process，BP）。

图10.3 GO功能聚类wego图

10.1.2 KEGG 功能分析

对于RNA项目，根据基因上调（up），下调(down)和上下调并集(all)分别进行分析。详细结果见./sup/10.gokegg/*_kegg。

KEGG功能聚类气泡图展示。纵坐标是KEGG pathway名称，横坐标是对应KEGG pathway中检出的基因占背景基因的个数，颜色代表显著性pvalue。

图10.4 功能聚类气泡图

KEGG功能聚类bar图展示。纵坐标是KEGG pathway名称，横坐标是对应KEGG pathway中检出的基因个数，颜色代表修正后的p.adjust。

图10.5 功能聚类气泡图

KEGG Pathway主要划分为7类：分别为Metabolism， Genetic information Processing等。其中每类又分为二、三、四级子条目。功能聚类分类图课展示pathway所属一级条目和二级条目的情况。

图10.6 功能聚类分类图

KEGG通路图具有强大的图形功能，将代谢、调控、通路、生化、疾病、药物等相关的分子相互作用整合到一起并形成关系网络，通过功能富集分析得到的KEGG通路示例图如下，部分富集结果中可能没有KEGG通路图。具体文件见./sup/10.gokegg/rpf/*_kegg 文件夹各组下面的map文件夹。

图10.7 KEGG通路示意图

10.2 翻译效率差异差异基因

我们对翻译效率差异基因进行GO与KEGG富集分析。基因注释结果请见./sup/10.gokegg/dteg/。

10.2.1 GO 功能分析

GO 功能聚类气泡图展示。纵坐标是GO Term 名称，横坐标是对应GO Term 中检出的基因占背景基因的个数，颜色代表显著性pvalue。

图10.8 GO功能聚类气泡图

GO 功能聚类bar图展示。纵坐标是GO Term 名称，横坐标是对应GO Term 中检出的基因个数，颜色代表显著性p.adjust。

图10.9 GO功能聚类bar图

GO 功能聚类wego图展示。横坐标坐标是GO Term 名称，纵坐标是对应GO Term 中检出的基因个数，取对数。根据GO分成了3大类，分别是分子功能（Molecular Function，MF）、细胞组分（Cellular Component，CC）、生物过程（Biological Process，BP）。

图10.10 GO功能聚类wego图

10.2.2 KEGG 功能分析

对于RNA项目，根据基因上调（up），下调(down)和上下调并集(all)分别进行分析。详细结果见./sup/10.gokegg/dteg/*_kegg。

KEGG功能聚类气泡图展示。纵坐标是KEGG pathway名称，横坐标是对应KEGG pathway中检出的基因占背景基因的个数，颜色代表显著性pvalue。

图10.11 功能聚类气泡图

KEGG功能聚类bar图展示。纵坐标是KEGG pathway名称，横坐标是对应KEGG pathway中检出的基因个数，颜色代表修正后的p.adjust。

图10.12 功能聚类气泡图

KEGG Pathway主要划分为7类：分别为Metabolism， Genetic information Processing等。其中每类又分为二、三、四级子条目。功能聚类分类图课展示pathway所属一级条目和二级条目的情况。

图10.13 功能聚类分类图

KEGG通路图具有强大的图形功能，将代谢、调控、通路、生化、疾病、药物等相关的分子相互作用整合到一起并形成关系网络，通过功能富集分析得到的KEGG通路示例图如下，部分富集结果中可能没有KEGG通路图。具体文件见./sup/10.gokegg/dteg/*_kegg 文件夹各组下面的map文件夹。

图10.14 KEGG通路示意图

10.3 转录翻译变化基因

我们对转录翻译变化基因进行GO与KEGG富集分析。基因注释结果请见./sup/10.gokegg/with_mrna/。其中“homo”代表转录和翻译同向变化，“opposite”表达转录和翻译反向变化。

10.3.1 GO 功能分析

GO 功能聚类气泡图展示。纵坐标是GO Term 名称，横坐标是对应GO Term 中检出的基因占背景基因的个数，颜色代表显著性pvalue。

图10.15 GO功能聚类气泡图

GO 功能聚类bar图展示。纵坐标是GO Term 名称，横坐标是对应GO Term 中检出的基因个数，颜色代表显著性p.adjust。

图10.16 GO功能聚类bar图

GO 功能聚类wego图展示。横坐标坐标是GO Term 名称，纵坐标是对应GO Term 中检出的基因个数，取对数。根据GO分成了3大类，分别是分子功能（Molecular Function，MF）、细胞组分（Cellular Component，CC）、生物过程（Biological Process，BP）。

图10.17 GO功能聚类wego图

10.3.2 KEGG 功能分析

KEGG详细结果见./sup/10.gokegg/with_mrna/*_kegg。

KEGG功能聚类气泡图展示。纵坐标是KEGG pathway名称，横坐标是对应KEGG pathway中检出的基因占背景基因的个数，颜色代表显著性pvalue。

图10.18 功能聚类气泡图

KEGG功能聚类bar图展示。纵坐标是KEGG pathway名称，横坐标是对应KEGG pathway中检出的基因个数，颜色代表修正后的p.adjust。

图10.19 功能聚类气泡图

KEGG Pathway主要划分为7类：分别为Metabolism， Genetic information Processing等。其中每类又分为二、三、四级子条目。功能聚类分类图课展示pathway所属一级条目和二级条目的情况。

图10.20 功能聚类分类图

KEGG通路图具有强大的图形功能，将代谢、调控、通路、生化、疾病、药物等相关的分子相互作用整合到一起并形成关系网络，通过功能富集分析得到的KEGG通路示例图如下，部分富集结果中可能没有KEGG通路图。具体文件见./sup/10.gokegg/with_mrna/*_kegg 文件夹各组下面的map文件夹。

图10.21 KEGG通路示意图

10.4 富集分析Q&A

11 备注

11.1 软件列表

此报告中涉及到的分析所用的软件及相关信息如下。

表11.1 软件说明

11.2 文件目录

报告的整体文件结构如下表，部分文件可能由于分析策略的不同和分析结果的差异，可能会没有。

12 参考文献