nanornaqtl 2.0.0

A comprehensive toolkit for multi-molecular phenotyping and QTL analysis based on Nanopore direct RNA sequencing data

nanornaqtl

A comprehensive toolkit for molecular phenotyping and QTL analysis using Nanopore direct RNA sequencing data

Introduction

nanornaqtl is an analysis tool specifically designed for Nanopore direct RNA sequencing data, supporting both standard Nanopore direct RNA-seq (poly(A)-enriched) and Nanopore direct total RNA-seq (covering both mRNA and non-coding RNA). It is used to identify multiple RNA molecular phenotypes and perform population-level QTL (Quantitative Trait Loci) analysis. This tool is particularly suitable for analyzing pooled samples (sequencing of multiple mixed individuals), with all analyses performed at the read level.

Main Features

• 12 Types of Molecular Phenotype Identification 📊:

  • RNA modification sites: m6A, m5C, pseudouridine (pseU), inosine, 2'-O-methylation (2OmeA / 2OmeC / 2OmeU / 2OmeG)
  • poly(A) tail length and alternative polyadenylation site (APA)
  • Intron retention rate
  • Transcript isoform

• 12 Types of QTL Analysis 🧪:

  • Modification QTL (m6A, m5C, pseU, inosine, 2'-O-methylation)
  • APA usage pattern QTL
  • Isoform usage pattern QTL
  • poly(A) tail length QTL
  • Intron retention rate QTL

Features ✨

• Support for both direct RNA-seq and direct total RNA-seq: Optional RNA type classification (mRNA / ncRNA / unknown) for total RNA data
• Read-level analysis: Suitable for population-level analysis of pooled samples
• Bayesian statistical methods: More robust for uneven coverage in Nanopore data
• 5' end-based read integrity assessment: Optional degradation detection to flag truncated reads
• Parallel processing: Supports multi-threading acceleration with chromosome-based parallel processing
• Flexible parameter settings: Customizable quality thresholds, coverage requirements, etc.

---

Installation

Install via pip 📦

pip install nanornaqtl

Install from source

git clone https://github.com/xinranxu0930/nanornaqtl.git
cd nanornaqtl
poetry install

Dependencies

External tools:

• bedtools (v2.30.0+)
• samtools (v1.15+)
• IsoQuant (for isoform analysis): https://github.com/ablab/IsoQuant

Python packages:

• pysam
• pymc
• arviz
• statsmodels
• pandas
• numpy
• scipy
• logomaker
• biopython

---

Data Preprocessing

Before using nanornaqtl, Nanopore sequencing data needs to be preprocessed.

Step 1: Data Format Conversion

If the sequencing data is in FAST5 format, it needs to be converted to POD5 format first (skip this step if already in POD5):

# Convert using pod5 tool
pod5 convert fast5 input.fast5 --output output.pod5

Step 2: Basecalling and Alignment 🔧

Use Dorado for basecalling and alignment:

Basic Alignment (Required)

dorado basecaller \
  model_path \
  pod5_directory/ \
  --reference fasta.fa \
  --mm2-opts "-x splice -k 14" \
  > output.bam

Parameter explanation:

• --mm2-opts "-x splice -k 14": RNA splice alignment parameters

Optional Features

If modification site analysis is needed:

dorado basecaller \
  model_path \
  pod5_directory/ \
  --reference fasta.fa \
  --mm2-opts "-x splice -k 14" \
  --modified-bases-models model_name \
  > output.bam

Download modification recognition models:

# Download all models
dorado download --model all

# Or download specific model
dorado download --model dna_r10.4.1_e8.2_400bps_hac@v4.1.0

If APA or polyA tail length analysis is needed:

dorado basecaller \
  model_path \
  pod5_directory/ \
  --reference fasta.fa \
  --mm2-opts "-x splice -k 14" \
  --estimate-poly-a \
  > output.bam

---

Detailed Usage Instructions

1. prep - Data Preprocessing

Function: Filter mapped reads from Dorado basecall BAM. Supports MAPQ / read length / sequencing quality filtering, optional RNA type classification (mRNA / ncRNA / unknown) based on gene annotation, and optional 5' end-based degradation detection.

Input:

• Raw BAM file output from Dorado (no preprocessing required)

Command:

nanornaqtl prep \
  -b <basecall_bam> \
  -p <output_prefix> \
  -t <threads> \
  -q <min_mapq> \
  -l <min_read_length> \
  --min_qscore <min_qscore> \
  [--RNA_classification -g <gtf>] \
  [--filter_degraded --tss_window 50 --tss_min_count 3]

Parameter explanation:

Parameter	Description	Default
-b, --basecall_bam	Input Dorado basecall BAM file	Required
-p, --dir_pre	Output prefix (can include path)	Required
-t, --threads	Number of threads	4
-q, --min_mapq	Minimum read mapping quality (MAPQ) score	0
-l, --min_read_length	Minimum read length	0
--min_qscore	Minimum read-level sequencing quality score (from Dorado qs tag)	0
--RNA_classification	Enable RNA type classification (requires -g)	off
-g, --gtf	Gene annotation file (GTF/GFF) for RNA type classification	None
--filter_degraded	[BETA] Enable 5' end-based degradation analysis	off
--tss_window	Window size for peak detection (bp)	50
--tss_min_count	Minimum fixed read count for peak calling	3

Output files:

Base output (always generated):

• <prefix>_map.sorted.bam / .bai: All mapped reads
• <prefix>_map[0/16].sorted.bam / .bai: Strand-separated BAM (flag=0 / flag=16)
• <prefix>_map.fastq: FASTQ of mapped reads

With --RNA_classification enabled (additional 7 BAMs):

• <prefix>_mRNA_map[0/16].sorted.bam / .bai
• <prefix>_mRNA_map.sorted.bam / .bai
• <prefix>_ncRNA_map[0/16].sorted.bam / .bai (with RT:Z:<biotype> tag)
• <prefix>_ncRNA_map.sorted.bam / .bai
• <prefix>_unknown_map.sorted.bam / .bai

With --filter_degraded enabled (additional 3 TSV):

• <prefix>_degradation_detail.tsv: Per-read 5' end label
• <prefix>_degradation_peaks.tsv: Peak summary
• <prefix>_degraded_read_list.txt: Non-intact read IDs

---

2. pheno - Molecular Phenotype Identification

2.1 m6A Modification Site Identification

Function: Identify m6A (N6-methyladenosine) modification sites.

Command:

nanornaqtl pheno m6A \
  -b <map_bam> \
  -o <output_dir> \
  -p <output_prefix> \
  -t <threads> \
  -f <fasta.fa> \
  [--motif] \
  [--metaPlotR]

Parameter explanation:

Parameter	Description	Default
-b, --bam	map.bam file path from prep step (*_map.sorted.bam)	Required
-o, --output_dir	Output file directory, will be created if not exists	Required
-p, --prefix	Output file prefix name	Required
-t, --threads	Number of threads	4
--mod_threshold	Modification probability threshold	0.75
--min_qscore	Minimum base quality score for modification site	10
--min_mapq	Minimum read mapping quality (MAPQ)	0
-f, --fasta	Reference genome fasta file path	Required
--min_rate	Minimum modification rate	0.1
--min_cov	Minimum coverage for modification site	5
--motif	Enable motif filtering (DRACH: [GAT][GA]AC[ATC])	False
--motifPaint	Enable motif logo plot generation (one plot per motif subtype when applicable)	False
--metaPlotR	Generate bed file for metaPlotR	False

Output files:

• <prefix>_m6A_sites_result.csv: m6A modification site information

Output column description:

Column name	Description
chrom	Chromosome
pos_1base	Position (1-based)
strand	Strand direction (+/-)
mod_num	Number of modified reads
cov	Total coverage
mod_rate	Modification rate
motif	Motif sequence (if --motif enabled)

---

2.2 m5C Modification Site Identification

Function: Identify m5C (5-methylcytosine) modification sites.

Command:

nanornaqtl pheno m5C \
  -b <map_bam> \
  -p <output_prefix> \
  -o <output_dir> \
  -t <threads> \
  -f <fasta.fa> \
  [--motif] \
  [--metaPlotR]

Parameters: Same as m6A

Motif classification (if --motif enabled):

• NSUN6: CTCCA motif
• NSUN2: GC-rich context surrounding the modification site

Output files:

• <prefix>_m5C_sites_result.csv: m5C modification site information

Output column description: Same as m6A, with additional motif_classification column (NSUN6/NSUN2)

---

2.3 pseudouridine (pseU) Modification Site Identification

Function: Identify pseudouridine modification sites.

Command:

nanornaqtl pheno pseU \
  -b <map_bam> \
  -p <output_prefix> \
  -o <output_dir> \
  -t <threads> \
  -f <fasta.fa> \
  [--motif] \
  [--metaPlotR]

Parameters: Same as m6A

Motif classification (if --motif enabled):

• TRUB1: GUΨCN
• PUS7_type1: UNΨAR (primary pattern)
• PUS7_type2: UGΨKG (secondary pattern)

Output files:

• <prefix>_pseU_sites_result.csv: pseU modification site information

Output column description: Same as m5C, includes motif classification (TRUB1/PUS7_type1/PUS7_type2)

---

2.4 inosine Modification Site Identification

Function: Identify inosine (A-to-I editing) sites.

Command:

nanornaqtl pheno inosine \
  -b <map_bam> \
  -p <output_prefix> \
  -o <output_dir> \
  -f <fasta.fa> \
  -t <threads> \
  [--motif] \
  [--metaPlotR]

Parameters: Same as m6A

Motif (if --motif enabled): HAG (H = A/C/U)

Output files:

• <prefix>_inosine_sites_result.csv: inosine site information

---

2.5 2'-O-methylation Modification Site Identification

Function: Identify 2'-O-methylation modification sites (2OmeA, 2OmeC, 2OmeU, 2OmeG).

Command:

nanornaqtl pheno <2OmeA|2OmeC|2OmeU|2OmeG> \
  -b <map_bam> \
  -p <output_prefix> \
  -o <output_dir> \
  -t <threads> \
  -f <fasta.fa> \
  [--metaPlotR]

Note: 2'-O-methylation only supports base mode. The --motif and --motifPaint options are not applicable.

Parameters: Same as m6A (excluding --motif and --motifPaint).

Output files:

• <prefix>_<modType>_sites_result.csv: 2'O-methylation modification site information

Output column description: chrom, pos_1base, strand, mod_num, cov, mod_rate (no motif column).

---

2.6 poly(A) Tail Length and APA Site Identification

Function: Identify poly(A) tail length and alternative polyadenylation (APA) cleavage sites for each read, based on Dorado's poly(A) signal estimation.

Prerequisites: BAM file must contain pt tag (generated by Dorado basecall with --estimate-poly-a).

Note: Results from this step are expected to be less stringent than those from curated APA databases. Dorado's poly(A) detection is signal-based and does not distinguish a true post-transcriptional poly(A) tail from a genome-encoded A-rich region (internal priming). This step applies a downstream A-stretch filter to remove obvious internal priming, but the hexamer detection rate is still expected to be lower than database-level estimates (~70%), since nanopore DRS data has not undergone the additional curation steps used by databases (oligo(dT) enrichment, 3' UTR restriction, multi-sample cross-validation, etc.).

