前面一篇文章讲了宏RNA测序时去除RNA的原因以及实验端去除让RNA的主流方法。这一章主要讲讲分析段如何去除rRNA。

实验时选择去除rRNA处理后，会大大降低rRNA的比例。但针对于痕量核酸样本（拭子类，脑脊液类），尤其是珍贵样本，不推荐实验时加去除rRNA处理。处理后回收核酸量太少，可能造成后续建库失败的风险。

本文主要讲拿到数据后如何将rRNA这部分数据剔除出来，后续的数据继续分析。那如何去除测序数据中的rRNA呢？

一、数据与软件

1. 宏RNA测序的下机数据
2. 参考数据集（rRNA）
3. seqkit 软件
4. Bowtie2 软件
5. Samtools 软件

二、 具体步骤：

1. 宏RNA测序的下机数据

目前宏RNA的测序主流平台还是二代测序平台，包括Illumina、BGIseq测序平台，还有少量来自于Thermo fisher 测序平台Ion S5（部分省市疾控会使用此平台针对新冠进行靶向测序）。测序策略一般为PE100或PE150。

2. 参考数据集（rRNA）

以人类rRNA为例。
登录NCBI网站，在搜索框输入“Homo sapiens”，选择Genomes 会跳转到Genomes-NCBI Datasets，截止写文前（2022-2-15）人类基因组有898条记录（里面有GCA 和GCF的重复）。但其实选一个参考序列即可，无需所有的基因组都下载。每个基因组3GB，全部下载下来比较浪费空间。

图1| homo sapiens.png

图2|人类基因组序列

image.png

点击RefSeq下载序列,

homo下载.png

image.png

另一种方式是去NCBI FTP Site（https://ftp.ncbi.nlm.nih.gov/）库里去搜索。
上述方法已经获取人类基因组参考序列的ID号为：GCF_000001405.39，将ID按照三个字符分隔，分别为/GCF/000/001/405，整个链接为https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405。

image.png

GRCh37和GRCh38是常用的两个人类参考序列版本，也即对应着UCSC的hg19和hg38，以GRCh38为例，选择了最新的p13版本，“GCF_000001405.39_GRCh38.p13”。每个版本有多次

GCF_000001405.39_GRCh38.p13网页下包含基因组序列、基因组注释、蛋白序列、编码区cds序列、RNA序列等。用的比较多的是基因组和蛋白序列以及注释信息(1基因组序列和注释信息)。2.RNA序列有两个文件，一个GCF_000001405.39_GRCh38.p13_rna.fna.gz 113MB ，一个GCF_000001405.39_GRCh38.p13_rna_from_genomic.fna.gz 80MB文件，

先将两个RNA序列下载下来。可浏览器直接点击下载（NCBI速度比较大，且中间可能会中断，多次尝试下载，要有耐心，哈哈哈哈），也可在Linux系统下使用wget下载。

wget -c https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405/GCF_000001405.39_GRCh38.p13/GCF_000001405.39_GRCh38.p13_rna_from_genomic.fna.gz
wget -c https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405/GCF_000001405.39_GRCh38.p13/GCF_000001405.39_GRCh38.p13_rna.fna.gz

下载好序列后，解压缩，提取其中的rRNA 序列，用作后续去除rRNA的参考数据集。
如果使用 GCF_000001405.39_GRCh38.p13_rna.fna文件提取rRNA，使用以下程序提取

grep 'ribosomal RNA' GCF_000001405.39_GRCh38.p13_rna.fna  |grep -v 'mRNA' |grep -v 'non-coding RNA'|grep -v 'misc_RNA' >GCF_000001405.39_GRCh38.p13_rna.ribosomal_RNA.acc

image.png

grep 'gbkey=rRNA' GCF_000001405.39_GRCh38.p13_rna_from_genomic.fna >GCF_000001405.39_GRCh38.p13_rna_from_genomic.ribosomal_RNA.acc
#提取并按照product列排序
grep 'gbkey=rRNA' GCF_000001405.39_GRCh38.p13_rna_from_genomic.fna |awk -F '[' '{print $4}' |sort

