graphld.blup

Best Linear Unbiased Prediction (BLUP) for effect size estimation.

Under the infinitesimal model with per-s.d. effect sizes $\beta \sim N(0, D)$, the BLUP effect sizes are:

$$E(\beta) = \sqrt{n} D (nD + R^{-1})^{-1} R^{-1}z$$

where $R^{-1}$ is approximated with the LDGM precision matrix. Public BLUP workflows take heritability for the analyzed variant scope and set $\mathrm{trace}(D)$ equal to that value across matched LDGM effect indices.

For usage examples, see the BLUP guide.

blup

BLUP

Bases: ParallelProcessor

Computes the best linear unbiased predictor using LDGMs and GWAS summary statistics.

prepare_block_data classmethod

prepare_block_data(metadata: DataFrame, **kwargs: Any) -> list[tuple]

Split summary statistics into blocks whose positions match the LDGMs.

Parameters:

Name	Type	Description	Default
metadata	DataFrame	DataFrame containing LDGM metadata	required
**kwargs	Any	Additional arguments from run(), including: sumstats: DataFrame containing summary statistics	{}

Returns:

Type	Description
list[tuple]	List of tuples containing a block-specific summary-statistics
list[tuple]	DataFrame and its offset in the concatenated output.

Source code in src/graphld/blup.py

@classmethod
def prepare_block_data(cls, metadata: pl.DataFrame, **kwargs: Any) -> list[tuple]:
    """Split summary statistics into blocks whose positions match the LDGMs.

    Args:
        metadata: DataFrame containing LDGM metadata
        **kwargs: Additional arguments from run(), including:
            sumstats: DataFrame containing summary statistics

    Returns:
        List of tuples containing a block-specific summary-statistics
        DataFrame and its offset in the concatenated output.
    """
    sumstats = kwargs.get('sumstats')

    # Partition summary statistics into blocks
    sumstats_blocks: list[pl.DataFrame] = partition_variants(metadata, sumstats)

    cumulative_num_variants = np.cumsum(np.array([len(df) for df in sumstats_blocks]))
    cumulative_num_variants = [0] + list(cumulative_num_variants[:-1])

    return list(zip(sumstats_blocks, cumulative_num_variants, strict=False))

create_shared_memory staticmethod

create_shared_memory(metadata: DataFrame, block_data: list[tuple], **kwargs: Any) -> SharedData

Create the BLUP output array for partitioned summary statistics.

Parameters:

Name	Type	Description	Default
metadata	DataFrame	Metadata DataFrame containing block information	required
block_data	list[tuple]	List of block-specific summary-statistics DataFrames and offsets	required
**kwargs	Any	Not used	{}

Returns:

Type	Description
SharedData	SharedData containing the BLUP weight array.

Source code in src/graphld/blup.py

@staticmethod
def create_shared_memory(metadata: pl.DataFrame, block_data: list[tuple], **kwargs: Any) -> SharedData:
    """Create the BLUP output array for partitioned summary statistics.

    Args:
        metadata: Metadata DataFrame containing block information
        block_data: List of block-specific summary-statistics DataFrames
            and offsets
        **kwargs: Not used

    Returns:
        SharedData containing the BLUP weight array.
    """
    total_variants = sum([len(df) for df, _ in block_data])
    return SharedData({
        'beta': total_variants,    # BLUP effect sizes
    })

process_block classmethod

process_block(ldgm: PrecisionOperator, flag: Value, shared_data: SharedData, block_offset: int, block_data: tuple, worker_params: tuple) -> None

Run BLUP on a single block.

Source code in src/graphld/blup.py

