miRNA-seq小RNA高通量测序pipeline：从raw

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

到此时，原始的单端小RNA测序已经完成了质量过滤和接头修剪、rfam+ensembl非编码小RNA的序列过滤工作。

5.与参考基因组外显子和内含子的比对

由于提取得到的短片段RNA中可能含有一些mRNA的降解片段，而目前已知miRNA可能存在于内含子中，所以接下来对通过非编码RNA过滤的reads进行外显子和内含子的比对。

策略为：reads比对exon genome→不能比对上exon的reads直接通过→能比对上exon的reads进行intron genome的比对→不能比对上intron的reads丢弃/能比对上则通过。

比对参数均为：最多允许20个比对位置，最大允许1个错配，开2个线程。-al：比对上的reads输出；-un：未必对上的reads输出。

首先在UCSC genome browser的table browser上下载hg38的全基因组外显子和内含子fa文件到~/hg.38.ref（如下图table browser选择各选项 > get output > genome > 外显子为CDS+UTR/内含子为intron），然后使用bowtie build index。

UCSC genome table browser

cd ~/miRNA/hg.38.ref/hg38.index
bowtie-build ../hg38.exon.fa hg.38.exon
bowtie-build ../hg38.intron.fa hg.38.intron

完成index后，编写脚本完成过滤步骤：脚本命名为exon_intron_filtering.sh，储存于~/miRNA/scripts路径。

#!/bin/bash
if [ -d ~/miRNA/exon.intron.filter ]
then
    echo 'directory for exon and intron filtering has existed.'
else
    mkdir ~/miRNA/exon.intron.filter
cd ~/miRNA/blastresult
for each in `ls *.fasta`
do
    cd ~/miRNA/exon.intron.filter
    bowtie -S -f -a -v 1 --best --strata -m 20 -p 2 -al ${each}_for_intron.fa -un ${each}_no_exon.fa \
    ~/hg.38.ref/hg38.index/hg.38.exon $each ${each}_exon_alignment.sam
    rm ${each}_exon_alignment.sam
    bowtie -S -f -a -v 1 --best --strata -m 20 -p 2 -al ${each}_intron_positive.fa \
    ~/hg.38.ref/hg38.index/hg.38.intron ${each}_for_intron.fa ${each}_intron_alignment.sam
    rm ${each}_intron_alignment.sam
    rm ${each}_for_intron.fa
    cat ${each}_no_exon.fa ${each}_intron_positive.fa > ${each}_exintr_filtered.fa
    rm ${each}_no_exon.fa
    rm ${each}_intron_positive.fa
done 

然后运行脚本：

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

可以获得6个完成所有过滤流程的fa文件。

6.mirdeep2软件mapper module进行基因组比对

使用mirdeep2中的mapper.pl可以对已经完成其他RNA类型过滤的样本进行基于Bowtie的参考基因组比对。由于预测novel miRNA的核心模组miRDeep2.pl需要mapper.pl的输出，因此首先进行这一步操作。

考虑到后续miRDeep2.pl对参考基因组的要求：其identifier列不允许出现空格（我领教过使用默认参考基因组比对好之后报错的奇景，奇葩），我们需要对下载的参考基因组进行去空格处理。然后基于处理后的fa文件进行index。

cd ~/hg.38.ref
sed -i 's/[ \t]*//g' Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa
cd ./hg38.index
bowtie-build ../Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa hg.38.whitespace.removed

接下来可以使用mirdeep2进行比对了。

输入：由于miRDeep2.pl的输入只能是一个fasta文件/config.txt（其中关系比较复杂，本pipeline不采用），所以我们合并各样本此前通过collapse_reads.pl得到并完成过滤步骤的unique_reads.fasta，从而将所有样本的reads同时输入进行预测，正好也不用调用-m参数（也就是封装了collapse_reads.pl而已）。此前我们已经通过collapse_reads在identifier列记录下了各unique reads的样本代号和reads出现次数，这种identifier的格式是适配于mirdeep2进行后续鉴定和定量分析的。（下面miRDeep2.pl会再讲到）

格式说明

输出：在该参数条件下，该module会对mirdeep2记录比对信息的arf格式文件。

这里用到的mapper.pl参数：-c：输入文件为fasta格式；-l：丢弃长度低于17nt的reads；-r：能容许的最高reads比对到基因组处数；-v：显式进程；-n：覆盖已存在的文件；-u：不去除临时文件；-q：容许一个mismatch；-p：参考基因组；-t：输出比对信息.arf。

