【原】miRNA測序數據的上游定量流程實戰演練

miRNA的上游分析流程跟mRNA的上游流程很相似：環境部署——數據下載——查看數據(非質控)——數據質控清洗——數據比對——數據定量https://www.bilibili.com/video/BV1zK411n7qr 

1.基于conda的環境部署/軟件安裝：

嘗試使用ARM架構(M1/M2芯片) 去安裝fastqc trim-galore hisat2 subread multiqc samtools salmon fastp，發現這些軟件中有幾個是不兼容的。所以需要改回原來的x86_64架構(Intel芯片)，如果非mac/M1/M2的不需要用這種方式。

CONDA_SUBDIR=osx-64 conda create -n miRNA_x86_64 python=3.9  
conda activate miRNA_x86_64  
conda install -y -c bioconda sra-tools hisat2 bowtie samtools fastp bowtie2 fastx_toolkit fastqc  
conda install -y -c hcc aspera-cli  
conda install -y sra-tools  

2.下載相應數據庫數據

miRbase是miRNA研究領域內最權威的數據庫之一，提供了miRNAs序列以及注釋，定位，發卡序列等信息，以及提供命名服務。

mkdir mirna  
mkdir reference  
cd ./reference  
#nohup wget ftp://mirbase.org/pub/mirbase/CURRENT/hairpin.fa &  
#nohup wget ftp://mirbase.org/pub/mirbase/CURRENT/mature.fa &  
#gunzip hairpin.fa.gz  
#gunzip mature.fa.gz  
  
#不知為何，筆者這邊一直出現連接失敗，實在沒辦法就直接從官網進行了點擊下載

1.  前體 miRNA（hairpin.fa）：

• 識別新 miRNA： 通過比對發夾狀序列，研究人員可以預測或識別新的 miRNA，因為新生成的 miRNA 在細胞內首先形成發夾結構。
• 結構分析： 發夾狀結構是 miRNA 特有的二級結構，通過分析它的結構和序列特征，可以更好地了解 miRNA 的生成機制。
• 處理和轉化為成熟 miRNA： 前體 miRNA 是成熟 miRNA 的來源，前體文件可以幫助模擬或研究 miRNA 在細胞內的生成過程，例如 Dicer 酶切割的具體位置和生成機制。

2.  成熟 miRNA（mature.fa）：

• 功能性分析： 成熟 miRNA 是直接調控基因表達的功能分子，可以結合到特定的 mRNA 靶位點。mature.fa 文件可以用于 miRNA 靶向基因預測、通路分析和功能研究。
• 表達定量： 在實際的表達定量分析（如 RNA-seq）中，比對成熟 miRNA 序列可以幫助準確識別 miRNA，并進行定量，從而用于下游的差異表達分析。
• 基因調控網絡： 成熟 miRNA 文件可用于構建 miRNA-mRNA 調控網絡，研究 miRNA 在特定生物學過程中的作用。

在這里這兩個文件的作用主要是進行序列比對。

3.Check 下載到本地的數據

打開hairpin.fa文件可以看到數據的格式

1. cel-let-7 是序列名稱。
2. MI0000001 是 miRBase 數據庫中對應的唯一 ID。
3. Caenorhabditis elegans 是該序列的來源物種。
4. let-7 stem-loop 描述了該序列是 let-7 miRNA 的發卡環結構。
5. 緊接的兩行是 let-7 的核苷酸序列。

cat hairpin.fa | grep '>'| awk '{print $3,$4}'| sort |uniq -c | wc -l  
# 265

