最近miRNA领域居然获得了诺贝尔奖，所以我也凑一下热闹想学习，所以就拜托神通广大的曾老师提供了一个miRNA测序数据案例，这个数据集是GSE181922，其中一共包括了40例样品的的miRNA数据。如下所示：
{jz:field.toptypename/}图片
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-tools2.下载相应数据库数据
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文件可以看到数据的格式
cel-let-7 是序列名称。MI0000001 是 miRBase 数据库中对应的唯一 ID。Caenorhabditis elegans 是该序列的来源物种。let-7 stem-loop 描述了该序列是 let-7 miRNA 的发卡环结构。紧接的两行是 let-7 的核苷酸序列。图片
cat hairpin.fa | grep '>'| awk '{print $3,$4}'| sort |uniq -c | wc -l # 265cat hairpin.fa: 读取 hairpin.fa 文件的全部内容，并输出到终端。grep '>': 筛选出包含 > 字符的行。FASTA 格式中，> 开头的行表示序列的注释信息，如 miRNA 名称和其他信息，而不是序列本身。awk '{print 图片
4} 表示输出第三个和第四个字段（即以空格或制表符分隔的第三和第四部分）。在 miRNA FASTA 文件中，第三个和第四个字段可能是与 miRNA 名称和种类相关的信息。sort对提取的第三和第四字段进行排序。uniq -c: 统计每个唯一的第三、四字段组合的出现次数。uniq -c 会对相同的行进行计数。例如，如果 miRNA_type1 出现了多次，则会输出类似 5 miRNA_type1。在使用 uniq 之前，必须先对内容进行 sort，否则无法识别相同的行。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 # 2656grep sapiens mature.fa | wc -l：grep sapiens mature.fa：从文件 mature.fa 中查找包含 "sapiens" 的行。| wc -l：将匹配的行数通过管道传递给 wc -l，统计这些行的数量。最终返回包含 "sapiens" 的总行数。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
接着检查一下前体和成熟体长度：
建站客服QQ：88888888前体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 -k24.构建索引构建 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.faperl -alne '...'： perl：调用 Perl 脚本。-a：启用自动分段模式，将每行分割成字段，保存在 @F 数组中（这里未用到 @F）。-l：自动处理换行符，使输出更整齐。-n：循环读取每一行，但不会自动打印输出。if(/^>/){if(/Homo/){$tmp=1}else{$tmp=0}};：/^>/ ：检测行是否以 > 开头，这通常表示 FASTA 格式中的序列 ID 行。**if(/Homo/)**：如果 ID 行中包含“Homo”（指代人类相关的序列），则将 $tmp 设置为 1（即允许输出该序列的标记）；否则设置为 0（排除该序列）。这一步确定每条序列是否为人类序列，仅处理包含 Homo 的序列。next if $tmp!=1;： next if 图片
Bowtie和Bowtie2的区别是什么：
Bowtie：采用基于 Burrows-Wheeler 变换（BWT）和 FM-index 的算法，适合对短序列（通常为小于 50bp 的短 RNA 或短读长 DNA）进行快速比对。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{02..25} 会生成一个从 02 到 25 的数字序列。sed 's/ /\n/g' 将生成的序列号中的空格替换为换行符，以便逐行读取数字。while read id; do ... done 形成一个循环，逐行读取序列号并存储在变量 id 中。其中 {id}：指定待转换的 SRA 文件，${id} 会插入循环读取的数字部分，生成类似 SRR15417902、SRR15417903 等文件名。图片
6.数据质控和清洗数据质控查看
# 对当前文件夹中所有以fastq.gz文件进行质量控制 fastqc -t 2 -o ./ *.fastq.gz # 对当前文件夹中所有以fastq.gz文件进行整合质量控制，输出到fastq_qc文件夹中 multiqc ./*zip -o ./fastq_qcfastqc：调用 FastQC 工具，用于分析测序数据的质量。-t 2：指定使用 2 个线程并行处理文件，以加快分析速度。-o ./：指定输出目录为当前目录（./），FastQC 生成的报告文件将保存在当前目录中。*.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*ls _clean.fq.gz: 列出所有以 _clean.fq.gz 结尾的文件，即预处理过的 miRNA 序列文件。while read id; do ... done: 使用 while 循环逐个读取文件名并将其赋值给变量 id，然后对每个文件执行循环内的命令。echo $id: 打印当前正在处理的文件名，用于检查进度。bowtie -p 2 -x 图片
对比结果中发现只有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 donels *_clean.fq.gz | while read id：列出当前目录下所有文件名以 _clean.fq.gz 结尾的文件。while read id 用来逐行读取这些文件名，并将文件名存储在变量 id 中。echo $id：打印当前正在处理的文件名，以便追踪进度。bowtie -p 2 -x 图片
这个对比结果情况就勉强能“让人接受”啦~
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>&1featureCounts：调用 featureCounts 工具进行基因计数。-T 2：指定使用 2 个线程。-F gff：指定输入文件的注释格式为 GFF。-M：允许多重比对的 reads（即与多个位置比对的 reads）。-t miRNA：在 GFF 文件中，选择类型为 miRNA 的条目，这样可以仅对成熟 miRNA 计数。-g Name：选择 GFF 文件中 Name 字段作为基因 ID。-a
图片
因为比对有问题，定量也很难保证，所以拿到了矩阵也很难进行后续分析：图片
后续的分析基本上等同于转录组测序表达量矩阵，就是差异分析等统计可视化：
图片
参考资料：Multiomic analysis of microRNA-mediated regulation reveals a proliferative axis involving miR-10b in fibrolamellar carcinoma. JCI Insight. 2022 Jun 8;7(11):e154743.DNAJB1-PRKACA fusion protein-regulated LINC00473 promotes tumor growth and alters mitochondrial fitness in fibrolamellar carcinoma. PLoS Genet 2024 Mar;20(3):e1011216.Chemical, Molecular, and Single-nucleus Analysis Reveal Chondroitin Sulfate Proteoglycan Aberrancy in Fibrolamellar Carcinoma. Cancer Res Commun 2022 Jul;2(7):663-678.生信技能树：后记我确实是看完了教学视频，以及配套的笔记，但是不知道为什么结果就大相径庭，一个人学习生信就是如此的枯燥和难受！
小RNA建库测序后的数据分析-实例讲解跟着生信技能树学习microRNA测序学习构建miRNA-seq数据分析环境miRNAseq数据分析这么多年了它的流程也没有固定这5个miRNA组成的肺鳞癌诊断基因集在tcga数据库能复现吗什么，给你了你这么多miRNA靶基因查询R包和网页工具你居然不知道怎么使用对miRNA进行go和kegg等功能数据库数据库注释使用miRNAtap数据源提取miRNA的预测靶基因结果谁说Windows下无法做生信分析（植物miRNA gene预测给你看）你希望这个探针注释到蛋白编码基因还是miRNA的基因呢如果miRNA的3p和5p功能不一样miRNA、LncRNA、CircRNA靠谱小结计算MiRNA–mRNA表达相关性使用多个网页工具预测MiRNA–mRNA相互作用一篇文章学会miRNA-seq分析致谢：感谢曾老师以及生信技能树团队全体成员。
注：若对内容有疑惑或者有发现明确错误的朋友，请联系后台(欢迎交流)。
本站仅提供存储服务，所有内容均由用户发布，如发现有害或侵权内容，请点击举报。
备案号: