RIsearch2使用方法-预测RNA-RNA互作（sRNA的靶基因）

Y大宽

发布于 2018-09-30 10:41:57

2.5K0

发布于 2018-09-30 10:41:57

文章被收录于专栏：Y大宽Y大宽

背景

非编码RNA经常和其它RNAs形成配对（双链）发挥其作用。这些RNA-RNA相互作用都是建立在碱基互补配对的基础上，两个RNA序列之间的高度互补是这种相互作用的强有力预测基础。RIsearch2是RNA-RNA相互作用预测工具，可以在给定的query和target序列之间形成互补定位。使用基于suffix arrays的seed-and-extend框架，RIsearch2可以发现RNA-RNA相互作用关系，这种发现可以基于基因组或转录组。类似之前的 RIsearch，RIsearch2也使用基于di-nucleotides to approximate nearest-neighbor energy parameters的修正Smith-Waterman-Gotoh algorithm算法。然而，不是执行整个序列比对，RIsearch2关注种子区域的完美互补并且向两端延伸。用户定义的seed and extension constraints 使得 RIsearch2 可应用于所有类型的RNA-RNA相互作用预测。

1下载安装,设置环境变量

cd biosoft
wget https://rth.dk/resources/risearch/RIsearch-2.1.tar.gz
tar -xzvf RIsearch-2.1.tar.gz
cd RIsearch-2.1
less README
less INSTALL

加入环境变量

cp /home/kelly/biosoft/RIsearch-2.1/bin/risearch2.x /home/kelly/bin/.

查看用法

risearch2.x

================================ RIsearch v2.1 ===============================
================ Energy based RNA-RNA interaction predictions ================

Usage: risearch2.x [options]

  -h,         --help
                 show this message
------------------------------------------------------------------------------
--------------------------- SUFFIX ARRAY CREATION ----------------------------
  -c <FILE>,  --create=FILE (.fa or .fa.gz)
                 create suffix array for target sequence(s) together with
                 their reverse complements, FASTA format, use '-' for stdin
  -o <FILE>,  --output=FILE
                 save created suffix array to given index file path
------------------------------------------------------------------------------
--------------------------- INTERACTION PREDICTION ---------------------------
  -q <FILE>,  --query=FILE (.fa or .fa.gz)
                 FASTA file for query sequence(s), use '-' for stdin
  -i <FILE>,  --index=FILE
                 pregenerated suffix array file for target sequence(s)
  -s n:m/l,   --seed=n:m/l
                 set seed length (-s l = length only; -s n:m = full interval;
                 -s n:m/l = length in interval; default -s 6)
  -l <int>,   --extension=L
                 max extension length(L) on the seed (do DP for max this length
                 up- and downstream of seed) (default L=20)
  -e <float>, --energy=dG
                 set deltaG energy threshold (in kcal/mol) to filter predictions
                 (default=-20)
  -z mat,     --matrix=mat
                 set energy matrix to t99 or t04 (default) for RNA-RNA duplexes
  -d <int>,   --penalty=dP
                 per-nucleotide extension penalty given in dacal/mol
                 (recommended: 30, default: 0)
  -t <int>,   --threads=N
                 set maximum number of threads to use (default=1)
  -p,         --report_alignment
                 report predictions in detailed format
  -p2,        --report_alignment=2
                 report predictions in a simple format together with CIGAR-like
                 string for interaction structure
  -p3,        --report_alignment=3
                 report predictions in a simple format together with
                 binding site (3'->5'), flanking 5'end (3'->5') and
                 flanking 3'end (5'->3') sequences of the target
                 (required for post-processing of CRISPR off-target predictions)
  --noGUseed     consider G-U wobble pairs as mismatch within the seed
                 (only for locating seeds, energy model is not affected)
  --verbose      verbose output
------------------------------------------------------------------------------
---------------------------- EXPERIMENTAL OPTIONS ----------------------------
  -m c:p,     --mismatch=c:p
                 introduce mismatched seeds
                 Set the max num of mismatches (c) allowed in the seed and
                 min num of consecutive matches required at seed start/end (p)
                 ! These seeds will not overlap with perfect complementary seeds
                 (default -m 0:0  (no mismatch);
                 if you set c>0, please also set p>0 to avoid overlaps)
  -x <float>, --seed_energy=F
                 set energy per length threshold that filters seeds (default=0)