cd ~/miRNA/exon.intron.filter
cat *.fa > all.fa
mkdir ~/miRNA/mirdeep_output cd ~/miRNA/mirdeep_output
mapper.pl ~/miRNA/exon.intron.filter/all.fa -v -u -c -l 17 -r 10 -n -q -p ~/hg.38.ref/hg38.index/hg.38.whitespace.removed -t all_aligned.arf

7.mirdeep2软件mirdeep2 module进行已知已知/novel miRNA的鉴定和预测

下载并提取miRbase数据库中人类的miRNA前体以及成熟miRNA的fasta文件。由于库中fasta文件为RNA序列，因此需要将U转换成T后再和reads进行比对。此外，还是因为miRDeep2.pl对fa文件的indentifier不允许有空格的原因，我们也对这两个文件进行去空格，导致identifier列所有信息挤在了一起（如：hsa-let-7c-5pMIMAT0000064Homosapienslet-7c-5p）。（做的精细一点可以只留miRNA的名字。如：hsa-miR-200a-3p。我这里偷了点懒，算啦）

mkdir ~/mirbase && cd ~/mirbase
wget ftp://mirbase.org/pub/mirbase/CURRENT/mature.fa.gz -P ~/mirbase
wget ftp://mirbase.org/pub/mirbase/CURRENT/hairpin.fa.gz -P ~/mirbase
gunzip *.gz
sed -n '/Homo sapiens/{p;n;p}' mature.fa | sed 's/U/T/g' > hsa_mature.fa 
sed -n '/Homo sapiens/{p;n;p}' hairpin.fa | sed 's/U/T/g' > hsa_hairpin.fa
sed -i 's/[ \t]*//g' hsa_mature.fa
sed -i 's/[ \t]*//g' hsa_hairpin.fa

接下来对所有样本中比对上的reads进行miRNA的识别或新miRNA的鉴定。

注：miRDeep2.pl对文件的位置有要求：miRDeep2 reads.fa genome.fa alignment_from_mapper.arf mature_miRNA_of_this_species.fa mature_miRNA_of_related_species.fa hairpin_miRNA_of_this_species.fa -t species
后三个文件若有空缺，需用none替代。

参数：-t：物种；-P：如果mirbase使用的miRNA id包含了-3p/5p时使用；-b：根据mirdeep score阈值对novel miRNA的预测结果进行过滤。这里以4为阈值是因为>4时信噪比>10，可以得到较高质量的novel miRNA预测结果。

cd ~/miRNA/mirdeep_output
miRDeep2.pl ~/miRNA/exon.intron.filter/all.fa ~/hg.38.ref/Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa all_aligned.arf ~/mirbase/hsa_mature.fa none ~/mirbase/hsa_hairpin.fa -t Human -b 4 -P 

mirdeep2核心组件开始运行了

经过一个多小时的运行，我们获得了一系列的结果：
（1）记录了识别的miRNA和预测的novel miRNA信息的csv和html：最重要的是包含了新预测miRNA的成熟片段、star片段（尽管它已经过时了）、hairpin前体序列（图上）
（2）对已知miRNA在样品间定量的表达量矩阵csv（图中）
（3）对所有已知和预测miRNA的样本来源和二级结构进行解读的pdf（图下）

识别到的miRNA信息 已知miRNA的定量矩阵 预测miRNA的二级结构及reads来源pdf

8. mirdeep2软件quantifier module进行novel miRNA的表达量定量

我们获得了novel miRNA的三套序列，但尚未对其定量。所以接下来使用mirdeep2软件的quantifier module。我们从获得的一系列已知和预测的novel miRNA汇总csv表格中获取包含新miRNA序列的hairpin、mature seq和star seq.fasta文件，作为新miRNA（mature和star序列分开定量）定量的参考。

cd ~/miRNA/mirdeep_output
cat result_26_05_2020_t_03_12_40.csv | sed -n '28,36p' > novel_miRNA_seq.txt
cat novel_miRNA_seq.txt | awk '{print ">"$1"\n"$16}' | sed 's/u/t/g' > novel_mature.fa
cat novel_miRNA_seq.txt | awk '{print ">""$1_star""\n"$17}' | sed 's/u/t/g' > novel_star.fa
cat novel_miRNA_seq.txt | awk '{print ">"$1"\n"$18}' | sed 's/u/t/g' > novel_hairpin.fa
quantifier.pl -p novel_hairpin.fa -m novel_mature.fa -s novel_star.fa -r ~/miRNA/exon.intron.filter/all.fa

最后将在运行miRDeep2.pl时与已经获得的已知miRNA的表达情况进行合并。

cat miRNAs_expressed_all_samples_26_05_2020_t_03_12_40.csv | awk '{print $1"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10}' > ../expression.matrix.txt
cat miRNAs_expressed_all_samples_1590465484.csv | sed -n '2,$p' | awk '{print $1"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10}' >> ../expression.matrix.txt