#输出
$ grep 'gbkey=rRNA' GCF_000001405.39_GRCh38.p13_rna_from_genomic.fna |awk -F '[' '{print $4}' |sort
product=l-rRNA]
product=RNA, 18S pre-ribosomal N4]
product=RNA, 18S ribosomal RNA N1]
product=RNA, 18S ribosomal RNA N1]
product=RNA, 18S ribosomal RNA N2]
product=RNA, 18S ribosomal RNA N3]
product=RNA, 18S ribosomal RNA N3]
product=RNA, 18S ribosomal RNA N5]
product=RNA, 28S ribosomal RNA N1]
product=RNA, 28S ribosomal RNA N1]
product=RNA, 28S ribosomal RNA N2]
product=RNA, 28S ribosomal RNA N3]
product=RNA, 28S ribosomal RNA N3]
product=RNA, 28S ribosomal RNA N4]
product=RNA, 28S ribosomal RNA N5]
product=RNA, 45S pre-ribosomal N1]
product=RNA, 45S pre-ribosomal N1]
product=RNA, 45S pre-ribosomal N2]
product=RNA, 45S pre-ribosomal N3]
product=RNA, 45S pre-ribosomal N3]
product=RNA, 45S pre-ribosomal N4]
product=RNA, 45S pre-ribosomal N5]
product=RNA, 5.8S ribosomal RNA N1]
product=RNA, 5.8S ribosomal RNA N1]
product=RNA, 5.8S ribosomal RNA N2]
product=RNA, 5.8S ribosomal RNA N3]
product=RNA, 5.8S ribosomal RNA N3]
product=RNA, 5.8S ribosomal RNA N4]
product=RNA, 5.8S ribosomal RNA N5]
product=RNA, 5S ribosomal 1]
product=RNA, 5S ribosomal 1]
product=RNA, 5S ribosomal 10]
product=RNA, 5S ribosomal 10]
product=RNA, 5S ribosomal 11]
product=RNA, 5S ribosomal 11]
product=RNA, 5S ribosomal 12]
product=RNA, 5S ribosomal 12]
product=RNA, 5S ribosomal 13]
product=RNA, 5S ribosomal 13]
product=RNA, 5S ribosomal 14]
product=RNA, 5S ribosomal 14]
product=RNA, 5S ribosomal 15]
product=RNA, 5S ribosomal 15]
product=RNA, 5S ribosomal 16]
product=RNA, 5S ribosomal 16]
product=RNA, 5S ribosomal 17]
product=RNA, 5S ribosomal 17]
product=RNA, 5S ribosomal 2]
product=RNA, 5S ribosomal 2]
product=RNA, 5S ribosomal 3]
product=RNA, 5S ribosomal 3]
product=RNA, 5S ribosomal 4]
product=RNA, 5S ribosomal 4]
product=RNA, 5S ribosomal 5]
product=RNA, 5S ribosomal 5]
product=RNA, 5S ribosomal 6]
product=RNA, 5S ribosomal 6]
product=RNA, 5S ribosomal 7]
product=RNA, 5S ribosomal 7]
product=RNA, 5S ribosomal 8]
product=RNA, 5S ribosomal 8]
product=RNA, 5S ribosomal 9]
product=RNA, 5S ribosomal 9]
product=s-rRNA]

除了有多copy的rRNA之外（5.8S、18S、28S、45S分别有7个copy，5S有34个copy），还有两个线类体rRNA（文后会补充线类体相关背景知识），为s-rRNA和l-rRNA，分别对应线类体的12S和18S。

image.png

至此，得到人类的rRNA有5.8S、18S、28S、45S、5S和线类体的s-rRNA和l-rRNA。至于45S rRNA前体，长度约为14000nt，后续会在细胞核的核仁中剪切生成80S核糖体的28S rRNA、5.8S rRNA和18S rRNA。

3. 提取序列