至此，已经完成安装，可以使用了。

2 如何使用：官方示例文件

和其它比对工具一样，RIsearch2也需要预先准备好的target 序列的index文件。所以先看RIsearch2如何产生index文件

2.1 为target序列产生index structure

目标序列只接受FASTA格式（或gzip压缩的FASTA文件），并且这些序列总是5'-3'格式。构建好的index会包括反向互补序列。下面这个命令展示如何产生目标序列(示例文件garget.fa为例)的的index文件

risearch2.x -c target.fa -o target.suf

这样会产生target.suf文件。当然，也可以支持多个文件

cat target*.fa.gz | risearch2.x -c - -o target.suf

2.2 RIsearch2进行相互作用预测

passing the query sequence(-q FILE),理解RIsearch2 输出 RIsearch2接受FASTA(或gzip压缩）格式序列文件（5'-3'），通过-q FILE option设定path。RIsearch2会为每个input文件（prefixed with risearch_seqID）产生gzipped格式结果文件。输入文件如果有重复ID，结果会被覆盖。
passing the index file as target(-i FILE) 前已述及，目标序列的file文件必须提前产生并且这个index 通过-i选项执行。下面是默认设置的最简单的运行方式

risearch2.x -q query.fa -i target.suf

这会产生risearch_query1.out.gz文件

seeds的最小长度（-s <int>） seeds定义为连续碱基配对的最大延伸和并且为种子确定任何序列互补性，无论是完美的还是接近完美的。所有种子由两个具有相同长度的序列组成，一个来自查询序列，另一个来自目标序列（或其反向互补）。如上所述，仅当种子由于位置限制或无效碱基配对而无法在任一端扩展时，种子才是最大的（除了摆动配对）。在RIsearch2中，种子最小长度| s | 由用户使用-s <int>选项设置; 排除比这更短的种子，并且不会急性下一步程序分析。例如，-s 7要求种子具有最少7个连续的碱基对而没有任何位置约束，这将在后面具体介绍。在默认设置下，它设置为-s 6。
预测相互作用的Energy threshold（能量阈值）(e<float>) RIsearch2仅在算法的第一步定位的种子区域周围预测相互作用。通过使用近似最近邻能量模型的动态编程方法在种子两端进行延伸来形成这些相互作用。对于每个定位的种子，发现具有最小自由能的相互作用。在某些情况下，这将仅仅是种子本身。但是，在报告产生之前可以通过应用能量阈值来过滤相互作用。此阈值由-e <float>选项设置，默认为-20，这意味着不会报告自由能变高于-20kcal / mol的预测。对于miRNA样相互作用，我们建议使用限制性较低的阈值，最高可达-10kcal / mol。下面是一个示例命令，能量阈值设置为7，最小种子长度设置为7。
动态编程表的大小限制（-l <int>） 这是种子延伸步骤最重要的设置。它决定了种子两端侧翼序列的长度，可以考虑延伸。此外，连同种子标准设置，对算法的运行时间影响最大。默认情况下，设置为20 。将其设置为0会使RIsearch2报告原始种子（假设它们足以满足能量阈值）。根据研究类型，建议使用10到30之间的值进行实际互作预测。但是，可以始终对small size的结果进行后续处理，以创建更长的互作预测。大小始终受序列边界限制，因此将其设置为非常高的值对应的是允许交互跨越整个查询序列。下面会详细解释。
报告不同输出格式的预测交互（-p，-p2，-p3） 如前所述，预测的query和target序列之间的互作结果会产生gzip压缩结果文件。预测的互作结果可以以默认的两种不同格式报告，这由用户设置参数-p或-p2控制。在默认设置中，交互以tsv格式报告。下面看实例

risearch2.x -c target.fa -o target.suf
risearch2.x -q query.fa -i target.suf -s 7 -e -13 -l 5
zcat risearch_query1.out.gz
zcat risearch_query2.out.gz

zcat risearch_query1.out.gz结果如下：

query1  9       21      ENST00000436685 451     463     +       -20.19

zcat risearch_query2.out.gz结果如下：

query2  9       22      ENST00000534717 281     294     -       -14.54
query2  9       22      ENST00000436685 156     169     -       -14.54

在上表中，

列分别表示查询序列ID
查询上的交互起始位置
查询上的交互结束位置
目标序列ID
目标上的交互起始位置
目标上的交互结束位置
相互作用的链
互作的自由能（以千卡/摩尔计） 当链为“—”时，代表在查询和反向互补靶序列之间发生实际预测相互作用，但是开始和结束位置总是指提供的靶序列上的位置，即正向链。 -p或-p2也将返回互作结构，这需要通过动态编程矩阵进行回溯。当-p被传递时，绑定站点额外可视化，并且每次相互作用输出需要四行：
（1）以5'到3'方向查询
（2）碱基对（GU之间“：”，CG和AU对之间“|”）
（3）3'到5'方向的目标序列
（4）预测的总结（如上所述的格式）。以下是与上面报告的相同交互的这种长输出格式的结果文件

risearch2.x -q query.fa -i target.suf -s 7 -e 13 -l 5 -p
zcat risearch_query1.out.gz

结果显示如下：

UUGAGAGGG
::||:||::
ggcuuucuu
query1  11      19      ENST00000534717 371     379     -       -5.31
UUGAGAGGG
::||:||::
ggcuuucuu
query1  11      19      ENST00000436685 246     254     -       -5.31
UUGAGAGGGCGA
:|||||:|  ||
gacucuucaccu
query1  11      22      ENST00000534717 110     121     -       -8.89
UUGAGAGGGCGA
:|||||:|  ||
gacucuucaccu
query1  11      22      ENST00000436685 10      21      -       -8.89
UGGGUUGAGAGGG
::|:::||: |||
gucuggcuugccc
query1  7       19      ENST00000534717 4       16      +       -8.46
UGGGUUG
:|::::|
gcuuggc
query1  7       13      ENST00000534717 35      41      -       0.54

对于更易于处理的格式（较小的文件大小和简单的解析），但仍提供有关结合位点本身的信息，可以用-p2，具有压缩互作结构的额外列添加到默认输出表中。它基本上是长格式第二行的记录，同时gap的信息使用字母编码如下：

P：规范碱基对
W：G-U摆动对
U：未配对
Q：查询中的凸起（查询中的核苷酸穿过靶中的间隙）
T：靶标中的凸起与输入序列一起，此信息足以重新创建以长格式（-p）给出的完整互作结果。使用-p2输出的示例(以前相同):

$ risearch2.x -c target.fa -o target.suf
$ risearch2.x -q query.fa -i target.suf -s 7 -e  -13 -l 5 -p2
$ zcat risearch_query1.out.gz
$ zcat risearch_query2.out.gz

zcat risearch_query1.out.gz结果如下：

query1  9       21      ENST00000436685 451     463     +       -20.19  PPUUPPPPPPPPP

zcat risearch_query2.out.gz结果如下：

query2  9       22      ENST00000534717 281     294     -       -14.54  PPUUPPPPWPPPWP
query2  9       22      ENST00000436685 156     169     -       -14.54  PPUUPPPPWPPPWP

更详细的-p3格式另外包括目标交互位点的序列以及其上游和下游区域的序列。当选择-p3时，程序比-p2格式输出预测增加三列：

目标结合位点（3'到5'方向）
5'末端上游（3'到5'方向）
目标序列相互作用位点3'末端（5'-3'方向）下游序列这是CRISPR-OFF网络服务器对CRISPR Cas9脱靶预测进行后续处理所需的格式，见下

$ risearch2.x -c target.fa -o target.suf
$ risearch2.x -q query.fa -i target.suf -s 7 -e -13 -l 5 -p3
$ zcat risearch_query1.out.gz
$ zcat risearch_query2.out.gz

zcat risearch_query1.out.gz结果如下：

query1  9       21      ENST00000436685 451     463     +       -20.19  PPUUPPPPPPPPP   ccuucucucccgc   ccgagaccucgacucuacuu    ugcagccugggaacuucagc

zcat risearch_query2.out.gz结果如下：

query2  9       22      ENST00000534717 281     294     -       -14.54  PPUUPPPPWPPPWP  acgccggagagagg  augggccugugacuacagua    gucgaucauagucuuccugc
query2  9       22      ENST00000436685 156     169     -       -14.54  PPUUPPPPWPPPWP  acgccggagagagg  augggccugugacuacagua    gucgaucauagucuuccugc