該教程分為2部分，第2部分在：miRNA-seq小RNA高通量測序pipeline：從raw reads巍举，鑒定已知miRNA-預(yù)測新miRNA脆荷，到表達(dá)矩陣【二】

因為之前碰到了一批小RNA測序的數(shù)據(jù)，所以很是琢磨了一番時間懊悯。把自己整理出來的心得記錄一下吧蜓谋，以后或許也還會有用。

miRNA-seq是一種建庫方法和數(shù)據(jù)處理上都稍有特殊的一種測序方法炭分。它的原理是首先提取短鏈RNA桃焕，直接在RNA兩端連接接頭后建庫，然后對所有小RNA進(jìn)行一起測序捧毛，通過數(shù)據(jù)處理提取出其中潛在的miRNA并進(jìn)行表達(dá)量分析观堂。由于小RNA的長度都很短让网，所以50bp長的單端測序已經(jīng)足夠可以獲得全長序列。

但是值得一提的是师痕，miRNA-seq的數(shù)據(jù)處理流程到目前都沒有固定溃睹。由于小RNA文庫中會存在一定比例的rRNA、tRNA胰坟、snRNA因篇、snoRNA和mRNA降解片段等其他類型的短RNA，而且有的研究需要發(fā)現(xiàn)新的miRNA笔横，而有的只需要基于已發(fā)現(xiàn)的miRNA竞滓，所以目前各家的處理流程（包括軟件、reads過濾流程吹缔、參考基因組類型）也是差異很大商佑。也有一些課題組在做一體化的miRNA-seq數(shù)據(jù)處理分析平臺的開發(fā)。我看過一些workflow涛菠，算是總結(jié)了一個比較適中的處理辦法吧莉御。

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

fastqc質(zhì)控→cutadapt去接頭→trimmomatic質(zhì)量過濾→clean data再次質(zhì)控→rfam過濾已注釋其他小RNA→bowtie過濾潛在mRNA降解片段→mirdeep2識別已知miRNA并鑒定新miRNA

系統(tǒng)路徑的構(gòu)架如下圖：

教程涉及的系統(tǒng)路徑

miRNA-seq數(shù)據(jù)處理需要的軟件：

fastQC、cutadapt俗冻、trimmomatic礁叔、blast+、bowtie迄薄、mirdeep2琅关。

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

它們都可以通過anaconda3直接一次性下載安裝添加環(huán)境變量：

conda install fastqc
conda install cutadapt
conda bowtie

注：很多小RNA測序是使用Bowtie進(jìn)行比對讥蔽，因為他們更適合短片段（<50bp）的比對涣易，且相比Bowtie2不支持gap比對。

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

在ncbi的ftp即可以下載安裝包冶伞。解壓后放到自己想要的路徑新症，最終設(shè)置環(huán)境變量即可。標(biāo)準(zhǔn)操作响禽。

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來調(diào)用。

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的運行依賴于數(shù)個軟件芋类。所以它的安裝步驟相對繁雜一些隆嗅。雖然官網(wǎng)上已經(jīng)寫出了全部一次性安裝的腳本，但是據(jù)說目前已經(jīng)不適用侯繁，所以保險起見胖喳，我們手動安裝。它的依賴軟件還包括:

（1）bowtie：用于比對短reads到參考基因組上贮竟。bowtie已經(jīng)通過anaconda安裝好了丽焊。
（2）ViennaRNA：預(yù)測RNA的二級結(jié)構(gòu)较剃。
（3）SQUID library：輔助（4）的安裝。個人尚不清楚它的具體作用粹懒。
（4）randfold：檢測RNA的折疊自由能是否屬于miRNA重付。
（5）PDF::API2：輔助miRdeep2輸出結(jié)果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

若每個命令都有相應(yīng)的軟件回應(yīng)凫乖，無報錯信息。則說明安裝成功弓颈。

數(shù)據(jù)處理

1.下機原始數(shù)據(jù)質(zhì)檢

我的原始數(shù)據(jù)已經(jīng)儲存于~/miRNA/rawdata路徑了帽芽。分別為con_1/2/3.fastq和treat_1/2/3.fastq。

我的數(shù)據(jù)

首先使用fastqc對原始數(shù)據(jù)進(jìn)行質(zhì)量控制翔冀，腳本代碼如下导街，并將其儲存在~/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格式的測序報告搬瑰，部分重要結(jié)果如下圖：

每個堿基位的平均質(zhì)量分?jǐn)?shù)

在測序數(shù)據(jù)中較高頻率出現(xiàn)的序列以及可能來源

接頭含量及可能來源

