<jupyter_start><jupyter_text>PIC-SURE python API use-case: Phenome-Wide analysis on BioData Catalyst studies This notebook is an illustration example about how to query data using the python **PIC-SURE API**. It takes as use-case a simple PheWAS analysis. This notebook is intentionally straightforward, and explanation provided are only aimed at guiding through the PheWAS analysis processus. For a more step-by-step introduction to the python PIC-SURE API, see the `PICSURE-API_101.ipynb` notebook. **Before running this notebook, please be sure to get a user-specific security token. For more information on how to proceed, see the `get_your_token.ipynb` notebook**  Environment set-up  System requirements- Python 3.6 or later- pip package manager- bash interpreter  Installation of external dependencies<jupyter_code>import sys
!{sys.executable} -m pip install -r requirements.txt<jupyter_output><empty_output><jupyter_text>Installing latest python PIC-SURE API libraries from github<jupyter_code>!{sys.executable} -m pip install --upgrade --force-reinstall git+https://github.com/hms-dbmi/pic-sure-python-adapter-hpds.git
!{sys.executable} -m pip install --upgrade --force-reinstall git+https://github.com/hms-dbmi/pic-sure-python-client.git
!{sys.executable} -m pip install --upgrade --force-reinstall git+https://github.com/hms-dbmi/pic-sure-biodatacatalyst-python-adapter-hpds.git
import json
from pprint import pprint

import pandas as pd
import numpy as np 
import matplotlib.pyplot as plt
from scipy import stats

import PicSureClient
import PicSureBdcAdapter

from python_lib.utils import get_multiIndex_variablesDict, joining_variablesDict_onCol
print("NB: This Jupyter Notebook has been written using PIC-SURE API following versions:\n- PicSureBdcAdapter: 1.0.0\n- PicSureClient: 1.1.0")
print("The installed PIC-SURE API libraries versions:\n- PicSureBdcAdapter: {0}\n- PicSureClient: {1}".format(PicSureBdcAdapter.__version__, PicSureClient.__version__))
# Pandas DataFrame display options
pd.set_option("max.rows", 100)

# Matplotlib display parameters
plt.rcParams["figure.figsize"] = (14,8)
font = {'weight' : 'bold',
        'size'   : 12}
plt.rc('font', **font)<jupyter_output><empty_output><jupyter_text>Connecting to a PIC-SURE network<jupyter_code>PICSURE_network_URL = "https://picsure.biodatacatalyst.nhlbi.nih.gov/picsure"
resource_id = "02e23f52-f354-4e8b-992c-d37c8b9ba140"
token_file = "token.<PASSWORD>"
with open(token_file, "r") as f:
    my_token = f.read()