1. cat hairpin.fa: 讀取 hairpin.fa 文件的全部內容，并輸出到終端。
2. grep '>': 篩選出包含 > 字符的行。FASTA 格式中，> 開頭的行表示序列的注釋信息，如 miRNA 名稱和其他信息，而不是序列本身。
3. awk '{print 4}': 這里 {print 4} 表示輸出第三個和第四個字段（即以空格或制表符分隔的第三和第四部分）。在 miRNA FASTA 文件中，第三個和第四個字段可能是與 miRNA 名稱和種類相關的信息。
4. sort對提取的第三和第四字段進行排序。
5. uniq -c: 統計每個唯一的第三、四字段組合的出現次數。uniq -c 會對相同的行進行計數。例如，如果 miRNA_type1 出現了多次，則會輸出類似 5 miRNA_type1。在使用 uniq 之前，必須先對內容進行 sort，否則無法識別相同的行。
6. wc -l：統計輸出的行數。wc -l 統計 uniq -c 輸出的總行數，即不同 miRNA 類型組合的數量。

cat mature.fa | grep '>'| awk '{print $3,$4}'| sort |uniq -c | wc -l  
# 265

接著觀察人類這個物種的miRNA的數量

grep sapiens mature.fa |wc -l  
# 2656  

1. grep sapiens mature.fa | wc -l：grep sapiens mature.fa：從文件 mature.fa 中查找包含 "sapiens" 的行。| wc -l：將匹配的行數通過管道傳遞給 wc -l，統計這些行的數量。最終返回包含 "sapiens" 的總行數。
2. grep sapiens hairpin.fa | wc：grep sapiens hairpin.fa：從文件 hairpin.fa 中查找包含 "sapiens" 的行。| wc：wc 會返回三個值：行數、單詞數、字節數。由于沒有加 -l 參數，結果中會包含所有這三個統計值，列出包含 "sapiens" 的行數、單詞總數以及字符總數。

接著觀察有多少序列，4行為一條序列

zless -S mature.fa | paste - - - - |wc -l  
# 24443  
zless -S hairpin.fa | paste - - - - |wc -l  
# 30029

接著檢查一下前體和成熟體長度：

前體miRNA和成熟體miRNA：前體miRNA長度一般是70-120堿基，通常是莖環(發卡，hairpin)結構。成熟之后一般是22個堿基。(曾老師的perl的示例代碼)

# 前體長度  
awk '/^>/ {printf("\n%s\t",$0);next;} {printf("%s",$0);} END {printf("\n");}' < hairpin.fa | tr "\t" "\n" | grep -v '>' | awk '{print length}' | uniq -c | sort -n -k2  
  
# 成熟體長度  
awk '/^>/ {printf("\n%s\t",$0);next;} {printf("%s",$0);} END {printf("\n");}' < mature.fa | tr "\t" "\n" | grep -v '>' | awk '{print length}' | uniq -c | sort -n -k2

4.構建索引

構建 miRNA 序列的索引庫（例如使用 bowtie 構建 hairpin.fa 和 mature.fa 的索引）可以顯著提升后續比對過程的速度和準確性，比如：1. 加速比對過程；2. 減少計算資源需求；3. 提升比對準確性；

U->T轉換

為什么要進行U-> T轉換：在 RNA 序列中，堿基用“U”（尿嘧啶）表示，而 DNA 序列中對應的是“T”（胸腺嘧啶）。大多數比對工具，如 Bowtie，主要是針對 DNA 序列設計的，因此它們默認識別“ATCG”四種堿基。在這種情況下，如果不將 RNA 中的“U”轉換為“T”，比對工具會無法正確識別和比對 RNA 序列。

perl -alne '{if(/^>/){if(/Homo/){$tmp=1}else{$tmp=0}};next if $tmp!=1;s/U/T/g if !/>/;print}' hairpin.fa > hairpin.human.fa  
perl -alne '{if(/^>/){if(/Homo/){$tmp=1}else{$tmp=0}};next if $tmp!=1;s/U/T/g if !/>/;print}' mature.fa > mature.human.fa

