<jupyter_start><jupyter_text>Description This notebooks will show the structure of the main data matrices in PhenoPLIER, and will guide you in analyzing gene associations for a particular trait: neutrophil counts, which is presented in the [manuscript](https://greenelab.github.io/phenoplier_manuscript/phenoplier-an-integration-framework-based-on-gene-co-expression-patterns) in Figure 1c.  Modules<jupyter_code>import tempfile

import numpy as np
from scipy import stats
import pandas as pd

from entity import Trait, Gene
import conf<jupyter_output><empty_output><jupyter_text>Load gene module-gene membership matrix (matrix Z) Here we load the gene module-gene membership matri, or "latent variables loadings matrix" (from the terminology of the [MultiPLIER article](https://doi.org/10.1016/j.cels.2019.04.003)).<jupyter_code>matrix_z = pd.read_pickle(conf.MULTIPLIER["MODEL_Z_MATRIX_FILE"])
matrix_z.shape
matrix_z.head()<jupyter_output><empty_output><jupyter_text>As you can see, this matrix Z contains the membership value for each gene across all LVs (or gene modules).A value of zero means that the gene does not belong to that LV, whereas a larger value represents how strongly that gene belongs to the LV.A group of genes that belong to the same LV represent a gene-set that has a similar expression profile across a set of tissues or cell types.We'll cover this in more detail in the next notebook (`02-LV_cell_types-...`).  Load information about LV alignment with pathways LV in matrix Z can represent a group of genes that align well with prior pathways (or prior knowledge) or be "novel" in the sense that the combination of genes do not represent a known unit but was found the PLIER when factorizing the recount2 data (see the MultiPLIER article for more details).Here we load that information, where for each LV and pathway, we have a p-value and area under the curve (AUC) that indicate how well the LV aligns to that pathway.<jupyter_code>lv_metadata = pd.read_pickle(conf.MULTIPLIER["MODEL_SUMMARY_FILE"])
lv_metadata.shape
lv_metadata.head()<jupyter_output><empty_output><jupyter_text>Load gene associations from PhenomeXcan Now we load the gene association across ~4,000 traits from [PhenomeXcan](https://doi.org/10.1126/sciadv.aba2083).The file we load here are the Summary-MultiXcan (or S-MultiXcan) results, essentially a p-value for each gene-trait pair.In the notebook I refer to these results generically as "TWAS results", meaning that we have gene-trait associations.All these TWAS results were derived solely from GWAS summary stats, so you can also generate yours relatively easily by using [S-MultiXcan](https://doi.org/10.1371/journal.pgen.1007889).<jupyter_code>phenomexcan_df = pd.read_pickle(conf.PHENOMEXCAN["SMULTIXCAN_MASHR_PVALUES_FILE"])
phenomexcan_df.shape
phenomexcan_df.head()<jupyter_output><empty_output><jupyter_text>Convert gene Ensembl IDs to symbols:<jupyter_code>phenomexcan_df = phenomexcan_df.rename(index=Gene.GENE_ID_TO_NAME_MAP)
phenomexcan_df = phenomexcan_df.loc[~phenomexcan_df.index.duplicated()]
phenomexcan_df.head()<jupyter_output><empty_output><jupyter_text>Keep genes present in our matrix Z only:<jupyter_code>common_genes = phenomexcan_df.index.intersection(matrix_z.index)
display(common_genes)
phenomexcan_df = phenomexcan_df.loc[common_genes]
phenomexcan_df.shape
matrix_z = matrix_z.loc[common_genes]
matrix_z.shape<jupyter_output><empty_output><jupyter_text>Take a look at genes associated with neutrophil counts Below I search the PhenomeXcan results to find traits related to "neutrophils".<jupyter_code>phenomexcan_df.columns[phenomexcan_df.columns.str.lower().str.contains("neutrophil")]<jupyter_output><empty_output><jupyter_text>For this demo, I select the the "neutrophil count" TWAS derived from the GWAS performed by [Astle et. al](https://doi.org/10.1016/j.cell.2016.10.042). Below you can see the sample size:<jupyter_code>trait_code = "Astle_et_al_2016_Neutrophil_count"
t = Trait.get_trait(full_code=trait_code)
display(f"{trait_code} - sample size: {t.n}")
traits_df = phenomexcan_df[[trait_code]].dropna()
traits_df.shape
traits_df.head()<jupyter_output><empty_output><jupyter_text>Here I quickly show the data summary for this trait's gene associations:<jupyter_code>traits_df.apply(lambda x: -np.log10(x)).describe()<jupyter_output><empty_output><jupyter_text>Make sure we don't have missing values or NaN.<jupyter_code>assert not traits_df.isna().any().any()
assert not np.isinf(traits_df).any().any()<jupyter_output><empty_output><jupyter_text>Convert p-values to z-scores This converts a p-value to a scalar > 0, where higher positive values mean stronger association, and lower values close to zero mean weaker associations. Check our manuscript for more details on this.<jupyter_code>traits_zscores = pd.DataFrame(
    data=np.abs(stats.norm.ppf(traits_df / 2)),
    index=traits_df.index.copy(),
    columns=traits_df.columns.copy(),
)
assert not traits_zscores.isna().any().any()
assert not np.isinf(traits_zscores).any().any()
traits_zscores.head()
traits_zscores.describe()<jupyter_output><empty_output><jupyter_text>Analysis of a neutrophil-termed LV Let's take as an example an LV that was previously analyzed in the MultiPLIER study, which we identify as `LV603`. This LV aligns well with pathways related to neutrophils, as you can see below.<jupyter_code>lv_metadata[
    (lv_metadata["LV index"] == "603") & (lv_metadata["FDR"] < 0.05)
].sort_values("FDR")<jupyter_output><empty_output><jupyter_text>Let's see which genes more strongly belong to LV603 (the numbers are the gene weights in this LV):<jupyter_code>lv603_top_genes = matrix_z["LV603"].sort_values(ascending=False)
display(lv603_top_genes.head(20))<jupyter_output><empty_output><jupyter_text>Are these top genes associated with our neutrophil count?<jupyter_code>traits_zscores.loc[lv603_top_genes.index].head(20)<jupyter_output><empty_output><jupyter_text>It seems so. But what about the rest of the genes? They might be also strongly associated.Let's take a random sample:<jupyter_code>traits_zscores.sample(n=20, random_state=0)<jupyter_output><empty_output><jupyter_text>They do not seem as high as those within the top genes in LV603.If we compute the correlation between LV603 gene weights (`lv603_top_genes`) and gene associations for neutrophil counts (`traits_zscores`) we get this:<jupyter_code>stats.pearsonr(
    traits_zscores.loc[lv603_top_genes.index].iloc[:, 0].to_numpy(),
    lv603_top_genes.to_numpy(),
)<jupyter_output><empty_output><jupyter_text>Although the correlation is significant (`8.16e-10`) and the slope positive (we are interested only in genes at the top of the LV), we need to account for correlated predicted expression from the TWAS models (for example, if the expression of two genes at the top of the LV is correlated that would invalidate our test).We have a class implemented in Python that computes this, as shown below.  Association between an LV and a trait You can use our class `gls.GLSPhenoplier` to compute an association between an LV and a trait, in a similar way we did it in the previous cells by correlating LV603 gene weight's and neutrophil count gene associations (transformed to z-scores). However, `gls.GLSPhenoplier` takes into account correlations between the predicted expression of genes by using a generalized least squares (GLS) model.<jupyter_code>from gls import GLSPhenoplier<jupyter_output><empty_output><jupyter_text>We need to save our trait gene association to a pickle file before moving on. Here I generate a temporary file path and save:<jupyter_code>with tempfile.NamedTemporaryFile(suffix=".pkl") as f:
    output_filepath = f.name
    display(output_filepath)

traits_zscores.to_pickle(output_filepath)<jupyter_output><empty_output><jupyter_text>Now select an LV id and the trait id, as shown below:<jupyter_code>lv_code = "LV603"

phenotype_code = traits_zscores.columns[0]
display(phenotype_code)<jupyter_output><empty_output><jupyter_text>Run the model:<jupyter_code>gls_model = GLSPhenoplier(
    smultixcan_result_set_filepath=output_filepath,
)

gls_model.fit_named(lv_code, phenotype_code)<jupyter_output><empty_output><jupyter_text>It will take a few seconds. You can show the entire table of results from the [statsmodels's GLS model](https://www.statsmodels.org/stable/generated/statsmodels.regression.linear_model.GLS.html):<jupyter_code>res = gls_model.results
print(res.summary())<jupyter_output>GLS Regression Results                            
==============================================================================
Dep. Variable:              phenotype   R-squared:                       0.005
Model:                            GLS   Adj. R-squared:                  0.005
Method:                 Least Squares   F-statistic:                     30.34
Date:                Fri, 08 Oct 2021   Prob (F-statistic):           3.76e-08
Time:                        15:58:08   Log-Likelihood:                -9329.5
No. Observations:                6450   AIC:                         1.866e+04
Df Residuals:                    6448   BIC:                         1.868e+04
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|[...]<jupyter_text>This p-value (`3.76e-08`) is a two-sided test on the LV coefficient (`0.0698`).You can see that the p-value is slightly less significant than the p-value from the Pearson correlation that we computed before. This is how you can access the model's estimated parameters.These are the coefficients:<jupyter_code>res.params<jupyter_output><empty_output><jupyter_text>The one-sided p-values (the ones we used in our manuscript, since we are only interested in the top genes of an LV):<jupyter_code>res.pvalues_onesided<jupyter_output><empty_output><jupyter_text>And the two-sided p-values (in case you need them or are interested in a different hypothesis):<jupyter_code>res.pvalues<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09732072498029945, 'fraction_numerical': 0.01901103230890465, 'mean_word_length': 4.91322073383809, 'pattern_counts': {'":': 0, '<': 97, '<?xml version=': 0, '>': 98, 'https://': 6, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '24', 'dir': 'jupyter-structured-clean-dedup', 'id': '592350', 'license': 'BSD-2-Clause-Patent', 'n_tokens_mistral': 3369, 'n_tokens_neox': 3045, 'n_words': 953, 'path': 'nbs/99_demo/01-LV_trait_association-neutrophil_counts.ipynb', 'repo_name': 'greenelab/phenoplier', 'size': '10152'}