RUISBIO ｜ 睿锶生物

RIC-Seq生信分析报告

1.项目简介

1.1 样本信息

1.2 实验原理及流程

RIC-seq(RNA in situ conformation sequencing)技术通过全局性、高分辨率的原位分析，揭示细胞内RNA分子的三维空间互作网络及其生物学功能（Cai et al., 2020）。RIC-seq通过生物素标记和邻近连接技术，生成全基因组范围的RNA-RNA互作图谱，定位RNA分子在细胞核内的空间分布及动态组装模式。RIC-seq可以用于发现非编码RNA（如增强子RNA、长链非编码RNA）的靶标及作用机制，揭示其在基因表达调控中的新功能，例如通过染色质环化调控增强子-启动子通讯。该技术主要包括以下步骤：
(1)​​细胞交联和通透处理​​：使用甲醛固定细胞，破坏细胞膜和核膜的完整性，使细胞处于可渗透状态;
(2)​​RNA随机切割和去磷酸化​​：使用微球菌核酸酶（MNase）随机切割RNA，并对3'突出端进行去磷酸化处理;
(3)​​RNA 3'端标记​​：使用生物素化的胞苷三磷酸（pCp-biotin）对RNA 3'端进行标记;
(4)​​原位邻近连接​​：在非变性条件下，对邻近的RNA片段进行连接，形成嵌合体。
RIC-Seq 实验原理如下图所示：

1.3 实验流程

RIC-Seq实验主要步骤有六步:
(1)细胞交联与通透处理​,固定细胞内RNA-蛋白质相互作用及RNA空间构象。首先用1%甲醛处理细胞10分钟，使RNA与结合蛋白及邻近RNA交联，之后加入0.125 M甘氨酸终止交联，用含NP-40的缓冲液通透细胞膜，暴露胞内RNA。
(2)RNA随机切割与标记​,生成短RNA片段并标记3'末端，便于后续连接。首先用微球菌核酸酶（MNase）随机切割RNA至约90 nt片段，然后使用碱性磷酸酶（FastAP）去除RNA 3'端磷酸基团，之后使用T4 RNA连接酶催化生物素化的胞苷三磷酸（pCp-biotin）连接至RNA 3'端。
(3)通过超声打断DNA成200-1000bp的片段。首先在温和条件下使用T4 RNA连接酶催化邻近RNA片段的5'端与生物素标记的3'端连接，之后用碱性磷酸酶去除未连接的pCp-biotin，T4激酶修复5'磷酸基团。
(4)利用抗体将与目的蛋白相结合的DNA片段沉淀下来。首先通过生物素-链霉亲和素结合捕获嵌合体RNA，之后通过反转录与扩增获取RIC-seq测序文库。
(5)高通量测序。
(6)生信分析，获取全基因组范围内RNA-RNA互作片段。

1.4 分析流程

获得测序原始数据（raw data）后，首先过滤原始数据中低质量序列，获得高质量的测序数据（clean data）。先将clean data与rRNA进行比对以去除rRNA数据，之后再与参考基因组比对，获取包含互作信息的reads。根据基因注释信息将互作reads区分为跨不同基因的分子间互作与同一基因内的分子内互作​。最后使用蒙特卡洛模拟筛选出显著性的分子间互作片段(Cao et al., 2021)。​

RIC-Seq 生物信息学分析流程

2. 数据质控

我们交付的原始数据为fastq（简称fq）格式文件的压缩包，文件名后缀通常为 “.fq.gz”。交付数据前我们会计算每个压缩文件的md5值。在您拿到数据之后，请您先校>验每个压缩文件的md5值，Linux下可以在数据目录使用“md5sum -c <*md5.txt>”命令进行校验，Windows下可使用hashmyfiles等校验工具，如发现压缩文件md5值与附在数据文件目录下的md5文档中的不一致则说明文件可能在传输的过程中被损坏。数据文件大小为文件占用磁盘空间的大小，文件的大小通常与磁盘格式、压缩比例等因素有关，与测序数据量（碱基数）的多少无对应关系，因此对应PE测序的 read1和read2两个文件大小也可能不相同。