@classmethod
def process_block(cls, ldgm: PrecisionOperator, flag: Value,
                 shared_data: SharedData, block_offset: int,
                 block_data: tuple,
                 worker_params: tuple) -> None:
    """Run BLUP on a single block."""
    per_effect_variance, sample_size, match_by_position = worker_params
    assert isinstance(per_effect_variance, float), "per-effect variance must be a float"
    assert isinstance(block_data, tuple), "block_data must be a tuple"
    sumstats, variant_offset = block_data
    num_variants = len(sumstats)

    # Merge annotations with LDGM variant info and get indices of merged variants
    ldgm, sumstat_indices = cls._merge_sumstats(ldgm, sumstats, match_by_position)

    # Keep only first occurrence of each index
    first_index_mask = ldgm.variant_info.select(pl.col('index').is_first_distinct()).to_numpy().flatten()
    ldgm.variant_info = ldgm.variant_info.filter(first_index_mask)
    sumstat_indices = sumstat_indices[first_index_mask]

    # Get Z-scores from the merged variant info
    z = ldgm.variant_info.select('Z').to_numpy()

    # Compute the BLUP for this block
    beta = ldgm @ z
    ldgm.update_matrix(np.full(ldgm.shape[0], sample_size * per_effect_variance))
    beta = np.sqrt(sample_size) * per_effect_variance * ldgm.solve(beta)
    ldgm.del_factor()

    # Store results for variants that were successfully merged
    beta_reshaped = np.zeros(num_variants)
    # Get indices of variants that were actually merged
    beta_reshaped[sumstat_indices] = beta.flatten()

    # Update the shared memory array
    block_slice = slice(variant_offset, variant_offset + num_variants)
    shared_data['beta', block_slice] = beta_reshaped

supervise classmethod

supervise(manager: WorkerManager, shared_data: Dict[str, Any], block_data: list, **kwargs: Any) -> pl.DataFrame

Supervise worker processes and collect results.

Parameters:

Name	Type	Description	Default
manager	WorkerManager	Worker manager	required
shared_data	Dict[str, Any]	Dictionary of shared memory arrays	required
**kwargs	Any	Additional arguments	{}

Returns:

Type	Description
DataFrame	DataFrame containing summary statistics with BLUP weights

Source code in src/graphld/blup.py

@classmethod
def supervise(cls, manager: WorkerManager, shared_data: Dict[str, Any], block_data: list, **kwargs: Any) -> pl.DataFrame:
    """Supervise worker processes and collect results.

    Args:
        manager: Worker manager
        shared_data: Dictionary of shared memory arrays
        **kwargs: Additional arguments

    Returns:
        DataFrame containing summary statistics with BLUP weights
    """

    manager.start_workers()
    manager.await_workers()
    beta = shared_data['beta']
    sumstats = pl.concat([df for df, _ in block_data])
    return sumstats.with_columns(pl.Series('weight', beta))

compute_blup classmethod

compute_blup(ldgm_metadata_path: str, sumstats: DataFrame, heritability: Optional[float] = None, sample_size: Optional[float] = None, populations: Optional[Union[str, List[str]]] = None, chromosomes: Optional[Union[int, List[int]]] = None, num_processes: Optional[int] = None, run_in_serial: bool = False, match_by_position: bool = False, verbose: bool = False, sigmasq: Optional[float] = None) -> pl.DataFrame

Compute BLUP weights for multiple LD blocks.

Parameters:

Name	Type	Description	Default
ldgm_metadata_path	str	Path to metadata CSV file	required
sumstats	DataFrame	Summary-statistics DataFrame containing Z scores	required
heritability	Optional[float]	Heritability for the analyzed variant scope. Internally, BLUP uses D = (heritability / m) I, where m is the number of unique matched LDGM effect indices.	None
sample_size	Optional[float]	GWAS sample size	None
populations	Optional[Union[str, List[str]]]	Optional population name(s)	None
chromosomes	Optional[Union[int, List[int]]]	Optional chromosome or list of chromosomes	None
num_processes	Optional[int]	Optional number of processes	None
run_in_serial	bool	Whether to run in serial mode	False
match_by_position	bool	Whether to match variants by position instead of ID	False
verbose	bool	Print additional information if True	False
sigmasq	Optional[float]	Removed compatibility keyword. Passing it raises ValueError.	None

Returns:

Type	Description
DataFrame	DataFrame with BLUP weights added to the partitioned summary statistics.

Source code in src/graphld/blup.py

