NanoRMS: 预测纳米孔RNA修饰计量学

从直接RNA测序数据集中的每读信息预测和可视化RNA修饰计量学

目录

• 一般描述
• 1. RNA修饰位点的预测
• 2. 使用Nanopolish重画波形图估计RNA修饰计量学（不推荐）
• 3. 使用Tombo重画波形图估计RNA修饰计量学（推荐）
• 4. 可视化个别位点的每读电流强度

依赖项和版本

• 引用
• 联系方式

一般描述

NanoRMS通过将读取分类为修饰/未修饰，基于它们的每读特征（电流强度、停留时间或追踪），来预测修饰计量学。

NanoRMS可以以单模式（1个样本）或配对模式（2个样本）运行。

NanoRMS可以运行无监督（例如K均值、聚合聚类、GMM）和有监督的机器学习算法（例如KNN、随机森林）。后者将需要成对样本，其中一个条件是敲除。

NanoRMS可以从Nanopolish重画波形图读取或从Tombo重画波形图读取中预测计量学。后者是推荐选项。

它的作用是什么？预测高度修饰和低度修饰RNA的计量学

NanoRMS可以使用无监督（K均值）或有监督（KNN）分类算法进行计量学预测。我们使用LC-MS/MS确认了伪尿嘧啶水平的合成分子来说明其定量能力。

NanoRMS将TRACE（TR）作为特征纳入预测每读修饰（从而计量学）。我们发现，将TRACE纳入极大地改善了RNA修饰计量学的预测。总体而言，我们发现最佳的特征组合是TRACE + 信号强度。

NanoRMS的计量学预测已在不同百分比修饰的RNA上的伪尿嘧啶和2-O-甲基化上进行了基准测试。

1. RNA修饰位点的预测

1.1. 使用Epinano-RMS提取基础呼叫特征

从参考fasta创建字典文件

请注意，我们使用的是专门为nanoRMS开发的EpiNano的修改版本（版本1.2）。EpiNano-RMS包括有关每种碱基（C_freq、A_freq、G_freq、T_freq）的不匹配比例的信息，以及总体不匹配频率。

首先，我们需要从参考fasta文件创建一个字典文件

java -jar epinano_RMS/picard.jar CreateSequenceDictionary REFERENCE=reference.fasta OUTPUT= reference.fasta.dict

运行Epinano-RMS

需要python3

python3 epinano_RMS/epinano_rms.py -R <reference_file> -b <bam_file> -s epinano_RMS/sam2tsv

使用测试数据的示例：

python3 epinano_RMS/epinano_rms.py -R test_data/yeast_rRNA_ref -b test_data/wt_sorted.bam -s epinano_RMS/sam2tsv

1.2. 预测RNA修饰

a) 单样本RNA修饰预测（即“de novo”预测）

单样本“de novo”RNA修饰预测已针对线粒体rRNA中的伪尿嘧啶RNA修饰进行了测试。使用基于CMC的探针测序（Nano-CMCseq）验证了新预测位点，所有3个生物重复中预测的2个位点均得到验证。

“de novo”伪尿嘧啶修饰位点的RNA修饰预测依赖于识别伪尿嘧啶基础呼叫错误的“签名”，这使我们能够在单个样本中de novo预测RNA修饰，只要修饰的计量学足够高（即能够与直接RNA测序的背景基础呼叫错误区分开来）。具体来说，伪尿嘧啶在修饰位置上引起强烈的不匹配签名，主要表现为C到U的不匹配（见下图）。

请注意，图5a和方法部分所示的阈值有错别字（但论文中的所有计算都是使用正确的阈值完成的）。正确的阈值如下：不匹配频率阈值：0.137，C频率阈值：0.578。

一般用法：

Rscript --vanilla Pseudou_prediction_singlecondition.R [options] -f <epinano_file1> (-s <epinano_file2> -t <epinano_file3>)

使用测试数据的示例（在线粒体核糖体RNA上预测伪尿嘧啶位点）：

