A Verified Workflow for MiRNA-Seq Data Processing and Bioinformatics Analysis Using R

Summary

Here, we present a protocol to analyze miRNA-Seq data using R. The workflow empowers researchers to explore miRNA regulated networks, and their significance in diverse biological and clinical questions. This work intends to serve as a practical guide for both novice and experienced researchers in the field of miRNA bioinformatics.

Abstract

MicroRNAs (miRNAs) are critical post-transcriptional regulators that influence a wide range of physiological and pathological processes. With the advancement of high-throughput sequencing technologies, miRNA-Seq has emerged as a powerful tool for profiling miRNA expression patterns. However, reliable interpretation of such data requires a standardized and reproducible analysis pipeline. Here, we present a verified workflow for miRNA-Seq data processing and bioinformatics analysis using R. This protocol encompasses all essential steps, including raw data preprocessing, quality control, alignment, quantification, normalization, differential expression analysis, target prediction, functional enrichment, and regulatory network construction. Designed for flexibility and transparency, the workflow integrates widely adopted R packages and supports species-specific annotation and modular customization. Additionally, users are guided to conduct downstream biological interpretation by leveraging curated databases and visualization tools such as Cytoscape. This protocol not only supports robust statistical analysis but also enables meaningful insights into miRNA-mRNA interactions and their roles in disease mechanisms. It is particularly well-suited for both novice and experienced researchers conducting miRNA biomarker discovery, disease modeling, or integrative multi-omics studies.

Introduction

MicroRNAs (miRNAs) are short non-coding RNA molecules that significantly influence gene expression by acting at the post-transcriptional stage¹. They typically function by binding to complementary sequences in the 3' untranslated regions (UTRs) of target messenger RNAs (mRNAs), leading to mRNA degradation or translational repression¹. Over the past two decades, miRNAs have been increasingly recognized as central regulators of various biological processes, including cell proliferation, differentiation, apoptosis, immune responses, and organ development². Moreover, dysregulation of miRNA expression has been implicated in the pathogenesis of numerous diseases, such as cancer, cardiovascular disease, neurological disorders, and kidney disease³. These findings highlight the potential of miRNAs not only as therapeutic targets but also as minimally invasive biomarkers in clinical diagnostics.

With the advent of next-generation sequencing (NGS) technologies, the study of miRNAs has entered a new era. Unlike microarray-based methods that are limited to known miRNAs, miRNA sequencing (miRNA-Seq) enables comprehensive, high-throughput, and unbiased profiling of both known and novel miRNAs across different sample types and conditions⁴. miRNA-Seq provides superior sensitivity, accuracy, and dynamic range, making it a preferred method for investigating miRNA expression patterns and discovering regulatory mechanisms in physiological and pathological settings⁵. However, the analysis of miRNA-Seq data presents specific computational challenges, including the handling of short read lengths, removal of adapter sequences, distinguishing between closely related miRNA family members, and managing high redundancy in read counts⁶. These characteristics necessitate a carefully designed and standardized analytical workflow.

Although various pipelines and software tools have been developed for miRNA-Seq data analysis, many of them rely on graphical user interfaces or fixed workflows that limit flexibility and reproducibility⁷. In contrast, the R programming environment provides a powerful and customizable platform for bioinformatics analysis⁸. R offers a rich ecosystem of packages for statistical modeling, data visualization, and integration with biological databases. This enables users to perform comprehensive and reproducible analyses in a transparent and script-based manner. Moreover, the modular nature of R workflows allows researchers to tailor each step according to specific experimental requirements, from raw data preprocessing to functional interpretation.

In this protocol, we present a verified and complete miRNA-Seq analysis workflow implemented entirely in R, with the goal of providing a reproducible and user-adaptable solution for researchers working with miRNA expression data. The workflow begins with quality control and adapter trimming of raw sequencing reads, followed by alignment to a reference genome or known miRNA sequences. Subsequent steps include quantification of read counts, normalization, differential expression analysis, target gene prediction, functional enrichment, and network visualization. The workflow incorporates several widely used and well-maintained R packages, ensuring both reliability and compatibility with future updates and extensions.

One of the core strengths of this protocol lies in its ability to move beyond differential expression results and provide meaningful biological interpretation. By integrating curated databases of validated and predicted miRNA-mRNA interactions, the workflow allows users to identify biologically relevant target genes. These targets can then be subjected to Gene Ontology and pathway enrichment analyses to uncover affected biological processes and molecular pathways. In the final step, miRNA-mRNA interaction networks can be visualized using external tools such as Cytoscape⁹, providing insight into the regulatory landscape and identifying key hub miRNAs with potential functional importance.

