<jupyter_start><jupyter_text>Compute perturbagen signatures from gold signature consensi<jupyter_code>import bz2
import json

import requests
import pandas
import sqlite3

import l1000<jupyter_output><empty_output><jupyter_text>Entrez Gene to Symbol mapping<jupyter_code># Read symbol to entrez_gene_id mapping
url = 'https://github.com/dhimmel/entrez-gene/raw/a7362748a34211e5df6f2d185bb3246279760546/data/symbol-map.json'
symbol_map = json.loads(requests.get(url).text)

# Read dataframe of entrez gene records
url = 'https://github.com/dhimmel/entrez-gene/raw/a7362748a34211e5df6f2d185bb3246279760546/data/genes-human.tsv'
entrez_gene_df = pandas.read_table(url)
renamer = {
    'GeneID': 'entrez_gene_id',
    'Symbol': 'symbol',
    'type_of_gene': 'type_of_gene',
    'description': 'description',
}
entrez_gene_df = entrez_gene_df.rename(columns=renamer)[list(renamer.values())]
entrez_gene_df.entrez_gene_id = entrez_gene_df.entrez_gene_id.astype(str)<jupyter_output><empty_output><jupyter_text>Read probe information<jupyter_code># construct probe_df
probe_df = pandas.read_table('data/geneinfo/geneinfo.tsv.gz')
probe_to_gene = dict(zip(probe_df.pr_id, probe_df.pr_gene_id))

# Landmark probes
landmark_probe_df = probe_df[probe_df.pr_pool_id.str.startswith('epsilon').fillna(False)]
print('landmark genes', len(landmark_probe_df))

# BIGS + Landmark probes
probe_df = probe_df[probe_df.is_bing.fillna(False)]
# If a gene is an epsilon landmark, do not include imputed probes
probe_df.query("pr_gene_id not in @landmark_probe_df.pr_gene_id or pr_id in @landmark_probe_df.pr_id")
print('best inferred gene set size', len(probe_df))

probe_df.head(2)
# Create a dataset of the genes represented by probes
gene_df = probe_df.groupby('pr_gene_id').apply(
    lambda df: pandas.Series({'status': 'measured' if df.pr_pool_id.iloc[0].startswith('epsilon') else 'imputed'})
)
gene_df.index.name = 'entrez_gene_id'
gene_df = gene_df.reset_index()
gene_df = gene_df.merge(entrez_gene_df).sort_values('entrez_gene_id')
gene_df.to_csv('data/consensi/genes.tsv', index=False, sep='\t')
gene_df.head(2)
# Number of genes after converting from probe to gene
gene_df.status.value_counts()
# Filter BING probes for valid entrez gene IDs
probe_df = probe_df.query("pr_gene_id in @gene_df.entrez_gene_id")

# Probes in BING
probe_df.pr_pool_id.value_counts()<jupyter_output><empty_output><jupyter_text>Preparations<jupyter_code>def run_consensi(sig_expr_df, pert_to_sigs, name):
    """Compute consensi signatures"""
    pert_expr_df = l1000.get_consensus_signatures(sig_expr_df, pert_to_sigs, weighting_subset=landmark_probe_df.pr_id)
    pert_expr_df = l1000.probes_to_genes(pert_expr_df, probe_to_gene)
    pert_expr_df = pert_expr_df.transpose()
    pert_expr_df.index.name = 'perturbagen'
    print(pert_expr_df.shape)
    path = 'data/consensi/consensi-{}.tsv.bz2'.format(name)
    with bz2.BZ2File(path, 'w') as write_file:
        pert_expr_df.reset_index().to_csv(write_file, sep='\t', index=False, float_format='%.3f')
    return pert_expr_df
# open database connection
connection = sqlite3.connect('data/l1000.db')
%%time

# get all gold signatures
query = "SELECT sigs.sig_id FROM sigs WHERE sigs.is_gold = 1"
sigs = pandas.read_sql(query, connection).sig_id.tolist()