Command:

nanornaqtl pheno polyA \
  -b <map_bam> \
  -p <output_prefix> \
  -o <output_dir> \
  -f <fasta.fa> \
  -t <threads> \
  [--min_mapq <min_mapq>]

Parameter explanation:

Parameter	Description	Default
-b, --bam	map.bam file path from prep step (*_map.sorted.bam)	Required
-p, --prefix	Output file prefix name	Required
-o, --output_dir	Output directory path, will be created if not exists	Required
-f, --fasta	Reference genome fasta file (must be indexed with .fai)	Required
-t, --threads	Number of threads (max 44)	4
--min_mapq	Minimum read mapping quality (MAPQ)	0

Analysis logic:

Starting from Dorado's pt tag, this step performs the following:

1. poly(A) length filtering. Each read is classified by its poly(A) length:
   • no_polyA (< 3 nt): essentially no poly(A) signal at the 3' end.
   • short_oligoA (3–14 nt): likely a short oligo(A) tag added by the TRAMP complex to mark RNA for exosome-mediated degradation, rather than a stable poly(A) tail.
   • Pass (≥ 15 nt): considered a true poly(A) tail. The 15 nt threshold is a balance — PABP requires roughly 20 nt to remain bound, while existing studies have reported tails as short as 8 nt.
