miRNA-seq小RNA高通量测序pipeline：从raw reads，鉴定已知miRNA-预测新miRNA，到表达矩阵【一】

该教程分为2部分，第2部分在：miRNA-seq小RNA高通量测序pipeline：从raw reads，鉴定已知miRNA-预测新miRNA，到表达矩阵【二】

因为之前碰到了一批小RNA测序的数据，所以很是琢磨了一番时间。把自己整理出来的心得记录一下吧，以后或许也还会有用。

miRNA-seq是一种建库方法和数据处理上都稍有特殊的一种测序方法。它的原理是首先提取短链RNA，直接在RNA两端连接接头后建库，然后对所有小RNA进行一起测序，通过数据处理提取出其中潜在的miRNA并进行表达量分析。由于小RNA的长度都很短，所以50bp长的单端测序已经足够可以获得全长序列。

但是值得一提的是，miRNA-seq的数据处理流程到目前都没有固定。由于小RNA文库中会存在一定比例的rRNA、tRNA、snRNA、snoRNA和mRNA降解片段等其他类型的短RNA，而且有的研究需要发现新的miRNA，而有的只需要基于已发现的miRNA，所以目前各家的处理流程（包括软件、reads过滤流程、参考基因组类型）也是差异很大。也有一些课题组在做一体化的miRNA-seq数据处理分析平台的开发。我看过一些workflow，算是总结了一个比较适中的处理办法吧。

本教程miRNA-seq的流程基本如下：

fastqc质控→cutadapt去接头→trimmomatic质量过滤→clean data再次质控→rfam过滤已注释其他小RNA→bowtie过滤潜在mRNA降解片段→mirdeep2识别已知miRNA并鉴定新miRNA

系统路径的构架如下图：

教程涉及的系统路径

miRNA-seq数据处理需要的软件：

fastQC、cutadapt、trimmomatic、blast+、bowtie、mirdeep2。

（1）fastQC、cutadapt和bowtie的安装：

它们都可以通过anaconda3直接一次性下载安装添加环境变量：

conda install fastqc
conda install cutadapt
conda bowtie

注：很多小RNA测序是使用Bowtie进行比对，因为他们更适合短片段（<50bp）的比对，且相比Bowtie2不支持gap比对。

（2）本地版blast+的安装：

在ncbi的ftp即可以下载安装包。解压后放到自己想要的路径，最终设置环境变量即可。标准操作。

wget ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/ncbi-blast-2.10.0+-x64-linux.tar.gz -P ~/Downloads
cd ~/Downloads && tar -zxvf ncbi-blast-2.10.0+-x64-linux.tar.gz
mv ncbi-blast-2.10.0+-x64-linux ~/tools/blast_2.10.0
echo 'export PATH=$PATH:/home/xuzihan/tools/blast_2.10.0/bin' >> ~/.bashrc
source ~/.bashrc

（3）trimmomatic的安装：

它是Java的jar文件，通过Java来调用。

wget http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.39.zip -P ~/Downloads
cd ~/Downloads && unzip Trimmomatic-0.39.zip
mv ~/Downloads/Trimmomatic-0.39 ~/tools/trimmomatic

（4）mirdeep2的安装：

mirdeep2的运行依赖于数个软件。所以它的安装步骤相对繁杂一些。虽然官网上已经写出了全部一次性安装的脚本，但是据说目前已经不适用，所以保险起见，我们手动安装。它的依赖软件还包括:

（1）bowtie：用于比对短reads到参考基因组上。bowtie已经通过anaconda安装好了。
（2）ViennaRNA：预测RNA的二级结构。
（3）SQUID library：辅助（4）的安装。个人尚不清楚它的具体作用。
（4）randfold：检测RNA的折叠自由能是否属于miRNA。
（5）PDF::API2：辅助miRdeep2输出结果pdf文件。

#all softwares are saved at /home/xuzihan/tools/ except for bowtie.
#we can directly download and install bowtie from bioconda.
conda install bowtie

