miRNA sequencing

Pre-alignment

Make directories

mkdir -p miRNA1224/mirge3/miRge3_Lib
mkdir -p miRNA1224/mirge3/alignment
mkdir -p miRNA1224/mirge3/scripts

Merge fastq.gz

Merge files from different lanes and runs for one sample (L1, L2, run2 etc).

#go to folder with raw files in it
cd home1/gdfenglab/akf010/miRNA1224/miRNA_rawfiles

#concatenate the files for each sample
cat Mock_1_L1_run1.fastq.gz Mock_1_L1_run2.fastq.gz Mock_1_L2_run2.fastq.gz > Mock1_merged.fastq.gz

cat Mock_2_L1_run1.fastq.gz Mock_2_L1_run2.fastq.gz Mock_2_L2_run2.fastq.gz > Mock2_merged.fastq.gz

cat TR_D7_1_L1_run1.fastq.gz TR_D7_1_L1_run2.fastq.gz  TR_D7_1_L2_run2.fastq.gz > TR1_merged.fastq.gz

cat TR_D7_2_L1_run1.fastq.gz TR_D7_2_L1_run2.fastq.gz  TR_D7_2_L2_run2.fastq.gz > TR2_merged.fastq.gz

Conda environment

name: mirge3_env
channels:
  - conda-forge
  - bioconda
  - defaults
dependencies:
  - python=3.8
  - mirge3
  - gffread
  - miranda
prefix: /home/gdfenglab/akf010/.conda/envs/mirge3_env

conda create --name mirge3_env python=3.8
conda activate mirge3_env
conda install mirge3 gffread miranda

Alignment

Setting up mirge3 library

cd miRge_Lib
wget -O human.tar.gz "https://sourceforge.net/projects/mirge3/files/miRge3_Lib/human.tar.gz/download"
tar -xzf human.tar.gz
cd ..

Running mirge3

Putting all four samples in the command so I can include differential gene expression in the output.

Made metadata file with the sample information

sample,group
Mock1_merged,control
Mock2_merged,control
TR1_merged,treatment
TR2_merged,treatment

Call the script

sbatch mirge3.all

#!/bin/bash

#SBATCH --mail-type=NONE
#SBATCH --mail-user=alison.francois@nationwidechildrens.org
#SBATCH --job-name=mirge3
#SBATCH --time=2:00:00
#SBATCH --ntasks=4
#SBATCH --output=R-%x.%j.out
#SBATCH --error=R-%x.%j.err

# conda environment set up before this
# load conda
# make sure correct directory

. /export/apps/opt/Miniconda3/23.10.0/etc/profile.d/conda.sh
conda activate mirge3_env
echo "Conda env is mirge3_env"

lib=/home/gdfenglab/akf010/miRNA1224/mirge3/miRge3_Lib
inDir=/home/gdfenglab/akf010/miRNA1224/miRNA_rawfiles

mkdir -p /home/gdfenglab/akf010/miRNA1224/mirge3/alignment/all_together
outDir=/home/gdfenglab/akf010/miRNA1224/mirge3/alignment/all_together

# ADAPTER IS QIAGEN, NOT ILLUMINA

miRge3.0 -s ${inDir}/Mock1_merged.fastq.gz,${inDir}/Mock2_merged.fastq.gz,${inDir}/TR1_merged.fastq.gz,${inDir}/TR2_merged.fastq.gz \
         -lib ${lib} \
         -on human \
         -db mirgenedb \
         -o "${outDir}" \
         -cpu 8 \
         -a AACTGTAGGCACCATCAAT \
         -dex \
         -mdt ${inDir}/metadata.csv

conda deactivate

echo "all done!"

Output appears in its own dated folder

tree miRNA1224/mirge3/alignment/all_together/miRge.2026-04-04_10-36-45

