25、RNA结构高通量分析及三维建模方法解析

RNA结构高通量分析及三维建模方法解析

一、RNA结构高通量分析

1.1 PROBer处理与归一化

在RNA结构高通量分析中，FPKM（每百万读取中每千碱基的片段数）信息是重要的数据。.beta和.gamma文件分别包含每个核苷酸的反应性谱和转录酶噪声。可以从.beta文件计算归一化的反应性谱。

以流感病毒的SHAPE - seq数据为例，该病毒基因组由八个分段的病毒RNA（vRNA）组成。使用PROBer处理这些数据，并通过90% Winsorization对概率进行归一化。在第5段vRNA中预测到一个假结结构，其长度约为1500 nt，因此选择第5段vRNA作为高通量测序（HTS）分析的示例。

1.2 BUMHMM分析步骤

1.2.1 数据准备

在开始BUMHMM分析之前，需要将每个核苷酸的覆盖度、终止计数和终止率计算并保存到以制表符分隔的文件中。以下是示例数据：
| | C1 | C2 | T1 | T2 |
| — | — | — | — | — |
| | 50258 | 10914 | 22580 | 10548 |
| | 61362 | 13558 | 24949 | 11371 |

其中，“C1”和“C2”表示对照实验的重复1和2，“T1”和“T2”表示处理实验的重复1和2。

在R中，将这些数据存储在SummarizedExperiment对象中，具体操作如下：

library(BUMHMM)
library(Biostrings)
library(SummarizedExperiment)
library(IRanges)
counts <- read.delim("(coverage).txt")
counts <- as.matrix(counts)
dropoff <- read.delim("(dropoff).txt")
dropoff <- as.matrix(dropoff)
rate <- read.delim("(dropoff_rate).txt")
rate <- as.matrix(rate)
colData <- DataFrame(replicate=c("control","control","treatment","treatment"), row.names=c("C1","C2","T1","T2"))
se_exp <- SummarizedExperiment(assays=list(coverage=counts, dropoff_count=dropoff, dropoff_rate=rate), colData=colData)

1.2.2 核苷酸对选择与处理

接下来，进行核苷酸对的选择以计算对数比率，跨重复缩放终止率，并计算用于隐马尔可夫模型（HMM）分析的核苷酸位置延伸。

Nc <- Nt <- 2
t <- 1
nuclSelection <- selectNuclPos(se_exp, Nc, Nt, t)
assay(se_exp, "dropoff_rate") <- scaleDOR(se_exp, nuclSelection, Nc, Nt)
stretches <- computeStretches(se_exp,t)

1.2.3 偏差校正

进行覆盖度偏差和序列偏差的校正：

varStab <- stabiliseVariance(se_exp, nuclSelection, Nc, Nt)
LDR_C <- varStab$LDR_C
LDR_CT <- varStab$LDR_CT

hist(LDR_C, breaks = 30)
nuclNum <- 3
patterns <- nuclPerm(nuclNum)
sequence <- DNAString(scan("(referece sequence).txt", what=""))
nuclPosition <- findPatternPos(patterns, sequence, '+')

1.2.4 后验概率计算与输出

计算HMM的后验概率并保存输出文件：

nuclPosition <- list()
nuclPosition[[1]] <- 1:nchar(sequence)
posteriors <- computeProbs(LDR_C, LDR_CT, Nc, Nt, '+', nuclPosition, nuclSelection$analysedC, nuclSelection$analysedCT, stretches)
head(posteriors)
shifted_posteriors <- matrix(, nrow=dim(posteriors)[1], ncol=1)
shifted_posteriors[1:(length(shifted_posteriors) - 1)] <- posteriors[2:dim(posteriors)[1],2]
plot(shifted_posteriors)
write.csv(shifted_posteriors,"(output probability file name).txt", quote=F, row.names = F)

输出文件包含范围从0到1的概率。

1.3 reactIDR分析步骤

1.3.1 RT终止和覆盖度计数

使用docker镜像中的rt_end_counter对每个核苷酸的RT终止和覆盖度进行计数。假设所有bam文件都存储在/bam目录中，具体操作如下：

$ docker run --name rtecount -it carushi/rt_end_counter -i
$ mkdir bam
$ exit
$ docker start rtecount
$ docker cp (local path of bam files)/. rtecount:/bam/
$ docker exec -it rtecount /bin/bash
$ bash count_and_cov.sh bam/(filename).bam
(convert all bam files to bed)
$ cd bam
$ mkdir bed
$ mv *.bed bed
$ exit
$ docker cp rtecount:/bam/bed (local path to save bed files)
$ docker stop rtecount
$ docker rm rtecount

1.3.2 bed文件转换与合并

将bed文件进行转换和合并：

