From 0d91e88ce7de4f19a01e459bae8ee4692cd0bf87 Mon Sep 17 00:00:00 2001 From: tn Date: Fri, 6 Feb 2026 09:46:31 -0500 Subject: [PATCH 1/4] scale factor --- src/pyrecount/accessor.py | 45 ++++++++++++++++++++++++++++++++++++--- 1 file changed, 42 insertions(+), 3 deletions(-) diff --git a/src/pyrecount/accessor.py b/src/pyrecount/accessor.py index 9cea22f..9a7bccb 100644 --- a/src/pyrecount/accessor.py +++ b/src/pyrecount/accessor.py @@ -86,7 +86,6 @@ def get_project_urls(self, dtype) -> List[str]: dtype=dtype, annotation=self.annotation, project_ids=self.project_ids, - # only needed for bigwig sample=self.sample, jxn_format=self.jxn_format, ) @@ -108,6 +107,7 @@ async def cache(self) -> None: makedirs(path.dirname(fpath), exist_ok=True) tasks.append(download_url_to_path(url=url, fpath=fpath)) + # launches all tasks, without rate limiting if tasks: await asyncio.gather(*tasks) @@ -118,7 +118,47 @@ def load(self, dtype) -> Union[pl.DataFrame, Tuple[pl.DataFrame, pl.DataFrame]]: case Dtype.JXN: return self._jxn_load() case Dtype.GENE: - return self._gene_load() + # XXX: + # TODO: force sample name alignment before scaling + + targetSize = 4e7 + L = 100 + + md = self._metadata_load() + mapped_reads = md["star.all_mapped_reads"].cast(pl.Float64) + avg_mapped_read_length = md["star.average_mapped_length"].cast( + pl.Float64 + ) + avg_read_length = md["avg_len"].cast(pl.Float64) + + ratio = (avg_mapped_read_length / avg_read_length).round(0) + ratio = ratio.set(~ratio.is_in([1, 2]), None) + + paired_end = ratio == 2 + paired_factor = pl.when(paired_end).then(2).otherwise(1) + + scaleFactor = md.select( + ( + targetSize + * L + * paired_factor + / (mapped_reads * avg_mapped_read_length**2) + ).alias("scaleFactor") + ).to_series() + + _, counts = self._gene_load() + + numeric_cols = counts.select(pl.selectors.numeric()).columns + scaled = counts.with_columns( + [ + (pl.col(c) * sf).round(0) + for c, sf in zip(numeric_cols, scaleFactor) + ] + ) + + print(scaled) + + return pl.DataFrame(), pl.DataFrame() case Dtype.EXON: return self._exon_load() case Dtype.BW: @@ -366,7 +406,6 @@ def _exon_load(self) -> pl.DataFrame: annotation = self._gtf_read(fpath) if url.endswith(f"{self.annotation.value}.gz"): counts = self._counts_read(fpath) - # TODO: extract first column (chromosome|start_1base|end_1ba…) exon_colname = counts.columns[0] exon_fields = ["chrom", "start", "end", "strand"] counts = ( From 71c0bb3a91aefea5ea02646c4abc6ebba3dd5322 Mon Sep 17 00:00:00 2001 From: tn Date: Sun, 8 Feb 2026 08:30:45 -0500 Subject: [PATCH 2/4] auc --- src/pyrecount/accessor.py | 95 ++++++++++++++++++++++----------------- 1 file changed, 54 insertions(+), 41 deletions(-) diff --git a/src/pyrecount/accessor.py b/src/pyrecount/accessor.py index 9a7bccb..472075d 100644 --- a/src/pyrecount/accessor.py +++ b/src/pyrecount/accessor.py @@ -111,6 +111,57 @@ async def cache(self) -> None: if tasks: await asyncio.gather(*tasks) + def _scale_mapped_reads(self, counts, target_size, L): + md = self._metadata_load() + + mapped_reads = pl.col("star.all_mapped_reads").cast(pl.Float64) + avg_mapped_len = pl.col("star.average_mapped_length").cast(pl.Float64) + avg_read_len = pl.col("avg_len").cast(pl.Float64) + + # paired-end detection + ratio = (avg_mapped_len / avg_read_len).round(0) + paired_end = ratio == 2 + paired_factor = pl.when(paired_end).then(2).otherwise(1) + + # scale factors + sf = md.select( + [ + pl.col("external_id"), + ( + target_size * L * paired_factor / (mapped_reads * avg_mapped_len**2) + ).alias("sf"), + ] + ) + + sf_map = dict(zip(sf["external_id"], sf["sf"])) + + return counts.with_columns( + [ + (pl.col(c) * sf_map[c]) + for c in counts.select(pl.selectors.numeric()).columns + ] + ) + + def _scale_auc(self, counts, target_size): + md = self._metadata_load() + + auc = pl.col("bc_auc.all_reads_all_bases").cast(pl.Float64) + + sf = md.select( + "external_id", + (target_size / auc).alias("sf"), + ) + + sf_map = dict(zip(sf["external_id"], sf["sf"])) + + return counts.with_columns( + [ + pl.col(c).mul(sf_map[c]).round(0).cast(pl.Int64).alias(c) + for c in counts.columns + if c != "gene_id" + ] + ) + def load(self, dtype) -> Union[pl.DataFrame, Tuple[pl.DataFrame, pl.DataFrame]]: match dtype: case Dtype.METADATA: @@ -118,47 +169,9 @@ def load(self, dtype) -> Union[pl.DataFrame, Tuple[pl.DataFrame, pl.DataFrame]]: case Dtype.JXN: return self._jxn_load() case Dtype.GENE: - # XXX: - # TODO: force sample name alignment before scaling - - targetSize = 4e7 - L = 100 - - md = self._metadata_load() - mapped_reads = md["star.all_mapped_reads"].cast(pl.Float64) - avg_mapped_read_length = md["star.average_mapped_length"].cast( - pl.Float64 - ) - avg_read_length = md["avg_len"].cast(pl.Float64) - - ratio = (avg_mapped_read_length / avg_read_length).round(0) - ratio = ratio.set(~ratio.is_in([1, 2]), None) - - paired_end = ratio == 2 - paired_factor = pl.when(paired_end).then(2).otherwise(1) - - scaleFactor = md.select( - ( - targetSize - * L - * paired_factor - / (mapped_reads * avg_mapped_read_length**2) - ).alias("scaleFactor") - ).to_series() - - _, counts = self._gene_load() - - numeric_cols = counts.select(pl.selectors.numeric()).columns - scaled = counts.with_columns( - [ - (pl.col(c) * sf).round(0) - for c, sf in zip(numeric_cols, scaleFactor) - ] - ) - - print(scaled) - - return pl.DataFrame(), pl.DataFrame() + annotation, counts = self._gene_load() + target_size = 4e7 + return annotation, self._scale_auc(counts, target_size) case Dtype.EXON: return self._exon_load() case Dtype.BW: From 809446ef59a999e6a809a5f5b021cf6622f2e912 Mon Sep 17 00:00:00 2001 From: drew Date: Sun, 8 Feb 2026 09:56:24 -0500 Subject: [PATCH 3/4] fix gene test --- src/pyrecount/accessor.py | 9 ++++++--- tests/test_accessor.py | 6 +++--- 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/src/pyrecount/accessor.py b/src/pyrecount/accessor.py index 472075d..433a173 100644 --- a/src/pyrecount/accessor.py +++ b/src/pyrecount/accessor.py @@ -147,6 +147,7 @@ def _scale_auc(self, counts, target_size): auc = pl.col("bc_auc.all_reads_all_bases").cast(pl.Float64) + # scale factors sf = md.select( "external_id", (target_size / auc).alias("sf"), @@ -368,9 +369,11 @@ def _gtf_read(self, rpath: str) -> pl.DataFrame: [ annotation_dataframe["attribute"] .map_elements( - lambda x: re.findall(rf'{field} "([^"]*)"', x)[0] - if rf'{field} "' in x - else "", + lambda x: ( + re.findall(rf'{field} "([^"]*)"', x)[0] + if rf'{field} "' in x + else "" + ), return_dtype=pl.Utf8, ) .alias(field) diff --git a/tests/test_accessor.py b/tests/test_accessor.py index 07359e6..8ac5cae 100644 --- a/tests/test_accessor.py +++ b/tests/test_accessor.py @@ -241,20 +241,20 @@ async def test_project_gene_accessor( (pl.col("project").is_in(project_ids)) ) - dtype = Dtype.GENE + dtype = [Dtype.METADATA, Dtype.GENE] gene = Project( metadata=project_meta_dataframe, dbase=dbase, organism=organism, - dtype=[dtype], + dtype=dtype, annotation=annotation, jxn_format=None, root_url=root_url, ) await gene.cache() - gene_annotation, gene_counts = gene.load(dtype) + gene_annotation, gene_counts = gene.load(Dtype.GENE) assert gene_annotation.shape == expected_annotation_shape assert gene_counts.shape == expected_counts_shape From f3dc0009611be350d4629b6f54cd84754b13b6ee Mon Sep 17 00:00:00 2001 From: tn Date: Mon, 9 Feb 2026 10:18:54 -0500 Subject: [PATCH 4/4] internal metadata state --- examples/example.py | 17 ++++++++++++++++- src/pyrecount/accessor.py | 23 ++++++++++++++--------- tests/test_accessor.py | 1 - 3 files changed, 30 insertions(+), 11 deletions(-) diff --git a/examples/example.py b/examples/example.py index 2ef849a..4a5b9ea 100755 --- a/examples/example.py +++ b/examples/example.py @@ -26,7 +26,7 @@ # subset project_meta_dataframe = recount_meta_dataframe.filter( - pl.col("project").is_in(["SRP009615"]) & pl.col("external_id").is_in(["SRR389077"]) + pl.col("project").is_in(["SRP009615"]) ) print(project_meta_dataframe) @@ -61,6 +61,21 @@ print(gene_annotation) print(gene_counts) +scaled_counts = project.scale_mapped_reads( + gene_counts, + target_size=4e7, + L=100, +) + +print(scaled_counts) + +scaled_counts = project.scale_auc( + gene_counts, + target_size=4e7, +) + +print(scaled_counts) + exon_annotation, exon_counts = project.load(Dtype.EXON) print(exon_annotation) diff --git a/src/pyrecount/accessor.py b/src/pyrecount/accessor.py index 433a173..42d7a28 100644 --- a/src/pyrecount/accessor.py +++ b/src/pyrecount/accessor.py @@ -47,7 +47,9 @@ class Project: project_ids: List[str] = field(init=False) sample: List[str] = field(init=False) - endpoints: "EndpointConnector" = field(init=False) + endpoints: EndpointConnector = field(init=False) + + _cached_metadata: Optional[pl.DataFrame] = field(init=False, default=None) def __post_init__(self): if not isinstance(self.metadata, pl.DataFrame): @@ -111,8 +113,8 @@ async def cache(self) -> None: if tasks: await asyncio.gather(*tasks) - def _scale_mapped_reads(self, counts, target_size, L): - md = self._metadata_load() + def scale_mapped_reads(self, counts: pl.DataFrame, target_size: float, L: int): + md = self.load_metadata() mapped_reads = pl.col("star.all_mapped_reads").cast(pl.Float64) avg_mapped_len = pl.col("star.average_mapped_length").cast(pl.Float64) @@ -142,8 +144,8 @@ def _scale_mapped_reads(self, counts, target_size, L): ] ) - def _scale_auc(self, counts, target_size): - md = self._metadata_load() + def scale_auc(self, counts: pl.DataFrame, target_size: float) -> pl.DataFrame: + md = self.load_metadata() auc = pl.col("bc_auc.all_reads_all_bases").cast(pl.Float64) @@ -163,16 +165,19 @@ def _scale_auc(self, counts, target_size): ] ) + def load_metadata(self) -> pl.DataFrame: + if self._cached_metadata is None: + self._cached_metadata = self._metadata_load() + return self._cached_metadata + def load(self, dtype) -> Union[pl.DataFrame, Tuple[pl.DataFrame, pl.DataFrame]]: match dtype: case Dtype.METADATA: - return self._metadata_load() + return self.load_metadata() case Dtype.JXN: return self._jxn_load() case Dtype.GENE: - annotation, counts = self._gene_load() - target_size = 4e7 - return annotation, self._scale_auc(counts, target_size) + return self._gene_load() case Dtype.EXON: return self._exon_load() case Dtype.BW: diff --git a/tests/test_accessor.py b/tests/test_accessor.py index 8ac5cae..15f32fe 100644 --- a/tests/test_accessor.py +++ b/tests/test_accessor.py @@ -5,7 +5,6 @@ from pyrecount.accessor import Metadata, Project from pyrecount.models import Dtype, Annotation -# TODO: transform raw counts # TODO: jxn dataframe headers, external_id not rail id # TODO: multi-project support for exon, gene dtypes # TODO: expand sra attributes