RNA-seq（转录组测序生信分析）去除rRNA的方法

      虽然湿实验中会有去除rRNA的步骤，但测出的原始序列中仍不可避免地存在有rRNA序列，所以我们在做转录组上游分析的时候可以添加去除rRNA的步骤，百度看到过很多软件可以去除rRNA，但是我做过测试，这些现成的软件运用都失败了，没有一个能去除rRNA，所以我用了自己的方法来，在linux操作系统下，以人的参考基因组hg19（GRCH37）为例，具体如下：

      首先，在NCBI上下载对应参考基因组的RNA序列，下载链接如下https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405/GCF_000001405.25_GRCh37.p13/GCF_000001405.25_GRCh37.p13_rna_from_genomic.fna.gz

      

        然后，将rRNA挑选出来：

     cat GCF_000001405.25_GRCh37.p13_rna_from_genomic.fna | grep "^>" | grep "gbkey=rRNA"  |  awk  '{print $1}'| sed 's/>//g' > id.list

            seqkit grep -f id.list GCF_000001405.25_GRCh37.p13_rna_from_genomic.fna > rRNA.fa

   

       再利用hisat2构建rRNA的索引：

             hisat2-build -p 16 rRNA.fa rRNA

            参数解释：

                                -p:是线程数；rRNA.fa是下载的rRNA序列，rRNA是所有索引的前缀名。

     

最后，利用hisat2比对，将没比对上的输出，即去除了rRNA：

             hisat2 --summary-file ${id}_duprRNA.txt -p 18 --un-gz ${id}_duprRNA.fq.gz -1 ${id}.fq.gz -2 ${id}_R2_clean.fq.gz -x ${ref_dir}/rRNA -S ${id}.sam  

           参数解释： 

                                -p：线程数

                                -un-gz：将未比对的单端读段输出到指定的文件路径，并对输出文件进行 gzip 压缩；--un-conc-gz，表示将未一致比对的成对读段输出到指定的文件路径，并对输出文件进行 gzip 压缩。

                                -x ：是对应的rRNA参考基因组，书写方式是——对应路径到前缀名。

                                -S：生成的sam文件，这个可以不写，但若不写，会在终端直接输出很长很多的sam文件，虽然不影响结果，但个人觉得眼花，所以我写了这个参数和指定输出文件名。

        

        上述步骤已经达到了去除rRNA的效果，对于后期分析提高精确度。