将高通量测序得到的原始图像数据经过Base Calling 转化为序列数据，即FASTQ格式，得到最原始的测序数据文件。FASTQ 格式文件可记录所测读段（read）的碱基及其质量分数。FASTQ 格式以测序读段为单位进行存储，每条读段占 4 行，第一行是序列标识（read ID）以及相关的描述信息，以“@” 开头；第二行即为碱基序列，长度由测序策略决定；第三行以“+”开头，后面是序列标示符、描述信息，或者什么也不加； 第四行是测序质量值（phred），与第二行一一对应，phred值以ASCII码标记，对应的 ASCII 值减去33，即为第二行对应碱基的测序质量值，示例如下：

@HWI-ST1276:71:C1162ACXX:1:1101:1208:2458 1:N:0:CGATGT
NAAGAACACGTTCGGTCACCTCAGCACACTTGTGAATGTCATGGGATCCAT
+
#55???BBBBB?BA@DEEFFCFFHHFFCFFHHHHHHHFAE0ECFFD/AEHH

测序错误率用e表示， 平台测得数据的碱基质量值用Qphred表示，则有：Qphred=-10log10(e)。软件中碱基识别正确率与Phred分值之间的简明对应关系见下表：

测序Reads的错误率往往会随着测序接近尾声而升高，这是由测序过程中化学试剂的消耗造成共有的特征。

2.1 原始数据质控

RIC-Seq实验基于第二代测序（NGS）平台完成，采用双端测序文库构建策略（插入片段~300 bp）。我们需要对原始测序数据进行质量评估与过滤，以确保后续分析的可靠性。首先，我们使用FastQC（version 0.12.1）(Andrews, 2010)对原始测序数据（raw data）进行全局质量分析，包括碱基质量分布（Phred score）、碱基组成平衡性（base content uniformity）、重复序列比例（duplication level）及GC含量偏差等指标，以全面评估测序质量。
我们使用Trimmomatic（version 0.39）(Bolger A. M. et al., 2014)对原始测序数据进行以下过滤操作。
接头序列去除：识别并切除双端reads中的接头序列；
修剪低质量碱基：修剪读段5'端质量低于25的碱基以及3'端质量低于20的碱基；
动态质量修剪：通过滑动窗口法（4 bp窗口步长）评估局部序列质量，当窗口平均质量值小于15时，执行3'端截断；
长度筛选：保留长度≥30 bp的reads。

在此之后我们根据reads的id信息去除PCR重复reads，再将去重reads使用Cutadapt(version 4.4)(Martin, 2011)去除低复杂度序列（如连续的AAAAA或CCCCC等低复杂度区域）以及过滤含高比例单一碱基的reads，得到过滤后数据clean data。

原始和过滤后质控结果请详见result/1.qc文件夹，raw为原始数据质控结果，clean为过滤后质控结果。

hela1_1hela1_2hela2_1hela2_2

❮

❯

图2.1 各个样本平均测序碱基质量分数，横坐标代表150 bp长度序列中各个位置，纵坐标为该位置平均的碱基质量值Q；盒形图中间的红线表示中位数（median value）；黄色部分代表四分位距（25-75%）；上下分割线代表 90%和 10%的上下临界值；蓝色的线代表碱基质量的平均值。

hela1_1hela1_2hela2_1hela2_2

❮

❯

图2.2 各个样本碱基平衡性，图中四条线代表A T C G在每个位置平均含量。理论上，A和T应该相等，G和C应该相等，且4种碱基平行且接近分布。正常情况下四种碱基的出现频率应该是接近的，而且没有位置差异。因此好的样本中四条线应该平行且接近。当部分位置碱基的比例出现 bias 时，即四条线波动较大时可能存在测序数据或者文库污染。如果所有位置的碱基比例一致的表现出bias 时，即四条线平行但分开，往往代表文库有 bias (建库过程或本身特点)，或者是测序中的系统误差。测序刚开始由于测序仪状态不稳定，在15bp之前很可能出现波动。

hela1_1hela1_2hela2_1hela2_2

❮

❯

