{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Identify the driver regulators of cell fate decisions\n",
"CEFCON is a computational tool for deciphering driver regulators of cell fate decisions from single-cell RNA-seq data. It takes a prior gene interaction network and expression profiles from scRNA-seq data associated with a given developmental trajectory as inputs, and consists of three main components, including cell-lineage-specific gene regulatory network (GRN) construction, driver regulator identification and regulon-like gene module (RGM) identification.\n",
"\n",
"Check out [(Wang et al., Nature Communications, 2023)](https://www.nature.com/articles/s41467-023-44103-3) for the detailed methods and applications.\n",
"\n",
"Code: [https://github.com/WPZgithub/CEFCON](https://github.com/WPZgithub/CEFCON)\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
" ____ _ _ __ \n",
" / __ \\____ ___ (_)___| | / /__ _____________ \n",
" / / / / __ `__ \\/ / ___/ | / / _ \\/ ___/ ___/ _ \\ \n",
"/ /_/ / / / / / / / /__ | |/ / __/ / (__ ) __/ \n",
"\\____/_/ /_/ /_/_/\\___/ |___/\\___/_/ /____/\\___/ \n",
"\n",
"Version: 1.5.6, Tutorials: https://omicverse.readthedocs.io/\n"
]
}
],
"source": [
"import omicverse as ov\n",
"#print(f\"omicverse version: {ov.__version__}\")\n",
"import scanpy as sc\n",
"#print(f\"scanpy version: {sc.__version__}\")\n",
"import pandas as pd\n",
"from tqdm.auto import tqdm\n",
"ov.plot_set()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Data loading and processing\n",
"Here, we use the mouse hematopoiesis data provided by [Nestorowa et al. (2016, Blood).](https://doi.org/10.1182/blood-2016-05-716480)\n",
"\n",
"**The scRNA-seq data requires processing to extract lineage information for the CEFCON analysis.** Please refer to the [original notebook](https://github.com/WPZgithub/CEFCON/blob/e74d2d248b88fb3349023d1a97d3cc8a52cc4060/notebooks/preprocessing_nestorowa16_data.ipynb) for detailed instructions on preprocessing scRNA-seq data."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Load mouse_hsc_nestorowa16_v0.h5ad\n"
]
},
{
"data": {
"text/plain": [
"AnnData object with n_obs × n_vars = 1645 × 3000\n",
" obs: 'E_pseudotime', 'GM_pseudotime', 'L_pseudotime', 'label_info', 'n_genes', 'leiden', 'cell_type_roughly', 'cell_type_finely'\n",
" var: 'n_cells', 'highly_variable', 'means', 'dispersions', 'dispersions_norm', 'E_pseudotime_logFC', 'GM_pseudotime_logFC', 'L_pseudotime_logFC'\n",
" uns: 'cell_type_finely_colors', 'cell_type_roughly_colors', 'draw_graph', 'hvg', 'leiden', 'leiden_colors', 'lineages', 'neighbors', 'pca', 'tsne', 'umap'\n",
" obsm: 'X_draw_graph_fa', 'X_pca'\n",
" varm: 'PCs'\n",
" layers: 'raw_count'\n",
" obsp: 'connectivities', 'distances'"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"adata = ov.single.mouse_hsc_nestorowa16()\n",
"adata"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"CEFCON fully exploit an available global and **context-free gene interaction network** as prior knowledge, from which we extract the cell-lineage-specific gene interactions according to the gene expression profiles derived from scRNA-seq data associated with a given developmental trajectory. \n",
"\n",
"You can download the prior network in the [zenodo](https://zenodo.org/records/8013900). **CEFCON only provides the prior network for human and mosue data anaylsis**. For other species, you should provide the prior network mannully.\n",
"\n",
"The author of CEFCON has provided several prior networks here; however, 'nichenet' yields the best results."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Load the prior gene interaction network: nichenet. #Genes: 25345, #Edges: 5290993\n"
]
}
],
"source": [
"prior_network = ov.single.load_human_prior_interaction_network(dataset='nichenet') "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**In the scRNA-seq analysis of human data, you should not run this step. Running it may change the gene symbol and result in errors.**\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Convert genes of the prior interaction network to mouse gene symbols:\n"
]
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "8a71aaf3892d4640872a7e28d51d604e",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"Processing: 0%| | 0/10 [00:00\n",
"\n",
"\n",
" \n",
" \n",
" | \n",
" from | \n",
" to | \n",
"
\n",
" \n",
" \n",
" \n",
" | 0 | \n",
" Klf2 | \n",
" Dlgap1 | \n",
"
\n",
" \n",
" | 2 | \n",
" Klf2 | \n",
" Bhlhe40 | \n",
"
\n",
" \n",
" | 3 | \n",
" Klf2 | \n",
" Rps6ka1 | \n",
"
\n",
" \n",
" | 4 | \n",
" Klf2 | \n",
" Pxn | \n",
"
\n",
" \n",
" | 5 | \n",
" Klf2 | \n",
" Ube2v1 | \n",
"
\n",
" \n",
" | ... | \n",
" ... | \n",
" ... | \n",
"
\n",
" \n",
" | 837982 | \n",
" Zranb1 | \n",
" Zfp141 | \n",
"
\n",
" \n",
" | 837983 | \n",
" Zranb1 | \n",
" Zfy1 | \n",
"
\n",
" \n",
" | 837984 | \n",
" Zranb1 | \n",
" Zfy2 | \n",
"
\n",
" \n",
" | 837987 | \n",
" Zscan21 | \n",
" Zfy1 | \n",
"
\n",
" \n",
" | 837988 | \n",
" Zscan21 | \n",
" Zfy2 | \n",
"
\n",
" \n",
"
\n",
"5029532 rows × 2 columns
\n",
""
],
"text/plain": [
" from to\n",
"0 Klf2 Dlgap1\n",
"2 Klf2 Bhlhe40\n",
"3 Klf2 Rps6ka1\n",
"4 Klf2 Pxn\n",
"5 Klf2 Ube2v1\n",
"... ... ...\n",
"837982 Zranb1 Zfp141\n",
"837983 Zranb1 Zfy1\n",
"837984 Zranb1 Zfy2\n",
"837987 Zscan21 Zfy1\n",
"837988 Zscan21 Zfy2\n",
"\n",
"[5029532 rows x 2 columns]"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Convert the gene symbols of the prior gene interaction network to the mouse gene symbols\n",
"prior_network = ov.single.convert_human_to_mouse_network(prior_network,server_name='asia') \n",
"prior_network"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"prior_network.to_csv('result/combined_network_Mouse.txt.gz',sep='\\t')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Alternatively, you can directly specify the file path of the input prior interaction network and import the specified file."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"#prior_network = './Reference_Networks/combined_network_Mouse.txt'\n",
"prior_network=ov.read('result/combined_network_Mouse.txt.gz',index_col=0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Training CEFCON model\n",
"\n",
"We recommend using GRUOBI to solve the integer linear programming (ILP) problem when identifying driver genes. GUROBI is a commercial solver that requires licenses to run. Thankfully, it provides free licenses in academia, as well as trial licenses outside academia. If there is no problem about the licenses, you need to install the `gurobipy` package.\n",
"\n",
"If difficulties arise while using GUROBI, the non-commercial solver, SCIP, will be employed as an alternative. But the use of SCIP does not come with a guarantee of achieving a successful solutio\n",
"\n",
"**By default, the program will verify the availability of GRUOBI. If GRUOBI is not accessible, it will automatically switch the solver to SCIP.**\n"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"tags": []
},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"CEFCON_obj = ov.single.pyCEFCON(adata, prior_network, repeats=5, solver='GUROBI')\n",
"CEFCON_obj"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Construct cell-lineage-specific GRNs"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Start data preparation\n",
"\n",
"[0] - Data loading and preprocessing...\n",
"Consider the input data with 3 lineages:\n",
" Lineage - E_pseudotime:\n",
" 2935 extra edges (Spearman correlation > 0.6) are added into the prior gene interaction network.\n",
" Total number of edges: 90036.\n",
" n_genes × n_cells = 2803 × 1065\n",
" Lineage - GM_pseudotime:\n",
" 100 extra edges (Spearman correlation > 0.6) are added into the prior gene interaction network.\n",
" Total number of edges: 90136.\n",
" n_genes × n_cells = 2803 × 882\n",
" Lineage - L_pseudotime:\n",
" 4 extra edges (Spearman correlation > 0.6) are added into the prior gene interaction network.\n",
" Total number of edges: 90140.\n",
" n_genes × n_cells = 2803 × 843\n"
]
}
],
"source": [
"CEFCON_obj.preprocess()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Lineage-by-lineage computation:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Start model training\n",
"\n",
"[1] - Constructing cell-lineage-specific GRN...\n",
" Lineage - E_pseudotime: \n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
" Iter: 1/5: 100%|██████████████| 350/350 [00:39<00:00, 8.89it/s, loss=0.000174, min_loss=0.000107]\n",
" Iter: 2/5: 100%|██████████████| 350/350 [00:39<00:00, 8.86it/s, loss=0.000228, min_loss=0.000144]\n",
" Iter: 3/5: 100%|██████████████| 350/350 [00:40<00:00, 8.70it/s, loss=0.000285, min_loss=0.000177]\n",
" Iter: 4/5: 100%|██████████████| 350/350 [00:40<00:00, 8.56it/s, loss=0.000265, min_loss=0.000151]\n",
" Iter: 5/5: 100%|██████████████| 350/350 [00:40<00:00, 8.54it/s, loss=0.000566, min_loss=0.000105]\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"[1] - Constructing cell-lineage-specific GRN...\n",
" Lineage - GM_pseudotime: \n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
" Iter: 1/5: 100%|███████████████| 350/350 [00:40<00:00, 8.56it/s, loss=0.00605, min_loss=0.000518]\n",
" Iter: 2/5: 100%|██████████████| 350/350 [00:40<00:00, 8.54it/s, loss=0.000623, min_loss=0.000159]\n",
" Iter: 3/5: 100%|██████████████| 350/350 [00:40<00:00, 8.56it/s, loss=0.000696, min_loss=0.000167]\n",
" Iter: 4/5: 100%|██████████████| 350/350 [00:40<00:00, 8.56it/s, loss=0.000175, min_loss=0.000149]\n",
" Iter: 5/5: 100%|██████████████| 350/350 [00:40<00:00, 8.57it/s, loss=0.000212, min_loss=0.000212]\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"[1] - Constructing cell-lineage-specific GRN...\n",
" Lineage - L_pseudotime: \n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
" Iter: 1/5: 100%|███████████████| 350/350 [00:40<00:00, 8.56it/s, loss=0.00046, min_loss=0.000136]\n",
" Iter: 2/5: 100%|██████████████| 350/350 [00:40<00:00, 8.56it/s, loss=0.000271, min_loss=0.000128]\n",
" Iter: 3/5: 100%|██████████████| 350/350 [00:40<00:00, 8.56it/s, loss=0.000186, min_loss=0.000186]\n",
" Iter: 4/5: 100%|██████████████| 350/350 [00:40<00:00, 8.57it/s, loss=0.000251, min_loss=0.000169]\n",
" Iter: 5/5: 100%|██████████████| 350/350 [00:40<00:00, 8.55it/s, loss=0.000875, min_loss=0.000139]\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Finish model training\n",
"\n"
]
}
],
"source": [
"CEFCON_obj.train()"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Start predict lineage - E_pseudotime:\n",
"Start calculate gene influence score - E_pseudotime:\n",
"Start calculate gene driver regulators - E_pseudotime:\n",
"[2] - Identifying driver regulators...\n",
" Solving MFVS problem...\n",
" 176 critical nodes are found.\n",
" 0 nodes left after graph reduction operation.\n",
" 176 MFVS driver genes are found.\n",
" Solving MDS problem...\n",
" 15 critical nodes are found.\n",
" 1124 nodes left after graph reduction operation.\n",
" Solving the Integer Linear Programming problem on the reduced graph...\n",
" Solving by GUROBI...(optimal value with GUROBI:116.0, status:optimal)\n",
" 131 MDS driver genes are found.\n",
"Start predict lineage - GM_pseudotime:\n",
"Start calculate gene influence score - GM_pseudotime:\n",
"Start calculate gene driver regulators - GM_pseudotime:\n",
"[2] - Identifying driver regulators...\n",
" Solving MFVS problem...\n",
" 232 critical nodes are found.\n",
" 0 nodes left after graph reduction operation.\n",
" 232 MFVS driver genes are found.\n",
" Solving MDS problem...\n",
" 18 critical nodes are found.\n",
" 1572 nodes left after graph reduction operation.\n",
" Solving the Integer Linear Programming problem on the reduced graph...\n",
" Solving by GUROBI...(optimal value with GUROBI:206.0, status:optimal)\n",
" 224 MDS driver genes are found.\n",
"Start predict lineage - L_pseudotime:\n",
"Start calculate gene influence score - L_pseudotime:\n",
"Start calculate gene driver regulators - L_pseudotime:\n",
"[2] - Identifying driver regulators...\n",
" Solving MFVS problem...\n",
" 180 critical nodes are found.\n",
" 4 nodes left after graph reduction operation.\n",
" Solving the Integer Linear Programming problem on the reduced graph...\n",
" Solving by GUROBI...(optimal value with GUROBI:2.0, status:optimal)\n",
" 182 MFVS driver nodes are found.\n",
" Solving MDS problem...\n",
" 18 critical nodes are found.\n",
" 1747 nodes left after graph reduction operation.\n",
" Solving the Integer Linear Programming problem on the reduced graph...\n",
" Solving by GUROBI...(optimal value with GUROBI:300.0, status:optimal)\n",
" 318 MDS driver genes are found.\n"
]
}
],
"source": [
"# Idenytify driver regulators for each lineage\n",
"CEFCON_obj.predicted_driver_regulators()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can find out the driver regulators identified by CEFCON."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"scrolled": true,
"tags": []
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" influence_score | \n",
" is_driver_regulator | \n",
" is_MFVS_driver | \n",
" is_MDS_driver | \n",
" is_TF | \n",
"
\n",
" \n",
" \n",
" \n",
" | JUN | \n",
" 7.352254 | \n",
" True | \n",
" True | \n",
" True | \n",
" True | \n",
"
\n",
" \n",
" | GATA1 | \n",
" 7.071392 | \n",
" True | \n",
" True | \n",
" True | \n",
" True | \n",
"
\n",
" \n",
" | FOS | \n",
" 6.930125 | \n",
" True | \n",
" True | \n",
" True | \n",
" True | \n",
"
\n",
" \n",
" | GATA2 | \n",
" 6.683559 | \n",
" True | \n",
" False | \n",
" True | \n",
" True | \n",
"
\n",
" \n",
" | MEIS1 | \n",
" 5.851068 | \n",
" True | \n",
" True | \n",
" True | \n",
" True | \n",
"
\n",
" \n",
"
\n",
"