<jupyter_start><jupyter_text>LINCS L1000 significant dysregulation<jupyter_code>import pandas
import scipy.stats
import numpy
from statsmodels.sandbox.stats import multicomp

%matplotlib inline
# Read gene info
gene_df = pandas.read_table('data/consensi/genes.tsv', dtype={'entrez_gene_id': str})
gene_df.head(2)
def process_matrix_df(z_matrix_df):
    """
    Take a perturbagen by gene dataframe and extract signficantly dysregulated pairs.
    Returns `signif_df` with a row per dysregulated pair and `summary_df` which counts the number
    of dysregulated genes per perturbation.
    """
    melt_df = pandas.melt(z_matrix_df.reset_index(), id_vars='perturbagen', var_name = 'entrez_gene_id', value_name = 'z_score')
    melt_df = melt_df.merge(gene_df[['entrez_gene_id', 'symbol', 'status']])
    signif_df = melt_df.groupby(['perturbagen', 'status']).apply(get_significance).reset_index(drop=True)
    signif_df = signif_df.sort_values(['perturbagen', 'symbol'])
    summary_df = signif_df.groupby(['perturbagen', 'direction', 'status']).apply(lambda df: pandas.Series({'count': len(df)})).reset_index()
    summary_df = summary_df.pivot_table('count', 'perturbagen', ['direction', 'status']).fillna(0).astype(int).reset_index()
    summary_df.columns = ['-'.join(col).strip('-') for col in summary_df.columns.values]
    return signif_df, summary_df
def get_significance(df):
    """
    Get signficant perturbagen-gene pairs using a Bonferroni adjustment.
    Different modes for measured and imputed genes.
    """
    p_values = 2 * scipy.stats.norm.cdf(-df.z_score.abs())
    if all(df.status == 'measured'):
        alpha = 0.05
        max_diffex = len(df)
    elif all(df.status == 'imputed'):
        alpha = 0.05
        max_diffex = 1000
    else:
        raise ValueError('Invalid status')
    reject, pvals_corrected, alpha_c_sidak, alpha_c_bonf = multicomp.multipletests(p_values, alpha=alpha, method='bonferroni')
    df['direction'] = df.z_score.map(lambda x: 'up' if x > 0 else 'down')
    df['nlog10_bonferroni_pval'] = -numpy.log10(pvals_corrected)
    df = df[reject]
    df = df.sort_values('nlog10_bonferroni_pval', ascending=False).iloc[:max_diffex, :]
    return df
# Iterate through all consensi and process them
for pert_kind in ['drugbank', 'knockdown', 'overexpression', 'pert_id']:
    print(pert_kind)
    path = 'data/consensi/consensi-{}.tsv.bz2'.format(pert_kind)
    z_matrix_df = pandas.read_table(path, index_col=0)
    signif_df, summary_df = process_matrix_df(z_matrix_df)
    path = 'data/consensi/signif/dysreg-{}.tsv'.format(pert_kind)
    signif_df.to_csv(path, index=False, sep='\t', float_format='%.3f')
    path = 'data/consensi/signif/dysreg-{}-summary.tsv'.format(pert_kind)
    summary_df.to_csv(path, index=False, sep='\t')
signif_df.head()
summary_df.head()
# compress `dysreg-pert_id.tsv` to reduce file size
! gzip data/consensi/signif/dysreg-pert_id.tsv<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1379544682296976, 'fraction_numerical': 0.009174311926605505, 'mean_word_length': 4.8297029702970296, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '606581', 'license': 'CC0-1.0', 'n_tokens_mistral': 1098, 'n_tokens_neox': 1038, 'n_words': 200, 'path': 'significance.ipynb', 'repo_name': 'dhimmel/lincs', 'size': '2943'}