#Download mirdeep2 package from github.
cd ~/tools
sudo apt-get install git
git clone https://github.com/rajewsky-lab/mirdeep2.git
echo 'export 
#Now a folder 'mirdeep2' has been created.
PATH=$PATH:/home/xuzihan/tools/mirdeep2/src' >> ~/.bashrc

#Download and install ViennaRNA software.
wget https://www.tbi.univie.ac.at/RNA/download/sourcecode/2_4_x/ViennaRNA-2.4.14.tar.gz -P ~/tools
tar -xzvf ViennaRNA-2.4.14.tar.gz
cd ./ViennaRNA-2.4.14 && ./configure
make
make install
echo 'export PATH=$PATH:/home/xuzihan/tools/ViennaRNA-2.4.14/src/bin' >> ~/.bashrc

#Download and install squid library.
cd ..
wget http://eddylab.org/software/squid/squid.tar.gz -P ~/tools
tar -xzvf squid.tar.gz
cd ./squid-1.9g && ./configure
make
make install

#Download and install randfold software.
cd ..
wget http://bioinformatics.psb.ugent.be/supplementary_data/erbon/nov2003/downloads/randfold-2.0.1.tar.gz -P ~/tools
tar -xzvf randfold-2.0.1.tar.gz
cd ./randfold-2.0.1
sudo apt-get install emacs25
emacs Makefile
#change line with INCLUDE=-I. to INCLUDE=-I. -I<your_path_to_squid-1.9g> -L<your_path_to_squid-1.9g>
#for me: INCLUDE=-I. -I/home/xuzihan/tools/squid-1.9g/ -L/home/xuzihan/tools/squid-1.9g/
#Then save it.
make
echo 'export PATH=$PATH:/home/xuzihan/tools/randfold-2.0.1' >> ~/.bashrc

#Download and install PDF generator.
cd ..
wget https://cpan.metacpan.org/authors/id/S/SS/SSIMMS/PDF-API2-2.037.tar.gz -P ~/tools
tar -xzvf PDF-API2-2.037.tar.gz
cd ./PDF-API2-2.037 && mkdir ../mirdeep2/lib
perl Makefile.PL INSTALL_BASE=/home/xuzihan/tools/mirdeep2 LIB=/home/xuzihan/tools/mirdeep2/lib
make
make test
make install
echo 'export PERL5LIB=PERL5LIB:/home/xuzihan/tools/mirdeep2/lib/perl5' >> ~/.bashrc

#complete the installation
cd ../mirdeep2
touch install_successful

#restart the terminal.
source ~/.bashrc

然后测试是否安装成功：

bowtie
RNAfold -h
randfold
make_html.pl

若每个命令都有相应的软件回应，无报错信息。则说明安装成功。

数据处理

1.下机原始数据质检

我的原始数据已经储存于~/miRNA/rawdata路径了。分别为con_1/2/3.fastq和treat_1/2/3.fastq。

我的数据

首先使用fastqc对原始数据进行质量控制，脚本代码如下，并将其储存在~/miRNA/scripts路径（raw.fastqc.sh）：

#!/bin/bash
if [ -d ~/miRNA/fastqc ]
then
    echo "fastqc directory existed."
else
    mkdir ~/miRNA/fastqc 
    echo "fastqc directory newly created."
fi
cd ~/miRNA/rawdata
all_miR_fastq=`ls *.fastq`
for each in ${all_miR_fastq}
do
fastqc $each -o ~/miRNA/fastqc
done

运行脚本：它会在miRNA目录下新建fastqc子目录，并储存6个样本的QC报告。

cd ~/miRNA/scripts && bash raw.fastqc.sh

接下来我们就获得了html格式的测序报告，部分重要结果如下图：

每个碱基位的平均质量分数

在测序数据中较高频率出现的序列以及可能来源

接头含量及可能来源

如图所示，本次小RNA测序数据的碱基质量不错，但由于小RNA测序建库的特殊性，且测序长度为50bp，常会测到3'端的接头序列，因此接下来我们需要进行接头的去除和质量的过滤。

2.进行质量过滤和接头去除

