技术服务描述

miRNA结题报告

生信部

2021年01月28日

1. 建库测序工作流程

1.1. 建库测序流程

  MicroRNA (miRNA) 是一类长度约为20-24个核苷酸长度的具有调控功能的非编码RNA，miRNA 是生物体内一类重要的特殊分子，诱导基因沉默，参与细胞生长、发育、基因转录和翻译等诸多生命活动的调控过程。

  从RNA样品到最终分析结果数据，需要经过样品检测、建库、测序和信息分析等过程，其中测序过程我们采用高通量测序illumina HiSeq/MiSeq测序平台；illumina测序得到的原始图像数据文件经碱基识别（Base Calling）分析转化为原始测序序列（Sequenced Reads），我们称之为Raw Data或Raw Reads，结果以FASTQ（简称为fq）文件格式存储，其中包含测序序列（reads）的序列信息以及其对应的测序质量信息。

  样品质量检测：利用NanoDrop 2000分光光度计对RNA样品进行初步定量，Aglient 2100对样品浓度精确定量。RNA样品总量、浓度、完整性（RIN值，RNA Integrity Number）、纯度（OD260/OD280比值）符合收样立项标准后，进入下一步的文库构建。

  文库构建：Total RNA PEG (polyethylene glycol)沉淀，连接3’接头，PAGE (polyacrylamide gel electrophoresis) 胶分离回收，连接5’接头，反转录，PCR扩增，文库切胶回收。

  文库质检：Q-PCR方法对文库进行精确定量，文库有效浓度>2nM即可上机测序。

  上机测序：将各文库按照有效浓度及目标下机数据量的需求pooling后，进行SE75测序。

图1: Small RNA 文库构建 Workflow

  Small RNA测序的接头信息：

  RNA 5‘Adapter (RA5), part:

  5’-GTTCAGAGTTCTACAGTCCGACGATC-3’

  RNA 3’Adapter (RA3), part:

  5’-AGATCGGAAGAGCACACGTCT-3’ (1ug流程)

  5’-CTGTAGGCACCATCAATCAGATCGGAAGAGCACACGTCTG-3’ (10ug流程)

2. 信息分析流程

  Raw Data通过去接头、去低质量、去污染、统计序列长度分布等过程完成初级分析；将初级分析得到的clean reads与ncRNA库比对，分类注释ncRNA，并去除reads中的rNRA、tRNA、snRNA、snoRNA等；然后再将reads比对到参考基因组、miRBase上，计算miRNA表达量，注释已知miRNA、预测新的miRNA；然后进行差异miRNA、靶基因预测、功能注释等后续分析。

miRNA生物信息分析流程图

2.1. miRNA-seq 测序数据质量评估

2.1.1. 数据质量评估 及 数据过滤

  对于测序得到的rawdata，我们首先使用fastqc进行数据质量评估，通过结果，我们可以了解到测序数据reads质量、长度、碱基分布，序列重复频次等基本信息。

  测序得到的raw data，里面含有带接头的、低质量的reads，为了保证后续信息分析的质量，必须对raw reads进行处理，得到clean reads。

  raw data的过滤步骤如下： 去除低质量 reads; 去除有5’接头污染的reads; 去除没有3’接头序列; 去除没有插入片段的reads; 去除长度小于16nt的reads;

  数据经过以上过滤后，我们对数据做一些基本统计。

结果：

2.2. 数据长度分布

  一般来说，小RNA的长度区间为18～30nt，数据过滤完成后，我们统计了clean reads数目及各自占总reads的百分比，并绘制reads长度分布图。reads长度分布图能帮助我们判断小RNA的种类，如miRNA集中在21或22nt，siRNA集中在24nt，piRNA集中在30nt。

Clean Data Reads 长度分布统计直方图：

Figure 2.2: Small RNA Clean Reads长度分布图

Clean Reads长度分布直方图，横坐标为Reads的长度，纵坐标为样品中对应该长度的total Reads数量。