2. Internal priming filtering. For reads passing the length filter, the genomic sequence within 20 nt downstream of the cleavage site is inspected. Reads with ≥ 6 consecutive A's downstream are flagged as internal_priming and excluded.
3. Cleavage site clustering. Cleavage sites within 24 nt of each other are merged into a single PAS cluster.
4. Motif characterization (per cluster, on the reference genome):
   • PAS hexamer (18 known motifs) is searched 0–60 nt upstream of the representative cleavage site, with priority and distance jointly scored (ideal distance ≈ 20 nt).
   • UGUA auxiliary motif is searched at −37 ± 5 nt.
   • U-rich auxiliary element is searched at −15 ± 3 nt using a 5 nt sliding window (U fraction ≥ 60%).
5. Confidence assignment:
   • high: hexamer found AND (UGUA or U-rich present)
   • medium: hexamer found only
   • low: no hexamer found

Output files:

• <prefix>_APAread.csv: Per-read poly(A) length and PAS cluster assignment
• <prefix>_APAres.csv: PAS cluster-level motif analysis result
• <prefix>_APAfailed.csv: Reads excluded from APA analysis with reason

Output column description:

<prefix>_APAread.csv:

Column	Description
read_id	Read identifier
cleavage_site	3' cleavage position on the genome (chrom_pos)
PAS_cluster_id	PAS cluster the read belongs to (chrom_start_end_strand)
polyA_length	poly(A) tail length from pt tag