上述步骤已经得到rRNA 序列的ID，根据ID 信息从fasta格式的GCF_000001405.39_GRCh38.p13_rna_from_genomic.fna文件中提取相应的序列。可以使用脚本语言，比如awk、awk、Perl等，也可以使用seqkit软件提取。
seqkit软件是一款序列梳理神器-小巧（6M大小），兼具统计、格式转换、长度筛选、质量值转换、翻译、反向互补、抽样、去重、滑窗、拆分等30项全能...。此处仅简要介绍怎么抽提序列。
seqkit网址：https://bioinf.shenwei.me/seqkit/
latest verison： releasev2.1.0：https://github.com/shenwei356/seqkit/releases/download/v2.1.0/seqkit_linux_amd64.tar.gz

tar -xvf seqkit_linux_amd64.tar.gz
解压缩后即为可执行文件。若为了使用方便可将软件路径加到环境变量里面。

第二步提取到的序列id进行格式化，只取id空格前的字符串，输出到文本。

awk '{print $1}' GCF_000001405.39_GRCh38.p13_rna_from_genomic.ribosomal_RNA.acc |sed 's/>//' >GCF_000001405.39_GRCh38.p13_rna_from_genomic.ribosomal_RNA.id

提取序列
seqkit grep -f GCF_000001405.39_GRCh38.p13_rna_from_genomic.ribosomal_RNA1.id GCF_000001405.39_GRCh38.p13_rna_from_genomic.fna > GCF_000001405.39_GRCh38.p13_rna_from_genomic_rRNA.fna

提取出序列后可选择去冗余，经alignment 多个copy的序列后，发现序列之间相似度很大，除部分Gap存在外，仅有少量SNP存在，如下表所示：

去冗余主要是为了减小参考数据集，加快筛查的速度。每条rRNA的copy之间差异度很小，故每个5.8S、5S、18S、28S、45S rRNA 随机选一条；加上2条线类体rRNA，作为去除人类rRNA的参考数据集；笔者也尝试了GRCh37版本中提取rRNA，步骤与此类似，copy数有差异（GRCh37 copy数少）。

4. bowtie2、samtools下载和参考数据集构建index

