sortmerna是一款檢測文庫中rRNA含量，并去除rRNA reads的軟件锅减。當(dāng)懷疑文庫中包含的rRNA比例較高時，這款軟件可以派上用場板壮。以下簡單介紹它的安裝與基本使用方法变秦。
Installation
conda install sortmerna
需要吐槽的是可以查到的官方manual還是2.1版的称龙，作者的github也停更了...而用conda 安裝的sortmerna都已經(jīng)是4.3.2了的烁，其中有許多更改的地方
$sortmerna --version
SortMeRNA version 4.3.2
Build Date: Apr 2 2021
在官方的流程中绎橘，使用了silva 的rRNA databases作為rRNA的參考序列。silva的rRNA databases包括了古生菌骑脱、原核生物和真核生物rRNA的保守序列菜枷。你可以在以下鏈接下載rRNA database的壓縮包，包括以下文件
rRNA database:
-
silva_ids_acc_tax.tar.gz
euk：Eukarya 真核生物
bac：Bacteria 原核生物
arc：Archaea 古生菌
Index rRNA database
為了避免程序運行時每次都index叁丧，我們可以先對參考的rRNA databases進行index啤誊。Index完成后，以后運行程序前將idx/目錄復(fù)制到指定的workdir/目錄下即可
sortmerna --index 1 --threads 32 \
--ref ~/sortmerna4.3/rRNA_databases/silva-bac-16s-id90.fasta \
--ref ~/sortmerna4.3/rRNA_databases/silva-bac-23s-id98.fasta \
--ref ~/sortmerna4.3/rRNA_databases/silva-arc-16s-id95.fasta \
--ref ~/sortmerna4.3/rRNA_databases/silva-arc-23s-id98.fasta \
--ref ~/sortmerna4.3/rRNA_databases/silva-euk-18s-id95.fasta \
--ref ~/sortmerna4.3/rRNA_databases/silva-euk-28s-id98.fasta \
--ref ~/sortmerna4.3/rRNA_databases/rfam-5s-database-id98.fasta \
--ref ~/sortmerna4.3/rRNA_databases/rfam-5.8s-database-id98.fasta \
--workdir ./
要注意的是sortmerna的運行需要指定工作目錄拥娄，默認是當(dāng)前目錄蚊锹。在工作目錄下，軟件會生成以下目錄
Default structure: WORKDIR/
idx/ (參考序列的index)
kvdb/ (以key-value形式存儲的比對結(jié)果)
out/ (輸出目錄)
readb/ (預(yù)處理的reads/index)
注意：在每次運行前要保證kvdb/下是空的稚瘾，否則程序會報錯
參數(shù)解釋：
--ref: 參考序列fasta文件牡昆，有一個參考fasta使用一次這個參數(shù)。
--index: 程序運行的index操作摊欠，不使用這個參數(shù)的話丢烘，程序會檢測workdir/idx/下是否存在已經(jīng)建好的index，沒有的話就根據(jù)--ref提供的fasta文件index些椒。
`--index 0` 則跳過index播瞳，但若在`workdir/idx/`下沒檢測到index則會終止程序，并提醒沒有index
`--index 1` 則只進行index
`--index 2` 和默認不使用該參數(shù)一樣
去除 rRNA
以下命令對雙端測序文件進行rRNA去除摊沉，并保留雙端測序中成對的序列到兩個文件中狐史。sortmerna會自動檢測輸入的文件類型痒给，并保持輸出的文件類型與輸入一致说墨，例如以下輸入.fq.gz，輸出也是.fq.gz苍柏。
sortmerna --index 0 --threads 48 \
--ref ~/sortmerna4.3/rRNA_databases/silva-bac-16s-id90.fasta \
--ref ~/sortmerna4.3/rRNA_databases/silva-bac-23s-id98.fasta \
--ref ~/sortmerna4.3/rRNA_databases/silva-arc-16s-id95.fasta \
--ref ~/sortmerna4.3/rRNA_databases/silva-arc-23s-id98.fasta \
--ref ~/sortmerna4.3/rRNA_databases/silva-euk-18s-id95.fasta \
--ref ~/sortmerna4.3/rRNA_databases/silva-euk-28s-id98.fasta \
--ref ~/sortmerna4.3/rRNA_databases/rfam-5s-database-id98.fasta \
--ref ~/sortmerna4.3/rRNA_databases/rfam-5.8s-database-id98.fasta \
--workdir ./ \
--reads NR1.fq.gz \
--reads NR2.fq.gz \
--aligned "NR/N_rRNA" --other "NR/N_clean" --paired_in --fastx --out2
參數(shù)解釋：
--reads: 輸入的序列文件(FASTA/FASTQ/FASTA.GZ/FASTQ.GZ)尼斧，如果是雙端測序就用兩次這個參數(shù)
--aligned: 比對上參考數(shù)據(jù)庫的序列，以rRNA數(shù)據(jù)庫為參考文件的話试吁，輸出的是rRNA棺棵。這里可以包含輸出的目錄和前綴，以下是一些可行的例子
'--aligned $MYDIR/dir_1/dir_2/1' -> $MYDIR/dir_1/dir_2/1.fasta
'--aligned dir_1/apfx' -> $PWD/dir_1/apfx.fasta
'--aligned dir_1/' -> $PWD/aligned.fasta
'--aligned apfx' -> $PWD/apfx.fasta
'--aligned (no argument)' -> WORKDIR/out/aligned.fasta
--other: 與參考文件比對不上的序列熄捍，在這里就是rRNA-free的序列烛恤。用法與--aligned一樣
--paired_in: flag，表明輸入的文件是雙端測序文件余耽，一端比對上就算成對的reads都比對上了缚柏。要與--fastx一起用，與--paired_out相互排斥碟贾。
--fastx: 輸出比對結(jié)果為fasta/fastq格式币喧。
--out2: 將雙端文件單獨輸出轨域。
--paired_out: flag，表明輸入的文件是雙端測序文件杀餐，一端比對不上的話干发，把成對的reads都輸出到non-aligned的文件中。要與--fastx一起用史翘，與--paired_in相互排斥枉长。
上面的命令輸出的結(jié)果為：
$ls NR/
N_clean_fwd.fq.gz N_clean_rev.fq.gz N_rRNA_fwd.fq.gz N_rRNA.log N_rRNA_rev.fq.gz
對于輸出的*rRNA.log我們可以關(guān)注其統(tǒng)計的比對結(jié)果
$cat NR/N_rRNA.log
...
Results:
Total reads passing E-value threshold = 181078446 (96.01) # 文庫中rRNA 含量
Total reads failing E-value threshold = 7531112 (3.99) # 文庫中非rRNA 含量
Minimum read length = 150
Maximum read length = 20
Mean read length = 9
Coverage by database:
~/sortmerna4.3/rRNA_databases/silva-bac-16s-id90.fasta 0.07
~/sortmerna4.3/rRNA_databases/silva-bac-23s-id98.fasta 8.82
~/sortmerna4.3/rRNA_databases/silva-arc-16s-id95.fasta 0.32
~/sortmerna4.3/rRNA_databases/silva-arc-23s-id98.fasta 4.92
~/sortmerna4.3/rRNA_databases/silva-euk-18s-id95.fasta 11.53
~/sortmerna4.3/rRNA_databases/silva-euk-28s-id98.fasta 70.22
~/sortmerna4.3/rRNA_databases/rfam-5s-database-id98.fasta 0.00
~/sortmerna4.3/rRNA_databases/rfam-5.8s-database-id98.fasta 0.00
其中 "Total reads passing E-value threshold = 181078446 (96.01)" 指的是輸入文件中的rRNA含量（參考序列是rRNA時），"Total reads failing E-value threshold = 7531112 (3.99)"則是非rRNA的含量琼讽。
實測如果有多個這樣的log文件也可以用multiqc生成報告搀暑。
對sortmerna軟件的簡介到此結(jié)束，使用最大的感受是這個軟件非常的費時跨琳。如果只是想大概檢測文庫中rRNA含量自点，可以下載相應(yīng)的rRNA fasta文件，并使用現(xiàn)在流行的比對軟件（如STAR脉让、Bowtie2等）進行Indexing和比對即可桂敛。這樣會快很多。以后有新的使用心得再做更新溅潜。
補充于：2021-12-02
實際上术唬，更加方便的rRNA remove手段應(yīng)該是將測序reads使用STAR、Bowtie2等方法比對到rRNA參考序列滚澜，保留unmapped reads即可粗仓。
完。
ref:
Filtering rRNA from RNAseq data: https://www.biostars.org/p/321714/
SILVA rRNA database: https://www.arb-silva.de/
https://github.com/biocore/sortmerna/
