Skip to content

Coloc

gentropy.method.colocalisation.Coloc

Calculate bayesian colocalisation based on overlapping signals from credible sets.

Based on the R COLOC package, which uses the Bayes factors from the credible set to estimate the posterior probability of colocalisation. This method makes the simplifying assumption that only one single causal variant exists for any given trait in any genomic region.

Hypothesis Description
H0 no association with either trait in the region
H1 association with trait 1 only
H2 association with trait 2 only
H3 both traits are associated, but have different single causal variants
H4 both traits are associated and share the same single causal variant

Bayes factors required

Coloc requires the availability of Bayes factors (BF) for each variant in the credible set (logBF column).

Attributes:

Name Type Description
PSEUDOCOUNT float

Pseudocount to avoid log(0). Defaults to 1e-10.
Source code in src/gentropy/method/colocalisation.py 
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
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
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
class Coloc:
    """Calculate bayesian colocalisation based on overlapping signals from credible sets.

    Based on the [R COLOC package](https://github.com/chr1swallace/coloc/blob/main/R/claudia.R), which uses the Bayes factors from the credible set to estimate the posterior probability of colocalisation. This method makes the simplifying assumption that **only one single causal variant** exists for any given trait in any genomic region.

    | Hypothesis    | Description                                                           |
    | ------------- | --------------------------------------------------------------------- |
    | H<sub>0</sub> | no association with either trait in the region                        |
    | H<sub>1</sub> | association with trait 1 only                                         |
    | H<sub>2</sub> | association with trait 2 only                                         |
    | H<sub>3</sub> | both traits are associated, but have different single causal variants |
    | H<sub>4</sub> | both traits are associated and share the same single causal variant   |

    !!! warning "Bayes factors required"

        Coloc requires the availability of Bayes factors (BF) for each variant in the credible set (`logBF` column).

    Attributes:
        PSEUDOCOUNT (float): Pseudocount to avoid log(0). Defaults to 1e-10.
    """

    METHOD_NAME: str = "COLOC"
    METHOD_METRIC: str = "llr"
    PSEUDOCOUNT: float = 1e-10

    @staticmethod
    def _get_posteriors(all_bfs: NDArray[np.float64]) -> DenseVector:
        """Calculate posterior probabilities for each hypothesis.

        Args:
            all_bfs (NDArray[np.float64]): h0-h4 bayes factors

        Returns:
            DenseVector: Posterior

        Example:
            >>> l = np.array([0.2, 0.1, 0.05, 0])
            >>> Coloc._get_posteriors(l)
            DenseVector([0.279, 0.2524, 0.2401, 0.2284])
        """
        diff = all_bfs - get_logsum(all_bfs)
        bfs_posteriors = np.exp(diff)
        return Vectors.dense(bfs_posteriors)

    @classmethod
    def colocalise(
        cls: type[Coloc],
        overlapping_signals: StudyLocusOverlap,
        priorc1: float = 1e-4,
        priorc2: float = 1e-4,
        priorc12: float = 1e-5,
    ) -> Colocalisation:
        """Calculate bayesian colocalisation based on overlapping signals.

        Args:
            overlapping_signals (StudyLocusOverlap): overlapping peaks

            priorc1 (float): Prior on variant being causal for trait 1. Defaults to 1e-4.
            priorc2 (float): Prior on variant being causal for trait 2. Defaults to 1e-4.
            priorc12 (float): Prior on variant being causal for traits 1 and 2. Defaults to 1e-5.

        Returns:
            Colocalisation: Colocalisation results
        """
        # register udfs
        logsum = f.udf(get_logsum, DoubleType())
        posteriors = f.udf(Coloc._get_posteriors, VectorUDT())
        return Colocalisation(
            _df=(
                overlapping_signals.df.select("*", "statistics.*")
                # Before summing log_BF columns nulls need to be filled with 0:
                .fillna(0, subset=["left_logBF", "right_logBF"])
                # Sum of log_BFs for each pair of signals
                .withColumn(
                    "sum_log_bf",
                    f.col("left_logBF") + f.col("right_logBF"),
                )
                # Group by overlapping peak and generating dense vectors of log_BF:
                .groupBy("chromosome", "leftStudyLocusId", "rightStudyLocusId")
                .agg(
                    f.count("*").alias("numberColocalisingVariants"),
                    fml.array_to_vector(f.collect_list(f.col("left_logBF"))).alias(
                        "left_logBF"
                    ),
                    fml.array_to_vector(f.collect_list(f.col("right_logBF"))).alias(
                        "right_logBF"
                    ),
                    fml.array_to_vector(f.collect_list(f.col("sum_log_bf"))).alias(
                        "sum_log_bf"
                    ),
                )
                .withColumn("logsum1", logsum(f.col("left_logBF")))
                .withColumn("logsum2", logsum(f.col("right_logBF")))
                .withColumn("logsum12", logsum(f.col("sum_log_bf")))
                .drop("left_logBF", "right_logBF", "sum_log_bf")
                # Add priors
                # priorc1 Prior on variant being causal for trait 1
                .withColumn("priorc1", f.lit(priorc1))
                # priorc2 Prior on variant being causal for trait 2
                .withColumn("priorc2", f.lit(priorc2))
                # priorc12 Prior on variant being causal for traits 1 and 2
                .withColumn("priorc12", f.lit(priorc12))
                # h0-h2
                .withColumn("lH0bf", f.lit(0))
                .withColumn("lH1bf", f.log(f.col("priorc1")) + f.col("logsum1"))
                .withColumn("lH2bf", f.log(f.col("priorc2")) + f.col("logsum2"))
                # h3
                .withColumn("sumlogsum", f.col("logsum1") + f.col("logsum2"))
                .withColumn("max", f.greatest("sumlogsum", "logsum12"))
                .withColumn(
                    "logdiff",
                    f.when(
                        f.col("sumlogsum") == f.col("logsum12"), Coloc.PSEUDOCOUNT
                    ).otherwise(
                        f.col("max")
                        + f.log(
                            f.exp(f.col("sumlogsum") - f.col("max"))
                            - f.exp(f.col("logsum12") - f.col("max"))
                        )
                    ),
                )
                .withColumn(
                    "lH3bf",
                    f.log(f.col("priorc1"))
                    + f.log(f.col("priorc2"))
                    + f.col("logdiff"),
                )
                .drop("right_logsum", "left_logsum", "sumlogsum", "max", "logdiff")
                # h4
                .withColumn("lH4bf", f.log(f.col("priorc12")) + f.col("logsum12"))
                # cleaning
                .drop(
                    "priorc1", "priorc2", "priorc12", "logsum1", "logsum2", "logsum12"
                )
                # posteriors
                .withColumn(
                    "allBF",
                    fml.array_to_vector(
                        f.array(
                            f.col("lH0bf"),
                            f.col("lH1bf"),
                            f.col("lH2bf"),
                            f.col("lH3bf"),
                            f.col("lH4bf"),
                        )
                    ),
                )
                .withColumn(
                    "posteriors", fml.vector_to_array(posteriors(f.col("allBF")))
                )
                .withColumn("h0", f.col("posteriors").getItem(0))
                .withColumn("h1", f.col("posteriors").getItem(1))
                .withColumn("h2", f.col("posteriors").getItem(2))
                .withColumn("h3", f.col("posteriors").getItem(3))
                .withColumn("h4", f.col("posteriors").getItem(4))
                .withColumn("h4h3", f.col("h4") / f.col("h3"))
                .withColumn("log2h4h3", f.log2(f.col("h4h3")))
                # clean up
                .drop(
                    "posteriors",
                    "allBF",
                    "h4h3",
                    "lH0bf",
                    "lH1bf",
                    "lH2bf",
                    "lH3bf",
                    "lH4bf",
                )
                .withColumn("colocalisationMethod", f.lit(cls.METHOD_NAME))
            ),
            _schema=Colocalisation.get_schema(),
        )