如圖所示，本次小RNA測序數(shù)據(jù)的堿基質(zhì)量不錯控硼，但由于小RNA測序建庫的特殊性泽论，且測序長度為50bp，常會測到3'端的接頭序列卡乾，因此接下來我們需要進(jìn)行接頭的去除和質(zhì)量的過濾翼悴。

2.進(jìn)行質(zhì)量過濾和接頭去除

編寫腳本：在miRNA總文件夾中建立trimmed_data路徑用于存儲過濾后數(shù)據(jù)，然后對~/miRNA/rawdata中的6個raw data.fastq使用cutadapt進(jìn)行3'端接頭的去除以及長度控制（去除3'端接頭序列后幔妨，僅保留17-35nt的reads作為潛在的miRNA序列）鹦赎，然后使用trimmomatic進(jìn)行堿基質(zhì)量的控制（reads開頭和末尾堿基質(zhì)量需>20，并以滑窗內(nèi)堿基平均質(zhì)量為20進(jìn)行滑窗質(zhì)檢误堡，在堿基修剪后僅保留>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數(shù)據(jù)進(jìn)行二次QC

編寫腳本：調(diào)用fastqc，將clean data的QC結(jié)果保存在~/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

對報告進(jìn)行簡要的解讀婚温，可見：

堿基質(zhì)量經(jīng)過過濾后是很好的

clean reads長度的分布：可見在22bp處出現(xiàn)峰。22bp正是miRNA的平均長度

現(xiàn)在檢測到的出現(xiàn)頻率較高的序列已經(jīng)不再是adapters/primers了

此前混雜在reads中的adapter序列已經(jīng)全部清除

4.過濾已經(jīng)在rfam和ensembl數(shù)據(jù)庫中注釋的rRNA媳否、tRNA栅螟、sRNA荆秦、snRNA

由于小RNA測序得到的reads中同時存在其他許多小RNA的種類，因此如果可以提前將這些可能造成混雜的序列去除力图，則可以加快后續(xù)比對的速度步绸，也減小假陽性。

rfam數(shù)據(jù)庫全面記錄了目前為止已經(jīng)注釋的小RNA序列吃媒。因此除開已經(jīng)注釋的miRNA序列瓤介，我們挑出庫中人類的rRNA、tRNA赘那、sRNA和snRNA（包含了snoRNA）對reads進(jìn)行比對刑桑，將比對上的reads去除。則剩下的reads為miRNA的可能性更大募舟。

此外為了進(jìn)一步的減小假陽性祠斧，我打算把rfam的這些ncRNA和ensembl數(shù)據(jù)庫中的ncRNA注釋合并進(jìn)行一起過濾。ensembl數(shù)據(jù)庫的大名就不多介紹了拱礁，它每一版release也包括了ncRNA序列琢锋，其類別可以分為snRNA、rRNA呢灶、scRNA吴超、sRNA、Mt_rRNA填抬、snoRNA烛芬、miRNA、Mt_tRNA飒责、scaRNA赘娄、lncRNA、ribozyme宏蛉、misc_RNA遣臼。

創(chuàng)建~/rfam路徑，在其中下載rfam數(shù)據(jù)庫中的全部fasta文件并解壓拾并；
創(chuàng)建~/ensembl路徑揍堰，下載ensembl數(shù)據(jù)庫中的Homo_sapiens.GRCh38.ncrna.fa.gz（隨基因組一起，目前最新發(fā)布版本是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文件一覽

我們通過復(fù)制rfam數(shù)據(jù)庫網(wǎng)頁中的taxonomy信息可以獲得庫中rRNA、tRNA之碗、sRNA和snRNA的對應(yīng)RF文件編號（圖中表格的第一列）蝙眶，然后通過對fasta文件中名字行（>xxx）中的物種信息可以挑出人類（Homo sapiens）的注釋。

在rfam-search中選擇相應(yīng)的taxonomy褪那，復(fù)制相應(yīng)的表格幽纷，用nano在終端中粘貼保存式塌，表格第一列可用于不同類型數(shù)據(jù)庫的提取

通過復(fù)制信息，我們得到并儲存了rfam.anno.rRNA.txt友浸、rfam.anno.tRNA.txt峰尝、rfam.anno.sRNA.txt和rfam.anno.snRNA.txt在rfam路徑中。

它們的第一列為相應(yīng)的RF號收恢。

rfam.anno.rRNA.txt的內(nèi)容

而ensembl下載的數(shù)據(jù)中僅包含人類的ncRNA fasta數(shù)據(jù)武学。其姓名列>xxx的第五列包含gene_biotype，可以提取出相應(yīng)類型的小RNA序列派诬。

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

