Metadata-Version: 2.4
Name: nanornaqtl
Version: 2.0.0
Summary: A comprehensive toolkit for multi-molecular phenotyping and QTL analysis based on Nanopore direct RNA sequencing data
Author-email: Xinran Xu <xinranxu0930@gmail.com>
License: MIT
Requires-Python: >=3.9
Description-Content-Type: text/markdown
License-File: LICENSE
Dynamic: license-file

# nanornaqtl

[![PyPI version](https://badge.fury.io/py/nanornaqtl.svg)](https://badge.fury.io/py/nanornaqtl)
[![Python versions](https://img.shields.io/pypi/pyversions/nanornaqtl.svg)](https://pypi.org/project/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 📦

```bash
pip install nanornaqtl
```

### Install from source

```bash
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](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):

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

### Step 2: Basecalling and Alignment 🔧

Use [Dorado](https://github.com/nanoporetech/dorado) for basecalling and alignment:

#### Basic Alignment (Required)

```bash
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**:

```bash
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:
```bash
# 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**:

```bash
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**:

```bash
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**:

```bash
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](https://github.com/olarerin/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**:

```bash
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**:

```bash
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**:

```bash
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**:

```bash
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**:

```bash
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**:

```bash
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](https://github.com/ablab/IsoQuant) tool for isoform analysis.

**Command example**:

```bash
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**:

```bash
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**:

```bash
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**:

```bash
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**:

```bash
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**:

```bash
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**:

```bash
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**:

```bash
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**:

```bash
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**:

```bash
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](https://github.com/xinranxu0930/nanornaqtl/issues)
- **Email**: xinranxu0930@gmail.com
- **GitHub**: [https://github.com/xinranxu0930/nanornaqtl](https://github.com/xinranxu0930/nanornaqtl)

---

## License

This project is licensed under the MIT License. See [LICENSE](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