编写脚本：在miRNA总文件夹中建立trimmed_data路径用于存储过滤后数据，然后对~/miRNA/rawdata中的6个raw data.fastq使用cutadapt进行3'端接头的去除以及长度控制（去除3'端接头序列后，仅保留17-35nt的reads作为潜在的miRNA序列），然后使用trimmomatic进行碱基质量的控制（reads开头和末尾碱基质量需>20，并以滑窗内碱基平均质量为20进行滑窗质检，在碱基修剪后仅保留>17nt的reads）。命名为read_filtering_cutadapt_trimmomatic.sh，脚本储存于~/miRNA/scripts路径中：

#!/bin/bash
if [ -d ~/miRNA/trimmed_data ]
then
    echo "trimmed_data directory has existed."
else
    mkdir ~/miRNA/trimmed_data
    echo "trimmed_data directory is newly created."
fi
cd ~/miRNA/rawdata
all_fastq=`ls *.fastq`
for each in ${all_fastq}
do
    cutadapt -a $1 --minimum-length=17 --maximum-length=35 \
     -o ../trimmed_data/"${each%.*}_trimmed.fastq" $each
    java -jar ~/trimmomatic/trimmomatic-0.39.jar SE -threads 4 \
     -phred33 "../trimmed_data/${each%.*}_trimmed.fastq" \
     "../trimmed_data/${each%.*}_clean.fastq" \
     LEADING:20 TRAILING:20 SLIDINGWINDOW:4:20 MINLEN:17 
    rm "../trimmed_data/${each%.*}_trimmed.fastq"
done

运行脚本，输入要去除的3'端接头序列即可：

cd ~/miRNA/scripts && bash read_filtering_cutadapt_trimmomatic.sh TGGAATTCTCGGGTGCCAAGGAACTC

脚本运行中

随后我们可获得~/miRNA/trimmed_data目录下的con_1/2/3_clean.fastq和treat_1/2/3_clean.fastq。

3. 对clean数据进行二次QC

编写脚本：调用fastqc，将clean data的QC结果保存在~/miRNA/trimmed_data/fastqc路径中。
脚本命名为second_fastqc.sh，保存于~/miRNA/scripts。

#!/bin/bash
if [ -d ~/miRNA/trimmed_data/fastqc ]
then
    echo "fastqc directory has existed."
else
    mkdir ~/miRNA/trimmed_data/fastqc
fi
cd ~/miRNA/trimmed_data
for each in `ls *.fastq`
do
fastqc $each -o ./fastqc
done

运行脚本。

cd ~/miRNA/scripts && bash second_fastqc.sh

对报告进行简要的解读，可见：

碱基质量经过过滤后是很好的

clean reads长度的分布：可见在22bp处出现峰。22bp正是miRNA的平均长度

现在检测到的出现频率较高的序列已经不再是adapters/primers了

此前混杂在reads中的adapter序列已经全部清除

4.过滤已经在rfam和ensembl数据库中注释的rRNA、tRNA、sRNA、snRNA

由于小RNA测序得到的reads中同时存在其他许多小RNA的种类，因此如果可以提前将这些可能造成混杂的序列去除，则可以加快后续比对的速度，也减小假阳性。

rfam数据库全面记录了目前为止已经注释的小RNA序列。因此除开已经注释的miRNA序列，我们挑出库中人类的rRNA、tRNA、sRNA和snRNA（包含了snoRNA）对reads进行比对，将比对上的reads去除。则剩下的reads为miRNA的可能性更大。

此外为了进一步的减小假阳性，我打算把rfam的这些ncRNA和ensembl数据库中的ncRNA注释合并进行一起过滤。ensembl数据库的大名就不多介绍了，它每一版release也包括了ncRNA序列，其类别可以分为snRNA、rRNA、scRNA、sRNA、Mt_rRNA、snoRNA、miRNA、Mt_tRNA、scaRNA、lncRNA、ribozyme、misc_RNA。

创建~/rfam路径，在其中下载rfam数据库中的全部fasta文件并解压；
创建~/ensembl路径，下载ensembl数据库中的Homo_sapiens.GRCh38.ncrna.fa.gz（随基因组一起，目前最新发布版本是release100）并解压：