可以使用如下命令查詢ensembl下載的ncRNA類別：

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

可得到：

ensembl收錄的ncRNA類型

編寫一個腳本劳淆，實現(xiàn)從rfam對應(yīng)特征表格的txt中獲取第一列文件信息，整合成分類的fasta文件并最后得到人類的rfam tRNA默赂、rRNA、snRNA fasta文件括勺。然后將ensembl庫中對應(yīng)類型的小RNA序列也添加到fasta數(shù)據(jù)庫中缆八。隨后，為了用于blastn比對序列疾捍，我們還需要在比對前使用makeblastdb進(jìn)行核酸比對數(shù)據(jù)庫的建立奈辰。腳本命名為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中

數(shù)據(jù)庫建立完畢乱豆，我們接下來即可利用blastdb對reads進(jìn)行比對了奖恰。

但！由于樣本中本身含有相同的小RNA序列宛裕，以及建庫時的PCR擴增瑟啃，測序獲得的reads中含有大量的相同序列。如果把相同的所有reads簡并到1條揩尸，在比對時將可以避免重復(fù)比對蛹屿，大大縮短時間。

mirdeep2軟件中的collapse_reads_md.pl腳本可以完成這一功能岩榆，輸入reads.fasta和樣本的三字符代號后错负，可將重復(fù)reads計數(shù)，以”xnumber“的形式加至序列名行（>開頭的那一行）勇边。（eg: >hsa_0_x328236犹撒，代表這個序列在hsa這個樣本中出現(xiàn)了328236次）

我們編寫如下腳本以獲得所有樣本的collapsed_reads.fasta文件。首先將儲存于~/miRNA/trimmed_data的各樣本clean data.fastq文件轉(zhuǎn)化為不含質(zhì)量信息的fasta文件粒褒，并將序列本身作為序列名字（eg: >AGCT）识颊，然后輸入collapse_reads_md.pl進(jìn)行折疊，以樣本組別前兩個字母（co/tr）+生物學(xué)重復(fù)編號（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

現(xiàn)在在trimmed_data路徑獲得了con_1/2/3_clean_collapsed.fasta和treat_1/2/3_clean_collapsed.fasta怕享。我們比較一下初始的fastq文件和collapsed_reads.fasta文件：

注意名稱列的suffix (>樣本編號_序列編號_x出現(xiàn)次數(shù))，和壓縮了不要太多的文件大小

接下來使用該腳本進(jìn)行collpased reads和之前建立的rfam人類數(shù)據(jù)庫的比對（blastn命令）镰踏，并從collapsed_reads.fasta中將這些比對上的序列去除函筋。比對的結(jié)果保存在新建立的~/miRNA/blastresult路徑中，比對結(jié)束后會出現(xiàn)con_1/2/3_clean_collapsed_rfam_filtered.fasta和treat_1/2/3_clean_collapsed_rfam_filtered.fasta奠伪，此外每個樣本含有一個clean_collapsed_blast_result文件夾含有更詳細(xì)的信息跌帐。
將該腳本命名為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循環(huán)+sed命令搜索未注釋fasta序列的速度實在太慢绊率。我的一個樣品需要處理9小時（......）谨敛。而且我發(fā)現(xiàn)直接用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路徑：根據(jù)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小時颈娜。期間就讓服務(wù)器跑吧，不用管它...（如果有更好效率更高的方式希望大家可以和我多交流Ｕ阋恕）

腳本運行結(jié)束后官辽，我們在~/miRNA/blastresult路徑中獲得了6個逃過了rfam注釋的collapsed reads.fasta。

Optional：查看過濾的數(shù)據(jù)情況

如果我們想查看本次blast比對后被過濾掉的reads數(shù)梆奈、保留的reads數(shù)以及比例野崇，可以使用如下的腳本完成。腳本命名為count_rfam_rate.sh亩钟，并保存于~/miRNA/scripts路徑乓梨。

計算方法：
對unique reads：直接查看annotated.txt和all_seq.txt的行數(shù)。
對total reads：計算annotated.txt和all_seq.txt中每行中的"xnnnn"的數(shù)字和清酥。

#!/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

運行腳本：

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

我的數(shù)據(jù)可以得到以下結(jié)果：

*用直接blastn比對的remove.rfam.annotated.sh的過濾結(jié)果 vs 用-task blastn-short比對的remove.rfam.annotated.v2.sh的過濾結(jié)果：

部分樣本的rfam過濾結(jié)果

對比結(jié)果顯然是說明blastn-short完勝扶镀。

未完待續(xù)！