Skip to content

Commit 83f383e

Browse files
ghar1821ghar1821
committed
update get obs and var function for R
1 parent dca6d2e commit 83f383e

4 files changed

Lines changed: 123 additions & 92 deletions

File tree

‎src/control_methods/perfect_integration/script.py‎

Lines changed: 10 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -3,28 +3,27 @@
33
## VIASH START
44
# The following code has been auto-generated by Viash.
55
par = {
6-
'input_unintegrated': 'resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/unintegrated.h5ad',
7-
'output_integrated_left': 'resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/integrated.h5ad',
8-
'output_integrated_right': 'resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/integrated.h5ad'
9-
}
10-
meta = {
11-
'name': 'perfect_integration'
6+
"input_unintegrated": "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/unintegrated.h5ad",
7+
"output_integrated_left": "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/integrated_left.h5ad",
8+
"output_integrated_right": "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/integrated_right.h5ad",
129
}
10+
meta = {"name": "perfect_integration"}
1311

1412
## VIASH END
1513

1614
print("Reading input files", flush=True)
1715
adata = ad.read_h5ad(par["input_unintegrated"])
1816

1917
print("Extracting and splitting unintegrated data", flush=True)
20-
#split 1
21-
adata_left = adata[(adata.obs.is_control>0) | (adata.obs.batch==1)]
18+
19+
# split 1
20+
adata_left = adata[(adata.obs.is_control > 0) | (adata.obs.batch == 1)]
2221
integrated_left = adata_left.layers["preprocessed"]
23-
#split 2 == split 1 in this case
2422

23+
# split 2 == split 1 in this case
2524

2625
print("Write output AnnData to file", flush=True)
27-
#split 1
26+
# split 1
2827
output_left = ad.AnnData(
2928
obs=adata_left.obs[[]],
3029
var=adata_left.var[[]],
@@ -35,7 +34,7 @@
3534
"parameters": {},
3635
},
3736
)
38-
#split 2
37+
# split 2
3938
output_right = ad.AnnData(
4039
obs=adata_left.obs[[]],
4140
var=adata_left.var[[]],

‎src/methods/harmonypy/script.py‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -4,8 +4,8 @@
44

55
## VIASH START
66
par = {
7-
"input": "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/unintegrated_censored.h5ad",
8-
"output": "output.h5ad",
7+
"input": "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/censored_right.h5ad",
8+
"output": "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/output_harmony_right.h5ad",
99
}
1010
meta = {"name": "harmonypy"}
1111
## VIASH END
@@ -30,7 +30,7 @@
3030
)
3131

3232
# have to add in the uncorrected markers as well
33-
uncorrected_data = adata[:, markers_not_correct].layers['preprocessed']
33+
uncorrected_data = adata[:, markers_not_correct].layers["preprocessed"]
3434
out_matrix = np.concatenate([out.Z_corr.transpose(), uncorrected_data], axis=1)
3535
out_var_idx = np.concatenate([markers_to_correct, markers_not_correct])
3636

‎src/utils/helper_functions.R‎

Lines changed: 83 additions & 62 deletions
Original file line numberDiff line numberDiff line change
@@ -1,61 +1,82 @@
11
# NOTE: These helper functions are ports of the original Python functions in 'helper_functions.py'
22

