From c4a449af840451804f5f647b5c6a5bb093c9c06b Mon Sep 17 00:00:00 2001 From: bruno-ariano Date: Wed, 28 Jan 2026 14:12:53 +0100 Subject: [PATCH 1/3] create modules part1 --- main.nf | 16 +- nextflow.config | 52 +- tdbsumstat/__init__.py | 2 + tdbsumstat/harmonize/__init__.py | 14 + tdbsumstat/harmonize/align_alleles.py | 47 ++ .../harmonize/align_alleles_with_pvar.py | 29 + tdbsumstat/harmonize/create_mapping.py | 26 + tdbsumstat/harmonize/exceptions.py | 2 + tdbsumstat/harmonize/fill_missing_columns.py | 118 ++++ tdbsumstat/harmonize/qc_sumstat.py | 70 ++ tdbsumstat/ingest/__init__.py | 16 + tdbsumstat/ingest/create_metadata.py | 107 +++ tdbsumstat/ingest/create_tiledb.py | 48 ++ tdbsumstat/ingest/export_metadata_to_csv.py | 95 +++ tdbsumstat/ingest/ingest_data.py | 46 ++ tdbsumstat/ingest/merge_metadata_files.py | 92 +++ .../ingest/write_individual_metadata.py | 25 + tdbsumstat/utils/compute_pval.py | 9 + tdbsumstat/utils/harmonize_ingest.py | 668 +----------------- 19 files changed, 791 insertions(+), 691 deletions(-) create mode 100644 tdbsumstat/harmonize/__init__.py create mode 100644 tdbsumstat/harmonize/align_alleles.py create mode 100644 tdbsumstat/harmonize/align_alleles_with_pvar.py create mode 100644 tdbsumstat/harmonize/create_mapping.py create mode 100644 tdbsumstat/harmonize/exceptions.py create mode 100644 tdbsumstat/harmonize/fill_missing_columns.py create mode 100644 tdbsumstat/harmonize/qc_sumstat.py create mode 100644 tdbsumstat/ingest/__init__.py create mode 100644 tdbsumstat/ingest/create_metadata.py create mode 100644 tdbsumstat/ingest/create_tiledb.py create mode 100644 tdbsumstat/ingest/export_metadata_to_csv.py create mode 100644 tdbsumstat/ingest/ingest_data.py create mode 100644 tdbsumstat/ingest/merge_metadata_files.py create mode 100644 tdbsumstat/ingest/write_individual_metadata.py create mode 100644 tdbsumstat/utils/compute_pval.py diff --git a/main.nf b/main.nf index cd24784..9dba474 100644 --- a/main.nf +++ b/main.nf @@ -9,11 +9,11 @@ include { RECOMPUTE_META } from "./modules/recompute_meta" workflow { if (params.ingestion){ - Channel.fromPath(params.file_path_ingestion, checkIfExists:true) + channel.fromPath(params.file_path_ingestion, checkIfExists:true) .splitText(by: params.ingestion_chunk_files, keepHeader: true, file: true) .set { list_files } - mapping_file = Channel.fromPath(params.mapping_file, checkIfExists:true) - create_tiledb = CREATE_TILEDB(mapping_file, Channel.of('dummy')) + mapping_file = channel.fromPath(params.mapping_file, checkIfExists:true) + create_tiledb = CREATE_TILEDB(mapping_file, channel.of('dummy')) // Pass the TileDB array through all ingestion steps ingestion_results = INGEST_DATA(create_tiledb.tiledb_storage, list_files, mapping_file, create_tiledb.dummy_file) // Collect all completion signals @@ -23,12 +23,12 @@ workflow { // Get one instance of the updated TileDB array (they should all be the same) updated_tiledb = ingestion_results.tiledb_updated.first() // After all ingestion is done, merge metadata using the updated TileDB - merged_metadata = MERGE_METADATA(updated_tiledb, mapping_file, all_metadata_parts, all_ingestion_done) + MERGE_METADATA(updated_tiledb, mapping_file, all_metadata_parts, all_ingestion_done) // The final output will be in merged_metadata.tiledb_final } if (params.export){ if (params.snp) { - Channel + channel .fromPath(params.snp, checkIfExists: true) .splitCsv(header: true) .map { row -> @@ -51,7 +51,7 @@ workflow { } } if (params.locusbreaker){ - Channel.fromPath(params.table_lb, checkIfExists:true) + channel.fromPath(params.table_lb, checkIfExists:true) .splitText(by: params.tiledb_batch_size, keepHeader: true, file: true) .map { batch_file -> def batch_index = (batch_file.name =~ /\.(\d+)\.csv$/)[0][1] @@ -61,7 +61,7 @@ workflow { EXPORT_LOCUSBREAKER(tiledb_metadata_batches) } if (params.export_traits){ - Channel.fromPath(params.list_traits, checkIfExists:true) + channel.fromPath(params.list_traits, checkIfExists:true) .splitText(by: params.tiledb_batch_size, keepHeader: true, file: true) .map { batch_file -> def batch_index = (batch_file.name =~ /\.(\d+)\.csv$/)[0][1] @@ -72,7 +72,7 @@ workflow { } if (params.recompute_meta){ - Channel.fromPath(params.list_traits, checkIfExists:true) + channel.fromPath(params.list_traits, checkIfExists:true) .splitText(by: params.tiledb_batch_size, keepHeader: true, file: true) .map { batch_file -> def batch_index = (batch_file.name =~ /\.(\d+)\.csv$/)[0][1] diff --git a/nextflow.config b/nextflow.config index d6d975f..9e81b66 100644 --- a/nextflow.config +++ b/nextflow.config @@ -43,24 +43,33 @@ params { config_profile_url = null config_profile_name = null } +// Execution reporting +timeline { + enabled = true + file = "${params.outdir}/pipeline_info/execution_timeline_${new java.util.Date().format('yyyy-MM-dd_HH-mm-ss')}.html" +} -// Load nf-core custom profiles from different Institutions -try { - includeConfig "${params.custom_config_base}/nfcore_custom.config" -} catch (Exception e) { - System.err.println("WARNING: Could not load nf-core/config profiles: ${params.custom_config_base}/nfcore_custom.config") +report { + enabled = true + file = "${params.outdir}/pipeline_info/execution_report_${new java.util.Date().format('yyyy-MM-dd_HH-mm-ss')}.html" } +trace { + enabled = true + file = "${params.outdir}/pipeline_info/execution_trace_${new java.util.Date().format('yyyy-MM-dd_HH-mm-ss')}.txt" +} + +dag { + enabled = true + file = "${params.outdir}/pipeline_info/pipeline_dag_${new java.util.Date().format('yyyy-MM-dd_HH-mm-ss')}.html" +} +// Load nf-core custom profiles from different Institutions +includeConfig "${params.custom_config_base}/nfcore_custom.config" + // Load base.config by default includeConfig 'conf/base.config' -// Load nf-core institutional configs (like your colleague) -try { - includeConfig "${params.custom_config_base}/nfcore_custom.config" -} catch (Exception e) { - System.err.println("WARNING: Could not load nf-core/config profiles: ${params.custom_config_base}/nfcore_custom.config") -} // Global process configuration for tdbsumstat process { @@ -143,28 +152,7 @@ profiles { test_recompute_meta { includeConfig 'conf/test.config' } } -// Execution reporting -def trace_timestamp = new java.util.Date().format('yyyy-MM-dd_HH-mm-ss') - -timeline { - enabled = true - file = "${params.outdir}/pipeline_info/execution_timeline_${trace_timestamp}.html" -} - -report { - enabled = true - file = "${params.outdir}/pipeline_info/execution_report_${trace_timestamp}.html" -} - -trace { - enabled = true - file = "${params.outdir}/pipeline_info/execution_trace_${trace_timestamp}.txt" -} -dag { - enabled = true - file = "${params.outdir}/pipeline_info/pipeline_dag_${trace_timestamp}.html" -} manifest { name = 'TileDB-Sumstat' diff --git a/tdbsumstat/__init__.py b/tdbsumstat/__init__.py index e69de29..f2ff7c9 100755 --- a/tdbsumstat/__init__.py +++ b/tdbsumstat/__init__.py @@ -0,0 +1,2 @@ +class HarmonizationError(Exception): + pass diff --git a/tdbsumstat/harmonize/__init__.py b/tdbsumstat/harmonize/__init__.py new file mode 100644 index 0000000..d04fa9a --- /dev/null +++ b/tdbsumstat/harmonize/__init__.py @@ -0,0 +1,14 @@ +from .create_mapping import _create_mapping +from .align_alleles_with_pvar import _align_alleles_with_pvar +from .align_alleles import _align_alleles +from .fill_missing_columns import _fill_missing_columns +from .qc_sumstat import _qc_sumstat + +__all__ = [ + "_create_mapping", + "_align_alleles_with_pvar", + "_align_alleles", + "_fill_missing_columns", + "_qc_sumstat" +] + diff --git a/tdbsumstat/harmonize/align_alleles.py b/tdbsumstat/harmonize/align_alleles.py new file mode 100644 index 0000000..1efd474 --- /dev/null +++ b/tdbsumstat/harmonize/align_alleles.py @@ -0,0 +1,47 @@ +import polars as pl +from pathlib import Path +import _align_alleles_with_pvar + +def _align_alleles( + chunk_pl:pl.DataFrame, + pvar_file: Path + ): + #Remove the current SNPID + if "SNPID" in chunk_pl.columns: + chunk_pl = chunk_pl.drop("SNPID") + swap = pl.col("A1") < pl.col("A2") + #A2 is the effect allel while A is the non effect one + chunk_pl = chunk_pl.with_columns([ + pl.when(swap).then(pl.col("A1")).otherwise(pl.col("A2")).alias("EA"), + # set NEA to the larger allele + pl.when(swap).then(pl.col("A2")).otherwise(pl.col("A1")).alias("NEA")] + ) + chunk_pl = chunk_pl.with_columns( + pl.concat_str( + [ + pl.col("CHR"), + pl.col("POS").cast(pl.Utf8), # cast POS if numeric + pl.col("EA"), # lexicographically smaller + pl.col("NEA") # lexicographically larger + ], + separator=":" + ).alias("SNPID") + ) + + if pvar_file: + _align_alleles_with_pvar(chunk_pl) + chunk_pl.drop(["REF","ALT"]) + else: + # flip the sign of BETA when swapping + chunk_pl = chunk_pl.with_columns([ + pl.when(swap).then(-pl.col("BETA")).otherwise(pl.col("BETA")).alias("BETA"), + # flip EAF to 1 - EAF when swapping + pl.when(swap).then(1.0 - pl.col("EAF")).otherwise(pl.col("EAF")).alias("EAF"), + # set EA to the smaller allele + ]) + chunk_pl = chunk_pl.drop(["A1","A2"]) + chunk_pl = chunk_pl.with_columns( + pl.concat_str( + pl.lit("chr"), + pl.col("SNPID")).alias("SNPID")) + return chunk_pl \ No newline at end of file diff --git a/tdbsumstat/harmonize/align_alleles_with_pvar.py b/tdbsumstat/harmonize/align_alleles_with_pvar.py new file mode 100644 index 0000000..8971bd8 --- /dev/null +++ b/tdbsumstat/harmonize/align_alleles_with_pvar.py @@ -0,0 +1,29 @@ +from .exceptions import HarmonizationError +from pathlib import Path +import polars as pl + +def _align_alleles_with_pvar( + pvar_file: Path | None = None): + """Align alleles based on pvar file.""" + if pvar_file is None: + raise HarmonizationError("pvar_file must be provided to verify alleles order") + if not pvar_file.is_file(): + raise FileNotFoundError(f"pvar_file {pvar_file} does not exist") + + #The ALT must correspond to the alternative allele + pvar_df = pl.read_csv( + pvar_file, + separator="\t", + has_header=True, + dtypes={"CHROM": pl.Utf8, "POS": pl.Utf8, "SNPID": pl.Utf8, + "REF": pl.Utf8, "ALT": pl.Utf8}, + ) + #Here we assume the SNPID is alphabetically sortedin both pvar and summary statistics + chunk_pl = chunk_pl.join(pvar_df, on="SNPID", how="inner", suffix="_pvar") + + swap = pl.col("ALT")< pl.col("REF") + + chunk_pl = chunk_pl.with_columns([ + pl.when(swap).then(pl.col("BETA")).otherwise(-pl.col("BETA")), + pl.when(swap).then(pl.col("EAF")).otherwise(1.0-pl.col("EAF")) + ]) diff --git a/tdbsumstat/harmonize/create_mapping.py b/tdbsumstat/harmonize/create_mapping.py new file mode 100644 index 0000000..ed460a4 --- /dev/null +++ b/tdbsumstat/harmonize/create_mapping.py @@ -0,0 +1,26 @@ +from .exceptions import HarmonizationError +import pandas as pd +from pathlib import Path + +def _create_mapping( + sumstat: pd.DataFrame, + mapping_file: Path | None = None + ): + if mapping_file is None: + raise HarmonizationError("No mapping file provided. Provide a mapping file") + if not mapping_file.is_file(): + raise FileNotFoundError(f"Mapping file {mapping_file} does not exist") + + df = pd.read_csv(mapping_file, header=None, names=["key", "value"]) + if df.empty: + raise HarmonizationError("Mapping file is empty or not formatted correctly.") + mapping_types = dict(zip(df["key"], df["value"])) + #check that "BETA", "SE" are in the vlaues of the mapping_types + if not all(col in mapping_types.values() for col in ["BETA", "SE"]): + raise HarmonizationError("Mapping file must contain BETA and SE columns.") + # Check if CHR and POS or SNP are present + if not all(col for col in ["CHR", "POS"] if col in mapping_types.values()): + if "SNPID" not in mapping_types.values(): + raise HarmonizationError("Mapping file must contain either CHR, and POS or SNPID columns.") + chunk_pl = sumstat.rename(mapping_types) + return chunk_pl \ No newline at end of file diff --git a/tdbsumstat/harmonize/exceptions.py b/tdbsumstat/harmonize/exceptions.py new file mode 100644 index 0000000..fb71ebc --- /dev/null +++ b/tdbsumstat/harmonize/exceptions.py @@ -0,0 +1,2 @@ +class HarmonizationError(Exception): + pass \ No newline at end of file diff --git a/tdbsumstat/harmonize/fill_missing_columns.py b/tdbsumstat/harmonize/fill_missing_columns.py new file mode 100644 index 0000000..25687ff --- /dev/null +++ b/tdbsumstat/harmonize/fill_missing_columns.py @@ -0,0 +1,118 @@ +import polars as pl +from tdbsumstat import HarmonizationError +import numpy as np +import _create_mapping +from pathlib import Path + + +def _fill_missing_columns( + chunk_pl: pl.DataFrame, + mapping_file:Path, + type_trait: str, + type_sumstat: str, + mac: int | None = 10, + trait: str | None = None, + cell: str | None = None, + gene: str | None = None, + pheno_var:int | None = None, + n: int | None = None, + n_controls: int | None = None, + n_cases: int | None = None): + + """Load and rename columns, and ensure CHR/POS exist.""" + chunk_pl = _create_mapping(sumstat = chunk_pl, mapping_file = mapping_file) + if "CHR" not in chunk_pl.columns or "POS" not in chunk_pl.columns: + chunk_pl = chunk_pl.with_columns( + pl.col("SNPID") + .str.split_exact(":", 4) + .struct.rename_fields(["CHR", "POS", "A1", "A2"]) + .alias("fields") + ).unnest("fields") + + if "EAF" not in chunk_pl.columns: + chunk_pl = chunk_pl.with_columns(pl.lit(0).alias("EAF")) + if "DIST" not in chunk_pl.columns: + chunk_pl = chunk_pl.with_columns(pl.lit(1).alias("DIST")) + #Check for removing double headers + + if type_trait == "quant": + if not "N" in chunk_pl.columns: + if n is not None: + chunk_pl = chunk_pl.with_columns( + pl.lit(int(n)).alias("N") + ) + else: + raise HarmonizationError("N column is missing and N parameter is not provided") + if mac is not None: + chunk_pl = chunk_pl.with_columns( + (2 * pl.col("N") * pl.min_horizontal(pl.col("EAF"), 1 - pl.col("EAF"))) + .alias("MAC") + ).filter(pl.col("MAC") >= mac) + + chunk_pl = chunk_pl.with_columns( + pl.lit(pheno_var).alias("PHENO_VAR") + ) + elif type_trait == "binary": + if not all(sample_size in chunk_pl.columns for sample_size in ["N_CASES", "N_CONTROLS"]): + if not None in [n_cases, n_controls]: + chunk_pl = chunk_pl.with_columns( + pl.lit(float(n_cases)).alias("N_CASES"), + pl.lit(float(n_controls)).alias("N_CONTROLS"), + pl.lit(float(n_cases) + float(n_controls)).alias("N"), + ) + else: + raise HarmonizationError("n_cases and n_controls columns are missing and were not provided") + if mac is not None: + chunk_pl = chunk_pl.with_columns( + (2 * pl.col("N") * pl.min_horizontal(pl.col("EAF"), 1 - pl.col("EAF"))) + .alias("MAC") + ).filter(pl.col("MAC") >= mac) + else: + raise HarmonizationError("Type of trait must be either binary or quant") + + if type_sumstat=="gwas": + tiledb_types = { + "CHR": np.uint16, + "TRAIT": str, + "POS": np.uint32, + "SNPID": str, + "RSID": str, + "EAF": np.float32, + "BETA": np.float32, + "SE": np.float32, + "P": np.float64, + } + if "TRAIT" not in chunk_pl.columns: + chunk_pl = chunk_pl.with_columns( + pl.lit(trait).alias("TRAIT") + ) + else: + tiledb_types = { + "CHR": np.uint16, + "CELL": str, + "GENE": str, + "POS": np.uint32, + "SNPID": str, + "RSID": str, + "DIST": np.int64, + "EAF": np.float32, + "BETA": np.float32, + "SE": np.float32, + "P": np.float64, + } + if "CELL" not in chunk_pl.columns: + chunk_pl = chunk_pl.with_columns( + pl.lit(cell).alias("CELL") + ) + if "GENE" not in chunk_pl.columns: + chunk_pl = chunk_pl.with_columns( + pl.lit(gene).alias("GENE") + ) + if "RSID" not in chunk_pl.columns: + chunk_pl = chunk_pl.with_columns( + pl.lit("None").alias("RSID") + ) + if "LOG10P" in chunk_pl.columns: + chunk_pl = chunk_pl.with_columns( + (10 ** (-pl.col("LOG10P"))).alias("P") + ) \ No newline at end of file diff --git a/tdbsumstat/harmonize/qc_sumstat.py b/tdbsumstat/harmonize/qc_sumstat.py new file mode 100644 index 0000000..a8812cc --- /dev/null +++ b/tdbsumstat/harmonize/qc_sumstat.py @@ -0,0 +1,70 @@ +import os +import gwaslab as gl +import polars as pl +import pandas as pd + +def _qc_sumstat( + uri: str = None, + type_sumstat: str = None, + type_trait:str = None, + file_path:str = None): + directory = uri + "_logs" + filename = os.path.basename(file_path) # "test.csv.gz" + # Remove all extensions + file_name = filename.split('.')[0] # "test" + + if not os.path.isdir(directory): + os.mkdir(directory) + sumstat_preqc = chunk_pl.to_pandas() + if type_sumstat == "gwas": + if type_trait== "quant": + sumstat_gl =gl.Sumstats(sumstat_preqc, + snpid="SNPID", + chrom="CHR", + pos="POS", + eaf="EAF", + beta="BETA", + se="SE", + p="P", + n="N", + ea = "EA", + nea = "NEA", + other = ["TRAIT","RSID"]) + else: + sumstat_gl =gl.Sumstats(sumstat_preqc, + snpid="SNPID", + chrom="CHR", + pos="POS", + eaf="EAF", + beta="BETA", + se="SE", + p="P", + n="N", + ncase = "N_CASES", + ncontrol = "N_CONTROLS", + ea = "EA", + nea = "NEA", + other = ["TRAIT","RSID"]) + + else: + sumstat_gl =gl.Sumstats(sumstat_preqc, + snpid="SNPID", + chrom="CHR", + pos="POS", + eaf="EAF", + beta="BETA", + se="SE", + p="P", + n="N", + other = ["CELL","GENE","RSID","DIST","PHENO_VAR"]) + #sumstat_gl.fix_id() + sumstat_gl.fix_chr(remove=True) + sumstat_gl.fix_pos(remove=True) + sumstat_gl.fix_allele(remove=True) + sumstat_gl.check_sanity() + sumstat_gl.check_data_consistency() + #sumstat_gl.remove_dup(mode="m") + #sumstat_gl.basic_check(n_cores = 4, remove=True, remove_dup=True) + + sumstat_gl.log.save(directory + "/" + file_name) + chunk_pl = pl.from_pandas(sumstat_gl.data) \ No newline at end of file diff --git a/tdbsumstat/ingest/__init__.py b/tdbsumstat/ingest/__init__.py new file mode 100644 index 0000000..a9aaa97 --- /dev/null +++ b/tdbsumstat/ingest/__init__.py @@ -0,0 +1,16 @@ +from .create_metadata import _create_metadata +from .create_tiledb import _create_tiledb +from .export_metadata_to_csv import _export_metadata_to_csv +from .ingest_data import _ingest_data +from .merge_metadata_files import _merge_metadata_files +from .write_individual_metadata import _write_individual_metadata + + +__all__ = [ + _create_metadata, + _create_tiledb, + _export_metadata_to_csv, + _ingest_data, + _merge_metadata_files, + _write_individual_metadata +] diff --git a/tdbsumstat/ingest/create_metadata.py b/tdbsumstat/ingest/create_metadata.py new file mode 100644 index 0000000..8212e23 --- /dev/null +++ b/tdbsumstat/ingest/create_metadata.py @@ -0,0 +1,107 @@ +from tdbsumstat import HarmonizationError +import polars as pl +from tdbsumstat.utils import acat_optimized +from tdbsumstat.utils import compute_pheno_variance +import logging +import _write_individual_metadata + +logger = logging.getLogger(__name__) +logging.basicConfig(level=logging.INFO, format="%(asctime)s - %(levelname)s - %(message)s") + +def _create_metadata( + type_sumstat:str, + type_trait:str, + file_path: str, + ): + """Create and store metadata as individual JSON files.""" + metadata = { + "traits": [], + "CELL": [] + } + chunk_pl = chunk_pl.drop_nulls() + if type_sumstat == "qtl": + # Get unique cell types + celltypes = chunk_pl["CELL"].unique().to_list() + if not celltypes: + raise HarmonizationError("No cell types found in the data") + + # Update CELL list + metadata["CELL"] = celltypes + + # Process each cell type + for cell in celltypes: + # Filter by this cell type and compute ACAT per gene + df_cell = chunk_pl.filter(pl.col("CELL") == cell) + + # Group by CHR and GENE, compute ACAT for each group + chr_gene_agg = df_cell.group_by(["CHR", "GENE"]).agg([ + pl.col("P").map_batches( + lambda s: pl.Series([acat_optimized(s)]), + return_dtype=pl.Float64 + ).alias("ACAT_LIST"), + pl.col("P").min().alias("min_P"), + pl.col("N").first().alias("N"), + pl.col("PHENO_VAR").first().alias("PHENO_VAR") + ]) + chr_gene_agg = chr_gene_agg.with_columns( + pl.col("ACAT_LIST").list.first().alias("ACAT") + ) + + # Initialize cell structure if not exists + if cell not in metadata: + metadata[cell] = {} + + # Populate metadata with chromosome -> gene structure + for row in chr_gene_agg.iter_rows(named=True): + chrom = str(row["CHR"]) # Convert to string for consistency + gene = row["GENE"] + acat_val = row["ACAT"] + n_val = float(row["N"]) + min_p = float(row["min_P"]) + pheno_val = float(row["PHENO_VAR"]) + + if chrom not in metadata[cell]: + metadata[cell][chrom] = {} + + gene_metadata = { + "ACAT": float(acat_val), + "N": n_val, + "PHENO_VAR": pheno_val, + "MIN_P":min_p + } + + metadata[cell][chrom][gene] = gene_metadata + + else: # GWAS case + if "TRAIT" not in chunk_pl.columns: + raise HarmonizationError("TRAIT column is missing in the data") + + traits = chunk_pl["TRAIT"].unique().to_list() + metadata["traits"] = traits + + for trait in traits: + df_trait = chunk_pl.filter(pl.col("TRAIT") == trait) + pheno_var = compute_pheno_variance(df_trait, type_trait) + n_val = df_trait["N"].unique().to_list()[0] + min_p = df_trait["P"].min().to_list()[0] + + trait_metadata = { + "N": float(n_val), + "PHENO_VAR": float(pheno_var), + "MIN_P": float(min_p) + } + + if type_trait == "binary": + n_cases = df_trait["N_CASES"].unique().to_list()[0] + n_controls = df_trait["N_CONTROLS"].unique().to_list()[0] + trait_metadata.update({ + "N_CASES": float(n_cases), + "N_CONTROLS": float(n_controls) + }) + + metadata[trait] = trait_metadata + + # Write individual metadata file instead of updating TileDB + _write_individual_metadata(metadata, file_path) + + logger.info(f"Individual metadata file created for {file_path}") \ No newline at end of file diff --git a/tdbsumstat/ingest/create_tiledb.py b/tdbsumstat/ingest/create_tiledb.py new file mode 100644 index 0000000..6996b06 --- /dev/null +++ b/tdbsumstat/ingest/create_tiledb.py @@ -0,0 +1,48 @@ +import tiledb +import numpy as np +from .exceptions import HarmonizationError + +def _create_tiledb( + type_sumstat:str, + uri:str + ): + """Create the mapping and dtype definitions.""" + pos_domain = (1, 300000000) # Example range for genomic positions + chr_domain = (1, 24) # Example range for genomic positions + attrs=[ + tiledb.Attr(name="SNPID", dtype="ascii", filters=tiledb.FilterList([tiledb.ZstdFilter(level=5)])), + tiledb.Attr(name="RSID", dtype="ascii", filters=tiledb.FilterList([tiledb.ZstdFilter(level=5)])), + tiledb.Attr(name="EAF", dtype=np.float32, filters=tiledb.FilterList([tiledb.ZstdFilter(level=5)])), + tiledb.Attr(name="BETA", dtype=np.float32, filters=tiledb.FilterList([tiledb.ZstdFilter(level=5)])), + tiledb.Attr(name="SE", dtype=np.float32, filters=tiledb.FilterList([tiledb.ZstdFilter(level=5)])), + tiledb.Attr(name="P", dtype=np.float64, filters=tiledb.FilterList([tiledb.ZstdFilter(level=5)])) + ] + + if type_sumstat == "gwas": + + dimension_tiledb = ["CHR", "TRAIT", "POS"] + dom = tiledb.Domain( + tiledb.Dim(name="CHR", domain = chr_domain, dtype=np.uint16, filters=tiledb.FilterList([tiledb.ZstdFilter(level=5)])), + tiledb.Dim(name="TRAIT", dtype="ascii", var=True, filters=tiledb.FilterList([tiledb.ZstdFilter(level=5)]) ), + tiledb.Dim(name="POS", domain = pos_domain, dtype=np.uint32, filters=tiledb.FilterList([tiledb.ZstdFilter(level=5)])) + ) + else: + + dimension_tiledb = ["CHR", "CELL", "GENE", "POS"] + dom = tiledb.Domain( + tiledb.Dim(name="CHR", domain = chr_domain, dtype=np.uint16, filters=tiledb.FilterList([tiledb.ZstdFilter(level=5)])), + tiledb.Dim(name="CELL", dtype="ascii", var=True, filters=tiledb.FilterList([tiledb.ZstdFilter(level=5)]) ), + tiledb.Dim(name="GENE",dtype="ascii", var=True, filters=tiledb.FilterList([tiledb.ZstdFilter(level=5)]) ), + tiledb.Dim(name="POS", domain = pos_domain, dtype=np.uint32, filters=tiledb.FilterList([tiledb.ZstdFilter(level=5)])) + ) + + attrs = attrs + [ + tiledb.Attr(name="DIST", dtype=np.float32, filters=tiledb.FilterList([tiledb.ZstdFilter(level=5)])) + ] + schema = tiledb.ArraySchema( + domain=dom, + attrs=attrs, + sparse=True, + allows_duplicates=False + ) + tiledb.Array.create(uri, schema) \ No newline at end of file diff --git a/tdbsumstat/ingest/export_metadata_to_csv.py b/tdbsumstat/ingest/export_metadata_to_csv.py new file mode 100644 index 0000000..cb127d6 --- /dev/null +++ b/tdbsumstat/ingest/export_metadata_to_csv.py @@ -0,0 +1,95 @@ +import tiledb +import logging +import json +import pandas as pd + +logger = logging.getLogger(__name__) +logging.basicConfig(level=logging.INFO, format="%(asctime)s - %(levelname)s - %(message)s") + +def _export_metadata_to_csv( + uri: str = None, + output_path: str = None): + + """Export metadata to CSV format.""" + if output_path is None: + output_path = f"{uri}_metadata.csv" + + # Get the merged metadata from TileDB + with tiledb.open(uri, "r") as array: + merged_metadata_json = array.meta.get("merged_metadata", "{}") + + if not merged_metadata_json: + logger.warning("No merged metadata found in TileDB array") + return + + merged_metadata = json.loads(merged_metadata_json) + + # Create rows for CSV + rows = [] + + # Process QTL data (cell types) + if "CELL" in merged_metadata and merged_metadata["CELL"]: + for cell_type in merged_metadata["CELL"]: + if cell_type in merged_metadata: + cell_data = merged_metadata[cell_type] + for chrom, genes in cell_data.items(): + for gene_id, gene_metadata in genes.items(): + row = { + "CHR": chrom, + "CELL": cell_type, + "GENE": gene_id, + "ACAT": gene_metadata.get("ACAT", ""), + "N": gene_metadata.get("N", ""), + "PHENO_VAR": gene_metadata.get("PHENO_VAR", ""), + "MIN_P": gene_metadata.get("MIN_P", "") + } + rows.append(row) + + # Process GWAS data (traits) + if "traits" in merged_metadata and merged_metadata["traits"]: + for trait in merged_metadata["traits"]: + if trait in merged_metadata: + trait_data = merged_metadata[trait] + row = { + "TRAIT": trait, + "N": trait_data.get("N", ""), + "PHENO_VAR": trait_data.get("PHENO_VAR", ""), + "MIN_P": trait_data.get("MIN_P", ""), + "N_CASES": trait_data.get("N_CASES", ""), + "N_CONTROLS": trait_data.get("N_CONTROLS", "") + } + rows.append(row) + + # Create DataFrame and save to CSV + if rows: + df = pd.DataFrame(rows) + + # Reorder columns for better readability + if "CELL" in df.columns: + # QTL format + column_order = ["CHR", "CELL", "GENE", "ACAT", "MIN_P", "N", "PHENO_VAR"] + # Only include columns that exist in the DataFrame + column_order = [col for col in column_order if col in df.columns] + df = df[column_order] + else: + # GWAS format + column_order = ["TRAIT", "N", "PHENO_VAR", "MIN_P", "N_CASES", "N_CONTROLS"] + column_order = [col for col in column_order if col in df.columns] + df = df[column_order] + + df.to_csv(output_path, index=False) + logger.info(f"Metadata exported to {output_path}") + + # Print summary + if "CELL" in df.columns: + logger.info(f"Exported {len(df)} gene-cell type combinations") + logger.info(f"Cell types: {df['CELL'].nunique()}") + logger.info(f"Genes: {df['GENE'].nunique()}") + logger.info(f"Chromosomes: {df['CHR'].nunique()}") + else: + logger.info(f"Exported {len(df)} traits") + + return df + else: + logger.warning("No metadata found to export") + return pd.DataFrame() \ No newline at end of file diff --git a/tdbsumstat/ingest/ingest_data.py b/tdbsumstat/ingest/ingest_data.py new file mode 100644 index 0000000..91c3694 --- /dev/null +++ b/tdbsumstat/ingest/ingest_data.py @@ -0,0 +1,46 @@ +import polars as pl +import tiledb +import logging + +logger = logging.getLogger(__name__) +logging.basicConfig(level=logging.INFO, format="%(asctime)s - %(levelname)s - %(message)s") + +def _ingest_data( + tiledb_types: str = None, + uri: str = None, + dimension_tiledb: int = None, + file_path: str = None, + type_sumstat: str = None): + """Append harmonized data to TileDB.""" + pl.Config.set_tbl_cols(-1) + if type_sumstat == "gwas": + dedup_keys = ["CHR", "POS", "TRAIT"] + else: + dedup_keys = ["CHR", "POS", "CELL", "GENE"] + + # Option A (recommended): window count -> keep only rows whose group count == 1 + chunk_pl = ( + chunk_pl + .with_columns(pl.count().over(dedup_keys).alias("_grp_count")) + .filter(pl.col("_grp_count") == 1) + .drop("_grp_count") + ) + chunk_pl = chunk_pl.with_columns([ + pl.col("CHR").cast(pl.UInt16), + pl.col("POS").cast(pl.UInt32) + ]) + chunk_pl_ingest = chunk_pl.select(tiledb_types.keys()) + chunk_pl_ingest = chunk_pl_ingest.drop_nulls() + try: + tiledb.from_pandas( + uri=uri, + dataframe=chunk_pl_ingest.to_pandas(), + index_dims=dimension_tiledb, + column_types=tiledb_types, + allows_duplicates = False, + mode="append" + ) + logger.info(f"Successfully appended chunk to TileDB for file {file_path}") + except Exception as e: + logger.error(f"Failed to append chunk to TileDB for file {file_path}: {e}") + raise \ No newline at end of file diff --git a/tdbsumstat/ingest/merge_metadata_files.py b/tdbsumstat/ingest/merge_metadata_files.py new file mode 100644 index 0000000..a52e4a9 --- /dev/null +++ b/tdbsumstat/ingest/merge_metadata_files.py @@ -0,0 +1,92 @@ +import os +import logging +from pathlib import Path +import tiledb +import json +import _export_metadata_to_csv +logger = logging.getLogger(__name__) +logging.basicConfig(level=logging.INFO, format="%(asctime)s - %(levelname)s - %(message)s") + + +def _merge_metadata_files( + uri: str = None + ): + """Merge all individual metadata files into final TileDB metadata.""" + metadata_dir = f"{uri}_metadata_parts" + + if not os.path.exists(metadata_dir): + logger.warning(f"No metadata directory found at {metadata_dir}") + return + + merged_metadata = { + "traits": [], + "CELL": [] + } + + # Process all individual metadata files + metadata_files = list(Path(metadata_dir).glob("*.json")) + logger.info(f"Found {len(metadata_files)} metadata files to merge") + + for metadata_file in metadata_files: + try: + with open(metadata_file, 'r') as f: + file_metadata = json.load(f) + + # Merge traits (for GWAS) + if "traits" in file_metadata and file_metadata["traits"]: + current_traits = set(merged_metadata.get("traits", [])) + new_traits = set(file_metadata["traits"]) + merged_metadata["traits"] = list(current_traits.union(new_traits)) + + for trait in file_metadata["traits"]: + if trait in file_metadata: + if trait not in merged_metadata: + merged_metadata[trait] = file_metadata[trait] + else: + logger.warning(f"Trait {trait} already exists in metadata, overwriting") + merged_metadata[trait] = file_metadata[trait] + + # Merge CELL and cell metadata (for QTL) + if "CELL" in file_metadata and file_metadata["CELL"]: + current_cells = set(merged_metadata.get("CELL", [])) + new_cells = set(file_metadata["CELL"]) + merged_metadata["CELL"] = list(current_cells.union(new_cells)) + + for cell in file_metadata["CELL"]: + if cell in file_metadata: + if cell not in merged_metadata: + merged_metadata[cell] = {} + + for chrom, genes in file_metadata[cell].items(): + if chrom not in merged_metadata[cell]: + merged_metadata[cell][chrom] = {} + + for gene, gene_data in genes.items(): + if gene in merged_metadata[cell][chrom]: + logger.warning(f"Gene {gene} already exists in cell {cell} chromosome {chrom}, overwriting") + merged_metadata[cell][chrom][gene] = gene_data + + logger.info(f"Processed {metadata_file.name}") + + except Exception as e: + logger.error(f"Error processing metadata file {metadata_file}: {e}") + continue + + # Store the final merged metadata in TileDB + with tiledb.open(uri, "w") as array: + array.meta["merged_metadata"] = json.dumps(merged_metadata) + + # Export to CSV + _export_metadata_to_csv() + + # Log summary + cell_count = len(merged_metadata.get("CELL", [])) + trait_count = len(merged_metadata.get("traits", [])) + + total_genes = 0 + for cell_type in merged_metadata.get("CELL", []): + if cell_type in merged_metadata: + for chrom_data in merged_metadata[cell_type].values(): + total_genes += len(chrom_data) + + logger.info(f"Final merged metadata: {cell_count} cell types, {trait_count} traits, {total_genes} total genes") diff --git a/tdbsumstat/ingest/write_individual_metadata.py b/tdbsumstat/ingest/write_individual_metadata.py new file mode 100644 index 0000000..ac132ff --- /dev/null +++ b/tdbsumstat/ingest/write_individual_metadata.py @@ -0,0 +1,25 @@ +import json +import os +from pathlib import Path +import logging + +logger = logging.getLogger(__name__) +logging.basicConfig(level=logging.INFO, format="%(asctime)s - %(levelname)s - %(message)s") + + +def _write_individual_metadata( + metadata, + file_path, + uri,): + """Write metadata to individual JSON file.""" + metadata_dir = f"{uri}_metadata_parts" + os.makedirs(metadata_dir, exist_ok=True) + + # Create a safe filename from the original file path + file_stem = Path(file_path).stem + metadata_file = os.path.join(metadata_dir, f"{file_stem}.json") + + with open(metadata_file, 'w') as f: + json.dump(metadata, f, indent=2) + + logger.info(f"Metadata written to {metadata_file}") \ No newline at end of file diff --git a/tdbsumstat/utils/compute_pval.py b/tdbsumstat/utils/compute_pval.py new file mode 100644 index 0000000..48317ff --- /dev/null +++ b/tdbsumstat/utils/compute_pval.py @@ -0,0 +1,9 @@ +from scipy import stats +import polars as pl +def _compute_pval(): + chunk_pl = chunk_pl.with_columns( + (pl.col("BETA") / pl.col("SE")).pow(2).map_batches( + lambda x: pl.Series(stats.chi2.sf(x.to_numpy(), df=1)), + return_dtype=pl.Float64 + ).alias('P') + ) \ No newline at end of file diff --git a/tdbsumstat/utils/harmonize_ingest.py b/tdbsumstat/utils/harmonize_ingest.py index c95bd07..78d2f85 100755 --- a/tdbsumstat/utils/harmonize_ingest.py +++ b/tdbsumstat/utils/harmonize_ingest.py @@ -1,653 +1,19 @@ -import logging -import json +from tdbsumstat.harmonize import _create_mapping, _align_alleles, _fill_missing_columns,_qc_sumstat +from tdbsumstat.ingest import _create_metadata, _create_tiledb, _export_metadata_to_csv,_ingest_data,_merge_metadata_files, _write_individual_metadata from pathlib import Path -import os -import pandas as pd -import polars as pl -import numpy as np -from scipy import stats -import tiledb -import gwaslab as gl -from collections import defaultdict -from tdbsumstat.utils import acat_optimized, compute_pheno_variance - - -logger = logging.getLogger(__name__) -logging.basicConfig(level=logging.INFO, format="%(asctime)s - %(levelname)s - %(message)s") - - -class HarmonizationError(Exception): - """Custom exception for harmonization errors.""" - pass - - class Harmonize: - def __init__(self, mapping_file: str, uri: str, type_sumstat: str, pvar_file: str, type_trait: str, mac: int): - self.mapping_file = mapping_file - self.uri = uri - self.pvar_file = pvar_file - self.type_sumstat = type_sumstat - self.type_trait = type_trait - self.mapping_types = {} - self.tiledb_types = {} - self.dimension_tiledb = [] - self.mac = mac - - def create_mapping(self): - df = pd.read_csv(self.mapping_file, header=None, names=["key", "value"]) - if df.empty: - raise HarmonizationError("Mapping file is empty or not formatted correctly.") - self.mapping_types = dict(zip(df["key"], df["value"])) - #check that "BETA", "SE" are in the vlaues of the mapping_types - if not all(col in self.mapping_types.values() for col in ["BETA", "SE"]): - raise HarmonizationError("Mapping file must contain BETA and SE columns.") - # Check if CHR and POS or SNP are present - if not all(col for col in ["CHR", "POS"] if col in self.mapping_types.values()): - if "SNPID" not in self.mapping_types.values(): - raise HarmonizationError("Mapping file must contain either CHR, and POS or SNPID columns.") - - def create_tiledb(self): - """Create the mapping and dtype definitions.""" - pos_domain = (1, 300000000) # Example range for genomic positions - chr_domain = (1, 24) # Example range for genomic positions - attrs=[ - tiledb.Attr(name="SNPID", dtype="ascii", filters=tiledb.FilterList([tiledb.ZstdFilter(level=5)])), - tiledb.Attr(name="RSID", dtype="ascii", filters=tiledb.FilterList([tiledb.ZstdFilter(level=5)])), - tiledb.Attr(name="EAF", dtype=np.float32, filters=tiledb.FilterList([tiledb.ZstdFilter(level=5)])), - tiledb.Attr(name="BETA", dtype=np.float32, filters=tiledb.FilterList([tiledb.ZstdFilter(level=5)])), - tiledb.Attr(name="SE", dtype=np.float32, filters=tiledb.FilterList([tiledb.ZstdFilter(level=5)])), - tiledb.Attr(name="P", dtype=np.float64, filters=tiledb.FilterList([tiledb.ZstdFilter(level=5)])) - ] - - if self.type_sumstat == "gwas": - - self.dimension_tiledb = ["CHR", "TRAIT", "POS"] - dom = tiledb.Domain( - tiledb.Dim(name="CHR", domain = chr_domain, dtype=np.uint16, filters=tiledb.FilterList([tiledb.ZstdFilter(level=5)])), - tiledb.Dim(name="TRAIT", dtype="ascii", var=True, filters=tiledb.FilterList([tiledb.ZstdFilter(level=5)]) ), - tiledb.Dim(name="POS", domain = pos_domain, dtype=np.uint32, filters=tiledb.FilterList([tiledb.ZstdFilter(level=5)])) - ) - else: - - self.dimension_tiledb = ["CHR", "CELL", "GENE", "POS"] - dom = tiledb.Domain( - tiledb.Dim(name="CHR", domain = chr_domain, dtype=np.uint16, filters=tiledb.FilterList([tiledb.ZstdFilter(level=5)])), - tiledb.Dim(name="CELL", dtype="ascii", var=True, filters=tiledb.FilterList([tiledb.ZstdFilter(level=5)]) ), - tiledb.Dim(name="GENE",dtype="ascii", var=True, filters=tiledb.FilterList([tiledb.ZstdFilter(level=5)]) ), - tiledb.Dim(name="POS", domain = pos_domain, dtype=np.uint32, filters=tiledb.FilterList([tiledb.ZstdFilter(level=5)])) - ) - - attrs = attrs + [ - tiledb.Attr(name="DIST", dtype=np.float32, filters=tiledb.FilterList([tiledb.ZstdFilter(level=5)])) - ] - schema = tiledb.ArraySchema( - domain=dom, - attrs=attrs, - sparse=True, - allows_duplicates=False - ) - tiledb.Array.create(self.uri, schema) - - def align_alleles(self): - """Align alleles based on pvar file.""" - if not self.pvar_file: - raise HarmonizationError("pvar_file must be provided to verify alleles order") - if not Path(self.pvar_file).is_file(): - raise FileNotFoundError(f"pvar_file {self.pvar_file} does not exist") - - #The ALT must correspond to the alternative allele - pvar_df = pl.read_csv( - self.pvar_file, - separator="\t", - has_header=True, - dtypes={"CHROM": pl.Utf8, "POS": pl.Utf8, "SNPID": pl.Utf8, - "REF": pl.Utf8, "ALT": pl.Utf8}, - ) - #Here we assume the SNPID is alphabetically sortedin both pvar and summary statistics - self.chunk_pl = self.chunk_pl.join(pvar_df, on="SNPID", how="inner", suffix="_pvar") - - swap = pl.col("ALT")< pl.col("REF") - - self.chunk_pl = self.chunk_pl.with_columns([ - pl.when(swap).then(pl.col("BETA")).otherwise(-pl.col("BETA")), - pl.when(swap).then(pl.col("EAF")).otherwise(1.0-pl.col("EAF")) - ]) - - def harmonize(self, - sumstat, - mac: int, - trait: str = None, - cell: str = None, - gene: str = None, - pheno_var:int = None, - n: int = None, - n_controls: int = None, - n_cases: int = None): - """Load and rename columns, and ensure CHR/POS exist.""" - self.chunk_pl = sumstat.rename(self.mapping_types) - # If CHR/POS missing, extract from SNP ID - if "CHR" not in self.chunk_pl.columns or "POS" not in self.chunk_pl.columns: - self.chunk_pl = self.chunk_pl.with_columns( - pl.col("SNPID") - .str.split_exact(":", 4) - .struct.rename_fields(["CHR", "POS", "A1", "A2"]) - .alias("fields") - ).unnest("fields") - - if "SNPID" in self.chunk_pl.columns: - self.chunk_pl = self.chunk_pl.drop("SNPID") - if "EAF" not in self.chunk_pl.columns: - self.chunk_pl = self.chunk_pl.with_columns(pl.lit(0).alias("EAF")) - if "DIST" not in self.chunk_pl.columns: - self.chunk_pl = self.chunk_pl.with_columns(pl.lit(1).alias("DIST")) - #Check for removing double headers - - if self.type_trait == "quant": - if not "N" in self.chunk_pl.columns: - if n is not None: - self.chunk_pl = self.chunk_pl.with_columns( - pl.lit(int(n)).alias("N") - ) - else: - raise HarmonizationError("N column is missing and N parameter is not provided") - if self.mac is not None: - self.chunk_pl = self.chunk_pl.with_columns( - (2 * pl.col("N") * pl.min_horizontal(pl.col("EAF"), 1 - pl.col("EAF"))) - .alias("MAC") - ).filter(pl.col("MAC") >= self.mac) - - self.chunk_pl = self.chunk_pl.with_columns( - pl.lit(pheno_var).alias("PHENO_VAR") - ) - elif self.type_trait == "binary": - if not all(sample_size in self.chunk_pl.columns for sample_size in ["N_CASES", "N_CONTROLS"]): - if not None in [n_cases, n_controls]: - self.chunk_pl = self.chunk_pl.with_columns( - pl.lit(float(n_cases)).alias("N_CASES"), - pl.lit(float(n_controls)).alias("N_CONTROLS"), - pl.lit(float(n_cases) + float(n_controls)).alias("N"), - ) - else: - raise HarmonizationError("n_cases and n_controls columns are missing and were not provided") - if self.mac is not None: - self.chunk_pl = self.chunk_pl.with_columns( - (2 * pl.col("N") * pl.min_horizontal(pl.col("EAF"), 1 - pl.col("EAF"))) - .alias("MAC") - ).filter(pl.col("MAC") >= self.mac) - else: - raise HarmonizationError("Type of trait must be either binary or quant") - - - - if "SNPID" in self.chunk_pl.columns: - self.chunk_pl = self.chunk_pl.drop("SNPID") - - - #Here we always assumbe that the SNPs are in REF=A1 and ALT=A2 - swap = pl.col("A1") < pl.col("A2") - - #Start by creating new SNPID aligned - self.chunk_pl = self.chunk_pl.with_columns([ - pl.when(swap).then(pl.col("A1")).otherwise(pl.col("A2")).alias("EA"), - # set NEA to the larger allele - pl.when(swap).then(pl.col("A2")).otherwise(pl.col("A1")).alias("NEA")] - ) - - self.chunk_pl = self.chunk_pl.with_columns( - pl.concat_str( - [ - pl.col("CHR"), - pl.col("POS").cast(pl.Utf8), # cast POS if numeric - pl.col("EA"), # lexicographically smaller - pl.col("NEA") # lexicographically larger - ], - separator=":" - ).alias("SNPID") - ) - - if self.pvar_file: - self.align_alleles() - self.chunk_pl.drop(["REF","ALT"]) - else: - # flip the sign of BETA when swapping - self.chunk_pl = self.chunk_pl.with_columns([ - pl.when(swap).then(-pl.col("BETA")).otherwise(pl.col("BETA")).alias("BETA"), - # flip EAF to 1 - EAF when swapping - pl.when(swap).then(1.0 - pl.col("EAF")).otherwise(pl.col("EAF")).alias("EAF"), - # set EA to the smaller allele - ]) - self.chunk_pl = self.chunk_pl.drop(["A1","A2"]) - - self.chunk_pl = self.chunk_pl.with_columns( - pl.concat_str( - pl.lit("chr"), - pl.col("SNPID")).alias("SNPID")) - - if self.type_sumstat=="gwas": - self.tiledb_types = { - "CHR": np.uint16, - "TRAIT": str, - "POS": np.uint32, - "SNPID": str, - "RSID": str, - "EAF": np.float32, - "BETA": np.float32, - "SE": np.float32, - "P": np.float64, - } - if "TRAIT" not in self.chunk_pl.columns: - self.chunk_pl = self.chunk_pl.with_columns( - pl.lit(trait).alias("TRAIT") - ) - else: - self.tiledb_types = { - "CHR": np.uint16, - "CELL": str, - "GENE": str, - "POS": np.uint32, - "SNPID": str, - "RSID": str, - "DIST": np.int64, - "EAF": np.float32, - "BETA": np.float32, - "SE": np.float32, - "P": np.float64, - } - if "CELL" not in self.chunk_pl.columns: - self.chunk_pl = self.chunk_pl.with_columns( - pl.lit(cell).alias("CELL") - ) - if "GENE" not in self.chunk_pl.columns: - self.chunk_pl = self.chunk_pl.with_columns( - pl.lit(gene).alias("GENE") - ) - if "RSID" not in self.chunk_pl.columns: - self.chunk_pl = self.chunk_pl.with_columns( - pl.lit("None").alias("RSID") - ) - if "LOG10P" in self.chunk_pl.columns: - self.chunk_pl = self.chunk_pl.with_columns( - (10 ** (-pl.col("LOG10P"))).alias("P") - ) - - #Calculate p-value from z-score - self.chunk_pl = self.chunk_pl.drop('P') - self.chunk_pl = self.chunk_pl.with_columns( - (pl.col("BETA") / pl.col("SE")).pow(2).map_batches( - lambda x: pl.Series(stats.chi2.sf(x.to_numpy(), df=1)), - return_dtype=pl.Float64 - ).alias('P') - ) - - def qc_sumstat(self, file_path:str): - directory = self.uri + "_logs" - filename = os.path.basename(file_path) # "test.csv.gz" - # Remove all extensions - file_name = filename.split('.')[0] # "test" - - if not os.path.isdir(directory): - os.mkdir(directory) - sumstat_preqc = self.chunk_pl.to_pandas() - if self.type_sumstat == "gwas": - if self.type_trait== "quant": - sumstat_gl =gl.Sumstats(sumstat_preqc, - snpid="SNPID", - chrom="CHR", - pos="POS", - eaf="EAF", - beta="BETA", - se="SE", - p="P", - n="N", - ea = "EA", - nea = "NEA", - other = ["TRAIT","RSID"]) - else: - sumstat_gl =gl.Sumstats(sumstat_preqc, - snpid="SNPID", - chrom="CHR", - pos="POS", - eaf="EAF", - beta="BETA", - se="SE", - p="P", - n="N", - ncase = "N_CASES", - ncontrol = "N_CONTROLS", - ea = "EA", - nea = "NEA", - other = ["TRAIT","RSID"]) - - else: - sumstat_gl =gl.Sumstats(sumstat_preqc, - snpid="SNPID", - chrom="CHR", - pos="POS", - eaf="EAF", - beta="BETA", - se="SE", - p="P", - n="N", - other = ["CELL","GENE","RSID","DIST","PHENO_VAR"]) - #sumstat_gl.fix_id() - sumstat_gl.fix_chr(remove=True) - sumstat_gl.fix_pos(remove=True) - sumstat_gl.fix_allele(remove=True) - sumstat_gl.check_sanity() - sumstat_gl.check_data_consistency() - #sumstat_gl.remove_dup(mode="m") - #sumstat_gl.basic_check(n_cores = 4, remove=True, remove_dup=True) - - sumstat_gl.log.save(directory + "/" + file_name) - self.chunk_pl = pl.from_pandas(sumstat_gl.data) - - def ingest_data(self, file_path): - """Append harmonized data to TileDB.""" - pl.Config.set_tbl_cols(-1) - if self.type_sumstat == "gwas": - dedup_keys = ["CHR", "POS", "TRAIT"] - else: - dedup_keys = ["CHR", "POS", "CELL", "GENE"] - - # Option A (recommended): window count -> keep only rows whose group count == 1 - self.chunk_pl = ( - self.chunk_pl - .with_columns(pl.count().over(dedup_keys).alias("_grp_count")) - .filter(pl.col("_grp_count") == 1) - .drop("_grp_count") - ) - self.chunk_pl = self.chunk_pl.with_columns([ - pl.col("CHR").cast(pl.UInt16), - pl.col("POS").cast(pl.UInt32) - ]) - chunk_pl_ingest = self.chunk_pl.select(self.tiledb_types.keys()) - chunk_pl_ingest = chunk_pl_ingest.drop_nulls() - try: - tiledb.from_pandas( - uri=self.uri, - dataframe=chunk_pl_ingest.to_pandas(), - index_dims=self.dimension_tiledb, - column_types=self.tiledb_types, - allows_duplicates = False, - mode="append" - ) - logger.info(f"Successfully appended chunk to TileDB for file {file_path}") - except Exception as e: - logger.error(f"Failed to append chunk to TileDB for file {file_path}: {e}") - raise - - - def create_metadata(self, file_path: str): - """Create and store metadata as individual JSON files.""" - metadata = { - "traits": [], - "CELL": [] - } - self.chunk_pl = self.chunk_pl.drop_nulls() - if self.type_sumstat == "qtl": - # Get unique cell types - celltypes = self.chunk_pl["CELL"].unique().to_list() - if not celltypes: - raise HarmonizationError("No cell types found in the data") - - # Update CELL list - metadata["CELL"] = celltypes - - # Process each cell type - for cell in celltypes: - # Filter by this cell type and compute ACAT per gene - df_cell = self.chunk_pl.filter(pl.col("CELL") == cell) - - # Group by CHR and GENE, compute ACAT for each group - chr_gene_agg = df_cell.group_by(["CHR", "GENE"]).agg([ - pl.col("P").map_batches( - lambda s: pl.Series([acat_optimized(s)]), - return_dtype=pl.Float64 - ).alias("ACAT_LIST"), - pl.col("P").min().alias("min_P"), - pl.col("N").first().alias("N"), - pl.col("PHENO_VAR").first().alias("PHENO_VAR") - ]) - chr_gene_agg = chr_gene_agg.with_columns( - pl.col("ACAT_LIST").list.first().alias("ACAT") - ) - - # Initialize cell structure if not exists - if cell not in metadata: - metadata[cell] = {} - - # Populate metadata with chromosome -> gene structure - for row in chr_gene_agg.iter_rows(named=True): - chrom = str(row["CHR"]) # Convert to string for consistency - gene = row["GENE"] - acat_val = row["ACAT"] - n_val = float(row["N"]) - min_p = float(row["min_P"]) - pheno_val = float(row["PHENO_VAR"]) - - if chrom not in metadata[cell]: - metadata[cell][chrom] = {} - - gene_metadata = { - "ACAT": float(acat_val), - "N": n_val, - "PHENO_VAR": pheno_val, - "MIN_P":min_p - } - - metadata[cell][chrom][gene] = gene_metadata - - else: # GWAS case - if "TRAIT" not in self.chunk_pl.columns: - raise HarmonizationError("TRAIT column is missing in the data") - - traits = self.chunk_pl["TRAIT"].unique().to_list() - metadata["traits"] = traits - - for trait in traits: - df_trait = self.chunk_pl.filter(pl.col("TRAIT") == trait) - pheno_var = compute_pheno_variance(df_trait, self.type_trait) - n_val = df_trait["N"].unique().to_list()[0] - min_p = df_trait["P"].min().to_list()[0] - - trait_metadata = { - "N": float(n_val), - "PHENO_VAR": float(pheno_var), - "MIN_P": float(min_p) - } - - if self.type_trait == "binary": - n_cases = df_trait["N_CASES"].unique().to_list()[0] - n_controls = df_trait["N_CONTROLS"].unique().to_list()[0] - trait_metadata.update({ - "N_CASES": float(n_cases), - "N_CONTROLS": float(n_controls) - }) - - metadata[trait] = trait_metadata - - # Write individual metadata file instead of updating TileDB - self._write_individual_metadata(metadata, file_path) - - logger.info(f"Individual metadata file created for {file_path}") - - def _write_individual_metadata(self, metadata, file_path): - """Write metadata to individual JSON file.""" - metadata_dir = f"{self.uri}_metadata_parts" - os.makedirs(metadata_dir, exist_ok=True) - - # Create a safe filename from the original file path - file_stem = Path(file_path).stem - metadata_file = os.path.join(metadata_dir, f"{file_stem}.json") - - with open(metadata_file, 'w') as f: - json.dump(metadata, f, indent=2) - - logger.info(f"Metadata written to {metadata_file}") - - def export_metadata_to_csv(self, output_path: str = None): - """Export metadata to CSV format.""" - if output_path is None: - output_path = f"{self.uri}_metadata.csv" - - # Get the merged metadata from TileDB - with tiledb.open(self.uri, "r") as array: - merged_metadata_json = array.meta.get("merged_metadata", "{}") - - if not merged_metadata_json: - logger.warning("No merged metadata found in TileDB array") - return - - merged_metadata = json.loads(merged_metadata_json) - - # Create rows for CSV - rows = [] - - # Process QTL data (cell types) - if "CELL" in merged_metadata and merged_metadata["CELL"]: - for cell_type in merged_metadata["CELL"]: - if cell_type in merged_metadata: - cell_data = merged_metadata[cell_type] - for chrom, genes in cell_data.items(): - for gene_id, gene_metadata in genes.items(): - row = { - "CHR": chrom, - "CELL": cell_type, - "GENE": gene_id, - "ACAT": gene_metadata.get("ACAT", ""), - "N": gene_metadata.get("N", ""), - "PHENO_VAR": gene_metadata.get("PHENO_VAR", ""), - "MIN_P": gene_metadata.get("MIN_P", "") - } - rows.append(row) - - # Process GWAS data (traits) - if "traits" in merged_metadata and merged_metadata["traits"]: - for trait in merged_metadata["traits"]: - if trait in merged_metadata: - trait_data = merged_metadata[trait] - row = { - "TRAIT": trait, - "N": trait_data.get("N", ""), - "PHENO_VAR": trait_data.get("PHENO_VAR", ""), - "MIN_P": trait_data.get("MIN_P", ""), - "N_CASES": trait_data.get("N_CASES", ""), - "N_CONTROLS": trait_data.get("N_CONTROLS", "") - } - rows.append(row) - - # Create DataFrame and save to CSV - if rows: - df = pd.DataFrame(rows) - - # Reorder columns for better readability - if "CELL" in df.columns: - # QTL format - column_order = ["CHR", "CELL", "GENE", "ACAT", "MIN_P", "N", "PHENO_VAR"] - # Only include columns that exist in the DataFrame - column_order = [col for col in column_order if col in df.columns] - df = df[column_order] - else: - # GWAS format - column_order = ["TRAIT", "N", "PHENO_VAR", "MIN_P", "N_CASES", "N_CONTROLS"] - column_order = [col for col in column_order if col in df.columns] - df = df[column_order] - - df.to_csv(output_path, index=False) - logger.info(f"Metadata exported to {output_path}") - - # Print summary - if "CELL" in df.columns: - logger.info(f"Exported {len(df)} gene-cell type combinations") - logger.info(f"Cell types: {df['CELL'].nunique()}") - logger.info(f"Genes: {df['GENE'].nunique()}") - logger.info(f"Chromosomes: {df['CHR'].nunique()}") - else: - logger.info(f"Exported {len(df)} traits") - - return df - else: - logger.warning("No metadata found to export") - return pd.DataFrame() - - def merge_metadata_files(self): - """Merge all individual metadata files into final TileDB metadata.""" - metadata_dir = f"{self.uri}_metadata_parts" - - if not os.path.exists(metadata_dir): - logger.warning(f"No metadata directory found at {metadata_dir}") - return - - merged_metadata = { - "traits": [], - "CELL": [] - } - - # Process all individual metadata files - metadata_files = list(Path(metadata_dir).glob("*.json")) - logger.info(f"Found {len(metadata_files)} metadata files to merge") - - for metadata_file in metadata_files: - try: - with open(metadata_file, 'r') as f: - file_metadata = json.load(f) - - # Merge traits (for GWAS) - if "traits" in file_metadata and file_metadata["traits"]: - current_traits = set(merged_metadata.get("traits", [])) - new_traits = set(file_metadata["traits"]) - merged_metadata["traits"] = list(current_traits.union(new_traits)) - - for trait in file_metadata["traits"]: - if trait in file_metadata: - if trait not in merged_metadata: - merged_metadata[trait] = file_metadata[trait] - else: - logger.warning(f"Trait {trait} already exists in metadata, overwriting") - merged_metadata[trait] = file_metadata[trait] - - # Merge CELL and cell metadata (for QTL) - if "CELL" in file_metadata and file_metadata["CELL"]: - current_cells = set(merged_metadata.get("CELL", [])) - new_cells = set(file_metadata["CELL"]) - merged_metadata["CELL"] = list(current_cells.union(new_cells)) - - for cell in file_metadata["CELL"]: - if cell in file_metadata: - if cell not in merged_metadata: - merged_metadata[cell] = {} - - for chrom, genes in file_metadata[cell].items(): - if chrom not in merged_metadata[cell]: - merged_metadata[cell][chrom] = {} - - for gene, gene_data in genes.items(): - if gene in merged_metadata[cell][chrom]: - logger.warning(f"Gene {gene} already exists in cell {cell} chromosome {chrom}, overwriting") - merged_metadata[cell][chrom][gene] = gene_data - - logger.info(f"Processed {metadata_file.name}") - - except Exception as e: - logger.error(f"Error processing metadata file {metadata_file}: {e}") - continue - - # Store the final merged metadata in TileDB - with tiledb.open(self.uri, "w") as array: - array.meta["merged_metadata"] = json.dumps(merged_metadata) - - # Export to CSV - self.export_metadata_to_csv() - - # Log summary - cell_count = len(merged_metadata.get("CELL", [])) - trait_count = len(merged_metadata.get("traits", [])) - - total_genes = 0 - for cell_type in merged_metadata.get("CELL", []): - if cell_type in merged_metadata: - for chrom_data in merged_metadata[cell_type].values(): - total_genes += len(chrom_data) - - logger.info(f"Final merged metadata: {cell_count} cell types, {trait_count} traits, {total_genes} total genes") - + def __init__( + self, + mapping_file: Path, + uri: str, + type_sumstat: str, + type_trait: str, + mac: int, + pvar_file: Path | None = None): + mapping_file = mapping_file + uri = uri + pvar_file = pvar_file + type_sumstat = type_sumstat + type_trait = type_trait + mac = mac + From 257bd16db2edc8cd1d82666b4d65b17266b36133 Mon Sep 17 00:00:00 2001 From: bruno-ariano Date: Wed, 29 Jul 2026 10:50:25 +0200 Subject: [PATCH 2/3] change delimiter trait export --- tdbsumstat/cli/export.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tdbsumstat/cli/export.py b/tdbsumstat/cli/export.py index 4811834..b864810 100755 --- a/tdbsumstat/cli/export.py +++ b/tdbsumstat/cli/export.py @@ -267,7 +267,7 @@ def query_spec(uri_path, chrom:int, trait: str = None, cell: str = None, gene: s attrs=attr.split(",") ).df[:, trait_list_np , :] else: - trait_list_pd[['cell','gene']] = trait_list_pd['TRAIT'].str.split(':', expand = True) + trait_list_pd[['cell','gene']] = trait_list_pd['TRAIT'].str.split('~', expand = True) cells = trait_list_pd['cell'].to_list() gene = trait_list_pd['gene'].to_list() tiledb_iterator = A.query( From 283cb762e46e91cfe1b8a2716163bec35583bef5 Mon Sep 17 00:00:00 2001 From: bruno-ariano Date: Wed, 29 Jul 2026 11:49:34 +0200 Subject: [PATCH 3/3] ADD beta and SNP ID related to min pvalue --- tdbsumstat/cli/export.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/tdbsumstat/cli/export.py b/tdbsumstat/cli/export.py index b864810..fcb5487 100755 --- a/tdbsumstat/cli/export.py +++ b/tdbsumstat/cli/export.py @@ -249,7 +249,10 @@ def query_spec(uri_path, chrom:int, trait: str = None, cell: str = None, gene: s return_dtype=pl.Float64 ).alias("ACAT_LIST"), pl.col("N").first().alias("N"), - pl.min("P").alias("MIN_P") + pl.min("P").alias("MIN_P"), + pl.min("P").alias("MIN_BETA"), + pl.col("SNP").sort_by("P").first().alias("MIN_P_SNP"), + pl.col("BETA").sort_by("P").first().alias("MIN_P_BETA") ]) chr_gene_agg = chr_gene_agg.with_columns( pl.col("ACAT_LIST").list.first().alias("ACAT"), @@ -275,6 +278,7 @@ def query_spec(uri_path, chrom:int, trait: str = None, cell: str = None, gene: s attrs=attr.split(",") ).df[:, cells, gene , :] + for chunk in tiledb_iterator: chunk.to_csv(f"{out}_{batch_name}.csv", mode="a", index=False, header = False) print(f"Saved filtered summary statistics in {out}")