../mirge3/alignment/all_together/
└── miRge.2026-04-04_10-36-45
    ├── annotation.report.csv # gives you a table categorized by RNA type
    ├── annotation.report.html # same thing but opens in your browser
    ├── index_data.js
    ├── mapped.csv # all mapped reads 
    ├── miR.Counts.csv  # the raw count data
    ├── miR.DESeq2.RData
    ├── miR.DiffExpn.txt #DGE results
    ├── miRge3_visualization.html # Open this in browser to look through results
    ├── mirge_diffex.pca.pdf
    ├── mirge_diffex.volcano.pdf
    ├── miR.RPM.csv # This one is NORMALIZED, use for later analysis.
    ├── run.log # details all the steps mirge3 went through
    └── unmapped.csv # the unmapped reads (leftovers)

I pasted miR.RPM.csv into excel and made a sortable table of my data. I don’t have a paid license to download graphs directly from the html page but I took screenshots. I could still export a data table from each graph on the html page.

Aligning to CMV transcripts

Create FASTA files from miRge3 output

We have to convert the .csv files into fasta files we can then run against the HCMV transcriptome. We are also filtering out reads where none of the samples showed counts above 5 RPM.

Call the scripts:

sbatch convert.csv # for the unmapped reads CSV
sbatch convert.mapped # same process but for the mapped reads file

#!/bin/bash

#SBATCH --mail-type=NONE
#SBATCH --mail-user=alison.francois@nationwidechildrens.org
#SBATCH --job-name=convert_csv
#SBATCH --time=1:00:00
#SBATCH --cpus-per-task=1
#SBATCH --output=R-%x.%j.out
#SBATCH --error=R-%x.%j.err

# conda environment set up before this
# load conda
# make sure correct directory

. /export/apps/opt/Miniconda3/23.10.0/etc/profile.d/conda.sh
conda activate mirge3_env
echo "Conda env is mirge3_env"

python << 'EOF'
import pandas as pd

def convert_unmapped_to_fasta(input_csv, output_fasta, min_count=5):
    # Load miRge3 unmapped reads
    df = pd.read_csv(input_csv)

    # 1. Identify your sample columns (usually columns 1 through 4)
    # iloc[:, 1:5] selects the 2nd, 3rd, 4th, and 5th columns
    sample_cols = df.columns[11:14]

    # 2. Filter: Keep row if ANY sample column value >= min_count
    # .gt(min_count-1) checks if value is >= threshold
    mask = df[sample_cols].gt(min_count - 1).any(axis=1)
    filtered_df = df[mask]

    # 3. Write to FASTA
    with open(output_fasta, 'w') as f:
        for i, row in filtered_df.iterrows():
            # Convert RNA sequence to DNA for miRanda
            sequence = str(row.iloc[0]).replace('U', 'T')

            # Get the max count among samples for the header
            max_val = row[sample_cols].max()

            f.write(f">read_{i}_maxcount_{max_val}\n{sequence}\n")

    print(f"Success! Created {output_fasta} with {len(filtered_df)} reads")
# Run the function
# Change 'unmapped.csv' to your actual file path
convert_unmapped_to_fasta('/home/gdfenglab/akf010/miRNA1224/mirge3/alignment/all_together/miRge.2026-04-04_10-36-45/unmapped.csv', '/home/gdfenglab/akf010/miRNA1224/mirge3/alignment/all_together/miRge.2026-04-04_10-36-45/filtered_unmapped.fasta', min_count=5)

EOF

conda deactivate

#!/bin/bash

#SBATCH --mail-type=NONE
#SBATCH --mail-user=alison.francois@nationwidechildrens.org
#SBATCH --job-name=convert_csv
#SBATCH --time=1:00:00
#SBATCH --cpus-per-task=1
#SBATCH --output=R-%x.%j.out
#SBATCH --error=R-%x.%j.err

# conda environment set up before this
# load conda
# make sure correct directory

. /export/apps/opt/Miniconda3/23.10.0/etc/profile.d/conda.sh
conda activate mirge3_env
echo "Conda env is mirge3_env"

python << 'EOF'
import pandas as pd

