diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml new file mode 100644 index 0000000..c42bd9e --- /dev/null +++ b/.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 diff --git a/MANIFEST.in b/MANIFEST.in index d0f20f9..4d7f0f8 100644 --- a/MANIFEST.in +++ b/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 diff --git a/README.md b/README.md index 9697ec0..7e4cde8 100644 --- a/README.md +++ b/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 diff --git a/docs/source/cmd.rst b/docs/source/cmd.rst index 75f9c2f..0c8e37a 100644 --- a/docs/source/cmd.rst +++ b/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. ================================= ========== =================== ============================================================================================================ diff --git a/docs/source/conf.py b/docs/source/conf.py index c9a245c..d3fa679 100644 --- a/docs/source/conf.py +++ b/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 --------------------------------------------------- diff --git a/docs/source/configuration.rst b/docs/source/configuration.rst index 3c8befa..0f89d1a 100644 --- a/docs/source/configuration.rst +++ b/docs/source/configuration.rst @@ -21,7 +21,14 @@ Only the ``data`` and ``out`` sections are required, other sections are optional ... out: - + + # 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: + criteria: readcount_min: <15> readcount_max: <1000> diff --git a/docs/source/outputtable.rst b/docs/source/outputtable.rst index 4b876a4..9aa1665 100644 --- a/docs/source/outputtable.rst +++ b/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 ========================================== ======================================================================================================================================== diff --git a/docs/source/preparation.rst b/docs/source/preparation.rst index ac44f7b..ab4456c 100644 --- a/docs/source/preparation.rst +++ b/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 `_:: - - guppy_basecaller -i -s --flowcell --kit --device auto -q 0 -r - -2. Align to transcriptome:: - - minimap2 -ax map-ont -uf -t 3 --secondary=no > 2>> - samtools view -Sb | samtools sort -o - &>> - samtools index &>> - -3. Resquiggle using `nanopolish eventalign `_:: - - nanopolish index -d - nanopolish eventalign --reads \ - --bam \ - --genome \ - --threads 32 > +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 `_, 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 `_:: + dorado basecaller --emit-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 `_ for the available models and options. + +2. Align the basecalled reads with `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 > 2>> + samtools view -Sb | samtools sort -o - &>> + samtools index &>> + + **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 `):: + + minimap2 -ax splice -uf -k14 -t 3 --secondary=no > 2>> + samtools view -Sb | samtools sort -o - &>> + samtools index &>> + +3. Resquiggle (align the raw signal to the reference) to produce the eventalign file. We recommend `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 for FAST5, or --slow5 for SLOW5/BLOW5 + f5c index -d + f5c eventalign --reads \ + --bam \ + --genome \ + --rna \ + --signal-index \ + --scale-events \ + --threads 32 > + + For **RNA004** data, add ``--kmer-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 `_). + + ``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. diff --git a/docs/source/quickstart.rst b/docs/source/quickstart.rst index 0ac8a20..4d5ca44 100644 --- a/docs/source/quickstart.rst +++ b/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 ` 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 `` + | ``--out_dir `` + - | ``xpore diffmod`` + | ``--config `` + * - Transcriptome alignment (output genome coordinates) + - | ``xpore dataprep`` + | ``--eventalign `` + | ``--out_dir `` + | ``--genome`` + | ``--transcript_fasta `` + | ``--gtf_or_gff `` + - | ``xpore diffmod`` + | ``--config `` + * - Genome alignment + - | ``xpore dataprep`` + | ``--eventalign `` + | ``--out_dir `` + | ``--kmer_source model_kmer`` + - | ``xpore diffmod`` + | ``--config `` + +Running Example Demo Data +------------------------- + Download and extract the demo dataset from our `zenodo `_:: 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 ` for more details. +See the :ref:`Configuration file page ` 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``. diff --git a/setup.py b/setup.py index 1f47df9..294f341 100644 --- a/setup.py +++ b/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 diff --git a/tests/conftest.py b/tests/conftest.py new file mode 100644 index 0000000..234c340 --- /dev/null +++ b/tests/conftest.py @@ -0,0 +1,108 @@ +"""Shared helpers and fixtures for the xpore test suite. + +The functions here build the small synthetic inputs the tests need: + +* ``make_eventalign_line`` / ``make_eventalign_str`` -- eventalign text, the input to + ``dataprep.combine``. +* ``make_events_array`` -- a single read's event array (what ``combine`` returns and what + ``preprocess_tx`` consumes), built with fixed string widths so several reads can be + ``np.concatenate``-d together. +* ``read_json_output`` -- parse the ``data.json`` that ``preprocess_tx`` writes. + +The ``locks`` and ``out_paths`` fixtures provide the two arguments ``preprocess_tx`` writes +through, so the tests never touch multiprocessing. +""" +import json +import threading + +import numpy as np +import pytest + + +# Column order of a nanopolish/f5c eventalign.txt, exactly as combine() parses it. +EVENTALIGN_COLUMNS = [ + 'contig', 'position', 'reference_kmer', 'read_index', 'strand', + 'event_index', 'event_level_mean', 'event_stdv', 'event_length', + 'model_kmer', 'model_mean', 'model_stdv', 'standardized_level', + 'start_idx', 'end_idx', +] + +# Sensible defaults for one eventalign row; override any field per call. With start_idx=0, +# end_idx=5 the per-row length is 5, so a lone row's norm_mean == its event_level_mean. +_EVENTALIGN_DEFAULTS = { + 'contig': 'ENST1', + 'position': 100, + 'reference_kmer': 'GGACT', + 'read_index': 0, + 'strand': 't', + 'event_index': 1, + 'event_level_mean': 120.0, + 'event_stdv': 2.0, + 'event_length': 0.005, + 'model_kmer': 'GGACT', + 'model_mean': 120.0, + 'model_stdv': 1.5, + 'standardized_level': 0.5, + 'start_idx': 0, + 'end_idx': 5, +} + + + +#helper functions for the combine() tests (in test_combine.py) +def make_eventalign_line(**overrides): + """Return one tab-separated eventalign row (no trailing newline).""" + row = dict(_EVENTALIGN_DEFAULTS, **overrides) + return '\t'.join(str(row[col]) for col in EVENTALIGN_COLUMNS) + +def make_eventalign_str(rows): + """Join a list of row-override dicts into one eventalign string (what combine expects).""" + return '\n'.join(make_eventalign_line(**row) for row in rows) + + + +#helper functions for the preprocess_tx() tests (in test_preprocess_tx.py) +def _events_dtype(kmer_col): + # Matches combine()'s output fields. Fixed string widths let arrays for different reads + # be np.concatenate-d inside preprocess_tx. + return np.dtype([ + ('transcript_id', '1 position (size > 1) +# so the caller's inclusion rule keeps it. +READS = [ + {"read_index": 0, "events": "\n".join([ + _line(100, "GGACT", 0, 1, 120.0, 1000, 1005), + _line(101, "GACTA", 0, 2, 95.0, 1005, 1010), + _line(102, "ACTAG", 0, 3, 104.0, 1010, 1015), + ])}, + {"read_index": 1, "events": "\n".join([ + _line(100, "GGACT", 1, 1, 121.0, 2000, 2005), + _line(101, "GACTA", 1, 2, 96.0, 2005, 2010), + _line(102, "ACTAG", 1, 3, 105.0, 2010, 2015), + ])}, + {"read_index": 2, "events": "\n".join([ + _line(101, "GACTA", 2, 1, 94.0, 3000, 3005), + _line(102, "ACTAG", 2, 2, 106.0, 3005, 3010), + _line(103, "CTAGT", 2, 3, 110.0, 3010, 3015), + ])}, +] + + +def _run_combine(events): + result = combine(events) # v2.1 returns array; v2.2 returns (array, kmer_col) + return result[0] if isinstance(result, tuple) else result + + +def _call_preprocess_tx(tx_id, data_dict, kmer_col, out_paths, locks): + # v2.2 added kmer_col + readcount_max params; v2.1 has neither. readcount_max=None on v2.2 + # matches v2.1 (no per-site cap), so both reshape the same reads. + if "kmer_col" in inspect.signature(preprocess_tx).parameters: + preprocess_tx(tx_id, data_dict, kmer_col, None, out_paths, locks) + else: + preprocess_tx(tx_id, data_dict, out_paths, locks) + + +def main(): + data_dict, kmer_col = {}, "reference_kmer" + for read in READS: + np_events = _run_combine(read["events"]) + kmer_col = np_events.dtype.names[2] + if np_events.size > 1: # mirror the caller's inclusion rule + data_dict[read["read_index"]] = np_events + + with tempfile.TemporaryDirectory() as tmp: + out_paths = {name: os.path.join(tmp, "data.%s" % name) + for name in ("json", "index", "readcount", "log")} + locks = {name: threading.Lock() for name in out_paths} + _call_preprocess_tx(TX_ID, data_dict, kmer_col, out_paths, locks) + expected = {} + with open(out_paths["json"]) as f: + for line in f: + if line.strip(): + expected.update(json.loads(line)) + + fixture = { + "_comment": "xPore v2.1 data.json reference captured from upstream/master; regenerate via _regenerate_data_json_v2_1.py.", + "tx_id": TX_ID, + "reads": READS, + "expected": expected, + } + out_path = os.path.join(os.path.dirname(__file__), "data_json_transcriptome_v2_1.json") + with open(out_path, "w") as f: + json.dump(fixture, f, indent=2) + f.write("\n") + + print("imported xpore from:", xpore.__file__) + print("wrote:", out_path) + print(json.dumps(expected, indent=2)) + + +if __name__ == "__main__": + main() diff --git a/tests/fixtures/combine_transcriptome_v2_1.json b/tests/fixtures/combine_transcriptome_v2_1.json new file mode 100644 index 0000000..8818bb7 --- /dev/null +++ b/tests/fixtures/combine_transcriptome_v2_1.json @@ -0,0 +1,25 @@ +{ + "_comment": "xPore v2.1 reference output captured from upstream/master; regenerate via _regenerate_combine_v2_1.py.", + "input_eventalign": "ENST1\t100\tGGACT\t0\tt\t1\t120.5\t2.0\t0.01\tGGACT\t120.0\t1.5\t0.5\t1000\t1005\nENST1\t100\tGGACT\t0\tt\t2\t118.0\t2.1\t0.01\tGGACT\t120.0\t1.5\t0.4\t1005\t1008\nENST1\t101\tGACTA\t0\tt\t3\t95.0\t1.8\t0.01\tGACTA\t95.5\t1.4\t0.3\t1008\t1013\nENST1\t102\tACTAG\t0\tt\t4\t104.5\t1.9\t0.01\tACTAG\t104.0\t1.6\t0.2\t1013\t1018", + "kmer_col": "reference_kmer", + "expected": [ + [ + "ENST1", + 102, + "GGACT", + 119.6 + ], + [ + "ENST1", + 103, + "GACTA", + 95.0 + ], + [ + "ENST1", + 104, + "ACTAG", + 104.5 + ] + ] +} diff --git a/tests/fixtures/data_json_transcriptome_v2_1.json b/tests/fixtures/data_json_transcriptome_v2_1.json new file mode 100644 index 0000000..a137fde --- /dev/null +++ b/tests/fixtures/data_json_transcriptome_v2_1.json @@ -0,0 +1,47 @@ +{ + "_comment": "xPore v2.1 data.json reference captured from upstream/master; regenerate via _regenerate_data_json_v2_1.py.", + "tx_id": "ENST1", + "reads": [ + { + "read_index": 0, + "events": "ENST1\t100\tGGACT\t0\tt\t1\t120.0\t2.0\t0.01\tGGACT\t120.0\t1.5\t0.5\t1000\t1005\nENST1\t101\tGACTA\t0\tt\t2\t95.0\t2.0\t0.01\tGACTA\t95.0\t1.5\t0.5\t1005\t1010\nENST1\t102\tACTAG\t0\tt\t3\t104.0\t2.0\t0.01\tACTAG\t104.0\t1.5\t0.5\t1010\t1015" + }, + { + "read_index": 1, + "events": "ENST1\t100\tGGACT\t1\tt\t1\t121.0\t2.0\t0.01\tGGACT\t121.0\t1.5\t0.5\t2000\t2005\nENST1\t101\tGACTA\t1\tt\t2\t96.0\t2.0\t0.01\tGACTA\t96.0\t1.5\t0.5\t2005\t2010\nENST1\t102\tACTAG\t1\tt\t3\t105.0\t2.0\t0.01\tACTAG\t105.0\t1.5\t0.5\t2010\t2015" + }, + { + "read_index": 2, + "events": "ENST1\t101\tGACTA\t2\tt\t1\t94.0\t2.0\t0.01\tGACTA\t94.0\t1.5\t0.5\t3000\t3005\nENST1\t102\tACTAG\t2\tt\t2\t106.0\t2.0\t0.01\tACTAG\t106.0\t1.5\t0.5\t3005\t3010\nENST1\t103\tCTAGT\t2\tt\t3\t110.0\t2.0\t0.01\tCTAGT\t110.0\t1.5\t0.5\t3010\t3015" + } + ], + "expected": { + "ENST1": { + "102": { + "GGACT": [ + 120.0, + 121.0 + ] + }, + "103": { + "GACTA": [ + 95.0, + 96.0, + 94.0 + ] + }, + "104": { + "ACTAG": [ + 104.0, + 105.0, + 106.0 + ] + }, + "105": { + "CTAGT": [ + 110.0 + ] + } + } + } +} diff --git a/tests/test_combine.py b/tests/test_combine.py new file mode 100644 index 0000000..cb6b09f --- /dev/null +++ b/tests/test_combine.py @@ -0,0 +1,102 @@ +"""Tests for dataprep.combine -- per-read event merging and kmer-column selection. + +combine() filters events where the reference and model kmers agree, merges multiple events +at the same position into one length-weighted norm_mean, and (new in v2.2) chooses which kmer +column to keep based on --kmer_source, accepting reverse-complement matches for genome reads. +""" +from conftest import make_eventalign_line, make_eventalign_str + +from xpore.scripts.dataprep import combine + +def test_transcriptome_basic_match(): + # reference_kmer == model_kmer: the transcriptome case, works with the default source. + events_str = make_eventalign_line( + contig='ENST1', position=100, read_index=0, + reference_kmer='GGACT', model_kmer='GGACT', + event_level_mean=120.5, start_idx=0, end_idx=5, + ) + result = combine(events_str, kmer_source='reference_kmer') + + assert result is not None + np_events, kmer_col = result + assert kmer_col == 'reference_kmer' + assert np_events.size == 1 + assert np_events['transcript_id'][0] == 'ENST1' + # combine reports the middle base of the 5-mer: position + 2. + assert np_events['transcriptomic_position'][0] == 102 + assert np_events['reference_kmer'][0] == 'GGACT' + assert np_events['norm_mean'][0] == 120.5 + +def test_reverse_complement_dropped_with_reference_kmer(): + # The same reverse read under reference_kmer source: no exact match, so it is filtered out + # and combine returns None. This is the behaviour in v2.1. + events_str = make_eventalign_line(reference_kmer='GGACT', model_kmer='AGTCC') + assert combine(events_str, kmer_source='reference_kmer') is None + +def test_reverse_complement_kept_with_model_kmer(): + # A reverse-oriented read: reference_kmer is the reverse complement of model_kmer (from the read) + # (revcomp('GGACT') == 'AGTCC'). With --kmer_source model_kmer this read must be KEPT. + events_str = make_eventalign_line( + contig='chr1', position=100, read_index=5, + reference_kmer='GGACT', model_kmer='AGTCC', + event_level_mean=88.0, start_idx=2000, end_idx=2005, + ) + result = combine(events_str, kmer_source='model_kmer') + + assert result is not None + np_events, kmer_col = result + assert kmer_col == 'model_kmer' + assert np_events.size == 1 + # the kmer recorded is the model_kmer (the kmer from the read). + assert np_events['model_kmer'][0] == 'AGTCC' + assert np_events['norm_mean'][0] == 88.0 + + +def test_direct_match_with_model_kmer(): + # Test forward reads (reference == model) are kept when kmer_source='model_kmer'. + events_str = make_eventalign_line(reference_kmer='GGACT', model_kmer='GGACT') + result = combine(events_str, kmer_source='model_kmer') + + assert result is not None + np_events, kmer_col = result + assert kmer_col == 'model_kmer' + assert np_events['model_kmer'][0] == 'GGACT' + + +def test_no_match_returns_none(): + # Neither an exact nor a reverse-complement match -> nothing survives the filter. + events_str = make_eventalign_line(reference_kmer='GGACT', model_kmer='TTTTT') + assert combine(events_str, kmer_source='model_kmer') is None + + +def test_norm_mean_is_length_weighted_ref_kmer(): + # when kmer_source='reference_kmer' + # Two events at the same position within a **forward** read (read, contig, position, kmer) group with different (start,end) + # widths -> norm_mean is the length-weighted mean, not the plain mean: + # (100*1 + 130*3) / (1 + 3) == 122.5 (the plain mean would be 115.0) + events_str = make_eventalign_str([ + dict(read_index=0, position=100, reference_kmer='GGACT', model_kmer='GGACT', + event_level_mean=100.0, start_idx=0, end_idx=1), + dict(read_index=0, position=100, reference_kmer='GGACT', model_kmer='GGACT', + event_level_mean=130.0, start_idx=1, end_idx=4), + ]) + np_events, _ = combine(events_str, kmer_source='reference_kmer') + + assert np_events.size == 1 + assert np_events['norm_mean'][0] == 122.5 + +def test_norm_mean_is_length_weighted_model_kmer(): + # when kmer_source='model_kmer' + # Two events at the same position within a **reverse** read (same read, contig, position, kmer) group with different (start,end) + # widths -> norm_mean is the length-weighted mean, not the plain mean: + # (100*1 + 130*3) / (1 + 3) == 122.5 (the plain mean would be 115.0) + events_str = make_eventalign_str([ + dict(contig='chr1', read_index=0, position=100, reference_kmer='GGACT', model_kmer='AGTCC', + event_level_mean=100.0, start_idx=0, end_idx=1), + dict(contig='chr1', read_index=0, position=100, reference_kmer='GGACT', model_kmer='AGTCC', + event_level_mean=130.0, start_idx=1, end_idx=4), + ]) + np_events, _ = combine(events_str, kmer_source='model_kmer') + + assert np_events.size == 1 + assert np_events['norm_mean'][0] == 122.5 diff --git a/tests/test_load_data.py b/tests/test_load_data.py new file mode 100644 index 0000000..442160a --- /dev/null +++ b/tests/test_load_data.py @@ -0,0 +1,94 @@ +"""Tests for diffmod.io.load_data -- assembling per-site read data for modelling. + +load_data() collects, per (position, kmer) present in every run, the reads for each +condition, applies the min/max read-count rules, and returns model-ready arrays. These tests +cover the v2.2 change (over-cap runs are truncated, not dropped) alongside the min_count +filter, output shapes, the pooled round-robin selection (when users have multiple samples per condition), and the >=2-conditions requirement. +""" +import numpy as np # noqa: F401 (imported for parity with the module under test) + +from xpore.diffmod.io import load_data + +IDX = 'geneX' +POS = 100 +KMER = 'GGACT' + + +def _wrap(values): + """Nest a read-value list into the {idx: {pos: {kmer: values}}} shape load_data reads.""" + return {IDX: {POS: {KMER: list(values)}}} + + +def test_over_max_count_truncates_not_drops(): + # A site in a run has MORE than max_count reads. The v2.1 behaviour dropped the site. + # the new behaviour truncates each site to max_count. + data_dict = { + ('KO', 'r1'): _wrap([0.1, 0.2, 0.3, 0.4, 0.5]), # 5 > 3 + ('WT', 'r2'): _wrap([1.1, 1.2, 1.3, 1.4]), # 4 > 3 + } + data = load_data(IDX, data_dict, min_count=1, max_count=3) + + key = (IDX, POS, KMER) + assert key in data # site kept, not dropped + result = data[key] + assert len(result['y']) == 6 # 3 (KO) + 3 (WT) + assert (result['y_condition_names'] == 'KO').sum() == 3 + assert (result['y_condition_names'] == 'WT').sum() == 3 + + +def test_run_below_min_count_dropped_site_survives(): + # KO has one run+site below min_count and one above. The below-min site+run is dropped. + data_dict = { + ('KO', 'r1'): _wrap([0.1]), # 1 < min_count -> dropped + ('KO', 'r2'): _wrap([0.2, 0.3, 0.4]), # 3 >= min_count -> kept + ('WT', 'r3'): _wrap([1.1, 1.2, 1.3]), + } + data = load_data(IDX, data_dict, min_count=2, max_count=100) + + result = data[(IDX, POS, KMER)] + ko_values = result['y'][result['y_condition_names'] == 'KO'] + assert sorted(ko_values) == [0.2, 0.3, 0.4] # r1's single read excluded + assert 'r1' not in result['run_names'] + + +def test_basic_two_condition_shapes(): + data_dict = { + ('KO', 'r1'): _wrap([0.1, 0.2, 0.3]), + ('WT', 'r2'): _wrap([1.1, 1.2, 1.3, 1.4]), + } + data = load_data(IDX, data_dict, min_count=1, max_count=100) + + result = data[(IDX, POS, KMER)] + assert len(result['y']) == 7 + assert result['x'].shape == (7, 2) # one-hot over 2 conditions + assert result['r'].shape == (7, 2) # one-hot over 2 runs + assert result['condition_names'] == ['KO', 'WT'] + assert result['run_names'] == ['r1', 'r2'] + + +def test_pooling_over_cap_round_robin(): + # KO's two runs total 7 reads > max_count 4. Round-robin selection takes reads evenly + # across the runs (2 + 2), not the first 4 (which would be 3 from r1 + 1 from r2). + data_dict = { + ('KO', 'r1'): _wrap([1.0, 2.0, 3.0]), + ('KO', 'r2'): _wrap([10.0, 20.0, 30.0, 40.0]), + ('WT', 'r3'): _wrap([100.0, 200.0]), + } + data = load_data(IDX, data_dict, min_count=1, max_count=4, pooling=True) + + result = data[(IDX, POS, KMER)] + ko = result['y_condition_names'] == 'KO' + assert ko.sum() == 4 # capped at max_count + ko_runs = result['y_run_names'][ko] + assert (ko_runs == 'r1').sum() == 2 # evenly split -> round-robin + assert (ko_runs == 'r2').sum() == 2 + + +def test_single_condition_site_skipped(): + # Only one condition present -> fewer than 2 conditions -> the site is skipped entirely. + data_dict = { + ('KO', 'r1'): _wrap([0.1, 0.2, 0.3]), + ('KO', 'r2'): _wrap([0.4, 0.5, 0.6]), + } + data = load_data(IDX, data_dict, min_count=1, max_count=100) + assert len(data) == 0 diff --git a/tests/test_preprocess_tx.py b/tests/test_preprocess_tx.py new file mode 100644 index 0000000..3bc0bac --- /dev/null +++ b/tests/test_preprocess_tx.py @@ -0,0 +1,47 @@ +"""Tests for dataprep.preprocess_tx -- per-transcript reshaping into per-position signal. + +preprocess_tx() concatenates a transcript's reads, groups them by position (or by +(position, kmer) in model_kmer mode), caps reads per site, runs a kmer-consistency check, +and appends the result to data.json. These tests exercise the v2.2 changes: the per-site +read cap, the break->continue fix, and the model_kmer (position, kmer) grouping. +""" +from conftest import make_events_array, read_json_output + +from xpore.scripts.dataprep import preprocess_tx + + +def test_per_site_readcount_max_caps_reads(out_paths, locks): + # Five reads at one site; a per-site cap of 3 keeps only 3 values for that site. + data_dict = { + i: make_events_array([(102, 'GGACT', float(i))], kmer_col='reference_kmer') + for i in range(5) + } + preprocess_tx('tx1', data_dict, 'reference_kmer', 3, out_paths, locks) + + data = read_json_output(out_paths['json'])['tx1'] + assert len(data['102']['GGACT']) == 3 + + +def test_readcount_max_none_keeps_all_reads(out_paths, locks): + # readcount_max=None means no per-site cap -> every read is kept. + data_dict = { + i: make_events_array([(102, 'GGACT', float(i))], kmer_col='reference_kmer') + for i in range(5) + } + preprocess_tx('tx1', data_dict, 'reference_kmer', None, out_paths, locks) + + data = read_json_output(out_paths['json'])['tx1'] + assert len(data['102']['GGACT']) == 5 + + +def test_model_kmer_mode_keeps_both_orientations(out_paths, locks): + # Forward and reverse reads at the same position carry different model_kmers. In + # model_kmer mode reads are grouped by (position, kmer), so both orientations are kept. + data_dict = { + 0: make_events_array([(102, 'GGACT', 100.0)], kmer_col='model_kmer'), + 1: make_events_array([(102, 'AGTCC', 88.0)], kmer_col='model_kmer'), + } + preprocess_tx('tx1', data_dict, 'model_kmer', None, out_paths, locks) + + data = read_json_output(out_paths['json'])['tx1'] + assert set(data['102'].keys()) == {'GGACT', 'AGTCC'} diff --git a/tests/test_xpore_v2_1_back_compatibility.py b/tests/test_xpore_v2_1_back_compatibility.py new file mode 100644 index 0000000..ebb1621 --- /dev/null +++ b/tests/test_xpore_v2_1_back_compatibility.py @@ -0,0 +1,67 @@ +"""Back-compatibility tests: the transcriptome path still matches xPore v2.1. + +Freezes outputs captured from xPore v2.1 and asserts the current (v2.2) +code reproduces them -- i.e. the v2.2 changes did not alter transcriptome (reference_kmer) +behaviour. + +Two levels are pinned: + * combine() -- the per-read event merge (fixtures/_regenerate_combine_v2_1.py) + * preprocess_tx() -- the per-position reshape to data.json + (combine + preprocess_tx; fixtures/_regenerate_data_json_v2_1.py) + +The data.json comparison is order-insensitive within each site: preprocess_tx orders the +reads at a position via np.argsort, whose tie-breaking is not guaranteed stable across numpy +versions, and read order at a site is not meaningful to diffmod. What must be preserved is the +SET of values at each (position, kmer). +""" +import json +from pathlib import Path + +from xpore.scripts.dataprep import combine, preprocess_tx + +FIXTURES = Path(__file__).parent / "fixtures" + + +def test_combine_transcriptome_matches_v2_1(): + reference = json.loads((FIXTURES / "combine_transcriptome_v2_1.json").read_text()) + + np_events, kmer_col = combine(reference["input_eventalign"], kmer_source="reference_kmer") + + assert kmer_col == reference["kmer_col"] + actual = [ + [str(row["transcript_id"]), int(row["transcriptomic_position"]), + str(row[kmer_col]), float(row["norm_mean"])] + for row in np_events + ] + assert actual == reference["expected"] + + +def _sort_site_values(data_json): + """Sort the read values at each (position, kmer) so comparison ignores within-site order.""" + return { + tx_id: {pos: {kmer: sorted(vals) for kmer, vals in kmers.items()} + for pos, kmers in positions.items()} + for tx_id, positions in data_json.items() + } + + +def test_preprocess_tx_data_json_matches_v2_1(out_paths, locks): + reference = json.loads((FIXTURES / "data_json_transcriptome_v2_1.json").read_text()) + + # Rebuild the per-read arrays the pipeline feeds preprocess_tx, using the current combine(). + data_dict, kmer_col = {}, "reference_kmer" + for read in reference["reads"]: + np_events, kmer_col = combine(read["events"], kmer_source="reference_kmer") + if np_events.size > 1: # mirror the caller's inclusion rule + data_dict[read["read_index"]] = np_events + + # readcount_max=None matches v2.1 (no per-site cap), so both reshape the same reads. + preprocess_tx(reference["tx_id"], data_dict, kmer_col, None, out_paths, locks) + + actual = {} + with open(out_paths["json"]) as f: + for line in f: + if line.strip(): + actual.update(json.loads(line)) + + assert _sort_site_values(actual) == _sort_site_values(reference["expected"]) diff --git a/xpore/__init__.py b/xpore/__init__.py index ab02866..d3c8880 100644 --- a/xpore/__init__.py +++ b/xpore/__init__.py @@ -1 +1 @@ -__version__ = "2.1" +__version__ = "2.2" diff --git a/xpore/diffmod/model_kmer.csv b/xpore/diffmod/RNA002_5mer_model.csv similarity index 100% rename from xpore/diffmod/model_kmer.csv rename to xpore/diffmod/RNA002_5mer_model.csv diff --git a/xpore/diffmod/RNA004_5mer_model.csv b/xpore/diffmod/RNA004_5mer_model.csv new file mode 100644 index 0000000..64aafa0 --- /dev/null +++ b/xpore/diffmod/RNA004_5mer_model.csv @@ -0,0 +1,1025 @@ +model_kmer,model_mean,model_stdv +AAAAA,95.82620936892796,3.0 +AAAAC,93.94534245544321,3.0 +AAAAG,96.75244217476782,3.0 +AAAAT,100.11003605201408,3.0 +AAACA,80.52107965758883,3.0 +AAACC,78.0608562659247,3.0 +AAACG,74.00812263706015,3.0 +AAACT,84.49851284450727,3.0 +AAAGA,108.77035360010254,3.0 +AAAGC,95.18957240215708,3.0 +AAAGG,97.60129573899219,3.0 +AAAGT,101.5521989303327,3.0 +AAATA,94.09832310112735,3.0 +AAATC,95.44718824854327,3.0 +AAATG,78.18692355388767,3.0 +AAATT,95.18957240215708,3.0 +AACAA,74.44112929936631,3.0 +AACAC,69.86647833170001,3.0 +AACAG,80.3248699230055,3.0 +AACAT,73.63390162534412,3.0 +AACCA,69.46561259687189,3.0 +AACCC,64.69771688982348,3.0 +AACCG,72.75552179698519,3.0 +AACCT,66.27001768740493,3.0 +AACGA,74.65703773399467,3.0 +AACGC,69.7031637746776,3.0 +AACGG,76.90415569727837,3.0 +AACGT,71.48891031799577,3.0 +AACTA,79.75890617143183,3.0 +AACTC,76.2959133559782,3.0 +AACTG,80.93571164235745,3.0 +AACTT,73.80230108940098,3.0 +AAGAA,109.93333677501815,3.0 +AAGAC,111.01234044909071,3.0 +AAGAG,110.6043524711086,3.0 +AAGAT,115.84039467189294,3.0 +AAGCA,100.78068665803413,3.0 +AAGCC,98.40977100906343,3.0 +AAGCG,99.70462988278295,3.0 +AAGCT,103.63028142087096,3.0 +AAGGA,104.44476240503103,3.0 +AAGGC,99.6576847110738,3.0 +AAGGG,102.55951159961332,3.0 +AAGGT,103.65252854799547,3.0 +AAGTA,108.46386821649175,3.0 +AAGTC,109.156603082218,3.0 +AAGTG,100.90005810977325,3.0 +AAGTT,107.8127511590893,3.0 +AATAA,78.36647407746253,3.0 +AATAC,73.22631984959625,3.0 +AATAG,85.59338036734157,3.0 +AATAT,78.07074391702294,3.0 +AATCA,72.80660084947708,3.0 +AATCC,67.96438524315963,3.0 +AATCG,74.93763932091056,3.0 +AATCT,71.43499083664527,3.0 +AATGA,80.08297649250969,3.0 +AATGC,73.33906082291048,3.0 +AATGG,82.79538558676504,3.0 +AATGT,76.21186832164322,3.0 +AATTA,77.61352043328077,3.0 +AATTC,73.66341014448489,3.0 +AATTG,79.69973171221524,3.0 +AATTT,72.51660390327108,3.0 +ACAAA,92.09052976087104,3.0 +ACAAC,81.02069839121474,3.0 +ACAAG,93.29942432983735,3.0 +ACAAT,90.80404072635456,3.0 +ACACA,72.26175803284734,3.0 +ACACC,68.72027981260051,3.0 +ACACG,70.2446803712819,3.0 +ACACT,75.89331107146006,3.0 +ACAGA,110.60889659506475,3.0 +ACAGC,95.70757424658491,3.0 +ACAGG,99.8293701614573,3.0 +ACAGT,104.37615516391364,3.0 +ACATA,77.91133448848697,3.0 +ACATC,77.60091814951757,3.0 +ACATG,73.22913199188059,3.0 +ACATT,79.01580960776492,3.0 +ACCAA,77.09862782801363,3.0 +ACCAC,68.43134809310136,3.0 +ACCAG,85.20038480080204,3.0 +ACCAT,75.75448724127297,3.0 +ACCCA,65.90900334339925,3.0 +ACCCC,60.411426990766685,3.0 +ACCCG,70.73268298599793,3.0 +ACCCT,65.27006268990887,3.0 +ACCGA,85.08943589616831,3.0 +ACCGC,77.92930770750148,3.0 +ACCGG,85.74661752481722,3.0 +ACCGT,80.64078000663874,3.0 +ACCTA,72.51660390327108,3.0 +ACCTC,67.75348001831705,3.0 +ACCTG,78.98454205022914,3.0 +ACCTT,68.75502749723665,3.0 +ACGAA,93.94534245544321,3.0 +ACGAC,93.4139617641573,3.0 +ACGAG,95.5996548164863,3.0 +ACGAT,97.26465423181352,3.0 +ACGCA,83.87508236975374,3.0 +ACGCC,81.40712216079402,3.0 +ACGCG,79.66821502199426,3.0 +ACGCT,86.45611051671538,3.0 +ACGGA,107.32407774846914,3.0 +ACGGC,94.8851173269455,3.0 +ACGGG,99.62747096019936,3.0 +ACGGT,103.529498992139,3.0 +ACGTA,86.50609535315411,3.0 +ACGTC,86.06678289310165,3.0 +ACGTG,79.61713596950237,3.0 +ACGTT,85.3337726893249,3.0 +ACTAA,82.87767966458334,3.0 +ACTAC,73.0982519453524,3.0 +ACTAG,89.24316032567219,3.0 +ACTAT,81.98169241005236,3.0 +ACTCA,71.07672485836207,3.0 +ACTCC,65.12462682906272,3.0 +ACTCG,73.88338235838582,3.0 +ACTCT,71.6906577204949,3.0 +ACTGA,95.0303715212215,3.0 +ACTGC,83.48838750586327,3.0 +ACTGG,92.99983735648848,3.0 +ACTGT,88.879463676515,3.0 +ACTTA,74.83791702378822,3.0 +ACTTC,69.96814555041973,3.0 +ACTTG,79.54477399311328,3.0 +ACTTT,71.70409542451063,3.0 +AGAAA,96.52008272534536,3.0 +AGAAC,88.5396853463734,3.0 +AGAAG,94.51404631516799,3.0 +AGAAT,97.1662028979103,3.0 +AGACA,76.42198064394117,3.0 +AGACC,74.93763932091056,3.0 +AGACG,70.0045521284364,3.0 +AGACT,81.91305863929077,3.0 +AGAGA,110.2930529233885,3.0 +AGAGC,96.68062203561509,3.0 +AGAGG,103.07261881253424,3.0 +AGAGT,105.79819486788014,3.0 +AGATA,85.64116412177741,3.0 +AGATC,88.95415973451395,3.0 +AGATG,71.87533426326806,3.0 +AGATT,88.72748114280871,3.0 +AGCAA,76.78885083589994,3.0 +AGCAC,70.07202298717874,3.0 +AGCAG,85.78789343454767,3.0 +AGCAT,75.40572344362867,3.0 +AGCCA,68.9127515027479,3.0 +AGCCC,64.50261279133882,3.0 +AGCCG,72.14442795816693,3.0 +AGCCT,66.64844338418587,3.0 +AGCGA,81.66473818287321,3.0 +AGCGC,73.65591455368912,3.0 +AGCGG,82.95281759049715,3.0 +AGCGT,76.08842746845527,3.0 +AGCTA,78.85884951545201,3.0 +AGCTC,75.75448724127297,3.0 +AGCTG,79.2752192797617,3.0 +AGCTT,73.90063365483088,3.0 +AGGAA,106.26705572938617,3.0 +AGGAC,105.87800324102278,3.0 +AGGAG,105.83769592684486,3.0 +AGGAT,111.96036691840172,3.0 +AGGCA,95.50221301460643,3.0 +AGGCC,93.77268507134005,3.0 +AGGCG,95.45272837616118,3.0 +AGGCT,98.86691297124995,3.0 +AGGGA,104.23296360555031,3.0 +AGGGC,96.77039729740251,3.0 +AGGGG,101.0975635802898,3.0 +AGGGT,102.46098875865587,3.0 +AGGTA,105.59895302336659,3.0 +AGGTC,104.14770031589859,3.0 +AGGTG,96.44610700867878,3.0 +AGGTT,104.54786398146717,3.0 +AGTAA,80.33006850341606,3.0 +AGTAC,72.39198151460504,3.0 +AGTAG,88.3627512430721,3.0 +AGTAT,80.11604718738549,3.0 +AGTCA,71.76279674245697,3.0 +AGTCC,66.89792394153456,3.0 +AGTCG,73.88338235838582,3.0 +AGTCT,71.624643532481,3.0 +AGTGA,87.72955856202796,3.0 +AGTGC,78.06600178704689,3.0 +AGTGG,88.1519746255113,3.0 +AGTGT,82.56639364522486,3.0 +AGTTA,75.73957266183288,3.0 +AGTTC,72.26175803284734,3.0 +AGTTG,76.92624997148602,3.0 +AGTTT,72.14442795816693,3.0 +ATAAA,91.19504551542164,3.0 +ATAAC,84.112089074928,3.0 +ATAAG,92.18594090296568,3.0 +ATAAT,91.65320878106304,3.0 +ATACA,74.13286291573449,3.0 +ATACC,70.9711279141348,3.0 +ATACG,70.62722958295811,3.0 +ATACT,77.80049978430839,3.0 +ATAGA,107.26157302828793,3.0 +ATAGC,92.67066157229425,3.0 +ATAGG,98.50177669358607,3.0 +ATAGT,101.75292274536783,3.0 +ATATA,81.88143846588814,3.0 +ATATC,82.13703384268355,3.0 +ATATG,74.19387723212975,3.0 +ATATT,83.33074853576787,3.0 +ATCAA,76.46671419206696,3.0 +ATCAC,69.29917527232718,3.0 +ATCAG,83.48175228372718,3.0 +ATCAT,76.168460552955,3.0 +ATCCA,67.16484150283765,3.0 +ATCCC,61.488729780829836,3.0 +ATCCG,69.94617619394069,3.0 +ATCCT,65.9000013607538,3.0 +ATCGA,83.51358750537146,3.0 +ATCGC,76.63372480289642,3.0 +ATCGG,82.90287158221317,3.0 +ATCGT,79.69973171221524,3.0 +ATCTA,74.70910910059513,3.0 +ATCTC,70.35041435605532,3.0 +ATCTG,77.23889708457665,3.0 +ATCTT,71.07672485836207,3.0 +ATGAA,100.16737627622827,3.0 +ATGAC,102.01925595172668,3.0 +ATGAG,99.6576847110738,3.0 +ATGAT,105.68235449421366,3.0 +ATGCA,89.64224347924223,3.0 +ATGCC,87.32155371751709,3.0 +ATGCG,85.9659341402592,3.0 +ATGCT,92.2671555229435,3.0 +ATGGA,104.6032865165002,3.0 +ATGGC,93.10956255932223,3.0 +ATGGG,98.74362022726764,3.0 +ATGGT,101.582164778373,3.0 +ATGTA,95.83219663525095,3.0 +ATGTC,95.13565274511357,3.0 +ATGTG,82.8465198068614,3.0 +ATGTT,94.15001655206783,3.0 +ATTAA,80.37665859955627,3.0 +ATTAC,72.38267979968812,3.0 +ATTAG,86.14699448169739,3.0 +ATTAT,80.65019469215972,3.0 +ATTCA,70.69907695452707,3.0 +ATTCC,65.55854041867562,3.0 +ATTCG,73.0102000562794,3.0 +ATTCT,71.2423747268454,3.0 +ATTGA,90.1606956248493,3.0 +ATTGC,80.65019469215972,3.0 +ATTGG,89.54929976975045,3.0 +ATTGT,85.46993951415361,3.0 +ATTTA,74.32268515532284,3.0 +ATTTC,70.4778372912672,3.0 +ATTTG,76.64122812418454,3.0 +ATTTT,71.64455323982851,3.0 +CAAAA,94.58179898549498,3.0 +CAAAC,91.89613001565503,3.0 +CAAAG,94.85145155985197,3.0 +CAAAT,98.82733531013379,3.0 +CAACA,80.2027014385807,3.0 +CAACC,77.3404713616952,3.0 +CAACG,73.93626138574356,3.0 +CAACT,83.57343767989641,3.0 +CAAGA,107.94017409430117,3.0 +CAAGC,94.82781399826297,3.0 +CAAGG,97.84770008746567,3.0 +CAAGT,101.0411334458565,3.0 +CAATA,91.72124906084396,3.0 +CAATC,92.7422215101023,3.0 +CAATG,78.12034872948558,3.0 +CAATT,92.83705813606113,3.0 +CACAA,73.79083307969675,3.0 +CACAC,69.07618517962966,3.0 +CACAG,79.0151725449183,3.0 +CACAT,73.42355527995966,3.0 +CACCA,68.87594908838443,3.0 +CACCC,63.563110771077795,3.0 +CACCG,72.13958972411555,3.0 +CACCT,65.78683485897446,3.0 +CACGA,75.56030377414967,3.0 +CACGC,71.1844552427842,3.0 +CACGG,76.94980109211512,3.0 +CACGT,72.7504314434675,3.0 +CACTA,78.04912893334082,3.0 +CACTC,74.58535235137897,3.0 +CACTG,79.17445108216339,3.0 +CACTT,72.92750697964016,3.0 +CAGAA,111.13084380167784,3.0 +CAGAC,110.74347814188324,3.0 +CAGAG,110.48863227145951,3.0 +CAGAT,117.29301602789243,3.0 +CAGCA,99.9569137977656,3.0 +CAGCC,97.76573139309797,3.0 +CAGCG,98.62269846042685,3.0 +CAGCT,102.19943651042783,3.0 +CAGGA,104.20612860691067,3.0 +CAGGC,99.71698039846584,3.0 +CAGGG,102.980577286638,3.0 +CAGGT,103.19735909120858,3.0 +CAGTA,110.04340176812916,3.0 +CAGTC,107.6810808454116,3.0 +CAGTG,101.95269729108132,3.0 +CAGTT,108.31702048663416,3.0 +CATAA,77.87507268150881,3.0 +CATAC,71.77308462222717,3.0 +CATAG,83.25691056241946,3.0 +CATAT,77.86136493733994,3.0 +CATCA,72.74055925335396,3.0 +CATCC,67.32412004008273,3.0 +CATCG,75.19500831862065,3.0 +CATCT,71.07672485836207,3.0 +CATGA,82.47756062503876,3.0 +CATGC,74.63182420612486,3.0 +CATGG,82.91374100583951,3.0 +CATGT,77.78424238320898,3.0 +CATTA,77.15645999264991,3.0 +CATTC,72.50863991491627,3.0 +CATTG,78.79153221239858,3.0 +CATTT,71.63142124164652,3.0 +CCAAA,84.47485437545035,3.0 +CCAAC,80.2027014385807,3.0 +CCAAG,84.59577614229113,3.0 +CCAAT,84.50275934450806,3.0 +CCACA,73.29722392540583,3.0 +CCACC,69.47344481547077,3.0 +CCACG,70.30647116292555,3.0 +CCACT,76.33530935052418,3.0 +CCAGA,102.28259657908205,3.0 +CCAGC,92.22021035123434,3.0 +CCAGG,90.3426754556443,3.0 +CCAGT,93.40646634905454,3.0 +CCATA,78.56512559346956,3.0 +CCATC,77.8683663037045,3.0 +CCATG,74.13286291573449,3.0 +CCATT,79.1855469740769,3.0 +CCCAA,77.17838472310495,3.0 +CCCAC,68.19546793149372,3.0 +CCCAG,84.83761950027969,3.0 +CCCAT,75.79164385274034,3.0 +CCCCA,66.81864545850372,3.0 +CCCCC,61.29389853358644,3.0 +CCCCG,71.28441121160702,3.0 +CCCCT,64.36153727956257,3.0 +CCCGA,81.76438194522102,3.0 +CCCGC,77.4759036140579,3.0 +CCCGG,84.48631008664402,3.0 +CCCGT,79.76696766940601,3.0 +CCCTA,73.06712037420472,3.0 +CCCTC,67.73074903264879,3.0 +CCCTG,79.58108007472944,3.0 +CCCTT,68.43907067926332,3.0 +CCGAA,95.48840688235315,3.0 +CCGAC,96.08334188384943,3.0 +CCGAG,96.4317400643443,3.0 +CCGAT,99.41876664166666,3.0 +CCGCA,84.89465805959932,3.0 +CCGCC,82.67568523771513,3.0 +CCGCG,80.84303217514953,3.0 +CCGCT,88.02327361151352,3.0 +CCGGA,103.14825922045063,3.0 +CCGGC,90.99656705702739,3.0 +CCGGG,93.87118524789948,3.0 +CCGGT,103.20428912621193,3.0 +CCGTA,88.72055110780536,3.0 +CCGTC,88.46527197842416,3.0 +CCGTG,80.60535853931731,3.0 +CCGTT,86.67762461245358,3.0 +CCTAA,80.24513903047934,3.0 +CCTAC,72.26175803284734,3.0 +CCTAG,87.23745105587551,3.0 +CCTAT,79.97022585608002,3.0 +CCTCA,71.12715191410166,3.0 +CCTCC,65.25349782604565,3.0 +CCTCG,73.8337406503915,3.0 +CCTCT,70.98752885642166,3.0 +CCTGA,89.85400204502413,3.0 +CCTGC,80.31998319768395,3.0 +CCTGG,88.14901507679333,3.0 +CCTGT,83.99590943806326,3.0 +CCTTA,74.84301827027241,3.0 +CCTTC,69.79293903646503,3.0 +CCTTG,79.09848739911249,3.0 +CCTTT,71.25321920206463,3.0 +CGAAA,95.4415263656709,3.0 +CGAAC,88.02017693440288,3.0 +CGAAG,93.75500543534174,3.0 +CGAAT,97.4607378188172,3.0 +CGACA,76.42890734077737,3.0 +CGACC,74.5309704460317,3.0 +CGACG,69.81440714079255,3.0 +CGACT,81.18004843551402,3.0 +CGAGA,108.74218878108616,3.0 +CGAGC,97.09702623575433,3.0 +CGAGG,101.00511409522922,3.0 +CGAGT,104.34144209080004,3.0 +CGATA,82.8465198068614,3.0 +CGATC,87.04281605772194,3.0 +CGATG,71.49722059726913,3.0 +CGATT,86.9308028926929,3.0 +CGCAA,76.6029683363206,3.0 +CGCAC,69.20754908470158,3.0 +CGCAG,83.84578942593366,3.0 +CGCAT,74.90991399795341,3.0 +CGCCA,68.26571614530484,3.0 +CGCCC,62.781493164044214,3.0 +CGCCG,71.48891031799577,3.0 +CGCCT,65.34906498353129,3.0 +CGCGA,81.61035258797277,3.0 +CGCGC,73.79083307969675,3.0 +CGCGG,82.28896120878535,3.0 +CGCGT,76.29868403471268,3.0 +CGCTA,78.25063493364709,3.0 +CGCTC,73.89628648273657,3.0 +CGCTG,78.70215542235304,3.0 +CGCTT,72.63597904456337,3.0 +CGGAA,107.03141061992918,3.0 +CGGAC,106.31586728791818,3.0 +CGGAG,106.41109904745176,3.0 +CGGAT,112.60331056138395,3.0 +CGGCA,96.20426365069022,3.0 +CGGCC,94.03030091772624,3.0 +CGGCG,95.61388103072146,3.0 +CGGCT,98.73590185773787,3.0 +CGGGA,103.81970483733139,3.0 +CGGGC,98.02334723948415,3.0 +CGGGG,101.14693999184219,3.0 +CGGGT,102.07871995302027,3.0 +CGGTA,107.35028815846573,3.0 +CGGTC,104.91320213039712,3.0 +CGGTG,94.0913152341215,3.0 +CGGTT,104.4269422146264,3.0 +CGTAA,79.69973171221524,3.0 +CGTAC,71.50535395368355,3.0 +CGTAG,87.17453415989777,3.0 +CGTAT,79.85129592086973,3.0 +CGTCA,71.88753702113131,3.0 +CGTCC,65.90900334339925,3.0 +CGTCG,73.33906082291048,3.0 +CGTCT,71.1734623069733,3.0 +CGTGA,87.01421095202777,3.0 +CGTGC,78.18692355388767,3.0 +CGTGG,87.90419495217203,3.0 +CGTGT,82.79538558676504,3.0 +CGTTA,75.38712001379486,3.0 +CGTTC,71.26500117517396,3.0 +CGTTG,76.90894666991062,3.0 +CGTTT,71.6106409754449,3.0 +CTAAA,91.39510503519394,3.0 +CTAAC,84.36693494535173,3.0 +CTAAG,91.41520009866477,3.0 +CTAAT,92.90427019900284,3.0 +CTACA,75.39460418454512,3.0 +CTACC,72.49186158404845,3.0 +CTACG,71.61771824118887,3.0 +CTACT,78.93608240479193,3.0 +CTAGA,107.29185740775154,3.0 +CTAGC,92.65535098081438,3.0 +CTAGG,96.08984305222052,3.0 +CTAGT,102.44891706777652,3.0 +CTATA,83.34755163934979,3.0 +CTATC,83.42311632357976,3.0 +CTATG,75.42184626388405,3.0 +CTATT,83.9804187612503,3.0 +CTCAA,77.50085184548578,3.0 +CTCAC,69.72240145626992,3.0 +CTCAG,84.112089074928,3.0 +CTCAT,76.75240982205372,3.0 +CTCCA,68.03787762836149,3.0 +CTCCC,61.658830304589124,3.0 +CTCCG,70.64013388300188,3.0 +CTCCT,66.63584338443178,3.0 +CTCGA,82.56639364522486,3.0 +CTCGC,77.23125162764516,3.0 +CTCGG,82.56639364522486,3.0 +CTCGT,80.33006850341606,3.0 +CTCTA,74.85401788241761,3.0 +CTCTC,70.87501716248273,3.0 +CTCTG,77.7593342097869,3.0 +CTCTT,71.15078947569066,3.0 +CTGAA,101.65969458889627,3.0 +CTGAC,104.55971904286781,3.0 +CTGAG,101.48359713569855,3.0 +CTGAT,107.47091458536025,3.0 +CTGCA,92.5586015850786,3.0 +CTGCC,90.91722971683892,3.0 +CTGCG,87.74798612316087,3.0 +CTGCT,94.9430368110067,3.0 +CTGGA,104.07732068371757,3.0 +CTGGC,93.29028671218981,3.0 +CTGGG,99.04350623011594,3.0 +CTGGT,101.5320757559806,3.0 +CTGTA,97.48771020938901,3.0 +CTGTC,98.09281678188428,3.0 +CTGTG,85.94461440867684,3.0 +CTGTT,95.56297081921016,3.0 +CTTAA,80.16197878613211,3.0 +CTTAC,72.3878253208103,3.0 +CTTAG,85.96370327829148,3.0 +CTTAT,80.20538655482034,3.0 +CTTCA,71.76279674245697,3.0 +CTTCC,66.03117165213104,3.0 +CTTCG,73.99270522423149,3.0 +CTTCT,71.66325646329081,3.0 +CTTGA,88.32554016677227,3.0 +CTTGC,79.34471447334099,3.0 +CTTGG,88.1412016573463,3.0 +CTTGT,82.75423863570238,3.0 +CTTTA,74.39481767664363,3.0 +CTTTC,70.73646794046698,3.0 +CTTTG,76.93007024024966,3.0 +CTTTT,71.75748870530248,3.0 +GAAAA,93.39688457948306,3.0 +GAAAC,92.88133400390136,3.0 +GAAAG,93.32901314068272,3.0 +GAAAT,98.52795565453424,3.0 +GAACA,81.05393915692218,3.0 +GAACC,77.91512699775537,3.0 +GAACG,75.02683532285096,3.0 +GAACT,83.58540764452421,3.0 +GAAGA,102.81702343573883,3.0 +GAAGC,91.19905166738764,3.0 +GAAGG,95.96242011700865,3.0 +GAAGT,95.57599634742937,3.0 +GAATA,96.23800918240934,3.0 +GAATC,95.59209720605875,3.0 +GAATG,84.112089074928,3.0 +GAATT,96.36186941480715,3.0 +GACAA,70.98752885642166,3.0 +GACAC,67.64636596496662,3.0 +GACAG,75.64756697731025,3.0 +GACAT,70.95822168146795,3.0 +GACCA,68.49119809193331,3.0 +GACCC,63.91266114631842,3.0 +GACCG,70.91790610971864,3.0 +GACCT,64.93851979381928,3.0 +GACGA,69.51594969779143,3.0 +GACGC,67.52162568629228,3.0 +GACGG,73.163413670892,3.0 +GACGT,68.39226219671431,3.0 +GACTA,78.87299543798257,3.0 +GACTC,75.66702550501348,3.0 +GACTG,78.75665100107643,3.0 +GACTT,72.63996007242923,3.0 +GAGAA,111.35422937933613,3.0 +GAGAC,111.09657786726935,3.0 +GAGAG,112.6496209542556,3.0 +GAGAT,117.43381201209311,3.0 +GAGCA,99.76623013535954,3.0 +GAGCC,96.29940387416669,3.0 +GAGCG,100.30363042048191,3.0 +GAGCT,102.13446752011791,3.0 +GAGGA,101.0411334458565,3.0 +GAGGC,92.69753293938652,3.0 +GAGGG,101.62996979166279,3.0 +GAGGT,99.45514932543428,3.0 +GAGTA,111.34263697898231,3.0 +GAGTC,109.78780639133375,3.0 +GAGTG,105.21696427238558,3.0 +GAGTT,109.60825586775888,3.0 +GATAA,73.43376884158754,3.0 +GATAC,69.50586702745444,3.0 +GATAG,80.03455585092024,3.0 +GATAT,73.5784615210103,3.0 +GATCA,71.6106409754449,3.0 +GATCC,66.84086201504488,3.0 +GATCG,73.95466241723229,3.0 +GATCT,69.39273056917943,3.0 +GATGA,74.94219012120101,3.0 +GATGC,68.35237057219305,3.0 +GATGG,76.69066971784285,3.0 +GATGT,70.95822168146795,3.0 +GATTA,75.95702262691249,3.0 +GATTC,72.0148761507784,3.0 +GATTG,77.17838472310495,3.0 +GATTT,70.39766101696604,3.0 +GCAAA,88.2234284447426,3.0 +GCAAC,81.6317873106288,3.0 +GCAAG,89.33642854014956,3.0 +GCAAT,89.07269277921941,3.0 +GCACA,73.11693654535587,3.0 +GCACC,69.54605715439611,3.0 +GCACG,69.8916918595698,3.0 +GCACT,75.98995206460995,3.0 +GCAGA,106.44060756659252,3.0 +GCAGC,91.49059699607221,3.0 +GCAGG,93.66880763458104,3.0 +GCAGT,98.55467996657892,3.0 +GCATA,79.34619503831851,3.0 +GCATC,79.16359922783784,3.0 +GCATG,73.36127333851246,3.0 +GCATT,79.71814609637256,3.0 +GCCAA,74.80111460942476,3.0 +GCCAC,68.19546793149372,3.0 +GCCAG,82.15739525675434,3.0 +GCCAT,74.32268515532284,3.0 +GCCCA,65.47250583671925,3.0 +GCCCC,60.53486801964765,3.0 +GCCCG,69.56316772074184,3.0 +GCCCT,63.92263611684158,3.0 +GCCGA,80.97523677126424,3.0 +GCCGC,75.25689670773018,3.0 +GCCGG,81.87369321532816,3.0 +GCCGT,77.83141367181929,3.0 +GCCTA,72.73577126750284,3.0 +GCCTC,68.15041866303274,3.0 +GCCTG,77.47069484345289,3.0 +GCCTT,68.6551354806687,3.0 +GCGAA,96.36186941480715,3.0 +GCGAC,97.91577128921597,3.0 +GCGAG,97.13443092574884,3.0 +GCGAT,101.13842678714471,3.0 +GCGCA,86.75950352769425,3.0 +GCGCC,84.74242586612881,3.0 +GCGCG,82.86468576110555,3.0 +GCGCT,89.29660894976153,3.0 +GCGGA,105.26429280623796,3.0 +GCGGC,91.19785396815283,3.0 +GCGGG,98.82284477254504,3.0 +GCGGT,101.43180371543657,3.0 +GCGTA,90.03725459596832,3.0 +GCGTC,90.03725459596832,3.0 +GCGTG,80.61935986650236,3.0 +GCGTT,88.37820818682755,3.0 +GCTAA,79.32152967263335,3.0 +GCTAC,71.89143512189744,3.0 +GCTAG,85.35949238875044,3.0 +GCTAT,79.06939509670742,3.0 +GCTCA,69.76695158089547,3.0 +GCTCC,64.4044733125246,3.0 +GCTCG,72.6245233333697,3.0 +GCTCT,70.1411724169185,3.0 +GCTGA,86.43063204378969,3.0 +GCTGC,79.69973171221524,3.0 +GCTGG,87.27999933436912,3.0 +GCTGT,83.2389067728216,3.0 +GCTTA,73.55238815155988,3.0 +GCTTC,69.46686862618532,3.0 +GCTTG,77.3404713616952,3.0 +GCTTT,70.55140100783028,3.0 +GGAAA,95.32530561664231,3.0 +GGAAC,87.67993073378128,3.0 +GGAAG,93.56956480367725,3.0 +GGAAT,96.22871976600298,3.0 +GGACA,75.52664521046947,3.0 +GGACC,73.74304950095393,3.0 +GGACG,69.59876224567604,3.0 +GGACT,80.53463068359464,3.0 +GGAGA,107.9415590822824,3.0 +GGAGC,95.22177447080186,3.0 +GGAGG,101.32625420832143,3.0 +GGAGT,103.90466329961441,3.0 +GGATA,84.67639059495399,3.0 +GGATC,88.1025066779018,3.0 +GGATG,71.21997070586482,3.0 +GGATT,87.35533755031193,3.0 +GGCAA,74.47555616857005,3.0 +GGCAC,68.14532725535699,3.0 +GGCAG,82.61448828843817,3.0 +GGCAT,73.4915092944737,3.0 +GGCCA,67.37502357525973,3.0 +GGCCC,62.832462408406165,3.0 +GGCCG,70.7040790223083,3.0 +GGCCT,64.64368812742207,3.0 +GGCGA,79.03337574608014,3.0 +GGCGC,71.6106409754449,3.0 +GGCGG,79.81471628524028,3.0 +GGCGT,73.8337406503915,3.0 +GGCTA,77.40887146075627,3.0 +GGCTC,73.77456619117491,3.0 +GGCTG,77.3214499581956,3.0 +GGCTT,72.13295010965425,3.0 +GGGAA,106.43480565639285,3.0 +GGGAC,105.89754610136981,3.0 +GGGAG,106.11984677470436,3.0 +GGGAT,111.47780200228001,3.0 +GGGCA,94.72226923485904,3.0 +GGGCC,92.99983735648848,3.0 +GGGCG,94.83975708186055,3.0 +GGGCT,97.95648199454003,3.0 +GGGGA,102.3712725297191,3.0 +GGGGC,93.61704214950731,3.0 +GGGGG,99.8293701614573,3.0 +GGGGT,100.47069813153894,3.0 +GGGTA,105.41874488086319,3.0 +GGGTC,103.90466329961441,3.0 +GGGTG,96.09264359876634,3.0 +GGGTT,103.66513522408387,3.0 +GGTAA,78.64270547650007,3.0 +GGTAC,70.53358468267183,3.0 +GGTAG,86.51037031542216,3.0 +GGTAT,78.37805786885896,3.0 +GGTCA,70.6228824108638,3.0 +GGTCC,65.63576645598825,3.0 +GGTCG,72.75552179698519,3.0 +GGTCT,69.94321506398565,3.0 +GGTGA,85.24669519367367,3.0 +GGTGC,76.0956886847785,3.0 +GGTGG,85.8402624251035,3.0 +GGTGT,79.95186646383416,3.0 +GGTTA,74.41636361297364,3.0 +GGTTC,70.79573762529385,3.0 +GGTTG,75.72365241561562,3.0 +GGTTT,70.91790610971864,3.0 +GTAAA,89.00679103473057,3.0 +GTAAC,83.37907887115125,3.0 +GTAAG,87.76804218278396,3.0 +GTAAT,89.97067977156624,3.0 +GTACA,74.48369514716072,3.0 +GTACC,71.12715191410166,3.0 +GTACG,71.27423050457162,3.0 +GTACT,77.4759036140579,3.0 +GTAGA,104.97897632176223,3.0 +GTAGC,88.51015012189545,3.0 +GTAGG,93.47385059114399,3.0 +GTAGT,97.70878472122146,3.0 +GTATA,82.86468576110555,3.0 +GTATC,83.44634012813518,3.0 +GTATG,75.15850563652866,3.0 +GTATT,83.32206701716882,3.0 +GTCAA,74.83791702378822,3.0 +GTCAC,68.51954826707303,3.0 +GTCAG,80.84303217514953,3.0 +GTCAT,74.96672494698132,3.0 +GTCCA,66.75554355778868,3.0 +GTCCC,61.29389853358644,3.0 +GTCCG,69.19931294787756,3.0 +GTCCT,65.10218310146236,3.0 +GTCGA,79.95836482111713,3.0 +GTCGC,75.09553287017441,3.0 +GTCGG,80.67549307975233,3.0 +GTCGT,77.77353986794013,3.0 +GTCTA,74.1527682307568,3.0 +GTCTC,70.45844675675802,3.0 +GTCTG,75.94852769428783,3.0 +GTCTT,70.51114095507147,3.0 +GTGAA,101.21565967648463,3.0 +GTGAC,102.83561210735999,3.0 +GTGAG,100.63073232441096,3.0 +GTGAT,107.18241436788847,3.0 +GTGCA,90.61065771657522,3.0 +GTGCC,88.10513539668591,3.0 +GTGCG,87.86066314422024,3.0 +GTGCT,93.44871348973263,3.0 +GTGGA,102.12942917173055,3.0 +GTGGC,90.29487044235451,3.0 +GTGGG,96.69543032078539,3.0 +GTGGT,98.95755879189149,3.0 +GTGTA,98.23672604618486,3.0 +GTGTC,97.42706783509146,3.0 +GTGTG,86.20332394403717,3.0 +GTGTT,96.31792999907702,3.0 +GTTAA,78.41928300331013,3.0 +GTTAC,71.25321920206463,3.0 +GTTAG,83.01257376926287,3.0 +GTTAT,78.31573147708076,3.0 +GTTCA,70.11861308331895,3.0 +GTTCC,65.07591892940536,3.0 +GTTCG,72.39056595604043,3.0 +GTTCT,70.210739293245,3.0 +GTTGA,83.8487484475726,3.0 +GTTGC,78.18692355388767,3.0 +GTTGG,85.24069211497996,3.0 +GTTGT,82.13703384268355,3.0 +GTTTA,73.27029686363005,3.0 +GTTTC,70.01643213824416,3.0 +GTTTG,75.2848016767879,3.0 +GTTTT,70.90390759362172,3.0 +TAAAA,93.10199546147226,3.0 +TAAAC,90.99292652221082,3.0 +TAAAG,92.75295547441955,3.0 +TAAAT,97.67277714202571,3.0 +TAACA,78.9810178241838,3.0 +TAACC,76.53765059538999,3.0 +TAACG,74.16870780320494,3.0 +TAACT,82.41918469054305,3.0 +TAAGA,108.06223171643799,3.0 +TAAGC,93.37016026743838,3.0 +TAAGG,97.8335249998958,3.0 +TAAGT,101.13118155888519,3.0 +TAATA,92.20785983555146,3.0 +TAATC,92.61644429129537,3.0 +TAATG,77.08002457387282,3.0 +TAATT,94.139292250866,3.0 +TACAA,73.1442295756671,3.0 +TACAC,68.55852962612039,3.0 +TACAG,77.25137111244409,3.0 +TACAT,73.22913199188059,3.0 +TACCA,67.96438524315963,3.0 +TACCC,62.91316347772191,3.0 +TACCG,71.36010239480268,3.0 +TACCT,65.47111891611495,3.0 +TACGA,74.2817008959855,3.0 +TACGC,69.91773659105995,3.0 +TACGG,75.51734349555255,3.0 +TACGT,71.90265716139062,3.0 +TACTA,77.85404908049125,3.0 +TACTC,73.91529260094448,3.0 +TACTG,78.81726034508848,3.0 +TACTT,72.00131124502497,3.0 +TAGAA,110.11026437767816,3.0 +TAGAC,109.63587902523194,3.0 +TAGAG,109.70880409771132,3.0 +TAGAT,116.29651093829936,3.0 +TAGCA,99.49230593690164,3.0 +TAGCC,96.44610700867878,3.0 +TAGCG,97.89453931629106,3.0 +TAGCT,101.07960845765511,3.0 +TAGGA,103.58049002243392,3.0 +TAGGC,97.42706783509146,3.0 +TAGGG,101.3489342430174,3.0 +TAGGT,101.88758563804898,3.0 +TAGTA,108.45679095074777,3.0 +TAGTC,108.17551645961166,3.0 +TAGTG,97.6741345462047,3.0 +TAGTT,106.51051704428443,3.0 +TATAA,77.4377647029834,3.0 +TATAC,70.74894987451977,3.0 +TATAG,83.85164228710562,3.0 +TATAT,77.60091814951757,3.0 +TATCA,71.44449389574923,3.0 +TATCC,66.20492219812931,3.0 +TATCG,74.19387723212975,3.0 +TATCT,70.1550324869252,3.0 +TATGA,80.72628013046508,3.0 +TATGC,73.0102000562794,3.0 +TATGG,82.25321365524128,3.0 +TATGT,76.03231797376137,3.0 +TATTA,76.34807906972391,3.0 +TATTC,71.74086445720259,3.0 +TATTG,78.66366108432356,3.0 +TATTT,71.36010239480268,3.0 +TCAAA,83.97742618224684,3.0 +TCAAC,79.9716118982193,3.0 +TCAAG,83.7338870353461,3.0 +TCAAT,84.37542970228338,3.0 +TCACA,72.61726791491571,3.0 +TCACC,69.65728962754457,3.0 +TCACG,70.1550324869252,3.0 +TCACT,75.50500650823128,3.0 +TCAGA,104.88071051967535,3.0 +TCAGC,90.40757750691822,3.0 +TCAGG,88.62036726515129,3.0 +TCAGT,93.06287108831647,3.0 +TCATA,79.1031863086086,3.0 +TCATC,78.76032667449456,3.0 +TCATG,74.5309704460317,3.0 +TCATT,79.39888923663196,3.0 +TCCAA,74.68568939833592,3.0 +TCCAC,67.7415286771481,3.0 +TCCAG,82.01359298949559,3.0 +TCCAT,73.95466241723229,3.0 +TCCCA,64.73186019051246,3.0 +TCCCC,60.1292763186002,3.0 +TCCCG,69.29917527232718,3.0 +TCCCT,63.91266114631842,3.0 +TCCGA,79.81471628524028,3.0 +TCCGC,75.84154558638633,3.0 +TCCGG,81.83317928617295,3.0 +TCCGT,78.31573147708076,3.0 +TCCTA,71.86230873504884,3.0 +TCCTC,67.54711013278025,3.0 +TCCTG,76.88338720250827,3.0 +TCCTT,67.92937893841585,3.0 +TCGAA,96.21726598743238,3.0 +TCGAC,96.71882419617243,3.0 +TCGAG,96.72695755258685,3.0 +TCGAT,100.77478381220142,3.0 +TCGCA,85.20177470818749,3.0 +TCGCC,84.22714076268204,3.0 +TCGCG,81.05393915692218,3.0 +TCGCT,88.4213004108736,3.0 +TCGGA,105.08717492805039,3.0 +TCGGC,92.82556781334489,3.0 +TCGGG,94.8851173269455,3.0 +TCGGT,102.19300649772495,3.0 +TCGTA,90.90134127105608,3.0 +TCGTC,90.1606956248493,3.0 +TCGTG,81.74202290165046,3.0 +TCGTT,88.90264847722264,3.0 +TCTAA,78.51979802728408,3.0 +TCTAC,72.52220482066971,3.0 +TCTAG,84.77783802172081,3.0 +TCTAT,77.83838446760048,3.0 +TCTCA,70.08238676632837,3.0 +TCTCC,64.93851979381928,3.0 +TCTCG,73.42355527995966,3.0 +TCTCT,70.08238676632837,3.0 +TCTGA,87.63357492426711,3.0 +TCTGC,78.83096334554614,3.0 +TCTGG,86.88557081565722,3.0 +TCTGT,82.31337814455736,3.0 +TCTTA,74.30491855128562,3.0 +TCTTC,69.83782666735875,3.0 +TCTTG,78.65294785178125,3.0 +TCTTT,70.7634403310388,3.0 +TGAAA,94.61856819387995,3.0 +TGAAC,87.27999933436912,3.0 +TGAAG,93.76338968137142,3.0 +TGAAT,95.83499718179678,3.0 +TGACA,75.41344181315844,3.0 +TGACC,73.47097552556217,3.0 +TGACG,70.26591287128585,3.0 +TGACT,80.40886066820104,3.0 +TGAGA,107.93245150813921,3.0 +TGAGC,95.96242011700865,3.0 +TGAGG,100.8154127324303,3.0 +TGAGT,104.24337429473304,3.0 +TGATA,83.30411189453413,3.0 +TGATC,86.42767091383465,3.0 +TGATG,71.41530584065487,3.0 +TGATT,86.37563538860779,3.0 +TGCAA,75.3159671566861,3.0 +TGCAC,68.35367017337245,3.0 +TGCAG,84.35897095699693,3.0 +TGCAT,74.11337223621082,3.0 +TGCCA,67.56855768102585,3.0 +TGCCC,62.2703024104534,3.0 +TGCCG,71.41530584065487,3.0 +TGCCT,65.01988902364405,3.0 +TGCGA,80.37665859955627,3.0 +TGCGC,72.75552179698519,3.0 +TGCGG,81.50762066962523,3.0 +TGCGT,75.01486518253017,3.0 +TGCTA,76.68935430429275,3.0 +TGCTC,72.6282631347357,3.0 +TGCTG,77.77626855604572,3.0 +TGCTT,71.74086445720259,3.0 +TGGAA,105.58348097001246,3.0 +TGGAC,104.15679787554033,3.0 +TGGAG,104.59255255218295,3.0 +TGGAT,109.93333677501815,3.0 +TGGCA,94.21972433418667,3.0 +TGGCC,92.32652306877577,3.0 +TGGCG,93.94534245544321,3.0 +TGGCT,97.39882395422151,3.0 +TGGGA,102.47359104241907,3.0 +TGGGC,95.2294099132319,3.0 +TGGGG,99.41876664166666,3.0 +TGGGT,100.7766170808923,3.0 +TGGTA,103.86264087029342,3.0 +TGGTC,102.39185496559384,3.0 +TGGTG,96.22871976600298,3.0 +TGGTT,103.04818501023348,3.0 +TGTAA,79.45660193002288,3.0 +TGTAC,71.11124045253479,3.0 +TGTAG,87.34399762081046,3.0 +TGTAT,79.16359922783784,3.0 +TGTCA,70.71606260314421,3.0 +TGTCC,65.15155933732174,3.0 +TGTCG,73.20439793022935,3.0 +TGTCT,70.35480826248974,3.0 +TGTGA,86.2666946550543,3.0 +TGTGC,76.87022006572471,3.0 +TGTGG,86.84678078629545,3.0 +TGTGT,81.4518107315098,3.0 +TGTTA,74.63182420612486,3.0 +TGTTC,70.6323375057766,3.0 +TGTTG,76.25625171360429,3.0 +TGTTT,71.02160715069773,3.0 +TTAAA,88.98389191085376,3.0 +TTAAC,83.99116730808721,3.0 +TTAAG,88.31704523414761,3.0 +TTAAT,90.72092756773344,3.0 +TTACA,74.46078846848457,3.0 +TTACC,71.30828314766855,3.0 +TTACG,71.31100735560244,3.0 +TTACT,78.06600178704689,3.0 +TTAGA,106.71323585342418,3.0 +TTAGC,90.71056888368105,3.0 +TTAGG,91.60367758397071,3.0 +TTAGT,99.91253040580453,3.0 +TTATA,83.37144342872121,3.0 +TTATC,83.20052522970317,3.0 +TTATG,75.37246985233166,3.0 +TTATT,84.112089074928,3.0 +TTCAA,74.94946205479772,3.0 +TTCAC,69.29917527232718,3.0 +TTCAG,81.64327007854568,3.0 +TTCAT,75.2243407933675,3.0 +TTCCA,66.92259791617711,3.0 +TTCCC,61.78357058326347,3.0 +TTCCG,69.713299679996,3.0 +TTCCT,65.49014031961454,3.0 +TTCGA,80.96040529461689,3.0 +TTCGC,76.08444538643136,3.0 +TTCGG,81.4243850529776,3.0 +TTCGT,78.68068731802552,3.0 +TTCTA,74.33268964227133,3.0 +TTCTC,70.70565402227757,3.0 +TTCTG,76.2959133559782,3.0 +TTCTT,70.74601808529705,3.0 +TTGAA,102.01925595172668,3.0 +TTGAC,103.65252854799547,3.0 +TTGAG,101.09349119205983,3.0 +TTGAT,107.15439941500782,3.0 +TTGCA,91.80219750591482,3.0 +TTGCC,89.68871550968113,3.0 +TTGCG,87.8013430355001,3.0 +TTGCT,94.2695157326237,3.0 +TTGGA,103.4332808920591,3.0 +TTGGC,92.35579756483001,3.0 +TTGGG,96.3856463009514,3.0 +TTGGT,99.95164037213367,3.0 +TTGTA,98.23163235450001,3.0 +TTGTC,97.74634085858878,3.0 +TTGTG,86.57601923767265,3.0 +TTGTT,96.20675691016524,3.0 +TTTAA,77.1335608687731,3.0 +TTTAC,72.14442795816693,3.0 +TTTAG,83.01257376926287,3.0 +TTTAT,77.81660064293777,3.0 +TTTCA,70.68171391732899,3.0 +TTTCC,65.73741856510927,3.0 +TTTCG,73.39636415159302,3.0 +TTTCT,70.87501716248273,3.0 +TTTGA,85.61040677673654,3.0 +TTTGC,78.05226066120652,3.0 +TTTGG,85.35949238875044,3.0 +TTTGT,81.02069839121474,3.0 +TTTTA,73.75864190401846,3.0 +TTTTC,70.37074783493787,3.0 +TTTTG,76.25344853166334,3.0 +TTTTT,71.44850004771521,3.0 diff --git a/xpore/diffmod/configurator.py b/xpore/diffmod/configurator.py index 5dc5cf9..e18f494 100644 --- a/xpore/diffmod/configurator.py +++ b/xpore/diffmod/configurator.py @@ -19,7 +19,7 @@ def get_paths(self): if 'prior' in self.yaml: paths['model_kmer'] = os.path.abspath(self.yaml['prior']) else: - paths['model_kmer'] = os.path.join(os.path.dirname(__file__),'model_kmer.csv') + paths['model_kmer'] = os.path.join(os.path.dirname(__file__),'RNA004_5mer_model.csv') paths['out_dir'] = os.path.join(os.path.abspath(self.yaml['out'])) paths.update(misc.makedirs(paths['out_dir'],sub_dirs=['models'])) diff --git a/xpore/diffmod/io.py b/xpore/diffmod/io.py index 8558233..c471a4d 100644 --- a/xpore/diffmod/io.py +++ b/xpore/diffmod/io.py @@ -43,19 +43,58 @@ def load_data(idx, data_dict, min_count, max_count, pooling=False): y, read_ids, condition_labels, run_labels = [], [], [], [] n_reads = defaultdict(list) - for (condition_name,run_name), d_dict in data_dict.items(): - if d_dict is not None: - norm_means = d_dict[idx][pos][kmer]#['norm_means'] - n_reads_per_run = len(norm_means) - # In case of pooling==False, if not enough reads, don't include them. - if (not pooling) and ((n_reads_per_run < min_count) or (n_reads_per_run > max_count)): - continue - # - n_reads[condition_name] += [n_reads_per_run] - y += norm_means - # read_ids += list(data_dict[run_name][idx][pos][kmer]['read_ids'][:]) - condition_labels += [condition_name]*n_reads_per_run - run_labels += [run_name]*n_reads_per_run + if not pooling: + for (condition_name,run_name), d_dict in data_dict.items(): + if d_dict is not None: + norm_means = d_dict[idx][pos][kmer]#['norm_means'] + n_reads_per_run = len(norm_means) + # Drop runs with too few reads (can't be modelled). For too many reads, + # truncate to the first max_count (file order, deterministic) rather than + # dropping the run, so a user can further subset at the diffmod stage on + # top of any cap already applied during dataprep. + if n_reads_per_run < min_count: + continue + if (max_count is not None) and (n_reads_per_run > max_count): + norm_means = norm_means[:max_count] + n_reads_per_run = max_count + n_reads[condition_name] += [n_reads_per_run] + y += norm_means + # read_ids += list(data_dict[run_name][idx][pos][kmer]['read_ids'][:]) + condition_labels += [condition_name]*n_reads_per_run + run_labels += [run_name]*n_reads_per_run + else: + # Pooling: reads from all runs of a condition are combined into one group. + # Gather the runs present at this (pos,kmer), preserving data_dict order. + runs_by_condition = defaultdict(list) + for (condition_name,run_name), d_dict in data_dict.items(): + if d_dict is not None: + runs_by_condition[condition_name] += [(run_name, d_dict[idx][pos][kmer])] + for condition_name, runs in runs_by_condition.items(): + total = sum(len(norm_means) for _, norm_means in runs) + if (max_count is None) or (total <= max_count): + # Under the cap: keep every read (per-run order). Read order does not + # affect the model output, so this matches the original pooled behaviour. + for run_name, norm_means in runs: + n_reads[condition_name] += [len(norm_means)] + y += norm_means + condition_labels += [condition_name]*len(norm_means) + run_labels += [run_name]*len(norm_means) + else: + # Over the cap: select reads round-robin across the runs (1st read of + # each run, then 2nd, ...) up to max_count, so the subset is evenly + # distributed across runs and deterministic, instead of front-loading + # the first run. + run_streams = [[(value, run_name) for value in norm_means] + for run_name, norm_means in runs] + selected = [read for group in itertools.zip_longest(*run_streams) + for read in group if read is not None][:max_count] + per_run = defaultdict(int) + for value, run_name in selected: + y += [value] + condition_labels += [condition_name] + run_labels += [run_name] + per_run[run_name] += 1 + n_reads[condition_name] += list(per_run.values()) y = np.array(y) # read_ids = np.array(read_ids) diff --git a/xpore/scripts/dataprep.py b/xpore/scripts/dataprep.py index 5b2d803..3be699d 100644 --- a/xpore/scripts/dataprep.py +++ b/xpore/scripts/dataprep.py @@ -108,18 +108,39 @@ def t2g(gene_id,fasta_dict,annotation_dict,g2t_mapping,df_eventalign_index,readc return n_reads, tx_ids, t2g_dict -def combine(events_str): +def combine(events_str,kmer_source='reference_kmer'): f_string = StringIO(events_str) eventalign_result = pd.read_csv(f_string,delimiter='\t',names=['contig','position','reference_kmer','read_index', 'strand','event_index','event_level_mean','event_stdv','event_length','model_kmer', 'model_mean', 'model_stdv', 'standardized_level', 'start_idx', 'end_idx']) f_string.close() - cond_successfully_eventaligned = eventalign_result['reference_kmer'] == eventalign_result['model_kmer'] + + # Reverse complement helper for matching genome-aligned reads where + # reference_kmer may be the reverse complement of model_kmer + def _revcomp(seq): + comp = str.maketrans('ACGT','TGCA') + return seq.translate(comp)[::-1] + + # When using model_kmer, accept both direct and reverse-complement matches + # (genome alignments can have reverse-oriented reads). + # When using reference_kmer, only accept exact matches (same as xPore v2.1 behaviour) + if kmer_source == 'model_kmer': + cond_successfully_eventaligned = ( + (eventalign_result['reference_kmer'] == eventalign_result['model_kmer']) | + (eventalign_result['reference_kmer'].apply(_revcomp) == eventalign_result['model_kmer']) + ) + else: + cond_successfully_eventaligned = ( + eventalign_result['reference_kmer'] == eventalign_result['model_kmer'] + ) if cond_successfully_eventaligned.sum() != 0: eventalign_result = eventalign_result[cond_successfully_eventaligned] - keys = ['read_index','contig','position','reference_kmer'] # for groupby + # Select kmer column based on kmer_source flag: model_kmer for genome alignments, + # reference_kmer for transcriptome alignments + kmer_col = 'model_kmer' if kmer_source == 'model_kmer' else 'reference_kmer' + keys = ['read_index','contig','position',kmer_col] # for groupby eventalign_result['length'] = pd.to_numeric(eventalign_result['end_idx'])-pd.to_numeric(eventalign_result['start_idx']) eventalign_result['sum_norm_mean'] = pd.to_numeric(eventalign_result['event_level_mean']) * eventalign_result['length'] @@ -151,11 +172,12 @@ def combine(events_str): # df_events = eventalign_result[['read_index']+features] # # print(df_events.head()) - features = ['transcript_id','transcriptomic_position','reference_kmer','norm_mean'] + # Use the selected kmer column in the output features + features = ['transcript_id','transcriptomic_position',kmer_col,'norm_mean'] # np_events = eventalign_result[features].reset_index().values.ravel().view(dtype=[('transcript_id', 'S15'), ('transcriptomic_position', ' 1): - data_dict[read_index] = data - readcount += 1 - if readcount > readcount_max: - break - - if readcount > readcount_max: - break + if (result is not None): + data, kmer_col = result + if data.size > 1: + data_dict[read_index] = data + readcount += 1 + if len(data_dict)>=readcount_min: # print(gene_id,len(data_dict)) #len(data_dict) is the number of reads to be processed. - task_queue.put((gene_id,data_dict,t2g_mapping,out_paths)) # Blocked if necessary until a free slot is available. + task_queue.put((gene_id,data_dict,kmer_col,t2g_mapping,readcount_max,out_paths)) # Blocked if necessary until a free slot is available. gene_ids_processed += [gene_id] @@ -375,7 +400,7 @@ def parallel_preprocess_gene(eventalign_filepath,fasta_dict,annotation_dict,is_g f.write('Total %d genes.\n' %len(gene_ids_processed)) f.write(helper.decor_message('successfully finished')) -def preprocess_gene(gene_id,data_dict,t2g_mapping,out_paths,locks): +def preprocess_gene(gene_id,data_dict,kmer_col,t2g_mapping,readcount_max,out_paths,locks): """ Convert transcriptomic to genomic coordinates for a gene. @@ -423,7 +448,8 @@ def preprocess_gene(gene_id,data_dict,t2g_mapping,out_paths,locks): # Based on Ensembl, remove transcript version. events_per_read['transcript_id'] = tx_ids - events_per_read = np.array(events_per_read,dtype=np.dtype([('transcript_id', 'S15'), ('transcriptomic_position', ' readcount_max: + y_array = y_array[:readcount_max] + g_kmer_array = g_kmer_array[:readcount_max] + g_positions_array = g_positions_array[:readcount_max] if (len(set(g_kmer_array)) == 1) and ('XXXXX' in set(g_kmer_array)) or (len(y_array) == 0): continue - + if 'XXXXX' in set(g_kmer_array): - y_array = y_array[g_kmer_array != 'XXXXX'] + y_array = y_array[g_kmer_array != 'XXXXX'] assert len(y_array) == len(g_kmer_array) - (g_kmer_array=='XXXXX').sum() - g_kmer_array = g_kmer_array[g_kmer_array != 'XXXXX'] - + g_kmer_array = g_kmer_array[g_kmer_array != 'XXXXX'] + try: assert len(set(g_kmer_array)) == 1 assert list(set(g_kmer_array))[0].count('N') == 0 ##to weed out the mapped kmers from tx_seq that contain 'N', which is not in diffmod's model_kmer assert {position} == set(g_positions_array) except: asserted = False - break + continue # skip this position, don't stop processing remaining positions kmer = set(g_kmer_array).pop() data[position] = {kmer: list(y_array)} #,'read_ids': [read_id.decode('UTF-8') for read_id in read_id_array]} @@ -502,7 +533,7 @@ def preprocess_gene(gene_id,data_dict,t2g_mapping,out_paths,locks): f.write(log_str + '\n') -def parallel_preprocess_tx(eventalign_filepath,out_dir,n_processes,readcount_min,readcount_max,resume): +def parallel_preprocess_tx(eventalign_filepath,out_dir,n_processes,readcount_min,readcount_max,resume,kmer_source='reference_kmer'): # Create output paths and locks. out_paths,locks = dict(),dict() @@ -539,6 +570,9 @@ def parallel_preprocess_tx(eventalign_filepath,out_dir,n_processes,readcount_min tx_ids = df_eventalign_index['transcript_id'].values.tolist() tx_ids = list(dict.fromkeys(tx_ids)) df_eventalign_index.set_index('transcript_id',inplace=True) + # kmer_col == kmer_source (argparse restricts it to a valid column name); bind it up front so it + # is always defined even if no read in this transcript returns a result from combine() (avoids NameError). + kmer_col = kmer_source with open(eventalign_filepath,'r') as eventalign_result: for tx_id in tx_ids: data_dict = dict() @@ -547,14 +581,15 @@ def parallel_preprocess_tx(eventalign_filepath,out_dir,n_processes,readcount_min read_index,pos_start,pos_end = row['read_index'],row['pos_start'],row['pos_end'] eventalign_result.seek(pos_start,0) events_str = eventalign_result.read(pos_end-pos_start) - data = combine(events_str) - if (data is not None) and (data.size > 1): - data_dict[read_index] = data - readcount += 1 - if readcount > readcount_max: - break + # Pass kmer_source to combine so it selects the correct kmer column + result = combine(events_str,kmer_source) + if (result is not None): + data, kmer_col = result + if data.size > 1: + data_dict[read_index] = data + readcount += 1 if readcount>=readcount_min: - task_queue.put((tx_id,data_dict,out_paths)) # Blocked if necessary until a free slot is available. + task_queue.put((tx_id,data_dict,kmer_col,readcount_max,out_paths)) # Blocked if necessary until a free slot is available. tx_ids_processed += [tx_id] # Put the stop task into task_queue. @@ -567,23 +602,25 @@ def parallel_preprocess_tx(eventalign_filepath,out_dir,n_processes,readcount_min f.write('Total %d transcripts.\n' %len(tx_ids_processed)) f.write(helper.decor_message('successfully finished')) -def preprocess_tx(tx_id,data_dict,out_paths,locks): +def preprocess_tx(tx_id,data_dict,kmer_col,readcount_max,out_paths,locks): """ - Convert transcriptomic to genomic coordinates for a gene. - + Reshape one aligned reference's per-read events into per-position signal arrays. + + Transcriptome mode (no --genome) (or genome alignment mode): positions stay in original reference coordinates; + no conversion to genomic coordinates for tx alignments is done here (that's preprocess_gene). Groups reads at each + position, caps them at readcount_max, checks k-mer consistency, and appends the + result to the data.json / .index / .readcount / .log files. Returns None. + Parameters ---------- - tx_id: str - Transcript ID. - data_dict: {read_id:events_array} - Events for each read. - features: [str] # todo - A list of features to collect from the reads that are aligned to each genomic coordinate in the output. - Returns - ------- - dict - A dict of all specified features collected for each genomic coordinate. + tx_id : str Reference (for tx alignments, Transcript ID). + data_dict : {read_index: events_array} Per-read arrays from combine(). + kmer_col : str 'reference_kmer' or 'model_kmer'. + readcount_max : int or None Per-site read cap; None means no cap. + out_paths, locks : dict Output file paths and their write locks. + """ + # features = ['read_id','transcript_id','transcriptomic_position','reference_kmer','norm_mean','start_idx','end_idx'] # columns in the eventalign file per read. @@ -603,39 +640,53 @@ def preprocess_tx(tx_id,data_dict,out_paths,locks): events = np.concatenate(events) - # Sort and split - idx_sorted = np.argsort(events['transcriptomic_position']) - unique_positions, index = np.unique(events['transcriptomic_position'][idx_sorted],return_index = True) + # Sort and split + if kmer_col == 'model_kmer': + # Group by (position, kmer) so forward/reverse reads at the same position + # produce separate entries with their respective kmers + idx_sorted = np.lexsort((events[kmer_col], events['transcriptomic_position'])) + positions_sorted = events['transcriptomic_position'][idx_sorted] + kmers_sorted = events[kmer_col][idx_sorted] + compound_keys = np.array(list(zip(positions_sorted, kmers_sorted)), + dtype=[('pos', positions_sorted.dtype), ('kmer', kmers_sorted.dtype)]) + unique_pairs, index = np.unique(compound_keys, return_index=True) + unique_positions = unique_pairs['pos'] + else: + idx_sorted = np.argsort(events['transcriptomic_position']) + unique_positions, index = np.unique(events['transcriptomic_position'][idx_sorted], return_index=True) y_arrays = np.split(events['norm_mean'][idx_sorted], index[1:]) -# read_id_arrays = np.split(events['read_id'][idx_sorted], index[1:]) - reference_kmer_arrays = np.split(events['reference_kmer'][idx_sorted], index[1:]) + kmer_arrays = np.split(events[kmer_col][idx_sorted], index[1:]) # Prepare # print('Reformating the data for each genomic position ...') data = defaultdict(dict) # for each position, make it ready for json dump asserted = True -# for key_tuple,y_array,reference_kmer_array in zip(key_tuples,y_arrays,reference_kmer_arrays): - for position,y_array,reference_kmer_array in zip(unique_positions,y_arrays,reference_kmer_arrays): - +# for key_tuple,y_array,kmer_array in zip(key_tuples,y_arrays,kmer_arrays): + for position,y_array,kmer_array in zip(unique_positions,y_arrays,kmer_arrays): + position = int(position) - if (len(set(reference_kmer_array)) == 1) and ('XXXXX' in set(reference_kmer_array)) or (len(y_array) == 0): + # Cap reads per site (not per transcript) to readcount_max; None means no limit + if readcount_max is not None and len(y_array) > readcount_max: + y_array = y_array[:readcount_max] + kmer_array = kmer_array[:readcount_max] + if (len(set(kmer_array)) == 1) and ('XXXXX' in set(kmer_array)) or (len(y_array) == 0): continue - - if 'XXXXX' in set(reference_kmer_array): - y_array = y_array[reference_kmer_array != 'XXXXX'] - assert len(y_array) == len(reference_kmer_array) - (reference_kmer_array=='XXXXX').sum() - reference_kmer_array = reference_kmer_array[reference_kmer_array != 'XXXXX'] - + + if 'XXXXX' in set(kmer_array): + y_array = y_array[kmer_array != 'XXXXX'] + assert len(y_array) == len(kmer_array) - (kmer_array=='XXXXX').sum() + kmer_array = kmer_array[kmer_array != 'XXXXX'] + try: - assert len(set(reference_kmer_array)) == 1 - assert list(set(reference_kmer_array))[0].count('N') == 0 ##to weed out the mapped kmers from tx_seq that contain 'N', which is not in diffmod's model_kmer + assert len(set(kmer_array)) == 1 + assert list(set(kmer_array))[0].count('N') == 0 ##to weed out the mapped kmers from tx_seq that contain 'N', which is not in diffmod's model_kmer except: asserted = False - break - kmer = set(reference_kmer_array).pop() + continue # skip this position, don't stop processing remaining positions + kmer = set(kmer_array).pop() - data[position] = {kmer: list(np.around(y_array,decimals=2))} + data[position][kmer] = list(np.around(y_array,decimals=2)) # write to file. log_str = '%s: %s.' %(tx_id,asserted) @@ -728,6 +779,20 @@ def dataprep(args): readcount_max = args.readcount_max resume = args.resume genome = args.genome + kmer_source = args.kmer_source + + # Check for likely genome alignment when user has selected reference_kmer + if kmer_source == 'reference_kmer': + with open(eventalign_filepath,'r') as f: + first_line = f.readline() + if first_line: + contig = first_line.split('\t')[0] + if contig.startswith('chr'): + raise ValueError( + "It looks like you are using a genome alignment (reference name '%s' starts with 'chr'), " + "but --kmer_source is set to 'reference_kmer'. Genome alignments can contain reverse-oriented " + "reads, using reference_kmer will cause all reverse reads to be filtered. Please use '--kmer_source model_kmer' instead." % contig + ) if genome and (None in [args.gtf_or_gff,args.transcript_fasta]): print('please provide the following') @@ -750,6 +815,6 @@ def dataprep(args): # gtf_or_gff = mergeGTFtxIDversion(gtf_or_gff,out_dir) annotation_dict,is_gff = readAnnotation(gtf_or_gff) fasta_dict = readFasta(transcript_fasta,is_gff) - parallel_preprocess_gene(eventalign_filepath,fasta_dict,annotation_dict,is_gff,out_dir,n_processes,readcount_min,readcount_max,resume) + parallel_preprocess_gene(eventalign_filepath,fasta_dict,annotation_dict,is_gff,out_dir,n_processes,readcount_min,readcount_max,resume,kmer_source) else: - parallel_preprocess_tx(eventalign_filepath,out_dir,n_processes,readcount_min,readcount_max,resume) + parallel_preprocess_tx(eventalign_filepath,out_dir,n_processes,readcount_min,readcount_max,resume,kmer_source) diff --git a/xpore/scripts/xpore.py b/xpore/scripts/xpore.py index 32ef87c..640e037 100644 --- a/xpore/scripts/xpore.py +++ b/xpore/scripts/xpore.py @@ -28,11 +28,12 @@ def parse_options(argv): # Optional arguments optional_dataprep.add_argument('--skip_eventalign_indexing', dest='skip_eventalign_indexing', help='skip indexing the eventalign nanopolish output.',default=False,action='store_true') # parser.add_argument('--features', dest='features', help='Signal features to extract.',type=list,default=['norm_mean']) - optional_dataprep.add_argument('--genome', dest='genome', help='to run on Genomic coordinates. Without this argument, the program will run on transcriptomic coordinates',default=False,action='store_true') + optional_dataprep.add_argument('--genome', dest='genome', help='to run on Genomic coordinates. Without this argument, the program will run on transcriptomic coordinates',default=False,action='store_true') + optional_dataprep.add_argument('--kmer_source', dest='kmer_source', help='kmer column to use from eventalign: reference_kmer (default, for transcriptome alignments) or model_kmer (for genome alignments).',type=str, choices=['reference_kmer','model_kmer'], default='reference_kmer') optional_dataprep.add_argument('--n_processes', dest='n_processes', help='number of processes to run.',type=int, default=1) optional_dataprep.add_argument('--chunk_size', dest='chunk_size', help='number of lines from nanopolish eventalign.txt for processing.',type=int, default=1000000) optional_dataprep.add_argument('--readcount_min', dest='readcount_min', help='minimum read counts per gene.',type=int, default=1) - optional_dataprep.add_argument('--readcount_max', dest='readcount_max', help='maximum read counts per gene.',type=int, default=1000) + optional_dataprep.add_argument('--readcount_max', dest='readcount_max', help='maximum read counts per site. Use None for no limit.',type=lambda x: None if x == 'None' else int(x), default=1000) optional_dataprep.add_argument('--resume', dest='resume', help='with this argument, the program will resume from the previous run.',default=False,action='store_true') #todo parser_dataprep._action_groups.append(optional_dataprep) parser_dataprep.set_defaults(func=dataprep) @@ -46,7 +47,7 @@ def parse_options(argv): # Optional arguments optional_diffmod.add_argument('--n_processes', dest='n_processes', help='number of processes to run.',type=int,default=1) optional_diffmod.add_argument('--save_models', dest='save_models', help='with this argument, the program will save the model parameters for each id.',default=False,action='store_true') # todo - optional_diffmod.add_argument('--resume', dest='resume', help='with this argument, the program will resume from the previous run.',default=False,action='store_true') + optional_diffmod.add_argument('--resume', dest='resume', help='with this argument, the program will resume from the previous run.',default=False,action='store_true') optional_diffmod.add_argument('--ids', dest='ids', help='gene or transcript ids to model.',default=[],nargs='*') parser_diffmod._action_groups.append(optional_diffmod) parser_diffmod.set_defaults(func=diffmod)