This method has been successfully applied in clinical research contexts, including studies on kidney disease, where circulating miRNAs serve as promising biomarkers for diagnosis and prognosis¹⁰. However, the modular and flexible design of the workflow makes it suitable for a wide range of applications, including disease modeling, drug response studies, developmental biology, and comparative genomics. Researchers can easily adapt the workflow to accommodate species-specific annotations, experimental conditions, or additional layers of omics data.

By offering an open-source, script-based solution, this R-centered pipeline addresses several of the major limitations associated with existing miRNA-Seq tools, including limited customization, reliance on non-transparent graphical interfaces, lack of support for non-model organisms, poor reproducibility due to the absence of version control, and difficulty in integrating with downstream statistical and functional analysis frameworks. It enables full control over data processing parameters, encourages reproducibility through version-controlled code, and promotes transparency in bioinformatics research. As the importance of miRNAs continues to grow in the context of systems biology and translational medicine, having access to a reliable and adaptable analysis framework becomes increasingly essential.

Protocol

NOTE: Materials with software links are listed in the Table of Materials.

1. Prepare RNA samples and sequence libraries

NOTE: Perform RNA extraction and sequencing outside this computational workflow. There is more than one way to analyze miRNA-sequencing data. This section provides the context of one practical.

1. Extract total RNA: Extract total RNA from the biological samples using a kit optimized for small RNA isolation (e.g., miRNA isolation kit). Follow the manufacturer's protocol carefully. Ensure that you use RNase-free consumables and keep samples on ice to minimize degradation.
2. Assess RNA integrity and quantity: Run 1-2 µL of the extracted RNA on a Bioanalyzer or equivalent device. Check the RNA integrity number (RIN) and ensure it is ≥ 7.0 for reliable sequencing. Record the concentration using a spectrophotometer or fluorometer.
3. Construct small RNA libraries: Prepare sequencing libraries from 1 µg of total RNA using a commercial small RNA-seq library preparation kit. Follow the kit protocol to ligate adapters, reverse transcribe, and amplify the cDNA. Purify the PCR products with size selection (e.g., 18-30 nt insert) to enrich for miRNA fragments.
4. Sequence libraries: Load the libraries onto a high-throughput sequencing platform. Set the run configuration for single-end sequencing with a read length of ~50 bp. Ensure that each sample generates approximately 10 million raw reads to achieve sufficient depth.
5. Export sequencing output: After sequencing, export the raw data as FASTQ files using the instrument's data output software. Verify that the output directory contains both the sequence reads and the associated quality score files. Store the FASTQ files in a structured directory for downstream analysis.

2. Preprocess raw reads and perform quality control

1. Trim adapter sequences
   1. Install and configure Cutadapt or fastp.
   2. Run adapter trimming on each FASTQ file using the following command:
      cutadapt -a XXXX -o trimmed_reads.fastq raw_reads.fastq
      NOTE: The adapter sequence followed after '-a' should apply with the specific sRNA library preparation kit. '-o' defines the output filename with the input file followed.
2. Assess read quality
   1. Use FastQC to generate quality control reports:
      fastqc trimmed_reads.fastq
   2. Review per-base quality scores, read length distribution, and adapter contamination: Open the FastQC HTML report generated for each FASTQ file in a web browser. Examine the following modules step by step:
      1. Per-base sequence quality: Ensure that most bases fall within the green zone (Phred score ≥30). Look for any drop in quality at the 3′ end that might indicate sequencing errors.
      2. Read length distribution: Confirm that the distribution corresponds to the expected insert size (e.g., 18-30 nt for miRNAs). Check that no unexpected peaks are present.
      3. Adapter content: Verify that adapter sequences have been effectively trimmed. Confirm that the adapter contamination percentage is close to zero after trimming.
   3. Save the FastQC summary report and flag any samples with poor quality metrics for re-trimming or exclusion from further analysis.

3. Map reads and generate count matrices

1. Align reads to reference
   1. Download the FASTA file of the reference genome or mature miRNA sequences (e.g., from miRBase)¹¹. Example:
      wget ftp://mirbase.org/pub/mirbase/CURRENT/mature.fa
   2. Index the reference genome using Bowtie.
      1. Open a terminal and run the following command to build the index:
         bowtie-build reference.fa reference_index
      2. Replace reference.fa with the actual FASTA file name.
      3. Replace reference_index with the desired prefix for the index.
      4. Ensure that Bowtie generates multiple index files (e.g., .ebwt). Verify that these files are present in the working directory, as they are required for alignment.
   3. Align reads using Bowtie with suitable parameters for short reads. Example:
      bowtie -v 0 -a --best --strata reference_index trimmed_reads.fastq > aligned_reads.sam
      NOTE: The input file is trimmed_reads.fastq, and the output file is aligned_reads.sam. '-v 0' means we allow no mistake in the entire read. '-a -best -strata' means we throw away any alignment that has more mismatches than the best.