<prefix>_APAres.csv:

Column	Description
PAS_cluster_id	PAS cluster ID
rep_position	Representative cleavage position of the cluster
read_count	Number of reads supporting this cluster
hexamer_motif	Best-matched PAS hexamer (none if no match)
hexamer_distance	Distance from hexamer start to cleavage site
UGUA_support	Whether UGUA auxiliary motif is present
U_rich_support	Whether U-rich auxiliary element is present
confidence	high / medium / low

<prefix>_APAfailed.csv:

Column	Description
read_id	Read identifier
fail_reason	no_polyA / short_oligoA / internal_priming
polyA_length	Detected poly(A) length
cleavage_site	3' end position (chrom_pos)

---

2.7 Intron Retention Rate Identification

Function: Calculate intron retention rate for each read.

Command:

nanornaqtl pheno intron_retention \
  -g <gtf_file> \
  -b <map_bam> \
  -o <output_dir> \
  -p <output_prefix>

Parameter explanation:

Parameter	Description
-g, --gtf	GTF annotation file path
-b, --bam	map.bam file path from prep step (*_map.sorted.bam)
-p, --prefix	Output file prefix name
-o, --output_dir	Output directory path, will be created if not exists

Output files:

• <prefix>_intronRetention_result.csv: Intron retention rate information

Output column description:

Column name	Description
readID	Read identifier
exon_len	Exon length
intron_len	Intron length
len	Total length
IntronRetentionRate	Intron retention rate

---

2.8 isoform Identification

Function: Identify and quantify transcript isoforms.

Tool: Uses IsoQuant tool for isoform analysis.

Command example:

isoquant.py \
  --reference hg38.fa \
  --genedb annotation.gtf \
  --bam map.bam \
  --data_type nanopore \
  -o ./work_dir \
  -t 60 \
  --complete_genedb

Parameter explanation:

Parameter	Description
--reference	Reference genome FASTA file
--genedb	Gene annotation GTF file
--bam	Input BAM file (*_map.sorted.bam)
--data_type	Data type (nanopore/pacbio_ccs)
-o	Output directory
-t	Number of threads
--complete_genedb	Use complete gene database mode

Output files:

IsoQuant will generate multiple output files, the key file is:

• OUT.transcript_model_reads.tsv.gz: Read-to-isoform mapping relationship

This file can be directly used for subsequent QTL analysis.

---

3. qtl - QTL Analysis

All QTL analysis modules require a SNP information file in the following format:

SNP information file format (snp_info.txt):

chrom    pos    A1    A2    AF
chr1    123456    A    G    0.35
chr1    789012    C    T    0.42

Column description:

Column name	Description
chrom	Chromosome name (must match BAM file)
pos	SNP position (1-based)
A1	Alternative allele
A2	Reference allele
AF	Alternative allele frequency

Genome size file format (genome.chrom.sizes):

chr1    249250621
chr2    243199373
chr3    198022430

Format description: Tab-separated, two columns: chromosome name and chromosome length

---

3.1 m6A QTL Analysis

Function: Identify genetic variants associated with m6A modification.

Command:

nanornaqtl qtl m6A \
  -b <map_bam> \
  --snp_info <snp_file> \
  -p <output_prefix> \
  -o <output_dir> \
  --modification_csv <m6A_sites_csv> \
  --read_mod_dict <m6A_reads_pkl> \
  --geno_size <genome_size_file> \
  -t <threads>

Parameter explanation:

Parameter	Description	Default
-b, --bam	map.bam file path from prep step (*_map.sorted.bam)	Required
--snp_info	SNP information file path	Required
-p, --prefix	Output file prefix name	Required
-o, --output_dir	Output directory path, will be created if not exists	Required
--modification_csv	m6A site result file (*_m6A_sites_result.csv)	Required
--read_mod_dict	m6A read dictionary file (*_m6A_reads_final.pkl)	Required
--geno_size	Genome size file	Required
-q, --min_qscore	Minimum base quality score for variant site	10
-c, --min_coverage	Minimum coverage for variant site	8
--mcmc_samples	Number of MCMC samples	2000
-t, --threads	Number of threads	4
--keep_tmp	Keep temporary files	False

Output files:

• <prefix>_m6A_QTLs_result.csv: m6A QTL results

Output column description:

Column name	Description
chrom	Chromosome
SNP	Variant ID
snp_pos_1base	SNP position
mod_pos_1base	Modification site position
A1, A2	Alleles
MAF	Minor allele frequency
A1_mod, A1_unmod	Number of modified/unmodified reads for A1
A2_mod, A2_unmod	Number of modified/unmodified reads for A2
log_OR_beta	Log odds ratio (posterior mean)
log_OR_se	Log odds ratio standard error
bayes_factor	Bayes Factor (H1 vs H0)
fisher_OR	Fisher's exact test odds ratio
fisher_pvalue	Fisher's exact test p-value
fisher_FDR	FDR-corrected Fisher p-value

Statistical method explanation:

1. Bayesian method (primary): Bayes Factor for significance, log_OR for effect size
2. Frequentist method (supplementary): Fisher's exact test with FDR correction

---

3.2 m5C QTL Analysis

Function: Identify genetic variants associated with m5C modification.

Command:

nanornaqtl qtl m5C \
  -b <map_bam> \
  --snp_info <snp_file> \
  -p <output_prefix> \
  -o <output_dir> \
  --modification_csv <m5C_sites_csv> \
  --read_mod_dict <m5C_reads_pkl> \
  --geno_size <genome_size_file> \
  -t <threads>

Parameters: Same as m6A QTL

Output files:

• <prefix>_m5C_QTLs_result.csv: m5C QTL results

Output columns: Same as m6A QTL

---

3.3 pseudouridine (pseU) QTL Analysis

Function: Identify genetic variants associated with pseudouridine modification.

Command:

nanornaqtl qtl pseU \
  -b <map_bam> \
  --snp_info <snp_file> \
  -p <output_prefix> \
  -o <output_dir> \
  --modification_csv <pseU_sites_csv> \
  --read_mod_dict <pseU_reads_pkl> \
  --geno_size <genome_size_file> \
  -t <threads>

Parameters: Same as m6A QTL

Output files:

• <prefix>_pseU_QTLs_result.csv: pseU QTL results

---

3.4 inosine QTL Analysis

Function: Identify genetic variants associated with inosine modification.

Command:

nanornaqtl qtl inosine \
  -b <map_bam> \
  --snp_info <snp_file> \
  -p <output_prefix> \
  -o <output_dir> \
  --modification_csv <inosine_sites_csv> \
  --read_mod_dict <inosine_reads_pkl> \
  --geno_size <genome_size_file> \
  -t <threads>

Parameters: Same as m6A QTL

Output files:

• <prefix>_inosine_QTLs_result.csv: inosine QTL results

---

3.5 2'O-methylation QTL Analysis