$ python ../script/bed_to_pars_format.py --offset -1 --fasta (reference fasta file).fa (filename)_l15q0filt_ctss.bed
(convert all ctss bed files)
$ python ../script/bed_to_pars_format.py --offset -0 --fasta (reference fasta file).fa (filename)_l15q0filt_cov.bed
(convert all cov bed files)
$ python ../reactIDR/score_converter.py --merge --output (name)_ctss (filename treatment rep1)_l15q0filt_ctss.bed.tab,(filename treatment rep2)_l15q0filt_ctss.bed.tab (filename control rep1)_l15q0filt_ctss.bed.tab,(filename control rep2)_l15q0filt_ctss.bed.tab
$ python ../reactIDR/score_converter.py --merge --output (name)_cov (filename treatment rep1)_l15q0filt_cov.bed.tab,(filename treatment rep2)_l15q0filt_cov.bed.tab (filename control rep1)_l15q0filt_cov.bed.tab,(filename control rep2)_l15q0filt_cov.bed.tab

1.3.3 参数评估（可选）

可以从预测的二级结构评估参数。使用RNA二级结构预测程序（如RNA fold和CentroidFold）预测RNA二级结构，并获得点括号格式的文件。在该文件所在目录运行以下脚本：

$ python ../reactIDR/IDR_hmm.py --train --time 10 --core 4 --grid --output_param (filename)_train.param.txt --output (filename).csv --ref (reference predicted secondary structure file).fa --case (filename)_ctss_case.tab --cont (filename)_ctss_cont.tab

1.3.4 概率计算

计算概率，如果上一步未创建训练参数文件，则- -param和- -ref选项不是必需的：

$ python ../reactIDR/IDR_hmm.py --test --time 10 --core 4 --param (filename)_train.param.txt --output (filename).csv --ref (reference predicted secondary structure file).fa --case (filename)_ctss_case.tab --cont (filename)_ctss_cont.tab

输出文件为test_(filename).csv，“IDR - HMM - final (transcript name)+case”项表示经过IDR和HMM后的后验概率。

1.3.5 结果可视化

reactIDR包包含用于可视化输出csv文件的Python程序：

$ python ../reactIDR/plot_bargraph.py --window 10 --ignore --output (output file name) (input csv file 1).csv (input csv file 2).csv

该脚本输出每个转录本的PDF文件。

1.4 注意事项

• 对于非生物信息学家，Python版本控制可能很复杂，建议使用Homebrew安装Python 3，并修改.bash_profile中的路径。
• 如果使用scipy时出现错误，检查系统中安装的scipy版本，可以使用以下命令安装scipy 1.2.2：

$ pip install scipy==1.2.2

• 在转录组范围分析中，根据使用的软件，需要使用不同的参考序列。PROBer和reactIDR可以一次计算每个转录本的概率，而BUMHMM一次读取和输出一个序列。
• 运行BUMHMM时，需要分别计算每个核苷酸的覆盖度、终止计数和终止率。覆盖度可以使用igvtools计算，终止计数可以从sam文件计算，终止率可以从覆盖度和终止计数计算。
• 可以使用R进行归一化，当使用90% Winsorization进行归一化时，在R中计算第5和第95百分位数，并基于这些值对概率进行归一化。
• 绘制图表需要安装seaborn和pandas Python包，可以使用以下命令安装：

$ pip install seaborn pandas

二、RNA 3D建模

2.1 引言

非编码RNA分子具有多种细胞功能，从催化作用到小分子检测再到翻译本身。它们通过采用复杂的三维折叠来执行这些功能。随着测序技术的发展，越来越多的RNA分子需要研究，实验方法和结构预测方法都在不断发展。

FARFAR2算法能够对相当复杂的RNA结构进行接近天然结构的预测，包括自动选择最终候选模型和估计模型准确性。可以使用公开可用的网络服务器（https://rosie.rosettacommons.org/farfar2）进行RNA建模。

2.2 FARFAR2 ROSIE服务器接口

2.2.1 基本接口

基本接口允许用户仅提供RNA结构预测任务中最常见的两个数据：序列和点括号二级结构。以甜菜西部黄化病毒（BWYV）的假结 - 1移码元件为例，用户可以指定RNA序列，大小写均可，链边界用逗号分隔。对于BWYV移码元件片段（PDB代码：1L2X），序列输入为：

gcgcggcaccguccgcggaacaaacgg

2.2.2 高级接口

高级接口允许用户指定更多选项，几乎可以指定影响Rosetta的FARFAR2算法命令行执行的每个参数。用户可以创建账户以获得更高优先级和电子邮件通知，也可以作为访客提交任务。如果涉及敏感数据，用户可以将任务设为私有。

2.3 示例分析

• 简单模型假结 ：使用基本接口对甜菜西部黄化病毒的移码元件中的简单模型假结进行建模。
• RNA - Puzzle 20复制 ：使用高级接口对宏基因组扭结姐妹核酶（RNA - Puzzle 20）进行建模。
• 包含实验数据的建模 ：
• 用MOHCA - seq实验的低分辨率约束对c - di - GMP核糖开关进行建模。
• 用1H NMR化学位移对串联GA基序进行建模。