Rscript --vanilla Pseudou_prediction_singlecondition.R -f WT_rRNA_Epinano.csv -s sn34KO_rRNA_Epinano.csv -t sn36KO_rRNA_Epinano.csv

b) 配对样本RNA修饰预测（即基于“差异错误”的预测）

伪尿嘧啶并不总是在高计量学中存在（例如rRNAs），但也可以以低计量学存在（例如在mRNAs中）。请注意，由于背景纳米孔“错误”与伪尿嘧啶引起的“错误”过于相似，因此不能使用“de novo”模式准确预测mRNA中的伪尿嘧啶。

对于这种情况，我们可以通过识别两个条件之间显示伪尿嘧啶差异错误签名的位点来预测差异性伪尿嘧啶化位点，如下所示。这种成对比较可以用于WT-KO，或两个条件之间（例如正常-热应激）。下图中，展示了使用此脚本识别的一些热响应位点的示例。差异错误基础预测可以应用于任何类型的RNA（mRNA、snoRNAs、snRNAs、rRNA等）。

对于转录组映射读取（来自一个链的读取）

一般用法：

Rscript --vanilla Pseudou_prediction_pairedcondition_transcript.R [options] -f <epinano_file1> -s <epinano_file2>

使用测试数据的示例：

Rscript --vanilla Pseudou_prediction_pairedcondition_transcript.R -f WT_ncRNA_Normal_Rep1_Epinano.csv -s WT_ncRNA_HeatShock_Rep1_Epinano.csv

对于基因组映射读取（来自两个链的读取）

我们首先需要处理Epinano输出和GTF文件，将它们转换为BED文件。此外，我们将使用Bedtools的相交功能将两者相交。

将Epinano输出转换为BED

一般用法：

./Epinano_to_BED.sh <epinano_file1>

使用测试数据的示例：

./Epinano_to_BED.sh WT_mRNA_Normal_Epinano.csv

将GTF（仅限CDS）输出转换为BED

一般用法：

Rscript --vanilla GTF_to_BED.R <GTF_File>

使用测试数据的示例：

Rscript --vanilla GTF_to_BED.R Saccer64.gtf

相交Epinano_BED文件和GTF_BED文件

所需工具：bedtools intersect

一般用法：

bedtools intersect -a <Epinano_Bed> -b <GTF_Bed> -wa -wb > output.bed

使用测试数据的示例：

bedtools intersect -a WT_mRNA_Normal_Epinano.csv.bed -b Saccer3.bed -wa -wb > WT_mRNA_Normal_Epinano_final.bed

一般用法：

Rscript --vanilla Pseudou_prediction_pairedcondition_genome.R [options] -f <epinano_file1> -s <epinano_file2>

使用测试数据的示例：

Rscript --vanilla Pseudou_prediction_pairedcondition_genome.R -f WT_mRNA_Normal_Epinano.csv -s WT_mRNA_HeatShock_Epinano.csv

2. 使用Nanopolish重画波形图估计RNA修饰计量学

这个版本已经过时。如果您仍想使用它，可以在这里找到详细信息和代码

3. 使用Tombo重画波形图估计RNA修饰计量学

要使用这个版本，您可以在这里找到安装详细信息

首先下载测试数据

这个数据集在per_read目录中有更详细的描述。