Function: Identify genetic variants associated with 2'-O-methylation modification. Supports four subtypes: 2OmeA, 2OmeC, 2OmeU, 2OmeG.

Command:

nanornaqtl qtl 2OmeA \
  -b <map_bam> \
  --snp_info <snp_file> \
  -p <output_prefix> \
  -o <output_dir> \
  --modification_csv <2OmeA_sites_csv> \
  --read_mod_dict <2OmeA_reads_pkl> \
  --geno_size <genome_size_file> \
  -t <threads>

Replace 2OmeA with 2OmeC, 2OmeU, or 2OmeG for other subtypes.

Parameters: Same as m6A QTL

Output files:

• <prefix>_2OmeA_QTLs_result.csv: 2OmeA QTL results (similarly for other subtypes)

Output columns: Same as m6A QTL

---

3.6 APA QTL Analysis

Function: Identify genetic variants affecting alternative polyadenylation site usage.

Command:

nanornaqtl qtl APA \
  -b <map_bam> \
  --snp_info <snp_file> \
  -p <output_prefix> \
  -o <output_dir> \
  --read_overlap_file <apa_result_csv> \
  --geno_size <genome_size_file> \
  -t <threads>

Parameter explanation:

Parameter	Description	Default
-b, --bam	map.bam file path from prep step (*_map.sorted.bam)	Required
--snp_info	SNP information file path	Required
-p, --prefix	Output file prefix name	Required
-o, --output_dir	Output directory path, will be created if not exists	Required
--read_overlap_file	APA result file (*_APA_result.csv)	Required
--geno_size	Genome size file	Required
-q, --min_qscore	Minimum base quality score for variant site	10
-c, --min_coverage	Minimum total coverage at variant site. Each allele must additionally have ≥ min_coverage/2 reads	8
--n_samples	Number of Dirichlet posterior samples for TVD effect size estimation	2000
-t, --threads	Number of threads	4
--keep_tmp	Keep temporary files	False

Output files:

• <prefix>_APA_QTLs_result.csv: APA QTL results

Output column description:

Column name	Description
chrom	Chromosome
strand	Strand (+/-)
snp_pos_1base	SNP position (1-based)
SNP	Variant ID
A1, A2	Alleles
MAF	Minor allele frequency
all_types	List of APA types observed at this SNP
A1_counts	Read counts per APA type carrying A1 allele
A2_counts	Read counts per APA type carrying A2 allele
bayes_factor	Bayes Factor (H1: A1 and A2 have different APA distributions vs H0: same)
TVD_beta	Posterior mean of Total Variation Distance (effect size, range 0-1)
TVD_se	Posterior standard deviation of TVD
dominant_shift	Main APA type usage change direction and magnitude
g_pvalue	G-test p-value (Monte Carlo permutation when expected counts < 5 or zero counts present)
g_FDR	Benjamini-Hochberg FDR-adjusted g_pvalue

Statistical method explanation:

1. Bayesian method (recommended): Bayes Factor with uniform Dirichlet prior. BF > 10 indicates strong evidence, BF > 100 indicates decisive evidence (Kass & Raftery 1995).
2. Frequentist method (reference): G-test (likelihood-ratio chi-square). Automatically falls back to Monte Carlo permutation when asymptotic approximation is unreliable (any expected cell count < 5 or any observed count = 0). FDR correction by Benjamini-Hochberg.

Effect size indicator explanation:

• TVD_beta (Total Variation Distance, posterior mean): Measures the degree of difference in APA usage distribution between A1 and A2 alleles, computed as the posterior mean of 0.5 × Σ|π₁ₖ − π₂ₖ| from Dirichlet posterior sampling.

  • 0: Two distributions are completely identical
  • 1: Two distributions are completely different

• TVD_se: Posterior standard deviation of TVD, indicating the uncertainty of the effect size estimate.

• dominant_shift: Describes the change in usage proportion of each APA type between A1 and A2.

  • Format example: type1:A1↓(-0.56) indicates reads carrying A1 allele use this APA type less (absolute proportion decrease of 0.56)

---

3.7 isoform QTL Analysis

Function: Identify genetic variants affecting transcript isoform usage.

Command:

nanornaqtl qtl isoform \
  -b <map_bam> \
  --snp_info <snp_file> \
  -p <output_prefix> \
  -o <output_dir> \
  --read_overlap_file <isoquant_output> \
  --geno_size <genome_size_file> \
  -t <threads>

Parameter explanation:

Parameter	Description
--read_overlap_file	IsoQuant output file (OUT.transcript_model_reads.tsv.gz)
Other parameters	Same as APA QTL

Output files:

• <prefix>_isoform_QTLs_result.csv: isoform QTL results

Output columns: Same as APA QTL. The all_types, A1_counts, A2_counts columns describe transcript IDs and their read counts per allele instead of APA types.

---

3.8 polyA Tail Length QTL Analysis

Function: Identify genetic variants affecting polyA tail length using JZS Bayesian two-sample t-test (Rouder et al., 2009) and Mann-Whitney U test.

Command:

nanornaqtl qtl polyA_tail \
  -b <map_bam> \
  --snp_info <snp_file> \
  -p <output_prefix> \
  -o <output_dir> \
  --polya_csv <polyA_csv> \
  --geno_size <genome_size_file> \
  -t <threads>

Parameter explanation:

Parameter	Description	Default
-b, --bam	map.bam file path from prep step (*_map.sorted.bam)	Required
--snp_info	SNP information file path	Required
-p, --prefix	Output file prefix name	Required
--polya_csv	polyA tail length result CSV file path (*_APAread.csv)	Required
--geno_size	Genome size file	Required
-q, --min_qscore	Minimum base quality score for variant site	10
-c, --min_coverage	Minimum coverage for variant site	8
-t, --threads	Number of threads	4
--keep_tmp	Keep temporary files	False

Output files:

• <prefix>_polyA_tail_length_QTLs_result.csv: polyA tail length QTL results

Output column description:

Column name	Description
chrom	Chromosome
strand	Strand
SNP	Variant ID
snp_pos_1base	SNP position (1-based)
A1, A2	Alleles
MAF	Minor allele frequency
A1_len	polyA tail length list for A1 allele
A2_len	polyA tail length list for A2 allele
beta	Effect size (log ratio of mean lengths on original scale)
SE	Standard error of beta
cohens_d	Cohen's d effect size (on log-transformed scale)
bayes_factor	JZS Bayes Factor (BF>3: moderate evidence, BF>10: strong evidence)
p_mw	Mann-Whitney U test p-value
p_mw_fdr	FDR-corrected Mann-Whitney U p-value (Benjamini-Hochberg)

Statistical methods:

• Bayesian: JZS (Jeffreys-Zellner-Siow) Bayes Factor based on two-sample t-test with log-transformed polyA tail lengths. Provides evidence strength for association without dependence on sample size.
• Frequentist: Mann-Whitney U test on original polyA tail lengths with FDR correction. Non-parametric, robust to non-normality and outliers.

---

3.9 Intron Retention Rate QTL Analysis

Function: Identify genetic variants affecting intron retention rate.

Command:

nanornaqtl qtl intron_retention \
  -b <map_bam> \
  --snp_info <snp_file> \
  -p <output_prefix> \
  -o <output_dir> \
  --ir_csv <ir_csv> \
  --geno_size <genome_size_file> \
  -t <threads>

Parameter explanation:

Parameter	Description	Default
--ir_csv	Intron retention rate result file (*_intronRetention_result.csv)	Required
Other parameters	Same as polyA tail length QTL

Output files:

• <prefix>_intron_retention_QTLs_result.csv: Intron retention rate QTL results

Output column description:

Column name	Description
chrom	Chromosome
strand	Strand
snp_pos_1base	SNP position (1-based)
SNP	Variant ID
A1, A2	Alleles
MAF	Minor allele frequency
beta	Posterior mean of logit-difference in overall IR mean (A1 − A2)
SE	Posterior standard deviation of the logit-difference
bayes_factor	Total Bayes Factor (BF_zero × BF_intensity)
BF_zero	Bayes Factor for the zero-inflation layer (difference in non-zero proportion between A1 and A2)
BF_intensity	Bayes Factor for the intensity layer (difference in IR distribution conditional on IR > 0)
MW_pvalue	Mann-Whitney U test p-value (two-sided)

