Merge pull request #238 from ClareRobin/clare_dev_pr

Add genome-alignment and RNA004 support (xPore v2.2)

Author: ClareRobin

.github/workflows/test.yml

@@ -0,0 +1,26 @@
+name: tests
+
+on:
+  push:
+  pull_request:
+
+jobs:
+  test:
+    # ubuntu-22.04 (not ubuntu-latest) so that Python 3.8 is reliably available via setup-python.
+    runs-on: ubuntu-22.04
+    strategy:
+      fail-fast: false
+      matrix:
+        python-version: ["3.8", "3.9", "3.10"]
+    steps:
+      - uses: actions/checkout@v4
+      - name: Set up Python ${{ matrix.python-version }}
+        uses: actions/setup-python@v5
+        with:
+          python-version: ${{ matrix.python-version }}
+      - name: Install package with test dependencies
+        run: |
+          python -m pip install --upgrade pip
+          pip install .[test]
+      - name: Run tests
+        run: pytest -v

MANIFEST.in

@@ -1 +1,2 @@
-include xpore/diffmod/model_kmer.csv
+include xpore/diffmod/RNA002_5mer_model.csv
+include xpore/diffmod/RNA004_5mer_model.csv

README.md

@@ -33,7 +33,7 @@ xPore is described in details in [Pratanwanich et al. *Nat Biotechnol* (2021)](h
 
 ### Release History
 
-The current release is xPore v2.1. 
+The current release is xPore v2.2, which adds support for genome-aligned eventalign output (via the new `--kmer_source model_kmer` flag) and RNA004 data (now the default `xpore-diffmod` prior). 
 
 Please refer to the github release history for previous releases: https://github.com/GoekeLab/xpore/releases
 

docs/source/cmd.rst

@@ -17,12 +17,13 @@ Argument name                       Required    Default value         Descriptio
 =================================   ==========  ===================  ============================================================================================================
 --eventalign=FILE                   Yes         NA                    Eventalign filepath, the output from nanopolish.         
 --out_dir=DIR                       Yes         NA                    Output directory.
---gtf_path_or_url                   No          NA                    GTF file path or url used for mapping transcriptomic to genomic coordinates.
---transcript_fasta_paths_or_urls    No          NA                    Transcript FASTA paths or urls used for mapping transcriptomic to genomic coordinates.
+--gtf_or_gff=FILE                   No          NA                    GTF or GFF file path used for mapping transcriptomic to genomic coordinates.
+--transcript_fasta=FILE             No          NA                    Transcript FASTA path used for mapping transcriptomic to genomic coordinates.
 --skip_eventalign_indexing          No          False                 To skip indexing the eventalign nanopolish output.
 --genome                            No          False                 To run on Genomic coordinates. Without this argument, the program will run on transcriptomic coordinates.
+--kmer_source=STR                   No          reference_kmer        Which kmer column to use from the eventalign file: ``reference_kmer`` (default, for transcriptome alignments) or ``model_kmer`` (for genome alignments, which contain reverse-strand reads).
 --n_processes=NUM                   No          1                     Number of processes to run.
---readcount_max=NUM                 No          1000                  Maximum read counts per gene.
+--readcount_max=NUM                 No          1000                  Maximum read counts per site. Use None for no limit.
 --readcount_min=NUM                 No          1                     Minimum read counts per gene.
 --resume                            No          False                 With this argument, the program will resume from the previous run.
 =================================   ==========  ===================  ============================================================================================================

docs/source/conf.py

@@ -26,9 +26,9 @@
 author = 'Ploy N. Pratanwanich'
 
 # The short X.Y version
-version = '2.0'
+version = '2.2'
 # The full version, including alpha/beta/rc tags
-release = '2.0'
+release = '2.2'
 
 
 # -- General configuration ---------------------------------------------------

docs/source/configuration.rst

@@ -21,7 +21,14 @@ Only the ``data`` and ``out`` sections are required, other sections are optional
         ...
 
     out: <DIR_PATH_FOR_OUTPUTS>
-    
+
+    # Optional: path to a CSV of unmodified-RNA kmer signal priors with columns
+    # model_kmer,model_mean,model_stdv. If omitted, xpore-diffmod defaults to the
+    # bundled RNA004 5mer model (xpore/diffmod/RNA004_5mer_model.csv).
+    # For RNA002 data, set this to the bundled RNA002 model, e.g.:
+    #     prior: /path/to/xpore/diffmod/RNA002_5mer_model.csv
+    prior: <PATH_TO_PRIOR_CSV>
+
     criteria:
         readcount_min: <15>
         readcount_max: <1000>

docs/source/outputtable.rst

@@ -18,6 +18,6 @@ mu_mod                                      inferred mean of the modified RNAs d
 sigma2_unmod                                inferred sigma^2 of the unmodified RNAs distribution
 sigma2_mod                                  inferred sigma^2 of the modified RNAs distribution
 conf_mu_unmod                               confidence level of mu_unmod compared to the unmodified reference signal
-conf_mu_mod                                 confidence level of mu_unmod compared to the unmodified reference signal
+conf_mu_mod                                 confidence level of mu_mod compared to the unmodified reference signal
 mod_assignment                              lower if mu_mod < mu_unmod and higher if mu_mod > mu_unmod
 ==========================================  ========================================================================================================================================

docs/source/preparation.rst

@@ -3,24 +3,38 @@
 Data preparation from raw reads
 ===================================
 
-1. After obtaining fast5 files, the first step is to basecall them. Below is an example script to run Guppy basecaller. You can find more detail about basecalling at `Oxford nanopore Technologies <https://nanoporetech.com>`_::
-
-    guppy_basecaller -i </PATH/TO/FAST5> -s </PATH/TO/FASTQ> --flowcell <FLOWCELL_ID> --kit <KIT_ID> --device auto -q 0 -r   
-
-2. Align to transcriptome::
-   
-    minimap2 -ax map-ont -uf -t 3 --secondary=no <MMI> <PATH/TO/FASTQ.GZ> > <PATH/TO/SAM> 2>> <PATH/TO/SAM_LOG>
-    samtools view -Sb <PATH/TO/SAM> | samtools sort -o <PATH/TO/BAM> - &>> <PATH/TO/BAM_LOG>
-    samtools index <PATH/TO/BAM> &>> <PATH/TO/BAM_INDEX_LOG>
-
-3. Resquiggle using `nanopolish eventalign <https://nanopolish.readthedocs.io/en/latest/quickstart_eventalign.html>`_::
-   
-    nanopolish index -d <PATH/TO/FAST5_DIR> <PATH/TO/FASTQ_FILE>
-    nanopolish eventalign --reads <PATH/TO/FASTQ_FILE> \
-    --bam <PATH/TO/BAM_FILE> \
-    --genome <PATH/TO/FASTA_FILE \
-    --signal-index \
-    --scale-events \
-    --summary <PATH/TO/summary.txt> \
-    --threads 32 > <PATH/TO/eventalign.txt>
+1. After obtaining the raw signal files (POD5, or FAST5 for older runs), the first step is to basecall them. Below is an example using `Dorado <https://github.com/nanoporetech/dorado>`_, Oxford Nanopore's current basecaller (it replaces the older Guppy/Albacore basecallers and supports both RNA002 and RNA004 chemistries — select the model that matches your chemistry). You can find more detail about basecalling at `Oxford Nanopore Technologies <https://nanoporetech.com>`_::
 
+    dorado basecaller <MODEL> </PATH/TO/POD5_DIR> --emit-fastq > <PATH/TO/FASTQ>
+
+   For RNA004 data use an ``rna004`` model (e.g. ``rna004_130bps_sup@v5.1.0``); for RNA002 data use an ``rna002`` model. See the `Dorado documentation <https://github.com/nanoporetech/dorado>`_ for the available models and options.
+
+2. Align the basecalled reads with `minimap2 <https://github.com/lh3/minimap2>`_. xPore supports both **transcriptome** and **genome** alignments — choose one depending on which coordinate system you want in the output.
+
+   **Transcriptome alignment** (align to a transcriptome reference). xPore reports transcriptomic coordinates, or genomic coordinates if you also pass ``--genome`` together with ``--gtf_or_gff`` and ``--transcript_fasta`` to ``xpore dataprep``::
+
+       minimap2 -ax map-ont -uf -t 3 --secondary=no <TRANSCRIPTOME.MMI> <PATH/TO/FASTQ.GZ> > <PATH/TO/SAM> 2>> <PATH/TO/SAM_LOG>
+       samtools view -Sb <PATH/TO/SAM> | samtools sort -o <PATH/TO/BAM> - &>> <PATH/TO/BAM_LOG>
+       samtools index <PATH/TO/BAM> &>> <PATH/TO/BAM_INDEX_LOG>
+
+   **Genome alignment** (align directly to a genome reference; use spliced alignment so that reads spanning introns map correctly). Genome alignments contain reverse-strand reads, so run ``xpore dataprep`` with ``--kmer_source model_kmer`` for these (see :ref:`Command line arguments <cmd>`)::
+
+       minimap2 -ax splice -uf -k14 -t 3 --secondary=no <GENOME.MMI> <PATH/TO/FASTQ.GZ> > <PATH/TO/SAM> 2>> <PATH/TO/SAM_LOG>
+       samtools view -Sb <PATH/TO/SAM> | samtools sort -o <PATH/TO/BAM> - &>> <PATH/TO/BAM_LOG>
+       samtools index <PATH/TO/BAM> &>> <PATH/TO/BAM_INDEX_LOG>
+
+3. Resquiggle (align the raw signal to the reference) to produce the eventalign file. We recommend `f5c <https://github.com/hasindu2008/f5c>`_, an optimised, CPU/GPU-accelerated re-implementation of ``nanopolish eventalign`` that produces equivalent output much faster on large datasets::
+
+       # index the raw signal: use -d <FAST5_DIR> for FAST5, or --slow5 <FILE.blow5> for SLOW5/BLOW5
+       f5c index -d <PATH/TO/FAST5_DIR> <PATH/TO/FASTQ_FILE>
+       f5c eventalign --reads <PATH/TO/FASTQ_FILE> \
+       --bam <PATH/TO/BAM_FILE> \
+       --genome <PATH/TO/REFERENCE_FASTA> \
+       --rna \
+       --signal-index \
+       --scale-events \
+       --threads 32 > <PATH/TO/eventalign.txt>
+
+   For **RNA004** data, add ``--kmer-model <PATH/TO/5-mer-model>``: recent versions of f5c auto-select the 9-mer model for RNA004, which xPore cannot use — xPore requires the 5-mer model (see `xPore issue #215 <https://github.com/GoekeLab/xpore/issues/215>`_).
+
+   ``nanopolish eventalign`` can be used instead with the same arguments. Note that the ``--genome`` argument here refers to the **alignment reference** (the transcriptome or genome FASTA used in step 2), not xPore's ``--genome`` flag.

docs/source/quickstart.rst

@@ -3,6 +3,44 @@
 Quickstart - Detection of differential RNA modifications
 =========================================================
 
+.. note::
+   **Updates in xPore v2.2:** xPore is now compatible with genome alignments and RNA004 data — see the table below and the :ref:`Data preparation from raw reads <preparation>` section for more information.
+
+xPore is now compatible with genome alignment. See below for the minimal commands to run xPore on transcriptome- or genome-aligned data:
+
+.. list-table::
+   :header-rows: 1
+   :widths: 22 48 30
+
+   * -
+     - **Dataprep**
+     - **Diffmod**
+   * - Transcriptome alignment (output transcriptome coordinates)
+     - | ``xpore dataprep``
+       | ``--eventalign <eventalign.txt>``
+       | ``--out_dir <out_dir>``
+     - | ``xpore diffmod``
+       | ``--config <config.yml>``
+   * - Transcriptome alignment (output genome coordinates)
+     - | ``xpore dataprep``
+       | ``--eventalign <eventalign.txt>``
+       | ``--out_dir <out_dir>``
+       | ``--genome``
+       | ``--transcript_fasta <transcript.fa>``
+       | ``--gtf_or_gff <annotation.gtf>``
+     - | ``xpore diffmod``
+       | ``--config <config.yml>``
+   * - Genome alignment
+     - | ``xpore dataprep``
+       | ``--eventalign <eventalign.txt>``
+       | ``--out_dir <out_dir>``
+       | ``--kmer_source model_kmer``
+     - | ``xpore diffmod``
+       | ``--config <config.yml>``
+
+Running Example Demo Data
+-------------------------
+
 Download and extract the demo dataset from our `zenodo <https://zenodo.org/record/5162402/files/demo.tar.gz>`_::
 
     wget https://zenodo.org/record/5162402/files/demo.tar.gz
@@ -61,8 +99,12 @@ Below is how it looks like::
 
     out: ./out # output dir
 
+    # The demo data is RNA002. Since v2.2 the default prior is the RNA004 model,
+    # so point xpore-diffmod at the bundled RNA002 model to reproduce the demo:
+    prior: /path/to/xpore/diffmod/RNA002_5mer_model.csv
+
 
-See the :ref:`Configuration file page <configuration>` for more details.
+See the :ref:`Configuration file page <configuration>` for more details. Note that since xPore v2.2 the default unmodified-signal prior is the RNA004 model; for RNA002 data (like this demo) set ``prior:`` to the bundled ``RNA002_5mer_model.csv`` as shown above.
 
 3. Now that we have the data and the configuration file ready for modelling differential modifications using ``xpore-diffmod``. 
 

setup.py

@@ -14,7 +14,7 @@
     name=__pkg_name__,
     license="MIT",
     description='xpore is a python package for Nanopore data analysis of differential RNA modifications.',
-    version='v2.1',
+    version='v2.2',
     long_description=README,
     long_description_content_type='text/markdown',
     url='https://github.com/GoekeLab/xpore',
@@ -30,6 +30,7 @@
             'ujson>=4.0.1'
             ],
     python_requires=">=3.8",
+    extras_require={'test': ['pytest']},
     entry_points={'console_scripts': ["xpore={}.scripts.xpore:main".format(__pkg_name__)]},
     classifiers=[
         # Trove classifiers