2.3. miRNA基因组比对分析

  我们将Clean Reads 与 mature miRNA, miRNA hairpin, mRNA, mature & primary tRNA, snoRNA, rRNA, other non-coding RNA, and (optional) known RNA 使用bowtie进行比对搜索，然后统计reads的比对情况，分析样本中各类ncRNA的数目及比例情况，统计结果如下。过滤出待分析的miRNA以便后续分析，将clean reads中ncRNA尽可能地去除，便于后续 miRNA 检测及预测。

Reads类型分布图

图显示的是各类型的RNA reads和miRNA reads占总distinct reads的比例。

与参考基因组比对质控结果：

2.4. 已知 miRNA 的检测以及 novel miRNA 预测分析

  首先用去除了ncRNA后的Reads比对miRbase数据库，比对上miRbase中 reads 用于检测已知成熟的miRNA，并统计其表达量和RPM值。我们采用miRge软件来进行 miRNA reads 检测注释、novel miRNA的预测、miRNA表达量的统计。

miRNA分析结果汇总：

2.4.1. miRNA异构体分析

  随着高通量测序技术的发展，早期认为的miRNA loci仅产生1条miRNA成熟体序列这一观点被颠覆。对于某一miRNA来说，它并不是单一的序列，而是由一系列长度/序列及表达不同的异构体 (isomiRs) 组成。这些isomiRs表达多样且序列多样，甚至引入多样的5'端及种子区域。特定miRNA位点在疾病组织中可具有异常的表达模式，现已证实，部分isomiRs具有重要的生物学功能。

  IsomiRs目前主要分为三类：5′isomiRs, 3′isomiRs, polymorphic isomiRs。IsomiRs的产生机制主要有：miRNA加工和成熟过程中Darsha和Dicer酶的不精确或选择性剪切；3'端核苷酸添加；RNA编辑和单核苷酸多态性( single nucleotide polymorphism，SNP)。主要表现为：5'端修剪；3' 端修剪；3'端核苷酸附加和碱基替换。其中，5'端修剪和碱基替换可发生在种子区域内部，产生“种子转移”（seed shifting）。

  我们通过序列与miRbase比对的情况，将各个miRNA Reads进行isomiR识别，并同时汇总至miRNA reads 注释结果中，分类统计，并选取表达量较高的前20个mirna进行各类IsomiRs表达丰度统计。

图2.4.1 所有样本的变体类型分布（isomiRs）

图2.4.2 各个样本的变体分布TOP20的丰度统计

2.4.2. miRNA 碱基编辑

  成熟miRNA序列的第2-8个碱基被称作“种子”序列，保守性很高。若在这一区域发生碱基突变，则可能改变miRNA的靶基因作用位点。我们通过将miRNA序列与miRBase中已知miRNA前体以及成熟miRNA序列进行比对(只允许一个位点的错配)找出发生了碱基突变的miRNA，统计汇总在该miRNA的碱基突变位点，以及发生碱基编辑的Reads数量及百分比，并对统计概况结果进行可视化。

图2.4.3 鉴定为已存在mirna碱基编辑的在各个样本中的相应readCount占比统计图

miRNA碱基编辑分析结果详细说明：

1. A2I在各个样本中统计汇总表

• 结果文件：02.miRNA/2.a2IEditing/a2IEditing.report.csv

