Rename tpm to cpm throughout the API (breaking) (2.1.0)

CHANGELOG.md

@@ -6,6 +6,10 @@
 * planned new feature: during import of long reads, (optionally) correct for short exon alignment issues.
 * separate new read import and classification of isoforms.
 
+## [2.1.0]
+
+* **breaking**: renamed `tpm` to `cpm` throughout the API (`transcript_table(cpm=..., cpm_pseudocount=...)`, `Gene.cpm()`, `estimate_cpm_threshold`, `_cpm`/`_sum_cpm` columns). The values were always counts per million (no transcript-length normalization) despite the "tpm" name -- correct for full-length long reads, where read count is already a direct proxy for molecule count, but mislabeled. No deprecated alias; update `tpm=` to `cpm=` in existing code.
+
 ## [2.0.7]
 
 * fixed: `add_sample_from_csv` crashed with a confusing `AttributeError` when one transcript_id/gene_id from the coverage csv wasn't found in the transcripts file; the warning now also suggests `infer_genes=True` when the file has no gene annotations (#25)

VERSION.txt

@@ -1 +1 @@
-2.0.7
+2.1.0

docs/notebooks/03_transcriptome_reconstruction.ipynb

@@ -659,7 +659,7 @@
    ],
    "source": [
     "# export transcript table with the same filter criteria:\n",
-    "transcript_tab=isoseq.transcript_table( groups=isoseq.groups(),tpm=True,coverage=True,\n",
+    "transcript_tab=isoseq.transcript_table( groups=isoseq.groups(),cpm=True,coverage=True,\n",
     "                                       min_coverage=5, progress_bar=True,\n",
     "                                       query=query_string)\n",
     "# write to csv file\n",

docs/notebooks/03b_transcriptome_import.ipynb

@@ -22,7 +22,7 @@
     "* corresponding table with **number of long reads per transcripts** for the samples to be added. This is a csv file, with one column for each sample. The sample names are specified in the header (first line). Row names (first column) must correspond to transcript names from the gtf file.\n",
     "\n",
     "For demonstration, we use the gtf file and transcript table exported from isotools in the [previous tutorial](03_transcriptome_reconstruction.html). \n",
-    "Note that in this table we exported not only the read counts, but also TPM values and further information on the transcripts. Therefore, we specify the suffix of the columns with the read counts.\n",
+    "Note that in this table we exported not only the read counts, but also CPM values and further information on the transcripts. Therefore, we specify the suffix of the columns with the read counts.\n",
     "Remember that we exported the sum coverage, so instead of individual replicates, we are adding the pooled samples. \n",
     "Also we filtered the transcripts, hence only a subset of transcripts will be imported. \n",
     "All files are assumed to be stored in a subfolder called 'demonstration_dataset'.\n",

docs/notebooks/04_saturation_analysis.ipynb

@@ -10,7 +10,7 @@
     "* Saturation analysis ([plot_saturation](../isotoolsAPI.html?highlight=plot_saturation#isotools.plots.plot_saturation)):\n",
     "    * Models the propability of observing a transcript with at least n reads. \n",
     "    * Depends on the expression level of the transcript\n",
-    "    * E.g.: \"With the given seq depth, the probability of observing a 1 TPM transcript is > 80%\"    \n",
+    "    * E.g.: \"With the given seq depth, the probability of observing a 1 CPM transcript is > 80%\"    \n",
     "* Rarefaction analysis ([plot_rarefaction](../isotoolsAPI.html?highlight=plot_rarefaction#isotools.plots.plot_rarefaction)):\n",
     "    * Subsamples the reads and counts the number of observed transcripts.\n",
     "    * Helps estimating, how much more transcripts would be observed with deeper sequencing.\n",
@@ -73,7 +73,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "The left figure shows the probability of observing 5 reads of a hypothetical transcript expressed at different levels, from 100 TPM to 5000 TPM (note that subsampling reads from one chromosome only scales the expression level), depending on the sequencing depth. With 20k reads, the probability of observing a 200 TPM transcript would be about 40% (orange line).\n",
+    "The left figure shows the probability of observing 5 reads of a hypothetical transcript expressed at different levels, from 100 CPM to 5000 CPM (note that subsampling reads from one chromosome only scales the expression level), depending on the sequencing depth. With 20k reads, the probability of observing a 200 CPM transcript would be about 40% (orange line).\n",
     "\n",
     "The right figure shows the number of different full splice match (FSM) transcripts, depending on the number of reads. Generaly, the deeper the sequencing, the less new transcripts are discovered, indicating beginning saturation. Lower curves, such as with GM12878_c (green) indicate lower library complexity (e.g. less different templates in the sample)."
    ]