mkdir ~/rfam && cd ~/rfam
curl -C -O ftp://ftp.ebi.ac.uk/pub/databases/Rfam/14.2/fasta_files/RF0[0001-3125].fa.gz
gunzip *.gz

mkdir ~/ensembl && cd ~/ensembl
wget ftp://ftp.ensembl.org/pub/release-100/fasta/homo_sapiens/ncrna/Homo_sapiens.GRCh38.ncrna.fa.gz -P ~/ensembl
gunzip *.gz

rfam14.2版本中的文件编号从RF00001至RF03125，每个文件都代表了一个小RNA注释的种类。

下载fa文件一览

我们通过复制rfam数据库网页中的taxonomy信息可以获得库中rRNA、tRNA、sRNA和snRNA的对应RF文件编号（图中表格的第一列），然后通过对fasta文件中名字行（>xxx）中的物种信息可以挑出人类（Homo sapiens）的注释。

在rfam-search中选择相应的taxonomy，复制相应的表格，用nano在终端中粘贴保存，表格第一列可用于不同类型数据库的提取

通过复制信息，我们得到并储存了rfam.anno.rRNA.txt、rfam.anno.tRNA.txt、rfam.anno.sRNA.txt和rfam.anno.snRNA.txt在rfam路径中。

它们的第一列为相应的RF号。

rfam.anno.rRNA.txt的内容

而ensembl下载的数据中仅包含人类的ncRNA fasta数据。其姓名列>xxx的第五列包含gene_biotype，可以提取出相应类型的小RNA序列。

Homo_sapiens.GRCh38.ncrna.fa中序列特征概览

可以使用如下命令查询ensembl下载的ncRNA类别：

cat Homo_sapiens.GRCh38.ncrna.fa | sed -n '/^>/p' | awk '{print $5}' | uniq

可得到：

ensembl收录的ncRNA类型

编写一个脚本，实现从rfam对应特征表格的txt中获取第一列文件信息，整合成分类的fasta文件并最后得到人类的rfam tRNA、rRNA、snRNA fasta文件。然后将ensembl库中对应类型的小RNA序列也添加到fasta数据库中。随后，为了用于blastn比对序列，我们还需要在比对前使用makeblastdb进行核酸比对数据库的建立。脚本命名为rfam.sRNA.hg.db.establish.sh，储存于~/miRNA/scripts：

#!/bin/bash
cd ~/rfam
snrna_no=`awk '{print $1}' rfam.anno.snRNA.txt|sed 's/$/&\.fa/g'`
cat $snrna_no|sed -n '/Homo sapiens/{p;n;p}' > rfam.hg.snrna.fa
rrna_no=`awk '{print $1}' rfam.anno.rRNA.txt|sed 's/$/&\.fa/g'`
cat $rrna_no|sed -n '/Homo sapiens/{p;n;p}' > rfam.hg.rrna.fa
trna_no=`awk '{print $1}' rfam.anno.tRNA.txt|sed 's/$/&\.fa/g'`
cat $trna_no|sed -n '/Homo sapiens/{p;n;p}' > rfam.hg.trna.fa
srna_no=`awk '{print $1}' rfam.anno.sRNA.txt|sed 's/$/&\.fa/g'`
cat $srna_no|sed -n '/Homo sapiens/{p;n;p}' > rfam.hg.srna.fa
cd ~/ensembl
cat Homo_sapiens.GRCh38.ncrna.fa | sed -n '/snRNA/{p;n;p}' >>  ~/rfam/rfam.hg.snrna.fa
cat Homo_sapiens.GRCh38.ncrna.fa | sed -n '/snoRNA/{p;n;p}' >>  ~/rfam/rfam.hg.snrna.fa
cat Homo_sapiens.GRCh38.ncrna.fa | sed -n '/rRNA/{p;n;p}' >>  ~/rfam/rfam.hg.rrna.fa
cat Homo_sapiens.GRCh38.ncrna.fa | sed -n '/tRNA/{p;n;p}' >>  ~/rfam/rfam.hg.trna.fa
cat Homo_sapiens.GRCh38.ncrna.fa | sed -n '/sRNA/{p;n;p}' >>  ~/rfam/rfam.hg.srna.fa
mkdir ~/miRNA/blastdb cd cd ~/miRNA/blastdb
makeblastdb -in ~/rfam/rfam.hg.trna.fa -dbtype nucl -title rfam.hg.trna -out rfam.hg.trna
makeblastdb -in ~/rfam/rfam.hg.rrna.fa -dbtype nucl -title rfam.hg.rrna -out rfam.hg.rrna
makeblastdb -in ~/rfam/rfam.hg.srna.fa -dbtype nucl -title rfam.hg.srna -out rfam.hg.srna
makeblastdb -in ~/rfam/rfam.hg.snrna.fa -dbtype nucl -title rfam.hg.snrna -out rfam.hg.snrna