1. perl -alne '...'： perl：調用 Perl 腳本。-a：啟用自動分段模式，將每行分割成字段，保存在 @F 數組中（這里未用到 @F）。-l：自動處理換行符，使輸出更整齊。-n：循環讀取每一行，但不會自動打印輸出。
2. if(/^>/){if(/Homo/){$tmp=1}else{$tmp=0}};：/^>/ ：檢測行是否以 > 開頭，這通常表示 FASTA 格式中的序列 ID 行。**if(/Homo/)**：如果 ID 行中包含“Homo”（指代人類相關的序列），則將 $tmp 設置為 1（即允許輸出該序列的標記）；否則設置為 0（排除該序列）。這一步確定每條序列是否為人類序列，僅處理包含 Homo 的序列。
3. next if $tmp!=1;： next if tmp 不是 1，跳過該行，即非人類序列直接跳過不處理。
4. s/U/T/g if !/>/;：s/U/T/g：將行中的“U”替換為“T”，全局替換（即一行中所有“U”均替換為“T”）。if !/>/：僅當行中不含 > 符號時執行替換，即應用在實際序列行，而非序列 ID 行。
5. print：print：打印處理后的行。對符合條件的序列和序列 ID 均輸出到指定文件。
6. > hairpin.human.fa 和 > mature.human.fa：>將標準輸出重定向到文件，分別保存到 hairpin.human.fa 和 mature.human.fa。

Bowtie和Bowtie2的區別是什么：

1. Bowtie：采用基于 Burrows-Wheeler 變換（BWT）和 FM-index 的算法，適合對短序列（通常為小于 50bp 的短 RNA 或短讀長 DNA）進行快速比對。
2. Bowtie2：采用了更復雜的比對算法，使用 Burrows-Wheeler 變換和動態規劃算法來支持長片段的局部和全局比對，因此適合較長的讀長（一般在 50bp 以上），包括 DNA 和 RNA-seq 數據。

bowtie-build hairpin.human.fa hsa-hairpin-bowtie-index  
bowtie-build mature.human.fa hsa-mature-bowtie-index  
  
# check  
ls -lh

5.下載數據

勾選想要下載的數據，并點擊accession list，并把list放入mirna文件夾中

cd ../  
cd ./mirna  
  
# check  
ls -lh

使用prefetech下載數據，這里需要把SRRlist和SRA toolkit軟件安裝好。除了這種方式，我們也可以選擇aspera下載方式

nohup bash -c 'cat SRR_Acc_List.txt | while read id; do  
echo "Downloading $id"  
prefetch $id  
done' &> download.log &  

把sra數據批量轉換為fastq數據

# 首先需要知道fastq-dump工具的位置  
which fastq-dump  
# /opt/anaconda3/envs/miRNA_x86_64/bin/fastq-dump  
  
# 查看文件夾中的數據是怎么樣的  
ls  
  
# 要明確一下echo和SRR的ID情況  
# 輸入進終端的時候一定要再三明確代碼情況  
dump=/opt/anaconda3/envs/miRNA_x86_64/bin/fastq-dump  
echo {02..25} |sed 's/ /\n/g' |while read id; \  
do ( $dump -O ./ --gzip --split-3 SRR154179${id}) ;\  
done  
  
# 數據有點大 筆者就分開下載了  
dump=/opt/anaconda3/envs/miRNA_x86_64/bin/fastq-dump  
echo {35..50} | sed 's/ /\n/g' | while read id; do  
($dump -O ./ --gzip --split-3 SRR154179${id})  
done

1. {02..25} 會生成一個從 02 到 25 的數字序列。
2. sed 's/ /\n/g' 將生成的序列號中的空格替換為換行符，以便逐行讀取數字。
3. while read id; do ... done 形成一個循環，逐行讀取序列號并存儲在變量 id 中。
4. 其中 {id}：指定待轉換的 SRA 文件，${id} 會插入循環讀取的數字部分，生成類似 SRR15417902、SRR15417903 等文件名。

6.數據質控和清洗

數據質控查看