2. Quantify miRNA expression
   1. Convert SAM files to BAM format using SAMtools.
      samtools view -S -b aligned_reads.sam > aligned_reads.bam
      NOTE: The input file aligned_reads.sam is the result from the last command. The output file aligned_reads.bam is prepared for subsequent analysis.
      1. Use SAMtools to compress and sort the alignment file:
         samtools sort aligned_reads.bam -o aligned_reads_sorted.bam
samtools index aligned_reads_sorted.bam
      2. Ensure the first command converts the SAM file into BAM format.
      3. Ensure the second command sorts the BAM file by genomic coordinates.
      4. Ensure the third command generates an index file (.bai), which is required for downstream analyses.
      5. Confirm that both the sorted BAM file and its index are successfully created before proceeding to quantification.
   2. Use featureCounts or HTSeq-count to generate a count matrix using miRNA annotation GTF:
      featureCounts -a miRNA.gtf -o counts.txt aligned_reads.bam
      ​NOTE: featureCounts quantifies reads in the input file aligned_reads.bam based on miRNA.gtf, and outputs counts.txt.

4. Perform Differential Expression Analysis in R

1. Load Count Data
   1. Import count matrix and sample metadata into R:
      library(DESeq2)
countData <- read.csv("counts.csv", row.names=1)
colData <- read.csv("metadata.csv", row.names=1)
dds <- DESeqDataSetFromMatrix(countData = countData, colData = colData, design = ~ condition)
      NOTE: Count data (counts.csv) and group of samples (metadata.csv) need to provided to DESeq2. 'condition' here clarifies the group of samples provided. For individual requirements, reference Manual of DESeq2¹².
2. Normalize and transform data
   1. Normalize count data using DESeq2's default method:
      dds <- DESeq(dds)
   2. Perform variance-stabilizing transformation:
      vsd <- vst(dds, blind=FALSE)
   3. Visualize sample clustering using PCA:
      plotPCA(vsd, intgroup="condition")
3. Identify differentially expressed miRNAs
   1. Extract and sort differential expression results:
      res <- results(dds)
resOrdered <- res[order(res$pvalue), ]
      NOTE: We reorder the results file 'res' based on the value of pvalue.
      summary(res)
   2. Filter for significantly differentially expressed miRNAs (p-value < 0.05, |log2FC| > 1).
      sig_miRNA <- subset(res, pvalue < 0.05 & abs(log2FC) > 1)
      NOTE: There are multiple thresholds for filtering significantly changed miRNAs. The requirements 'p-value < 0.05, |log2FC| > 1' are widely applied. The thresholds can be adjusted for individual data.
4. Visualize expression changes
   1. Install and load EnhancedVolcano package.
   2. Create the volcano plot:
      library(EnhancedVolcano)
EnhancedVolcano(res,
lab = rownames(res),
x = 'log2FoldChange',
y = 'pvalue',
title = 'Differentially Expressed miRNAs')

5. Predict target genes of miRNAs

1. Query databases
   1. Use online resources such as TargetScan, miRDB, and miRTarBase to search for specific microRNA and retrieve target genes.
   2. Focus on experimentally validated targets for greater confidence.
2. Automate prediction in R
   1. Load multiMiR package and query validated targets:
      library(multiMiR)
target_results <- get_multimir(mirna = c("hsa-miR-21-5p"), table = "validated")
      NOTE: Here we take "hsa-miR-21-5p" for example and get validated targets.
   2. Extract unique target gene symbols for enrichment analysis:
      genes <- unique(target_results@data$target_symbol)

6. Conduct functional enrichment analysis

1. Perform GO enrichment
   1. Load enrichment tools:
      library(clusterProfiler)
library(org.Hs.eg.db)
      NOTE: Here, we load a database containing human genome annotations, which are useful for converting common gene identifiers.
   2. Run GO enrichment analysis for Biological Processes:
      ego <- enrichGO(gene = genes,
OrgDb = org.Hs.eg.db,
keyType = "SYMBOL",
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = 0.05)
dotplot(ego)
      NOTE: For the usage of enrichGO, we provide a list of genes, and clarify they are nominated as 'SYMBOL' here. We conducted biological process enrichment, which corresponds to parameter 'ont ="BP" '. We perform multiple testing correction, so we designate pAdjustMethod = "BH". For a significant threshold, we choose pvalueCutoff = 0.05.dotplot presents a visualized result. For more individual options, reference the Manual of clusterProfiler¹³.
