sortmerna是一款检测文库中rRNA含量，并去除rRNA reads的软件。当怀疑文库中包含的rRNA比例较高时，这款软件可以派上用场。以下简单介绍它的安装与基本使用方法。
Installation
conda install sortmerna
需要吐槽的是可以查到的官方manual还是2.1版的，作者的github也停更了...而用conda 安装的sortmerna都已经是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/目录复制到指定的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的运行需要指定工作目录，默认是当前目录。在工作目录下，软件会生成以下目录
Default structure: WORKDIR/
idx/ (参考序列的index)
kvdb/ (以key-value形式存储的比对结果)
out/ (输出目录)
readb/ (预处理的reads/index)
注意：在每次运行前要保证kvdb/下是空的，否则程序会报错
参数解释：
--ref: 参考序列fasta文件，有一个参考fasta使用一次这个参数。
--index: 程序运行的index操作，不使用这个参数的话，程序会检测workdir/idx/下是否存在已经建好的index，没有的话就根据--ref提供的fasta文件index。
`--index 0` 则跳过index，但若在`workdir/idx/`下没检测到index则会终止程序，并提醒没有index
`--index 1` 则只进行index
`--index 2` 和默认不使用该参数一样
去除 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
参数解释：
--reads: 输入的序列文件(FASTA/FASTQ/FASTA.GZ/FASTQ.GZ)，如果是双端测序就用两次这个参数
--aligned: 比对上参考数据库的序列，以rRNA数据库为参考文件的话，输出的是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: 输出比对结果为fasta/fastq格式。
--out2: 将双端文件单独输出。
--paired_out: flag，表明输入的文件是双端测序文件，一端比对不上的话，把成对的reads都输出到non-aligned的文件中。要与--fastx一起用，与--paired_in相互排斥。
上面的命令输出的结果为：
$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我们可以关注其统计的比对结果
$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软件的简介到此结束，使用最大的感受是这个软件非常的费时。如果只是想大概检测文库中rRNA含量，可以下载相应的rRNA fasta文件，并使用现在流行的比对软件（如STAR、Bowtie2等）进行Indexing和比对即可。这样会快很多。以后有新的使用心得再做更新。
补充于：2021-12-02
实际上，更加方便的rRNA remove手段应该是将测序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/
