Code for the simulated data

[DruncBread]

Hi, it's a great work!

I find your method for simulating partially overlapped slices very inspiring for my project, but I failed to find the code to generate simulated data.

I wonder if you can provide the code for the simulated data?

Best regards.

[x-h-liu]

Hi, here's the code:

def drop_spots_by_neighbors(slice_,eps=1e-1,power=1):
    dist = scipy.spatial.distance_matrix(slice_.obsm['spatial'],slice_.obsm['spatial'])
    neighbors = (((dist - dist[dist>0].min())/dist[dist>0].min())**2)<eps
    np.fill_diagonal(neighbors,0)
    neighbors = scipy.sparse.coo_matrix(neighbors)
    neighbors_ = neighbors.copy()
    for i in range(power-1): neighbors = neighbors.dot(neighbors_)
    neighbors = neighbors.toarray()
    prob = neighbors.sum(axis=0)/neighbors.sum(axis=0).max()
    select = np.array(np.random.binomial(1,prob),dtype=bool)
    temp_slice = slice_[slice_.obs.index[select]]
    return temp_slice

def rotate_coordinates(slice_,theta):
    R = np.array([[np.cos(theta),-np.sin(theta)],[np.sin(theta),np.cos(theta)]])
    temp_slice = slice_.copy()
    temp_slice.obsm['spatial'] = slice_.obsm['spatial'].dot(R)   # is multiplying R on the right of the coordinate doing the right thing? OK it might be clockwise rotation, and becomes counterclockwise when flipped along y axis.
    return temp_slice

def resample_counts(slice_, delta):
    temp_slice = slice_.copy()
    temp_mat = np.zeros(temp_slice.X.shape)
    for i in range(len(slice_)):
        v = temp_slice.X[i].copy().toarray().flatten()
        sumv = v.sum()
        if delta == 0:
            v = v.astype("float64")
            v = v / v.sum()
            v = np.random.multinomial(sumv, v)
            temp_mat[i] = v
        else:
            v = v + delta * 0.0002 * sumv
            v = v/v.sum()
            v = np.random.multinomial(sumv,v)
            temp_mat[i] = v
    temp_slice.X = scipy.sparse.csr.csr_matrix(temp_mat)
    sc.pp.calculate_qc_metrics(temp_slice, inplace=True)
    return temp_slice

def split_slice_overlap_single_resample(slice_, nrow, ncol, overlap, drop_edge_spots=False, power=1, rotate=False, resample=False, delta=0.1, seed=0):
    (min_x, min_y), (max_x, max_y) = slice_.obsm['spatial'].min(axis=0), slice_.obsm['spatial'].max(axis=0)
    len_x, len_y = max_x-min_x, max_y-min_y

    if ncol > 1:
        stride_x = len_x / ((overlap / (1 - overlap)) + ncol)
        overlap_length_x = (overlap / (1 - overlap)) * stride_x
        len_block_x = stride_x + overlap_length_x
    else:
        len_block_x, stride_x = len_x, 0
    if nrow > 1:
        stride_y = len_y / ((overlap / (1 - overlap)) + nrow)
        overlap_length_y = (overlap / (1 - overlap)) * stride_y
        len_block_y = stride_y + overlap_length_y
    else:
        len_block_y, stride_y = len_y, 0

    np.random.seed(seed)
    split_slices = {}
    for i in range(nrow):
        for j in range(ncol):
            if (i, j) == (0, 0) or (i, j) == (1, 0):
                idx = np.logical_and(np.logical_and(slice_.obsm['spatial'][:, 0] >= min_x + j * stride_x,
                                                    slice_.obsm['spatial'][:, 0] <= min_x + j * stride_x + len_block_x),
                                     np.logical_and(slice_.obsm['spatial'][:, 1] >= min_y + i * stride_y,
                                                    slice_.obsm['spatial'][:, 1] <= min_y + i * stride_y + len_block_y))
                temp_slice = slice_[slice_.obs.index[idx]]
                if drop_edge_spots:
                    temp_slice = drop_spots_by_neighbors(temp_slice, power=power)
                if rotate:
                    temp_slice = rotate_coordinates(temp_slice, theta=2*np.pi*(np.random.rand()-0.5))
                if resample and (i, j) == (1, 0):
                    temp_slice = resample_counts(temp_slice, delta=float(delta))
                plot_slice(temp_slice)
                split_slices[i, j] = temp_slice
    plt.show()
    return split_slices


delta_to_sim=[0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3.0, 3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 3.7, 3.8, 3.9, 4.0, 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9, 5.0]
overlap_to_sim=[0.9, 0.7, 0.5, 0.3]


for overlap in overlap_to_sim:
    for delta in delta_to_sim:
        print(overlap, delta)
        slice_name,nrow,ncol,drop_edge_spots,rotate,resample,seed = "151674",2,1,False,False,True,0
        split_slices = split_slice_overlap_single_resample(slices[slice_name], nrow, ncol, overlap, drop_edge_spots=drop_edge_spots, rotate=rotate, resample=resample, delta=delta, seed=seed)
        save=True
        if save:
            for i in range(nrow):
                for j in range(ncol):
                    if (i, j) == (0, 0) or (i, j) == (1, 0):
                        split_slices[i,j].write('SAVE_PATH/{0}_overlap{1}_drop{2}_rotate{3}_resample{4}_delta{5}_row{6}_col{7}.h5ad'.format(slice_name,overlap,drop_edge_spots, rotate,resample,delta,i,j))

slices is a dictionary mapping slice name (e,g, "151674") to the AnnData object of that slice.

Hope this is useful! Let me know if you have any questions, or if you need code to generate the AnnData object of each slice as well.

[DruncBread]

Thank you very much!