該教程分為2部分森枪，第1部分在：miRNA-seq小RNA高通量測序pipeline：從raw reads，鑒定已知miRNA-預測新miRNA审孽，到表達矩陣【一】
到此時县袱，原始的單端小RNA測序已經(jīng)完成了質量過濾和接頭修剪、rfam+ensembl非編碼小RNA的序列過濾工作佑力。
5.與參考基因組外顯子和內(nèi)含子的比對
由于提取得到的短片段RNA中可能含有一些mRNA的降解片段式散，而目前已知miRNA可能存在于內(nèi)含子中，所以接下來對通過非編碼RNA過濾的reads進行外顯子和內(nèi)含子的比對打颤。
策略為：reads比對exon genome→不能比對上exon的reads直接通過→能比對上exon的reads進行intron genome的比對→不能比對上intron的reads丟棄/能比對上則通過暴拄。
比對參數(shù)均為：最多允許20個比對位置，最大允許1個錯配瘸洛，開2個線程揍移。-al：比對上的reads輸出；-un：未必對上的reads輸出反肋。
首先在UCSC genome browser的table browser上下載hg38的全基因組外顯子和內(nèi)含子fa文件到~/hg.38.ref（如下圖table browser選擇各選項 > get output > genome > 外顯子為CDS+UTR/內(nèi)含子為intron），然后使用bowtie build index踏施。
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可以對已經(jīng)完成其他RNA類型過濾的樣本進行基于Bowtie的參考基因組比對。由于預測novel miRNA的核心模組miRDeep2.pl需要mapper.pl的輸出日熬，因此首先進行這一步操作棍厌。
考慮到后續(xù)miRDeep2.pl對參考基因組的要求：其identifier列不允許出現(xiàn)空格（我領教過使用默認參考基因組比對好之后報錯的奇景，奇葩）竖席，我們需要對下載的參考基因組進行去空格處理耘纱。然后基于處理后的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參數(shù)（也就是封裝了collapse_reads.pl而已）弄慰。此前我們已經(jīng)通過collapse_reads在identifier列記錄下了各unique reads的樣本代號和reads出現(xiàn)次數(shù)，這種identifier的格式是適配于mirdeep2進行后續(xù)鑒定和定量分析的蝶锋。（下面miRDeep2.pl會再講到）
輸出：在該參數(shù)條件下陆爽，該module會對mirdeep2記錄比對信息的arf格式文件。
這里用到的mapper.pl參數(shù)：-c：輸入文件為fasta格式扳缕；-l：丟棄長度低于17nt的reads慌闭；-r：能容許的最高reads比對到基因組處數(shù)；-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數(shù)據(jù)庫中人類的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替代卖怜。
參數(shù)：-t：物種；-P：如果mirbase使用的miRNA id包含了-3p/5p時使用阐枣；-b：根據(jù)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
經(jīng)過一個多小時的運行甩鳄，我們獲得了一系列的結果：
（1）記錄了識別的miRNA和預測的novel miRNA信息的csv和html：最重要的是包含了新預測miRNA的成熟片段、star片段（盡管它已經(jīng)過時了）宪哩、hairpin前體序列（圖上）
（2）對已知miRNA在樣品間定量的表達量矩陣csv（圖中）
（3）對所有已知和預測miRNA的樣本來源和二級結構進行解讀的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時與已經(jīng)獲得的已知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:
我仔細觀察了一下跑程序的進程信息甜熔，發(fā)現(xiàn)miRDeep2.pl對已知miRNA的定量也是依賴于quantifier.pl來計算。而且由于設置了參數(shù)-P（適配于加了臂名的新版mirbase的miRNA名稱突倍，如miR-88-3p/5p）腔稀，有一些缺少臂名的成熟miRNA序列（如miR-888）沒有辦法和相應的precursor配對，所以miRNA的定量可能潛在地存在問題（我不確定）羽历。
所以這次我的策略是：
使用之前mapper.pl獲得的arf焊虏，只用miRDeep2.pl得到的novel miRNA序列，丟棄它的known miRNA定量結果秕磷，將mature和star序列都并入到hsa_mature.fa中诵闭，將novel hairpin并入到hsa_hairpin.fa，然后不加參數(shù)-P澎嚣，使用quantifier.pl對known和novel miRNA進行統(tǒng)一的定量疏尿。
有點強迫癥，既然要對fa文件的identifier去空格易桃，還是把從mirbase下載的mature和hairpin數(shù)據(jù)庫的identifier格式精簡一下（之前直接去空格褥琐，可以看到都擠在一起了）。畢竟輸出miRNA的名稱就是identifier晤郑，現(xiàn)在規(guī)范好可以免了后期的麻煩敌呈。
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
好的，現(xiàn)在mirbase下載的miRNA序列庫中的identifier就清清爽爽的了造寝，最后輸出的miRNA名字也不需要再修整：
接下來就是從之前miRDeep2.pl的輸出中找出驱富，并把預測的novel miRNA序列并入到mirbase下載-修整過的種子庫中：（我把novel miRNA的identifier也做了一下修飾：前體為hsa-xxx-stem-loop；mature序列為hsa-xxx-mature匹舞；star序列為hsa-xxx-star。前面的hsa用于適配軟件的-t參數(shù)线脚，我們既然輸入了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參數(shù)（識別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只能用后面這種方法獲得。（前面的方法是拼接了兩次定量的表達矩陣伶选，所以標準化是不準的史飞，已經(jīng)通過代碼刪掉沒有輸出了尖昏。但是也可以自己算一下RPM即reads per million，也很簡單）
我推薦用這個option构资，獲得更為穩(wěn)妥的表達矩陣抽诉。
上游分析到此結束。
最后：由于存在不同前體對應同一mature序列的情況吐绵，所以表達矩陣中的miRNA可能存在重復迹淌，但如下圖所示，reads數(shù)差異并不大己单，所以隨機去重應該沒有問題唉窃。接下來就可以使用基于R的edgeR或DESeq2包進行差異表達分析了。若要對數(shù)據(jù)進行建模纹笼，可使用RPM（reads per million）方法進行標準化后使用（因為miRNA長度差異不大纹份，所以可以忽略基因長度差異對測序reads數(shù)的影響：RPKM/FPKM/TPM考慮了基因長度）。
全文完允乐！