• 表头说明：

  结果说明	结果文件
  基本注释结果：
  miRNA reads 变体类型注释结果	02.miRNA/1.isomiR/sample_miRge3.gff
  miRNA reads 变体类型可视化统计	02.miRNA/1.isomiR/*svg
  miRNA碱基编辑分析结果：
  A2I在各个样本中统计汇总表	02.miRNA/2.a2IEditing/a2IEditing.report.csv
  A2I在各个样本中统计汇总表的筛选及简化信息 ( 即汇总表的*.RPM.mismatch、*.AtoI.percentage列 )	02.miRNA/2.a2IEditing/a2IEditing.report.newform.csv
  A2I分析的序列及统计详情信息	02.miRNA/2.a2IEditing/a2IEditing.detail.txt
  二级结构预测结果：
  novel miRNA 二级结构预测结果	02.miRNA/3.novel_predict/*.R1_novel_miRNAs/
  novel miRNA 二级结构预测详细信息表格汇总	02.miRNA/3.novel_predict/*.R1_novel_miRNAs/*_novel_miRNAs_report.csv

2. A2I在各个样本中统计汇总表的筛选及简化信息
   ( 即汇总表的*.RPM.mismatch、*.AtoI.percentage列，该表对应上面的可视化结果 )

• 结果文件：02.miRNA/2.a2IEditing/a2IEditing.report.newform.csv

• 表头说明：




结果说明	结果文件
基本注释结果：
miRNA reads 变体类型注释结果	02.miRNA/1.isomiR/sample_miRge3.gff
miRNA reads 变体类型可视化统计	02.miRNA/1.isomiR/*svg
miRNA碱基编辑分析结果：
A2I在各个样本中统计汇总表	02.miRNA/2.a2IEditing/a2IEditing.report.csv
A2I在各个样本中统计汇总表的筛选及简化信息 ( 即汇总表的*.RPM.mismatch、*.AtoI.percentage列 )	02.miRNA/2.a2IEditing/a2IEditing.report.newform.csv
A2I分析的序列及统计详情信息	02.miRNA/2.a2IEditing/a2IEditing.detail.txt
二级结构预测结果：
novel miRNA 二级结构预测结果	02.miRNA/3.novel_predict/*.R1_novel_miRNAs/
novel miRNA 二级结构预测详细信息表格汇总	02.miRNA/3.novel_predict/*.R1_novel_miRNAs/*_novel_miRNAs_report.csv




3. A2I分析的序列及统计详情信息

• 结果文件：02.miRNA/2.a2IEditing/a2IEditing.detail.txt

• 结果说明：

  ①每一个miRNA块，包含3小块内容：原始序列、筛选后的序列统计、筛选后的序列。重复的miRNA块代表不同样本，顺序为 A2I在各个样本中统计汇总表 的样本信息。

  ②每一行序列包含3个信息， miRNA sequence（reads序列），*.readCount（Count计数），Filter（是否筛选到后续进行碱基编辑分析）

2.4.3. novel miRNA 预测及二级结构分析

  miRge根据序列特征，碱基配对原则，最小自由能等特征模型，进行新颖的miRNA及二级结构预测。预测示意图如下。结果见：02.miRNA/3.novel_predict/

图2.4.4 新miRNA的二级结构预测示意图

2.5. miRNA 定量分析

2.5.1. 样本间相关性分析

  生物学重复通常是任何生物学实验所必须的，目前主流期刊也基本要求生物学重复。生物学重复主要有两个用途：一个是证明所涉及的生物学实验操作不是偶然，而是可重复的。另一个是为了确保后续的差异分析得到更可靠的结果。样品间基因表达水平相关性是检验实验可靠性和样本选择是否合理的重要指标。相关系数越接近1，表明样品之间表达模式的相似度越高。Encode计划建议皮尔逊相关系数的平方(R2)大于0.92(理想的取样和实验条件下)。具体的项目操作中，我们要求生物学重复样品间R2至少要大于0.8，否则需要对样品做出合适的解释，或者重新进行实验。根据各样本所有基因的表达值计算组内及组间样本的相关性系数，绘制成热图，可直观显示组间样本差异及组内样本重复情况。样本间相关性系数越高，其表达模式越为接近，样本相关性热图如下图所示。

图 2.5.1. 样本间相关性热图

图中横纵坐标为各样本相关系数的平方

结果文件：

样本间相关性热图结果：04.DE/Quant/cor_pheatmap*

2.5.2. 主成分分析

  主成分分析（PCA）也常用来评估组间差异及组内样本重复情况，PCA采用线性代数的计算方法，对数以万计的基因变量进行降维及主成分提取。我们对所有样本的基因表达值进行PCA分析，如下图所示。理想条件下，PCA图中，组间样本应该分散，组内样本应该聚在一起。

图 2.5.2. 主成分分析结果图

图中横坐标为第一主成分，纵坐标为第二主成分

结果文件：

主成分分析结果：04.DE/Quant/pca*

2.6. miRNA 差异分析

  基因表达定量完成后，需要对其表达数据进行统计学分析，筛选样本在不同状态下表达水平显著差异的基因。差异分析主要分为三个步骤。

• 首先对原始的readcount进行标准化（normalization），主要是对测序深度的校正。

• 然后统计学模型进行假设检验概率（pvalue）的计算

• 最后进行多重假设检验校正，得到FDR值（错误发现率，padj是其常见形式)[1-2]。

  针对不同的实验情况，我们选用合适的软件进行基因表达差异显著性分析，具体如下表所示。

  若按照以上标准筛选得到的差异基因过少（低于100），很有可能导致后面的功能富集分析没有显著性结果，所以，我们会根据项目的具体情况，适当地降低筛选差异基因的阈值标准。若项目实验只关注某几个基因的表达情况（如基因敲除），不在意富集结果，从下面的差异分析表格中筛选关注的那几个基因即可。

  一般来说，如果一个基因在两组样品中的表达量差异达到两倍以上，我们认为这样的基因是具有表达差异的。为了判断两个样品之间的表达量差异究竟是由于各种误差导致的还是本质差异，我们需要对所有基因在这两个样本中的表达量数据进行假设检验。而转录组分析是针对成千上万个基因进行的，这样会导致假阳性的累积，基因数目越多，假设检验的假阳性累积程度会越高，所以引入padj对假设检验的P-value进行校正，从而控制假阳性的比例[3]。

  差异基因的筛选标准是非常重要的，我们给出的标准|log2(FoldChange)| > 1 & padj< 0.05是常用的经验值，在实际项目中可以根据情况灵活选择。例如，差异倍数可以选择1.5倍，也可以选择3倍，padj常用的阈值包括0.01、0.05、0.1等。若按照以上标准筛选得到的差异基因过少，很有可能导致后⾯的功能富集分析没有显著性结果。若项目实验只关注某几个基因的表达情况（如基因敲除），不在意富集结果，从下面的差异分析表格中筛选关注的那几个基因即可。反之，如果得到的差异基因数目过多，不利于后续目标基因的筛选，这个时候可使用更严格的阈值标准进行筛选，则可以使用更严格的阈值标准进行筛选。

2.6.1. 差异miRNA的筛选

  通过Deseq2进行差异分析，我们通常采用 |log2FC|>1 & padj < 0.05 进行差异miRNA的筛选，差异计算结果如下。

结果文件：

04.DE/Allgene_anno.xls表头

2.6.2. 差异miRNA的热图聚类

  将所有比较组的差异miRNA取并集之后作为差异miRNA集。两组以上的实验，可对差异miRNA集进行聚类分析，将表达模式相近的基因聚在一起。我们采用主流的层次聚类对基因的表达值进行聚类分析，对行（row）进行均一化处理（Z-score）。热图中表达模式相近的基因或样本会被聚集在一起，每个方格中的颜色反映的不是基因表达值，而是表达数据的行进行均一化处理后得到的数值（一般在-1到1之间），所以热图中的颜色只能横向比较（同一基因在不同样本中的表达情况），不能纵向比较（同一样本不同基因的表达情况）。结果文件中既有组间的聚类，也有样品间的聚类。结题报告展示了样品间的聚类，具体如下图所示。

图 2.6.2. 差异表达基因聚类热图

图中横坐标为样品名，纵坐标为差异miRNA归一化后的数值，颜色越红，表达量越高，越蓝，表达量越低。

结果文件：

差异miRNA的热图结果：04.DE/heatmap/

2.6.3. 差异miRNA的火山图分布

  火山图可直观展示每个比较组合的差异miRNA分布情况，如下图所示。图中横坐标表示基因在处理和对照两组中的表达倍数变化(log2FoldChange)，纵坐标表示基因在处理和对照两组中表达差异的显著性水平(-log10padj或-log10pvalue)。为上调基因用红色点表示，下调基因用蓝色点表示。

图 2.6.3. 差异miRNA火山图

图中横坐标为log2FoldChange值，纵坐标为-log10padj或-log10pvalue，蓝色的虚线表示差异基因筛选标准的阈值线

结果文件：

差异基因的火山图结果：04.DE/volcano/

2.7. 靶基因预测

  miRNA是一类在生物体内起到重要调控作用的的小片段非编码RNA。一般认为miRNA通过和mRNA的结合，可以抑制mRNA的表达，从而影响到Gene的表达。

all

down

up

图 2.7.1. TargetScan vs miRDB 的靶基因交集veen图