# 對當前文件夾中所有以fastq.gz文件進行質量控制  
fastqc -t 2 -o ./ *.fastq.gz  
# 對當前文件夾中所有以fastq.gz文件進行整合質量控制，輸出到fastq_qc文件夾中  
multiqc ./*zip -o ./fastq_qc

1. fastqc：調用 FastQC 工具，用于分析測序數據的質量。
2. -t 2：指定使用 2 個線程并行處理文件，以加快分析速度。
3. -o ./：指定輸出目錄為當前目錄（./），FastQC 生成的報告文件將保存在當前目錄中。
4. *.fastq.gz：匹配當前目錄下所有以 .fastq.gz 結尾的文件，作為輸入文件進行質量控制分析。

正式數據清洗

# 檢查文件壓縮格式  
file /Users/zaneflying/Desktop/miRNA_Z/mirna/SRR15417902.fastq.gz  
  
# trim+clean  
# 文章用了miRquant 2.0這個工具進行trim，但筆者還是按照曾老師的流程進行處理  
ls /Users/zaneflying/Desktop/miRNA_Z/mirna/*.gz | while read id; do  
echo $id  
gzip -cd $id | fastq_quality_filter -v -q 20 -p 80 -Q33 -i - -o tmp  
fastx_trimmer -v -f 1 -l 27 -m 15 -i tmp -Q33 -z -o ${id%%.*}_clean.fq.gz  
fastqc -t 2 -o ./ ${id%%.*}_clean.fq.gz  
done  

# check一下  
ls -lh *_clean.fq.gz

7.數據比對

根據miRBase數據庫下載的序列進行比對和定量。

# mature清洗和定量  
index=/Users/zaneflying/Desktop/miRNA_Z/reference/hsa-mature-bowtie-index  
  
ls *_clean.fq.gz | while read id; do  
echo $id  
bowtie -p 2 -x $index $id -S tmp  
samtools view -bS -@ 2 tmp -o ${id}_mature.bam  
done  
  
# hairpin清洗和定量  
index=/Users/zaneflying/Desktop/miRNA_Z/reference/hsa-hairpin-bowtie-index.  
  
ls *_clean.fq.gz | while read id; do  
echo $id  
bowtie -p 2 -x $index $id -S tmp  
samtools view -bS -@ 2 tmp -o ${id}_hairpin.bam  
done

1. *ls _clean.fq.gz: 列出所有以 _clean.fq.gz 結尾的文件，即預處理過的 miRNA 序列文件。
2. while read id; do ... done: 使用 while 循環逐個讀取文件名并將其賦值給變量 id，然后對每個文件執行循環內的命令。
3. echo $id: 打印當前正在處理的文件名，用于檢查進度。
4. bowtie -p 2 -x id -S tmp: -p 2：指定使用 2 個線程來加速處理。-x index 應該是您之前創建的 miRNA 參考序列的索引（在您的例子中應該是 /Users/zaneflying/Desktop/miRNA_Z/reference/mature）。$id：輸入文件，即當前的 _clean.fq.gz 文件。-S tmp**：指定輸出文件 tmp，這將是一個 SAM 格式的中間文件。
5. samtools view -bS -@ 2 tmp -o {id}_mature.bam：輸出文件，將 BAM 文件命名為輸入文件名加 _mature.bam 后綴。

對比結果中發現只有1507條reads對應上，也就是說幾乎所有都沒有比對上的情況，很費解。應該是我沒有學好：

然后嘗試更換一下參考基因組，文章中提到的是hg19

筆者這里使用GRCh38進行對比，不過這個并不是重點哈。下載流GRCh38程可見轉錄組上游分析流程推文。

# 生成索引  
bowtie2-build Homo_sapiens.GRCh38.dna.primary_assembly.fa GRCh38.dna  

# index索引條目  
index=/Users/zaneflying/Desktop/miRNA_Z/GRCh38.113/GRCh38.dna  
# bowtie開始轉化  
ls *_clean.fq.gz | while read id  
do  
echo $id  
bowtie -p 2 -x $index $id -S tmp ;  
samtools view -bS -@ 2 tmp -o ${id}_genome.bam  
done

1. ls *_clean.fq.gz | while read id：列出當前目錄下所有文件名以 _clean.fq.gz 結尾的文件。while read id 用來逐行讀取這些文件名，并將文件名存儲在變量 id 中。
2. echo $id：打印當前正在處理的文件名，以便追蹤進度。
3. bowtie -p 2 -x id -S tmp：使用 bowtie 對文件 index 指定索引文件路徑，即 id 是輸入的 FASTQ 文件。-S tmp 將比對結果輸出為 SAM 格式，臨時保存在文件 tmp 中。
4. samtools view -bS -@ 2 tmp -o {id}_genome.bam 指定輸出文件名為 ${id}_genome.bam，即與輸入文件同名，但以 _genome.bam 結尾。

這個對比結果情況就勉強能“讓人接受”啦~

8.數據定量

文章中用的是miRquant 2.0

筆者使用featurecounts去定量， 需要先去miRBase上下載hsa.gff3

# 下載hsa.gff文件,放到reference文件夾中  
nohup wget ftp://mirbase.org/pub/mirbase/CURRENT/genomes/hsa.gff3 &  
# 如果網絡不好，就直接手動下載  
  
# 安裝subread  
conda install -c bioconda subread  
  
# 開始比對  
gtf=/Users/zaneflying/Desktop/miRNA_Z/reference/hsa.gff3  
featureCounts -T 2 -F gff -M -t miRNA -g Name -a $gtf -o all.counts.mature.txt *genome* 1>counts.mature.log 2>&1  
featureCounts -T 2 -F gff -M -t miRNA_primary_transcript -g Name -a $gtf -o all.counts.hairpin.txt *genome* 1>counts.hairpin.log 2>&1  

1. featureCounts：調用 featureCounts 工具進行基因計數。-T 2：指定使用 2 個線程。-F gff：指定輸入文件的注釋格式為 GFF。-M：允許多重比對的 reads（即與多個位置比對的 reads）。-t miRNA：在 GFF 文件中，選擇類型為 miRNA 的條目，這樣可以僅對成熟 miRNA 計數。-g Name：選擇 GFF 文件中 Name 字段作為基因 ID。-a gtf 變量。-o all.counts.mature.txt：輸出文件名稱，包含成熟 miRNA 的計數。genome：指定輸入的 BAM 文件，genome 是通配符，表示所有名稱包含“genome”的 BAM 文件。1>counts.mature.log 2>&1：將標準輸出和標準錯誤重定向到日志文件 counts.mature.log。
2. 第二行命令與第一行基本相同，唯一的區別是 -t 選項為 miRNA_primary_transcript，用于選擇前體 miRNA 條目，從而統計前體 miRNA 的計數。輸出文件為 all.counts.hairpin.txt，并將日志信息輸出到 counts.hairpin.log。

因為比對有問題，定量也很難保證，所以拿到了矩陣也很難進行后續分析：

后續的分析基本上等同于轉錄組測序表達量矩陣，就是差異分析等統計可視化：

參考資料：

1. Multiomic analysis of microRNA-mediated regulation reveals a proliferative axis involving miR-10b in fibrolamellar carcinoma. JCI Insight. 2022 Jun 8;7(11):e154743.
2. DNAJB1-PRKACA fusion protein-regulated LINC00473 promotes tumor growth and alters mitochondrial fitness in fibrolamellar carcinoma. PLoS Genet 2024 Mar;20(3):e1011216.
3. Chemical, Molecular, and Single-nucleus Analysis Reveal Chondroitin Sulfate Proteoglycan Aberrancy in Fibrolamellar Carcinoma. Cancer Res Commun 2022 Jul;2(7):663-678.
4. 生信技能樹：

致謝：感謝曾老師以及生信技能樹團隊全體成員。

注：若對內容有疑惑或者有發現明確錯誤的朋友，請聯系后臺(歡迎交流)。