def convert_mapped_to_fasta(input_csv, output_fasta, min_count=5):
    # Load miRge3 mapped reads
    df = pd.read_csv(input_csv)

    # 1. Identify your sample columns (usually columns 1 through 4)
    # iloc[:, 1:5] selects the 2nd, 3rd, 4th, and 5th columns
    sample_cols = df.columns[11:14]

    # 2. Filter: Keep row if ANY sample column value >= min_count
    # .gt(min_count-1) checks if value is >= threshold
    mask = df[sample_cols].gt(min_count - 1).any(axis=1)
    filtered_df = df[mask]

    # 3. Write to FASTA
    with open(output_fasta, 'w') as f:
        for i, row in filtered_df.iterrows():
            # Convert RNA sequence to DNA for miRanda
            sequence = str(row.iloc[0]).replace('U', 'T')

            # Get the max count among samples for the header
            max_val = row[sample_cols].max()

            f.write(f">read_{i}_maxcount_{max_val}\n{sequence}\n")

    print(f"Success! Created {output_fasta} with {len(filtered_df)} reads")
# Run the function
# Change 'mapped.csv' to your actual file path
convert_mapped_to_fasta('/home/gdfenglab/akf010/miRNA1224/mirge3/alignment/all_together/miRge.2026-04-04_10-36-45/mapped.csv', '/home/gdfenglab/akf010/miRNA1224/mirge3/alignment/miRanda/filtered_mapped.fasta', min_count=5)

EOF

conda deactivate

We now have fasta files with the unmapped and mapped reads, filtered, separately.

Getting HCMV transcript fasta file

This is the website I used to get a TR strain genome assembly.

https://www.ncbi.nlm.nih.gov/datasets/genome/GCA_027927695.1/

In the reference folder I downloaded the following files:

cd home1/gdfenglab/akf010/miRNA1224/reference

#download from genome assembly
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/027/927/695/GCA_027927695.1_ASM2792769v1/GCA_027927695.1_ASM2792769v1_rna_from_genomic.fna.gz

#unzip the file
gunzip GCA_027927695.1_ASM2792769v1_rna_from_genomic.fna.gz 

#rename it
mv GCA_027927695.1_ASM2792769v1_rna_from_genomic.fna TR_transcripts.fasta

Miranda alignment

MiRanda lets you align a fasta of your collapsed reads (unmapped or mapped) to the viral transcripts and see if anything could be a match.

Call the script:

sbatch miranda.batch

#!/bin/bash

#SBATCH --mail-type=NONE
#SBATCH --mail-user=alison.francois@nationwidechildrens.org
#SBATCH --job-name=miRanda
#SBATCH --time=4:00:00
#SBATCH --ntasks=8
#SBATCH --nodes=2-4
#SBATCH --output=R-%x.%j.out
#SBATCH --error=R-%x.%j.err

# conda environment set up before this
# load conda
# make sure correct directory

. /export/apps/opt/Miniconda3/23.10.0/etc/profile.d/conda.sh
conda activate mirge3_env
echo "Conda env is mirge3_env"

genomeDir=/home/gdfenglab/akf010/miRNA1224/reference
inDir=/home/gdfenglab/akf010/miRNA1224/mirge3/alignment/miRanda
outDir=/home/gdfenglab/akf010/miRNA1224/mirge3/alignment/miRanda

miranda ${inDir}/filtered_unmapped.fasta ${genomeDir}/TR_transcripts.fasta -sc 140 -en -20 -strict -out ${outDir}/unmapped_miranda_results.txt
#get a summary
grep -e ">>" ${inDir}/unmapped_miranda_results.txt | awk '{print $1"\t"$2"\t"$3"\t"$4}' > ${outDir}/unmapped_miranda_summary.tsv

#######

miranda ${inDir}/filtered_mapped.fasta ${genomeDir}/TR_transcripts.fasta -sc 140 -en -20 -strict -out ${outDir}/mapped_miranda_results.txt
#get a summary
grep -e ">>" ${inDir}/mapped_miranda_results.txt | awk '{print $1"\t"$2"\t"$3"\t"$4}' > ${outDir}/mapped_miranda_summary.tsv

conda deactivate