mecfs_bio.build_system.task.gwaslab.ldsc_diagnostic_plot_task

Task that draws the LD-score-regression diagnostic plot for a heritability analysis.

LD score regression fits E[chi^2_j] = intercept + (N * h2 / M) * ld_j but, because there are millions of variants, the fit is never shown -- standard regression diagnostics are essentially never performed (see the discussion in gwaslab_snp_heritability_by_ldsc_task). This task supplies one: it bins the analysis's variants into equal-count LD-score bins and plots each bin's mean chi-square against the fitted line, so model misfit (bins drifting off the line) and confounding (an intercept away from one, or a low-LD bin sitting below one) are visible by eye.

It is a sibling of SNPHeritabilityByLDSCTask rather than a reader of its output: it takes that task purely to reach a single, consistent definition of the analysis (which summary statistics, which LD reference, the preprocessing pipe, the sample size), and re-runs LD-score regression on that definition to draw the fitted line.

The fit is re-run with gwaslab -- the same estimator the heritability task uses -- so the line and the annotated intercept/heritability match what that task reports, rather than an independent solver that differs in its treatment of extreme chi-square (gwaslab keeps all variants and uses a two-step intercept estimator with a cutoff at chi-square 30; the repo's GenomicSEM port cuts variants above max(0.001*N, 80) and fits a single-step intercept). The re-run passes NaN prevalences so gwaslab reports observed-scale heritability, sidestepping the liability-scale rescaling it applies to binary traits, which would otherwise make the fitted slope hard to reconstruct. Consistently, the binned points keep every merged variant -- no high-chi-square cut -- to match the variant set gwaslab's slope regression uses.

The left y-axis is chi-square; the right y-axis is the same points rescaled to chi^2 * M / N, in which the fitted slope equals the observed-scale heritability and the axis is comparable across traits of different sample size. The two axes are one affine relabeling of a single set of marks, not two independent measures.

Classes:

• LdscDiagnosticPlotConfig –

  How the diagnostic is binned and framed.

• LdscDiagnosticPlotTask –

  Draw the LD-score-regression diagnostic plot for a SNPHeritabilityByLDSCTask analysis.

• LdscFit –

  The observed-scale LD-score-regression fit needed to draw the line and annotate it.

• LdscSumstats –

  The slice of a gwaslab Sumstats the diagnostic uses: its mutable variant table (which the

• MergedDiagnosticData –

  Per-variant chi-square and LD score for the variants shared by the summary statistics and

Functions:

• build_diagnostic_figure –

  Assemble the interactive figure: binned mean chi-square, the fitted line, a chi^2 = 1

• compute_chi2 –

  The per-variant chi-square. Prefer an existing Z column -- the same precedence

• gwaslab_observed_fit –

  Run gwaslab's LD-score regression on the observed scale and return its intercept and

• liability_h2 –

  The liability-scale heritability for a binary phenotype, or None for a quantitative one

• merge_chi2_with_ld_scores –

  Inner-join the summary statistics to the LD reference on rsID, pairing each variant's

• resolve_sample_size –

  The scalar N the fit and the right-axis rescaling use, mirroring how the heritability task

LdscDiagnosticPlotConfig

How the diagnostic is binned and framed.

x_range and chi2_m_over_n_range are left unset for per-trait autoscaling; pin them to a shared window to make a set of plots comparable across traits (equal LD-score and heritability-unit axes). chi2_m_over_n_range is given in the right-axis units (chi^2 * M / N).

Attributes:

• chi2_m_over_n_range (tuple[float, float] | None) –
• n_bins (int) –
• show_error_bars (bool) –
• title (str | None) –
• x_range (tuple[float, float] | None) –

chi2_m_over_n_range class-attribute instance-attribute

chi2_m_over_n_range: tuple[float, float] | None = None

n_bins class-attribute instance-attribute

n_bins: int = 20

show_error_bars class-attribute instance-attribute

show_error_bars: bool = False

title class-attribute instance-attribute

title: str | None = None

x_range class-attribute instance-attribute

x_range: tuple[float, float] | None = None

LdscDiagnosticPlotTask

Bases: Task

Draw the LD-score-regression diagnostic plot for a SNPHeritabilityByLDSCTask analysis.

Methods:

• create –
• execute –

Attributes:

• config (LdscDiagnosticPlotConfig) –
• deps (list[Task]) –
• fit_estimator (Callable[[LdscSumstats, str, GenomeBuild], LdscFit]) –
• ldsc_task (SNPHeritabilityByLDSCTask) –
• meta (Meta) –
• sumstats_reader (Callable[[Asset], LdscSumstats]) –

config instance-attribute

config: LdscDiagnosticPlotConfig

deps property

deps: list[Task]

fit_estimator class-attribute instance-attribute

fit_estimator: Callable[
    [LdscSumstats, str, GenomeBuild], LdscFit
] = field(default=gwaslab_observed_fit)

ldsc_task instance-attribute

ldsc_task: SNPHeritabilityByLDSCTask

meta instance-attribute

meta: Meta

sumstats_reader class-attribute instance-attribute

sumstats_reader: Callable[[Asset], LdscSumstats] = field(
    default=_read_sumstats
)

create classmethod

create(
    asset_id: str,
    ldsc_task: SNPHeritabilityByLDSCTask,
    config: LdscDiagnosticPlotConfig = LdscDiagnosticPlotConfig(),
    sumstats_reader: Callable[
        [Asset], LdscSumstats
    ] = _read_sumstats,
    fit_estimator: Callable[
        [LdscSumstats, str, GenomeBuild], LdscFit
    ] = gwaslab_observed_fit,
) -> LdscDiagnosticPlotTask

Source code in mecfs_bio/build_system/task/gwaslab/ldsc_diagnostic_plot_task.py

@classmethod
def create(
    cls,
    asset_id: str,
    ldsc_task: SNPHeritabilityByLDSCTask,
    config: LdscDiagnosticPlotConfig = LdscDiagnosticPlotConfig(),
    sumstats_reader: Callable[[Asset], LdscSumstats] = _read_sumstats,
    fit_estimator: Callable[
        [LdscSumstats, str, GenomeBuild], LdscFit
    ] = gwaslab_observed_fit,
) -> LdscDiagnosticPlotTask:
    source_meta = ldsc_task.meta
    assert isinstance(source_meta, ResultTableMeta)
    meta = GWASPlotDirectoryMeta(
        trait=source_meta.trait,
        project=source_meta.project,
        id=AssetId(asset_id),
    )
    return cls(
        meta=meta,
        ldsc_task=ldsc_task,
        config=config,
        sumstats_reader=sumstats_reader,
        fit_estimator=fit_estimator,
    )

execute

execute(scratch_dir: Path, fetch: Fetch, wf: WF) -> Asset

Source code in mecfs_bio/build_system/task/gwaslab/ldsc_diagnostic_plot_task.py

def execute(self, scratch_dir: Path, fetch: Fetch, wf: WF) -> Asset:
    sumstats = self.sumstats_reader(fetch(self._source_sumstats_id))
    sumstats.data = self.ldsc_task.pipe.process_pandas(sumstats.data)
    sumstats.data = drop_variants_with_degenerate_z(sumstats.data)
    _ensure_sample_size_column(
        sumstats.data, self.ldsc_task.set_N, self.ldsc_task.phenotype_info
    )
    n = resolve_sample_size(
        sumstats.data, self.ldsc_task.set_N, self.ldsc_task.phenotype_info
    )

    ld_asset = fetch(self._ld_ref_id)
    assert isinstance(ld_asset, DirectoryAsset)
    ld_frame = read_ld_scores(ld_asset.path).collect()
    total_m = total_m_5_50(ld_frame)

    # Bin from the processed sumstats before the fit: gwaslab's estimate mutates the Sumstats,
    # and the bins must reflect the variants as they enter the regression.
    merged = merge_chi2_with_ld_scores(sumstats.data, ld_frame.to_pandas())
    bins = bin_by_ld_score(merged.chi2, merged.ld, self.config.n_bins)

    ref_ld_chr = str(ld_asset.path) + self.ldsc_task.ld_file_filename_pattern
    fit = self.fit_estimator(sumstats, ref_ld_chr, self.ldsc_task.build)

    fig = build_diagnostic_figure(
        bins=bins,
        fit=fit,
        n=n,
        m=total_m,
        config=self.config,
        h2_liability=liability_h2(fit.h2_obs, self.ldsc_task.phenotype_info),
    )
    out_dir = scratch_dir / "ldsc_diagnostic_plot"
    write_plots_to_dir(out_dir, {"ldsc_diagnostic": fig})
    return DirectoryAsset(out_dir)

LdscFit

The observed-scale LD-score-regression fit needed to draw the line and annotate it.

Attributes:

• h2_obs (float) –
• intercept (float) –

h2_obs instance-attribute

h2_obs: float

intercept instance-attribute

intercept: float

LdscSumstats

Bases: Protocol

The slice of a gwaslab Sumstats the diagnostic uses: its mutable variant table (which the task reads and rewrites through the preprocessing pipe) plus enough to run observed-scale LD-score regression and read back its summary. gwaslab's Sumstats satisfies this; a test can inject a stand-in. ldsc_h2 is declared read-only (a property) so gwaslab's plain DataFrame attribute satisfies the pd.DataFrame | None slot covariantly.

Methods:

• estimate_h2_by_ldsc –
• infer_build –

Attributes:

• data (DataFrame) –
• ldsc_h2 (DataFrame | None) –
• meta (dict) –

data instance-attribute

data: DataFrame

ldsc_h2 property

ldsc_h2: DataFrame | None

meta instance-attribute

meta: dict

estimate_h2_by_ldsc

estimate_h2_by_ldsc(
    *,
    ref_ld_chr: str,
    w_ld_chr: str,
    samp_prev: float | None,
    pop_prev: float | None,
) -> None

Source code in mecfs_bio/build_system/task/gwaslab/ldsc_diagnostic_plot_task.py

def estimate_h2_by_ldsc(
    self,
    *,
    ref_ld_chr: str,
    w_ld_chr: str,
    samp_prev: float | None,
    pop_prev: float | None,
) -> None: ...

infer_build

infer_build() -> None

Source code in mecfs_bio/build_system/task/gwaslab/ldsc_diagnostic_plot_task.py

def infer_build(self) -> None: ...

MergedDiagnosticData

Per-variant chi-square and LD score for the variants shared by the summary statistics and the LD reference. Both arrays are 1-D of the same length, aligned row-for-row.

Attributes:

• chi2 (ndarray) –
• ld (ndarray) –

chi2 instance-attribute

chi2: ndarray

ld instance-attribute

ld: ndarray

build_diagnostic_figure

build_diagnostic_figure(
    bins: LdscDiagnosticBins,
    fit: LdscFit,
    n: float,
    m: float,
    config: LdscDiagnosticPlotConfig,
    h2_liability: float | None = None,
) -> go.Figure

Assemble the interactive figure: binned mean chi-square, the fitted line, a chi^2 = 1 reference, and a right-hand axis in chi^2 * M / N units whose slope reads as heritability.

h2_liability, when given (a binary phenotype), is annotated as a second heritability line.

Source code in mecfs_bio/build_system/task/gwaslab/ldsc_diagnostic_plot_task.py

def build_diagnostic_figure(
    bins: LdscDiagnosticBins,
    fit: LdscFit,
    n: float,
    m: float,
    config: LdscDiagnosticPlotConfig,
    h2_liability: float | None = None,
) -> go.Figure:
    """Assemble the interactive figure: binned mean chi-square, the fitted line, a chi^2 = 1
    reference, and a right-hand axis in chi^2 * M / N units whose slope reads as heritability.

    h2_liability, when given (a binary phenotype), is annotated as a second heritability line."""
    ld_range = bins.mean_ld.max() - bins.mean_ld.min()
    x_min = (
        config.x_range[0]
        if config.x_range
        else float(bins.mean_ld.min() - 0.02 * ld_range)
    )
    x_max = (
        config.x_range[1]
        if config.x_range
        else float(bins.mean_ld.max() + 0.02 * ld_range)
    )
    ld_line = np.array([x_min, x_max])
    y_fit = fit_line_chi2(fit.intercept, fit.h2_obs, n, m, ld_line)

    if config.chi2_m_over_n_range is not None:
        secondary_range = list(config.chi2_m_over_n_range)
        primary_range = [v * n / m for v in secondary_range]
    else:
        candidates = np.concatenate([bins.mean_chi2, y_fit, [1.0]])
        low, high = float(candidates.min()), float(candidates.max())
        pad = 0.05 * (high - low) if high > low else 0.05 * high
        primary_range = [low - pad, high + pad]
        secondary_range = [v * m / n for v in primary_range]

    fig = go.Figure()
    fig.add_trace(
        go.Scatter(
            x=bins.mean_ld,
            y=bins.mean_chi2,
            mode="markers",
            name="Binned mean χ²",
            marker=dict(size=9, color=_POINT_COLOR),
            error_y=dict(
                type="data", array=bins.se_chi2, visible=config.show_error_bars
            ),
            customdata=np.stack(
                [chi2_to_heritability_units(bins.mean_chi2, n, m), bins.count], axis=-1
            ),
            hovertemplate=(
                "LD score = %{x:.2f}<br>"
                "χ² = %{y:.4f}<br>"
                "χ²·M/N = %{customdata[0]:.3e}<br>"
                "variants in bin = %{customdata[1]}<extra></extra>"
            ),
        )
    )
    fig.add_trace(
        go.Scatter(
            x=ld_line,
            y=y_fit,
            mode="lines",
            name="LDSC fit",
            line=dict(color=_FIT_COLOR, width=2),
            hoverinfo="skip",
        )
    )
    fig.add_trace(
        go.Scatter(
            x=[x_min, x_max],
            y=[1.0, 1.0],
            mode="lines",
            name="χ² = 1",
            line=dict(color=_REFERENCE_COLOR, width=2, dash="dash"),
            hoverinfo="skip",
        )
    )
    # Plotly only draws a secondary axis that a trace references; this invisible trace spans the
    # right-axis range so the chi^2 * M / N axis renders without adding a visible mark.
    fig.add_trace(
        go.Scatter(
            x=[x_min, x_max],
            y=secondary_range,
            yaxis="y2",
            mode="markers",
            marker=dict(opacity=0),
            showlegend=False,
            hoverinfo="skip",
        )
    )
    annotation_lines = [f"h² (observed) = {fit.h2_obs:.4f}"]
    if h2_liability is not None:
        annotation_lines.append(f"h² (liability) = {h2_liability:.4f}")
    annotation_lines.append(f"intercept = {fit.intercept:.4f}")

    fig.update_layout(
        template="plotly_white",
        title=config.title,
        xaxis=dict(title="LD score", range=[x_min, x_max]),
        yaxis=dict(title="mean χ²", range=primary_range),
        yaxis2=dict(
            title="χ² · M / N",
            overlaying="y",
            side="right",
            range=secondary_range,
            showgrid=False,
        ),
        legend=dict(yanchor="top", y=0.99, xanchor="left", x=0.01),
        annotations=[
            # Nudged off the exact bottom-right corner and drawn as an opaque bordered box so it
            # cannot run into the lowest right-axis (chi^2 * M / N) tick label when the plot is small.
            dict(
                xref="paper",
                yref="paper",
                x=0.99,
                y=0.02,
                xanchor="right",
                yanchor="bottom",
                showarrow=False,
                align="right",
                bgcolor="white",
                bordercolor=_REFERENCE_COLOR,
                borderwidth=1,
                borderpad=5,
                text="<br>".join(annotation_lines),
            )
        ],
    )
    return fig

compute_chi2

compute_chi2(sumstats: DataFrame) -> np.ndarray

The per-variant chi-square. Prefer an existing Z column -- the same precedence drop_variants_with_degenerate_z uses -- and otherwise form Z from BETA / SE.

Source code in mecfs_bio/build_system/task/gwaslab/ldsc_diagnostic_plot_task.py

def compute_chi2(sumstats: pd.DataFrame) -> np.ndarray:
    """The per-variant chi-square. Prefer an existing Z column -- the same precedence
    drop_variants_with_degenerate_z uses -- and otherwise form Z from BETA / SE."""
    if LDSC_Z_COL in sumstats.columns:
        z = sumstats[LDSC_Z_COL].to_numpy(dtype=float)
    elif GWASLAB_BETA_COL in sumstats.columns and GWASLAB_SE_COL in sumstats.columns:
        z = sumstats[GWASLAB_BETA_COL].to_numpy(dtype=float) / sumstats[
            GWASLAB_SE_COL
        ].to_numpy(dtype=float)
    else:
        raise AssertionError(
            f"cannot compute chi-square: need a {LDSC_Z_COL} column, or both "
            f"{GWASLAB_BETA_COL} and {GWASLAB_SE_COL}"
        )
    return z * z

gwaslab_observed_fit

gwaslab_observed_fit(
    sumstats: LdscSumstats,
    ref_ld_chr: str,
    build: GenomeBuild,
) -> LdscFit

Run gwaslab's LD-score regression on the observed scale and return its intercept and heritability -- the same estimator the heritability task uses, so the drawn line matches what that task reports. None (not NaN) prevalences keep gwaslab on the observed scale: its summary labels the result 'Liability' whenever both prevalences are not None, which renames the h2 column, whereas None reports plain observed-scale h2 -- the scale in which the chi-square-units slope is reconstructable.

sumstats is a gwaslab Sumstats already carrying the analysis's preprocessing and sample size.

Source code in mecfs_bio/build_system/task/gwaslab/ldsc_diagnostic_plot_task.py

def gwaslab_observed_fit(
    sumstats: LdscSumstats, ref_ld_chr: str, build: GenomeBuild
) -> LdscFit:
    """Run gwaslab's LD-score regression on the observed scale and return its intercept and
    heritability -- the same estimator the heritability task uses, so the drawn line matches what
    that task reports. None (not NaN) prevalences keep gwaslab on the observed scale: its summary
    labels the result 'Liability' whenever both prevalences are not None, which renames the h2
    column, whereas None reports plain observed-scale h2 -- the scale in which the chi-square-units
    slope is reconstructable.

    sumstats is a gwaslab Sumstats already carrying the analysis's preprocessing and sample size."""
    sumstats.infer_build()
    assert sumstats.meta["gwaslab"]["genome_build"] == build, (
        f"sumstats build {sumstats.meta['gwaslab']['genome_build']} does not match {build}"
    )
    sumstats.estimate_h2_by_ldsc(
        ref_ld_chr=ref_ld_chr,
        w_ld_chr=ref_ld_chr,
        samp_prev=None,
        pop_prev=None,
    )
    ldsc_h2 = unwrap(sumstats.ldsc_h2)
    assert isinstance(ldsc_h2, pd.DataFrame)
    return LdscFit(
        intercept=float(ldsc_h2[_GWASLAB_LDSC_INTERCEPT_COL].iloc[0]),
        h2_obs=float(ldsc_h2[_GWASLAB_LDSC_H2_OBS_COL].iloc[0]),
    )

liability_h2

liability_h2(
    h2_obs: float, phenotype_info: PhenotypeInfo
) -> float | None

The liability-scale heritability for a binary phenotype, or None for a quantitative one (which has no liability scale). Uses gwaslab's own observed-to-liability conversion so the value matches the liability heritability the heritability task reports, with the sample prevalence as P and the estimated population prevalence as K.

Source code in mecfs_bio/build_system/task/gwaslab/ldsc_diagnostic_plot_task.py

def liability_h2(h2_obs: float, phenotype_info: PhenotypeInfo) -> float | None:
    """The liability-scale heritability for a binary phenotype, or None for a quantitative one
    (which has no liability scale). Uses gwaslab's own observed-to-liability conversion so the value
    matches the liability heritability the heritability task reports, with the sample prevalence as
    P and the estimated population prevalence as K."""
    if isinstance(phenotype_info, BinaryPhenotypeSampleInfo):
        return float(
            h2_obs_to_liab(
                h2_obs,
                phenotype_info.sample_prevalence,
                phenotype_info.estimated_population_prevalence,
            )
        )
    return None

merge_chi2_with_ld_scores

merge_chi2_with_ld_scores(
    sumstats: DataFrame, ld_scores: DataFrame
) -> MergedDiagnosticData

Inner-join the summary statistics to the LD reference on rsID, pairing each variant's chi-square with its LD score. Variants missing from either side drop out, which restricts the plot to the LD reference's variant set (HapMap3, in the standard reference).

Source code in mecfs_bio/build_system/task/gwaslab/ldsc_diagnostic_plot_task.py

def merge_chi2_with_ld_scores(
    sumstats: pd.DataFrame, ld_scores: pd.DataFrame
) -> MergedDiagnosticData:
    """Inner-join the summary statistics to the LD reference on rsID, pairing each variant's
    chi-square with its LD score. Variants missing from either side drop out, which restricts the
    plot to the LD reference's variant set (HapMap3, in the standard reference)."""
    joined = pd.DataFrame(
        {
            GWASLAB_RSID_COL: sumstats[GWASLAB_RSID_COL].to_numpy(),
            "chi2": compute_chi2(sumstats),
        }
    ).merge(
        ld_scores[[LD_SCORE_RSID_COL, LD_SCORE_LD_SCORE_COL]],
        left_on=GWASLAB_RSID_COL,
        right_on=LD_SCORE_RSID_COL,
        how="inner",
    )
    return MergedDiagnosticData(
        chi2=joined["chi2"].to_numpy(dtype=float),
        ld=joined[LD_SCORE_LD_SCORE_COL].to_numpy(dtype=float),
    )

resolve_sample_size

resolve_sample_size(
    sumstats: DataFrame,
    set_n: int | None,
    phenotype_info: PhenotypeInfo,
) -> float

The scalar N the fit and the right-axis rescaling use, mirroring how the heritability task sources it: an explicit override, then a per-variant N column (averaged), then the phenotype's total sample size.

Source code in mecfs_bio/build_system/task/gwaslab/ldsc_diagnostic_plot_task.py

def resolve_sample_size(
    sumstats: pd.DataFrame, set_n: int | None, phenotype_info: PhenotypeInfo
) -> float:
    """The scalar N the fit and the right-axis rescaling use, mirroring how the heritability task
    sources it: an explicit override, then a per-variant N column (averaged), then the phenotype's
    total sample size."""
    if set_n is not None:
        return float(set_n)
    if GWASLAB_SAMPLE_SIZE_COLUMN in sumstats.columns:
        return float(sumstats[GWASLAB_SAMPLE_SIZE_COLUMN].mean())
    assert phenotype_info.total_sample_size is not None, (
        "no sample size available: set a sample size on the heritability task, provide an "
        f"{GWASLAB_SAMPLE_SIZE_COLUMN} column, or a phenotype total_sample_size"
    )
    return float(phenotype_info.total_sample_size)