3-
#' Adds annotations (.var and .obs) from the unintegrated dataset to the
4-
#' integrated dataset. In the case of the control method "perfect_integration",
5-
#' the function will fetch annotations from the validation dataset instead.
3+
library(dplyr)
4+
requireNamespace("anndataR", quietly = TRUE)
5+
6+
#' Adds annotations (.var and .obs) from the unintegrated data to the
7+
#' integrated dataset.
8+
#' In the case of the control method "perfect_integration",
9+
#' the function will fetch the batch label from the unintegrated data
10+
#' based on the split.
11+
#' i.e., if in split 1, donor 3-5 is from batch 2, then the batch label for that split
12+
#' will be changed from batch 1 to batch 2.
613
#'
7-
#' @param i_adata AnnData object, batch-integrated dataset
8-
#' @param v_adata AnnData object, validation dataset
14+
#' @param s1_adata AnnData object, integrated data from split 1
15+
#' @param s2_adata AnnData object, integrated data from split 2
916
#' @param u_adata AnnData object, unintegrated dataset
1017
#' @return AnnData object with .var and .obs added
11-
get_obs_var_for_integrated <- function(i_adata, v_adata, u_adata) {
12-
if (i_adata$uns$method_id == "perfect_integration_horizontal") {
13-
if (i_adata$n_obs != v_adata$n_obs) {
14-
stop(
15-
"The number of cells in the integrated (perfect_integration_horizontal) ",
16-
"and validation datasets do not match"
17-
)
18-
}
19-
i_adata$obs <- v_adata$obs[rownames(i_adata), , drop = FALSE]
20-
i_adata$var <- v_adata$var[colnames(i_adata), , drop = FALSE]
21-
} else if (i_adata$uns$method_id == "perfect_integration_vertical") {
22-
comb_adata <- anndata::concat(list(v_adata, u_adata))
23-
# subset to just batch 1
24-
# Check if 'batch' column exists
25-
if (!"batch" %in% colnames(comb_adata$obs)) {
26-
stop(
27-
"Column 'batch' not found in comb_adata$obs for ",
28-
"perfect_integration_vertical."
29-
)
30-
}
31-
comb_adata <- comb_adata[comb_adata$obs$batch == 1, ]
18+
#'
19+
get_obs_var_for_integrated <- function(s1_adata, s2_adata, u_adata) {
20+
21+
s1_adata$obs <- u_adata$obs[s1_adata$obs_names, ]
22+
s2_adata$obs <- u_adata$obs[s2_adata$obs_names, ]
23+
s1_adata$var <- u_adata$var[s1_adata$var_names, ]
24+
s2_adata$var <- u_adata$var[s2_adata$var_names, ]
3225

33-
if (i_adata$n_obs != comb_adata$n_obs) {
34-
stop(
35-
"The number of cells in the integrated (perfect_integration_vertical) ",
36-
"and validation + unintegrated datasets do not match."
37-
)
38-
}
39-
i_adata$obs <- comb_adata$obs[rownames(i_adata), , drop = FALSE]
40-
i_adata$var <- v_adata$var[colnames(i_adata), , drop = FALSE]
41-
} else {
42-
if (i_adata$n_obs != u_adata$n_obs) {
43-
stop(
44-
"The number of cells in the integrated and unintegrated datasets do not match"
45-
)
46-
}
47-
# Compare obs_names for ordering
48-
if (!all(rownames(i_adata) == rownames(u_adata))) {
49-
warning(
50-
"The cell ordering in the integrated and unintegrated datasets do not match"
51-
)
26+
# if integrated data came from perfect integration, change the batch labels of the samples
27+
# everything is from batch 1, but some samples need to be labelled to come from batch 2
28+
if (s1_adata$uns["method_id"] == "perfect_integration") {
29+
cat(
30+
"Control method 'perfect_integration' detected. Changing batch labels for split 2.\n"
31+
)
32+
33+
cat("Computing new batch labels\n")
34+
# mutate is needed as donors that are used for controls, we won't have the mapping
35+
s1_adata_new_batch_labels <- get_batch_label_perfect_integration(
36+
u_adata = u_adata,
37+
i_adata = s1_adata,
38+
split_id = 1
39+
)
40+
41+
s2_adata_new_batch_labels <- get_batch_label_perfect_integration(
42+
u_adata = u_adata,
43+
i_adata = s2_adata,
44+
split_id = 2
45+
)
46+
47+
cat("Attaching new batch labels\n")
48+
s1_adata$obs$batch <- s1_adata_new_batch_labels$new_batch_label
49+
s2_adata$obs$batch <- s2_adata_new_batch_labels$new_batch_label
5250
}
5351

54-
i_adata$obs <- u_adata$obs[rownames(i_adata), , drop = FALSE]
55-
i_adata$var <- u_adata$var[colnames(i_adata), , drop = FALSE]
56-
}
52+
return(list(
53+
"s1_adata" = s1_adata,
54+
"s2_adata" = s2_adata
55+
))
56+
}
57+
58+
#' Helper function to get the batch label for perfect integration.
59+
#' First, get donor batch map for a given split.
60+
#' Then apply the map to the integrated data.
61+
#'
62+
#' @param u_adata AnnData object, unintegrated dataset.
63+
#' @param i_adata AnnData object, integrated data.
64+
#' @param split_id numeric, split id of the integrated data.
65+
#'
66+
#' @return a dataframe with donor and new batch label
67+
#'
68+
get_batch_label_perfect_integration <- function(u_adata, i_adata, split_id) {
69+
actual_donor_batch_map <- unique(
70+
u_adata$obs[(u_adata$obs$split == split_id),
71+
c("donor", "batch")]
72+
)
73+
# mutate is needed as donors that are used for controls, we won't have the mapping
74+
i_adata_new_batch_labels <- i_adata$obs[, c("donor", "batch")] %>%
75+
left_join(actual_donor_batch_map, by="donor", suffix = c("_old", "_new")) %>%
76+
mutate(new_batch_label = ifelse(is.na(batch_new), batch_old, batch_new)) %>%
77+
select(donor, new_batch_label)
5778

58-
i_adata
79+
return(i_adata_new_batch_labels)
5980
}
6081

6182
#' Subsets the anndata object to remove the control cells.
@@ -64,12 +85,12 @@ get_obs_var_for_integrated <- function(i_adata, v_adata, u_adata) {
6485
#' @param adata AnnData object
6586
#' @return AnnData object with cells from control samples removed
6687
subset_nocontrols <- function(adata) {
67-
if (!"is_control" %in% colnames(adata$obs)) {
68-
stop("The column 'is_control' is not present in the adata object.")
69-
}
88+
if (!"is_control" %in% colnames(adata$obs)) {
89+
stop("The column 'is_control' is not present in the adata object.")
90+
}
7091

71-
# Subset the adata to remove cells where is_control != 0
72-
adata[adata$obs$is_control == 0, ]
92+
# Subset the adata to remove cells where is_control != 0
93+
adata[adata$obs$is_control == 0, ]
7394
}
7495

7596
#' Subsets the anndata object to only include markers that need to be
@@ -79,7 +100,7 @@ subset_nocontrols <- function(adata) {
79100
#' @param adata AnnData object
80101
#' @return AnnData object with only the markers to correct
81102
subset_markers_tocorrect <- function(adata) {
82-
adata[, adata$var$to_correct]
103+
adata[, adata$var$to_correct]
83104
}
84105

85106
#' Subsets the anndata object to remove all cells where the marker is not
@@ -89,12 +110,12 @@ subset_markers_tocorrect <- function(adata) {
89110
#' @param adata AnnData object
90111
#' @return AnnData object with only the labeled cells
91112
remove_unlabelled <- function(adata) {
92-
if (!"cell_type" %in% colnames(adata$obs)) {
93-
stop("The column 'cell_type' is not present in the adata object.")
94-
}
113+
if (!"cell_type" %in% colnames(adata$obs)) {
114+
stop("The column 'cell_type' is not present in the adata object.")
115+
}
95116

96-
# Convert to lowercase and filter out "unlabelled" and "unlabeled"
97-
is_unlabelled <- tolower(adata$obs$cell_type) %in%
98-
c("unlabelled", "unlabeled")
99-
adata[!is_unlabelled, ]
117+
# Convert to lowercase and filter out "unlabelled" and "unlabeled"
118+
is_unlabelled <- tolower(adata$obs$cell_type) %in%
119+
c("unlabelled", "unlabeled")
120+
adata[!is_unlabelled, ]
100121
}

‎src/utils/helper_functions.py‎

Lines changed: 27 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,10 @@
11
import anndata as ad
2+
import numpy as np
23

34

45
def get_obs_var_for_integrated(
56
s1_adata: ad.AnnData, s2_adata: ad.AnnData, u_adata: ad.AnnData
6-
) :
7+
):
78
"""
89
Fetch annotations (.var and .obs) from the unintegrated dataset to the integrated datasets (left and right).
910
In the case of the control method 'perfect_integration', annotations are fetched only from batch 1.
@@ -14,33 +15,37 @@ def get_obs_var_for_integrated(
1415
u_adata: AnnData object, unintegrated dataset
1516
1617
Outputs:
17-
s1_adata: AnnData object, split == 1 dataset with annotations
18-
s2_adata: AnnData object, split == 2 dataset with annotations
18+
s1_adata: AnnData object, split == 1 dataset with annotations
19+
s2_adata: AnnData object, split == 2 dataset with annotations
1920
"""
2021

2122
s1_adata.obs = u_adata.obs.loc[s1_adata.obs_names]
2223
s1_adata.var = u_adata.var.loc[s1_adata.var_names]
2324
s2_adata.obs = u_adata.obs.loc[s2_adata.obs_names]
2425
s2_adata.var = u_adata.var.loc[s2_adata.var_names]
2526

26-
if s1_adata.uns['method_id'] == 'perfect_integration':
27+
if s1_adata.uns["method_id"] == "perfect_integration":
2728
print(
2829
"Control method 'perfect_integration' detected. Changing batch labels for split 2"
2930
)
3031

31-
#Apply mapping to all non-control cells of split 1 and split 2 (+ change the split label)
32+
# Apply mapping to all non-control cells of split 1 and split 2 (+ change the split label)
3233
split_dict_s1 = get_donor_batch_map(u_adata, split_of_interest=1)
33-
s1_adata.obs.loc[(u_adata.obs.is_control==0), 'batch'] = s1_adata.obs['donor'].map(split_dict_s1)
34-
s1_adata.obs.loc[(u_adata.obs.is_control==0), 'split'] = 1
34+
s1_adata.obs.loc[(u_adata.obs.is_control == 0), "batch"] = s1_adata.obs[
35+
"donor"
36+
].map(split_dict_s1)
37+
s1_adata.obs.loc[(u_adata.obs.is_control == 0), "split"] = 1
3538

3639
split_dict_s2 = get_donor_batch_map(u_adata, split_of_interest=2)
37-
s2_adata.obs.loc[(u_adata.obs.is_control==0), 'batch'] = s2_adata.obs['donor'].map(split_dict_s2)
38-
s2_adata.obs.loc[(u_adata.obs.is_control==0), 'split'] = 2
40+
s2_adata.obs.loc[(u_adata.obs.is_control == 0), "batch"] = s2_adata.obs[
41+
"donor"
42+
].map(split_dict_s2)
43+
s2_adata.obs.loc[(u_adata.obs.is_control == 0), "split"] = 2
3944

4045
return s1_adata, s2_adata
4146

42-
def get_donor_batch_map(
43-
u_adata: ad.AnnData, split_of_interest: int) -> dict:
47+
48+
def get_donor_batch_map(u_adata: ad.AnnData, split_of_interest: int) -> dict:
4449
"""
4550
Create a dictionary that represent the correct donor/batch mapping for a split of interest.
4651
Note: This helper function is only meant to be used for 'perfect_integration' control method.
@@ -55,17 +60,23 @@ def get_donor_batch_map(
5560
import pandas as pd
5661

5762
split_dict = (
58-
u_adata[(u_adata.obs.is_control==0) & (u_adata.obs.split==split_of_interest)].obs
59-
.groupby('donor')['batch']
63+
u_adata[
64+
(u_adata.obs.is_control == 0) & (u_adata.obs.split == split_of_interest)
65+
]
66+
.obs.groupby("donor")["batch"]
6067
.apply(pd.Series.unique)
6168
.apply(list)
6269
.to_dict()
6370
)
6471

65-
#Safeguard to ensure that each donor has a unique batch in the split
66-
assert False not in [True if np.unique(el).size == 1 else False for el in split_dict.values()], "Donor/Batch mapping is ambiguous. Some donors have multiple batches in the same split."
72+
# Safeguard to ensure that each donor has a unique batch in the split
73+
assert False not in [
74+
True if np.unique(el).size == 1 else False for el in split_dict.values()
75+
], (
76+
"Donor/Batch mapping is ambiguous. Some donors have multiple batches in the same split."
77+
)
6778

68-
#transform the elements of the dictionary in integers
79+
# transform the elements of the dictionary in integers
6980
split_dict = {k: int(v[0]) for k, v in split_dict.items()}
7081

7182
return split_dict

0 commit comments

Comments
 (0)