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. 生信技能树：

致谢：感谢曾老师以及生信技能树团队全体成员。

注：若对内容有疑惑或者有发现明确错误的朋友，请联系后台(欢迎交流)。