Help with the simulations from the paper

[hiraksarkar]

Hi,

Great paper and package. Thanks alot for making the code open source. I am trying to simulate the dataset for benchmarking cell-cell interaction tools following the Figure 1 f,g,h of the paper. I read the an adapted version of SVCA is being used for simulation. However I am unable to find the source code (or scripts for running SVCA) to generate or reproduce the ROC plots from the paper. It would be great to get some pointers to simulate the data as described.

Thanks again,
Best

[leeyoyohku]

Hi @hiraksarkar ,

Thanks for your question and patience. I didn't include simulation codes because it was mostly identical to SVCA. Please see the adapted codes below and refer to SVCA when needed. Note that you need to stack the simulated counts for each scenario afterwards. Also you need to turn the negatively correlated simulations to positive manually. We didn't explore other possibilities to simulate spatial interaction data, so it may not be the best approach especially after a year of submission.

import numpy as np
import scipy as s
import pandas as pd
from svca.models.model1 import Model1
from svca.simulations.from_real import FromRealSimulation
from svca.models.io import *
from svca.util_functions import utils
import sys
from limix.utils.preprocess import covar_rescaling_factor_efficient
from copy import deepcopy


def run(data_dir, protein_index, output_dir, bootstrap_index,
        normalisation='standard', permute=False):
    # reading all data
    ####################################################################
    expression_file = data_dir + 'exp.txt'
    position_file = data_dir+'pos.txt'
    # protein_name, phenotype, X = pr names, pr expr values, positions
    db = pd.read_csv('/CellChatDB_db2.csv',
                     index_col=0)
    protein_names, phenotypes, X = utils.read_data(expression_file,
                                                   position_file)

  # import pdb; pdb.set_trace()
  # protein_name = protein_names[protein_index, :]
  protein_name = db.index[protein_index]
  print(protein_name)
  phenotype = phenotypes[:, (protein_names==protein_name.split("_")[0])[:,0]]
  #import pdb; pdb.set_trace()
  sel = np.arange(phenotypes.shape[1])
  sel =list(sel[pd.Series(protein_names[:,0]).isin(db.loc[protein_name,
                                                ['R1', 'R2']].dropna().values).values])
  kin_from = phenotypes[:, sel] # remove the selected pr
  N_samples = X.shape[0]

  # permuting cells
  if permute:
      perm = np.random.permutation(X.shape[0])
      X = X[perm, :]

  # do simulations
  ####################################################################
  sim = FromRealSimulation(X, phenotype[:,0], kin_from)
  file_prefix = protein_name
  Y_sim = sim.simulate(interactions_size = None)
  np.save(output_dir+ 'nul/'+file_prefix, Y_sim)
  Y_sim = sim.simulate(interactions_size=0.25)
  np.save(output_dir+ 'rescale25/' + file_prefix, Y_sim)
  Y_sim = sim.simulate(interactions_size=0.99)
  np.save(output_dir+ 'full/' + file_prefix, Y_sim)
  Y_sim = sim.simulate(interactions_size=0.75)
  np.save(output_dir+ 'rescale75/' + file_prefix, Y_sim)
  Y_sim = sim.simulate(interactions_size=0.5)
  np.save(output_dir+ 'half/' + file_prefix, Y_sim)
  Y_sim = sim.simulate(interactions_size = 0.99)
  np.save(output_dir+ 'full/'+file_prefix, Y_sim)

  # run model on simulated data
  ####################################################################
  # intrinsic and environmental term
  ####################################################################
  cterms = ['intrinsic', 'environmental']
  model = Model1(Y_sim, X, norm=normalisation, oos_predictions=0., cov_terms=cterms, kin_from=kin_from)
  model.reset_params()
  model.train_gp(grid_size=10)

  file_prefix = protein_name[0] + '_' + str(bootstrap_index) + '_local'
  write_variance_explained(model, output_dir, file_prefix) # local_effects
  write_LL(model, output_dir, file_prefix) # LL

  int_param = model.intrinsic_cov.getParams()
  env_param = model.environmental_cov.getParams()
  noise_param = model.noise_cov.getParams()

  ####################################################################
  # add cell-cell interactions
  ####################################################################
  model.add_cov(['interactions'])

  LL = np.Inf
  for i in range(5):
      if i == 0:
          int_bk = int_param
          local_bk = env_param
          noise_bk = noise_param
          scale_interactions = True
      else:
          int_bk = int_param * s.random.uniform(0.8, 1.2, len(int_param))
          local_bk = env_param * s.random.uniform(0.8, 1.2, len(env_param))
          noise_bk = noise_param * s.random.uniform(0.8, 1.2, len(noise_param))
          scale_interactions = False
      model.set_initCovs({'intrinsic': int_bk,
                      'noise': noise_bk,
                      'environmental':local_bk})
      if scale_interactions:
          model.set_scale_down(['interactions'])
      else:
          model.use_scale_down = False

      model.reset_params()
      model.train_gp(grid_size=10)
      if model.gp.LML() < LL:
          LL = model.gp.LML()
          saved_params = model.gp.getParams()

  model.gp.setParams(saved_params)
  file_prefix = protein_name[0] + '_' + str(bootstrap_index) + '_interactions'
  write_variance_explained(model, output_dir, file_prefix)
  write_r2(model, output_dir, file_prefix)
  write_LL(model, output_dir, file_prefix)
  # write_Ks(model, output_dir, file_prefix)


  if __name__ == '__main__':
    data_dir = sys.argv[1]
    output_dir = sys.argv[2]
    protein_index = int(0)
    bootstrap_index = 1
    normalisation =  'quantile'

    run(data_dir, protein_index, output_dir, bootstrap_index, normalisation)

[hiraksarkar]

Hi @leeyoyohku

Thanks for replying. Assuming "CellChatDB_db2.csv" is a standard cellchat database I assume the R1 and R2 are possible receptors? Am I right? I am actually trying to reproduce the plot in Figure 1f. Can you please share the script for that? Precisely I want to check the ROC score for a few different methods for benchamrking. Thanks!

[yujias424]

Hi @leeyoyohku,

Thank you for presenting this excellent package. I share a similar question to @hiraksarkar regarding the "CellChatDB_db2.csv" file. Should we consider this file as the standard CellChat database, and are R1 and R2 intended to represent possible receptors?

Would it be possible for you to provide the "CellChatDB_db2.csv" file as a reference for further clarification?

Thank you very much for your time and assistance!

Best，
Larry

[yujias424]

Hi @leeyoyohku,

Furthermore, another problem is about the following code snippet,

sel = np.arange(phenotypes.shape[1])
sel = list(sel[pd.Series(protein_names[:,0]).isin(db.loc[protein_name,['R1', 'R2']].dropna().values).values])
kin_from = phenotypes[:, sel] # remove the selected pr

If protein_name can retrieve multiple rows, it seems that the code won't work since the db.loc[protein_name,['R1', 'R2']].dropna().values will return a two dimension matrix?

Could you please provide some clarification or guidance on this?

Thank you again for your help!

Best regards,
Larry