运行：

cd ~/miRNA/scripts && bash rfam.sRNA.hg.db.establish.sh

建立blastdb中

数据库建立完毕，我们接下来即可利用blastdb对reads进行比对了。

但！由于样本中本身含有相同的小RNA序列，以及建库时的PCR扩增，测序获得的reads中含有大量的相同序列。如果把相同的所有reads简并到1条，在比对时将可以避免重复比对，大大缩短时间。

mirdeep2软件中的collapse_reads_md.pl脚本可以完成这一功能，输入reads.fasta和样本的三字符代号后，可将重复reads计数，以”xnumber“的形式加至序列名行（>开头的那一行）。（eg: >hsa_0_x328236，代表这个序列在hsa这个样本中出现了328236次）

我们编写如下脚本以获得所有样本的collapsed_reads.fasta文件。首先将储存于~/miRNA/trimmed_data的各样本clean data.fastq文件转化为不含质量信息的fasta文件，并将序列本身作为序列名字（eg: >AGCT），然后输入collapse_reads_md.pl进行折叠，以样本组别前两个字母（co/tr）+生物学重复编号（1/2/3）作为三字符样本编号。
脚本命名为get.collapsed.reads.sh，储存于~/miRNA/scripts路径：

#!/bin/bash
cd ~/miRNA/trimmed_data
for each in *.fastq
do
    series=`echo $each | grep -o "^[a-z]*\_[1-3]"`
    prefix=`echo $series | grep -o "^[a-z][a-z]"`
    suffix=`echo $series | grep -o "[0-9]"`
    awk '{if(NR%4 == 2){print">"$0"\n"$0}}' $each > ${each%.*}.fasta
    collapse_reads_md.pl ${each%.*}.fasta "${prefix}${suffix}" > ${each%.*}_collapsed.fasta
done

运行脚本：

cd ~/miRNA/scripts && bash get.collapsed.reads.sh

现在在trimmed_data路径获得了con_1/2/3_clean_collapsed.fasta和treat_1/2/3_clean_collapsed.fasta。我们比较一下初始的fastq文件和collapsed_reads.fasta文件：

注意名称列的suffix (>样本编号_序列编号_x出现次数)，和压缩了不要太多的文件大小

接下来使用该脚本进行collpased reads和之前建立的rfam人类数据库的比对（blastn命令），并从collapsed_reads.fasta中将这些比对上的序列去除。比对的结果保存在新建立的~/miRNA/blastresult路径中，比对结束后会出现con_1/2/3_clean_collapsed_rfam_filtered.fasta和treat_1/2/3_clean_collapsed_rfam_filtered.fasta，此外每个样本含有一个clean_collapsed_blast_result文件夹含有更详细的信息。
将该脚本命名为remove.rfam.annotated.sh，储存于~/miRNA/scripts路径：

#!/bin/bash
if [ -d ~/miRNA/blastresult ]
then
    echo "the directory for blast_result has existed."
else
    mkdir ~/miRNA/blastresult 
    echo "the directory for blast_result is newly created."
