Skip to content

Commit

Permalink
positive mask added
Browse files Browse the repository at this point in the history
  • Loading branch information
@leeyoyohku
leeyoyohku committed Dec 12, 2022
1 parent 5a855d4 commit aa17624
Show file tree
Hide file tree
Showing 3 changed files with 25 additions and 25 deletions.
3 changes: 2 additions & 1 deletion spatialdm/plottings.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,8 @@
from sklearn import linear_model
import scipy.stats as stats
import seaborn as sns
from utils import compute_pathway
#from utils import compute_pathway
from .utils import *
import holoviews as hv
from holoviews import opts, dim
from bokeh.io import output_file, show
Expand Down
45 changes: 22 additions & 23 deletions spatialdm/utils.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -149,7 +149,6 @@ def pair_selection_matrix(adata, n_perm, sel_ind, method):
(adata.obsp['weight'] @ L_mat_use[_idx, n_short_lri:]) * R_mat_use[_idx, n_short_lri:], axis=0)



def spot_selection_matrix(adata, ligand, receptor, ind, n_perm, method):
# local variables (only live in this function scope)
# ligand = adata.uns['ligand'].loc[sel_ind]
Expand All @@ -167,71 +166,71 @@ def spot_selection_matrix(adata, ligand, receptor, ind, n_perm, method):
if adata.uns['mean'] == 'geometric':
from scipy.stats.mstats import gmean
L1 = [pd.Series(x[0]).dropna().values for x in ligand.values]
L_mat0 = [adata[:, L1[l]].X.astype(np.float)[:,0] for l in range(len(L1))]
for i,k in enumerate(ligand.index):
L_mat0 = [adata.raw[:, L1[l]].X.astype(np.float)[:, 0] for l in range(len(L1))]
for i, k in enumerate(ligand.index):
if len(ligand.loc[k].dropna()) > 1:
if adata.uns['mean'] == 'geometric':
L_mat0[i] = gmean(adata[:, ligand.loc[k].dropna()].X, axis=1)
L_mat0[i] = gmean(adata.raw[:, ligand.loc[k].dropna()].X, axis=1)
else:
L_mat0[i] = adata[:, ligand.loc[k].dropna()].X.mean(1)
L_mat0[i] = adata.raw[:, ligand.loc[k].dropna()].X.mean(1)

# averaged receptor values
R1 = [pd.Series(x[0]).dropna().values for x in receptor.values]
R_mat0 = [adata[:, R1[r]].X.astype(np.float)[:, 0] for r in range(len(R1))]
R_mat0 = [adata.raw[:, R1[r]].X.astype(np.float)[:, 0] for r in range(len(R1))]
for i, k in enumerate(receptor.index):
if len(receptor.loc[k].dropna()) > 1:
if adata.uns['mean'] == 'geometric':
R_mat0[i] = gmean(adata[:, receptor.loc[k].dropna()].X, axis=1)
R_mat0[i] = gmean(adata.raw[:, receptor.loc[k].dropna()].X, axis=1)
else:
R_mat0[i] = adata[:, receptor.loc[k].dropna()].X.mean(1)
R_mat0[i] = adata.raw[:, receptor.loc[k].dropna()].X.mean(1)
n_short_lri = (adata.uns['geneInter'].loc[ligand.index, 'annotation'] \
!='Secreted Signaling').sum()
!= 'Secreted Signaling').sum()
ranges = [np.arange(n_short_lri), np.arange(n_short_lri, len(L1))]
weight_matrices = [adata.obsp['nearest_neighbors'], adata.obsp['weight']]
N = adata.shape[0]
L_mat0 = np.array(L_mat0)
R_mat0 = np.array(R_mat0)
pos = np.zeros((N, len(ligand)))

for r,weight_matrix in zip(ranges, weight_matrices):
for r, weight_matrix in zip(ranges, weight_matrices):
R_mat = R_mat0[r].T
L_mat = L_mat0[r].T
pos = (L_mat > 0) + (R_mat > 0)
pos[:, r] = (L_mat > 0) + (R_mat > 0)
R_mat_use = _standardise(R_mat, Local=True, axis=0)
L_mat_use = _standardise(L_mat, Local=True, axis=0)
wij_sq = (weight_matrix ** 2).sum(1)
rbf_d = csc_matrix(weight_matrix)

