Skip to content

mecfs_bio.build_system.task.sldsc.sldsc_regression_inputs

Turn gwaslab sumstats plus prepared S-LDSC references into kernel inputs, reproducing the preprocessing of ldsc's cell_type_specific as wrapped by gwaslab's estimate_h2_cts_by_ldsc: rows restricted to gwaslab's bundled HapMap3 list (rsID present, and EA/NEA equal to the HapMap3 A1/A2 in either orientation; the list depends on the sumstats' genome build), Z from BETA/SE, rows with missing SNP/Z/alleles dropped, an inner join onto the baseline SNP order (the prepared baseline is already restricted to SNPs with a weight LD score, matching ldsc's inner join with w_ld), N cast to int64 by truncation, and the chi-square filter Z^2 < max(0.001 * max(N), 80).

Fail-fast checks with no gwaslab counterpart: fewer than min_hapmap3_snps sumstats SNPs in the baseline's HapMap3 regression-SNP universe raises; a regression SNP with missing N raises (gwaslab fills with the median N); a SNP appearing more than once raises (gwaslab's merge would duplicate the row).

M_5_50 (passed through to the kernel with each LD-score column): for each annotation, the sum of its values over reference-panel SNPs with minor allele frequency above 5% (for a binary annotation, the number of common SNPs it contains). In the LDSC model E[chi2_j] = N_j * sum_c tau_c * l(j, c) + N_j * a + 1, with l(j, c) the LD score of SNP j with annotation c and tau_c the per-SNP heritability contributed by annotation c, M_5_50 converts per-SNP effects into the heritability attributed to an annotation (h2_c = tau_c * M_5_50_c).

Classes:

Functions:

AlignedSLDSCRegression

Attributes:

inputs instance-attribute 

inputs: SLDSCRegressionInputs

models instance-attribute 

models: tuple[SLDSCModel, ...]

build_regression_inputs 

build_regression_inputs(
    sumstats: DataFrame,
    baseline: PreparedBaseline,
    annotation_set: PreparedAnnotationSet,
    min_hapmap3_snps: int,
) -> AlignedSLDSCRegression
Source code in mecfs_bio/build_system/task/sldsc/sldsc_regression_inputs.py 
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
def build_regression_inputs(
    sumstats: pl.DataFrame,
    baseline: PreparedBaseline,
    annotation_set: PreparedAnnotationSet,
    min_hapmap3_snps: int,
) -> AlignedSLDSCRegression:
    joined = _join_baseline(sumstats=sumstats, baseline=baseline)
    _require_min_snps(joined=joined, min_hapmap3_snps=min_hapmap3_snps)
    regression = _filter_chisq(joined=_cast_n(joined=joined))
    baseline_rows = regression[_ROW_COL].to_numpy()
    annotation_rows = _annotation_rows(regression=regression, ld=annotation_set.ld)
    z = regression[GWASLAB_Z_COL].to_numpy()
    return AlignedSLDSCRegression(
        inputs=SLDSCRegressionInputs(
            chi2=z * z,
            n=regression[GWASLAB_SAMPLE_SIZE_COLUMN].to_numpy().astype(np.float64),
            w_ld=baseline.w_ld[W_LD_COL].to_numpy()[baseline_rows].astype(np.float64),
            baseline_ld=baseline.ld.matrix[baseline_rows].astype(np.float64),
            baseline_m=baseline.ld.m_5_50,
            annotation_ld=annotation_set.ld.matrix[annotation_rows],
            annotation_m=annotation_set.ld.m_5_50,
        ),
        models=tuple(
            SLDSCModel(
                label=m.label,
                target_index=annotation_set.ld.column_index(m.target_column),
                control_indices=tuple(
                    annotation_set.ld.column_index(c) for c in m.control_columns
                ),
            )
            for m in annotation_set.models
        ),
    )

load_sldsc_sumstats 

load_sldsc_sumstats(
    asset: Asset, pre_pipe: DataProcessingPipe
) -> pl.DataFrame

SNP, Z, N (Float64, possibly null) for every usable HapMap3 row of the sumstats.

Source code in mecfs_bio/build_system/task/sldsc/sldsc_regression_inputs.py 
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
def load_sldsc_sumstats(asset: Asset, pre_pipe: DataProcessingPipe) -> pl.DataFrame:
    """SNP, Z, N (Float64, possibly null) for every usable HapMap3 row of the sumstats."""
    sumstats = read_sumstats(asset)
    data = (
        pre_pipe.process(narwhals.from_native(sumstats.data).lazy())
        .collect()
        .to_polars()
    )
    required = (
        GWASLAB_RSID_COL,
        GWASLAB_SAMPLE_SIZE_COLUMN,
        GWASLAB_EFFECT_ALLELE_COL,
        GWASLAB_NON_EFFECT_ALLELE_COL,
    )
    missing = [c for c in required if c not in data.columns]
    assert not missing, f"sumstats need columns {missing} for S-LDSC"
    z = (
        pl.col(GWASLAB_Z_COL)
        if GWASLAB_Z_COL in data.columns
        else pl.col(GWASLAB_BETA_COL) / pl.col(GWASLAB_SE_COL)
    )
    selected = data.select(
        pl.col(GWASLAB_RSID_COL).cast(pl.String).alias(SNP_COL),
        z.cast(pl.Float64).alias(GWASLAB_Z_COL),
        pl.col(GWASLAB_SAMPLE_SIZE_COLUMN).cast(pl.Float64).fill_nan(None),
        pl.col(GWASLAB_EFFECT_ALLELE_COL).cast(pl.String),
        pl.col(GWASLAB_NON_EFFECT_ALLELE_COL).cast(pl.String),
    )
    return (
        _restrict_to_hapmap3(sumstats=selected, build=sumstats.build)
        .drop_nulls(
            subset=[
                SNP_COL,
                GWASLAB_Z_COL,
                GWASLAB_EFFECT_ALLELE_COL,
                GWASLAB_NON_EFFECT_ALLELE_COL,
            ]
        )
        .filter(pl.col(GWASLAB_Z_COL).is_not_nan())
        .select(SNP_COL, GWASLAB_Z_COL, GWASLAB_SAMPLE_SIZE_COLUMN)
    )