colocalise(overlapping_signals: StudyLocusOverlap, priorc1: float = 0.0001, priorc2: float = 0.0001, priorc12: float = 1e-05) -> Colocalisation classmethod

Calculate bayesian colocalisation based on overlapping signals.

Parameters:

Name Type Description Default
overlapping_signals StudyLocusOverlap

overlapping peaks

 required
priorc1 float

Prior on variant being causal for trait 1. Defaults to 1e-4.

 0.0001
priorc2 float

Prior on variant being causal for trait 2. Defaults to 1e-4.

 0.0001
priorc12 float

Prior on variant being causal for traits 1 and 2. Defaults to 1e-5.

 1e-05

Returns:

Name Type Description
Colocalisation Colocalisation

Colocalisation results
Source code in src/gentropy/method/colocalisation.py 
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
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
@classmethod
def colocalise(
    cls: type[Coloc],
    overlapping_signals: StudyLocusOverlap,
    priorc1: float = 1e-4,
    priorc2: float = 1e-4,
    priorc12: float = 1e-5,
) -> Colocalisation:
    """Calculate bayesian colocalisation based on overlapping signals.

    Args:
        overlapping_signals (StudyLocusOverlap): overlapping peaks

        priorc1 (float): Prior on variant being causal for trait 1. Defaults to 1e-4.
        priorc2 (float): Prior on variant being causal for trait 2. Defaults to 1e-4.
        priorc12 (float): Prior on variant being causal for traits 1 and 2. Defaults to 1e-5.

    Returns:
        Colocalisation: Colocalisation results
    """
    # register udfs
    logsum = f.udf(get_logsum, DoubleType())
    posteriors = f.udf(Coloc._get_posteriors, VectorUDT())
    return Colocalisation(
        _df=(
            overlapping_signals.df.select("*", "statistics.*")
            # Before summing log_BF columns nulls need to be filled with 0:
            .fillna(0, subset=["left_logBF", "right_logBF"])
            # Sum of log_BFs for each pair of signals
            .withColumn(
                "sum_log_bf",
                f.col("left_logBF") + f.col("right_logBF"),
            )
            # Group by overlapping peak and generating dense vectors of log_BF:
            .groupBy("chromosome", "leftStudyLocusId", "rightStudyLocusId")
            .agg(
                f.count("*").alias("numberColocalisingVariants"),
                fml.array_to_vector(f.collect_list(f.col("left_logBF"))).alias(
                    "left_logBF"
                ),
                fml.array_to_vector(f.collect_list(f.col("right_logBF"))).alias(
                    "right_logBF"
                ),
                fml.array_to_vector(f.collect_list(f.col("sum_log_bf"))).alias(
                    "sum_log_bf"
                ),
            )
            .withColumn("logsum1", logsum(f.col("left_logBF")))
            .withColumn("logsum2", logsum(f.col("right_logBF")))
            .withColumn("logsum12", logsum(f.col("sum_log_bf")))
            .drop("left_logBF", "right_logBF", "sum_log_bf")
            # Add priors
            # priorc1 Prior on variant being causal for trait 1
            .withColumn("priorc1", f.lit(priorc1))
            # priorc2 Prior on variant being causal for trait 2
            .withColumn("priorc2", f.lit(priorc2))
            # priorc12 Prior on variant being causal for traits 1 and 2
            .withColumn("priorc12", f.lit(priorc12))
            # h0-h2
            .withColumn("lH0bf", f.lit(0))
            .withColumn("lH1bf", f.log(f.col("priorc1")) + f.col("logsum1"))
            .withColumn("lH2bf", f.log(f.col("priorc2")) + f.col("logsum2"))
            # h3
            .withColumn("sumlogsum", f.col("logsum1") + f.col("logsum2"))
            .withColumn("max", f.greatest("sumlogsum", "logsum12"))
            .withColumn(
                "logdiff",
                f.when(
                    f.col("sumlogsum") == f.col("logsum12"), Coloc.PSEUDOCOUNT
                ).otherwise(
                    f.col("max")
                    + f.log(
                        f.exp(f.col("sumlogsum") - f.col("max"))
                        - f.exp(f.col("logsum12") - f.col("max"))
                    )
                ),
            )
            .withColumn(
                "lH3bf",
                f.log(f.col("priorc1"))
                + f.log(f.col("priorc2"))
                + f.col("logdiff"),
            )
            .drop("right_logsum", "left_logsum", "sumlogsum", "max", "logdiff")
            # h4
            .withColumn("lH4bf", f.log(f.col("priorc12")) + f.col("logsum12"))
            # cleaning
            .drop(
                "priorc1", "priorc2", "priorc12", "logsum1", "logsum2", "logsum12"
            )
            # posteriors
            .withColumn(
                "allBF",
                fml.array_to_vector(
                    f.array(
                        f.col("lH0bf"),
                        f.col("lH1bf"),
                        f.col("lH2bf"),
                        f.col("lH3bf"),
                        f.col("lH4bf"),
                    )
                ),
            )
            .withColumn(
                "posteriors", fml.vector_to_array(posteriors(f.col("allBF")))
            )
            .withColumn("h0", f.col("posteriors").getItem(0))
            .withColumn("h1", f.col("posteriors").getItem(1))
            .withColumn("h2", f.col("posteriors").getItem(2))
            .withColumn("h3", f.col("posteriors").getItem(3))
            .withColumn("h4", f.col("posteriors").getItem(4))
            .withColumn("h4h3", f.col("h4") / f.col("h3"))
            .withColumn("log2h4h3", f.log2(f.col("h4h3")))
            # clean up
            .drop(
                "posteriors",
                "allBF",
                "h4h3",
                "lH0bf",
                "lH1bf",
                "lH2bf",
                "lH3bf",
                "lH4bf",
            )
            .withColumn("colocalisationMethod", f.lit(cls.METHOD_NAME))
        ),
        _schema=Colocalisation.get_schema(),
    )
  • 2023-01-20
  • 2024-01-18
  • Contributors