(cd per_read && wget https://public-docs.crg.es/enovoa/public/lpryszcz/src/nanoRMS/per_read/guppy3.0.3.hac -q --show-progress -r -c -nc -np -nH --cut-dirs=6 --reject="index.html*")

从所有样本中检索每读特征。

ref=per_read/guppy3.0.3.hac/Saccharomyces_cerevisiae.R64-1-1_firstcolumn.ncrna.fa
per_read/get_features.py --rna -f $ref -t 6 -i per
_read/guppy3.0.3.hac/*WT??C/workspace/*.fast5

您的Fast5文件必须已经基础呼叫并包含FastQ条目，

移动和追踪表（这可以通过h5ls -r batch0.fast5 | less进行检查）。

注意，当通过MinKNOW执行基础呼叫时，默认情况下不存储移动和追踪表。

在这种情况下，您需要重新基础呼叫您的Fast5文件

指定guppy_basecaller的--fast5_out参数。

估计两个样本之间的修饰频率差异

注意，您需要提供可能被修饰的候选位置。这些之前已经确定过——请参阅上面的第1.2节。预测RNA修饰。所以我们这里只是从现有的候选文件生成BED文件。

# 准备 BED - 这不再是需要的步骤！

f=per_read/results/predictions_ncRNA_WT30C_WT45C.tsv.gz
zgrep -v X.Ref $f |awk -F'\t' 'BEGIN {OFS = FS} {print $1,$2-1,$2,".",100,"+"}' > $f.bed

# 计算对照组和样本之间的修饰频率差异

per_read/get_freq.py -f $ref -b $f.bed -o $f.bed.tsv.gz -1 per_read/guppy3.0.3.hac/WT30C/workspace/.fast5.bam -2 per_read/guppy3.0.3.hac/WT45C/workspace/.fast5.bam

注意，候选位置文件（-b $f.bed）必须仅引用一些读取对齐的参考位置——

否则get_freq.py将报告有关低覆盖率的警告。

您可以使用--mincov更改所需的最小读取数，但至少为5，由于KNN的要求。

请注意，KMEANS不准确分配计量学变化的方向性，而KNN则准确（因为KMEANS随机分配一个聚类为“修饰”，另一个为“未修饰”。因此，要知道KMEANS计量学预测的变化方向，您需要从该给定位置的不匹配错误的方向性推断。如果您不关心变化的方向性，而只是关心变化的大小，您可以直接取预测计量学变化的绝对值。

4. 可视化个别位点的每读电流强度

4.1. 数据预处理

首先，通过折叠所有对于同一位置的多个观察结果，生成一个折叠的Nanopolish事件对齐输出。

python3 per_read_mean.py <event_align_file>

使用测试数据的示例：

python3 per_read_mean.py test_data/data1_eventalign_output.txt

其次，创建以感兴趣位置为中心的每读电流强度的15-mer窗口

上一步生成的Nanopolish事件对齐的输出在此脚本中用作输入。

Rscript --vanilla nanopolish_window.R positions_file <input_table> <label>

使用测试数据的示例：

Rscript --vanilla nanopolish_window.R test_data/positions test_data/data1_eventalign_output.txt_processed_perpos_mean.csv data1

4.2. 可视化电流强度信息：

4.2.1 密度图

Rscript --vanilla density_nanopolish.R <window_file1> <window_file2> <window_file3(optional)> <window_file4(optional)>

使用测试数据的示例：

Rscript --vanilla nanopolish_density_plot.R test_data/sn34_window_file.tsv test_data/wt_window_file.tsv

4.2.2 以修饰位点为中心的平均电流强度图

Rscript --vanilla nanopolish_meanlineplot.R <window_file1> <window_file2> <window_file3(optional)> <window_file4(optional)>

使用测试数据的示例：

Rscript --vanilla nanopolish_meanlineplot.R test_data/sn34_window_file.tsv test_data/wt_window_file.tsv

4.2.3 以修饰位点为中心的每读电流强度图

Rscript --vanilla nanopolish_perreadlineplot.R <window_file1> <window_file2> <window_file3(optional)> <window_file4(optional)>

使用测试数据的示例：

Rscript --vanilla nanopolish_perreadlineplot.R test_data/sn34_window_file.tsv test_data/wt_window_file.tsv

4.2.4 来自每读15-mer电流强度数据的主成分分析图

Rscript --vanilla nanopolish_pca.R <window_file1.tsv> <window_file2.tsv> <window_file3.tsv(optional)> <window_file4.tsv(optional)>

使用测试数据的示例：

Rscript --vanilla nanopolish_pca.R test_data/sn34_window_file.tsv test_data/wt_window_file.tsv

依赖项和版本

如果您遇到无法将浮点数NAN转换为整数的错误，

确保将numpy降级到1.20以下的版本

pip install "numpy<1.20"