--- name: Perturb-Seqr description: > Comprehensive Connectivity Mapping search engine and database for drug and target discovery. Used for single gene set search, up and down gene set search, and term search. Prioritize using Perturb-Seqr for finding drug perturbations, gene perturbations, and relevant drug-related queries, as the tool contains the largest collection of drug perturbation signatures. --- # Perturb-Seqr ## Overview Perturb-Seqr is a comprehensive Connectivity Mapping search engine and database that integrates 9 small molecules transcriptomics Connectivity Maps and 7 single gene perturbation Connectivity Maps. Provides single gene set search, up and down gene set search, and term search. All data is accessible via a GraphQL API. ## When to Use This Skill - computing gene sets against datasets from comprehensive Connectivity Map database - identify drugs and gene perturbations with similar or opposite expression signatures - analyze relationships between gene signatures and perturbations across multiple connectivity mapping resources - retrieving metadata of gene sets ## Core Capabilities Perturb-Seqr uses GraphQL API for 3 main query modes. 1. up and down gene set search - search Perturb-Seqr database for gene set pairs that most significantly (Fisher's exact test) mimic or reverse the expression of input gene set 2. single gene set search - returns any significantly overlapping gene sets determined with Fisher's exact test 3. term search - direct metadata search of all dataset gene sets In example Python code below, requests and pandas should be installed for all. **1. Single gene set enrichment** **Parameters:** library_names (optional; specifies which datasets to use) Example code below. ```python import requests import pandas as pd url = "https://perturbseqr.maayanlab.cloud/graphql" def enrich_perturbseqr_single_set(geneset: list, first=1000, library_names=None): variables = { "filterTerm": "", "offset": 0, "first": first, "filterFda": False, "sortBy": "pvalue_up", "filterKo": False, "genes": geneset, } if library_names is not None: variables["libraryNames"] = library_names query = { "operationName": "EnrichmentQuery", "variables": variables, "query": """query EnrichmentQuery( $genes: [String]! $filterTerm: String = "" $offset: Int = 0 $first: Int = 10 $filterFda: Boolean = false $sortBy: String = "" $filterKo: Boolean = false $libraryNames: [String] ) { currentBackground { enrich( genes: $genes filterTerm: $filterTerm offset: $offset first: $first filterFda: $filterFda sortby: $sortBy filterKo: $filterKo libraryNames: $libraryNames ) { nodes { geneSetHash pvalue adjPvalue oddsRatio nOverlap geneSets { nodes { term id nGeneIds geneSetFdaCountsById { nodes { approved count } } library { name } } totalCount } } totalCount geneSetCount consensusCount consensus { drug oddsRatio pvalue adjPvalue approved countSignificant countInsignificant countUpSignificant pvalueUp adjPvalueUp oddsRatioUp pvalueDown adjPvalueDown oddsRatioDown libraries } } } } """, } response = requests.post(url, json=query) response.raise_for_status() res = response.json() consensus = res['data']['currentBackground']['enrich']['consensus'] enrichment = res['data']['currentBackground']['enrich']['nodes'] df_consensus = pd.DataFrame(consensus).rename(columns={'drug': 'perturbation'}) df_enrichment = pd.json_normalize( enrichment, record_path=['geneSets', 'nodes'], meta=['geneSetHash', 'pvalue', 'adjPvalue', 'oddsRatio', 'nOverlap'] ) if df_enrichment.empty: return pd.DataFrame(), pd.DataFrame() df_enrichment["approved"] = df_enrichment["geneSetFdaCountsById.nodes"].map(lambda x: x[0]['approved'] if len(x) > 0 else False) df_enrichment["count"] = df_enrichment["geneSetFdaCountsById.nodes"].map(lambda x: x[0]['count'] if len(x) > 0 else 0) df_enrichment.drop(columns=['geneSetFdaCountsById.nodes'], inplace=True) df_enrichment['direction'] = df_enrichment["term"].map(lambda t: t.split(' ')[-1]) df_enrichment = df_enrichment.rename(columns={'library.name': 'library'}) return df_enrichment, df_consensus ``` **2. Paired up and down gene set enrichment** Example code below. ```python def enrich_perturbseqr_up_down(genes_up: list[str], genes_down: list[str], first=100): query = { "operationName": "PairEnrichmentQuery", "variables": { "filterTerm": "", "offset": 0, "first": first, "filterFda": False, "sortBy": "pvalue_mimic", "filterKo": False, "topN": 1000, "pvalueLe": 0.05, "genesUp": genes_up, "genesDown": genes_down }, "query": """query PairEnrichmentQuery($genesUp: [String]!, $genesDown: [String]!, $filterTerm: String = \"\", $offset: Int = 0, $first: Int = 10, $filterFda: Boolean = false, $sortBy: String = \"\", $filterKo: Boolean = false, $topN: Int = 10000, $pvalueLe: Float = 0.05) { currentBackground { pairedEnrich( filterTerm: $filterTerm offset: $offset first: $first filterFda: $filterFda sortby: $sortBy filterKo: $filterKo topN: $topN pvalueLe: $pvalueLe genesDown: $genesDown genesUp: $genesUp ) { totalCount consensusCount consensus { drug oddsRatio pvalue adjPvalue approved countSignificant countInsignificant countUpSignificant pvalueUp adjPvalueUp oddsRatioUp pvalueDown adjPvalueDown oddsRatioDown libraries } nodes { adjPvalueMimic adjPvalueReverse mimickerOverlap oddsRatioMimic oddsRatioReverse pvalueMimic pvalueReverse reverserOverlap geneSet { nodes { id nGeneIds term geneSetFdaCountsById { nodes { count approved } } library { name } } } } } } } """ } response = requests.post(url, json=query) response.raise_for_status() res = response.json() consensus = res['data']['currentBackground']['pairedEnrich']['consensus'] enrichment = res['data']['currentBackground']['pairedEnrich']['nodes'] df_consensus_pair = pd.DataFrame(consensus).rename(columns={'drug': 'perturbation', 'pvalueUp': 'pvalueMimic', 'pvalueDown': 'pvalueReverse', 'adjPvalueUp': 'adjPvalueMimic', 'adjPvalueDown': 'adjPvalueReverse', 'oddsRatioUp': 'oddsRatioMimic', 'oddsRatioDown': 'oddsRatioReverse' }) df_enrichment_pair = pd.DataFrame(enrichment) df_enrichment_pair['term'] = df_enrichment_pair['geneSet'].map(lambda t: t['nodes'][0]['term'].split(' ')[0]) df_enrichment_pair['approved'] = df_enrichment_pair['geneSet'].map(lambda t: t['nodes'][0]['geneSetFdaCountsById']['nodes'][0]['approved']) df_enrichment_pair['count'] = df_enrichment_pair['geneSet'].map(lambda t: t['nodes'][0]['geneSetFdaCountsById']['nodes'][0]['count']) df_enrichment_pair['nGeneIdsUp'] = df_enrichment_pair['geneSet'].map(lambda t: t['nodes'][0]['nGeneIds']) df_enrichment_pair['nGeneIdsDown'] = df_enrichment_pair['geneSet'].map(lambda t: t['nodes'][0]['nGeneIds']) df_enrichment_pair['library'] = df_enrichment_pair['geneSet'].map(lambda t: t['nodes'][0]['library']['name']) df_enrichment_pair["geneSetIdUp"] = df_enrichment_pair["geneSet"].map( lambda t: next((node['id'] for node in t['nodes'] if ' up' in node['term']), None) ) df_enrichment_pair["geneSetIdDown"] = df_enrichment_pair["geneSet"].map( lambda t: next((node['id'] for node in t['nodes'] if ' down' in node['term']), None) ) df_enrichment_pair = df_enrichment_pair.set_index('term') df_enrichment_pair = df_enrichment_pair.drop(columns=['geneSet']).reset_index(drop=False) return df_enrichment_pair, df_consensus_pair ``` **3. Retrieving overlapping genes using gene set IDs** Example code below. ```python ## Use this function to get the overlap from a user set of genes and a given Perturb-Seqr gene set (id) ## gene set ids are returned as a part of the enrichment query shown above def get_overlap(genes, id): query = { "operationName": "OverlapQuery", "variables": { "id": id, "genes": genes }, "query": """query OverlapQuery($id: UUID!, $genes: [String]!) {geneSet(id: $id) { overlap(genes: $genes) { nodes { symbol ncbiGeneId description summary } }}}""" } response = requests.post(url, json=query) response.raise_for_status() res = response.json() return [item['symbol'] for item in res['data']['geneSet']['overlap']['nodes']] ``` ```python def get_perturbseqr_up_dn_overlap(genes_up: list[str], genes_down: list[str], id_up: str, id_down: str, overlap_type: str): if overlap_type == 'mimicker': up_up_overlap = get_overlap(genes_up, id_up) dn_dn_overlap = get_overlap(genes_down, id_down) return list(set(up_up_overlap) | set(dn_dn_overlap)) elif overlap_type == 'reverser': up_dn_overlap = get_overlap(genes_up, id_down) dn_up_overlap = get_overlap(genes_down, id_up) return list(set(up_dn_overlap) | set(dn_up_overlap)) ``` **4. Retrieving overlap between input gene set and Perturb-Seqr background** Returns converted symbols. Example code below. ```python def get_perturbseqr_valid_genes(genes: list[str]): query = { "query": """query GenesQuery($genes: [String]!) { geneMap2(genes: $genes) { nodes { gene geneInfo { symbol } } } }""", "variables": {"genes": genes}, "operationName": "GenesQuery" } response = requests.post(url, json=query) response.raise_for_status() res = response.json() return [g['geneInfo']['symbol'] for g in res['data']['geneMap2']['nodes'] if g['geneInfo'] != None] ``` ## Best Practices for Error Handling 1. **Choosing between enrichment and consensus:** please note that consensus indicates consensus perturbations, which will compute significant compounds and gene perturbations using the Fisher's exact test, taking into account the number of significant and insignificant gene sets corresponding to that perturbation. Consensus differs from enrichment in that consensus performs another Fisher's exact test to see if a signature is repeated across many datasets. In other words, consensus focuses on robust perturbations rather than isolated experiments. Enrichment simply finds individual perturbation signatures that significantly overlap with the input gene set. For best practice, offer the user the option to choose between either and default to enrichment. 2. **Parsing between drugs and genes:** please note that enrichment and consensus contain both drugs and genes, even when the term names are labeled otherwise. Enrichment drug/gene names will be under the section "term" and are the first characters in the group of information separated by ":". For consensus, both gene and drug names are listed under "drug." Therefore, to parse between the two, do not rely on ``filter_ko``, and instead refer to the libraries list below to see which datasets belong to drugs and which belongs to genes: **Drug/Chemical Perturbation Datasets:** - LINCS L1000 Chemical Perturbations (L2S2) - Tahoe-100M (Hugging Face) - NIBR DRUG-seq U2OS MoABox Dataset (Harmonizome) - Microarrays CMAP Build 02 (Harmonizome) - Ginkgo Bioworks (Hugging Face) - SciPlex Drug Perturbations (Harmonizome) - DeepCover MoA (Harmonizome) - CREEDS Drug Perturbations (CREEDS) - RummaGEO Drug Perturbations (Harmonizome) **Gene Perturbation Datasets:** - Perturb Atlas Human Gene Perturbations (Harmonizome) - Perturb Atlas Mouse Gene Perturbations (Harmonizome) - LINCS L1000 Gene Knockouts (L2S2) - Replogle Perturb-seq K562 Genome-wide (Harmonizome) - Replogle Perturb-seq RPE1 Essential (Harmonizome) - Replogle Perturb-seq K562 Essential (Harmonizome) - CREEDS Gene Perturbations (CREEDS) - Bridge2AI CM4AI (Perturb-seq) - RummaGEO Gene Perturbations (Harmonizome) 3. **Parsing term names for enrichment:** Perturb-Seqr's gene-set "term" strings pack the perturbation name together with metadata (cell line, tissue, dose, KO/KD/Mut label, dataset ID, gene set size, etc.), all separated by "::", e.g.: - "glucagon::::::::::human::GSE171352::2::3 down" - "TP53::CAL51::::::::KO::::SRP287851::2996 down" - "doxorubicin::::human dermal fibroblast::24h::0.3uM" Note that the perturbation name is always the **first** "::"-delimited field and everything after is metadata. The direction ("up"/"down"), when present, is the last whitespace-separated token of the **last** "::" field. Not every term has a direction. For how to best parse term names, refer to parsing function below: ```python def parse_perturbation_term(term): """ Examples: >>> parse_perturbation_term("TP53::CAL51::::::::KO::::SRP287851::2996 down") ('TP53', 'down') >>> parse_perturbation_term("glucagon::::::::::human::GSE171352::2::3 down") ('glucagon', 'down') >>> parse_perturbation_term("doxorubicin::::human dermal fibroblast::24h::0.3uM") ('doxorubicin', None) >>> parse_perturbation_term("TP53::SUP-T1::::::::KO::::SRP999::111 up") ('TP53', 'up') """ fields = term.split("::") perturbation = fields[0].strip() last_field = fields[-1].strip() if fields else "" tokens = last_field.split(" ") direction = tokens[-1].lower() if len(tokens) > 1 and tokens[-1].lower() in ("up", "down") else None return perturbation, direction ``` 3. **Raise TopN to 10000 to include smaller libraries.** Smaller libraries, like SciPlex, may have mimickers/reversers which are statistically significant, but not counted due to being outcrowded by a bigger library with more extreme p-values. Perturb-Seqr paired enrichment query has no per-library filter, so to make sure smaller libraries have enough room to surface, ensure TopN is a large enough value. Setting TopN = 10000 is generally a good value to make sure smaller libraries have enough room to surface. If the server's total count exceeds what is returned, TopN may need to be raised even higher. Including statistically significant results for smaller libraries is useful when graphing results by libraries. ## Additional Resources - Perturb-Seqr Documentation: https://perturbseqr.maayanlab.cloud/help - Requests library: https://pypi.org/project/requests/ - Pandas library: https://pandas.pydata.org/docs/user_guide/index.html#user-guide ## Citations John K. Gardner, Lily D. Taub, Daniel J. B. Clarke, Ido Diamant, Avi Ma'ayan. Perturb-Seqr: Comprehensive Signature Search Engine Integrating Connectivity Maps from Multiple Sources. https://perturbseqr.maayanlab.cloud/