adata.uns['local_stat']['local_I'][:,r] = (rbf_d @ R_mat_use) * L_mat_use
adata.uns['local_stat']['local_I_R'][:,r] = (rbf_d @ L_mat_use) * R_mat_use
adata.uns['local_stat']['local_I'][:, r] = (rbf_d @ R_mat_use) * L_mat_use
adata.uns['local_stat']['local_I_R'][:, r] = (rbf_d @ L_mat_use) * R_mat_use
## Calculate p values
if method in ['both', 'z-score']:
norm_res1 = [stats.norm.fit(L_mat_use[:, i]) for i in range(L_mat_use.shape[1])]
norm_res2 = [stats.norm.fit(R_mat_use[:, i]) for i in range(R_mat_use.shape[1])]
norm_res1 = np.array(norm_res1)
norm_res2 = np.array(norm_res2)
mu1_ls, std_L_ls = norm_res1[:,0], norm_res1[:,1]
mu2_ls, std_R_ls = norm_res2[:,0], norm_res2[:,1]
mu1_ls, std_L_ls = norm_res1[:, 0], norm_res1[:, 1]
mu2_ls, std_R_ls = norm_res2[:, 0], norm_res2[:, 1]
sigma1_sq_ls = [(std1 * N / (N - 1)) for std1 in std_L_ls]
sigma2_sq_ls = [(std2 * N / (N - 1)) for std2 in std_R_ls]
std_ls = [compute_var_local(sigma1_sq, sigma2_sq, wij_sq, N) \
for (sigma1_sq, sigma2_sq) in zip(sigma1_sq_ls, sigma2_sq_ls)]
adata.uns['local_z'][r] = (adata.uns['local_stat']['local_I'][:,r] + \
adata.uns['local_stat']['local_I_R'][:,r]).T / std_ls
for (sigma1_sq, sigma2_sq) in zip(sigma1_sq_ls, sigma2_sq_ls)]
adata.uns['local_z'][r] = (adata.uns['local_stat']['local_I'][:, r] + \
adata.uns['local_stat']['local_I_R'][:, r]).T / std_ls
adata.uns['local_z_p'][r] = stats.norm.sf(adata.uns['local_z'][r])
#adata.uns['local_z_p'] = np.where(pos.T == False, 1, adata.uns['local_z_p'])

if method in ['both', 'permutation']:
for i in tqdm(range(n_perm)):
_idx = np.random.permutation(L_mat.shape[0])
adata.uns['local_stat']['local_permI'][r, i,:] = ((rbf_d @ R_mat_use[_idx, :]) * L_mat_use).T
adata.uns['local_stat']['local_permI_R'][r, i,:] = ((rbf_d @ L_mat_use[_idx, :]) * R_mat_use).T
#adata.uns['local_perm_p'] = np.where(pos.T == False, 1, adata.uns['local_perm_p'])
adata.uns['local_stat']['local_permI'][r, i, :] = ((rbf_d @ R_mat_use[_idx, :]) * L_mat_use).T
adata.uns['local_stat']['local_permI_R'][r, i, :] = ((rbf_d @ L_mat_use[_idx, :]) * R_mat_use).T
try:
adata.uns['local_z_p'] = np.where(pos.T == False, 1, adata.uns['local_z_p'])
adata.uns['local_z_p'] = pd.DataFrame(adata.uns['local_z_p'], index=ind)
adata.uns['local_perm_p'] = (np.expand_dims(adata.uns['local_stat']['local_I'].T + \
adata.uns['local_stat']['local_I_R'].T, 1) <= \
(adata.uns['local_stat']['local_permI'] + adata.uns['local_stat'][
'local_permI_R'])).sum(1) / n_perm
adata.uns['local_perm_p'] = np.where(pos.T == False, 1, adata.uns['local_perm_p'])
adata.uns['local_perm_p'] = pd.DataFrame(adata.uns['local_perm_p'], index=ind)
except OSError as e:
if e.errno != e.errno:
Expand Down
2 changes: 1 addition & 1 deletion spatialdm/version.py
View file
Original file line number Diff line number Diff line change
@@ -1 +1 @@
__version__ = '0.0.2'
__version__ = '0.0.5'

0 comments on commit aa17624

Please sign in to comment.