Skip to content

Commit 2cac736

Browse files
janursajanursa
committed
raw scores for sem and vc
1 parent c57cc7f commit 2cac736

9 files changed

Lines changed: 118 additions & 210 deletions

File tree

‎scripts/configs/labels_tw.config‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -46,11 +46,11 @@ process {
4646

4747

4848
// Resource labels
49-
withLabel: {lowtime: 1.h}
50-
withLabel: {midtime: 4.h}
51-
withLabel: {hightime: 8.h}
52-
withLabel: {onedaytime: 24.h}
53-
withLabel: {twodaytime: 48.h}
49+
withLabel: lowtime {time = 1.h}
50+
withLabel: midtime {time = 4.h}
51+
withLabel: hightime {time = 8.h}
52+
withLabel: onedaytime {time = 24.h}
53+
withLabel: twodaytime {time = 48.h}
5454

5555
withLabel: lowcpu { cpus = 5 }
5656
withLabel: midcpu { cpus = 15 }

‎scripts/run_all.sh‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
set -e
22

33
# datasets=( 'replogle' 'op' 'nakatake' 'adamson' 'norman' 'xaira_HEK293T' 'xaira_HCT116' 'parsebioscience' 'ibd_uc' 'ibd_cd' '300BCG' ) #'replogle' 'op' 'nakatake' 'adamson' 'norman' 'xaira_HEK293T' 'xaira_HCT116' 'parsebioscience' 'ibd_uc' 'ibd_cd' '300BCG') #
4-
datasets=( 'op' ) #'replogle' 'op' 'nakatake' 'adamson' 'norman' 'xaira_HEK293T' 'xaira_HCT116' 'parsebioscience' 'ibd_uc' 'ibd_cd' '300BCG') #
4+
datasets=( 'op' 'replogle' ) #'replogle' 'op' 'nakatake' 'adamson' 'norman' 'xaira_HEK293T' 'xaira_HCT116' 'parsebioscience' 'ibd_uc' 'ibd_cd' '300BCG') #
55
run_local=false # set to true to run locally, false to run on AWS
66

77
run_grn_inference=false

‎src/metrics/sem/helper.py‎

Lines changed: 61 additions & 44 deletions
Original file line numberDiff line numberDiff line change
@@ -383,7 +383,9 @@ def main(par):
383383
gene_mask = np.logical_or(np.any(A, axis=1), np.any(A, axis=0))
384384
in_degrees = np.sum(A != 0, axis=0)
385385
out_degrees = np.sum(A != 0, axis=1)
386-
idx = np.argsort(np.maximum(out_degrees, in_degrees))[:-par['n_top_genes']]
386+
# n_genes = par['n_top_genes']
387+
n_genes = 3000
388+
idx = np.argsort(np.maximum(out_degrees, in_degrees))[:-n_genes]
387389
gene_mask[idx] = False
388390
X = X[:, gene_mask]
389391
X = X.toarray() if isinstance(X, csr_matrix) else X
@@ -445,52 +447,67 @@ def main(par):
445447
# Evaluate inferred GRN
446448
print("\n======== Evaluate inferred GRN ========")
447449
scores = evaluate_grn(X_controls, delta_X, is_train, is_reporter, A, signed=use_signs)
448-
449-
# Evaluate baseline GRN
450-
print("\n======== Evaluate shuffled GRN ========")
451-
scores_baseline = evaluate_grn(X_controls, delta_X, is_train, is_reporter, A_baseline, signed=use_signs)
452-
453-
# Keep only the genes for which both GRNs got a score
454-
mask = ~np.logical_or(np.isnan(scores), np.isnan(scores_baseline))
455-
scores = scores[mask]
456-
scores_baseline = scores_baseline[mask]
457-
458-
rr_all = {}
459-
# Perform rank test between actual scores and baseline
460-
rr_all['spearman'] = float(np.mean(scores))
461-
rr_all['spearman_shuffled'] = float(np.mean(scores_baseline))
462-
if len(scores) == 0:
450+
451+
# Keep only valid scores (non-NaN)
452+
valid_scores = scores[~np.isnan(scores)]
453+
454+
if len(valid_scores) == 0:
463455
# No valid genes to evaluate
464-
df_results = pd.DataFrame({'sem_precision': [np.nan], 'sem_balanced': [np.nan]})
465-
elif np.all(scores - scores_baseline == 0):
466-
# Identical performance (suspicious - likely an error)
467-
print("WARNING: Identical scores detected - possible evaluation error!")
468-
df_results = pd.DataFrame({'sem_precision': [1.0], 'sem_balanced': [0.0]})
456+
print("WARNING: No valid genes to evaluate!")
457+
results = {'sem': [0.0]}
469458
else:
470-
res = wilcoxon(scores - scores_baseline, zero_method='wilcox', alternative='greater')
471-
rr_all['Wilcoxon pvalue'] = float(res.pvalue)
472-
473-
print(rr_all)
459+
# Final score is mean of valid R² scores
460+
final_score = float(np.mean(valid_scores))
474461

475-
eps = 1e-300 # very small number to avoid log(0)
476-
pval_clipped = max(res.pvalue, eps)
462+
print(f"\nMethod: {method_id}")
463+
print(f"SEM score (mean R²): {final_score:.4f}")
464+
print(f"Valid genes evaluated: {len(valid_scores)}/{len(scores)}")
465+
print(f"SEM score (min): {np.min(valid_scores):.4f}")
466+
print(f"SEM score (max): {np.max(valid_scores):.4f}")
477467

478-
# Set to 0 if not significant (p >= 0.05)
479-
if res.pvalue >= 0.05:
480-
score = 0.0
481-
print(f"p-value: {res.pvalue:.6f} (not significant, p >= 0.05)")
482-
print(f"SEM score set to 0")
468+
results = {'sem': [float(final_score)]}
469+
470+
# Evaluate baseline GRN
471+
if False:
472+
print("\n======== Evaluate shuffled GRN ========")
473+
scores_baseline = evaluate_grn(X_controls, delta_X, is_train, is_reporter, A_baseline, signed=use_signs)
474+
475+
# Keep only the genes for which both GRNs got a score
476+
mask = ~np.logical_or(np.isnan(scores), np.isnan(scores_baseline))
477+
scores = scores[mask]
478+
scores_baseline = scores_baseline[mask]
479+
480+
rr_all = {}
481+
# Perform rank test between actual scores and baseline
482+
rr_all['spearman'] = float(np.mean(scores))
483+
rr_all['spearman_shuffled'] = float(np.mean(scores_baseline))
484+
if len(scores) == 0:
485+
raise ValueError("No valid scores to compare between inferred GRN and baseline GRN.")
486+
elif np.all(scores - scores_baseline == 0):
487+
# Identical performance (suspicious - likely an error)
488+
raise ValueError("Identical performance between inferred GRN and baseline GRN - likely an error.")
483489
else:
484-
# Compute final score
485-
score = -np.log10(pval_clipped)
486-
print(f"p-value: {res.pvalue:.6f} (significant)")
487-
488-
print(f"Final score: {score}")
489-
490-
results = {
491-
'sem_precision': [float(np.log2(np.mean(scores) / (np.mean(scores_baseline) + 1e-6)))],
492-
'sem': [float(score)]
493-
}
494-
495-
df_results = pd.DataFrame(results)
490+
res = wilcoxon(scores - scores_baseline, zero_method='wilcox', alternative='greater')
491+
rr_all['Wilcoxon pvalue'] = float(res.pvalue)
492+
493+
print(rr_all)
494+
495+
eps = 1e-300 # very small number to avoid log(0)
496+
pval_clipped = max(res.pvalue, eps)
497+
498+
# Set to 0 if not significant (p >= 0.05)
499+
if res.pvalue >= 0.05:
500+
score = 0.0
501+
print(f"p-value: {res.pvalue:.6f} (not significant, p >= 0.05)")
502+
print(f"SEM score set to 0")
503+
else:
504+
# Compute final score
505+
score = -np.log10(pval_clipped)
506+
print(f"p-value: {res.pvalue:.6f} (significant)")
507+
508+
print(f"Final score: {score}")
509+
results['sem_precision'] = [float(np.log2(np.mean(scores) / (np.mean(scores_baseline) + 1e-6)))]
510+
results['sem_n'] = [float(score)]
511+
512+
df_results = pd.DataFrame(results)
496513
return df_results

‎src/metrics/sem/script.py‎

Lines changed: 5 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -18,18 +18,11 @@
1818
'genes_n': 5000
1919
}
2020
## VIASH END
21-
22-
import argparse
23-
parser = argparse.ArgumentParser()
24-
parser.add_argument('--prediction', type=str, help='Path to the predicted GRN in h5ad format')
25-
parser.add_argument('--evaluation_data', type=str, help='Path to the evaluation data in h5ad format')
26-
parser.add_argument('--score', type=str)
27-
parser.add_argument('--layer', type=str, default='lognorm', help='Layer in the h5ad file to use')
28-
parser.add_argument('--num_workers', type=int, default=20, help='Number of workers to use')
29-
21+
run_local = False
3022
try:
3123
sys.path.append(meta["resources_dir"])
3224
except:
25+
run_local=True
3326
meta = {
3427
"resources_dir":'src/metrics/sem/',
3528
"util_dir": 'src/utils',
@@ -39,14 +32,11 @@
3932
sys.path.append(meta["util_dir"])
4033
sys.path.append(meta["helper_dir"])
4134
from helper import main as main_sem
42-
from util import format_save_score
43-
35+
from util import format_save_score, parse_args
4436

37+
if run_local:
38+
par = parse_args(par)
4539

46-
args = parser.parse_args()
47-
for key, value in vars(args).items():
48-
if value is not None:
49-
par[key] = value
5040

5141
if __name__ == "__main__":
5242

‎src/metrics/vc/helper.py‎

Lines changed: 36 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -486,45 +486,42 @@ def main(par):
486486
j = gene_dict[target]
487487
A[i, j] = float(weight)
488488

489-
# Only consider the genes that are actually present in the inferred GRN
490-
gene_mask = np.logical_or(np.any(A, axis=1), np.any(A, axis=0))
491-
X = X[:, gene_mask]
492-
A = A[gene_mask, :][:, gene_mask]
493-
gene_names = gene_names[gene_mask]
494-
495-
print(f"Using {len(gene_names)} genes present in the GRN")
496-
497-
# Additional memory-aware gene filtering for very large GRNs
498-
MAX_GENES_FOR_MEMORY = par['n_top_genes'] # Reduced further to avoid memory issues
499-
if len(gene_names) > MAX_GENES_FOR_MEMORY:
500-
print(f"Too many genes ({len(gene_names)}) for memory. Selecting top {MAX_GENES_FOR_MEMORY} by GRN connectivity.")
489+
# Gene filtering based on n_top_genes parameter
490+
if par['n_top_genes'] == -1:
491+
# Use all genes from evaluation data
492+
# A already has zeros for genes without GRN connections
493+
print(f"Using all {len(gene_names)} genes from evaluation data (including those without GRN connections)")
494+
else:
495+
# Filter to genes present in the inferred GRN
496+
gene_mask = np.logical_or(np.any(A, axis=1), np.any(A, axis=0))
501497

502-
# Select genes with highest connectivity in the GRN
503-
gene_connectivity = np.sum(np.abs(A), axis=0) + np.sum(np.abs(A), axis=1)
504-
top_gene_indices = np.argsort(gene_connectivity)[-MAX_GENES_FOR_MEMORY:]
498+
# Additional memory-aware gene filtering for very large GRNs
499+
MAX_GENES_FOR_MEMORY = par['n_top_genes']
500+
if np.sum(gene_mask) > MAX_GENES_FOR_MEMORY:
501+
print(f"Too many genes with GRN connections ({np.sum(gene_mask)}) for memory. Selecting top {MAX_GENES_FOR_MEMORY} by GRN connectivity.")
502+
503+
# Select genes with highest connectivity in the GRN
504+
gene_connectivity = np.sum(np.abs(A), axis=0) + np.sum(np.abs(A), axis=1)
505+
# Set connectivity to 0 for genes not in mask
506+
gene_connectivity[~gene_mask] = 0
507+
top_gene_indices = np.argsort(gene_connectivity)[-MAX_GENES_FOR_MEMORY:]
508+
gene_mask = np.zeros(len(gene_names), dtype=bool)
509+
gene_mask[top_gene_indices] = True
505510

506-
X = X[:, top_gene_indices]
507-
A = A[top_gene_indices, :][:, top_gene_indices]
508-
gene_names = gene_names[top_gene_indices]
511+
# Apply the gene mask
512+
X = X[:, gene_mask]
513+
A = A[gene_mask, :][:, gene_mask]
514+
gene_names = gene_names[gene_mask]
509515

510-
print(f"Final: Using {len(gene_names)} most connected genes for evaluation")
516+
print(f"Using {len(gene_names)} genes (filtered by GRN connectivity)")
511517

512518
# Remove self-regulations
513519
np.fill_diagonal(A, 0)
514520

515-
# Create baseline model
516-
A_baseline = np.copy(A)
517-
for j in range(A.shape[1]):
518-
np.random.shuffle(A[:j, j])
519-
np.random.shuffle(A[j+1:, j])
520-
assert np.any(A_baseline != A)
521-
522521
# Mapping between gene expression profiles and their matched negative controls
523-
524522
control_map, _ = create_control_matching(are_controls, match_groups)
525523
loose_control_map, _ = create_control_matching(are_controls, loose_match_groups)
526524

527-
528525
ss_res = 0
529526
ss_tot = 0
530527
cv = GroupKFold(n_splits=5)
@@ -561,22 +558,26 @@ def main(par):
561558
test_data_loader = torch.utils.data.DataLoader(test_dataset, batch_size=512)
562559

563560
# Evaluate inferred GRN
561+
print(f"\n======== Fold {i+1}: Evaluate inferred GRN ========")
564562
res = evaluate(A, train_data_loader, test_data_loader, n_perturbations)
565563
ss_res = ss_res + res[0]
566564
ss_tot = ss_tot + res[1]
567565

568-
# Evaluate baseline GRN (shuffled target genes)
569-
#ss_tot = ss_tot + evaluate(A_baseline, train_data_loader, test_data_loader, n_perturbations)
570-
571566
r2 = 1 - ss_res / ss_tot
572567

573-
final_score = np.mean(np.clip(r2, 0, 1))
574-
print(f"Method: {method_id}")
575-
print(f"R2: {final_score}")
568+
# Compute scores per gene
569+
r2_per_gene = np.clip(r2, 0, 1)
570+
571+
# Final score is mean R2 across genes
572+
final_score = np.mean(r2_per_gene)
573+
574+
print(f"\nMethod: {method_id}")
575+
print(f"R2 (mean): {final_score:.4f}")
576+
print(f"R2 (min): {np.min(r2_per_gene):.4f}")
577+
print(f"R2 (max): {np.max(r2_per_gene):.4f}")
576578

577579
results = {
578580
'vc': [float(final_score)]
579-
580581
}
581582

582583
df_results = pd.DataFrame(results)

‎src/metrics/vc/run_local.sh‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -17,7 +17,7 @@ mkdir -p "$save_dir"
1717

1818
# datasets to process
1919
datasets=( "replogle" "xaira_HEK293T" "xaira_HCT116" "nakatake" "norman" "adamson" 'parsebioscience' 'op' "300BCG" 'ibd_uc' 'ibd_cd') #"300BCG" "ibd" 'parsebioscience', 'xaira_HEK293T'
20-
# datasets=( "op" "300BCG" "parsebioscience" "ibd" )
20+
datasets=( "op" "300BCG" "parsebioscience" "ibd" )
2121
# methods to process
2222
methods=( "pearson_corr" "positive_control" "negative_control" "ppcor" "portia" "scenic" "grnboost" "scprint" "scenicplus" "celloracle" "scglue" "figr" "granie")
2323
# methods=( "grnboost")

‎src/utils/config.py‎

Lines changed: 5 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -93,7 +93,6 @@
9393
'norman': ['regression', 'ws_distance', 'tf_binding', 'sem', 'gs_recovery'],
9494
'nakatake': ['regression', 'sem', 'gs_recovery'],
9595
'op': ['regression', 'vc', 'rc_tf_act', 'tf_binding', 'sem', 'gs_recovery'],
96-
# 'op': ['regression', 'tf_binding', 'gs_recovery'],
9796
'300BCG': ['regression', 'vc', 'rc_tf_act', 'tf_binding', 'sem', 'gs_recovery'],
9897
'ibd_uc': ['regression', 'vc', 'tf_binding', 'sem', 'gs_recovery'],
9998
'ibd_cd': ['regression', 'vc', 'tf_binding', 'sem', 'gs_recovery'],
@@ -111,28 +110,21 @@
111110
METRICS = [
112111
'r_precision', 'r_recall', 'r_f1',
113112
'ws_precision', 'ws_recall', 'ws_f1',
114-
'vc',
115-
'sem',
113+
'vc', 'vc_raw', 'vc_precision',
114+
'sem', 'sem_precision', 'sem_raw',
116115
't_rec_precision', 't_rec_recall', 't_rec_f1',
117-
118-
'rc_tf_act',
119-
120-
'anchor_regression_raw',
121-
116+
'rc_tf_act',
122117
'tfb_precision', 'tfb_recall', 'tfb_f1',
123118
'gs_precision', 'gs_recall', 'gs_f1',
124119
]
125120

126121
FINAL_METRICS = [
127122
'r_precision', 'r_recall',
128-
'ws_precision', 'ws_recall',
129123
'vc',
130-
'sem',
124+
'sem',
125+
'ws_precision', 'ws_recall',
131126
't_rec_precision', 't_rec_recall',
132-
133127
'rc_tf_act',
134-
135-
'anchor_regression_raw',
136128
'tfb_f1',
137129
'gs_f1',
138130
]

‎src/utils/dataset_config.env‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -19,7 +19,7 @@ METRICS_replogle="regression,ws_distance,tf_recovery,tf_binding,sem,gs_recovery"
1919
METRICS_adamson="regression,tf_binding,sem,gs_recovery"
2020
METRICS_norman="regression,ws_distance,tf_binding,sem,gs_recovery"
2121
METRICS_nakatake="regression,sem,gs_recovery"
22-
METRICS_op="tf_binding,gs_recovery"
22+
METRICS_op="regression,vc,rc_tf_act,tf_binding,sem,gs_recovery"
2323
METRICS_300BCG="regression,vc,rc_tf_act,tf_binding,sem,gs_recovery"
2424
METRICS_ibd_uc="regression,vc,tf_binding,sem,gs_recovery"
2525
METRICS_ibd_cd="regression,vc,tf_binding,sem,gs_recovery"

‎test.ipynb‎

Lines changed: 3 additions & 95 deletions
Large diffs are not rendered by default.

0 commit comments

Comments
 (0)