diff --git a/docker/Dockerfile b/docker/Dockerfile index 985e7ab..5f22e1b 100644 --- a/docker/Dockerfile +++ b/docker/Dockerfile @@ -3,12 +3,16 @@ FROM eu.gcr.io/finngen-refinery-dev/bioinformatics:0.8 RUN apt-get update && apt-get upgrade --yes RUN curl https://clickhouse.com/ | sh && ./clickhouse install RUN apt install pigz -ADD ./docker/requirements.txt requirements.txt -RUN pip3 install --upgrade scipy -RUN pip3 install -r requirements.txt -ADD finngen_qc finngen_qc -ADD core core -ADD sb_release sb_release -ADD qc_scripts qc_scripts +# Install kanta package +RUN curl -LsSf https://astral.sh/uv/install.sh | sh +ENV PATH="/root/.local/bin:$PATH" +RUN uv python install 3.14 +ADD pyproject.toml pyproject.toml +ADD README.md README.md +ADD uv.lock uv.lock +ADD src src +ENV UV_PROJECT_ENVIRONMENT="/opt/kanta-venv" +RUN uv sync --python 3.14 --locked +ENV PATH="/opt/kanta-venv/bin:$PATH" diff --git a/docker/requirements.txt b/docker/requirements.txt deleted file mode 100644 index d37a3f4..0000000 --- a/docker/requirements.txt +++ /dev/null @@ -1,17 +0,0 @@ -# Progress bars for loops -tqdm -psutil - -# TODO(Vincent 2024-05-03) Upgrade to numpy 2.0.0, should be out soon. -# Note: 2.2.3 is already out and specified here! -numpy==2.2.3 - -# TODO(Vincent 2024-05-03) -# Could add the "compression" extra if we want to use the faster Zstandard compression. -# More info: https://pandas.pydata.org/pandas-docs/stable/getting_started/install.html#compression -pandas[pyarrow]==2.2.2 -duckdb==1.1.3 - -# FOR PLOTTING -matplotlib>=3.9.0 -seaborn>=0.13.2 diff --git a/src/kanta/intake/__main__.py b/src/kanta/intake/__main__.py index c1e23e3..d1d6f5e 100644 --- a/src/kanta/intake/__main__.py +++ b/src/kanta/intake/__main__.py @@ -53,6 +53,9 @@ ) output.check_safe_write(output_file_assemble_stage) output.check_safe_write(output_file_tidyup_stage) + output.check_safe_write( + output_file_tidyup_stage.with_name(f"{output_file_tidyup_stage.stem}_duplicates.parquet") + ) tmp_dir = output.create_tmp_dir() diff --git a/src/kanta/intake/tidyup.py b/src/kanta/intake/tidyup.py index 935016e..2448798 100644 --- a/src/kanta/intake/tidyup.py +++ b/src/kanta/intake/tidyup.py @@ -5,7 +5,6 @@ Differences from the WDL implementation ======================================= -- No logging of duplicates/err lines. - Outputs to a single parquet file, no .txt.gz, as this is very slow. @@ -18,19 +17,26 @@ Lowest tested working spec: 8 CPUs / 8 GB RAM, 32 buckets. Runs in 6-12 min. -If failing due to OOM in the sort+dedup stage, try increasing the bucket count. +If failing due to OOM, try increasing the bucket count — both the sort+dedup +stage and the final join+sanitize stage process one bucket at a time. The GCP VM type appears to matter. N2D is about 2x faster than E2. """ from argparse import ArgumentParser +from datetime import datetime from pathlib import Path import polars as pl +import pyarrow.parquet as pq from kanta import output +def log_step(message: str) -> None: + print(f"[{datetime.now():%Y-%m-%d %H:%M:%S}] {message}", flush=True) + + COLUMNS_UNIQUENESS_SORT = [ "FINNGENID", "APPROX_EVENT_DAY", @@ -40,6 +46,8 @@ "tutkimusvastauksentila", "tutkimustulosarvo", "tutkimustulosyksikko", + "tuloksenpoikkeavuus", + "tutkimustulosteksti" ] @@ -66,24 +74,41 @@ def main( tmp_dir_sort_dedup = tmp_dir / "sort_dedup" tmp_dir_sort_dedup.mkdir() - print("# Consolidate") + tmp_dir_duplicates = tmp_dir / "duplicates" + tmp_dir_duplicates.mkdir() + + log_step("# Consolidate: start") consolidated_file = consolidate_columns(assembled_file, tmp_file_consolidate) + log_step("# Consolidate: done") - print("# Partition") + log_step("# Partition: start") partition(consolidated_file, tmp_dir_partition, partition_n_buckets) - - print("# Sort + Dedup") - for bucket_file in tmp_dir_partition.glob("bucket_id__*.parquet"): - ( - pl.scan_parquet(bucket_file) - .pipe(sort_dedup) - .sink_parquet(tmp_dir_sort_dedup / bucket_file.name) - ) - - print("# Concatenate + join SEX") - bucket_files = [] - for bucket_id in range(partition_n_buckets): - bucket_files.append(tmp_dir_sort_dedup / f"bucket_id__{bucket_id}.parquet") + log_step("# Partition: done") + + log_step("# Sort + Dedup: start") + bucket_files_to_process = sorted(tmp_dir_partition.glob("bucket_id__*.parquet")) + total_n_before = 0 + total_n_after = 0 + bucket_row_counts = {} + for i, bucket_file in enumerate(bucket_files_to_process, start=1): + log_step(f"# Sort + Dedup: bucket {i}/{len(bucket_files_to_process)} ({bucket_file.name})") + n_before = pl.scan_parquet(bucket_file).select(pl.len()).collect().item() + output_bucket_file = tmp_dir_sort_dedup / bucket_file.name + kept, dropped = sort_dedup(pl.scan_parquet(bucket_file)) + kept.sink_parquet(output_bucket_file) + dropped.sink_parquet(tmp_dir_duplicates / bucket_file.name) + n_after = pl.scan_parquet(output_bucket_file).select(pl.len()).collect().item() + total_n_before += n_before + total_n_after += n_after + bucket_row_counts[bucket_file.name] = n_after + n_removed = total_n_before - total_n_after + pct_removed = 100 * n_removed / total_n_before if total_n_before else 0 + log_step(f"# Sort + Dedup: done, removed {n_removed} duplicate rows total ({pct_removed:.2f}% of {total_n_before})") + + duplicates_output_file = output_file.with_name(f"{output_file.stem}_duplicates.parquet") + log_step(f"# Sort + Dedup: writing duplicate rows to {duplicates_output_file}") + duplicate_bucket_files = sorted(tmp_dir_duplicates.glob("bucket_id__*.parquet")) + concatenate_parquet_files(duplicate_bucket_files, duplicates_output_file) df_pheno = pl.scan_csv( phenotype_file, @@ -91,20 +116,6 @@ def main( separator="\t", ).select("FINNGENID", "SEX") - df_concat = ( - pl.scan_parquet(bucket_files) - # Join SEX - .join( - df_pheno, - left_on="FINNGENID", - right_on="FINNGENID", - how="left", - maintain_order="left", - ) - .with_row_index(name="_rowid", offset=1) - ) - - print("# Sanitize text fields") # Unicode "SYMBOL FOR NEWLINE", displayed as: ␤ unicode_newline = "\u2424" # Unicode "SYMBOL FOR HORIZONTAL TABULATION", displayed as: ␉ @@ -120,42 +131,56 @@ def main( "load_id_pseudo", "file_name_pseudo", "laboratoriotutkimusoid", - "_rowid", + "ROWID", "_rowid_source", "SEX", ] - ( - df_concat.with_columns( - pl.selectors.exclude(*trusted_columns) - .str.replace_all(pattern="\r\n|\r|\n", value=unicode_newline) - .str.replace_all(pattern="\t", value=unicode_tab, literal=True) + + tmp_dir_final = tmp_dir / "final" + tmp_dir_final.mkdir() + + # Process one bucket at a time (rather than one join + with_row_index + + # sanitize + sink over the whole dataset) because maintain_order="left" on + # the join and the global counter in with_row_index both disable Polars' + # streaming engine, forcing full materialization of the joined/indexed/ + # sanitized frame in memory. Per-bucket, each chunk is small enough to + # materialize safely; ROWID stays globally monotonic via a running offset. + log_step("# Join SEX + Sanitize: start (per bucket)") + row_index_offset = 1 + final_bucket_files = [] + for bucket_id in range(partition_n_buckets): + bucket_name = f"bucket_id__{bucket_id}.parquet" + bucket_file = tmp_dir_sort_dedup / bucket_name + n_rows = bucket_row_counts[bucket_name] + log_step( + f"# Join SEX + Sanitize: bucket {bucket_id + 1}/{partition_n_buckets} " + f"({bucket_name}), {n_rows} rows, ROWID offset {row_index_offset}" ) - # Re-order column to be somewhat backward compatible with previous implementation - .select( - "_rowid", - "_rowid_source", - "FINNGENID", - "EVENT_AGE", - "APPROX_EVENT_DAY", - "TIME", - "laboratoriotutkimusnimike", - "paikallinentutkimusnimike_koodi", - "paikallinentutkimusnimike_selite", - "tutkimuskoodistonjarjestelma", - "tutkimusvastauksentila", - "tutkimustulosarvo", - "tutkimustulosyksikko", - "tuloksenpoikkeavuus", - "viitearvoryhma", - "viitevalialkuarvo", - "viitevalialkuyksikko", - "viitevaliloppuarvo", - "viitevaliloppuyksikko", - "tutkimustulosteksti", - "SEX", + output_bucket_file = tmp_dir_final / bucket_name + ( + pl.scan_parquet(bucket_file) + .join( + df_pheno, + left_on="FINNGENID", + right_on="FINNGENID", + how="left", + maintain_order="left", + ) + .with_row_index(name="ROWID", offset=row_index_offset) + .with_columns( + pl.selectors.exclude(*trusted_columns) + .str.replace_all(pattern="\r\n|\r|\n", value=unicode_newline) + .str.replace_all(pattern="\t", value=unicode_tab, literal=True) + ) + .sink_parquet(output_bucket_file) ) - .sink_parquet(output_file) - ) + final_bucket_files.append(output_bucket_file) + row_index_offset += n_rows + log_step("# Join SEX + Sanitize: done") + + log_step("# Concatenate final bucket files: start") + concatenate_parquet_files(final_bucket_files, output_file) + log_step("# Concatenate final bucket files: done") def init_cli(): @@ -196,27 +221,16 @@ def init_cli(): def consolidate_columns(assembled_file: Path, output_file: Path) -> Path: - """Remove unnecessary columns form the assembled file and rename the ones we will keep.""" + """Keep all main.* columns (stripped of prefix) plus freetext text field and _rowid_source.""" + schema = pl.scan_parquet(assembled_file).collect_schema() + + internal_columns = {"main._rowid", "main._filename"} rename_columns = { - "main.FINNGENID": "FINNGENID", - "main.EVENT_AGE": "EVENT_AGE", - "main.APPROX_EVENT_DAY": "APPROX_EVENT_DAY", - "main.TIME": "TIME", - "main.laboratoriotutkimusnimike": "laboratoriotutkimusnimike", - "main.paikallinentutkimusnimike_koodi": "paikallinentutkimusnimike_koodi", - "main.paikallinentutkimusnimike_selite": "paikallinentutkimusnimike_selite", - "main.tutkimuskoodistonjarjestelma": "tutkimuskoodistonjarjestelma", - "main.tutkimusvastauksentila": "tutkimusvastauksentila", - "main.tutkimustulosarvo": "tutkimustulosarvo", - "main.tutkimustulosyksikko": "tutkimustulosyksikko", - "main.tuloksenpoikkeavuus": "tuloksenpoikkeavuus", - "main.viitearvoryhma": "viitearvoryhma", - "main.viitevalialkuarvo": "viitevalialkuarvo", - "main.viitevalialkuyksikko": "viitevalialkuyksikko", - "main.viitevaliloppuarvo": "viitevaliloppuarvo", - "main.viitevaliloppuyksikko": "viitevaliloppuyksikko", - "freetext.tutkimustulosteksti": "tutkimustulosteksti", + col: col.removeprefix("main.") + for col in schema.names() + if col.startswith("main.") and col not in internal_columns } + rename_columns["freetext.tutkimustulosteksti"] = "tutkimustulosteksti" out_columns = list(rename_columns.keys()) + ["_rowid_source"] @@ -232,6 +246,7 @@ def consolidate_columns(assembled_file: Path, output_file: Path) -> Path: def partition(assembled_file: Path, tmp_dir: Path, n_buckets): for bucket_id in range(n_buckets): + log_step(f"# Partition: bucket {bucket_id + 1}/{n_buckets}") ( pl.scan_parquet(assembled_file) .filter(pl.col("FINNGENID").hash() % n_buckets == bucket_id) @@ -239,7 +254,29 @@ def partition(assembled_file: Path, tmp_dir: Path, n_buckets): ) +def concatenate_parquet_files(input_files: list[Path], output_file: Path) -> None: + """Concatenate parquet files by copying row batches directly, one at a time. + + Polars' `sink_parquet` over a multi-file `scan_parquet` is not reliably + streaming in all versions/engines, and has been observed to materialize + the whole dataset in memory. Using pyarrow directly guarantees a bounded + memory footprint regardless of total row count. + """ + writer = None + try: + for input_file in input_files: + reader = pq.ParquetFile(input_file) + if writer is None: + writer = pq.ParquetWriter(output_file, reader.schema_arrow, compression="zstd") + for batch in reader.iter_batches(): + writer.write_batch(batch) + finally: + if writer is not None: + writer.close() + + def sort_dedup(frame: pl.LazyFrame | pl.DataFrame): + """Sort by the full column order, then split into kept rows (first per duplicate key) and dropped duplicates.""" all_columns = frame.collect_schema().names() sort_subset_columns = set(COLUMNS_UNIQUENESS_SORT) other_columns = [] @@ -249,15 +286,19 @@ def sort_dedup(frame: pl.LazyFrame | pl.DataFrame): sort_full_columns = COLUMNS_UNIQUENESS_SORT + other_columns - return frame.sort(by=sort_full_columns).unique( - subset=COLUMNS_UNIQUENESS_SORT, keep="first", maintain_order=True - ) + sorted_frame = frame.sort(by=sort_full_columns) + is_first = pl.struct(COLUMNS_UNIQUENESS_SORT).is_first_distinct() + + kept = sorted_frame.filter(is_first) + dropped = sorted_frame.filter(~is_first) + return kept, dropped if __name__ == "__main__": args = init_cli() output.check_safe_write(args.output_file) + output.check_safe_write(args.output_file.with_name(f"{args.output_file.stem}_duplicates.parquet")) tmp_dir = output.create_tmp_dir() main( diff --git a/wdl/analysis/kanta_ca.json b/wdl/analysis/kanta_ca.json deleted file mode 100644 index 922e738..0000000 --- a/wdl/analysis/kanta_ca.json +++ /dev/null @@ -1,47 +0,0 @@ -{ - "kanta_ca.test": true, # runs only 10 regions - "kanta_ca.prefix": "kanta_test", - "kanta_ca.pheno_list_file": "gs://r13-data/kanta/cond/inputs/kanta_test_phenos.txt", - "kanta_ca.pheno_file": "gs://r13-data/kanta/cond/inputs/kanta_pheno_cov.txt.gz", - "kanta_ca.ss_root": "gs://r12-data/kanta-lab/release_v2/summary_stats/PHENO.gz", - "kanta_ca.covariates": ["SEX_IMPUTED","AGE_AT_DEATH_OR_END_OF_FOLLOWUP","PC1","PC2","PC3","PC4","PC5","PC6","PC7","PC8","PC9","PC10","IS_FINNGEN2_CHIP","BATCH_DS1_BOTNIA_Dgi_norm","BATCH_DS10_FINRISK_Palotie_norm","BATCH_DS11_FINRISK_PredictCVD_COROGENE_Tarto_norm","BATCH_DS12_FINRISK_Summit_norm","BATCH_DS13_FINRISK_Bf_norm","BATCH_DS14_GENERISK_norm","BATCH_DS15_H2000_Broad_norm","BATCH_DS16_H2000_Fimm_norm","BATCH_DS17_H2000_Genmets_norm_relift","BATCH_DS18_MIGRAINE_1_norm_relift","BATCH_DS19_MIGRAINE_2_norm","BATCH_DS2_BOTNIA_T2dgo_norm","BATCH_DS20_SUPER_1_norm_relift","BATCH_DS21_SUPER_2_norm_relift","BATCH_DS22_TWINS_1_norm","BATCH_DS23_TWINS_2_norm_nosymmetric","BATCH_DS24_SUPER_3_norm","BATCH_DS25_BOTNIA_Regeneron_norm","BATCH_DS26_DIREVA_norm","BATCH_DS27_NFBC66_norm","BATCH_DS28_NFBC86_norm","BATCH_DS3_COROGENE_Sanger_norm","BATCH_DS4_FINRISK_Corogene_norm","BATCH_DS5_FINRISK_Engage_norm","BATCH_DS6_FINRISK_FR02_Broad_norm_relift","BATCH_DS7_FINRISK_FR12_norm","BATCH_DS8_FINRISK_Finpcga_norm","BATCH_DS9_FINRISK_Mrpred_norm"], - - # DOCKER (Used for the task via the workflow) - "kanta_ca.extract_regions.docker": "eu.gcr.io/finngen-sandbox-v3-containers/finemap-suite:56205d3", - "kanta_ca.docker": "eu.gcr.io/finngen-refinery-dev/regenie:3.3_r12_cond", - - # CONDITIONAL PARAMS - "kanta_ca.regenie_conditional.regenie_params": " ' --qt --bsize 200 --ref-first ' ", - "kanta_ca.regenie_conditional.null_root": "gs://r12-data/kanta-lab/release_v2/locos/R12_GRM_V0_LD_0.2.PHENO.loco.gz", - "kanta_ca.regenie_conditional.bgen_root": "gs://finngen-production-library-red/finngen_R13/bgen_1.0/data/chrom/finngen_R13_CHROM.bgen", - "kanta_ca.regenie_conditional.conditioning_mlogp_threshold": 6, - "kanta_ca.regenie_conditional.max_steps": 10, - - # WORKFLOW-LEVEL SUMSTATS COLUMN DEFINITIONS - "kanta_ca.rsid_col": "", - "kanta_ca.chromosome_col": "#chrom", - "kanta_ca.position_col": "pos", - "kanta_ca.allele1_col": "ref", - "kanta_ca.allele2_col": "alt", - "kanta_ca.freq_col": "af_alt", - "kanta_ca.beta_col": "beta", - "kanta_ca.se_col": "sebeta", - "kanta_ca.p_col": "pval", - "kanta_ca.mlogp_col": "mlogp", - "kanta_ca.delimiter": "TAB", - - # TASK-LEVEL PARAMETERS (Still under kanta_ca.extract_regions) - "kanta_ca.extract_regions.window": 1500000, - "kanta_ca.extract_regions.max_region_width": -1, - "kanta_ca.extract_regions.window_shrink_ratio": 0.9, - "kanta_ca.extract_regions.p_threshold": 0.00000005, - "kanta_ca.extract_regions.exclude_MHC": false, - "kanta_ca.extract_regions.scale_se_by_pval": false, - "kanta_ca.extract_regions.x_chromosome": true, - "kanta_ca.extract_regions.set_variant_id": true, - "kanta_ca.extract_regions.set_variant_id_map_chr": "23=X", - "kanta_ca.extract_regions.mem": 32, - - "kanta_ca.filter_covariates.threshold_cov_count": 10 - -} diff --git a/wdl/analysis/kanta_ca.wdl b/wdl/analysis/kanta_ca.wdl deleted file mode 100644 index 0d32723..0000000 --- a/wdl/analysis/kanta_ca.wdl +++ /dev/null @@ -1,353 +0,0 @@ -version 1.0 - -workflow kanta_ca { - input { - String prefix - File pheno_list_file - File pheno_file - String ss_root - String docker - Boolean test - Array[String] covariates - - String rsid_col - String chromosome_col - String position_col - String allele1_col - String allele2_col - String freq_col - String beta_col - String se_col - String p_col - String mlogp_col - String delimiter - } - - call filter_covariates { - input: - pheno_file = pheno_file, - pheno_list = pheno_list_file, - covariates = covariates - } - - scatter (pheno in read_lines(pheno_list_file)) { - call extract_regions { - input: - sumstats = sub(ss_root, "PHENO", pheno), - pheno = pheno, - # Pass workflow inputs to task - rsid_col = rsid_col, - chromosome_col = chromosome_col, - position_col = position_col, - allele1_col = allele1_col, - allele2_col = allele2_col, - freq_col = freq_col, - beta_col = beta_col, - se_col = se_col, - p_col = p_col, - delimiter = delimiter - } - } - - call merge_regions { - input: - hits = extract_regions.region, - test = test - } - - Map[String, String] cov_map = read_map(filter_covariates.cov_pheno_map) - Array[Array[String]] all_regions = read_tsv(merge_regions.regions) - - scatter (region in all_regions) { - String pheno = region[0] - String chrom = region[1] - String region_limits = region[2] - String locus = region[3] - - call regenie_conditional { - input: - docker = docker, - prefix = prefix, - locus = locus, - region = region_limits, - pheno = pheno, - chrom = chrom, - covariates = cov_map[pheno], - mlogp_col = mlogp_col, - chr_col = chromosome_col, - pos_col = position_col, - ref_col = allele1_col, - alt_col = allele2_col, - sumstats_root = ss_root, - beta_col = beta_col, - se_col = se_col, - pheno_file = pheno_file - } - } - - output { - Array[File] results = flatten(regenie_conditional.conditional_chains) - } -} - -task regenie_conditional { - input { - # GENERAL PARAMS - String docker - String prefix - # hit info - String locus - String region - String pheno - String chrom - # files to localize - File pheno_file - String bgen_root - String null_root - String sumstats_root - # column names and stuff - String chr_col - String pos_col - String ref_col - String alt_col - String mlogp_col - String beta_col - String se_col - # Script parameters/options - Float conditioning_mlogp_threshold - Int max_steps - String covariates - String? regenie_params - } - - # localize all files based on roots and pheno/chrom info - File sumstats = sub(sumstats_root, "PHENO", pheno) - File sum_tabix = sumstats + ".tbi" - File null = sub(null_root, "PHENO", pheno) - File bgen = sub(bgen_root, 'CHROM', chrom) - File bgen_sample = bgen + ".sample" - File bgen_index = bgen + ".bgi" - - # runtime params based on file sizes - Int disk_size = 120 - - command <<< - tabix -h ~{sumstats} ~{region} > region_sumstats.txt - - python3 /scripts/regenie_conditional.py \ - --out ./~{prefix} \ - --bgen ~{bgen} \ - --null-file ~{null} \ - --sumstats region_sumstats.txt \ - --pheno-file ~{pheno_file} \ - --pheno ~{pheno} \ - --locus-region ~{locus} \ - ~{region} \ - --pval-threshold ~{conditioning_mlogp_threshold} \ - --max-steps ~{max_steps} \ - --chr-col '~{chr_col}' \ - --pos-col '~{pos_col}' \ - --ref-col '~{ref_col}' \ - --alt-col '~{alt_col}' \ - --mlogp-col '~{mlogp_col}' \ - --beta-col '~{beta_col}' \ - --sebeta-col '~{se_col}' \ - --covariates ~{covariates} \ - ~{if defined(regenie_params) then " --regenie-params " + regenie_params else ""} \ - --log info - >>> - - output { - Array[File] conditional_chains = glob("./${prefix}*.snps") - Array[File] logs = glob("./${prefix}*.log") - Array[File] regenie_output = glob("./${prefix}*.conditional") - } - - runtime { - docker: "~{docker}" - disks: "local-disk ${disk_size} HDD" - } -} - -task merge_regions { - input { - Array[File] hits - Boolean test - } - - String outfile = "regions.txt" - - command <<< - while read f; do cat $f >> tmp.txt; done < ~{write_lines(hits)} - cat tmp.txt ~{if test then " | shuf | head -n 10" else ""} > ~{outfile} - >>> - - output {File regions = outfile} -} - -task filter_covariates { - - input { - File pheno_file - Array[String] covariates - File pheno_list - Int threshold_cov_count - } - - String outfile = "./pheno_cov_map_" + threshold_cov_count + ".txt" - Int disk_size = ceil(size(pheno_file,'GB')) + 2 * 2 - - command <<< - - set -euxo pipefail - - python3 <0).sum().to_frame(pheno) - df = pd.concat([df,m],axis =1) - - print(f"{i+1}/{len(phenos_groups)} {pheno_name}") - #If it's a group of phenos the min function will return the lowest count across all phenos - tmp_df = df[pheno_list].min(axis =1) - covs = tmp_df.index[tmp_df >= ~{threshold_cov_count}].tolist() - missing_covs = [elem for elem in covariates if elem not in covs] - if missing_covs:tmp_err.write(f"{pheno_name}\t{','.join(missing_covs)}\n") - o.write(f"{pheno_name}\t{','.join(covs)}\n") - - CODE - - >>> - output { - File cov_pheno_map = outfile - File cov_pheno_warning = sub(outfile,'.txt','.err.txt') - } - - runtime { - memory: "64 GB" - disks: "local-disk ${disk_size} HDD" - } - -} - - -task extract_regions { - input { - File sumstats - String pheno - # Passed from Workflow - String rsid_col - String chromosome_col - String position_col - String allele1_col - String allele2_col - String freq_col - String beta_col - String se_col - String p_col - String delimiter - - # Task-specific parameters remain - Int window - Int max_region_width - Float window_shrink_ratio - Float p_threshold - Float? minimum_pval - - Boolean scale_se_by_pval - Boolean exclude_MHC - Boolean x_chromosome - Boolean set_variant_id - String? set_variant_id_map_chr - - String docker - Int mem - } - - command <<< - set -euo pipefail - - make_finemap_inputs.py \ - --sumstats ~{sumstats} \ - --rsid-col '~{rsid_col}' \ - --chromosome-col '~{chromosome_col}' \ - --position-col '~{position_col}' \ - --allele1-col '~{allele1_col}' \ - --allele2-col '~{allele2_col}' \ - --freq-col '~{freq_col}' \ - --beta-col '~{beta_col}' \ - --se-col '~{se_col}' \ - --p-col '~{p_col}' \ - --delimiter '~{delimiter}' \ - --grch38 \ - ~{true='--exclude-MHC ' false='' exclude_MHC} \ - --prefix ~{pheno} \ - --out ~{pheno} \ - --window ~{window} \ - ~{if (max_region_width < 0) then "" else "--max-region-width " + max_region_width } \ - --window-shrink-ratio ~{window_shrink_ratio} \ - ~{true='--scale-se-by-pval ' false='' scale_se_by_pval} \ - ~{true='--x-chromosome ' false='' x_chromosome} \ - ~{true='--set-variant-id ' false='' set_variant_id} \ - ~{if defined(set_variant_id_map_chr) then '--set-variant-id-map-chr ' + set_variant_id_map_chr else ''} \ - --p-threshold ~{p_threshold} \ - ~{if defined(minimum_pval) then '--min-p-threshold ' + minimum_pval else ''} \ - --wdl - - # Define the res variable by reading the file created by make_finemap_inputs.py - res=$(cat ~{pheno}_had_results) - - # Handle case where no results were found (res == "False") - if [ "$res" == "False" ]; then - touch ~{pheno}.lead_snps.txt - touch ~{pheno}.bed - touch ~{pheno}.formatted.txt - else - # Existing logic for successful runs - python3 -c "import pandas as pd; snps = pd.read_csv('~{pheno}.lead_snps.txt', sep='\t'); bed = pd.read_csv('~{pheno}.bed', sep='\t', names=['chr', 'start', 'end']); bed['chr'] = bed['chr'].astype(str); snps['chromosome'] = snps['chromosome'].astype(str); snps_with_regions = snps.apply(lambda row: pd.Series({'region': '{}:{}-{}'.format(row['chromosome'], m.iloc[0]['start'], m.iloc[0]['end']), 'chr': row['chromosome'], 'start': m.iloc[0]['start'], 'end': m.iloc[0]['end'], 'rsid': row['rsid'], 'mlogp': row['mlogp']}) if len(m := bed[(bed['chr'] == row['chromosome']) & (bed['start'] <= row['position']) & (bed['end'] >= row['position'])]) > 0 else pd.Series({'region': None, 'chr': None, 'start': None, 'end': None, 'rsid': None, 'mlogp': None}), axis=1).dropna(); top_snps = snps_with_regions.loc[snps_with_regions.groupby('region')['mlogp'].idxmax()]; top_snps['chr_num'] = pd.to_numeric(top_snps['chr'], errors='coerce'); top_snps = top_snps.sort_values(['chr_num', 'start']); open('~{pheno}.formatted.txt', 'w').write('\n'.join(['~{pheno}\t{}\t{}\t{}\t{:.5g}'.format(row['chr'], row['region'], row['rsid'], row['mlogp']) for _, row in top_snps.iterrows()]) + '\n')" - fi - >>> - - output { - File region = pheno + ".formatted.txt" - File leadsnps = pheno + ".lead_snps.txt" - File bed = pheno + ".bed" - File log = pheno + ".log" - Boolean had_results = read_boolean("~{pheno}_had_results") - } - - runtime { - docker: "~{docker}" - memory: "${mem} GB" - disks: "local-disk 20 HDD" - preemptible: 2 - noAddress: true - } -} diff --git a/wdl/analysis/kanta_regenie.json b/wdl/analysis/kanta_regenie.json deleted file mode 100644 index 9602940..0000000 --- a/wdl/analysis/kanta_regenie.json +++ /dev/null @@ -1,27 +0,0 @@ -{ "kanta_regenie.docker": "eu.gcr.io/finngen-refinery-dev/regenie:3.3_r12", - #"kanta_regenie.kanta_parquet": "gs://fg-3/kanta_v3/dev/release/kanta_dev_core.parquet", - "kanta_regenie.kanta_parquet": "gs://fg-3/kanta_v2/core/release/data/finngen_R13_kanta_lab_2.0.parquet", - "kanta_regenie.measurement_col_name": "MEASUREMENT_VALUE_MERGED", - "kanta_regenie.cov_file": "gs://r12-data/pheno/R12_COV_V2.FID.txt.gz", - "kanta_regenie.is_binary": false, - "kanta_regenie.bgen_list": "gs://fg-3/kanta_v2/inputs/bgen_chunks.txt", - - # PHENO INFO - "kanta_regenie.create_pheno_file.docker": "eu.gcr.io/finngen-refinery-dev/kanta:kanta_analysis_v2.2", - "kanta_regenie.pheno_omop": '3009542', - "kanta_regenie.create_pheno_file.min_count": 1, - #"kanta_regenie.validate_hits": "gs://fg-3/kanta_v2/inputs/3026361_hits.txt", - - #STEP 2 - "kanta_regenie.step2.test": "additive", - "kanta_regenie.step2.options": "--firth --approx --pThresh 0.01 --firth-se", - "kanta_regenie.step2.bsize": 400, - - # STEP1 - "kanta_regenie.step1.covariates": "SEX_IMPUTED,EVENT_AGE,PC{1:10},IS_FINNGEN2_CHIP,BATCH_DS1_BOTNIA_Dgi_norm,BATCH_DS10_FINRISK_Palotie_norm,BATCH_DS11_FINRISK_PredictCVD_COROGENE_Tarto_norm,BATCH_DS12_FINRISK_Summit_norm,BATCH_DS13_FINRISK_Bf_norm,BATCH_DS14_GENERISK_norm,BATCH_DS15_H2000_Broad_norm,BATCH_DS16_H2000_Fimm_norm,BATCH_DS17_H2000_Genmets_norm_relift,BATCH_DS18_MIGRAINE_1_norm_relift,BATCH_DS19_MIGRAINE_2_norm,BATCH_DS2_BOTNIA_T2dgo_norm,BATCH_DS20_SUPER_1_norm_relift,BATCH_DS21_SUPER_2_norm_relift,BATCH_DS22_TWINS_1_norm,BATCH_DS23_TWINS_2_norm_nosymmetric,BATCH_DS24_SUPER_3_norm,BATCH_DS25_BOTNIA_Regeneron_norm,BATCH_DS26_DIREVA_norm,BATCH_DS27_NFBC66_norm,BATCH_DS28_NFBC86_norm,BATCH_DS3_COROGENE_Sanger_norm,BATCH_DS4_FINRISK_Corogene_norm,BATCH_DS5_FINRISK_Engage_norm,BATCH_DS6_FINRISK_FR02_Broad_norm_relift,BATCH_DS7_FINRISK_FR12_norm,BATCH_DS8_FINRISK_Finpcga_norm,BATCH_DS9_FINRISK_Mrpred_norm", - "kanta_regenie.step1.grm_bed": "gs://r12-data/grm/R12_GRM_V0_LD_0.2.bed", - "kanta_regenie.step1.bsize": 1000, - "kanta_regenie.step1.options": "", - "kanta_regenie.step1.covariate_inclusion_threshold": 10, -} - diff --git a/wdl/analysis/kanta_regenie.wdl b/wdl/analysis/kanta_regenie.wdl deleted file mode 100644 index 56156a8..0000000 --- a/wdl/analysis/kanta_regenie.wdl +++ /dev/null @@ -1,370 +0,0 @@ -version 1.0 - -workflow kanta_regenie { - input { - String docker - File kanta_parquet - String pheno_omop - File cov_file - Boolean is_binary - String measurement_col_name - File bgen_list - File? validate_hits - } - - # EXTRACTS FROM HARMONIZED KANTA DATA THE ENTRIES WITH MATCHING OMOP - # CALCULATES AVERAGE AGE AND MEASUREMENT_UNIT AND MERGES WITH COV FILE - call create_pheno_file { - input : - pheno_omop = pheno_omop, - cov_file=cov_file, - parquet_file = kanta_parquet, - mes_col = measurement_col_name, - } - - # REGENIE STEP1 - call step1 { - input : - docker = docker, - phenolist=[pheno_omop], - is_binary = is_binary, - cov_pheno = create_pheno_file.pheno_file, - } - - # REGENIE STEP2 - scatter (bgen in read_lines(bgen_list)) { - call step2 { - input: - docker = docker, - bgen = bgen, - validate_hits=validate_hits, - phenolist=[pheno_omop], - is_binary=is_binary, - cov_pheno=create_pheno_file.pheno_file, - covariates=step1.covars_used, - pred=step1.pred, - loco=step1.loco, - nulls=step1.nulls, - firth_list=step1.firth_list, - } - } - - Array[File] results = flatten(step2.regenie) - call gather { input : files = results } - -} - - -task gather { - input { - Array[File] files - } - command <<< - pheno=`basename ~{files[0]} .regenie.gz | awk -F "." '{sub(/[^_]*_/, "", $NF); print $NF}'` - mkdir regenie munged - - echo -e "`date`\tconcatenating result pieces into regenie/$pheno.regenie.gz, sorting by chr pos just in case" - cat <(zcat ~{files[0]} | head -1) <(for file in ~{sep=" " files};do zcat $file | tail -n+2 ;done | sort -k1,1g -k2,2g) | bgzip > regenie/$pheno.regenie.gz - echo -e "`date`\tconverting to munged/$pheno.gz to a format used for importing to pheweb" - zcat regenie/$pheno.regenie.gz | awk 'BEGIN {FS=" "; OFS="\t"; split("CHROM GENPOS ALLELE0 ALLELE1 LOG10P BETA SE A1FREQ", REQUIRED_FIELDS)} NR==1 {for(i=1;i<=NF;i++) h[$i]=i;for(i in REQUIRED_FIELDS) if (!(REQUIRED_FIELDS[i] in h)) {print REQUIRED_FIELDS[i]" expected in regenie header">>"/dev/stderr"; exit 1} print "#chrom","pos","ref","alt","pval","mlogp","beta","sebeta","af_alt"} NR>1 {print $h["CHROM"],$h["GENPOS"],$h["ALLELE0"],$h["ALLELE1"],10^-$h["LOG10P"],$h["LOG10P"],$h["BETA"],$h["SE"],$h["A1FREQ"]}' | bgzip > munged/$pheno.gz - - echo -e "`date`\ttabixing munged/$pheno.gz" - tabix -s1 -b2 -e2 munged/$pheno.gz - echo -e "`date`\tdone" - >>> - output { - File regenie = glob("regenie/*.regenie.gz")[0] - File pheweb = glob("munged/*.gz")[0] - File pheweb_tbi = glob("munged/*.gz.tbi")[0] - } - - runtime { - memory: "8 GB" - disks: "local-disk 200 HDD" - preemptible: 2 - noAddress: true - } -} - -task step2 { - - input { - Array[String] phenolist - File cov_pheno - String covariates - String test - Boolean is_binary - File bgen - File? validate_hits - File pred - File firth_list - Array[File] loco - Array[File] nulls - Int bsize - String options - String docker - } - String test_cmd = if test == "additive" then "" else "--test "+ test - File bgi = bgen + ".bgi" - File sample = bgen + ".sample" - String prefix = sub(basename(pred), "_pred.list", "") + "." + basename(bgen) - - command <<< - ## continue statement exits with 1.... switch to if statement below in case want to pipefail back - ##set -euxo pipefail - n_cpu=`grep -c ^processor /proc/cpuinfo` - # move loco files to /cromwell_root as pred file paths point there - for file in ~{sep=" " loco}; do - mv $file . - done - - # move null files to /cromwell_root as firth_list file paths point there - for file in ~{sep=" " nulls}; do - mv $file . - done - - regenie \ - --step 2 \ - ~{test_cmd} \ - ~{if is_binary then "--bt --af-cc" else ""} \ - --bgen ~{bgen} \ - --ref-first \ - --sample ~{sample} \ - --covarFile ~{cov_pheno} \ - --covarColList ~{covariates} \ - --phenoFile ~{cov_pheno} \ - --phenoColList ~{sep="," phenolist} \ - --pred ~{pred} \ - ~{if is_binary then "--use-null-firth ~{firth_list}" else ""} \ - ~{if defined(validate_hits) then "--extract ~{validate_hits}" else ""} \ - --bsize ~{bsize} \ - --threads $n_cpu \ - --gz \ - --out ~{prefix} \ - ~{options} - - >>> - - output { - Array[File] log = glob("*.log") - Array[File] regenie = glob("*.regenie.gz") - } - - runtime { - docker: "~{docker}" - memory:"8 GB" - disks: "local-disk " + (ceil(size(bgen, "G")) + 5) + " HDD" - zones: "europe-west1-b europe-west1-c europe-west1-d" - preemptible: 2 - noAddress: true - } -} - - -task step1 { - input { - Array[String] phenolist - Boolean is_binary - File grm_bed - String prefix = basename(grm_bed, ".bed") - File cov_pheno - String covariates - Int bsize - String options - String docker - Int covariate_inclusion_threshold - } - File grm_bim = sub(grm_bed, ".bed", ".bim") - File grm_fam = sub(grm_bed, ".bed", ".fam") - - command <<< - set -eux - n_cpu=`grep -c ^processor /proc/cpuinfo` - - #filter out covariates with too few observations - covars=~{covariates} - COVARFILE=~{cov_pheno} - PHENO="~{sep=',' phenolist}" - THRESHOLD=~{covariate_inclusion_threshold} - # Filter binary covariates that don't have enough covariate values in them - # Inputs: covariate file, comma-separated phenolist, comma-separated covariate list, threshold for excluding a covariate - # If a covariate is quantitative (=having values different from 0,1,NA), it is masked and will be passed through. - # If a binary covariate has value NA, it will not be counted towards 0 or 1 for that covariate. - # If a covariate does not seem to exist (e.g. PC{1:10}), it will be passed through. - # If any of the phenotypes is not NA on that row, that row will be counted. This is in line with the step1 mean-imputation for multiple phenotypes. - zcat -f $COVARFILE | awk -v covariates=$covars -v phenos=$PHENO -v th=$THRESHOLD > new_covars ' - BEGIN{FS="\t"} - NR == 1 { - covlen = split(covariates,covars,",") - phlen = split(phenos,phenoarr,",") - for (i=1; i<=NF; i++){ - h[$i] = i - mask[$i] = 0 - } - } - NR > 1 { - #if any of the phenotypes is not NA, then take the row into account - process_row=0 - for (ph in phenoarr){ - if ($h[phenoarr[ph]] != "NA"){ - process_row = 1 - } - } - if (process_row == 1){ - for (co in covars){ - if($h[covars[co]] == 0) { - zerovals[covars[co]] +=1 - } - else if($h[covars[co]] == 1) { - onevals[covars[co]] +=1 - } - else if($h[covars[co]] == "NA"){ - #no-op - na=0; - } - else { - #mask this covariate to be included, no matter the counts - #includes both covariate not found in header and quantitative covars - mask[covars[co]] =1 - } - } - } - - } - END{ - SEP="" - for (co in covars){ - if( ( zerovals[covars[co]] > th && onevals[covars[co]] > th ) || mask[covars[co]] == 1 ){ - printf("%s%s" ,SEP,covars[co]) - SEP="," - } - printf "Covariate %s zero count: %d one count: %d mask: %d\n",covars[co],zerovals[covars[co]],onevals[covars[co]],mask[covars[co]] >> "/dev/stderr"; - } - }' - - NEWCOVARS=$(cat new_covars) - # fid needs to be the same as iid in fam - awk '{$1=$2} 1' ~{grm_fam} > tempfam && mv tempfam ~{grm_fam} - - regenie \ - --step 1 \ - ~{if is_binary then "--bt" else ""} \ - --bed ~{sub(grm_bed, "\\.bed$", "")} \ - --covarFile ~{cov_pheno} \ - --covarColList $NEWCOVARS \ - --phenoFile ~{cov_pheno} \ - --phenoColList ~{sep="," phenolist} \ - --bsize ~{bsize} \ - --lowmem \ - --lowmem-prefix tmp_rg \ - --gz \ - --threads $n_cpu \ - --out ~{prefix} \ - ~{if is_binary then "--write-null-firth" else ""} \ - ~{options} - - # rename loco files with phenotype names and update pred.list accordingly giving it a unique name - awk '{orig=$2; sub(/_[0-9]+.loco.gz/, "."$1".loco.gz", $2); print "mv "orig" "$2} ' ~{prefix}_pred.list | bash - phenohash=`echo ~{sep="," phenolist} | md5sum | awk '{print $1}'` - awk '{sub(/_[0-9]+.loco.gz/, "."$1".loco.gz", $2)} 1' ~{prefix}_pred.list > ~{prefix}.$phenohash.pred.list - loco_n=$(wc -l ~{prefix}.$phenohash.pred.list|cut -d " " -f 1) - - #check that loco predictions were created for every pheno - if [[ $loco_n -ne ~{length(phenolist)} ]]; then - echo "The model did not converge. This job will abort." - exit 1 - fi - - if [[ "~{is_binary}" == "true" ]] - then - # rename firth files with phenotype names and update firth.list accordingly giving it a unique name - awk '{orig=$2; sub(/_[0-9]+.firth.gz/, "."$1".firth.gz", $2); print "mv "orig" "$2} ' ~{prefix}_firth.list | bash - awk '{sub(/_[0-9]+.firth.gz/, "."$1".firth.gz", $2)} 1' ~{prefix}_firth.list > ~{prefix}.$phenohash.firth.list - - #check if there is a firth approx per every null - firth_n=$(wc -l ~{prefix}.$phenohash.firth.list|cut -d " " -f 1) - if [[ $loco_n -ne $firth_n ]]; then - echo "fitting firth null approximations FAILED. This job will abort." - exit 1 - fi - else # touch files to have quant phenos not fail output globbing - touch ~{prefix}.$phenohash.firth.list - touch get_globbed.firth.gz - fi - >>> - runtime { - docker: "~{docker}" - disks: "local-disk 200 HDD" - memory:"16 GB" - zones: "europe-west1-b europe-west1-c europe-west1-d" - preemptible: 2 - noAddress: true - } - output { - File log = prefix + ".log" - Array[File] loco = glob("*.loco.gz") - File pred = glob("*.pred.list")[0] - Array[File] nulls = glob("*.firth.gz") - File firth_list = glob("*.firth.list")[0] - String covars_used = read_string("new_covars") - File covariatelist = "new_covars" - } -} - -task create_pheno_file { - input { - String pheno_omop - File parquet_file - File cov_file - Int min_count - String mes_col - String docker - } - String out_file = pheno_omop + "_median_RINT.txt.gz" - command <<< - set -eux - PHENO="~{pheno_omop}" - PARQUET="~{parquet_file}" - MEASUREMENT_COLUMN="~{mes_col}" - # EXTRACT IDs and dump to tsv - clickhouse --query " - SELECT - FINNGENID, - median(EVENT_AGE) AS EVENT_AGE, - median($MEASUREMENT_COLUMN) AS median_value - FROM ( - SELECT - FINNGENID, - EVENT_AGE, - $MEASUREMENT_COLUMN, - FROM - file('$PARQUET', Parquet) - WHERE - OMOP_CONCEPT_ID = '$PHENO' - AND QC_PASS == 1 - AND $MEASUREMENT_COLUMN IS NOT NULL - - ) - GROUP BY FINNGENID - HAVING count(*) >= ~{min_count} - ORDER BY FINNGENID ASC - INTO OUTFILE 'median.tsv' - FORMAT TSVWithNames; - " - head median.tsv - # with pandas, calculate medians and RINT the median measurement column into PHENO column - python3 -c "import sys, pandas as pd, numpy as np, scipy.stats as stats; df = pd.read_csv('./median.tsv',sep='\t');print(df);new_col = '$PHENO'; col = 'median_value';print(new_col,col); df[new_col] = stats.norm.ppf((stats.rankdata(df[col]) - 0.375) / (len(df[col]) + 0.25)); df.to_csv('./transformed.tsv', sep='\t',index=False)" - head transformed.tsv - # sort cov file and merge - zcat -f ~{cov_file}| (sed -u 1q; sort) > cov.txt - join -t $'\t' --header cov.txt transformed.tsv | bgzip -c > ~{out_file} - >>> - runtime { - docker : "~{docker}" - disks: "local-disk ~{ceil(size(parquet_file,'GB')) + ceil(size(cov_file,'GB')) + 10} HDD" - } - output { - File pheno_file = out_file - } -} - diff --git a/wdl/google_inputs.json b/wdl/google_inputs.json index 41aebc0..6e07414 100644 --- a/wdl/google_inputs.json +++ b/wdl/google_inputs.json @@ -7,8 +7,8 @@ "default_runtime_attributes": { "docker": "eu.gcr.io/finngen-refinery-dev/bioinformatics:0.8", "zones": "europe-west1-b europe-west1-c europe-west1-d", - "cpu":4, - "disk":"2 HDD", + "cpu":2, + "disk":"20 HDD", "memory" : "4 GB", "preemptible": "1" } diff --git a/wdl/intake.json b/wdl/intake.json new file mode 100644 index 0000000..cc9c32e --- /dev/null +++ b/wdl/intake.json @@ -0,0 +1,9 @@ +{ + "intake.test": false, + "intake.source_list_file": "gs://fg-3/kanta_v4/inputs/kanta_inputs.txt", + "intake.phenotype_file": "gs://finngen-production-library-red/finngen_R14/phenotype_1.0/data/finngen_R14_minimum_1.0.txt.gz", + "intake.assemble_docker": "europe-west1-docker.pkg.dev/finngen-refinery-dev/fg-refinery-registry/kanta:v4-intake.assemble", + "intake.tidyup_docker": "europe-west1-docker.pkg.dev/finngen-refinery-dev/fg-refinery-registry/kanta:v4-intake.latest", + "intake.partition_n_buckets": 36, + "intake.prefix": "finngen_R14_kanta_laboratory_responses_internal_1.0" +} diff --git a/wdl/intake.wdl b/wdl/intake.wdl new file mode 100644 index 0000000..4e3f5a5 --- /dev/null +++ b/wdl/intake.wdl @@ -0,0 +1,133 @@ +version 1.0 + +workflow intake { + input { + File source_list_file + File phenotype_file + String assemble_docker + String tidyup_docker + Int partition_n_buckets = 24 + String prefix = "finngen_R14_kanta_laboratory_responses_internal_1.0" + Boolean test = false + } + + Array[Array[File]] all_source_pairs = read_tsv(source_list_file) + Array[Array[File]] source_pairs = if test then [all_source_pairs[0]] else all_source_pairs + + call assemble { + input: + source_pairs = source_pairs, + docker = assemble_docker, + prefix = prefix + } + + call tidyup { + input: + assembled_file = assemble.assembled_file, + phenotype_file = phenotype_file, + docker = tidyup_docker, + partition_n_buckets = partition_n_buckets, + prefix = prefix + } + + call export_tsv { + input: + docker = tidyup_docker, + tidied_parquet = tidyup.tidied_parquet, + prefix = prefix + } + + output { + File assembled = assemble.assembled_file + File tidied_parquet = tidyup.tidied_parquet + File tidied_tsv_gz = export_tsv.tidied_tsv_gz + File tidied_duplicates_parquet = tidyup.tidied_duplicates_parquet + } +} + + +task assemble { + input { + Array[Array[File]] source_pairs + String docker + String prefix + } + + command <<< + set -euxo pipefail + python3 -m kanta.intake.assemble \ + --source-list-file ~{write_tsv(source_pairs)} \ + --output-file ~{prefix}_assembled.parquet + >>> + + output { + File assembled_file = "~{prefix}_assembled.parquet" + } + + runtime { + docker: docker + disks: "local-disk ~{ceil(size(flatten(source_pairs), 'GB') * 3)} SSD" + } +} + + +task tidyup { + input { + File assembled_file + File phenotype_file + String docker + Int partition_n_buckets + String prefix + } + + command <<< + set -euxo pipefail + echo "cpus: $(nproc)" + python3 -m kanta.intake.tidyup \ + --assembled-file ~{assembled_file} \ + --phenotype-file ~{phenotype_file} \ + --output-file ~{prefix}.parquet \ + --partition-n-buckets ~{partition_n_buckets} + >>> + + output { + File tidied_parquet = "~{prefix}.parquet" + File tidied_duplicates_parquet = "~{prefix}_duplicates.parquet" + } + + runtime { + docker: docker + predefinedMachineType: "n2d-highcpu-32" + disks: "local-disk ~{ceil(size(assembled_file, 'GB')) * 3 + 1} SSD" + } +} + + +task export_tsv { + input { + File tidied_parquet + String prefix + String docker + } + + Int disk_size = ceil(size(tidied_parquet, 'GB') * 2 + 10) + String out = prefix + ".txt.gz" + command <<< + set -euxo pipefail + clickhouse --query "SELECT * FROM '~{tidied_parquet}'" \ + --format TSVWithNames \ + --max_threads "$(nproc)" \ + --input_format_parquet_preserve_order 1 \ + | pigz > ~{out} + >>> + + output { + File tidied_tsv_gz = out + } + + runtime { + docker: docker + disks: "local-disk ~{disk_size} HDD" + cpu: 16 + } +} diff --git a/wdl/kanta_core.json b/wdl/kanta_core.json deleted file mode 100644 index 5bace5e..0000000 --- a/wdl/kanta_core.json +++ /dev/null @@ -1,16 +0,0 @@ -{ - "kanta_core.test": false, - "kanta_core.prefix": "kanta_dev", - # DATA INPUTS - "kanta_core.kanta_munged_data":"gs://fg-3/kanta_v3/munged/kanta_v3_harmonized_2026_03_06.txt.gz", - "kanta_core.validate_outputs.inclusion_list": "gs://finngen-production-library-red/all_allowed_ids_to_sb/finngen_R14_finngenid_actual_inclusion_list.txt", - # RUNTIME - "kanta_core.kanta_docker":"europe-west1-docker.pkg.dev/finngen-refinery-dev/fg-refinery-registry/kanta:v3.core.qc_fix", - #"kanta_core.release_docker":"eu.gcr.io/finngen-sandbox-v3-containers/kanta:parquet", - #"kanta_core.analysis_docker":"europe-west1-docker.pkg.dev/finngen-refinery-dev/fg-refinery-registry/kanta:dev.analysis.3", - "kanta_core.split.n_chunks": 32, - "kanta_core.munge.cpus": 8, - #ANALYSIS - "kanta_core.compare_versions.old_parquet": "gs://fg-3/kanta_v3/dev/release/kanta_dev_core.parquet", - -} diff --git a/wdl/kanta_core.wdl b/wdl/kanta_core.wdl deleted file mode 100644 index 12b1654..0000000 --- a/wdl/kanta_core.wdl +++ /dev/null @@ -1,279 +0,0 @@ -version 1.0 - -workflow kanta_core { - input { - String prefix - String kanta_docker - String? release_docker - String? analysis_docker - # test mode will use only 100k lines and 4 cpus - Boolean test - File kanta_munged_data - } - - # splits input in chunks - call split { input:test = test,kanta_data = kanta_munged_data} - scatter (i in range(length(split.chunks))) { - call munge { - input: docker = kanta_docker, prefix = i,chunk = split.chunks[i] - } - } - String base_prefix = if test then prefix + "_test" else prefix - # merge chunks & logs - call merge { input: prefix = base_prefix,munged_chunks = munge.munged_chunk,docker=kanta_docker} - call merge_logs {input: prefix = base_prefix,logs = flatten(munge.logs)} - # build parquet and release file - call release { input: docker = select_first([release_docker, kanta_docker]), mem = if test then 4 else 64, prefix = prefix, munged_data = merge.merged_file} - call validate_outputs {input : parquet_file = release.core_files[1],docker=kanta_docker} - # CHECKS AND PLOTS - call compare_versions {input: new_parquet=release.core_files[1],docker=select_first([analysis_docker, kanta_docker]),prefix=prefix} - # builds updated pos/neg tables - call build_pos_tables{input:merged_file = merge.merged_file,docker=select_first([analysis_docker, kanta_docker])} - # checks extraction automatically - call qc_extracted {input: core_parquet=release.core_files[1],docker=select_first([analysis_docker, kanta_docker]),prefix=prefix} - -} - - -task qc_extracted { - input { - File core_parquet - String docker - String prefix - } - - String dist_summary = prefix + "_extraction_summary.tsv" - command <<< - # makes KS comparison between extracted and harmonized data for all OMOP IDs with both type of entries - python3 /qc_scripts/omop_extracted_dist.py --full --file_path ~{core_parquet} --summary-file ~{dist_summary} - >>> - output { - File summary = dist_summary - Array[File] plots = glob("./plots/*png") - } - runtime { - disks: "local-disk ~{2*ceil(size(core_parquet,'GB')) + 10} HDD" - docker : "~{docker}" - memory: "16 GB" - } -} - -task compare_versions { - input { - String docker - String prefix - File new_parquet - File old_parquet - } - command <<< - # extracs counts from new release - python3 /qc_scripts/counts.py ~{new_parquet} -c FINNGENID,MEASUREMENT_VALUE_HARMONIZED,MEASUREMENT_VALUE_EXTRACTED,TEST_OUTCOME,TEST_OUTCOME_TEXT_EXTRACTED,OUTCOME_POS_EXTRACTED -b QC_PASS -o . --prefix ~{prefix} - # extracs counts from previous release - python3 /qc_scripts/counts.py ~{old_parquet} -c FINNGENID,MEASUREMENT_VALUE_HARMONIZED,MEASUREMENT_VALUE_EXTRACTED,TEST_OUTCOME,TEST_OUTCOME_TEXT_EXTRACTED,OUTCOME_POS_EXTRACTED -b QC_PASS -o . --prefix old - # extracts omop names - python3 -c "import csv, sys; reader = csv.DictReader(sys.stdin); print('conceptId\tconceptName'); [print(f\"{row['conceptId']}\t{row['conceptName']}\") for row in reader]" < /finngen_qc/data/LABfi_ALL.usagi.csv | (sed -u 1q;sort -u) > omop_name_table.tsv - ls - # builds summary table and plots comparisons - python3 /qc_scripts/count_plot.py --new ~{prefix}_omop_analysis.tsv --old ./old_omop_analysis.tsv --names ./omop_name_table.tsv --new_suffix ~{prefix} --old_suffix ~{basename(old_parquet,'.parquet')} --out_tsv ~{prefix}_omop_comparison.tsv --out_img ~{prefix}_omop_comparison.png - ls - >>> - output { - File figure = "~{prefix}_omop_comparison.png" - File comparison = "~{prefix}_omop_comparison.tsv" - File summary_table = "./~{prefix}_omop_analysis.tsv" - File summary_md = "./~{prefix}_omop_analysis.md" - } - runtime { - disks: "local-disk ~{2*ceil(size(new_parquet,'GB')) + 10} HDD" - docker : "~{docker}" - } -} - -task build_pos_tables{ - input { - String docker - File merged_file - } - - command <<< - python3 /qc_scripts/extract_pos_counts.py ~{merged_file} --map /finngen_qc/data/LABfi_ALL.usagi.csv --pn_orig /core/data/negpos_mapping.tsv --plus_orig /core/data/kanta_plusplus_abnormality.tsv - >>> - output{ - File plus_summary = "./plusplus_summary.tsv" - File posneg_summary = "./pos_neg_summary.tsv" - Array[File] pasteable = glob("./*pasteable*") - } - - runtime { - disks: "local-disk ~{ceil(size(merged_file,'GB')) + 10} HDD" - docker : "~{docker}" - } - -} - -task validate_outputs { - input { - String docker - File parquet_file - File inclusion_list - } - command <<< - RELEASE_SAMPLES='./release_samples.txt' - INCLUSION_SAMPLES='./inclusion_samples.txt' - - # Check if RELEASE_SAMPLES exists - clickhouse --query="SELECT DISTINCT FINNGENID FROM file('~{parquet_file}', Parquet) ORDER BY FINNGENID " | sort > $RELEASE_SAMPLES - echo "N samples in release file: $(wc -l "$RELEASE_SAMPLES" | awk '{print $1}')" - - #Inclusion list - zcat -f ~{inclusion_list} | sort | uniq > "$INCLUSION_SAMPLES" - echo "N samples in inclusion file: $(wc -l "$INCLUSION_SAMPLES" | awk '{print $1}')" - - # Calculate EXTRA_SAMPLES regardless of file creation - comm -23 "$RELEASE_SAMPLES" "$INCLUSION_SAMPLES" > extra_samples.txt - EXTRA_SAMPLES=$(cat extra_samples.txt | wc -l) - echo "Release samples that are not in exclusion list: $EXTRA_SAMPLES" - - >>> - runtime { - disks: "local-disk ~{ceil(size(parquet_file,'GB')) + 10} HDD" - docker : "~{docker}" - } - output{ File extra_samples = "./extra_samples.txt"} -} - -task release { - input { - String docker - File munged_data - Int mem - String prefix - } - String core_prefix = prefix + "_core" - String meta_prefix = prefix + "_extended_columns" - command <<< - echo ~{mem} - set -euxo pipefail - awk '/^MemTotal:/{print $2/1024/1024}' /proc/meminfo - /usr/bin/time -v bash /sb_release/run.sh ~{munged_data} . ~{core_prefix} core 2> tmp.txt - cat tmp.txt && cat tmp.txt | awk '/Maximum resident set size/ {print "Max memory usage (GB):", $6/1024/1024}' - /usr/bin/time -v bash /sb_release/run.sh ~{munged_data} . ~{meta_prefix} metadata 2> tmp.txt - cat tmp.txt && cat tmp.txt | awk '/Maximum resident set size/ {print "Max memory usage (GB):", $6/1024/1024}' - - >>> - runtime { - docker : "~{docker}" - disks: "local-disk ~{ceil(size(munged_data,'GB')) * 4 + 10} HDD" - memory: "~{mem} GB" - cpu : mem/4 - } - output { - Array[File] core_files = ["~{core_prefix}.txt.gz","~{core_prefix}.parquet","~{core_prefix}.log","~{core_prefix}_schema.json"] - Array[File] meta_files = ["~{meta_prefix}.txt.gz","~{meta_prefix}.parquet","~{meta_prefix}.log","~{meta_prefix}_schema.json"] - } -} - -task merge { - input { - Array[File] munged_chunks - String prefix - String docker - } - String out_file = prefix + ".txt.gz" - String dup_file = prefix +"_duplicates.txt.gz" - command <<< - # write header to reports file - zcat ~{munged_chunks[0]} | head -n1 | bgzip -c > tmp.txt.gz - # merge files including reports - while read f; do echo $f && date +%Y-%m-%dT%H:%M:%S && zcat $f | sed -E 1d | bgzip -c >> tmp.txt.gz ; done < <(cat ~{write_lines(munged_chunks)} | sort -V ) - - python3 /core/duplicates.py --input tmp.txt.gz --prefix ~{prefix} - - >>> - runtime { - disks: "local-disk ~{ceil(size(munged_chunks,'GB')) * 4 + 10} HDD" - docker : "~{docker}" - } - output { - File merged_file = out_file - File duplicates = dup_file - } -} - -task merge_logs { - input { - Array[File] logs - String prefix - } - command <<< - # merge all warn,abbr,unit files - cat ~{write_lines(logs)} > logs.txt - # write headers - for f in {err,warn} ; do cat logs.txt | grep $f | head -n1 | xargs head -n1 > ~{prefix}"_"$f".txt"; done - for f in {err,warn,log} ;do while read i ;do cat $i | sed -E 1d >> ~{prefix}"_"$f".txt"; done < <(cat logs.txt | grep $f | sort -V);done - >>> - runtime {disks: "local-disk ~{ceil(size(logs,'GB')) * 4 + 10} HDD"} - output { - File out_log = "~{prefix}_log.txt" - File out_err = "~{prefix}_err.txt" - File out_warn = "~{prefix}_warn.txt" - } -} - -task munge { - input { - String docker - File chunk - String prefix - Int cpus - } - - command <<< - set -euxo pipefail - python3 /core/main.py --gz --mp --raw-data ~{chunk} --prefix ~{prefix} - ls - >>> - runtime { - docker : "~{docker}" - disks: "local-disk ~{ceil(size(chunk,'GB')) * 3 + 10} HDD" - mem: "~{cpus} GB" - cpu : "~{cpus}" - } - output { - File munged_chunk = "~{prefix}.txt.gz" - Array[File] logs = glob("./~{prefix}*txt") - } -} - -task split{ - input { - File kanta_data - Int n_chunks - Boolean test - } - - Int chunks = if test then 4 else n_chunks - command <<< - FILE=~{kanta_data} - TEST=~{if test then 1 else 0} - CHUNKS=$(( TEST ? 4 : ~{n_chunks} )) - LINES=$(( TEST ? 4000000 : 250000000 )) - - echo "$([[ $TEST -eq 1 ]] && echo TEST || echo FULL) MODE: CHUNKS=$CHUNKS, LINES=$LINES" - CHUNK_SIZE=$(( (LINES + CHUNKS - 1) / CHUNKS )) - echo "Splitting $LINES lines into $CHUNKS chunks of ~$CHUNK_SIZE each" - gzip -dc "$FILE" | head -n1 > header.txt && echo "Header saved to header.txt" - cmd="gzip -dc \"$FILE\" | tail -n +2" - [[ $TEST -eq 1 ]] && cmd+=" | head -n $LINES" - eval "$cmd" | split -l "$CHUNK_SIZE" -d --verbose --filter='{ cat header.txt; cat; } | bgzip -c > ${FILE}.gz' - kanta - >>> - - runtime { - disks: "local-disk ~{ceil(size(kanta_data,'GB')) * 10 + 20} HDD" - } - - output { - Array[File] chunks = glob("./kanta*gz") - File header = "header.txt" - } -} diff --git a/wdl/kanta_munge.json b/wdl/kanta_munge.json deleted file mode 100644 index 5dfcd60..0000000 --- a/wdl/kanta_munge.json +++ /dev/null @@ -1,11 +0,0 @@ -{ - "kanta_munge.test": false, - "kanta_munge.kanta_docker": "europe-west1-docker.pkg.dev/finngen-refinery-dev/fg-refinery-registry/kanta:v3.release", - #"kanta_munge.analysis_docker": "europe-west1-docker.pkg.dev/finngen-refinery-dev/fg-refinery-registry/kanta:dev.34", - "kanta_munge.kanta_data": "gs://fg-3/kanta_v3/pre-munged/finngen_R14_kanta_laboratory_responses_internal_1.0_2025_12_23_unique.tsv.gz", - "kanta_munge.prefix": "kanta_v3_harmonized", - "kanta_munge.munge.cpus": 8, - "kanta_munge.munge.harmonization_branch": "development", - "kanta_munge.split.n_chunks": 32, -} - diff --git a/wdl/kanta_munge.wdl b/wdl/kanta_munge.wdl deleted file mode 100644 index 54f2055..0000000 --- a/wdl/kanta_munge.wdl +++ /dev/null @@ -1,216 +0,0 @@ -version 1.0 - -workflow kanta_munge { - input { - File kanta_data - String kanta_docker - String? analysis_docker - String prefix - # test mode will use only 100k lines and 4 cpus - Boolean test - } - - # splits input in chunks - call split { - input: - test = test, - kanta_data = kanta_data - } - - # PROPER MUNGING ACTION - scatter (i in range(length(split.chunks))) { - call munge { - input: - docker = kanta_docker, - prefix = i, - chunk = split.chunks[i] - } - } - - call GetCurrentDate{} - # MERGE CHUNKS AND LOGS - String base_prefix = prefix + (if test then "_test" else "") + "_" + GetCurrentDate.date_string - call merge_logs { - input: - prefix = base_prefix, - logs = flatten(munge.logs) - } - call merge { - input: - docker = select_first([analysis_docker,kanta_docker]), - prefix = base_prefix, - munged_chunks = munge.munged_chunk - } - call qc { - input: - docker = select_first([analysis_docker,kanta_docker]), - merged_parquet=merge.merged_parquet, - prefix = base_prefix - } -} - - -task qc { - input { - File merged_parquet - String prefix - String docker - } - - String unmap = prefix+ "_unmapped_entries.txt" - String injection = prefix+ "_candidate_injections.txt" - String injection_issues = prefix+ "_injection_check.tsv" - - command <<< - # this steps checks the quality of the injected entries - time python3 /qc_scripts/injection_check.py ~{merged_parquet} -o ~{injection_issues} - # this step creates a candidate injection based on KS values for unharmonized data with source values - # it also returns the counts of TEST_NAME,UNIT(cleaned) that do not have a mapping - time python3 /qc_scripts/unharmonized.py ~{merged_parquet} --min_count 500 --ks-n 100000 -a ~{injection} -u ~{unmap} - # this step creates the table of most common unit per OMOP_ID - time python3 /qc_scripts/create_harmonization_table.py --input ~{merged_parquet} - #Abnormality estimates - time python3 /qc_scripts/abnormality.py --parquet_file ~{merged_parquet} --min-count 1000 - >>> - runtime { - disks: "local-disk ~{ceil(size(merged_parquet,'GB')) + 20} HDD" - docker : "~{docker}" - memory: "16 GB" - } - - output { - File umapped_entries = "~{unmap}" - File injection_candidates = "~{injection}" - File injection_mismathces = "~{injection_issues}" - File harmonization_counts = "harmonization_counts.tsv" - File harmonization_diffs = "harmonization_diffs.tsv" - File ab_table = "abnormality_estimation.table.tsv" - File ab_counts = "abnormality_estimation.txt" - - } -} - - -task merge { - input { - Array[File] munged_chunks - String prefix - String docker - } - String parquet_prefix = prefix + "_formatted" - String out_file = prefix +".txt.gz" - - command <<< - # Get the exact first line to use as a filter - FIRST_LINE=$(zcat ~{munged_chunks[0]} | head -n1) - # Use that string to delete duplicates in the stream - sort -V ~{write_lines(munged_chunks)} | xargs pigz -dc | sed "1b; /^$FIRST_LINE$/d" | pigz -c > ~{out_file} - bash /sb_release/run.sh ~{out_file} . ~{parquet_prefix} munged - >>> - runtime { - disks: "local-disk ~{ceil(size(munged_chunks,'GB')) * 4 + 10} HDD" - docker : "~{docker}" - memory: "16 GB" - - } - - output { - File merged = out_file - File merged_parquet = "~{parquet_prefix}.parquet" - } -} - -task merge_logs { - input { - Array[File] logs - String prefix - } - command <<< - # merge all warn,abbr,unit files - cat ~{write_lines(logs)} > logs.txt - # write headers - for f in {err,warn,abbr,unit} ; do cat logs.txt | grep $f | head -n1 | xargs head -n1 > ~{prefix}"_"$f".txt"; done - for f in {err,warn,abbr,unit,log} ;do while read i ;do cat $i | sed -E 1d >> ~{prefix}"_"$f".txt"; done < <(cat logs.txt | grep $f | sort -V);done - >>> - runtime { - disks: "local-disk ~{ceil(size(logs,'GB')) * 4 + 10} HDD" - } - - output { - File out_log = "~{prefix}_log.txt" - File out_err = "~{prefix}_err.txt" - File out_abbr = "~{prefix}_abbr.txt" - File out_unit = "~{prefix}_unit.txt" - File out_warn = "~{prefix}_warn.txt" - } -} - -task munge { - input { - String docker - File chunk - String prefix - Int cpus - String harmonization_branch - } - String out_chunk = "~{prefix}_munged.txt.gz" - command <<< - set -euxo pipefail - python3 /finngen_qc/main.py --gz --out . --raw-data ~{chunk} --log info --mp --prefix ~{prefix} --harmonization-gh-branch ~{harmonization_branch} - zcat ~{out_chunk} | wc -l - >>> - runtime { - docker : "~{docker}" - disks: "local-disk ~{ceil(size(chunk,'GB')) * 4 + 8} HDD" - mem: "~{cpus} GB" - cpu : "~{cpus}" - } - - output { - File munged_chunk = out_chunk - Array[File] logs = glob("./~{prefix}*txt") - Array[File] problematic = glob("./*duplicates*gz") - } -} - -task split{ - input { - File kanta_data - Int n_chunks - Boolean test - } - Int chunks = if test then 4 else n_chunks - command <<< - FILE=~{kanta_data} - TEST=~{if test then 1 else 0} - CHUNKS=$(( TEST ? 4 : ~{n_chunks} )) - LINES=$(( TEST ? 4000000 : 250000000 )) - - echo "$([[ $TEST -eq 1 ]] && echo TEST || echo FULL) MODE: CHUNKS=$CHUNKS, LINES=$LINES" - CHUNK_SIZE=$(( (LINES + CHUNKS - 1) / CHUNKS )) - echo "Splitting $LINES lines into $CHUNKS chunks of ~$CHUNK_SIZE each" - gzip -dc "$FILE" | head -n1 > header.txt && echo "Header saved to header.txt" - cmd="gzip -dc \"$FILE\" | tail -n +2" - [[ $TEST -eq 1 ]] && cmd+=" | head -n $LINES" - eval "$cmd" | split -l "$CHUNK_SIZE" -d --verbose --filter='{ cat header.txt; cat; } | bgzip -c > ${FILE}.gz' - kanta - - >>> - runtime {disks: "local-disk ~{ceil(size(kanta_data,'GB')) * 10 + 20} HDD"} - output { - Array[File] chunks = glob("./kanta*gz") - File header = "header.txt" - } -} - - -task GetCurrentDate { - command <<< - date +%Y_%m_%d | tr -d '\n' - >>> - output { - String date_string = read_string(stdout()) - } - meta { - volatile: true - } -} diff --git a/wdl/release/analysis_v2.json b/wdl/release/analysis_v2.json deleted file mode 100644 index 9724fdb..0000000 --- a/wdl/release/analysis_v2.json +++ /dev/null @@ -1,12 +0,0 @@ -{ - "kanta_analysis.test": true, - "kanta_analysis.prefix": "finngen_R13_kanta_analysis_1.0", - # DATA INPUTS - "kanta_analysis.kanta_munged_data": "gs://fg-3/kanta_v2/munged/kanta_2025_03_18_munged.txt.gz", - "kanta_analysis.validate_outputs.inclusion_list": "gs://finngen-production-library-red/all_allowed_ids_to_sb/finngen_R13_finngenid_actual_inclusion_list.txt", - # RUNTIME - "kanta_analysis.kanta_docker": "eu.gcr.io/finngen-refinery-dev/kanta:kanta_analysis_v2.2", - "kanta_analysis.split.n_chunks": 16, - "kanta_analysis.analysis.cpus": 8, - -} diff --git a/wdl/release/kanta_munge_v2.json b/wdl/release/kanta_munge_v2.json deleted file mode 100644 index a7986c3..0000000 --- a/wdl/release/kanta_munge_v2.json +++ /dev/null @@ -1,16 +0,0 @@ -{ - "kanta_munge.test": false, - "kanta_munge.prefix": "finngen_R13_kanta_lab_1.0", - "kanta_munge.kanta_data": "gs://fg-3/kanta_v2/pre-munged/finngen_R12_kanta_laboratory_responses_reports_internal_2.0_2024_11_20_unique.tsv.gz", - "kanta_munge.merge.reports_columns": ["STATEMENT_ID","STATEMENT_TEXT"], - # PHENO/SAMPLE DATA - "kanta_munge.sex_map.min_pheno": "gs://finngen-production-library-red/finngen_R13/phenotype_1.0/data/finngen_R13_minimum_1.0.txt.gz", - "kanta_munge.validate_outputs.inclusion_list": "gs://finngen-production-library-red/all_allowed_ids_to_sb/finngen_R13_finngenid_actual_inclusion_list.txt", - - # RUNTIME OPTIONS - "kanta_munge.kanta_docker": "eu.gcr.io/finngen-refinery-dev/kanta:kanta_analysis_v2.2", - "kanta_munge.munge.cpus": 8, - "kanta_munge.split.n_chunks": 32, - -} - diff --git a/wdl/v2/google_inputs.json b/wdl/v2/google_inputs.json deleted file mode 100644 index 41aebc0..0000000 --- a/wdl/v2/google_inputs.json +++ /dev/null @@ -1,15 +0,0 @@ -{ - "google_labels":{ - "project":"kanta", - "product":"core-analysis" - }, - - "default_runtime_attributes": { - "docker": "eu.gcr.io/finngen-refinery-dev/bioinformatics:0.8", - "zones": "europe-west1-b europe-west1-c europe-west1-d", - "cpu":4, - "disk":"2 HDD", - "memory" : "4 GB", - "preemptible": "1" - } -} diff --git a/wdl/v2/kanta_core.json b/wdl/v2/kanta_core.json deleted file mode 100644 index c0036be..0000000 --- a/wdl/v2/kanta_core.json +++ /dev/null @@ -1,14 +0,0 @@ -{ - "kanta_core.test": false, - "kanta_core.prefix": "kanta_dev", - # DATA INPUTS - "kanta_core.kanta_munged_data":"gs://fg-3/kanta_v3/dev/munged/kanta_dev_munged.txt.gz", - #"kanta_core.kanta_munged_data":"gs://fg-3/kanta_v2/munged/kanta_2025_04_30_munged.txt.gz", - - "kanta_core.validate_outputs.inclusion_list": "gs://finngen-production-library-red/all_allowed_ids_to_sb/finngen_R13_finngenid_actual_inclusion_list.txt", - # RUNTIME - "kanta_core.kanta_docker": "europe-west1-docker.pkg.dev/finngen-refinery-dev/fg-refinery-registry/kanta:v5.dev", - "kanta_core.split.n_chunks": 16, - "kanta_core.munge.cpus": 8, - -} diff --git a/wdl/v2/kanta_core.wdl b/wdl/v2/kanta_core.wdl deleted file mode 100644 index 3e33518..0000000 --- a/wdl/v2/kanta_core.wdl +++ /dev/null @@ -1,185 +0,0 @@ -version 1.0 - -workflow kanta_core { - input { - String prefix - String kanta_docker - # test mode will use only 100k lines and 4 cpus - Boolean test - File kanta_munged_data - } - - # splits input in chunks - call split { input:test = test,kanta_data = kanta_munged_data} - scatter (i in range(length(split.chunks))) { - call munge { - input: docker = kanta_docker, prefix = i,chunk = split.chunks[i] - } - } - String base_prefix = if test then prefix + "_test" else prefix - # merge chunks & logs - call merge { input: prefix = base_prefix,munged_chunks = munge.munged_chunk,docker=kanta_docker} - call merge_logs {input: prefix = base_prefix,logs = flatten(munge.logs)} - # build parquet and release file - call release { input: docker = kanta_docker, mem = if test then 4 else 64, prefix = prefix, munged_data = merge.merged_file} - call validate_outputs {input : parquet_file = release.core_files[1],docker=kanta_docker} -} - -task validate_outputs { - input { - String docker - File parquet_file - File inclusion_list - } - command <<< - RELEASE_SAMPLES='./release_samples.txt' - INCLUSION_SAMPLES='./inclusion_samples.txt' - - # Check if RELEASE_SAMPLES exists - clickhouse --query="SELECT DISTINCT FINNGENID FROM file('~{parquet_file}', Parquet) ORDER BY FINNGENID " | sort > $RELEASE_SAMPLES - echo "N samples in release file: $(wc -l "$RELEASE_SAMPLES" | awk '{print $1}')" - - #Inclusion list - zcat -f ~{inclusion_list} | sort | uniq > "$INCLUSION_SAMPLES" - echo "N samples in inclusion file: $(wc -l "$INCLUSION_SAMPLES" | awk '{print $1}')" - - # Calculate EXTRA_SAMPLES regardless of file creation - comm -23 "$RELEASE_SAMPLES" "$INCLUSION_SAMPLES" > extra_samples.txt - EXTRA_SAMPLES=$(cat extra_samples.txt | wc -l) - echo "Release samples that are not in exclusion list: $EXTRA_SAMPLES" - - >>> - runtime { - disks: "local-disk ~{ceil(size(parquet_file,'GB')) + 10} HDD" - docker : "~{docker}" - } - output{ File extra_samples = "./extra_samples.txt"} -} - -task release { - input { - String docker - File munged_data - Int mem - String prefix - } - String core_prefix = prefix + "_core" - String meta_prefix = prefix + "_metadata_columns" - command <<< - echo ~{mem} - set -euxo pipefail - awk '/^MemTotal:/{print $2/1024/1024}' /proc/meminfo - /usr/bin/time -v bash /sb_release/run.sh ~{munged_data} . ~{core_prefix} core 2> tmp.txt - cat tmp.txt && cat tmp.txt | awk '/Maximum resident set size/ {print "Max memory usage (GB):", $6/1024/1024}' - /usr/bin/time -v bash /sb_release/run.sh ~{munged_data} . ~{meta_prefix} metadata 2> tmp.txt - cat tmp.txt && cat tmp.txt | awk '/Maximum resident set size/ {print "Max memory usage (GB):", $6/1024/1024}' - - >>> - runtime { - docker : "~{docker}" - disks: "local-disk ~{ceil(size(munged_data,'GB')) * 4 + 10} HDD" - memory: "~{mem} GB" - cpu : mem/4 - } - output { - Array[File] core_files = ["~{core_prefix}.txt.gz","~{core_prefix}.parquet","~{core_prefix}.log","~{core_prefix}_schema.json"] - Array[File] meta_files = ["~{meta_prefix}.txt.gz","~{meta_prefix}.parquet","~{meta_prefix}.log","~{meta_prefix}_schema.json"] - } -} - -task merge { - input { - Array[File] munged_chunks - String prefix - String docker - } - String out_file = prefix + ".txt.gz" - String dup_file = prefix +"_duplicates.txt.gz" - command <<< - # write header to reports file - zcat ~{munged_chunks[0]} | head -n1 | bgzip -c > tmp.txt.gz - # merge files including reports - while read f; do echo $f && date +%Y-%m-%dT%H:%M:%S && zcat $f | sed -E 1d | bgzip -c >> tmp.txt.gz ; done < <(cat ~{write_lines(munged_chunks)} | sort -V ) - - python3 /core/duplicates.py --input tmp.txt.gz --prefix ~{prefix} - - >>> - runtime { - disks: "local-disk ~{ceil(size(munged_chunks,'GB')) * 4 + 10} HDD" - docker : "~{docker}" - } - output { - File merged_file = out_file - File duplicates = dup_file - } -} - -task merge_logs { - input { - Array[File] logs - String prefix - } - command <<< - # merge all warn,abbr,unit files - cat ~{write_lines(logs)} > logs.txt - # write headers - for f in {err,warn} ; do cat logs.txt | grep $f | head -n1 | xargs head -n1 > ~{prefix}"_"$f".txt"; done - for f in {err,warn,log} ;do while read i ;do cat $i | sed -E 1d >> ~{prefix}"_"$f".txt"; done < <(cat logs.txt | grep $f | sort -V);done - >>> - runtime {disks: "local-disk ~{ceil(size(logs,'GB')) * 4 + 10} HDD"} - output { - File out_log = "~{prefix}_log.txt" - File out_err = "~{prefix}_err.txt" - File out_warn = "~{prefix}_warn.txt" - } -} - -task munge { - input { - String docker - File chunk - String prefix - Int cpus - } - - command <<< - set -euxo pipefail - python3 /core/main.py --gz --mp --raw-data ~{chunk} --prefix ~{prefix} - ls - >>> - runtime { - docker : "~{docker}" - disks: "local-disk ~{ceil(size(chunk,'GB')) * 3 + 10} HDD" - mem: "~{cpus} GB" - cpu : "~{cpus}" - } - output { - File munged_chunk = "~{prefix}.txt.gz" - Array[File] logs = glob("./~{prefix}*txt") - } -} - -task split{ - input { - File kanta_data - Int n_chunks - Boolean test - } - - Int chunks = if test then 4 else n_chunks - command <<< - zcat ~{kanta_data} | head -n1 > header.txt - zcat ~{kanta_data} | sed -E 1d ~{if test then " | head -n 10000 " else ""} > tmp.tsv - for f in {00..~{chunks-1}}; do cat header.txt | bgzip -c > kanta$f.gz; done - split tmp.tsv -n l/~{chunks} -d kanta --filter='gzip >> $FILE.gz' - >>> - - runtime { - disks: "local-disk ~{ceil(size(kanta_data,'GB')) * 10 + 20} HDD" - } - - output { - Array[File] chunks = glob("./kanta*gz") - File header = "header.txt" - } -} diff --git a/wdl/v2/kanta_munge.json b/wdl/v2/kanta_munge.json deleted file mode 100644 index 13b23c0..0000000 --- a/wdl/v2/kanta_munge.json +++ /dev/null @@ -1,15 +0,0 @@ -{ - "kanta_munge.test": false, - "kanta_munge.kanta_docker": "europe-west1-docker.pkg.dev/finngen-refinery-dev/fg-refinery-registry/kanta:v4.dev", - #"kanta_munge.kanta_data": "gs://fg-3/kanta_v2/pre-munged/finngen_R12_kanta_laboratory_responses_reports_internal_2.0_2025_04_24_unique.tsv.gz", - #"kanta_munge.prefix": "finngen_R13_kanta_lab_1.0.txt.gz", - "kanta_munge.munge.cpus": 8, - "kanta_munge.munge.harmonization_branch": "development", - - "kanta_munge.split.n_chunks": 32, - # TEST INPUTS - "kanta_munge.prefix": "_dev", - "kanta_munge.kanta_data": "gs://fg-3/kanta_v3/dev/finngen_R12_kanta_laboratory_responses_reports_internal_2.0_2025_11_21_unique.tsv.gz", - -} - diff --git a/wdl/v2/kanta_munge.wdl b/wdl/v2/kanta_munge.wdl deleted file mode 100644 index 73cb2fe..0000000 --- a/wdl/v2/kanta_munge.wdl +++ /dev/null @@ -1,142 +0,0 @@ -version 1.0 - -workflow kanta_munge { - input { - File kanta_data - String kanta_docker - String prefix - # test mode will use only 100k lines and 4 cpus - Boolean test - } - - # splits input in chunks - call split { - input: - test = test, - kanta_data = kanta_data - } - - # PROPER MUNGING ACTION - scatter (i in range(length(split.chunks))) { - call munge { - input: - docker = kanta_docker, - prefix = i, - chunk = split.chunks[i] - } - } - # MERGE CHUNKS AND LOGS - String base_prefix = "kanta" + if test then "_test" else prefix - call merge_logs { - input: - prefix = base_prefix, - logs = flatten(munge.logs) - } - call merge { - input: - docker = kanta_docker, - prefix = base_prefix, - munged_chunks = munge.munged_chunk - } -} - -task merge { - input { - Array[File] munged_chunks - String prefix - String docker - - } - String out_file = prefix +"_munged.txt.gz" - - command <<< - # write header - zcat ~{munged_chunks[0]} | head -n1 | bgzip -c > ~{out_file} - # merge files without headers - while read f; do echo $f && date +%Y-%m-%dT%H:%M:%S && zcat $f | sed -E 1d | bgzip -c >> ~{out_file} ; done < <(cat ~{write_lines(munged_chunks)} | sort -V ) - >>> - runtime { - disks: "local-disk ~{ceil(size(munged_chunks,'GB')) * 4 + 10} HDD" - docker : "~{docker}" - } - - output { - File munged = out_file - } -} - -task merge_logs { - input { - Array[File] logs - String prefix - } - command <<< - # merge all warn,abbr,unit files - cat ~{write_lines(logs)} > logs.txt - # write headers - for f in {err,warn,abbr,unit} ; do cat logs.txt | grep $f | head -n1 | xargs head -n1 > ~{prefix}"_"$f".txt"; done - for f in {err,warn,abbr,unit,log} ;do while read i ;do cat $i | sed -E 1d >> ~{prefix}"_"$f".txt"; done < <(cat logs.txt | grep $f | sort -V);done - >>> - runtime { - disks: "local-disk ~{ceil(size(logs,'GB')) * 4 + 10} HDD" - } - - output { - File out_log = "~{prefix}_log.txt" - File out_err = "~{prefix}_err.txt" - File out_abbr = "~{prefix}_abbr.txt" - File out_unit = "~{prefix}_unit.txt" - File out_warn = "~{prefix}_warn.txt" - } -} - -task munge { - input { - String docker - File chunk - String prefix - Int cpus - String harmonization_branch - } - String out_chunk = "~{prefix}_munged.txt.gz" - command <<< - set -euxo pipefail - python3 /finngen_qc/main.py --out . --raw-data ~{chunk} --log info --mp --harmonization --gz --prefix ~{prefix} --harmonization-gh-branch ~{harmonization_branch} - zcat ~{out_chunk} | wc -l - >>> - runtime { - docker : "~{docker}" - disks: "local-disk ~{ceil(size(chunk,'GB')) * 4 + 10} HDD" - mem: "~{cpus} GB" - cpu : "~{cpus}" - } - - output { - File munged_chunk = out_chunk - Array[File] logs = glob("./~{prefix}*txt") - Array[File] problematic = glob("./*duplicates*gz") - } -} - -task split{ - input { - File kanta_data - Int n_chunks - Boolean test - } - Int chunks = if test then 4 else n_chunks - command <<< - zcat ~{kanta_data} | head -n1 > header.txt - zcat ~{kanta_data} | sed -E 1d ~{if test then " | head -n 4000000 " else ""} > tmp.tsv - for f in {00..~{chunks-1}}; do cat header.txt | bgzip -c > kanta$f.gz; done - split tmp.tsv -n l/~{chunks} -d kanta --filter='gzip >> $FILE.gz' - >>> - runtime {disks: "local-disk ~{ceil(size(kanta_data,'GB')) * 10 + 20} HDD"} - output { - Array[File] chunks = glob("./kanta*gz") - File header = "header.txt" - } -} - - - diff --git a/wdl/v2/pre-merge.json b/wdl/v2/pre-merge.json deleted file mode 100644 index 91b72f2..0000000 --- a/wdl/v2/pre-merge.json +++ /dev/null @@ -1,10 +0,0 @@ -{ - "pre_merge.test":false, - "pre_merge.kanta_list": "gs://fg-3/kanta_v2/inputs/kanta_file_list.txt", - "pre_merge.prefix": "finngen_R12_kanta_laboratory_responses_reports_internal_VERSION.txt.gz", - "pre_merge.version": "2.0", - "pre_merge.reports": "gs://fg-3/kanta_v2/raw_thl/finngen_R12_kanta_laboratory_reports_internal_1.0.txt.gz", - "pre_merge.sort_merge_cols": "asiakirjaoid_pseudo,merkintaoid_pseudo,entryoid_pseudo,load_id_pseudo,file_name_pseudo", - "pre_merge.sort_report.out_cols": "antaja_organisaatioid,lausunnontilaid,lausuntoteksti", - -} diff --git a/wdl/v2/pre-merge.wdl b/wdl/v2/pre-merge.wdl deleted file mode 100644 index ad1d2ea..0000000 --- a/wdl/v2/pre-merge.wdl +++ /dev/null @@ -1,218 +0,0 @@ -version 1.0 - -workflow pre_merge { - input { - Boolean test - File kanta_list - File reports - String sort_merge_cols - String prefix - String version - } - ############# - #--REPORTS--# - ############# - # sort report file based on new multi column key - call sort_file as sort_report { - input: - input_file = reports, - sort_merge_cols = sort_merge_cols, - test = false # i want to use all reports - } - # merge the reports file where ids match so free text is just one row - call join_reports { input:reports= sort_report.sorted_file } - - ############### - #--RESPONSES--# - ############### - # merge the various chunks into single files for each year - call merge_responses {input:kanta_list=kanta_list} - scatter (responses_file in merge_responses.merged_responses) { - # sort it as in the other task - call sort_file as sort_responses { - input : - test = test, # work with 10k lines only - input_file= responses_file, - sort_merge_cols = sort_merge_cols, - out_cols = merge_responses.responses_cols - } - # join each response with the pre sorted reports file - call merge_reports_responses {input:reports =join_reports.joined_reports ,responses = sort_responses.sorted_file} - } - ############# - #--MERGING--# - ############# - call merge_files {input:rr_files = merge_reports_responses.merged_file,out_file = sub(prefix,"VERSION",if test then version +"_test" else version) } -} - -task merge_files { - input { - Array[File] rr_files - String out_file - } - command <<< - zcat ~{rr_files[0]} | head -n1 | bgzip -c > ~{out_file} - while read f; do echo $f && zcat $f | sed -E 1d | bgzip -c >> ~{out_file}; done < ~{write_lines(rr_files)} - zcat ~{out_file} | wc -l - >>> - runtime {disks: "local-disk ~{ceil(size(rr_files,'GB'))*3 + 10} HDD"} - output { File merged_file = out_file} -} - - -task join_reports{ - input { - File reports - } - String out_file = sub(basename(reports),".txt.gz","_joined_ids.txt.gz") - command <<< - # Process the file - zcat -f ~{reports} | awk ' - BEGIN { FS=OFS="\t"; seen_key = 0 } - { - if (NR == 1) { - # Initialize with first row - prev_key = $1 SUBSEP $2 SUBSEP $3 - prev_value = $4 - seen_key = 1 - next - } - current_key = $1 SUBSEP $2 SUBSEP $3 - if (current_key == prev_key) { - # Same key, append value - prev_value = prev_value "::" $4 - seen_key++ - } else { - # Different key, output previous key - split(prev_key, cols, SUBSEP) - prefix = seen_key > 1 ? "HUOM::" seen_key "::" : "" - print cols[1], cols[2], cols[3], prefix prev_value - # Reset for new key - prev_key = current_key - prev_value = $4 - seen_key=1 - } - } - END { - # Output last key - split(prev_key, cols, SUBSEP) - prefix = seen_key > 1 ? "HUOM::" seen_key "::" : "" - print cols[1], cols[2], cols[3], prefix prev_value - }' | bgzip -c > ~{out_file} - - zcat ~{reports} | wc -l - zcat ~{out_file} | wc -l - >>> - runtime { disks: "local-disk ~{ceil(size(reports,'GB')*3) + 10} HDD" } - output { - File joined_reports = out_file - } -} - -task merge_reports_responses { - input { - File reports - File responses - } - String out_file = sub(basename(responses),'responses','responses_reports') - command <<< - echo ~{out_file} - zcat ~{responses} | wc -l - join -t $'\t' -a 1 -o auto -e NA --header <(zcat ~{responses}) <(zcat ~{reports} ) | cut -f 2- | bgzip -c > ~{out_file} - wc -l ~{out_file} - >>> - runtime { - disks:"local-disk ~{ceil(size(responses,'GB') + size(reports,'GB')*3) + 2} HDD" - } - output { - File merged_file = out_file - } -} - -task sort_file { - input { - Boolean test - File input_file - String sort_merge_cols - String out_cols - } - String out_file = basename(input_file,".txt.gz") + "_sorted.txt.gz" - command <<< - echo ~{input_file} ~{out_file} - # MERGE COLS - echo ~{sort_merge_cols} - merge_cols=$(zcat -f ~{input_file} | head -n1 | tr '\t' '\n' | nl | grep -wf <(echo ~{sort_merge_cols} | tr ',' '\n') | awk '{print $1}' | tr '\n' ',' | rev | cut -c 2- | rev) - echo $merge_cols - # OUT_COLS - echo ~{out_cols} - out_cols=$(zcat -f ~{input_file} | head -n1 | tr '\t' '\n' | nl | grep -wf <(echo ~{out_cols} | tr ',' '\n') | awk '{print $1}' | tr '\n' ',' | rev | cut -c 2- | rev) - echo $out_cols - # OUT FILE - echo "MERGE COLUMNS" - CMD="paste <(zcat ~{input_file} | cut -f $merge_cols | tr '\t' '-') <(zcat ~{input_file} | cut -f $out_cols ) ~{if test then ' | head -n 10000 ' else ''} > tmp.txt" - echo $CMD - eval $CMD - echo "SORT ID COLUMNS" - head -n1 tmp.txt |bgzip -c > ~{out_file} - cat tmp.txt| sed -E 1d | sort | bgzip -c >> ~{out_file} - >>> - runtime { - disks: "local-disk ~{ceil(size(input_file,'GB'))*10 + 2} HDD" - } - output { - File sorted_file = out_file - } -} - -task merge_responses { - - input { - File kanta_list - } - - Array[File] kanta_files = read_lines(kanta_list) - Int disk = ceil(size(kanta_files,"GB"))*4 + 20 - - command <<< - echo "MERGE PART FILES" - merged_dir="merged" - mkdir $merged_dir - # merge partial gz files - cat ~{write_lines(kanta_files)} > file_list.txt - current=0 - total=$(wc -l < file_list.txt) - while read file; - do - ((current++)) - echo $current/$total $(du -sb $file | awk '{print $1/2^30"GB"}') $(basename $file) - # Check if the file has multiple parts - if [[ $file == *".part"* ]]; then - # Get the base filename without the part number - base_filename=$(echo "$file" | sed 's/\.part[0-9][0-9]*//') - # Concatenate the parts into a single file - cat "$file" >> "$merged_dir/$(basename "$base_filename")" - else - # Copy the single-part file to the merged directory - cp "$file" "$merged_dir" - fi - done < file_list.txt - - echo "CHECK SIZE" - for file in $merged_dir/*; - do - echo $(du -sb $file | awk '{print $1/2^30"GB"}') $(basename $file) - zcat $file | head -n1 | tr '\t' ',' > header.txt - done - cat header.txt - >>> - runtime { - disks: "local-disk ~{disk} HDD" - } - - output { - Array[File] merged_responses = glob("merged/*gz") - String responses_cols = read_string("header.txt") - } - -} - diff --git a/wdl/v2/sort_dup.json b/wdl/v2/sort_dup.json deleted file mode 100644 index 74b0c40..0000000 --- a/wdl/v2/sort_dup.json +++ /dev/null @@ -1,8 +0,0 @@ -{ - "kanta_sort_dup.test":false, - "kanta_sort_dup.kanta_data":"gs://fg-3/kanta_v2/finngen_R12_kanta_laboratory_responses_reports_internal_2.0.txt.gz", - "kanta_sort_dup.sex_map.min_pheno": "gs://finngen-production-library-red/finngen_R12/phenotype_1.0/data/finngen_R12_minimum_1.0.txt.gz", - "kanta_sort_dup.get_cols.branch": "dev", - "kanta_sort_dup.split.n_chunks": 15, #N.B it should be < 16 - -} diff --git a/wdl/v2/sort_dup.wdl b/wdl/v2/sort_dup.wdl deleted file mode 100644 index 0f09211..0000000 --- a/wdl/v2/sort_dup.wdl +++ /dev/null @@ -1,222 +0,0 @@ -version 1.0 - -workflow kanta_sort_dup{ - input { - # works with 100k lines - Boolean test - File kanta_data - } - - call get_cols {} # metadata - # split input in chunks - call split { - input: - test = test, - kanta_data = kanta_data, - cols = get_cols.cols, - s_cols = get_cols.s_cols - } - - # builds sex dictionary mapping from pheno file - call sex_map {} - - # extract columns sort and extract duplicates/errs - scatter (i in range(length(split.chunks))) { - call sort { - input : - index = i, - chunk = split.chunks[i], - sort_cols = split.sort_cols, - sex_map = sex_map.sex_map - - } - } - # merge chunks (unique/dup/err) - String prefix = basename(kanta_data,'.txt.gz') - call merge { - input : - sorted_chunks = sort.sorted_chunk, - sort_cols = split.sort_cols, - header = split.header, - prefix = if test then prefix+ "_test" else prefix - } -} - -task merge { - input { - File header - Array[File] sorted_chunks - Array[String] sort_cols - String prefix - } - - Int chunk_size = ceil(size(sorted_chunks,"GB")) - command <<< - # CONCAT PRE-SORTED FILES - echo "SORT FILES" - /usr/bin/time -v sort -t $'\t' -m -k ~{sep=" -k " sort_cols} ~{sep=" " sorted_chunks} > sorted.txt - # REMOVE DUPS - python3 <>> - runtime { - disks: "local-disk ~{chunk_size*4+10} HDD" - } - output { - Array[File] kanta_files = glob("~{prefix}*gz") - } -} - -task sort { - input { - File chunk - Array[String] sort_cols - Int index - File sex_map - } - String out_file = "kanta_sorted_" + index - command <<< - zcat ~{chunk} | sort -t $'\t' -k ~{sep=" -k " sort_cols} > tmp.txt - - #add sex - awk -F'\t' 'BEGIN {OFS="\t"} NR==FNR {sex[$1]=$2; next} NR==1 {print $0, "SEX"; next} {print $0, (sex[$1] ? sex[$1] : "NA")}' \ - ~{sex_map} tmp.txt > ~{out_file} - - # check file size - count_tmp=$(wc -l < tmp.txt) - count_out=$(wc -l < ~{out_file}) - - # Perform the assertion - if [[ "$count_tmp" -ne "$count_out" ]]; then - echo "❌ Assertion Failed: Line counts do not match!" >&2 - echo "tmp.txt has $count_tmp lines." >&2 - echo "~{out_file} has $count_out lines." >&2 - exit 1 # Exit with a non-zero status to signal an error - else - echo "✅ Assertion Passed: Both files have $count_tmp lines." - fi - >>> - - runtime { - disks: "local-disk ~{ceil(size(chunk,'GB'))*3 + 10} HDD" - } - - output { - File sorted_chunk = out_file - } -} - -task get_cols { - input {String branch} - - command <<< - # get required columns to cut from git repository - curl -s https://raw.githubusercontent.com/FINNGEN/kanta_lab_preprocessing/~{branch}/finngen_qc/magic_config.py > config.py - python3 -c "import config;o= open('./columns.txt','wt') ;o.write('\n'.join(list(config.config['rename_cols'].keys())) + '\n');o.write('\n'.join(config.config['other_cols'])+ '\n')" - python3 -c "import config;o= open('./sort_columns.txt','wt') ;o.write('\n'.join(config.config['sort_cols'])+ '\n')" - >>> - runtime { - disks: "local-disk 10 HDD" - } - output { - Array[String] cols = read_lines("columns.txt") - Array[String] s_cols = read_lines("sort_columns.txt") - } -} - -task split { - input { - Boolean test - File kanta_data - Int n_chunks - Array[String] cols - Array[String] s_cols - } - - Int disk_size = ceil(size(kanta_data,"GB"))*10*n_chunks - - command <<< - echo "SORT KANTA" - cat ~{write_lines(cols)} > columns.txt - cat ~{write_lines(s_cols)} > sort_columns.txt - COLS=$(zcat ~{kanta_data} | head -n1 | tr '\t' '\n' | grep -wnf columns.txt | cut -f 1 -d ':' | tr '\n' ',' | rev | cut -c2- | rev) - echo $COLS - - # uncompress and split new header from body - zcat ~{kanta_data} | cut -f $COLS | head -n1 > header.txt - zcat ~{kanta_data} | cut -f $COLS | sed -E 1d ~{if test then " | head -n 10000 " else ""}> tmp.tsv - - # GET SORT COLS AND KEEP ORDER - echo "COLS" - while read f; - do - cat header.txt | head -n1 | tr '\t' '\n'| grep -wn $f | cut -f 1 -d ':' >> sort_cols.txt - done < sort_columns.txt - cat sort_cols.txt - - # SPLIT INTO N FILES - split tmp.tsv -n l/~{n_chunks} -d kanta_chunk --filter='gzip > $FILE.gz' - >>> - - runtime { - disks: "local-disk ~{disk_size} HDD" - } - - output { - Array[File] chunks = glob("./kanta_chunk*gz") - File header = "header.txt" - Array[String] sort_cols = read_lines("sort_cols.txt") - } -} - -task sex_map { - input {File min_pheno} - String sex_file = "sex_map.txt" - command <<< - # get sex col - sexcol=$(awk '{for(i=1;i<=NF;i++){if($i=="SEX"){print i; exit}}}' <(zcat ~{min_pheno} | head -n1)) - # extract sex only and sort - zcat ~{min_pheno} | cut -f 1,$sexcol | (sed -u 1q ; sort )>> ~{sex_file} - >>> - runtime {disks: "local-disk ~{ceil(size(min_pheno,'GB')) * 3} HDD"} - output {File sex_map = sex_file} -}