图2.3 各个样本重复序列水平，测序深度越高，越容易产生一定程度的重复（duplication），这属于正常的现象。但如果duplication 的程度很高，就提示我们可能有 bias 的存在（如建库过程中由于 PCR 扩增引起的duplication）。横坐标为 reads 重复的次数，纵坐标为重复次数对应的 reads 占 unique reads 的比例，以unique reads 的总数作为 100%。这里，我们仅对文件前 2000000 个reads 进行统计：对长度小于75bp 的reads 将其截短为 50bp，用于统计重复。

2.2 过滤后数据质控

这里展示Fastp过滤后的数据质控结果，图片内容与上面raw data类似。

hela1_1.dedup.clean.pair.cuthela1_2.dedup.clean.pair.cuthela2_1.dedup.clean.pair.cuthela2_2.dedup.clean.pair.cut

❮

❯

图2.4 各个样本平均测序碱基质量分数，横坐标代表150 bp长度序列中各个位置，纵坐标为该位置平均的碱基质量值Q；盒形图中间的红线表示中位数（median value）；黄色部分代表四分位距（25-75%）；上下分割线代表 90%和 10%的上下临界值；蓝色的线代表碱基质量的平均值。

hela1_1.dedup.clean.pair.cuthela1_2.dedup.clean.pair.cuthela2_1.dedup.clean.pair.cuthela2_2.dedup.clean.pair.cut

❮

❯

图2.5 各个样本碱基平衡性，图中四条线代表A T C G在每个位置平均含量。理论上，A和T应该相等，G和C应该相等，且4种碱基平行且接近分布。正常情况下四种碱基的出现频率应该是接近的，而且没有位置差异。因此好的样本中四条线应该平行且接近。当部分位置碱基的比例出现 bias 时，即四条线在某些位置波动较大时，可能测序数据或者文库存在污染。当所有位置的碱基比例一致的表现出bias 时，即四条线平行但分开，往往代表文库有 bias (建库过程或本身特点)，或者是测序中的系统误差。一般测序的时候，刚开始测序仪状态不稳定，在15bp之前很可能出现波动。

hela1_1.dedup.clean.pair.cuthela1_2.dedup.clean.pair.cuthela2_1.dedup.clean.pair.cuthela2_2.dedup.clean.pair.cut

❮

❯

图2.6 各个样本重复序列水平，测序深度越高，越容易产生一定程度的重复（duplication），这属于正常的现象。但如果duplication 的程度很高，就提示我们可能有 bias 的存在（如建库过程中由于 PCR 扩增引起的duplication）。横坐标为 reads 重复的次数，纵坐标为重复次数对应的 reads 占 unique reads 的比例，以unique reads 的总数作为 100%。这里，我们仅对文件前 2000000 个reads 进行统计：对长度小于75bp 的reads 将其截短为 50bp，用于统计重复。

2.3 数据过滤结果统计

我们对数据过滤结果进行统计，如下表所示：

表 2.1数据过滤结果统计：
Sample：样品名称；
Raw_Total_Reads/Clean_Total_Reads：过滤前后样本总reads数量，单位为百万；
Raw_Total_Bases/Clean_Total_Bases：过滤前后样本总碱基数量，单位为百万；
Raw_GC_Content/Clean_GC_Content:过滤前后样本GC含量。

3. 比对参考基因组

我们将各样品过滤后的clean data的reads首先与物种rRNA进行比对以去除来自rRNA的序列。之后与参考基因组进行比对，获取reads在基因组上的位置信息。这里使用的软件是STAR(version 2.7.11b)(Dobin A. et al., 2013)。

3.1 去除rRNA情况

表 3.1 去除rRNA统计：
Sample：样品名称；
Input_reads：clean后reads总数；
Mapped_to_rRNA：比对到rRNA序列的reads数；
rRNA_percent：rRNA序列百分比；
After_rRNA：去除rRNA后剩余的reads数量

3.2 比对参考基因组情况

获得去除rRNA的reads后，将其与参考基因组比对，标记为唯一比对reads和嵌合体reads将作为可用reads，在后续分析中用于提取RNA-RNA相互作用信息。