client = PicSureClient.Client()
connection = client.connect(PICSURE_network_URL, my_token)
adapter = PicSureBdcAdapter.Adapter(connection)
resource = adapter.useResource(resource_id)<jupyter_output><empty_output><jupyter_text>PheWAS analysis In a nutshell, this PheWAS analysis follows the subsequent steps:- Retrieving the variables dictionary, using the PIC-SURE API dedicated methods- From the info provided by the dictionary, retrieving the data in an exploitable format through PIC-SURE API calls- Data management- Running univariate tests against every phenotype variables- Accounting for multiple hypotheses testing issue- Plotting resultsThis analysis is conducted using individuals enrolled in the COPDGene Study. Overall goal of this cohort is to detect underlying genetic factors to develop Chronic Obstructive Pulmonary Disease (COPD), and currently includes more than 10,000 individuals ([more information on COPDGene Study](http://www.copdgene.org)).  1. Retrieving variables dictionary from the Database Retrieving variables dictionary for COPDGene related variables.<jupyter_code>copdgene_dic = resource.dictionary().find("Genetic Epidemiology of COPD (COPDGene)").DataFrame()
consent_dic = resource.list_consents()['consent']
# line needed bc dictionaries doesn't return the same column type for categorical and numerical variables
plain_variablesDict = copdgene_dic.transpose().join(consent_dic.transpose(), how="outer").transpose()<jupyter_output><empty_output><jupyter_text>Here we retrieved information about all the available variables of the COPDGene study, as well as information about the individual consents. The Consent variable will be eventually used to filter the patients pertaining to the COPDGene study. We are then combining both dictionaries into a single one.<jupyter_code>plain_variablesDict.head()
resource.dictionary().help()
# Display the variables tree hierarchy from the variables name
variablesDict = get_multiIndex_variablesDict(plain_variablesDict)
variablesDict.iloc[10:20,:]
dic = resource.dictionary().find().DataFrame()<jupyter_output><empty_output><jupyter_text>2. Selecting variables and retrieving data from the database Subseting to keep only the phenotype variables + the variable "affection status" that will be used as the dependent variable for this illustration use-case.<jupyter_code>mask_pheno = variablesDict.index.get_level_values(1) == 'Subject Phenotype'
mask_status = variablesDict.index.get_level_values(2) == 'Affection status'
mask_drop = variablesDict["simplified_name"].isin(["Dbgap_id", "De-identified site code", "A1AD: phenotype/genotype"])
mask_vars = mask_pheno | mask_status
variablesDict = variablesDict.loc[mask_vars & ~ mask_drop,:]
selected_vars = variablesDict.loc[:, "name"].tolist()
pprint(selected_vars[0:5])<jupyter_output><empty_output><jupyter_text>Subsetting only COPDGene records To subset only the subjects from the COPDGene study, we are using the different COPDGene study accession codes, which all begins with 'phs000179...'.<jupyter_code>phs = pd.Series(consent_dic)
mask_phs_copdgene = phs.str.startswith("phs000179").tolist()
phs_copdgene = phs[mask_phs_copdgene]<jupyter_output><empty_output><jupyter_text>Retrieving the data:<jupyter_code>%%capture
query = resource.query()
query.filter().add(consent_dic.index[0], phs_copdgene)
query.select().add(selected_vars)
facts = query.getResultsDataFrame(low_memory=False)
facts.head()
status_var = variablesDict.loc[variablesDict.index.get_level_values(2) == 'Affection status', "name"]
facts = facts.dropna(subset=status_var)\
.set_index("Patient ID")
mask_to_drop = variablesDict["simplified_name"]\
.isin(["Dbgap_id", "De-identified site code", "A1AD: phenotype/genotype"])
variablesDict = variablesDict.loc[~mask_to_drop, :]
var_to_keep = variablesDict.loc[:, "name"]
facts = facts.loc[:, var_to_keep]
print("{0} rows, {1} columns".format(*facts.shape))<jupyter_output><empty_output><jupyter_text>Here is a sample of our dataset, one row per patient and one column par variable:<jupyter_code>facts.tail(5)<jupyter_output><empty_output><jupyter_text>3. Data-management  Selecting variables regarding their types One important step in a PheWAS is to get the distinction between categorical and numerical variables. This distinction is straightforward using the variable dictionary.<jupyter_code>mask_categories = variablesDict.loc[:, "categorical"] == True
categorical_names = variablesDict.loc[mask_categories, "name"].tolist()
continuous_names = variablesDict.loc[~mask_categories, "name"].tolist()<jupyter_output><empty_output><jupyter_text>Selecting the dependent variable to studyMost of PheWAS studies use a genetic variant as the dependent variable to separate the population between cases and controls. However the population doesn't have to be dichotomized using a genetic variant, and any phenotypic variable could be used to run a PheWAS analysis (see for example [*Neuraz et al.*, 2013](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1003405)). Here we will use the **COPD status** as the case-control variable to dichotomize the population in our analysis, and keep only the population subset containing relevant values for the COPD status (i.e. keeping "Case" and "Control" individuals, thus discarding "Other", "Control, Exclusionary Disease", and null values).<jupyter_code>dependent_var_name = variablesDict.loc[variablesDict["simplified_name"] == "Affection status", "name"].values[0]
categorical_names.remove(dependent_var_name)
mask_dependent_var_name = facts[dependent_var_name].isin(["Case", "Control"])
facts = facts.loc[mask_dependent_var_name,:]\
             .astype({dependent_var_name: "category"})
print("Control: {0} individuals\nCase: {1} individuals".format(*facts[dependent_var_name].value_counts().tolist()))<jupyter_output><empty_output><jupyter_text>4. Univariate statistical tests From this point, each variable present in the `facts_dummies` dataset will be tested again the selected dependent variable, (ie presence or absence of COPD). Two different association tests will be carried out according to variables data types: - Mann-Whitney U test for continuous ones- Fisher exact test for categorical ones  LLR result<jupyter_code>from statsmodels.discrete.discrete_model import Logit
independent_names = variablesDict["name"].tolist()
independent_names.remove(dependent_var_name)
dependent_var = facts[dependent_var_name].astype("category").cat.codes
dic_pvalues = {}
simple_index_variablesDict = variablesDict.set_index("name", drop=True)
from scipy.linalg import LinAlgError
from statsmodels.tools.sm_exceptions import PerfectSeparationError
from tqdm import tqdm
for independent_name in tqdm(independent_names, position=0, leave=True):
    matrix = facts.loc[:, [dependent_var_name, independent_name]]\
                  .dropna(how="any")
    if matrix.shape[0] == 0:
        dic_pvalues[independent_name] = np.NaN
        continue
    if simple_index_variablesDict.loc[independent_name, "categorical"]:
        matrix = pd.get_dummies(matrix,
                                columns=[independent_name],
                                drop_first=False)\
                    .iloc[:, 0:-1]
    dependent_var = matrix[dependent_var_name].cat.codes
    independent_var = matrix.drop(dependent_var_name, axis=1)\
                            .assign(intercept = 1)
    model = Logit(dependent_var, independent_var)
    try:
        results = model.fit(disp=0)
        dic_pvalues[independent_name] = results.llr_pvalue
    except (LinAlgError, PerfectSeparationError) as e:
        dic_pvalues[independent_name] = np.NaN<jupyter_output><empty_output><jupyter_text>p-values distribution (univariate tests)<jupyter_code>pd.Series([v for v in dic_pvalues.values()]).plot.hist(bins=30)
plt.suptitle("Distribution of individual p-values",
             weight="bold",
            fontsize=15)<jupyter_output><empty_output><jupyter_text>5. Multiple hypotheses testing correction: Bonferroni Method In order to handle the multiple testing problem (increase in the probability of getting statistically significant associations), we will use the Bonferroni correction method. Although many other multiple-comparison corrections exist, Bonferroni is the most straightforward to use, because it doesn't require assumptions about variables correlation. Other PheWAS analysis also use False Discovery Rate controlling procedures ([see](https://en.wikipedia.org/wiki/False_discovery_rate)). In a nutshell, Bonferonni allows to calculate a corrected "statistical significant threshold" according to the number of test performed. Every p-value below this threshold will be deemed statistically significant.<jupyter_code>%%capture
# Merging pvalues from different tests
df_pvalues = pd.DataFrame.from_dict(dic_pvalues, orient="index", columns=["pvalues"])\
.rename_axis("name")\
.reset_index(drop=False)

# Adding pvalues results as a new column to variablesDict
variablesDict = joining_variablesDict_onCol(variablesDict,
                                              df_pvalues,
                                              left_col="name",
                                              right_col="name")

adjusted_alpha = 0.05/len(variablesDict["pvalues"])
variablesDict["p_adj"] = variablesDict["pvalues"] / len(variablesDict["pvalues"])
variablesDict['log_p'] = -np.log10(variablesDict['pvalues'])
variablesDict = variablesDict.sort_index()
variablesDict["group"] = variablesDict.reset_index(level=2)["level_2"].values
print("Bonferonni adjusted significance threshold: {0:.2E}".format(adjusted_alpha))<jupyter_output><empty_output><jupyter_text>6. Result visualisations: Manhattan plot Manhattan plot is the classical way to plot the results of a PheWAS analysis. It plots every tested phenotypical variables on the X-axis, against their *-log(pvalue)* on the Y-axis. The horizontal line represents the adjusted significance level threshold.<jupyter_code>mask = variablesDict["pvalues"].isna()
df_results = variablesDict.loc[~mask,:].copy().replace([np.inf, -np.inf], np.nan)
df_results = df_results.loc[~df_results["log_p"].isna().values,:]

#### Specific adjustment to make this specific plot looks nicer
####### to adapt when changing data or dependent variable
df_results = df_results.replace({"TLC": "Spirometry",
                                 "New Gold Classification": "Quantitative Analysis", 
                  "Other": "Demographics",
                                "Eligibility Form": "Sociodemography and Administration"})
group_order={'6MinWalk': 0,
 'CT Acquisition Parameters': 1,
 'CT Assessment Scoresheet': 2,
 'Demographics and Physical Characteristics': 3,
 'Longitudinal Analysis': 5,
 'Medical History': 4,
 'Medication History': 13,
 'Quantitative Analysis': 9,
 'Respiratory Disease': 6,
 'SF-36 Health Survey': 11,
 'Sociodemography and Administration': 12,
 'Spirometry': 7,
 'VIDA': 15}
df_results["group_order"] = df_results["group"].replace(group_order)
df_results = df_results.sort_values("group_order", ascending=True)
df_results["simplified_name"] = df_results["simplified_name"].str.replace("[0-9]+[A-z]*", "").to_frame()
###


fig = plt.figure()
ax = fig.add_subplot(111)
colors = plt.get_cmap('Set1')
x_labels = []
x_labels_pos = []

y_lims = (0, df_results["log_p"].max(skipna=True) + 50)
threshold_top_values = df_results["log_p"].sort_values(ascending=False)[0:6][-1]

df_results["ind"] = np.arange(1, len(df_results)+1)
df_grouped = df_results.groupby(('group'))
for num, (name, group) in enumerate(df_grouped):
    group.plot(kind='scatter', x='ind', y='log_p',color=colors.colors[num % len(colors.colors)], ax=ax, s=20)
    x_labels.append(name)
    x_labels_pos.append((group['ind'].iloc[-1] - (group['ind'].iloc[-1] - group['ind'].iloc[0])/2)) # Set label in the middle
    for n, row in group.iterrows():
        if row["log_p"] > threshold_top_values:
            ax.text(row['ind'] + 3, row["log_p"] + 0.05, row["simplified_name"], rotation=0, alpha=1, size=8, color="black")
                
ax.set_xticks(x_labels_pos)
ax.set_xticklabels(x_labels)
ax.set_xlim([0, len(df_results) +1])
ax.set_ylim(y_lims)
ax.set_ylabel('-log(p-values)', style="italic")
ax.set_xlabel('Phenotypes', fontsize=15)
ax.axhline(y=-np.log10(adjusted_alpha), linestyle=":", color="black", label="Adjusted Threshold")
plt.xticks(fontsize = 9,rotation=90)
plt.yticks(fontsize = 8)
plt.title("Statistical Association Between COPD Status and Phenotypes", 
          loc="center",
          style="oblique", 
          fontsize = 20,
         y=1)
xticks = ax.xaxis.get_major_ticks()
xticks[0].set_visible(False)
handles, labels = ax.get_legend_handles_labels()
plt.legend(handles = handles, labels = labels, loc = "upper left")
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10864585378624256, 'fraction_numerical': 0.01191177433078081, 'mean_word_length': 5.513213981244672, 'pattern_counts': {'":': 5, '<': 62, '<?xml version=': 0, '>': 63, 'https://': 6, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '459397', 'license': 'Apache-2.0', 'n_tokens_mistral': 4768, 'n_tokens_neox': 4489, 'n_words': 1239, 'path': 'NHLBI_BioData_Catalyst/python/3_PheWAS.ipynb', 'repo_name': 'hms-dbmi/Access-to-Data-using-PIC-SURE-API', 'size': '15272'}