fi
cd ~/miRNA/trimmed_data
all_collapsed=`ls *collapsed.fasta`
cd ~/miRNA/blastresult
for each in $all_collapsed
    do
    mkdir ~/miRNA/blastresult/${each%.*}_blast_result && cd ~/miRNA/blastresult/${each%.*}_blast_result
    blastn -query ~/miRNA/trimmed_data/$each -out ${each%.*}_srna.blast -db ~/miRNA/blastdb/rfam.hg.srna -outfmt 6 -evalue 0.01
    blastn -query ~/miRNA/trimmed_data/$each -out ${each%.*}_rrna.blast -db ~/miRNA/blastdb/rfam.hg.rrna -outfmt 6 -evalue 0.01
    blastn -query ~/miRNA/trimmed_data/$each -out ${each%.*}_trna.blast -db ~/miRNA/blastdb/rfam.hg.trna -outfmt 6 -evalue 0.01
    blastn -query ~/miRNA/trimmed_data/$each -out ${each%.*}_snrna.blast -db ~/miRNA/blastdb/rfam.hg.snrna -outfmt 6 -evalue 0.01
    echo "blastn on $each against rfam rRNA/tRNA/snRNA completed."
    cat *.blast|awk '{print $1}'|sort|uniq > annotated.txt
    sed -n '/^>/p' ~/miRNA/trimmed_data/$each|sed 's/^>//g'|sort > all_seq.txt
    sort all_seq.txt annotated.txt annotated.txt|uniq -u > unannotated.txt
    touch ../${each%.*}_rfam_filtered.fasta
    for each_unanno_uniq in `cat unannotated.txt`
    do
        sed -n '/'$each_unanno_uniq'/''{p;n;p}' ~/miRNA/trimmed_data/$each >> ../${each%.*}_rfam_filtered.fasta
    done
done

但是这个脚本中最后使用单线程for循环+sed命令搜索未注释fasta序列的速度实在太慢。我的一个样品需要处理9小时（......）。而且我发现直接用blastn的过滤效果并不是太好，过滤的比例很低。

所以我改用了blastn-short，它更适用于<30bp的短序列的比对。此外，用一个python脚本替换了直接用shell搜寻未注释的fasta序列。将下面的脚本命名为remove.rfam.annotated.v2.sh，储存在~/miRNA/scripts路径：

#!/bin/bash
if [ -d ~/miRNA/blastresult ]
then
    echo "the directory for blast_result has existed."
else
    mkdir ~/miRNA/blastresult
    echo "the directory for blast_result is newly created."
fi
cd ~/miRNA/trimmed_data
all_collapsed=`ls *collapsed.fasta`
cd ~/miRNA/blastresult
for each in $all_collapsed
    do
    mkdir ~/miRNA/blastresult/${each%.*}_blast_result && cd ~/miRNA/blastresult\
/${each%.*}_blast_result
    blastn -task blastn-short -query ~/miRNA/trimmed_data/$each -out ${each%.*}_srna.blast \
 -db ~/miRNA/blastdb/rfam.hg.srna -outfmt 6 -evalue 0.01 -num_threads 1
    blastn -task blastn-short -query ~/miRNA/trimmed_data/$each -out ${each%.*}_rrna.blast \
-db ~/miRNA/blastdb/rfam.hg.rrna -outfmt 6 -evalue 0.01 -num_threads 1
    blastn -task blastn-short -query ~/miRNA/trimmed_data/$each -out ${each%.*}_trna.blast \
-db ~/miRNA/blastdb/rfam.hg.trna -outfmt 6 -evalue 0.01 -num_threads 1
    blastn -task blastn-short -query ~/miRNA/trimmed_data/$each -out ${each%.*}_snrna.blast \
-db ~/miRNA/blastdb/rfam.hg.snrna -outfmt 6 -evalue 0.01 -num_threads 1
echo "blastn on $each against rfam rRNA/tRNA/snRNA completed."
    cat *.blast|awk '{print $1}'|sort|uniq > annotated.txt
    sed -n '/^>/p' ~/miRNA/trimmed_data/$each|sed 's/^>//g'|sort > all_seq.txt
    sort all_seq.txt annotated.txt annotated.txt|uniq -u > unannotated.txt
    python3 ~/miRNA/scripts/find_fasta.py ~/miRNA/trimmed_data/$each \