表 3.2 比对参考基因组统计：
Sample：样品名称；
Input_reads：去除rRNA后reads总数；
Uniquely_mapped：唯一比对reads数；
Multi_mapped：比对到多个区域reads数；
Unmapped: 没有比对上基因组的reads数；
Chimeric：鉴定为嵌合reads数；
valid_percent：可用reads百分比，它等于唯一比对reads数加上嵌合reads数占总数的百分比。

3.3 比对可视化

客户可以结合物种参考基因组和注释文件使用IGV (Integrative Genomics Viewer) 浏览器对bam、bw、bed等文件进行可视化浏览。IGV浏览器使用方法可参考我们提供的使用说明文档IGV快速上手

4. 鉴定相互作用

获取可用的reads后，筛选出含有空隙（skipped bases on the reference，CIGAR标识符为N）的reads，根据空隙将reads划分为左右两个标签，之后根据参考基因组注释信息区分为分子内互作片段和分子间互作片段。如果存在重复样本，那么会将同组的reads结果进行合并。我们重点关注分子间互作片段，它反映了RNA-RNA相互作用。为了获取可信的互作片段，我们使用蒙特卡洛模拟随机生成大量模拟互作事件（1万次模拟），假设无真实互作时，互作reads的分布服从随机概率。计算模拟背景下互作的期望值和方差，之后比较实验数据与模拟结果差异的显著性，最后保留P值低于阈值（0.05）的显著性互作。本节结果请详见位于report/result/3.pair文件夹中。
名为*.significant.interMolecular.interaction.list.gz的压缩包是全部显著性分子间互作结果文件，它的“GeneA”和“GeneB”代表发生互作的两个RNA对应的转录本id。
名为*.hic的文件是用于Juicebox的可视化文件，如同Hi-C记录了全基因组DNA-DNA互作信息一样，这里使用.hic格式文件记录RNA之间的互作。
名为*.bedpe.gz的文件是压缩的互作片段坐标文件，解压后可以用文本编辑器打开，也可在IGV中查看。

hela_count_heatmap

❮

❯

图 4.1 互作基因类型统计：横轴坐标为发生互作的RNA类型，颜色代表log10转化后的互作片段总数。

5. RNA互作图谱

获取全基因组范围的RNA-RNA相互作用片段后，我们绘制全基因组以及各个染色体的RNA互作热图。网页中只展示部分染色体结果，全部染色体结果请详见位于report/result/3.heatmap文件夹中。

hela_genome

❮

❯

图 5.1 全基因组RNA互作图谱：颜色代表该区域相互作用片段数量。

hela_chr1hela_chr11hela_chr21

❮

❯

图 5.2 部分染色体RNA互作图谱：颜色代表该区域相互作用片段数量。

参考文献

Andrews S. FastQC: a quality control tool for high throughput sequence data.https://www.bioinformatics.babraham.ac.uk/projects/fastqc/, 2010.
Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014 Aug 1;30(15):2114-20.
Cai Z, Cao C, Ji L, Ye R, Wang D, Xia C, Wang S, Du Z, Hu N, Yu X, Chen J, Wang L, Yang X, He S, Xue Y. RIC-seq for global in situ profiling of RNA-RNA spatial interactions. Nature. 2020 Jun;582(7812):432-437.
Cao C, Cai Z, Ye R, Su R, Hu N, Zhao H, Xue Y. Global in situ profiling of RNA-RNA spatial interactions with RIC-seq. Nat Protoc. 2021 Jun;16(6):2916-2946.
Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, Batut P, Chaisson M, Gingeras TR. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013 Jan 1;29(1):15-21. doi: 10.1093/bioinformatics/bts635. Epub 2012 Oct 25. PMID: 23104886; PMCID: PMC3530905.
Martin, M., Cutadapt Removes Adapter Sequences from High-Throughput Sequencing Reads, EMBnet.Journal, vol. 17, no. 1, pp. 10–12, from https://journal.embnet.org/index.php/embnetjournal/article/view/200, May 2, 2011.

联系我们

官网：武汉睿锶生物技术有限公司 咨询热线：400-689-7068

邮箱：sales@ruisbio.com 地址：湖北省武汉市江夏区东湖健康城

报告文件结构