Dean's blog
0%

sortmerna-去除rRNA

Posted on
Symbols count in article: 4.4k Reading time ≈ 4 mins.

sortmerna是一款检测文库中rRNA含量，并去除rRNA reads的软件。当怀疑文库中包含的rRNA比例较高时，这款软件可以派上用场。以下简单介绍它的安装与基本使用方法。

Installation

使用conda安装sortmerna

1
conda install sortmerna

需要吐槽的是可以查到的官方manual还是2.1版的，作者的github也停更了…而用conda安装的sortmerna都已经是4.3.2了，其中有许多更改的地方

1
2
3
$ sortmerna --version
SortMeRNA version 4.3.2
Build Date: Apr  2 2021

在官方的流程中，使用了silva 的rRNA databases作为rRNA的参考序列。silva的rRNA databases包括了古生菌、原核生物和真核生物rRNA的保守序列。你可以在以下链接下载rRNA database的压缩包，包括以下文件

rRNA database:

Index rRNA database

为了避免程序运行时每次都index，我们可以先对参考的rRNA databases进行index。Index完成后，以后运行程序前将idx/目录复制到指定的workdir/目录下即可

1
2
3
4
5
6
7
8
9
10
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的运行需要指定工作目录，默认是当前目录。在工作目录下，软件会生成以下目录

1
2
3
4
5
Default structure: WORKDIR/
                              idx/   # 参考序列的index
                              kvdb/  # 以key-value形式存储的比对结果
                              out/   # 输出目录
                              readb/ # 预处理的reads/index

注意：在每次运行前要保证kvdb/下是空的，否则程序会报错

参数解释：
--ref: 参考序列fasta文件，有一个参考fasta使用一次这个参数。
--index: 程序运行的index操作，不使用这个参数的话，程序会检测workdir/idx/下是否存在已经建好的index，没有的话就根据--ref提供的fasta文件index。

1
2
3
`--index 0` # 则跳过index，但若在`workdir/idx/`下没检测到index则会终止程序，并提醒没有index
`--index 1` # 则只进行index
`--index 2` # 和默认不使用该参数一样

去除 rRNA

以下命令对双端测序文件进行rRNA去除，并保留双端测序中成对的序列到两个文件中。sortmerna会自动检测输入的文件类型，并保持输出的文件类型与输入一致，例如以下输入.fq.gz，输出也是.fq.gz。

1
2
3
4
5
6
7
8
9
10
11
12
13
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。这里可以包含输出的目录和前缀，以下是一些可行的例子

1
2
3
4
5
'--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相互排斥。

上面的命令输出的结果为：

1
2
$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我们可以关注其统计的比对结果

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
$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/

完。