Network Construction and Feature Selection with BoBa-T
This tutorial demonstrates how to:
Build gene regulatory networks from scATAC-seq data
Perform LASSO-based feature selection to prune networks
Adaptively select regularization parameters for optimal network topology
Compare and visualize network properties
Setup and Imports
[1]:
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import networkx as nx
import bobaT as bb
# Set plotting style
sns.set_style('whitegrid')
plt.rcParams['figure.dpi'] = 100
# Suppress warnings for cleaner output
import warnings
warnings.filterwarnings('ignore')
/Users/smgroves/opt/anaconda3/lib/python3.8/site-packages/graph_tool/all.py:40: RuntimeWarning: Error importing draw module, proceeding nevertheless: dlopen(/Users/smgroves/opt/anaconda3/lib/python3.8/site-packages/graph_tool/draw/libgraph_tool_draw.so, 0x0009): Library not loaded: @rpath/libsigc-2.0.0.dylib
Referenced from: <455D357A-072E-3747-8F45-BCDCA894EEE5> /Users/smgroves/opt/anaconda3/lib/libcairomm-1.0.1.dylib
Reason: tried: '/Users/smgroves/opt/anaconda3/lib/libsigc-2.0.0.dylib' (no such file), '/System/Volumes/Preboot/Cryptexes/OS/Users/smgroves/opt/anaconda3/lib/libsigc-2.0.0.dylib' (no such file), '/Users/smgroves/opt/anaconda3/lib/libsigc-2.0.0.dylib' (no such file), '/Users/smgroves/opt/anaconda3/lib/libsigc-2.0.0.dylib' (no such file), '/System/Volumes/Preboot/Cryptexes/OS/Users/smgroves/opt/anaconda3/lib/libsigc-2.0.0.dylib' (no such file), '/Users/smgroves/opt/anaconda3/lib/python3.8/site-packages/graph_tool/draw/../../../../libsigc-2.0.0.dylib' (no such file), '/Users/smgroves/opt/anaconda3/lib/libsigc-2.0.0.dylib' (no such file), '/System/Volumes/Preboot/Cryptexes/OS/Users/smgroves/opt/anaconda3/lib/libsigc-2.0.0.dylib' (no such file), '/Users/smgroves/opt/anaconda3/lib/python3.8/site-packages/graph_tool/draw/../../../../libsigc-2.0.0.dylib' (no such file), '/Users/smgroves/opt/anaconda3/lib/libsigc-2.0.0.dylib' (no such file), '/System/Volumes/Preboot/Cryptexes/OS/Users/smgroves/opt/anaconda3/lib/libsigc-2.0.0.dylib' (no such file), '/Users/smgroves/opt/anaconda3/bin/../lib/libsigc-2.0.0.dylib' (no such file), '/Users/smgroves/opt/anaconda3/lib/libsigc-2.0.0.dylib' (no such file), '/System/Volumes/Preboot/Cryptexes/OS/Users/smgroves/opt/anaconda3/lib/libsigc-2.0.0.dylib' (no such file), '/Users/smgroves/opt/anaconda3/bin/../lib/libsigc-2.0.0.dylib' (no such file), '/usr/local/lib/libsigc-2.0.0.dylib' (no such file), '/usr/lib/libsigc-2.0.0.dylib' (no such file, not in dyld cache)
warnings.warn(msg, RuntimeWarning)
Configure Paths
[2]:
# Set your working directory
dir_prefix = '.'
# Input directories
DIRECT_NET_DIR = os.path.join(dir_prefix, 'test_data')
DATA_DIR = os.path.join(dir_prefix, 'test_data')
# Output directory
NETWORK_OUTDIR = os.path.join(dir_prefix, 'networks')
os.makedirs(NETWORK_OUTDIR, exist_ok=True)
print(f"Working directory: {dir_prefix}")
print(f"DIRECT-NET data: {DIRECT_NET_DIR}")
print(f"Output directory: {NETWORK_OUTDIR}")
Working directory: .
DIRECT-NET data: ./test_data
Output directory: ./networks
Part 1: Load and Explore DIRECT-NET Data
DIRECT-NET provides TF-target predictions with motif scores. We’ll explore the distribution of scores and set an initial threshold. BoBa-T works with any data file that has TFs, target nodes, and optimally some sort of score (like a motif score).
[3]:
# Load DIRECT-NET predictions
direct_net_file = os.path.join(DIRECT_NET_DIR, 'Direct_net_0.1.csv')
direct_net = pd.read_csv(direct_net_file, header=0, index_col=0)
# Standardize gene names to uppercase
direct_net['Target_gene'] = direct_net['Target_gene'].str.upper()
print(f"DIRECT-NET data loaded: {direct_net.shape[0]} TF-target interactions")
print(f"Columns: {list(direct_net.columns)}")
direct_net.head()
DIRECT-NET data loaded: 1298 TF-target interactions
Columns: ['regulatory regions', 'TF motif', 'motif score', 'Target_gene']
[3]:
| regulatory regions | TF motif | motif score | Target_gene | |
|---|---|---|---|---|
| 1 | chr5_136394557_136395437_Cux1 | PPARG | 9.432222 | CUX1 |
| 2 | chr6_72736650_72737562_Tcf7l1 | PPARG | 13.168407 | TCF7L1 |
| 3 | chr6_72740374_72741224_Tcf7l1 | PPARG | 9.306865 | TCF7L1 |
| 4 | chr9_63968292_63969164_Smad3 | PPARG | 9.969222 | SMAD3 |
| 5 | chr16_50594630_50595552_Bbx | PPARG | 12.514526 | BBX |
Visualize Motif Score Distribution
[4]:
plt.figure(figsize=(10, 5))
plt.hist(direct_net['motif score'], bins=50, edgecolor='black', alpha=0.7)
plt.xlabel('Motif Score')
plt.ylabel('Frequency')
plt.title('Distribution of DIRECT-NET Motif Scores')
plt.axvline(0, color='red', linestyle='--', label='Score = 0')
plt.legend()
plt.tight_layout()
plt.show()
print(f"\nMotif score statistics:")
print(direct_net['motif score'].describe())
Motif score statistics:
count 1298.000000
mean 12.524504
std 2.138108
min 3.607144
25% 11.659460
50% 12.551152
75% 13.576800
max 21.885676
Name: motif score, dtype: float64
Part 2: Build Initial Network
We’ll construct a directed graph from the DIRECT-NET data, applying a threshold to filter weak interactions.
[5]:
# Set threshold for motif score
# You can adjust this based on the distribution above
THRESHOLD = 0 # Change this value based on your data
print(f"Building network with motif score threshold: {THRESHOLD}")
# Build network
G = bb.net.build_network_from_directnet(
direct_net,
threshold=THRESHOLD,
threshold_var='motif score'
)
print(f"\nNetwork statistics:")
print(f" Nodes: {G.number_of_nodes()}")
print(f" Edges: {G.number_of_edges()}")
print(f" Average in-degree: {np.mean([d for n, d in G.in_degree()]):.2f}")
print(f" Average out-degree: {np.mean([d for n, d in G.out_degree()]):.2f}")
Building network with motif score threshold: 0
Network built with 72 nodes and 1004 edges
Network statistics:
Nodes: 72
Edges: 1004
Average in-degree: 13.94
Average out-degree: 13.94
Analyze Network Topology
[6]:
# Calculate degree distributions
in_degrees = [d for n, d in G.in_degree()]
out_degrees = [d for n, d in G.out_degree()]
# Plot degree distributions
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# In-degree (number of regulators per gene)
axes[0].hist(in_degrees, bins=range(max(in_degrees) + 2),
edgecolor='black', alpha=0.7)
axes[0].set_xlabel('Number of Regulators (In-degree)')
axes[0].set_ylabel('Number of Genes')
axes[0].set_title('In-Degree Distribution')
axes[0].axvline(np.mean(in_degrees), color='red',
linestyle='--', label=f'Mean = {np.mean(in_degrees):.1f}')
axes[0].legend()
# Out-degree (number of targets per TF)
axes[1].hist(out_degrees, bins=range(max(out_degrees) + 2),
edgecolor='black', alpha=0.7)
axes[1].set_xlabel('Number of Targets (Out-degree)')
axes[1].set_ylabel('Number of TFs')
axes[1].set_title('Out-Degree Distribution')
axes[1].axvline(np.mean(out_degrees), color='red',
linestyle='--', label=f'Mean = {np.mean(out_degrees):.1f}')
axes[1].legend()
plt.tight_layout()
plt.show()
# Identify genes with many regulators
high_indegree = [(n, d) for n, d in G.in_degree() if d > 10]
if high_indegree:
print(f"\nGenes with >10 regulators: {len(high_indegree)}")
print("Top 5:")
for gene, degree in sorted(high_indegree, key=lambda x: x[1], reverse=True)[:5]:
print(f" {gene}: {degree} regulators")
Genes with >10 regulators: 44
Top 5:
NFE2L2: 34 regulators
SMAD3: 33 regulators
CUX2: 31 regulators
EHF: 30 regulators
FOXO3: 30 regulators
Save Initial Network
[7]:
# Save network to CSV
network_name = "DIRECT-NET_network_2020db_0.1"
network_file = os.path.join(NETWORK_OUTDIR, f"{network_name}.csv")
bb.net.save_network(network_file, G, attributes=True, overwrite=True)
# Also save as DataFrame for easier manipulation
network_df = pd.DataFrame([
{
'source': edge[0],
'target': edge[1],
'weight': G[edge[0]][edge[1]].get('weight', np.nan),
'evidence': G[edge[0]][edge[1]].get('evidence', '')
}
for edge in G.edges()
])
print(f"\nNetwork saved as DataFrame:")
print(network_df.head())
Network saved to ./networks/DIRECT-NET_network_2020db_0.1.csv
Network saved as DataFrame:
source target weight evidence
0 NFKB1 PRDM16 5.254761 direct-net
1 NFKB1 NFIX 18.630681 direct-net
2 NFKB1 SMAD3 6.030182 direct-net
3 NFKB1 FOS 8.319178 direct-net
4 NFKB1 NFE2L2 9.749579 direct-net
Part 3: Load Expression Data
We need expression data to perform LASSO regression. Each column should be a gene, each row an observation (cell, sample, etc.).
[14]:
# Load expression data
# Adjust the filename as needed
data_path = os.path.join(DATA_DIR, './demo_dataset.csv')
# Check if file exists
if os.path.exists(data_path):
data = pd.read_csv(data_path, index_col=0, header=0)
# Standardize gene names to uppercase
data.columns = [col.upper() for col in data.columns]
print(f"Expression data loaded: {data.shape[0]} observations, {data.shape[1]} genes")
print(f"\nFirst few genes: {list(data.columns[:5])}")
# Verify overlap with network
network_genes = set(network_df['source']).union(set(network_df['target']))
data_genes = set(data.columns)
overlap = network_genes.intersection(data_genes)
print(f"\nGene overlap:")
print(f" Network genes: {len(network_genes)}")
print(f" Data genes: {len(data_genes)}")
print(f" Overlap: {len(overlap)} ({len(overlap)/len(network_genes)*100:.1f}%)")
if len(overlap) < len(network_genes) * 0.5:
print("\nWARNING: Less than 50% gene overlap! Check gene naming conventions.")
else:
print(f"ERROR: Expression data file not found at {data_path}")
print("Please update the data_path variable with the correct file location.")
data = None
Expression data loaded: 501 observations, 94 genes
First few genes: ['REST', 'ETS1', 'NFE2L2', 'PLAGL1', 'ICAM1']
Gene overlap:
Network genes: 72
Data genes: 94
Overlap: 72 (100.0%)
Part 4: LASSO Feature Selection
Now we’ll use LASSO regression to identify the most important regulators for each target gene. LASSO will:
Perform L1 regularization to drive weak coefficients to zero
Select only the most predictive TFs for each target
Test multiple alpha values to control regularization strength
Note: This step can take several minutes depending on network size.
[15]:
# Only run if expression data is loaded
if data is not None:
# Define alpha values to test
# Smaller alpha = less regularization = more features retained
# Larger alpha = more regularization = fewer features retained
alphas = [0.00001, 0.0001, 0.001, 0.01, 0.1, 1]
print(f"Running LASSO feature selection with alphas: {alphas}")
print("This may take several minutes...\n")
# Run LASSO feature selection
alpha_scores = bb.net.lasso_feature_selection(
network=network_df,
data=data,
network_name=network_name,
output_dir=NETWORK_OUTDIR,
alphas=alphas,
save_network=True,
plot=True
)
print("\nLASSO feature selection complete!")
print(f"Results saved to: {NETWORK_OUTDIR}/feature_selection/{network_name}/")
else:
print("Skipping LASSO: expression data not loaded")
Running LASSO feature selection with alphas: [1e-05, 0.0001, 0.001, 0.01, 0.1, 1]
This may take several minutes...
Running LASSO feature selection for 59 targets...
LASSO feature selection complete!
LASSO feature selection complete!
Results saved to: ./networks/feature_selection/DIRECT-NET_network_2020db_0.1/
Visualize LASSO Performance Across Alpha Values
[16]:
if data is not None and alpha_scores:
# Calculate mean scores for each alpha
mean_scores = {alpha: np.mean(scores) for alpha, scores in alpha_scores.items()}
plt.figure(figsize=(10, 5))
plt.plot(list(mean_scores.keys()), list(mean_scores.values()),
marker='o', linewidth=2, markersize=8)
plt.xscale('log')
plt.xlabel('Alpha (regularization strength)', fontsize=12)
plt.ylabel('Mean R² Score', fontsize=12)
plt.title('LASSO Performance vs Regularization Strength', fontsize=14)
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()
print("\nMean R² scores by alpha:")
for alpha, score in sorted(mean_scores.items()):
print(f" α = {alpha:7.5f}: R² = {score:.4f}")
Mean R² scores by alpha:
α = 0.00001: R² = 0.9097
α = 0.00010: R² = 0.9030
α = 0.00100: R² = 0.8505
α = 0.01000: R² = 0.6722
α = 0.10000: R² = 0.0776
α = 1.00000: R² = -0.0145
Part 5: Compare Networks Across Alpha Values
Let’s compare how different alpha values affect network topology.
[17]:
if data is not None:
comparison_results = []
for alpha in alphas:
# Load filtered network
filtered_net_path = os.path.join(
NETWORK_OUTDIR,
'feature_selection',
network_name,
f'Lasso_alpha_{alpha}',
f'{network_name}_Lasso_{alpha}.csv'
)
if os.path.exists(filtered_net_path):
filtered_net = pd.read_csv(filtered_net_path)
# Compare networks
results = bb.net.compare_networks(
original_network=network_df,
filtered_network=filtered_net,
alpha=alpha,
save=True,
save_dir=os.path.join(
NETWORK_OUTDIR,
'feature_selection',
network_name
)
)
comparison_results.append(results)
else:
print(f"Warning: Filtered network not found for alpha={alpha}")
# Summarize results
if comparison_results:
print("\n" + "="*70)
print("SUMMARY: Network Pruning Across Alpha Values")
print("="*70)
summary_df = pd.DataFrame([
{
'Alpha': r['alpha'],
'Edges Dropped': r['n_dropped'],
'Max Parents': r['max_parents'],
'Target (max parents)': r['max_parents_target']
}
for r in comparison_results
])
print(summary_df.to_string(index=False))
else:
print("Skipping comparison: expression data not loaded")
Network comparison (α=1e-05):
Dropped edges: 71
Target with most parents: NFE2L2 (31 parents)
Parents: ['ASCL1', 'BACH1', 'BACH2', 'CREB1', 'EGR1', 'EHF', 'FOS', 'FOXO3', 'GLIS3', 'GRHL2', 'HSF2', 'JUN', 'JUNB', 'JUND', 'LMX1B', 'MEIS2', 'NFIA', 'NFKB1', 'NFYC', 'NR3C2', 'PKNOX2', 'RBPJ', 'REST', 'RORA', 'SIX1', 'TCF4', 'TCF7L1', 'TCF7L2', 'THRB', 'ZBTB7A', 'ZEB1']
Network comparison (α=0.0001):
Dropped edges: 356
Target with most parents: EHF (19 parents)
Parents: ['BACH1', 'CREB1', 'EHF', 'ESR1', 'FOS', 'GLIS3', 'JUN', 'JUNB', 'JUND', 'NFATC2', 'NFIB', 'NFIX', 'NR6A1', 'RORA', 'RORB', 'STAT1', 'TCF7L2', 'TEAD1', 'THRB']
Network comparison (α=0.001):
Dropped edges: 700
Target with most parents: EHF (10 parents)
Parents: ['EHF', 'ESR1', 'FOS', 'JUN', 'JUND', 'NFATC2', 'NR6A1', 'RORB', 'RUNX1', 'STAT1']
Network comparison (α=0.01):
Dropped edges: 876
Target with most parents: MEIS2 (5 parents)
Parents: ['ASCL1', 'EGR1', 'NFIB', 'REST', 'RUNX1']
Network comparison (α=0.1):
Dropped edges: 975
Target with most parents: HES1 (2 parents)
Parents: ['HES1', 'RUNX1']
Network comparison (α=1):
Dropped edges: 989
Target with most parents: BACH1 (1 parents)
Parents: ['BACH1']
======================================================================
SUMMARY: Network Pruning Across Alpha Values
======================================================================
Alpha Edges Dropped Max Parents Target (max parents)
0.00001 71 31 NFE2L2
0.00010 356 19 EHF
0.00100 700 10 EHF
0.01000 876 5 MEIS2
0.10000 975 2 HES1
1.00000 989 1 BACH1
Part 6: Adaptive LASSO Pruning
Instead of using a single alpha for all targets, we can adaptively select the alpha that reduces each target’s regulators to a desired maximum (e.g., 8 regulators).
[18]:
if data is not None:
# Set maximum number of regulators per target
MAX_PARENTS = 8
print(f"Running adaptive LASSO pruning (max {MAX_PARENTS} parents per target)...\n")
# Perform adaptive pruning
pruned_network, tf_alpha_dict = bb.net.adaptive_lasso_pruning(
network=network_df,
output_dir=NETWORK_OUTDIR,
network_name=network_name,
max_parents=MAX_PARENTS,
alphas=alphas
)
print(f"\nAdaptive pruning complete!")
print(f"Original network: {len(network_df)} edges")
print(f"Pruned network: {len(pruned_network)} edges")
print(f"Reduction: {(1 - len(pruned_network)/len(network_df))*100:.1f}%")
# Show distribution of selected alphas
alpha_counts = pd.Series(list(tf_alpha_dict.values())).value_counts().sort_index()
print(f"\nAlpha distribution across targets:")
print(alpha_counts)
# Save pruned network
pruned_network_file = os.path.join(
NETWORK_OUTDIR,
'feature_selection',
network_name,
f'combined_{network_name}_Lasso_adaptive_max{MAX_PARENTS}.csv'
)
pruned_network.to_csv(pruned_network_file, index=False, header=False)
print(f"\nPruned network saved to: {pruned_network_file}")
else:
print("Skipping adaptive pruning: expression data not loaded")
Running adaptive LASSO pruning (max 8 parents per target)...
Adaptively pruning network to max 8 parents per target...
Adaptive pruning complete! Final network: 322 edges
Adaptive pruning complete!
Original network: 1004 edges
Pruned network: 322 edges
Reduction: 67.9%
Alpha distribution across targets:
0.00000 24
0.00001 1
0.00010 4
0.00100 39
0.01000 4
dtype: int64
Pruned network saved to: ./networks/feature_selection/DIRECT-NET_network_2020db_0.1/combined_DIRECT-NET_network_2020db_0.1_Lasso_adaptive_max8.csv
Visualize Alpha Selection
[19]:
if data is not None and 'tf_alpha_dict' in locals():
# Plot alpha distribution
alpha_values = [v for v in tf_alpha_dict.values() if v != 'NA' and v != 0]
if alpha_values:
plt.figure(figsize=(10, 5))
plt.hist(alpha_values, bins=len(set(alpha_values)),
edgecolor='black', alpha=0.7)
plt.xlabel('Selected Alpha Value', fontsize=12)
plt.ylabel('Number of Targets', fontsize=12)
plt.title('Distribution of Selected Alpha Values (Adaptive Pruning)',
fontsize=14)
plt.xscale('log')
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()
# Show some example targets
print("\nExample targets and their selected alphas:")
for i, (target, alpha) in enumerate(list(tf_alpha_dict.items())[:10]):
n_parents = len(pruned_network[pruned_network['target'] == target])
print(f" {target}: α = {alpha}, {n_parents} regulators")
Example targets and their selected alphas:
AHR: α = 0.001, 7 regulators
ASCL1: α = 0.001, 6 regulators
BACH1: α = 0.001, 6 regulators
BACH2: α = 0.001, 4 regulators
BBX: α = 0.001, 3 regulators
CREB1: α = 0, 0 regulators
CUX1: α = 0.0001, 6 regulators
CUX2: α = 0.001, 6 regulators
EGR1: α = 0, 8 regulators
EHF: α = 0.01, 2 regulators
Part 7: Visualize Network Changes
Let’s visualize a subnetwork to see the effects of pruning.
[22]:
if data is not None and 'pruned_network' in locals():
# Pick an interesting target gene (one that had many regulators)
target_of_interest = "ASCL1"
if target_of_interest:
# Get original and pruned edges
original_edges = network_df[network_df['target'] == target_of_interest]
pruned_edges = pruned_network[pruned_network['target'] == target_of_interest]
print(f"Example: {target_of_interest}")
print(f" Original regulators: {len(original_edges)}")
print(f" Pruned regulators: {len(pruned_edges)}")
print(f" Kept: {list(pruned_edges['source'].values)}")
# Create subgraphs
fig, axes = plt.subplots(1, 2, figsize=(16, 6))
# Original network
G_orig = nx.DiGraph()
G_orig.add_node(target_of_interest)
for _, row in original_edges.iterrows():
G_orig.add_edge(row['source'], row['target'])
pos_orig = nx.spring_layout(G_orig, k=0.5, seed=42)
nx.draw_networkx(
G_orig, pos_orig, ax=axes[0],
node_color=['#ff7f0e' if n == target_of_interest else '#1f77b4'
for n in G_orig.nodes()],
node_size=1000,
font_size=8,
arrows=True,
arrowsize=10,
edge_color='gray',
alpha=0.7
)
axes[0].set_title(f'Original Network: {target_of_interest}\n'
f'{len(original_edges)} regulators',
fontsize=12)
axes[0].axis('off')
# Pruned network
G_pruned = nx.DiGraph()
G_pruned.add_node(target_of_interest)
for _, row in pruned_edges.iterrows():
G_pruned.add_edge(row['source'], row['target'])
pos_pruned = nx.spring_layout(G_pruned, k=0.5, seed=42)
nx.draw_networkx(
G_pruned, pos_pruned, ax=axes[1],
node_color=['#ff7f0e' if n == target_of_interest else '#2ca02c'
for n in G_pruned.nodes()],
node_size=1000,
font_size=8,
arrows=True,
arrowsize=10,
edge_color='gray',
alpha=0.7
)
axes[1].set_title(f'Pruned Network: {target_of_interest}\n'
f'{len(pruned_edges)} regulators (LASSO)',
fontsize=12)
axes[1].axis('off')
plt.tight_layout()
plt.show()
Example: ASCL1
Original regulators: 13
Pruned regulators: 6
Kept: ['EGR1', 'HES1', 'JUN', 'JUND', 'SIX1', 'STAT1']
Bonus: Export for Boolean Modeling
The pruned network can now be used as input for network inference.
[21]:
if data is not None and 'pruned_network' in locals():
# Create a NetworkX graph from pruned network
G_pruned_final = nx.DiGraph()
for _, row in pruned_network.iterrows():
G_pruned_final.add_edge(row['source'], row['target'])
print("Pruned network summary for Boolean modeling:")
print(f" Nodes: {G_pruned_final.number_of_nodes()}")
print(f" Edges: {G_pruned_final.number_of_edges()}")
print(f" Avg in-degree: {np.mean([d for n, d in G_pruned_final.in_degree()]):.2f}")
print(f" Max in-degree: {max([d for n, d in G_pruned_final.in_degree()])}")
# Identify source nodes (potential initial conditions)
source_nodes = [n for n, d in G_pruned_final.in_degree() if d == 0]
print(f"\nSource nodes (no regulators): {len(source_nodes)}")
if source_nodes:
print(f" Examples: {source_nodes[:5]}")
# Ready for BooleaBayes!
print("\n✅ Network is ready for Boolean modeling with BoBa-T!")
Pruned network summary for Boolean modeling:
Nodes: 71
Edges: 322
Avg in-degree: 4.54
Max in-degree: 8
Source nodes (no regulators): 12
Examples: ['CREB1', 'STAT1', 'JUND', 'TCF4', 'RORA']
✅ Network is ready for Boolean modeling with BoBa-T!
[ ]: