aboutsummaryrefslogtreecommitdiff
path: root/gn2/wqflask/api/correlation.py
blob: 090d13ac59891ed6548dab72adda05fd4128efa9 (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
import collections
import scipy
import numpy

from gn2.base import data_set
from gn2.base.trait import create_trait, retrieve_sample_data
from gn2.utility import corr_result_helpers
from gn2.utility.tools import get_setting
from gn2.wqflask.correlation import correlation_functions
from gn2.wqflask.database import database_connection

def do_correlation(start_vars):
    if 'db' not in start_vars:
        raise ValueError("'db' not found!")
    if 'target_db' not in start_vars:
        raise ValueError("'target_db' not found!")
    if 'trait_id' not in start_vars:
        raise ValueError("'trait_id' not found!")

    this_dataset = data_set.create_dataset(dataset_name=start_vars['db'])
    target_dataset = data_set.create_dataset(
        dataset_name=start_vars['target_db'])
    this_trait = create_trait(dataset=this_dataset,
                              name=start_vars['trait_id'])
    this_trait = retrieve_sample_data(this_trait, this_dataset)

    corr_params = init_corr_params(start_vars)

    corr_results = calculate_results(
        this_trait, this_dataset, target_dataset, corr_params)

    final_results = []
    for _trait_counter, trait in enumerate(list(corr_results.keys())[:corr_params['return_count']]):
        if corr_params['type'] == "tissue":
            [sample_r, num_overlap, sample_p, symbol] = corr_results[trait]
            result_dict = {
                "trait": trait,
                "sample_r": sample_r,
                "#_strains": num_overlap,
                "p_value": sample_p,
                "symbol": symbol
            }
        elif corr_params['type'] == "literature" or corr_params['type'] == "lit":
            [gene_id, sample_r] = corr_results[trait]
            result_dict = {
                "trait": trait,
                "sample_r": sample_r,
                "gene_id": gene_id
            }
        else:
            [sample_r, sample_p, num_overlap] = corr_results[trait]
            result_dict = {
                "trait": trait,
                "sample_r": sample_r,
                "#_strains": num_overlap,
                "p_value": sample_p
            }
        final_results.append(result_dict)
    return final_results


def calculate_results(this_trait, this_dataset, target_dataset, corr_params):
    corr_results = {}

    target_dataset.get_trait_data()

    if corr_params['type'] == "tissue":
        trait_symbol_dict = this_dataset.retrieve_genes("Symbol")
        corr_results = do_tissue_correlation_for_all_traits(
            this_trait, trait_symbol_dict, corr_params)
        sorted_results = collections.OrderedDict(sorted(list(corr_results.items()),
                                                        key=lambda t: -abs(t[1][1])))
    # ZS: Just so a user can use either "lit" or "literature"
    elif corr_params['type'] == "literature" or corr_params['type'] == "lit":
        trait_geneid_dict = this_dataset.retrieve_genes("GeneId")
        corr_results = do_literature_correlation_for_all_traits(
            this_trait, this_dataset, trait_geneid_dict, corr_params)
        sorted_results = collections.OrderedDict(sorted(list(corr_results.items()),
                                                        key=lambda t: -abs(t[1][1])))
    else:
        for target_trait, target_vals in list(target_dataset.trait_data.items()):
            result = get_sample_r_and_p_values(
                this_trait, this_dataset, target_vals, target_dataset, corr_params['type'])
            if result is not None:
                corr_results[target_trait] = result

        sorted_results = collections.OrderedDict(
            sorted(list(corr_results.items()), key=lambda t: -abs(t[1][0])))

    return sorted_results


def do_tissue_correlation_for_all_traits(this_trait, trait_symbol_dict, corr_params, tissue_dataset_id=1):
    # Gets tissue expression values for the primary trait
    primary_trait_tissue_vals_dict = correlation_functions.get_trait_symbol_and_tissue_values(
        symbol_list=[this_trait.symbol])

    if this_trait.symbol.lower() in primary_trait_tissue_vals_dict:
        primary_trait_tissue_values = primary_trait_tissue_vals_dict[this_trait.symbol.lower(
        )]

        corr_result_tissue_vals_dict = correlation_functions.get_trait_symbol_and_tissue_values(
            symbol_list=list(trait_symbol_dict.values()))

        tissue_corr_data = {}
        for trait, symbol in list(trait_symbol_dict.items()):
            if symbol and symbol.lower() in corr_result_tissue_vals_dict:
                this_trait_tissue_values = corr_result_tissue_vals_dict[symbol.lower(
                )]

                result = correlation_functions.cal_zero_order_corr_for_tiss(primary_trait_tissue_values,
                                                                            this_trait_tissue_values,
                                                                            corr_params['method'])

                tissue_corr_data[trait] = [
                    result[0], result[1], result[2], symbol]

        return tissue_corr_data


def do_literature_correlation_for_all_traits(this_trait, target_dataset, trait_geneid_dict, corr_params):
    input_trait_mouse_gene_id = convert_to_mouse_gene_id(
        target_dataset.group.species.lower(), this_trait.geneid)

    lit_corr_data = {}
    for trait, gene_id in list(trait_geneid_dict.items()):
        mouse_gene_id = convert_to_mouse_gene_id(
            target_dataset.group.species.lower(), gene_id)

        if mouse_gene_id and str(mouse_gene_id).find(";") == -1:
            result = ""
            with database_connection(get_setting("SQL_URI")) as conn:
                with conn.cursor() as cursor:
                    cursor.execute(
                        ("SELECT value FROM LCorrRamin3 "
                         "WHERE GeneId1=%s AND GeneId2=%s"),
                        (mouse_gene_id,
                         input_trait_mouse_gene_id))
                    result = cursor.fetchone()
                    if not result:
                        cursor.execute(
                            ("SELECT value FROM LCorrRamin3 "
                             "WHERE GeneId2=%s AND GeneId1=%s"),
                            (mouse_gene_id,
                             input_trait_mouse_gene_id))
                        result = cursor.fetchone()
            if result:
                lit_corr = result[0]
                lit_corr_data[trait] = [gene_id, lit_corr]
            else:
                lit_corr_data[trait] = [gene_id, 0]
        else:
            lit_corr_data[trait] = [gene_id, 0]

    return lit_corr_data


def get_sample_r_and_p_values(this_trait, this_dataset, target_vals, target_dataset, type):
    """
    Calculates the sample r (or rho) and p-value

    Given a primary trait and a target trait's sample values,
    calculates either the pearson r or spearman rho and the p-value
    using the corresponding scipy functions.
    """

    this_trait_vals = []
    shared_target_vals = []
    for i, sample in enumerate(target_dataset.group.samplelist):
        if sample in this_trait.data:
            this_sample_value = this_trait.data[sample].value
            target_sample_value = target_vals[i]
            this_trait_vals.append(this_sample_value)
            shared_target_vals.append(target_sample_value)

    this_trait_vals, shared_target_vals, num_overlap = corr_result_helpers.normalize_values(
        this_trait_vals, shared_target_vals)

    if type == 'pearson':
        sample_r, sample_p = scipy.stats.pearsonr(
            this_trait_vals, shared_target_vals)
    else:
        sample_r, sample_p = scipy.stats.spearmanr(
            this_trait_vals, shared_target_vals)

    if num_overlap > 5:
        if numpy.isnan(sample_r):
            return None
        else:
            return [sample_r, sample_p, num_overlap]


def convert_to_mouse_gene_id(species=None, gene_id=None):
    """If the species is rat or human, translate the gene_id to the mouse geneid

    If there is no input gene_id or there's no corresponding mouse gene_id, return None

    """
    if not gene_id:
        return None

    mouse_gene_id = None
    with database_connection(get_setting("SQL_URI")) as conn:
        with conn.cursor() as cursor:
            if species == 'mouse':
                mouse_gene_id = gene_id
            elif species == 'rat':
                cursor.execute(
                    ("SELECT mouse FROM GeneIDXRef "
                     "WHERE rat=%s"), gene_id)
                result = cursor.fetchone()
                if result:
                    mouse_gene_id = result[0]
            elif species == 'human':
                cursor.execute(
                    "SELECT mouse FROM GeneIDXRef "
                    "WHERE human=%s", gene_id)
                result = cursor.fetchone()
                if result:
                    mouse_gene_id = result[0]
    return mouse_gene_id


def init_corr_params(start_vars):
    method = "pearson"
    if 'method' in start_vars:
        method = start_vars['method']

    type = "sample"
    if 'type' in start_vars:
        type = start_vars['type']

    return_count = 500
    if 'return_count' in start_vars:
        assert(start_vars['return_count'].isdigit())
        return_count = int(start_vars['return_count'])

    corr_params = {
        'method': method,
        'type': type,
        'return_count': return_count
    }

    return corr_params