2. Perform KEGG pathway enrichment
   1. Run KEGG enrichment:
      ekegg <- enrichKEGG(gene = genes, organism = 'hsa')
dotplot(ekegg)
      NOTE: For the usage of enrichKEGG, we provide a list of genes, and clarify the organism as human ("has"). dotplot presents a visualized result. For more individual options, reference the Manual of clusterProfiler¹³.

7. Construct and visualize miRNA - mRNA interaction network

1. Export data for network visualization
2. Create a data frame of miRNA-target gene pairs based on the targets generated from TargetScan, miRDB or miRTarBase.
3. Write the network table to CSV:
   write.csv(miRNA_target_pairs, "network.csv")
4. Import into Cytoscape
   1. Open Cytoscape and import the network table.
   2. Visualize the network using force-directed or circular layout.
   3. Analyze topological properties (e.g., degree centrality) to identify hub miRNAs.

Results

We downloaded the microRNA expression matrix from GSE133530 and conducted differential expression analysis directly. We provided an example analytical R script for the dataset in Supplementary File 1. The dataset performed global miRNA profiling on 16 renal cysts of different sizes (minimal cysts: less than 1-5 mL, n = 10; medium cysts: between 10-25 mL, n = 4; large cysts: greater than 50 mL, n = 4) and minimally cystic tissue (MCT, n = 7, including 1 replicate) from fou...

Discussion

The analysis of miRNA-Seq data presents distinct challenges due to the small size and redundancy of reads, making rigorous quality control and preprocessing critical. One of the most important steps in the workflow is adapter trimming. Because miRNAs are approximately 22 nucleotides long, adapter sequences can easily dominate the reads if not properly removed. Failure to perform accurate trimming can result in misalignment and inflation of false-positive reads. Similarly, quality filtering should be implemented before al...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
Agilent-021827 Human miRNA Microarray	Agilent	/	A commercial array for profiling microRNAs of human samples
Bowtie	Johns Hopkins University	http://bowtie-bio.sourceforge.net/index.shtml	A software tool for aligning sequencing reads to long reference sequences
clusterProfiler (R package)	Bioconductor	https://bioconductor.org/packages/clusterProfiler/	An R package designed for functional enrichment analysis and visualization of high-throughput biological data.
Cutadapt	Open Source	https://cutadapt.readthedocs.io	A command-line tool that removes adapter sequences, primers, poly-A tails and other unwanted fragments from high-throughput sequencing reads.
Cytoscape	Cytoscape Consortium	https://cytoscape.org/	An open-source software platform designed for the visualization and analysis of complex biological networks.
DESeq2 (R package)	Bioconductor	https://bioconductor.org/packages/DESeq2/	An R package designed for differential gene expression analysis of count data
EnhancedVolcano (R package)	Bioconductor	https://bioconductor.org/packages/EnhancedVolcano/	An R package designed to create publication-quality volcano plots.
FastQC	Babraham Bioinformatics	https://www.bioinformatics.babraham.ac.uk/projects/fastqc/	An open-source quality-control tool for high-throughput sequencing data.
featureCounts	Subread / SourceForge	http://subread.sourceforge.net/	A program used for counting reads mapped to genomic features
HTSeq-count	Python Package	https://htseq.readthedocs.io	A command-line tool that counts how many aligned high-throughput sequencing reads overlap genomic features such as genes or exons. I
Illumina Human v2 MicroRNA expression beadchip	Illumina	/	A commercial array for profiling microRNAs of human samples
multiMiR (R package)	Bioconductor	https://bioconductor.org/packages/multiMiR/	An R package that provides the largest integrated collection of predicted and experimentally validated microRNA–target interactions together with their associations to diseases and drugs.
org.Hs.eg.db (R package)	Bioconductor	https://bioconductor.org/packages/org.Hs.eg.db/	An annotation package designed for human (Homo sapiens) genomics research.
R Software	R Project	https://www.r-project.org/	An open-source Project for Statistical Computing
Rstudio	Posit PBC	/	An integrated development environment helps to be more productive with R and Python
SAMtools	Open Source	http://www.htslib.org/	A software package for manipulating next-generation sequencing (NGS) data.