bowtie2(http://bowtie-bio.sourceforge.net/bowtie2/index.shtml) 是将测序的reads与长参考序列比对的工具。bowtie2使用FM索引（基于Burrows-Wheeler transform 或BWT）对基因组进行索引，以此来保持其占用较小内存。bowtie2支持间隔、局部和双端对其模式。可以同时使用多个处理器来极大的提升比对速度。

Samtools(http://www.htslib.org/)，是一个用于操作sam和bam文件的工具合集。能够实现二进制查看、格式转换、排序及合并等功能，结合sam格式中的flag、tag等信息，还可以完成对比对结果的统计汇总。

（1）**下载bowtie2和samtools **，可以使用conda安装，就是可能不是最新版本。如果使用最新版本，可使用以下命令下载：

#下载bowtie2
wget -c --no-check-certificate  https://jaist.dl.sourceforge.net/project/bowtie-bio/bowtie2/2.4.5/bowtie2-2.4.5-linux-x86_64.zip
#解压后即可用
unzip bowtie2-2.4.5-linux-x86_64.zip 
#下载samtools
wget -c  --no-check-certificate https://jaist.dl.sourceforge.net/project/samtools/samtools/1.15/samtools-1.15.tar.bz2
tar -jxvf samtools-1.15.tar.bz2
#下载的为源码，需编译
cd samtools-1.15/
./configure
make
./samtools --help

构建参考数据集的索引index

$ bowtie2-build Del_human_rRNA_refseq.fasta human_rRNA

建好的人rRNA索引

（2）使用bowtie2去除rRNA

bowtie2 --very-sensitive-local --no-unal -I 1 -X 1000 -p 6 -x rRNA -1 R1.fq.gz -2 R2.fq.gz --un-conc-gz sample_rRNAremoved.fq.gz 2>sample_Map2rRNAStat.xls | samtools view -S -b -o sample_Map2rRNA.bam -

实例：
bowtie2 --very-sensitive-local --no-unal -I 1 -X 1000 -p 6 -x ../DB/human_rRNA -1 6-1-1_L01_72_1.fq.gz -2 6-1-1_L01_72_2.fq.gz --un-conc-gz 6-1-1_L01.fq.gz 2> 6-1-1_L01_Map2rRNAStat.xls |samtools view -S -b -o 6-1-1_L01_Map2rRNA.bam -

Bowtie2参数设置：
--very-sensitive-local --no-unal -I 1 -X 1000：这几个参数不描述，可以看看bowtie2的帮助页；
-p：线程数。
-x：参考序列的前缀（就是bowtie2-build rRNA.fa rRNA中的第二个参数）。
-1，-2：分别接测序数据的R1和R2，没有R2的话，只需要-1。
--un-conc-gz：比对不上的序列，以此前缀输出，压缩格式。
--un-conc <path> ：不能map的reads，fastq格式，非压缩格式，供后续分析直接使用。
sample_Map2rRNAStat.xls是bowtie2最后给出的比对结果，包括比对率、唯一比对、多重比对率等。
后面可以直接通过管道命令“|”，接samtools view -S -b -o sample_Map2rRNA.bam - 直接输出比对到rRNA的bam文件，用作其他分析。也可以不用samtools，直接 通过 ">"输出成 sam 文件。

samtools参数介绍：
-o 输出文件的名称
-@ 指使用的线程数
-b 该参数设置输出BAM格式，默认下输出是SAM格式文件
-S 默认下输入是BAM文件，若是输入为SAM文件，则最好加该参数，否则有时会报错

bowtie2跑完之后，将得到几个文件：

image.png

*_1.fq.gz 和 *_2.fq.gz 是原始文件
*rRNAremoved的R1和R2的fq，这个就是我们最后想拿到的，去除了rRNA的序列，用于后续分析；
*Map2rRNA.bam,比对上rRNA的bam；
*Map2rRNAStat.xls是比对到rRNA的情况。

image.png

最后，得到 *rRNAremoved.fq.1.gz 和 *rRNAremoved.fq.2.gz，就是去除rRNA序列之后的病毒宏基因组数据。

该样本为人肛拭子样本。比对人的rRNA的占比为69.41%。出乎我的医意料。但是和我的前一篇文章中的《宏RNA测序去除宿主核糖体rRNA（一）-实验端去除》(https://www.jianshu.com/p/3bfb4fa934d1)中描述的 在这转录组RNA中，核糖体RNA（rRNA）越占RNA总量的80%，是含量最多的一类RNA的观点相吻合。此次选用的原始文件左右双端均为1.2GB，耗时17分钟。故在时间要求比较高的样本中去宿主基因组的步骤可以改为先去除rRNA，因为rRNA的库很小，能快速扫完，把数据量降下来后再进行后续分析，有助于整个分析环节所需时长的缩短。

线类体是细胞核外唯一含有DNA的细胞器，线类体DNA（mitochondrial DNA,mtDNA）。mtDNA全长16569bp，双链，闭环分子，外环为重链（H链），内环为轻链（L链）。其基因结构在哺乳动物中保守，共编码37个基因：13种编码蛋白的基因，22种tRNA，2种rRNA，无内含子。唯一的非编码区：1000bpde D-loop环，包含重链复制起始点、轻重链转录启动子、4个高度保守序列。
异质性：每个细胞中含有数百个线类体，每个线类体含有2-10个拷贝的mtDNA,因此每个细胞中含有数千个mtDNA。
线类体广泛存在于除红细胞以外的组织细胞中。脑和骨骼肌作为高能量代谢器官，其细胞中含有大量线类体，因而mtDNA 缺陷产生极其广泛的病原谱，最主要的表现为线类体性脑疾病。

image.png

参考文献：

1. https://zhuanlan.zhihu.com/p/329967121
2. 焦老师讲遗传系列之9：线粒体遗传
3. 人类线粒体基因组简介