docs/notebooks/06_filtering.ipynb

@@ -427,7 +427,7 @@
     "isoseq.write_gtf(f'{path}/demonstration_dataset_transcripts_balanced.gtf', **balanced)\n",
     "\n",
     "transcript_tab=isoseq.transcript_table(groups=isoseq.groups(),\n",
-    "                                       tpm=True,\n",
+    "                                       cpm=True,\n",
     "                                       coverage=True,\n",
     "                                       progress_bar=True,\n",
     "                                       ** balanced)\n",

docs/notebooks/ENCODE_S2.ipynb

@@ -382,11 +382,11 @@
     "To estimate the saturation of discovered isoforms, IsoTools offers two complimentary figures: \n",
     "* The saturation plot models the expected probability of discovering a transcript \n",
     "    * at given level of confidence (number of supporting long reads) \n",
-    "    * which is present at a given abundance level in the sample (in transcripts per milion, TPM)\n",
+    "    * which is present at a given abundance level in the sample (in counts per million, CPM)\n",
     "    * depending on the total number of long reads\n",
     "* The rarefaction analysis depicts the number of different discovered transcripts when subsampling the long reads\n",
     "\n",
-    "The saturation analysis reveals that at the given sequencing depth, transcripts expressed at 2 TPM have > 80% probability of beeing covered by at least two reads in the individual samples, while transcripts expressed at 1 TPM have 40% to 75% probability. The slope of the rarefaction curve at the right end reflects the number of novel transcripts that could be discovered if sequencing depth would be increased. "
+    "The saturation analysis reveals that at the given sequencing depth, transcripts expressed at 2 CPM have > 80% probability of beeing covered by at least two reads in the individual samples, while transcripts expressed at 1 CPM have 40% to 75% probability. The slope of the rarefaction curve at the right end reflects the number of novel transcripts that could be discovered if sequencing depth would be increased. "
    ]
   },
   {

docs/notebooks/isotools_alzheimer.ipynb

@@ -139,11 +139,11 @@
     "To estimate the saturation of discovered isoforms, IsoTools offers two complimentary figures: \n",
     "* The saturation plot models the expected probability of discovering a transcript \n",
     "    * at given level of confidence (number of supporting long reads) \n",
-    "    * which is present at a given abundance level in the sample (in transcripts per milion, TPM)\n",
+    "    * which is present at a given abundance level in the sample (in counts per million, CPM)\n",
     "    * depending on the total number of long reads\n",
     "* The rarefaction analysis depicts the number of different discovered transcripts when subsampling the long reads\n",
     "\n",
-    "The saturation analysis reveals that at the given sequencing depth, transcripts expressed at 1 TPM have about 90% probability of beeing covered by at least two reads in the individual samples, while transcripts expressed at 0.5 TPM have 60% probability. The slope of the rarefaction curve at the right end reflects the number of novel transcripts that could be discovered if sequencing depth would be increased. "
+    "The saturation analysis reveals that at the given sequencing depth, transcripts expressed at 1 CPM have about 90% probability of beeing covered by at least two reads in the individual samples, while transcripts expressed at 0.5 CPM have 60% probability. The slope of the rarefaction curve at the right end reflects the number of novel transcripts that could be discovered if sequencing depth would be increased. "
    ]
   },
   {

src/isotools/__init__.py

@@ -30,7 +30,7 @@
 from .gene import Gene
 from .transcriptome import Transcriptome
 from .splice_graph import SegmentGraph, SegGraphNode
-from ._transcriptome_stats import estimate_tpm_threshold
+from ._transcriptome_stats import estimate_cpm_threshold
 
 from ._transcriptome_filter import (
     DEFAULT_GENE_FILTER,
@@ -44,7 +44,7 @@
     "Gene",
     "SegmentGraph",
     "SegGraphNode",
-    "estimate_tpm_threshold",
+    "estimate_cpm_threshold",
     "DEFAULT_GENE_FILTER",
     "DEFAULT_TRANSCRIPT_FILTER",
     "DEFAULT_REF_TRANSCRIPT_FILTER",

src/isotools/_transcriptome_io.py

@@ -2327,8 +2327,8 @@ def transcript_table(
     samples=None,
     groups=None,
     coverage=False,
-    tpm=False,
-    tpm_pseudocount=0,
+    cpm=False,
+    cpm_pseudocount=0,
     extra_columns=None,
     **filter_args,
 ):
@@ -2339,8 +2339,8 @@ def transcript_table(
     :param samples: provide a list of samples for which coverage / expression information is added.
     :param groups: provide groups as a dict (as from Transcriptome.groups()), for which coverage / expression information is added.
     :param coverage: If set, coverage information is added for specified samples / groups.
-    :param tpm: If set, expression information (in tpm) is added for specified samples / groups.
-    :param tpm_pseudocount: This value is added to the coverage for each transcript, before calculating tpm.
+    :param cpm: If set, expression information (in cpm) is added for specified samples / groups.
+    :param cpm_pseudocount: This value is added to the coverage for each transcript, before calculating cpm.
     :param extra_columns: Specify the additional information added to the table.
         These can be any transcript property as defined by the key in the transcript dict.
     :param filter_args: Parameters (e.g. "region", "query", "min_coverage",...) are passed to Transcriptome.iter_transcripts.
@@ -2353,7 +2353,7 @@ def transcript_table(
             samples = []
     if groups is None:
         groups = {}
-    if coverage is False and tpm is False:
+    if coverage is False and cpm is False:
         samples = []
         groups = {}
     if extra_columns is None:
@@ -2448,13 +2448,13 @@ def transcript_table(
         if samples:
             if coverage:
                 df_list.append(cov[samples].add_suffix("_coverage"))
-            if tpm:
+            if cpm:
                 total = (
                     stab.loc[samples, "nonchimeric_reads"]
-                    + tpm_pseudocount * cov.shape[0]
+                    + cpm_pseudocount * cov.shape[0]
                 )
                 df_list.append(
-                    ((cov[samples] + tpm_pseudocount) / total * 1e6).add_suffix("_tpm")
+                    ((cov[samples] + cpm_pseudocount) / total * 1e6).add_suffix("_cpm")
                 )
         if groups:
             cov_gr = pd.DataFrame(
@@ -2465,14 +2465,14 @@ def transcript_table(
             )
             if coverage:
                 df_list.append(cov_gr.add_suffix("_sum_coverage"))
-            if tpm:
+            if cpm:
                 total = {
                     group_name: stab.loc[sample, "nonchimeric_reads"].sum()
-                    + tpm_pseudocount * cov.shape[0]
+                    + cpm_pseudocount * cov.shape[0]
                     for group_name, sample in groups.items()
                 }
                 df_list.append(
-                    ((cov_gr + tpm_pseudocount) / total * 1e6).add_suffix("_sum_tpm")
+                    ((cov_gr + cpm_pseudocount) / total * 1e6).add_suffix("_sum_cpm")
                 )
         df = pd.concat(df_list, axis=1)
 

src/isotools/_transcriptome_stats.py

@@ -551,21 +551,21 @@ def alternative_splicing_events(
 # function to optimize (inverse nbinom cdf)
 
 
-def _tpm_fun(tpm_th, n_reads, cov_th=2, p=0.8):
-    return (p - nbinom.cdf(n_reads - cov_th, n=cov_th, p=tpm_th * 1e-6)) ** 2
+def _cpm_fun(cpm_th, n_reads, cov_th=2, p=0.8):
+    return (p - nbinom.cdf(n_reads - cov_th, n=cov_th, p=cpm_th * 1e-6)) ** 2
 
 
-def estimate_tpm_threshold(n_reads, cov_th=2, p=0.8):
+def estimate_cpm_threshold(n_reads, cov_th=2, p=0.8):
     """Estimate the minimum expression level of observable transcripts at given coverage.
 
-    The function returns the expression level in transcripts per million (TPM), that can be observed
+    The function returns the expression level in counts per million (CPM), that can be observed
     at the given sequencing depth.
 
     :param n_reads: The sequencing depth (total number of reads) for the sample.
     :param cov_th: The requested minimum number of reads per transcripts.
     :param p: The probability of a transcript at threshold expression level to be observed.
     """
-    return minimize_scalar(_tpm_fun, bounds=(0.01, 1000), args=(n_reads, cov_th, p))[
+    return minimize_scalar(_cpm_fun, bounds=(0.01, 1000), args=(n_reads, cov_th, p))[
         "x"
     ]
 

src/isotools/gene.py

@@ -856,13 +856,13 @@ def _get_info(self, transcript_id, key, sample_i, group_i, **kwargs):
             )  # only the names of the subcategories
         elif key == "coverage":
             return self.coverage[sample_i, transcript_id]
-        elif key == "tpm":
-            return self.tpm(kwargs.get("pseudocount", 1))[sample_i, transcript_id]
+        elif key == "cpm":
+            return self.cpm(kwargs.get("pseudocount", 1))[sample_i, transcript_id]
         elif key == "group_coverage_sum":
             return tuple(self.coverage[si, transcript_id].sum() for si in group_i)
-        elif key == "group_tpm_mean":
+        elif key == "group_cpm_mean":
             return tuple(
-                self.tpm(kwargs.get("pseudocount", 1))[si, transcript_id].mean()
+                self.cpm(kwargs.get("pseudocount", 1))[si, transcript_id].mean()
                 for si in group_i
             )
         elif key in self.transcripts[transcript_id]:
@@ -893,10 +893,10 @@ def _set_coverage(self, force=False):
         self.data["coverage"] = cov
         self.data["segment_graph"] = None
 
-    def tpm(self, pseudocount=1):
-        """Returns the transcripts per million (TPM).
+    def cpm(self, pseudocount=1):
+        """Returns the counts per million (CPM).
 
-        TPM is returned as a numpy array, with samples in columns and transcript isoforms in the rows.
+        CPM is returned as a numpy array, with samples in columns and transcript isoforms in the rows.
         """
         return (
             (self.coverage + pseudocount)

src/isotools/plots.py

@@ -434,14 +434,14 @@ def plot_saturation(
     """Plots Negative Binomial model to analyze the saturation of LRTS data.
 
     Saturation (e.g. the probability to observe a transcript of interest in the sample) is dependent on the sequencing depth (number of reads),
-    the concentration of the transcripts of interest in the sample (in TPM),
+    the concentration of the transcripts of interest in the sample (in CPM),
     and the requested coverage of the transcript in the data (minimum number of reads per transcript).
     This function models the relation with a Negative Binomial distribution, to help estimate the required sequencing depth.
 
     :param isoseq: If provided, the sequencing depth of samples from this isotools.Transcriptome object are depicted as vertical lines.
     :param ax: The axis for the plot.
     :param cov_th: The requested coverage, e.g. the minimum number of reads per transcript.
-    :param expr_th: A list of transcript concentrations in TPM for transcripts of interest.
+    :param expr_th: A list of transcript concentrations in CPM for transcripts of interest.
     :param x_range: Specify the range of the x axis (e.g. the sequencing depth)
     :param legend: If set True, a legend is added to the plot.
     :param label: If set True, the sample names and sequencing depth from the isoseq parameter is printed in the plot.
@@ -467,11 +467,11 @@ def plot_saturation(
         if isoseq is not None
         else {}
     )
-    for tpm_th in expr_th:
+    for cpm_th in expr_th:
         chance = nbinom.cdf(
-            k - cov_th, n=cov_th, p=tpm_th * 1e-6
+            k - cov_th, n=cov_th, p=cpm_th * 1e-6
         )  # 0 to k-cov_th failiors
-        ax.plot(k / 1e6, chance, label=f"{tpm_th} TPM")
+        ax.plot(k / 1e6, chance, label=f"{cpm_th} CPM")
     for sample, cov in n_reads.items():
         ax.axvline(cov / 1e6, color="grey", ls="--")
         if label:

src/isotools/run_isotools.py

@@ -191,7 +191,7 @@ def main():
         df = isoseq.transcript_table(
             groups=groups,
             coverage=True,
-            tpm=True,
+            cpm=True,
             query=args.filter_query,
             progress_bar=args.progress_bar,
         )

tests/data_import_test.py

@@ -84,6 +84,29 @@ def test_import_bam():
     isoseq.save("tests/data/example_1_isotools.pkl")
 
 
+@pytest.mark.dependency(depends=["test_import_bam"])
+def test_cpm():
+    # regression test for the tpm -> cpm rename: values are counts per
+    # million (count / total_reads * 1e6), not length-normalized TPM --
+    # correct for full-length long reads, where read count already is a
+    # direct proxy for molecule count.
+    isoseq = Transcriptome.load("tests/data/example_1_isotools.pkl")
+    tab = isoseq.transcript_table(coverage=True, cpm=True)
+    stab = isoseq.sample_table.set_index("name")
+    for sample in isoseq.samples:
+        cov_col, cpm_col = f"{sample}_coverage", f"{sample}_cpm"
+        assert cov_col in tab.columns and cpm_col in tab.columns
+        total = stab.loc[sample, "nonchimeric_reads"]
+        expected = tab[cov_col] / total * 1e6
+        assert (
+            tab[cpm_col] - expected
+        ).abs().max() < 1e-6, "cpm should be count/total_reads*1e6"
+
+    gene = next(iter(isoseq.iter_genes(query="EXPRESSED")))
+    gene_cpm = gene.cpm()
+    assert gene_cpm.shape == (len(isoseq.samples), gene.n_transcripts)
+
+
 @pytest.mark.dependency(depends=["test_import_bam"])
 def test_fsm():
     isoseq = Transcriptome.load("tests/data/example_1_isotools.pkl")