通过这些方法和工具，研究人员可以更深入地了解RNA的结构和功能，为相关领域的研究提供有力支持。

graph LR
    A[RNA结构高通量分析] --> B[PROBer处理]
    A --> C[BUMHMM分析]
    A --> D[reactIDR分析]
    C --> C1[数据准备]
    C --> C2[核苷酸对选择与处理]
    C --> C3[偏差校正]
    C --> C4[后验概率计算与输出]
    D --> D1[RT终止和覆盖度计数]
    D --> D2[bed文件转换与合并]
    D --> D3[参数评估（可选）]
    D --> D4[概率计算]
    D --> D5[结果可视化]
    E[RNA 3D建模] --> F[FARFAR2 ROSIE服务器]
    F --> G[基本接口]
    F --> H[高级接口]
    G --> I[简单模型假结建模]
    H --> J[RNA - Puzzle 20复制]
    H --> K[包含实验数据的建模]

以上内容展示了RNA结构高通量分析和3D建模的详细方法和步骤，希望对相关研究人员有所帮助。

三、RNA结构高通量分析与3D建模的关联及应用拓展

3.1 高通量分析为3D建模提供基础数据

RNA结构高通量分析得到的结果，如核苷酸的反应性谱、后验概率等数据，能够为RNA 3D建模提供重要的约束信息。例如，在使用FARFAR2进行3D建模时，如果结合高通量分析得到的低分辨率约束，像MOHCA - seq实验数据，能够更准确地构建RNA的三维结构。这种关联可以用以下表格来展示：

高通量分析方法	提供的数据类型	对3D建模的作用
PROBer	归一化的反应性谱	辅助确定RNA的局部结构特征
BUMHMM	后验概率	作为结构建模的概率约束
reactIDR	经过IDR和HMM后的后验概率	帮助筛选更合理的3D模型

3.2 应用拓展案例

3.2.1 疾病相关RNA结构研究

在研究与疾病相关的RNA时，先通过高通量分析确定RNA的结构特征，再利用3D建模构建其三维结构。以流感病毒的RNA为例，通过PROBer、BUMHMM和reactIDR等高通量分析方法，得到RNA的反应性和概率信息，然后使用FARFAR2服务器进行3D建模。这样可以深入了解病毒RNA的结构与功能，为开发抗病毒药物提供靶点。

3.2.2 新型RNA分子的功能预测

对于新发现的非编码RNA分子，利用高通量分析和3D建模相结合的方法，可以预测其功能。首先进行高通量分析，确定RNA的二级结构特征，然后使用FARFAR2进行3D建模，通过分析其三维结构，推测其可能参与的细胞过程，如催化作用、小分子结合等。

3.3 未来发展趋势

3.3.1 算法优化

随着数据量的增加和计算能力的提升，高通量分析和3D建模的算法将不断优化。例如，FARFAR2算法可能会进一步提高预测的准确性和效率，能够处理更复杂的RNA结构。

3.3.2 多组学数据整合

未来可能会将RNA结构高通量分析与其他组学数据（如转录组学、蛋白质组学等）进行整合，从而更全面地了解RNA在细胞中的功能和调控机制。

3.3.3 临床应用拓展

在临床诊断和治疗方面，RNA结构分析和建模的应用将不断拓展。例如，通过分析患者体内特定RNA的结构变化，进行疾病的早期诊断和个性化治疗。

graph LR
    A[RNA结构高通量分析] --> B[提供数据] --> C[RNA 3D建模]
    C --> D[疾病相关RNA研究] --> E[抗病毒药物开发]
    C --> F[新型RNA功能预测] --> G[细胞过程推测]
    H[未来发展趋势] --> I[算法优化]
    H --> J[多组学数据整合]
    H --> K[临床应用拓展]

四、总结

RNA结构高通量分析和3D建模是研究RNA结构与功能的重要手段。高通量分析方法（如PROBer、BUMHMM和reactIDR）能够提供RNA的结构特征和概率信息，而FARFAR2算法和其对应的网络服务器则为RNA的3D建模提供了便利。这两种方法相互关联，为研究人员深入了解RNA的结构和功能提供了有力支持。

在实际应用中，它们可以用于疾病相关RNA研究、新型RNA分子的功能预测等多个领域。随着技术的不断发展，未来这些方法将在算法优化、多组学数据整合和临床应用拓展等方面取得更大的进展，为生命科学和医学领域带来更多的突破。

研究人员在使用这些方法时，需要注意一些细节，如Python版本控制、软件的参数设置等。同时，结合不同的实验数据和工具，能够更好地利用这些方法，推动RNA研究的发展。

通过本文的介绍，希望能够帮助研究人员更全面地了解RNA结构高通量分析和3D建模的方法和应用，为相关研究提供参考。