背景
非編碼RNA經(jīng)常和其它RNAs形成配對（雙鏈）發(fā)揮其作用。這些RNA-RNA相互作用都是建立在堿基互補配對的基礎上，兩個RNA序列之間的高度互補是這種相互作用的強有力預測基礎。RIsearch2是RNA-RNA相互作用預測工具，可以在給定的query和target序列之間形成互補定位。使用基于suffix arrays的seed-and-extend框架，RIsearch2可以發(fā)現(xiàn)RNA-RNA相互作用關系，這種發(fā)現(xiàn)可以基于基因組或轉錄組。類似之前的 RIsearch，RIsearch2也使用基于di-nucleotides to approximate nearest-neighbor energy parameters的修正Smith-Waterman-Gotoh algorithm算法。然而，不是執(zhí)行整個序列比對，RIsearch2關注種子區(qū)域的完美互補并且向兩端延伸。 用戶定義的seed and extension constraints 使得 RIsearch2 可應用于所有類型的RNA-RNA相互作用預測。
1下載安裝,設置環(huán)境變量
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
加入環(huán)境變量
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)
至此，已經(jīng)完成安裝，可以使用了。
2 如何使用：官方示例文件
和其它比對工具一樣，RIsearch2也需要預先準備好的target 序列的index文件。所以先看RIsearch2如何產(chǎn)生index文件
2.1 為target序列產(chǎn)生index structure
目標序列只接受FASTA格式（或gzip壓縮的FASTA文件），并且這些序列總是5'-3'格式。構建好的index會包括反向互補序列。下面這個命令展示如何產(chǎn)生目標序列(示例文件target.fa為例)的的index文件
risearch2.x -c target.fa -o target.suf
這樣會產(chǎn)生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）產(chǎn)生gzipped格式結果文件。輸入文件如果有重復ID，結果會被覆蓋。 -
passing the index file as target(-i FILE)
前已述及，目標序列的file文件必須提前產(chǎn)生并且這個index 通過-i選項執(zhí)行。下面是默認設置的最簡單的運行方式
risearch2.x -q query.fa -i target.suf
這會產(chǎn)生risearch_query1.out.gz文件
-
seeds的最小長度（
-s <int>）
seeds定義為連續(xù)堿基配對的最大延伸和并且為種子確定任何序列互補性，無論是完美的還是接近完美的。 所有種子由兩個具有相同長度的序列組成，一個來自查詢序列，另一個來自目標序列（或其反向互補）。如上所述，僅當種子由于位置限制或無效堿基配對而無法在任一端擴展時，種子才是最大的（除了擺動配對）。在RIsearch2中，種子最小長度| s | 由用戶使用-s <int>選項設置; 排除比這更短的種子，并且不會急性下一步程序分析。 例如，-s 7要求種子具有最少7個連續(xù)的堿基對而沒有任何位置約束，這將在后面具體介紹。 在默認設置下，它設置為-s 6。 -
預測相互作用的Energy threshold（能量閾值）(
e<float>)
RIsearch2僅在算法的第一步定位的種子區(qū)域周圍預測相互作用。 通過使用近似最近鄰能量模型的動態(tài)編程方法在種子兩端進行延伸來形成這些相互作用。 對于每個定位的種子，發(fā)現(xiàn)具有最小自由能的相互作用。 在某些情況下，這將僅僅是種子本身。 但是，在報告產(chǎn)生之前可以通過應用能量閾值來過濾相互作用。 此閾值由-e <float>選項設置，默認為-20，這意味著不會報告自由能變高于-20kcal / mol的預測。 對于miRNA樣相互作用，我們建議使用限制性較低的閾值，最高可達-10kcal / mol。 下面是一個示例命令，能量閾值設置為7，最小種子長度設置為7。 -
動態(tài)編程表的大小限制（
-l <int>）
這是種子延伸步驟最重要的設置。 它決定了種子兩端側翼序列的長度，可以考慮延伸。 此外，連同種子標準設置，對算法的運行時間影響最大。 默認情況下，設置為20 。將其設置為0會使RIsearch2報告原始種子（假設它們足以滿足能量閾值）。 根據(jù)研究類型，建議使用10到30之間的值進行實際互作預測。 但是，可以始終對small size的結果進行后續(xù)處理，以創(chuàng)建更長的互作預測。 大小始終受序列邊界限制，因此將其設置為非常高的值對應的是允許交互跨越整個查詢序列。下面會詳細解釋。 -
報告不同輸出格式的預測交互（
-p，-p2，-p3）
如前所述，預測的query和target序列之間的互作結果會產(chǎn)生gzip壓縮結果文件。 預測的互作結果可以以默認的兩種不同格式報告，這由用戶設置參數(shù)-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
- 目標上的互作起始位置
- 目標上的互作結束位置
- 相互作用的鏈
- 互作的自由能 （以千卡/摩爾計）
當鏈為“—”時，代表在查詢和反向互補靶序列之間發(fā)生實際預測相互作用，但是開始和結束位置總是指提供的靶序列上的位置，即正向鏈。
-p或-p2也將返回互作結構，這需要通過動態(tài)編程矩陣進行回溯。
當-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：規(guī)范堿基對
- W：G-U擺動對
- U：未配對
- Q：查詢中的凸起（查詢中的核苷酸穿過靶中的間隙）
- T：靶標中的凸起
與輸入序列一起，此信息足以重新創(chuàng)建以長格式（-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格式另外包括目標交互位點的序列以及其上游和下游區(qū)域的序列。 當選擇-p3時，程序比-p2格式輸出預測增加三列：
- 目標結合位點（3'到5'方向）
- 5'末端上游（3'到5'方向）
- 目標序列相互作用位點3'末端（5'-3'方向）下游序列
這是CRISPR-OFF網(wǎng)絡服務器對CRISPR Cas9脫靶預測進行后續(xù)處理所需的格式，見下
$ 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
