aboutsummaryrefslogtreecommitdiff
path: root/gn3/computations/partial_correlations.py
blob: f43c4d480866d5ea2c41a720415f1ec464fc873c (about) (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
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
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
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
"""
This module deals with partial correlations.

It is an attempt to migrate over the partial correlations feature from
GeneNetwork1.
"""

import math
from functools import reduce
from typing import Any, Tuple, Union, Sequence
from scipy.stats import pearsonr, spearmanr

import pandas
import pingouin

from gn3.settings import TEXTDIR
from gn3.data_helpers import parse_csv_line

def control_samples(controls: Sequence[dict], sampleslist: Sequence[str]):
    """
    Fetches data for the control traits.

    This migrates `web/webqtl/correlation/correlationFunction.controlStrain` in
    GN1, with a few modifications to the arguments passed in.

    PARAMETERS:
    controls: A map of sample names to trait data. Equivalent to the `cvals`
        value in the corresponding source function in GN1.
    sampleslist: A list of samples. Equivalent to `strainlst` in the
        corresponding source function in GN1
    """
    def __process_control__(trait_data):
        def __process_sample__(acc, sample):
            if sample in trait_data["data"].keys():
                sample_item = trait_data["data"][sample]
                val = sample_item["value"]
                if val is not None:
                    return (
                        acc[0] + (sample,),
                        acc[1] + (val,),
                        acc[2] + (sample_item["variance"],))
            return acc
        return reduce(
            __process_sample__, sampleslist, (tuple(), tuple(), tuple()))

    return reduce(
        lambda acc, item: (
            acc[0] + (item[0],),
            acc[1] + (item[1],),
            acc[2] + (item[2],),
            acc[3] + (len(item[0]),),
        ),
        [__process_control__(trait_data) for trait_data in controls],
        (tuple(), tuple(), tuple(), tuple()))

def dictify_by_samples(samples_vals_vars: Sequence[Sequence]) -> Sequence[dict]:
    """
    Build a sequence of dictionaries from a sequence of separate sequences of
    samples, values and variances.

    This is a partial migration of
    `web.webqtl.correlation.correlationFunction.fixStrains` function in GN1.
    This implementation extracts code that will find common use, and that will
    find use in more than one place.
    """
    return tuple(
        {
            sample: {"sample_name": sample, "value": val, "variance": var}
            for sample, val, var in zip(*trait_line)
        } for trait_line in zip(*(samples_vals_vars[0:3])))

def fix_samples(primary_trait: dict, control_traits: Sequence[dict]) -> Sequence[Sequence[Any]]:
    """
    Corrects sample_names, values and variance such that they all contain only
    those samples that are common to the reference trait and all control traits.

    This is a partial migration of the
    `web.webqtl.correlation.correlationFunction.fixStrain` function in GN1.
    """
    primary_samples = tuple(
        present[0] for present in
        ((sample, all(sample in control.keys() for control in control_traits))
         for sample in primary_trait.keys())
        if present[1])
    control_vals_vars: tuple = reduce(
        lambda acc, x: (acc[0] + (x[0],), acc[1] + (x[1],)),
        ((item["value"], item["variance"])
         for sublist in [tuple(control.values()) for control in control_traits]
         for item in sublist),
        (tuple(), tuple()))
    return (
        primary_samples,
        tuple(primary_trait[sample]["value"] for sample in primary_samples),
        control_vals_vars[0],
        tuple(primary_trait[sample]["variance"] for sample in primary_samples),
        control_vals_vars[1])

def find_identical_traits(
        primary_name: str, primary_value: float, control_names: Tuple[str, ...],
        control_values: Tuple[float, ...]) -> Tuple[str, ...]:
    """
    Find traits that have the same value when the values are considered to
    3 decimal places.

    This is a migration of the
    `web.webqtl.correlation.correlationFunction.findIdenticalTraits` function in
    GN1.
    """
    def __merge_identicals__(
            acc: Tuple[str, ...],
            ident: Tuple[str, Tuple[str, ...]]) -> Tuple[str, ...]:
        return acc + ident[1]

    def __dictify_controls__(acc, control_item):
        ckey = "{:.3f}".format(control_item[0])
        return {**acc, ckey: acc.get(ckey, tuple()) + (control_item[1],)}

    return (reduce(## for identical control traits
        __merge_identicals__,
        (item for item in reduce(# type: ignore[var-annotated]
            __dictify_controls__, zip(control_values, control_names),
            {}).items() if len(item[1]) > 1),
        tuple())
            or
            reduce(## If no identical control traits, try primary and controls
                __merge_identicals__,
                (item for item in reduce(# type: ignore[var-annotated]
                    __dictify_controls__,
                    zip((primary_value,) + control_values,
                        (primary_name,) + control_names), {}).items()
                 if len(item[1]) > 1),
                tuple()))

def tissue_correlation(
        primary_trait_values: Tuple[float, ...],
        target_trait_values: Tuple[float, ...],
        method: str) -> Tuple[float, float]:
    """
    Compute the correlation between the primary trait values, and the values of
    a single target value.

    This migrates the `cal_tissue_corr` function embedded in the larger
    `web.webqtl.correlation.correlationFunction.batchCalTissueCorr` function in
    GeneNetwork1.
    """
    def spearman_corr(*args):
        result = spearmanr(*args)
        return (result.correlation, result.pvalue)

    method_fns = {"pearson": pearsonr, "spearman": spearman_corr}

    assert len(primary_trait_values) == len(target_trait_values), (
        "The lengths of the `primary_trait_values` and `target_trait_values` "
        "must be equal")
    assert method in method_fns.keys(), (
        "Method must be one of: {}".format(",".join(method_fns.keys())))

    corr, pvalue = method_fns[method](primary_trait_values, target_trait_values)
    return (corr, pvalue)

def batch_computed_tissue_correlation(
        primary_trait_values: Tuple[float, ...], target_traits_dict: dict,
        method: str) -> Tuple[dict, dict]:
    """
    This is a migration of the
    `web.webqtl.correlation.correlationFunction.batchCalTissueCorr` function in
    GeneNetwork1
    """
    def __corr__(acc, target):
        corr = tissue_correlation(primary_trait_values, target[1], method)
        return ({**acc[0], target[0]: corr[0]}, {**acc[0], target[1]: corr[1]})
    return reduce(__corr__, target_traits_dict.items(), ({}, {}))

def correlations_of_all_tissue_traits(
        primary_trait_symbol_value_dict: dict, symbol_value_dict: dict,
        method: str) -> Tuple[dict, dict]:
    """
    Computes and returns the correlation of all tissue traits.

    This is a migration of the
    `web.webqtl.correlation.correlationFunction.calculateCorrOfAllTissueTrait`
    function in GeneNetwork1.
    """
    primary_trait_values = tuple(primary_trait_symbol_value_dict.values())[0]
    return batch_computed_tissue_correlation(
        primary_trait_values, symbol_value_dict, method)

def good_dataset_samples_indexes(
        samples: Tuple[str, ...],
        samples_from_file: Tuple[str, ...]) -> Tuple[int, ...]:
    """
    Return the indexes of the items in `samples_from_files` that are also found
    in `samples`.

    This is a partial migration of the
    `web.webqtl.correlation.PartialCorrDBPage.getPartialCorrelationsFast`
    function in GeneNetwork1.
    """
    return tuple(sorted(
        samples_from_file.index(good) for good in
        set(samples).intersection(set(samples_from_file))))

def partial_correlations_fast(# pylint: disable=[R0913, R0914]
        samples, primary_vals, control_vals, database_filename,
        fetched_correlations, method: str, correlation_type: str) -> Tuple[
            float, Tuple[float, ...]]:
    """
    Computes partial correlation coefficients using data from a CSV file.

    This is a partial migration of the
    `web.webqtl.correlation.PartialCorrDBPage.getPartialCorrelationsFast`
    function in GeneNetwork1.
    """
    assert method in ("spearman", "pearson")
    with open(f"{TEXTDIR}/{database_filename}", "r") as dataset_file:
        dataset = tuple(dataset_file.readlines())

    good_dataset_samples = good_dataset_samples_indexes(
        samples, parse_csv_line(dataset[0])[1:])

    def __process_trait_names_and_values__(acc, line):
        trait_line = parse_csv_line(line)
        trait_name = trait_line[0]
        trait_data = trait_line[1:]
        if trait_name in fetched_correlations.keys():
            return (
                acc[0] + (trait_name,),
                acc[1] + tuple(
                    trait_data[i] if i in good_dataset_samples else None
                    for i in range(len(trait_data))))
        return acc

    processed_trait_names_values: tuple = reduce(
        __process_trait_names_and_values__, dataset[1:], (tuple(), tuple()))
    all_target_trait_names: Tuple[str, ...] = processed_trait_names_values[0]
    all_target_trait_values: Tuple[float, ...] = processed_trait_names_values[1]

    all_correlations = compute_partial(
        primary_vals, control_vals, all_target_trait_names,
        all_target_trait_values, method)
    ## Line 772 to 779 in GN1 are the cause of the weird complexity in the
    ## return below. Once the surrounding code is successfully migrated and
    ## reworked, this complexity might go away, by getting rid of the
    ## `correlation_type` parameter
    return len(all_correlations), tuple(
        corr + (
            (fetched_correlations[corr[0]],) if correlation_type == "literature"
            else fetched_correlations[corr[0]][0:2])
        for idx, corr in enumerate(all_correlations))

def build_data_frame(
        xdata: Tuple[float, ...], ydata: Tuple[float, ...],
        zdata: Union[
            Tuple[float, ...],
            Tuple[Tuple[float, ...], ...]]) -> pandas.DataFrame:
    """
    Build a pandas DataFrame object from xdata, ydata and zdata
    """
    x_y_df = pandas.DataFrame({"x": xdata, "y": ydata})
    if isinstance(zdata[0], float):
        return x_y_df.join(pandas.DataFrame({"z": zdata}))
    interm_df = x_y_df.join(pandas.DataFrame(
        {"z{}".format(i): val for i, val in enumerate(zdata)}))
    if interm_df.shape[1] == 3:
        return interm_df.rename(columns={"z0": "z"})
    return interm_df

def compute_partial(
        primary_vals, control_vals, target_vals, target_names,
        method: str) -> Tuple[
            Union[
                Tuple[str, int, float, float, float, float], None],
            ...]:
    """
    Compute the partial correlations.

    This is a re-implementation of the
    `web.webqtl.correlation.correlationFunction.determinePartialsByR` function
    in GeneNetwork1.

    This implementation reworks the child function `compute_partial` which will
    then be used in the place of `determinPartialsByR`.

    TODO: moving forward, we might need to use the multiprocessing library to
    speed up the computations, in case they are found to be slow.
    """
    # replace the R code with `pingouin.partial_corr`
    def __compute_trait_info__(target):
        primary = [
            prim for targ, prim in zip(target, primary_vals)
            if targ is not None]
        datafrm = build_data_frame(
            primary,
            [targ for targ in target if targ is not None],
            [cont for i, cont in enumerate(control_vals)
             if target[i] is not None])
        covariates = "z" if datafrm.shape[1] == 3 else [
            col for col in datafrm.columns if col not in ("x", "y")]
        ppc = pingouin.partial_corr(
            data=datafrm, x="x", y="y", covar=covariates, method=method)
        pc_coeff = ppc["r"]

        zero_order_corr = pingouin.corr(
            datafrm["x"], datafrm["y"], method=method)

        if math.isnan(pc_coeff):
            return (
                target[1], len(primary), pc_coeff, 1, zero_order_corr["r"],
                zero_order_corr["p-val"])
        return (
            target[1], len(primary), pc_coeff,
            (ppc["p-val"] if not math.isnan(ppc["p-val"]) else (
                0 if (abs(pc_coeff - 1) < 0.0000001) else 1)),
            zero_order_corr["r"], zero_order_corr["p-val"])

    return tuple(
        __compute_trait_info__(target)
        for target in zip(target_vals, target_names))

def partial_correlations_normal(# pylint: disable=R0913
        primary_vals, control_vals, input_trait_gene_id, trait_database,
        data_start_pos: int, db_type: str, method: str) -> Tuple[
            float, Tuple[float, ...]]:
    """
    Computes the correlation coefficients.

    This is a migration of the
    `web.webqtl.correlation.PartialCorrDBPage.getPartialCorrelationsNormal`
    function in GeneNetwork1.
    """
    def __add_lit_and_tiss_corr__(item):
        if method.lower() == "sgo literature correlation":
            # if method is 'SGO Literature Correlation', `compute_partial`
            # would give us LitCorr in the [1] position
            return tuple(item) + trait_database[1]
        if method.lower() in (
                "tissue correlation, pearson's r",
                "tissue correlation, spearman's rho"):
            # if method is 'Tissue Correlation, *', `compute_partial` would give
            # us Tissue Corr in the [1] position and Tissue Corr P Value in the
            # [2] position
            return tuple(item) + (trait_database[1], trait_database[2])
        return item

    target_trait_names, target_trait_vals = reduce(
        lambda acc, item: (acc[0]+(item[0],), acc[1]+(item[data_start_pos:],)),
        trait_database, (tuple(), tuple()))

    all_correlations = compute_partial(
        primary_vals, control_vals, target_trait_vals, target_trait_names,
        method)

    if (input_trait_gene_id and db_type == "ProbeSet" and method.lower() in (
            "sgo literature correlation", "tissue correlation, pearson's r",
            "tissue correlation, spearman's rho")):
        return (
            len(trait_database),
            tuple(
                __add_lit_and_tiss_corr__(item)
                for idx, item in enumerate(all_correlations)))

    return len(trait_database), all_correlations