unannotated.txt ../${each%.*}_rfam_filtered.fasta
done

最后使用未注释名称匹配总fasta文件的步骤直接借用了其他大佬的python脚本。这个脚本命名为find_fasta.py，也储存在~/miRNA/scripts路径：根据ID从FASTA文件中批量提取序列【Python脚本】

# -*- coding: utf-8 -*-
"""
@author: gyw
@Date:  Wed Feb 28 2019
@E-mail: willgyw@126.com
@Description: A script to extract sequences from fasta file.
"""
import sys 

def usage():
    print('Usage: python3 script.py [fasta_file] [idlist_file] [outfile_name]')


def main():
    outf = open(sys.argv[3],'w')
    dict = {}
    with open(sys.argv[1], 'r') as fastaf:
        for line in fastaf:
            if line.startswith('>'):
                name = line.strip().split()[0][1:]
                dict[name] = ''
            else:
                dict[name] += line.replace('\n','')
     
    with open(sys.argv[2],'r') as listf:
        for row in listf:
            row = row.strip()
            for key in dict.keys():
                if key == row:
                    outf.write('>' + key+ '\n')
                    outf.write(dict[key] + '\n')
    outf.close()


try:
    main()
except IndexError:
    usage()

运行该脚本：

cd ~/miRNA/scripts && bash remove.rfam.annotated.v2.sh

py脚本使用了字典查询功能，所以速度快了很多，单线程一个样本处理约2小时。期间就让服务器跑吧，不用管它...（如果有更好效率更高的方式希望大家可以和我多交流！）

脚本运行结束后，我们在~/miRNA/blastresult路径中获得了6个逃过了rfam注释的collapsed reads.fasta。

Optional：查看过滤的数据情况

如果我们想查看本次blast比对后被过滤掉的reads数、保留的reads数以及比例，可以使用如下的脚本完成。脚本命名为count_rfam_rate.sh，并保存于~/miRNA/scripts路径。

计算方法：
对unique reads：直接查看annotated.txt和all_seq.txt的行数。
对total reads：计算annotated.txt和all_seq.txt中每行中的"xnnnn"的数字和。

#!/bin/bash
cd ~/miRNA/blastresult
dir=`ls *_blast_result | grep -o "^[a-z]*\_[1-3]\_[a-z]*\_[a-z]*\_[a-z]*\_[a-z]*"`
for each_dir in $dir
do
    cd ~/miRNA/blastresult/$each_dir
    filtered_uniq=`cat annotated.txt | wc -l`
    total_uniq=`cat all_seq.txt | wc -l`
    filtered_ratio=`echo "scale=4;$filtered_uniq/$total_uniq" | bc`
    uniq_ratio=`echo "scale=4;($total_uniq-$filtered_uniq)/$total_uniq" | bc`
    sample_name=`echo "$each_dir" | grep -o '^[a-z]*\_[1-3]'`
    echo -e "**For unique reads of $sample_name **""\n"\
    "total unique reads = $total_uniq""\n"\
    "rfam filtered unique reads = $filtered_uniq""\n"\
    "rate of alignment = $filtered_ratio""\n"\
    "pass rate = $uniq_ratio"
    total_read=0
    total_read_char=`cat all_seq.txt | grep -o '[0-9]*$'`
    for each_num in $total_read_char
    do
         total_read=$((total_read+${each_num}))
    done
    total_filtered_read=0
    total_filtered_read_char=`cat annotated.txt | grep -o '[0-9]*$'`
    for each_num in $total_filtered_read_char
    do
        total_filtered_read=$((total_filtered_read+${each_num}))
    done
    total_filtered_ratio=`echo "scale=4;$total_filtered_read/$total_read" | bc`
    total_ratio=`echo "scale=4;($total_read-$total_filtered_read)/$total_read" | bc`
    echo -e "**For total reads of $sample_name **""\n"\
    "total reads = $total_read""\n"\
    "rfam filtered total reads = $total_filtered_read""\n"\
    "rate of alignment = $total_filtered_ratio""\n"\
    "pass rate = $total_ratio"
done