# get probes and extract signatures
probes = probe_df.pr_id.tolist()
path = 'download/modzs.gctx'
sig_expr_df = l1000.extract_from_gctx(path, probes, sigs)<jupyter_output>CPU times: user 56min 42s, sys: 9min 9s, total: 1h 5min 52s
Wall time: 1h 6min 47s<jupyter_text>drugbank consensi<jupyter_code>query = """
SELECT unichem.resource_id AS drugbank_id, perts.pert_id, sigs.sig_id, sigs.is_gold
FROM unichem, perts, sigs
WHERE unichem.resource = 'drugbank'
AND unichem.pert_uid = perts.pert_uid
AND sigs.pert_id = perts.pert_id
"""
sig_df = pandas.read_sql(query, connection)
sig_df = sig_df.query("is_gold == 1")
%%time
pert_to_sigs = {k: g['sig_id'].tolist() for k, g in sig_df.groupby('drugbank_id')}
pert_expr_df = run_consensi(sig_expr_df, pert_to_sigs, name='drugbank')<jupyter_output>(1170, 7467)
CPU times: user 1h 54min 57s, sys: 1min 16s, total: 1h 56min 14s
Wall time: 1h 56min 28s<jupyter_text>knockdown consensi<jupyter_code>query = """
SELECT perts.pert_id, perts.pert_iname, perts.pert_type, sigs.sig_id
FROM perts, sigs
WHERE sigs.pert_id = perts.pert_id
AND pert_type = 'trt_sh'
AND sigs.is_gold = 1
"""
sig_df = pandas.read_sql(query, connection)
sig_df['pertubation_entrez_gene_id'] = sig_df.pert_iname.map(symbol_map.get)
sig_df.head(2)
%%time
# Condense to perturbagens
pert_to_sigs = {int(k): g['sig_id'].tolist() for k, g in sig_df.groupby('pertubation_entrez_gene_id')}
pert_expr_df = run_consensi(sig_expr_df, pert_to_sigs, name='knockdown')<jupyter_output>(4326, 7467)
CPU times: user 5h 38min 35s, sys: 1min 5s, total: 5h 39min 40s
Wall time: 5h 40min 16s<jupyter_text>overexpression consensi<jupyter_code>query = """
SELECT perts.pert_id, perts.pert_iname, perts.pert_type, sigs.sig_id
FROM perts, sigs
WHERE sigs.pert_id = perts.pert_id
AND pert_type = 'trt_oe'
AND sigs.is_gold = 1
"""
sig_df = pandas.read_sql(query, connection)
sig_df['pertubation_entrez_gene_id'] = sig_df.pert_iname.map(symbol_map.get)
sig_df.head(2)
%%time
# Condense to perturbagens
pert_to_sigs = {int(k): g['sig_id'].tolist() for k, g in sig_df.groupby('pertubation_entrez_gene_id')}
pert_expr_df = run_consensi(sig_expr_df, pert_to_sigs, name='overexpression')<jupyter_output>(2413, 7467)
CPU times: user 2h 25min 22s, sys: 6.9 s, total: 2h 25min 29s
Wall time: 2h 25min 42s<jupyter_text>consensi for all pert_ids<jupyter_code>query = """
SELECT perts.pert_id, perts.pert_iname, perts.pert_type, sigs.sig_id
FROM perts, sigs
WHERE sigs.pert_id = perts.pert_id
AND sigs.is_gold = 1
"""
sig_df = pandas.read_sql(query, connection)
sig_df.head(2)
%%time
# Condense to perturbagens
pert_to_sigs = {k: g['sig_id'].tolist() for k, g in sig_df.groupby('pert_id')}
pert_expr_df = run_consensi(sig_expr_df, pert_to_sigs, name='pert_id')
pert_expr_df.head(2)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1199009083402147, 'fraction_numerical': 0.031544178364987616, 'mean_word_length': 6.87516254876463, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 29, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '529108', 'license': 'CC0-1.0', 'n_tokens_mistral': 2641, 'n_tokens_neox': 2424, 'n_words': 453, 'path': 'consensi.ipynb', 'repo_name': 'dhimmel/lincs', 'size': '6055'}