@classmethod
def compute_blup(cls,
        ldgm_metadata_path: str,
        sumstats: pl.DataFrame,
        heritability: Optional[float] = None,
        sample_size: Optional[float] = None,
        populations: Optional[Union[str, List[str]]] = None,
        chromosomes: Optional[Union[int, List[int]]] = None,
        num_processes: Optional[int] = None,
        run_in_serial: bool = False,
        match_by_position: bool = False,
        verbose: bool = False,
        sigmasq: Optional[float] = None,
        ) -> pl.DataFrame:
    """Compute BLUP weights for multiple LD blocks.

    Args:
        ldgm_metadata_path: Path to metadata CSV file
        sumstats: Summary-statistics DataFrame containing Z scores
        heritability: Heritability for the analyzed variant scope. Internally, BLUP uses
            D = (heritability / m) I, where m is the number of unique
            matched LDGM effect indices.
        sample_size: GWAS sample size
        populations: Optional population name(s)
        chromosomes: Optional chromosome or list of chromosomes
        num_processes: Optional number of processes
        run_in_serial: Whether to run in serial mode
        match_by_position: Whether to match variants by position instead of ID
        verbose: Print additional information if True
        sigmasq: Removed compatibility keyword. Passing it raises ValueError.

    Returns:
        DataFrame with BLUP weights added to the partitioned summary statistics.
    """
    if sigmasq is not None:
        raise ValueError(
            "The BLUP sigmasq keyword is no longer supported. "
            "Pass total heritability for the analyzed variant scope with heritability=."
        )
    if heritability is None:
        raise ValueError("BLUP requires heritability for the analyzed variant scope")
    if sample_size is None:
        raise ValueError("BLUP requires sample_size")
    if not 0 <= heritability <= 1:
        raise ValueError(f"Heritability must be between 0 and 1, got {heritability}")

    matched_effects = cls._count_matched_effects(
        ldgm_metadata_path=ldgm_metadata_path,
        sumstats=sumstats,
        populations=populations,
        chromosomes=chromosomes,
        match_by_position=match_by_position,
    )
    if matched_effects == 0:
        raise ValueError("No summary statistics matched LDGM variants for BLUP")

    per_effect_variance = heritability / matched_effects
    run_fn = cls.run_serial if run_in_serial else cls.run
    result = run_fn(
        ldgm_metadata_path=ldgm_metadata_path,
        populations=populations,
        chromosomes=chromosomes,
        num_processes=num_processes,
        worker_params=(per_effect_variance, sample_size, match_by_position),
        sumstats=sumstats
    )

    if verbose:
        print(f"Number of variants in summary statistics: {len(result)}")
        print(f"Number of matched LDGM effect indices: {matched_effects}")
        print(f"Per-effect variance: {per_effect_variance}")
        nonzero_count = (result['weight'] != 0).sum()
        print(f"Number of variants with nonzero weights: {nonzero_count}")

    return result

run_blup

run_blup(*args: Any, **kwargs: Any) -> pl.DataFrame

Compute Best Linear Unbiased Prediction (BLUP) weights.

Positional and keyword arguments are forwarded to :meth:BLUP.compute_blup.

Parameters:

Name	Type	Description	Default
*args	Any	Positional arguments for :meth:BLUP.compute_blup.	()
**kwargs	Any	Keyword arguments for :meth:BLUP.compute_blup.	{}

Returns:

Type	Description
DataFrame	DataFrame with BLUP weights and associated statistics.

Source code in src/graphld/blup.py

def run_blup(*args: Any, **kwargs: Any) -> pl.DataFrame:
    """Compute Best Linear Unbiased Prediction (BLUP) weights.

    Positional and keyword arguments are forwarded to :meth:`BLUP.compute_blup`.

    Args:
        *args: Positional arguments for :meth:`BLUP.compute_blup`.
        **kwargs: Keyword arguments for :meth:`BLUP.compute_blup`.

    Returns:
        DataFrame with BLUP weights and associated statistics.
    """
    return BLUP.compute_blup(*args, **kwargs)