Note: SNPs are filtered by --min_coverage (default=8, total reads across A1 and A2). The model is a two-layer zero-inflated Beta: BF_zero tests for differential probability of triggering intron retention (IR > 0); BF_intensity tests for differential IR magnitude conditional on IR > 0.

---

Output File Description

File Naming Rules

• prep module: <prefix>_map[0/16].sorted.bam, <prefix>_map.fastq
• pheno module: <prefix>_<phenotype>_result.csv
• qtl module: <prefix>_<phenotype>_QTLs_result.csv

Common Output Format

All CSV output files are standard tab-separated or comma-separated files that can be opened and analyzed with Excel, R, Python, and other tools.

Key Statistical Indicator Interpretation 📈

Bayesian Methods

• BF (Bayes Factor):
  • BF > 3: Significant evidence supporting association
  • BF > 10: Strong evidence
  • BF > 30: Very strong evidence

Frequentist Methods

• p-value: Traditional significance test p-value, typically p < 0.05 is considered significant

Comparison Between BF and p-value in This Study

We compared SNPs called significant by BF > 3 versus those called by p < 0.05. The two sets are largely overlapping, but each method has a small fraction of method-specific hits, and these hits are systematically distributed along the coverage axis:

• p-value-only significant hits are enriched at very high coverage. With large sample sizes, the Mann-Whitney U test becomes overpowered: even biologically negligible differences in IR distribution can yield p < 0.05, inflating the false-positive rate.

• BF-only significant hits are enriched at low coverage. This reflects a fundamental difference in how the two methods quantify evidence. Frequentist tests rely on the sampling distribution under H0 and lose power rapidly when n is small (the discrete p-value is bounded away from zero by the limited number of rank permutations). In contrast, the Bayes Factor compares marginal likelihoods of H1 versus H0 directly, so when a small number of reads already show clear separation between A1 and A2 (e.g., A1 reads cluster near IR = 0 while A2 reads cluster at higher IR), the BF can still register strong evidence. The weakly informative priors used in the model (Beta(1,1) for the zero-inflation layer, Gamma(2, 0.1) for the intensity layer) also stabilize BF estimation at low coverage.

Recommendation: BF is preferred as the primary metric in this study, because it remains calibrated across the coverage spectrum—neither over-rejecting at high coverage nor under-rejecting at low coverage.

---

Citation

If you use nanornaqtl in your research, please cite our GitHub repository:

nanornaqtl: A comprehensive toolkit for molecular phenotyping and QTL analysis using Nanopore direct RNA sequencing
https://github.com/xinranxu0930/nanornaqtl

---

Contact 📧

• Issue reporting: GitHub Issues
• Email: [e-mail]
• GitHub: https://github.com/xinranxu0930/nanornaqtl

---

License

This project is licensed under the MIT License. See LICENSE file for details.

---

Changelog

v1.0.0 (2026-01-18)

• Initial release
• Support for 7 types of molecular phenotype identification
• Support for 8 types of QTL analysis
• Implementation of Bayesian statistical methods
• Support for multi-threading parallel processing

v1.0.3 (2026-01-20)

• Fixed coverage calculation
• Added logging system
• Optimized base mode support issues
• Added version query

v1.0.4 (2026-01-27)

• Added new de novo APA identification feature
• Fixed parameter description errors in README

v2.0.0 (2026-05-07)

• Added minimum read length filtering for mapped reads
• Added read-level sequencing quality (qscore) filtering based on Dorado qs tag
• Added RNA type classification (mRNA / ncRNA / unknown) based on user-provided gene annotation, with biotype tag (RT:Z:<biotype>) attached to ncRNA reads
• Added [BETA] 5' end-based degradation detection via peak calling, outputting per-read integrity labels (intact / ambiguous / truncated / undetermined) and a degraded read list
• Added identification of four 2'-O-methylation modifications (2OmeA, 2OmeC, 2OmeU, 2OmeG)
• Updated motif rules for m5C, pseU, and inosine based on the latest literature
• Added poly(A) length-based filtering (no_polyA / short_oligoA / pass)
• Added internal priming filtering and PAS hexamer / UGUA / U-rich motif characterization for APA site identification
• Merged polyA tail length and APA site identification into a single sub-command