美文网首页便捷工具
sortmerna 去除rRNA

sortmerna 去除rRNA

作者: 尘世中一个迷途小书僮 | 来源:发表于2021-06-08 20:01 被阅读0次

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:

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/

相关文章

网友评论

    本文标题:sortmerna 去除rRNA

    本文链接:https://www.haomeiwen.com/subject/galxeltx.html