所以我们最终获得了这次小RNA测序的已知+novel miRNA的全部表达矩阵：

表达矩阵

Optional:

我仔细观察了一下跑程序的进程信息，发现miRDeep2.pl对已知miRNA的定量也是依赖于quantifier.pl来计算。而且由于设置了参数-P（适配于加了臂名的新版mirbase的miRNA名称，如miR-88-3p/5p），有一些缺少臂名的成熟miRNA序列（如miR-888）没有办法和相应的precursor配对，所以miRNA的定量可能潜在地存在问题（我不确定）。

可能干扰miRNA定量的部分

所以这次我的策略是：
使用之前mapper.pl获得的arf，只用miRDeep2.pl得到的novel miRNA序列，丢弃它的known miRNA定量结果，将mature和star序列都并入到hsa_mature.fa中，将novel hairpin并入到hsa_hairpin.fa，然后不加参数-P，使用quantifier.pl对known和novel miRNA进行统一的定量。

有点强迫症，既然要对fa文件的identifier去空格，还是把从mirbase下载的mature和hairpin数据库的identifier格式精简一下（之前直接去空格，可以看到都挤在一起了）。毕竟输出miRNA的名称就是identifier，现在规范好可以免了后期的麻烦。

cd ~/mirbase
awk -F " " '{if($1~/^>.*/){print $1"-"$6}else{print $0}}' hairpin.fa | sed -n '/hsa/{p;n;p}' | sed 's/[uU]/T/g' > hsa_hairpin.fa
awk -F " " '{if($1~/^>.*/){print $1}else{print $0}}' mature.fa | sed -n '/hsa/{p;n;p}'  | sed 's/[uU]/T/g' > hsa_mature.fa

好的，现在mirbase下载的miRNA序列库中的identifier就清清爽爽的了，最后输出的miRNA名字也不需要再修整：

修整后的数据库identifier

接下来就是从之前miRDeep2.pl的输出中找出，并把预测的novel miRNA序列并入到mirbase下载-修整过的种子库中：（我把novel miRNA的identifier也做了一下修饰：前体为hsa-xxx-stem-loop；mature序列为hsa-xxx-mature；star序列为hsa-xxx-star。前面的hsa用于适配软件的-t参数，我们既然输入了Human，那么库中只有带hsa的序列才会被挑出定量）

cd ~/mirbase
cat ~/miRNA/mirdeep_output/novel_miRNA_seq.txt | awk '{print ">hsa-"$1"-stem-loop""\n"$18}' | sed 's/[uU]/t/g' > novel_hairpin.fa
cat ~/miRNA/mirdeep_output/novel_miRNA_seq.txt | awk '{print ">hsa-"$1"-star""\n"$17}' | sed 's/[uU]/t/g' > novel_star.fa
cat ~/miRNA/mirdeep_output/novel_miRNA_seq.txt | sed 's/[uU]/t/g' | awk '{print ">hsa-"$1"-mature""\n"$16}'  > novel_mature.fa
cat hsa_mature.fa novel_mature.fa novel_star.fa > known_novel_mature.fa
cat hsa_hairpin.fa novel_hairpin.fa > known_novel_hairpin.fa

然后就是用mirdeep2的quantifier.pl对所有已知和novel miRNA进行同时定量。不使用-P参数（识别miRNA时不会揪着-3/5p不放），不输入star序列，只输入mature序列。

cd ~/miRNA/mirdeep_output
quantifier.pl -p ~/mirbase/known_novel_hairpin.fa -m ~/mirbase/known_novel_mature.fa -r ../exon.intron.filter/all.fa -t Human -y now -k

获得结果如下：

获得的表达矩阵

我大致地和前面得到的结果做了一个比较，绝对read counts是基本一致的，但是normalized counts只能用后面这种方法获得。（前面的方法是拼接了两次定量的表达矩阵，所以标准化是不准的，已经通过代码删掉没有输出了。但是也可以自己算一下RPM即reads per million，也很简单）

我推荐用这个option，获得更为稳妥的表达矩阵。

上游分析到此结束。

最后：由于存在不同前体对应同一mature序列的情况，所以表达矩阵中的miRNA可能存在重复，但如下图所示，reads数差异并不大，所以随机去重应该没有问题。接下来就可以使用基于R的edgeR或DESeq2包进行差异表达分析了。若要对数据进行建模，可使用RPM（reads per million）方法进行标准化后使用（因为miRNA长度差异不大，所以可以忽略基因长度差异对测序reads数的影响：RPKM/FPKM/TPM考虑了基因长度）。

表达矩阵中重复miRNA的展示

全文完！