DEMO: Hepatic lobule

In the hepatic lobule dataset, we first integrate information from 10 Cosmx hepatic lobules to create landmark library and a spatial gene expression atlas of the lobules. We then introduce two Visium liver slices with different health states and use ARIEL to assist in lobule localization and analysis. By performing correlation analysis with the Cosmx atlas, we identify genes potentially associated with fatty liver disease.

%load_ext autoreload
%autoreload 2

import numpy as np
import math
import torch
import matplotlib.pyplot as plt
import scanpy as sc
import squidpy as sq
import anndata as ad
import pandas as pd
import anndata
from sklearn.cluster import KMeans
from scipy.spatial import ConvexHull, Delaunay
from ariel_srt.Landmark import normalization_spatial, Rasterization, alternative_landmark, screen_landmark
from ariel_srt.Alignment import rigid_alignment, non_rigid_alignment
from ariel_srt.Transfer import Appro_GP, Exact_GP
import random
random.seed(10)

c:\Anaconda3\envs\align\lib\site-packages\numba\core\decorators.py:246: RuntimeWarning: nopython is set for njit and is ignored
  warnings.warn('nopython is set for njit and is ignored', RuntimeWarning)

path = './data/hepatic lobule'
device = "cuda" if torch.cuda.is_available() else "cpu"
device

'cuda'

COSMX

for i in range(1,11):
    globals()[f"data_slice{i}"]=sc.read_h5ad(f'{path}/data_slice{i}.h5ad')

# Original information of ten Cosmx hepatic lobules

cl={}
sp={}
for i in range(1, 11):
    cl[f"data_slice{i}"]=globals()[f"data_slice{i}"].obsm['color']
    sp[f"data_slice{i}"]=globals()[f"data_slice{i}"].obsm['spatial']
num_plot = len(cl)

dx = []
dy = []
for j in range(10):
    dx.append(sp[f"data_slice{j+1}"][:,0].max() - sp[f"data_slice{j+1}"][:,0].min())
    dy.append(sp[f"data_slice{j+1}"][:,1].max() - sp[f"data_slice{j+1}"][:,1].min())
dx = max(dx)
dy = max(dy)

plt.figure(figsize=((50,20)))
for j in range(num_plot):
    plt.subplot(2,5,j+1)
    sp0 = sp[f'data_slice{j+1}']
    cl0 = cl[f'data_slice{j+1}']

    plt.scatter(sp0[:,0],sp0[:,1],s=5,c=cl0)
    if(j==0):
        plt.title('Slice1 (Reference)',fontsize=40)
    else:
        plt.title(f'Slice{j+1}',fontsize=40)
    plt.xticks(fontsize=30)
    plt.yticks(fontsize=30)

    dd = sp[f"data_slice{j+1}"][:,0].max() - sp[f"data_slice{j+1}"][:,0].min()
    plt.xlim(xmin=sp[f"data_slice{j+1}"][:,0].min() - (dx-dd)/2, xmax=sp[f"data_slice{j+1}"][:,0].max() + (dx-dd)/2)
    dd = sp[f"data_slice{j+1}"][:,1].max() - sp[f"data_slice{j+1}"][:,1].min()
    plt.ylim(ymin=sp[f"data_slice{j+1}"][:,1].min() - (dy-dd)/2, ymax=sp[f"data_slice{j+1}"][:,1].max() + (dy-dd)/2)

Select landmarks

spatial1 = normalization_spatial(data_slice1.obsm['spatial'])
data_slice1.obsm['spatial'] = spatial1

for i in range(2,11):
    globals()[f'lm1_{i}'], globals()[f'lm{i}'] = alternative_landmark(spatial1, globals()[f"data_slice{i}"].obsm['spatial'], data_slice1.obsm['pca'], globals()[f"data_slice{i}"].obsm['pca'], n = 100)

slice 2 to reference(slice 1)

In this experiment, we use landmarks to redefine the hepatic lobules, which imposes high requirements on the spatial distribution of the landmarks. Therefore, we further screen the alternative landmarks generated by ARIEL manually. Check interactive_lm_selection.py for detail.

# Here is a selection example. We use the interactive landmark selecting algorithm to find landmarks manually. 
from IPython.display import Image, display
display(Image(filename=f'{path}/select.gif'))

import matplotlib.image as mpimg
I = mpimg.imread(f'{path}/selected_lm.png')
plt.imshow(I)
plt.axis('off')

(-0.5, 900.5, 767.5, -0.5)

We demonstrate the process of manually selecting landmark pairs through animated graphs. The figure above shows the output results of the algorithm, and the numbers enclosed by red boxes are the selected landmark pairs in slice 2 - reference (slice 1) experiment.

Given the structural characteristics of liver lobules as polygonal prisms with symmetry, we sort the landmarks on two slices in a clockwise direction, then select the optimal matching set.

lm1_2=lm1_2[[38,3,48,49,30,24],:]
lm2=lm2[[24,3,49,30,48,38],:]
D=[]
for i in range(6):
    lm2_test = np.vstack((lm2[range(i,6),:],lm2[range(0,i),:])) 
    lm2_test = rigid_alignment(lm1_2, lm2_test)
    d=np.sum(np.sqrt(np.sum(np.power(lm2_test-lm1_2,2),axis=1)))
    D.append(d)
lm2 = np.vstack((lm2[range(np.argmin(D),6),:],lm2[range(0,np.argmin(D)),:])) 

plt.figure(figsize=((20,10)))
plt.subplot(1,2,1)
plt.scatter(spatial1[:,0],spatial1[:,1],s=1,c=data_slice1.obsm['color'])
plt.scatter(lm1_2[:,0],lm1_2[:,1],s = 1000,c=range(lm1_2.shape[0]),marker = "*",edgecolor = "black")
for i in range(lm1_2.shape[0]):
    plt.text(lm1_2[i,0],lm1_2[i,1],f'{i}', c='red',fontsize=30)

plt.subplot(1,2,2)
plt.scatter(data_slice2.obsm['spatial'][:,0],data_slice2.obsm['spatial'][:,1],s=1,c=data_slice2.obsm['color'])
plt.scatter(lm2[:,0],lm2[:,1],s = 1000,c=range(lm2.shape[0]),marker = "*",edgecolor = "black")
for i in range(lm2.shape[0]):
    plt.text(lm2[i,0],lm2[i,1],f'{i}', c='red',fontsize=30)

slice3 to reference(slice 1)

lm1_3=lm1_3[[13,12,5,22,1,46],:]
lm3=lm3[[12,1,22,5,46,13],:]
D=[]
for i in range(6):
    lm3_test = np.vstack((lm3[range(i,6),:],lm3[range(0,i),:])) 
    lm3_test = rigid_alignment(lm1_3, lm3_test)
    d=np.sum(np.sqrt(np.sum(np.power(lm3_test-lm1_3,2),axis=1)))
    D.append(d)
lm3=np.vstack((lm3[range(np.argmin(D),6),:],lm3[range(0,np.argmin(D)),:])) 

plt.figure(figsize=((20,10)))
plt.subplot(1,2,1)
plt.scatter(spatial1[:,0],spatial1[:,1],s=1,c=data_slice1.obsm['color'])
plt.scatter(lm1_3[:,0],lm1_3[:,1],s = 1000,c=range(lm1_3.shape[0]),marker = "*",edgecolor = "black")
for i in range(lm1_3.shape[0]):
    plt.text(lm1_3[i,0],lm1_3[i,1],f'{i}', c='red',fontsize=30)

plt.subplot(1,2,2)
plt.scatter(data_slice3.obsm['spatial'][:,0],data_slice3.obsm['spatial'][:,1],s=1,c=data_slice3.obsm['color'])
plt.scatter(lm3[:,0],lm3[:,1],s = 1000,c=range(lm3.shape[0]),marker = "*",edgecolor = "black")
for i in range(lm3.shape[0]):
    plt.text(lm3[i,0],lm3[i,1],f'{i}', c='red',fontsize=30)

slice4 to reference(slice 1)

lm1_4=lm1_4[[27,43,44,28,89,48],:]
lm4=lm4[[28,48,89,59,27,43],:]

D=[]
for i in range(6):
    lm4_test = np.vstack((lm4[range(i,6),:],lm4[range(0,i),:])) 
    lm4_test = rigid_alignment(lm1_4, lm4_test)
    d=np.sum(np.sqrt(np.sum(np.power(lm4_test-lm1_4,2),axis=1)))
    D.append(d)
lm4=np.vstack((lm4[range(np.argmin(D),6),:],lm4[range(0,np.argmin(D)),:])) 

plt.figure(figsize=((20,10)))
plt.subplot(1,2,1)
plt.scatter(spatial1[:,0],spatial1[:,1],s=1,c=data_slice1.obsm['color'])
plt.scatter(lm1_4[:,0],lm1_4[:,1],s = 1000,c=range(lm1_4.shape[0]),marker = "*",edgecolor = "black")
for i in range(lm1_4.shape[0]):
    plt.text(lm1_4[i,0],lm1_4[i,1],f'{i}', c='red',fontsize=30)

plt.subplot(1,2,2)
plt.scatter(data_slice4.obsm['spatial'][:,0],data_slice4.obsm['spatial'][:,1],s=1,c=data_slice4.obsm['color'])
plt.scatter(lm4[:,0],lm4[:,1],s = 1000,c=range(lm4.shape[0]),marker = "*",edgecolor = "black")
for i in range(lm4.shape[0]):
    plt.text(lm4[i,0],lm4[i,1],f'{i}', c='red',fontsize=30)

slice5 to reference(slice 1)

lm1_5=lm1_5[[33,16,26,45,30,52],:]
lm5=lm5[[33,52,30,16,45,26],:]

D=[]
for i in range(6):
    lm5_test = np.vstack((lm5[range(i,6),:],lm5[range(0,i),:])) 
    lm5_test = rigid_alignment(lm1_5, lm5_test)
    d=np.sum(np.sqrt(np.sum(np.power(lm5_test-lm1_5,2),axis=1)))
    D.append(d)
lm5=np.vstack((lm5[range(np.argmin(D),6),:],lm5[range(0,np.argmin(D)),:])) 

plt.figure(figsize=((20,10)))
plt.subplot(1,2,1)
plt.scatter(spatial1[:,0],spatial1[:,1],s=1,c=data_slice1.obsm['color'])
plt.scatter(lm1_5[:,0],lm1_5[:,1],s = 1000,c=range(lm1_5.shape[0]),marker = "*",edgecolor = "black")
for i in range(lm1_5.shape[0]):
    plt.text(lm1_5[i,0],lm1_5[i,1],f'{i}', c='red',fontsize=30)

plt.subplot(1,2,2)
plt.scatter(data_slice5.obsm['spatial'][:,0],data_slice5.obsm['spatial'][:,1],s=1,c=data_slice5.obsm['color'])
plt.scatter(lm5[:,0],lm5[:,1],s = 1000,c=range(lm5.shape[0]),marker = "*",edgecolor = "black")
for i in range(lm5.shape[0]):
    plt.text(lm5[i,0],lm5[i,1],f'{i}', c='red',fontsize=30)

slice6 to reference(slice 1)

lm1_6=lm1_6[[53,5,23,60,9,65],:]
lm6=lm6[[5,60,53,65,9,23],:]

D=[]
for i in range(6):
    lm6_test = np.vstack((lm6[range(i,6),:],lm6[range(0,i),:])) 
    lm6_test = rigid_alignment(lm1_6, lm6_test)
    d=np.sum(np.sqrt(np.sum(np.power(lm6_test-lm1_6,2),axis=1)))
    D.append(d)
lm6 = np.vstack((lm6[range(np.argmin(D),6),:],lm6[range(0,np.argmin(D)),:])).copy() 

plt.figure(figsize=((20,10)))
plt.subplot(1,2,1)
plt.scatter(spatial1[:,0],spatial1[:,1],s=1,c=data_slice1.obsm['color'])
plt.scatter(lm1_6[:,0],lm1_6[:,1],s = 1000,c=range(lm1_6.shape[0]),marker = "*",edgecolor = "black")
for i in range(lm1_6.shape[0]):
    plt.text(lm1_6[i,0],lm1_6[i,1],f'{i}', c='red',fontsize=30)

plt.subplot(1,2,2)
plt.scatter(data_slice6.obsm['spatial'][:,0],data_slice6.obsm['spatial'][:,1],s=1,c=data_slice6.obsm['color'])
plt.scatter(lm6[:,0],lm6[:,1],s = 1000,c=range(lm6.shape[0]),marker = "*",edgecolor = "black")
for i in range(lm6.shape[0]):
    plt.text(lm6[i,0],lm6[i,1],f'{i}', c='red',fontsize=30)

slice7 to reference(slice 1)

lm1_7=lm1_7[[5,37,18,10,7,0],:]
lm7=lm7[[10,18,5,0,37,7],:]

D=[]
for i in range(6):
    lm7_test = np.vstack((lm7[range(i,6),:],lm7[range(0,i),:])) 
    lm7_test = rigid_alignment(lm1_7, lm7_test)
    d=np.sum(np.sqrt(np.sum(np.power(lm7_test-lm1_7,2),axis=1)))
    D.append(d)
lm7=np.vstack((lm7[range(np.argmin(D),6),:],lm7[range(0,np.argmin(D)),:])) 

plt.figure(figsize=((20,10)))
plt.subplot(1,2,1)
plt.scatter(spatial1[:,0],spatial1[:,1],s=1,c=data_slice1.obsm['color'])
plt.scatter(lm1_7[:,0],lm1_7[:,1],s = 1000,c=range(lm1_7.shape[0]),marker = "*",edgecolor = "black")
for i in range(lm1_7.shape[0]):
    plt.text(lm1_7[i,0],lm1_7[i,1],f'{i}', c='red',fontsize=30)

plt.subplot(1,2,2)
plt.scatter(data_slice7.obsm['spatial'][:,0],data_slice7.obsm['spatial'][:,1],s=1,c=data_slice7.obsm['color'])
plt.scatter(lm7[:,0],lm7[:,1],s = 1000,c=range(lm7.shape[0]),marker = "*",edgecolor = "black")
for i in range(lm7.shape[0]):
    plt.text(lm7[i,0],lm7[i,1],f'{i}', c='red',fontsize=30)

slice8 to reference(slice 1)

lm1_8=lm1_8[[23,0,45,25,6,15],:]
lm8=lm8[[25,6,15,23,45,0],:]

D=[]
for i in range(6):
    lm8_test = np.vstack((lm8[range(i,6),:],lm8[range(0,i),:])) 
    lm8_test = rigid_alignment(lm1_8, lm8_test)
    d=np.sum(np.sqrt(np.sum(np.power(lm8_test-lm1_8,2),axis=1)))
    D.append(d)
lm8=np.vstack((lm8[range(np.argmin(D),6),:],lm8[range(0,np.argmin(D)),:])) 

plt.figure(figsize=((20,10)))
plt.subplot(1,2,1)
plt.scatter(spatial1[:,0],spatial1[:,1],s=1,c=data_slice1.obsm['color'])
plt.scatter(lm1_8[:,0],lm1_8[:,1],s = 1000,c=range(lm1_8.shape[0]),marker = "*",edgecolor = "black")
for i in range(lm1_8.shape[0]):
    plt.text(lm1_8[i,0],lm1_8[i,1],f'{i}', c='red',fontsize=30)

plt.subplot(1,2,2)
plt.scatter(data_slice8.obsm['spatial'][:,0],data_slice8.obsm['spatial'][:,1],s=1,c=data_slice8.obsm['color'])
plt.scatter(lm8[:,0],lm8[:,1],s = 1000,c=range(lm8.shape[0]),marker = "*",edgecolor = "black")
for i in range(lm8.shape[0]):
    plt.text(lm8[i,0],lm8[i,1],f'{i}', c='red',fontsize=30)

slice9 to reference(slice 1)

lm1_9=lm1_9[[29,9,45,58,32,12],:]
lm9=lm9[[9,29,58,45,12,32],:]

D=[]
for i in range(6):
    lm9_test = np.vstack((lm9[range(i,6),:],lm9[range(0,i),:])) 
    lm9_test = rigid_alignment(lm1_9, lm9_test)
    d=np.sum(np.sqrt(np.sum(np.power(lm9_test-lm1_9,2),axis=1)))
    D.append(d)
lm9=np.vstack((lm9[range(np.argmin(D),6),:],lm9[range(0,np.argmin(D)),:])) 

plt.figure(figsize=((20,10)))
plt.subplot(1,2,1)
plt.scatter(spatial1[:,0],spatial1[:,1],s=1,c=data_slice1.obsm['color'])
plt.scatter(lm1_9[:,0],lm1_9[:,1],s = 1000,c=range(lm1_9.shape[0]),marker = "*",edgecolor = "black")
for i in range(lm1_9.shape[0]):
    plt.text(lm1_9[i,0],lm1_9[i,1],f'{i}', c='red',fontsize=30)

plt.subplot(1,2,2)
plt.scatter(data_slice9.obsm['spatial'][:,0],data_slice9.obsm['spatial'][:,1],s=1,c=data_slice9.obsm['color'])
plt.scatter(lm9[:,0],lm9[:,1],s = 1000,c=range(lm9.shape[0]),marker = "*",edgecolor = "black")
for i in range(lm9.shape[0]):
    plt.text(lm9[i,0],lm9[i,1],f'{i}', c='red',fontsize=30)

slice10 to reference(slice 1)

lm1_10=lm1_10[[19,44,0,15,47,26],:]
lm10=lm10[[44,15,0,47,19,26],:]

D=[]
for i in range(6):
    lm10_test = np.vstack((lm10[range(i,6),:],lm10[range(0,i),:])) 
    lm10_test = rigid_alignment(lm1_10, lm10_test)
    d=np.sum(np.sqrt(np.sum(np.power(lm10_test-lm1_10,2),axis=1)))
    D.append(d)
lm10=np.vstack((lm10[range(np.argmin(D),6),:],lm10[range(0,np.argmin(D)),:])) 

plt.figure(figsize=((20,10)))
plt.subplot(1,2,1)
plt.scatter(spatial1[:,0],spatial1[:,1],s=1,c=data_slice1.obsm['color'])
plt.scatter(lm1_10[:,0],lm1_10[:,1],s = 1000,c=range(lm1_10.shape[0]),marker = "*",edgecolor = "black")
for i in range(lm1_10.shape[0]):
    plt.text(lm1_10[i,0],lm1_10[i,1],f'{i}', c='red',fontsize=30)

plt.subplot(1,2,2)
plt.scatter(data_slice10.obsm['spatial'][:,0],data_slice10.obsm['spatial'][:,1],s=1,c=data_slice10.obsm['color'])
plt.scatter(lm10[:,0],lm10[:,1],s = 1000,c=range(lm10.shape[0]),marker = "*",edgecolor = "black")
for i in range(lm10.shape[0]):
    plt.text(lm10[i,0],lm10[i,1],f'{i}', c='red',fontsize=30)

Landmark results

dx = []
dy = []
for j in range(10):
    dx.append(globals()[f"data_slice{j+1}"].obsm['spatial'][:,0].max() - globals()[f"data_slice{j+1}"].obsm['spatial'][:,0].min())
    dy.append(globals()[f"data_slice{j+1}"].obsm['spatial'][:,1].max() - globals()[f"data_slice{j+1}"].obsm['spatial'][:,1].min())
dx = max(dx)
dy = max(dy)

plt.figure(figsize=((50,20)))
for j in range(10):
    plt.subplot(2,5,j+1)

    plt.scatter(globals()[f'data_slice{j+1}'].obsm['spatial'][:,0],globals()[f'data_slice{j+1}'].obsm['spatial'][:,1],s=5,c=globals()[f'data_slice{j+1}'].obsm['color'])
    plt.xticks(fontsize=30)
    plt.yticks(fontsize=30)
    if(j==0):
        plt.title('Reference',fontsize=40)
        for i in range(2,11):
            plt.scatter(globals()[f"lm1_{i}"][:,0],globals()[f"lm1_{i}"][:,1],s=500,c=plt.cm.tab10(i-1))
    else:
        plt.title(f'Slice{j+1}',fontsize=40)
        plt.scatter(globals()[f"lm{j+1}"][:,0],globals()[f"lm{j+1}"][:,1],s=500,c=plt.cm.tab10(j))
    
        dd = globals()[f"data_slice{j+1}"].obsm['spatial'][:,0].max() - globals()[f"data_slice{j+1}"].obsm['spatial'][:,0].min()
        plt.xlim(xmin=globals()[f"data_slice{j+1}"].obsm['spatial'][:,0].min() - (dx-dd)/2, xmax=globals()[f"data_slice{j+1}"].obsm['spatial'][:,0].max() + (dx-dd)/2)
        dd = globals()[f"data_slice{j+1}"].obsm['spatial'][:,1].max() - globals()[f"data_slice{j+1}"].obsm['spatial'][:,1].min()
        plt.ylim(ymin=globals()[f"data_slice{j+1}"].obsm['spatial'][:,1].min() - (dy-dd)/2, ymax=globals()[f"data_slice{j+1}"].obsm['spatial'][:,1].max() + (dy-dd)/2)
     

C:\Users\Hanson Yin\AppData\Local\Temp\ipykernel_123848\457440794.py:11: UserWarning: *c* argument looks like a single numeric RGB or RGBA sequence, which should be avoided as value-mapping will have precedence in case its length matches with *x* & *y*.  Please use the *color* keyword-argument or provide a 2D array with a single row if you intend to specify the same RGB or RGBA value for all points.
  plt.scatter(globals()[f"lm1_{i}"][:,0],globals()[f"lm1_{i}"][:,1],s=500,c=plt.cm.tab10(i-1))
C:\Users\Hanson Yin\AppData\Local\Temp\ipykernel_123848\457440794.py:14: UserWarning: *c* argument looks like a single numeric RGB or RGBA sequence, which should be avoided as value-mapping will have precedence in case its length matches with *x* & *y*.  Please use the *color* keyword-argument or provide a 2D array with a single row if you intend to specify the same RGB or RGBA value for all points.
  plt.scatter(globals()[f"lm{j+1}"][:,0],globals()[f"lm{j+1}"][:,1],s=500,c=plt.cm.tab10(j))

Landmark library

We collect the expression information from all landmarks and their surrounding points on the reference, and output the result as a landmark library.

lm1 = np.empty((0,2))
for i in range(2,11):
    lm1 = np.vstack((lm1, globals()[f'lm1_{i}']))

spatial_new1, gene1=Rasterization(spatial1, data_slice1.X,nx=50,ny=50)
index_lmlibrary=[]
for i in range(lm1.shape[0]):
    k=np.argmin(np.power((spatial_new1[:,0]-lm1[i,0]),2)+np.power((spatial_new1[:,1]-lm1[i,1]),2))
    if not k in index_lmlibrary:
        index_lmlibrary.append(k)

lm_library = anndata.AnnData(X=gene1)
lm_library.obsm['spatial'] = spatial_new1
lm_library.var_names=data_slice1.var_names
lm_library.uns['index_lm'] = index_lmlibrary
print(lm_library)

AnnData object with n_obs × n_vars = 1992 × 1000
    uns: 'index_lm'
    obsm: 'spatial'

plt.figure(figsize=((10,10)))

plt.scatter(spatial1[:,0],spatial1[:,1],s=5)
plt.scatter(lm_library.obsm['spatial'][lm_library.uns['index_lm'],0],lm_library.obsm['spatial'][lm_library.uns['index_lm'],1],s=300)
plt.axis("off")

(-0.05, 1.05, -0.05, 1.05)

lm_library.write_h5ad(f'{path}/lm_library.h5ad')
gene = pd.DataFrame(lm_library.X)
gene.columns=data_slice1.var_names
gene.to_csv(f'{path}/lm_library_gene.csv', index=False)
pd.DataFrame(lm_library.obsm['spatial']).to_csv(f'{path}/lm_library_spatial.csv',header=False, index=False)

Redefine the hepatic lobule

Based on the landmarks selected in each slice, we redefine the lobule boundaries.

lm = lm_library.obsm['spatial'][lm_library.uns['index_lm'],:]
hull=ConvexHull(lm)
delaunay = Delaunay(lm[hull.vertices])
result = delaunay.find_simplex(spatial1)>=0
data_slice1 = data_slice1[result,:]
spatial1 = spatial1[result,:]

for i in range(2,11):
    hull=ConvexHull(globals()[f"lm{i}"])
    delaunay = Delaunay(globals()[f"lm{i}"][hull.vertices])
    result = delaunay.find_simplex(globals()[f'data_slice{i}'].obsm['spatial'])>=0
    globals()[f'data_slice{i}'] = globals()[f'data_slice{i}'][result,:]
    globals()[f'spatial{i}'] = globals()[f'data_slice{i}'].obsm['spatial']

plt.figure(figsize=((50,20)))
for j in range(10):
    plt.subplot(2,5,j+1)
    plt.scatter(globals()[f"data_slice{j+1}"].obsm['spatial'][:,0],globals()[f"data_slice{j+1}"].obsm['spatial'][:,1],s=5,c=globals()[f"data_slice{j+1}"].obsm['color'])
    plt.xticks(fontsize=30)
    plt.yticks(fontsize=30)
    if(j==0):
        plt.title('Reference',fontsize=40)
        for i in range(2,11):
            plt.scatter(globals()[f"lm1_{i}"][:,0],globals()[f"lm1_{i}"][:,1],s=500,c=plt.cm.tab10(i-1))
    else:
        plt.title(f'Slice{j+1}',fontsize=40)
        plt.scatter(globals()[f"lm{j+1}"][:,0],globals()[f"lm{j+1}"][:,1],s=500,c=plt.cm.tab10(j))
    
        dd = globals()[f"data_slice{j+1}"].obsm['spatial'][:,0].max() - globals()[f"data_slice{j+1}"].obsm['spatial'][:,0].min()
        plt.xlim(xmin=globals()[f"data_slice{j+1}"].obsm['spatial'][:,0].min() - (dx-dd)/2, xmax=globals()[f"data_slice{j+1}"].obsm['spatial'][:,0].max() + (dx-dd)/2)
        dd = globals()[f"data_slice{j+1}"].obsm['spatial'][:,1].max() - globals()[f"data_slice{j+1}"].obsm['spatial'][:,1].min()
        plt.ylim(ymin=globals()[f"data_slice{j+1}"].obsm['spatial'][:,1].min() - (dy-dd)/2, ymax=globals()[f"data_slice{j+1}"].obsm['spatial'][:,1].max() + (dy-dd)/2)

plt.tight_layout()

C:\Users\Hanson Yin\AppData\Local\Temp\ipykernel_123848\3322210828.py:10: UserWarning: *c* argument looks like a single numeric RGB or RGBA sequence, which should be avoided as value-mapping will have precedence in case its length matches with *x* & *y*.  Please use the *color* keyword-argument or provide a 2D array with a single row if you intend to specify the same RGB or RGBA value for all points.
  plt.scatter(globals()[f"lm1_{i}"][:,0],globals()[f"lm1_{i}"][:,1],s=500,c=plt.cm.tab10(i-1))
C:\Users\Hanson Yin\AppData\Local\Temp\ipykernel_123848\3322210828.py:13: UserWarning: *c* argument looks like a single numeric RGB or RGBA sequence, which should be avoided as value-mapping will have precedence in case its length matches with *x* & *y*.  Please use the *color* keyword-argument or provide a 2D array with a single row if you intend to specify the same RGB or RGBA value for all points.
  plt.scatter(globals()[f"lm{j+1}"][:,0],globals()[f"lm{j+1}"][:,1],s=500,c=plt.cm.tab10(j))

for i in range(1,11):
    globals()[f"data_slice{i}"].write_h5ad(f'{path}/data_slice{i}_redefine.h5ad')

Alignment

Rigid alignment

for i in range(2,11):
    globals()[f'lm{i}'], globals()[f'spatial{i}'] = rigid_alignment(globals()[f'lm1_{i}'], globals()[f'lm{i}'], globals()[f'spatial{i}'])

plt.figure(figsize=((50,20)))

for i in range(1,11):
    plt.subplot(2,5, i)
    
    plt.scatter(globals()[f"spatial{i}"][:,0],
                globals()[f"spatial{i}"][:,1],
                c=globals()[f'data_slice{i}'].obsm['color'],
                s=5)
    if(i==1):
        plt.title('Slice1 (Reference)',fontsize=40)
    else:
        plt.title(f'Slice{i}',fontsize=40)
        
    plt.xticks(fontsize=30)
    plt.yticks(fontsize=30)
    plt.xlim(xmin=0, xmax=1)
    plt.ylim(ymin=0, ymax=1)
    
plt.tight_layout()

Non-rigid alignment

for i in range(2,11):
    globals()[f'lm1_{i}'] = lm1[range(6*(i-2),6*(i-1)),:]
    globals()[f'lm{i}'], globals()[f'spatial{i}'] = non_rigid_alignment(globals()[f'lm1_{i}'], globals()[f'lm{i}'], globals()[f'spatial{i}'])

plt.figure(figsize=((50,20)))

for i in range(1,11):
    plt.subplot(2,5, i)
    
    plt.scatter(globals()[f"spatial{i}"][:,0],
                globals()[f"spatial{i}"][:,1],
                c=globals()[f'data_slice{i}'].obsm['color'],
                s=5)
    if(i==1):
        plt.title('Slice1 (Reference)',fontsize=40)
    else:
        plt.title(f'Slice{i}',fontsize=40)
    plt.xlim(xmin=0, xmax=1)
    plt.ylim(ymin=0, ymax=1)
    plt.xticks(fontsize=30)
    plt.yticks(fontsize=30)

plt.tight_layout()

plt.figure(figsize=((10,10)))
for i in range(1,11):
    plt.scatter(globals()[f"spatial{i}"][:,0],
                globals()[f"spatial{i}"][:,1],
                c=globals()[f'data_slice{i}'].obsm['color'],
                s=5)

plt.xticks(fontsize=30)
plt.yticks(fontsize=30)

(array([0. , 0.2, 0.4, 0.6, 0.8, 1. ]),
 [Text(0, 0.0, '0.0'),
  Text(0, 0.2, '0.2'),
  Text(0, 0.4, '0.4'),
  Text(0, 0.6000000000000001, '0.6'),
  Text(0, 0.8, '0.8'),
  Text(0, 1.0, '1.0')])

Information transfer

for i in range(1,11):
    x1 = torch.tensor(spatial1).float()
    x2 = torch.tensor(globals()[f'spatial{i}']).float()
    y2 = torch.tensor(globals()[f'data_slice{i}'].obsm['pca']).float()
    globals()[f'predictions{i}_pca'] = Appro_GP(x1, x2, y2, device = device, n_intro = 100, num_epochs = 500, lr = 0.01)
    if(device == 'cpu'): 
        globals()[f'predictions{i}_pca'] = globals()[f'predictions{i}_pca'].detach().numpy()
    else:
        globals()[f'predictions{i}_pca'] = globals()[f'predictions{i}_pca'].detach().cpu().numpy()
    globals()[f'y_predict{i}_pca_all'] = np.dot(globals()[f'predictions{i}_pca'],globals()[f'data_slice{i}'].varm['PCs'].T)
    pd.DataFrame(globals()[f'y_predict{i}_pca_all']).to_csv(f'{path}/data_slice{i}_pca_all_cut.csv',index=False)

c:\Anaconda3\envs\align\lib\site-packages\sklearn\cluster\_kmeans.py:1416: FutureWarning: The default value of `n_init` will change from 10 to 'auto' in 1.4. Set the value of `n_init` explicitly to suppress the warning
  super()._check_params_vs_input(X, default_n_init=10)
c:\Anaconda3\envs\align\lib\site-packages\sklearn\cluster\_kmeans.py:1416: FutureWarning: The default value of `n_init` will change from 10 to 'auto' in 1.4. Set the value of `n_init` explicitly to suppress the warning
  super()._check_params_vs_input(X, default_n_init=10)
c:\Anaconda3\envs\align\lib\site-packages\sklearn\cluster\_kmeans.py:1416: FutureWarning: The default value of `n_init` will change from 10 to 'auto' in 1.4. Set the value of `n_init` explicitly to suppress the warning
  super()._check_params_vs_input(X, default_n_init=10)
c:\Anaconda3\envs\align\lib\site-packages\sklearn\cluster\_kmeans.py:1416: FutureWarning: The default value of `n_init` will change from 10 to 'auto' in 1.4. Set the value of `n_init` explicitly to suppress the warning
  super()._check_params_vs_input(X, default_n_init=10)
c:\Anaconda3\envs\align\lib\site-packages\sklearn\cluster\_kmeans.py:1416: FutureWarning: The default value of `n_init` will change from 10 to 'auto' in 1.4. Set the value of `n_init` explicitly to suppress the warning
  super()._check_params_vs_input(X, default_n_init=10)
c:\Anaconda3\envs\align\lib\site-packages\sklearn\cluster\_kmeans.py:1416: FutureWarning: The default value of `n_init` will change from 10 to 'auto' in 1.4. Set the value of `n_init` explicitly to suppress the warning
  super()._check_params_vs_input(X, default_n_init=10)
c:\Anaconda3\envs\align\lib\site-packages\sklearn\cluster\_kmeans.py:1416: FutureWarning: The default value of `n_init` will change from 10 to 'auto' in 1.4. Set the value of `n_init` explicitly to suppress the warning
  super()._check_params_vs_input(X, default_n_init=10)
c:\Anaconda3\envs\align\lib\site-packages\sklearn\cluster\_kmeans.py:1416: FutureWarning: The default value of `n_init` will change from 10 to 'auto' in 1.4. Set the value of `n_init` explicitly to suppress the warning
  super()._check_params_vs_input(X, default_n_init=10)
c:\Anaconda3\envs\align\lib\site-packages\sklearn\cluster\_kmeans.py:1416: FutureWarning: The default value of `n_init` will change from 10 to 'auto' in 1.4. Set the value of `n_init` explicitly to suppress the warning
  super()._check_params_vs_input(X, default_n_init=10)
c:\Anaconda3\envs\align\lib\site-packages\sklearn\cluster\_kmeans.py:1416: FutureWarning: The default value of `n_init` will change from 10 to 'auto' in 1.4. Set the value of `n_init` explicitly to suppress the warning
  super()._check_params_vs_input(X, default_n_init=10)

for i in range(1,11):
    globals()[f'y_predict{i}_pca_all']=np.array(pd.read_csv(f'{path}/data_slice{i}_pca_all_cut.csv',header=0))
    globals()[f'n_x{i}']=globals()[f'y_predict{i}_pca_all'].copy()
    for j in range(globals()[f'n_x{i}'].shape[1]):
        globals()[f'n_x{i}'][:,j]=globals()[f'n_x{i}'][:,j]-min(globals()[f'n_x{i}'][:,j])

# Original spatial expression of GLUL
plt.figure(figsize=((50,15)))
for i in range(1,11):
    plt.subplot(2,5, i)
    
    plt.scatter(globals()[f'data_slice{i}'].obsm['spatial'][:,0],
                globals()[f'data_slice{i}'].obsm['spatial'][:,1],
                c=globals()[f'data_slice{i}'].X[:,np.where(data_slice1.var_names=='GLUL')[0][0]],
                cmap='PuBu',
                s=5)
    plt.colorbar().ax.tick_params(labelsize=30)
    plt.title(f'Slice{i}',fontsize=40)
        
    plt.xticks(fontsize=30)
    plt.yticks(fontsize=30)

plt.tight_layout()

# spatial expression of GLUL after transfer
plt.figure(figsize=((50,15)))
for i in range(1,11):
    plt.subplot(2,5, i)
    
    plt.scatter(spatial1[:,0],
                spatial1[:,1],
                c=globals()[f'n_x{i}'][:,np.where(data_slice1.var_names=='GLUL')[0][0]],
                cmap='PuBu',
                s=5)
    plt.colorbar().ax.tick_params(labelsize=30)
    plt.title(f'Slice{i}',fontsize=40)
        
    plt.xticks(fontsize=30)
    plt.yticks(fontsize=30)

plt.tight_layout()

# Original spatial expression of SERPINA1
plt.figure(figsize=((50,15)))
for i in range(1,11):
    plt.subplot(2,5, i)
    
    plt.scatter(globals()[f'data_slice{i}'].obsm['spatial'][:,0],
                globals()[f'data_slice{i}'].obsm['spatial'][:,1],
                c=globals()[f'data_slice{i}'].X[:,np.where(data_slice1.var_names=='SERPINA1')[0][0]],
                cmap='PuBu',
                s=5)
    plt.colorbar().ax.tick_params(labelsize=30)
    plt.title(f'Slice{i}',fontsize=40)
        
    plt.xticks(fontsize=30)
    plt.yticks(fontsize=30)

plt.tight_layout()

# spatial expression of SERPINA1 after transfer
plt.figure(figsize=((50,15)))
for i in range(1,11):
    plt.subplot(2,5, i)
    
    plt.scatter(spatial1[:,0],
                spatial1[:,1],
                c=globals()[f'n_x{i}'][:,np.where(data_slice1.var_names=='SERPINA1')[0][0]],
                cmap='PuBu',
                s=5)
    plt.colorbar().ax.tick_params(labelsize=30)
    plt.title(f'Slice{i}',fontsize=40)
        
    plt.xticks(fontsize=30)
    plt.yticks(fontsize=30)

plt.tight_layout()

Spatial composite expression profiles

for i in range(1,11):
    globals()[f'n_x{i}_to01']=globals()[f'n_x{i}']
    for j in range(globals()[f'n_x{i}_to01'].shape[1]):
        globals()[f'n_x{i}_to01'][:,j]=(globals()[f'n_x{i}'][:,j]-np.min(globals()[f'n_x{i}'][:,j]))/(np.max(globals()[f'n_x{i}'][:,j])-np.min(globals()[f'n_x{i}'][:,j]))

composite=np.zeros(globals()[f'n_x{i}_to01'].shape)
for i in range(1,11):
    composite=composite+globals()[f'n_x{i}_to01']

composite=composite/10
composite=anndata.AnnData(X=composite)
composite.obsm['spatial']=data_slice1.obsm['spatial']
composite.var_names=data_slice1.var_names
composite.write_h5ad(f'{path}/composite.h5ad')

plt.figure(figsize=((10,9)))
scatter1 = plt.scatter(spatial1[:,0],spatial1[:,1],c=composite.X[:,np.where(data_slice1.var_names=='GLUL')],s=5,cmap='PuBu')
cbar = plt.colorbar(scatter1)
plt.axis("off")
plt.title('GLUL',fontsize=60, style='italic')
cbar.ax.tick_params(labelsize=40)

plt.figure(figsize=((10,9)))
scatter1 = plt.scatter(spatial1[:,0],spatial1[:,1],c=composite.X[:,np.where(data_slice1.var_names=='SERPINA1')],s=5,cmap='PuBu')
cbar = plt.colorbar(scatter1)
plt.axis("off")
plt.title('SERPINA1',fontsize=60, style='italic')
cbar.ax.tick_params(labelsize=40)

Two new Visium liver (different health states)

fat2 = sc.read_h5ad(f'{path}/fat2.h5ad')
health1 = sc.read_h5ad(f'{path}/health1.h5ad')
composite = sc.read_h5ad(f'{path}/composite.h5ad')
lm_library = sc.read_h5ad(f'{path}/lm_library.h5ad')

embedding_fat2_ll=np.array(pd.read_csv(f'{path}/ll_fat2_ll.csv')).T
embedding_fat2_ll=embedding_fat2_ll[lm_library.uns['index_lm'],:]
embedding_fat2=np.array(pd.read_csv(f'{path}/ll_fat2_fat2.csv')).T

embedding_health1_ll=np.array(pd.read_csv(f'{path}/ll_health1_ll.csv')).T
embedding_health1_ll=embedding_health1_ll[lm_library.uns['index_lm'],:]
embedding_health1=np.array(pd.read_csv(f'{path}/ll_health1_health1.csv')).T

lm_library_spatial=lm_library.obsm['spatial']
lm_library_spatial=lm_library_spatial[lm_library.uns['index_lm'],:]

Select landmarks and find hepatic lobule in Visium slice

lm_fat2_reference, lm_fat2 = alternative_landmark(lm_library_spatial, fat2.obsm['spatial'], embedding_fat2_ll, embedding_fat2, n = 100, replace1 = True)
lm_fat2_plot = lm_fat2.copy()

lm_health1_reference, lm_health1 = alternative_landmark(lm_library_spatial, health1.obsm['spatial'], embedding_health1_ll, embedding_health1, n = 600, replace1 = True)
lm_health1_plot = lm_health1.copy()

# Manually select the landmarks
lm_fat2=lm_fat2[[77,94,50,46,85,53,97,60,3,5,95,23,21,0,62,48,47,52,33,58,25,20,86,69,34,24,40,87,38,98,67,6,84,43],:][[23,17,8,4,0,24],:]
lm_fat2_reference=lm_fat2_reference[[77,94,50,46,85,53,97,60,3,5,95,23,21,0,62,48,47,52,33,58,25,20,86,69,34,24,40,87,38,98,67,6,84,43],:][[23,0,17,8,4,24],:]
D=[]
for i in range(6):
    lm_fat2_test=np.vstack((lm_fat2[range(i,6),:],lm_fat2[range(0,i),:])) 
    lm_fat2_test=rigid_alignment(lm_fat2_reference, lm_fat2_test)
    d=np.sum(np.sqrt(np.sum(np.power(lm_fat2_test-lm_fat2_reference,2),axis=1)))
    D.append(d)
lm_fat2=np.vstack((lm_fat2[range(np.argmin(D),6),:],lm_fat2[range(0,np.argmin(D)),:])) 

lm_health1=lm_health1[[343,545,532,527,340,510,139,158,161,540,587,309,66,549,557,227,547,46,164,576,74,210,292,574,303,4,186,591,118,530,272,529],:][[1,30,26,22,15,5],:]
lm_health1_reference=lm_health1_reference[[343,545,532,527,340,510,139,158,161,540,587,309,66,549,557,227,547,46,164,576,74,210,292,574,303,4,186,591,118,530,272,529],:][[22,30,15,26,1,5],:]
D=[]
for i in range(6):
    lm_health1_test=np.vstack((lm_health1[range(i,6),:],lm_health1[range(0,i),:])) 
    lm_health1_test=rigid_alignment(lm_health1_reference, lm_health1_test)
    d=np.sum(np.sqrt(np.sum(np.power(lm_health1_test-lm_health1_reference,2),axis=1)))
    D.append(d)
lm_health1=np.vstack((lm_health1[range(np.argmin(D),6),:],lm_health1[range(0,np.argmin(D)),:])) 

# Defining liver lobules on Visium slices using convex sets formed by landmarks
hull=ConvexHull(lm_fat2)
delaunay = Delaunay(lm_fat2[hull.vertices])
result_fat2 = delaunay.find_simplex(fat2.obsm['spatial'])>=0
line_fat2 = np.concatenate((lm_fat2[hull.vertices],lm_fat2[hull.vertices][0,:].reshape((1,2))),axis=0)

hull=ConvexHull(lm_health1)
delaunay = Delaunay(lm_health1[hull.vertices])
result_health1 = delaunay.find_simplex(health1.obsm['spatial'])>=0
line_health1 = np.concatenate((lm_health1[hull.vertices],lm_health1[hull.vertices][0,:].reshape((1,2))),axis=0)

Visualize the landmark selection result

plt.figure(figsize=(20,10))
plt.subplot(1,2,1)
color_map = {'Mid': 'yellow', 'Central': 'red', 'Portal': 'blue','Periportal':'green','Unknown':'black'}
group=pd.DataFrame(health1.obsm['zonationGroup'])
group.columns=['zonationGroup']
c = group['zonationGroup'].map(color_map)
plt.scatter(health1.obsm['spatial'][:,0],health1.obsm['spatial'][:,1],c=c,s=70)
plt.scatter(lm_health1[:,0],lm_health1[:,1],s=800,marker = "*",edgecolor = "black",c='orange')
plt.plot(line_health1[:,0],line_health1[:,1],c='black', linewidth=3)

plt.axis('equal')
plt.axis("off")

plt.subplot(1,2,2)
color_map = {'Mid': 'yellow', 'Central': 'red', 'Portal': 'blue','Periportal':'green','Unknown':'black'}
group=pd.DataFrame(fat2.obsm['zonationGroup'])
group.columns=['zonationGroup']
c = group['zonationGroup'].map(color_map)
plt.scatter(fat2.obsm['spatial'][:,0],fat2.obsm['spatial'][:,1],c=c,s=70)
plt.scatter(lm_fat2[:,0],lm_fat2[:,1],s=800,marker = "*",edgecolor = "black",c='orange')
plt.plot(line_fat2[:,0],line_fat2[:,1],c='black', linewidth=3)

import matplotlib. patches as mpatches
patch1 = mpatches.Patch (color='yellow', label='Mid')
patch2 = mpatches.Patch (color='green', label='Periportal')   
patch3 = mpatches.Patch (color='blue', label='Portal')
patch4 = mpatches.Patch (color='red', label='Central') 
patch5 = mpatches.Patch (color='black', label='Unknown')
handles, labels = plt. gca (). get_legend_handles_labels ()
handles. extend ([patch1, patch2, patch3, patch4, patch5])
plt.legend (handles=handles,loc="lower left", fontsize=35)

plt.axis('equal')
plt.axis("off")
plt.tight_layout()

Alignment

fat2=fat2[result_fat2,:]
health1=health1[result_health1,:]

lm_fat2, spatial_fat2 = rigid_alignment(lm_fat2_reference, lm_fat2, fat2.obsm['spatial'])
lm_fat2, spatial_fat2 = non_rigid_alignment(lm_fat2_reference, lm_fat2, spatial_fat2)

lm_health1, spatial_health1 = rigid_alignment(lm_health1_reference, lm_health1, health1.obsm['spatial'])
lm_health1, spatial_health1 = non_rigid_alignment(lm_health1_reference, lm_health1, spatial_health1)

Information transfer

composite_1 = composite.copy()
fat2_index=[]
health1_index=[]
composite_1_index=[]

for i in range(fat2.shape[1]):
    if(fat2.var_names[i] in health1.var_names):
        if(fat2.var_names[i] in composite_1.var_names):
            fat2_index.append(i)
            health1_index.append(np.where(health1.var_names==fat2.var_names[i])[0][0])
            composite_1_index.append(np.where(composite_1.var_names==fat2.var_names[i])[0][0])

fat2=fat2[:,fat2_index]
health1=health1[:,health1_index]
composite_1=composite_1[:,composite_1_index]

fat2.X = fat2.layers['counts']
sc.pp.normalize_total(fat2, target_sum=1e4, inplace=True)
sc.pp.log1p(fat2)
sc.tl.pca(fat2, svd_solver='arpack')

health1.X = health1.layers['counts']
sc.pp.normalize_total(health1, target_sum=1e4, inplace=True)
sc.pp.log1p(health1)
sc.tl.pca(health1, svd_solver='arpack')

WARNING: adata.X seems to be already log-transformed.
WARNING: adata.X seems to be already log-transformed.

c:\Anaconda3\envs\align\lib\site-packages\scanpy\preprocessing\_normalization.py:169: UserWarning: Received a view of an AnnData. Making a copy.
  view_to_actual(adata)
c:\Anaconda3\envs\align\lib\site-packages\scanpy\preprocessing\_normalization.py:169: UserWarning: Received a view of an AnnData. Making a copy.
  view_to_actual(adata)

x1 = torch.tensor(composite.obsm['spatial']).float()
x2 = torch.tensor(spatial_health1).float()
y2 = torch.tensor(health1.obsm['X_pca'].copy()).float()
pre_health1_pca = Exact_GP(x1, x2, y2, device, num_epochs = 500)
if(device == 'cpu'):
    pre_health1_pca = pre_health1_pca.detach().numpy()
else:
    pre_health1_pca = pre_health1_pca.detach().cpu().numpy()
pre_health1_pca_all = np.dot(pre_health1_pca,health1.varm['PCs'].T)

x1 = torch.tensor(composite.obsm['spatial']).float()
x2 = torch.tensor(spatial_fat2.copy()).float()
y2 = torch.tensor(fat2.obsm['X_pca'].copy()).float()
pre_fat2_pca = Exact_GP(x1, x2, y2, device, num_epochs = 500)
if(device == 'cpu'):
    pre_fat2_pca = pre_fat2_pca.detach().numpy()
else:
    pre_fat2_pca = pre_fat2_pca.detach().cpu().numpy()
pre_fat2_pca_all = np.dot(pre_fat2_pca,fat2.varm['PCs'].T)

pre_health1_pca_all_to01=pre_health1_pca_all.copy()
for j in range(pre_health1_pca_all_to01.shape[1]):
    if(np.max(pre_health1_pca_all[:,j])!= np.min(pre_health1_pca_all[:,j])):
        pre_health1_pca_all_to01[:,j]=(pre_health1_pca_all[:,j]-np.min(pre_health1_pca_all[:,j]))/(np.max(pre_health1_pca_all[:,j])-np.min(pre_health1_pca_all[:,j]))
    else:
        pre_health1_pca_all_to01[:,j]=np.repeat(0.5, pre_health1_pca_all_to01.shape[0])
pre_fat2_pca_all_to01=pre_fat2_pca_all.copy()
for j in range(pre_fat2_pca_all_to01.shape[1]):
    if (np.max(pre_health1_pca_all[:,j])!= np.min(pre_health1_pca_all[:,j])):
        pre_fat2_pca_all_to01[:,j]=(pre_fat2_pca_all[:,j]-np.min(pre_fat2_pca_all[:,j]))/(np.max(pre_fat2_pca_all[:,j])-np.min(pre_fat2_pca_all[:,j]))
    else:
        pre_fat2_pca_all_to01[:,j]=np.repeat(0.5, pre_health1_pca_all_to01.shape[0])

Find suspicious genes

fat2_cor=[]
health1_cor=[]
for i in range(fat2.shape[1]):
    health1_cor.append(np.corrcoef(composite_1.X[:,i],pre_health1_pca_all_to01[:,i])[0,1])
    fat2_cor.append(np.corrcoef(composite_1.X[:,i],pre_fat2_pca_all_to01[:,i])[0,1])

c:\Anaconda3\envs\align\lib\site-packages\numpy\lib\function_base.py:2897: RuntimeWarning: invalid value encountered in divide
  c /= stddev[:, None]
c:\Anaconda3\envs\align\lib\site-packages\numpy\lib\function_base.py:2898: RuntimeWarning: invalid value encountered in divide
  c /= stddev[None, :]

d = []
for i in range(len(health1_cor)):
    d.append(health1_cor[i] - fat2_cor[i])

D = np.array(d)
nan_indices = np.where(~np.isnan(D))[0].tolist()
index = health1.var_names.tolist()

d = [item for i, item in enumerate(d) if i in nan_indices]
index = [item for i, item in enumerate(index) if i in nan_indices]
health1_cor = [item for i, item in enumerate(health1_cor) if i in nan_indices]
fat2_cor = [item for i, item in enumerate(fat2_cor) if i in nan_indices]

sq.gr.spatial_neighbors(data_slice1)
sq.gr.spatial_autocorr(data_slice1,mode="moran",)

c:\Anaconda3\envs\align\lib\site-packages\squidpy\gr\_utils.py:198: ImplicitModificationWarning: Setting element `.obsp['spatial_connectivities']` of view, initializing view as actual.
  obj[key] = data

def get_top_n_indices(lst, n):
    sorted_indices = sorted(range(len(lst)), key=lambda i: lst[i], reverse=True)
    return sorted_indices[:n]

for i in get_top_n_indices(d, 4):
    if data_slice1.uns['moranI'].loc[index[i]]['I']>0.1:
        print(index[i])
        print(health1_cor[i])
        print(fat2_cor[i])

FOS
0.5130907010640495
-0.44299015007410847
SERPINA1
0.6735148710636133
-0.14194582237725056
ADIRF
0.7588274789972049
0.20245839126783868

plt.figure(dpi=300,figsize=(60,30))

plt.subplot(3,6, 1)
plt.scatter(composite_1.obsm['spatial'][:,0],
            composite_1.obsm['spatial'][:,1],
            c = data_slice1.X[:,np.where(data_slice1.var_names=='ADIRF')[0][0]],
            vmin=data_slice1.X[:,np.where(data_slice1.var_names=='ADIRF')[0][0]].min(),
            vmax=np.percentile(data_slice1.X[:,np.where(data_slice1.var_names=='ADIRF')[0][0]], 99.5),
            cmap='PuBu',
            s=20)
plt.axis("off")

plt.subplot(3,6, 2)
plt.scatter(composite_1.obsm['spatial'][:,0], composite_1.obsm['spatial'][:,1], cmap='PuBu',c=composite_1.X[:, np.where(composite_1.var_names == 'ADIRF')[0][0]],s=20)
plt.axis("off")

plt.subplot(3,6,3)
plt.scatter(health1.obsm['spatial'][:,0], health1.obsm['spatial'][:,1], cmap='PuBu',c=health1.X[:, np.where(health1.var_names == 'ADIRF')[0][0]],s=1000)
plt.axis("off")

plt.subplot(3,6,4)
plt.scatter(composite_1.obsm['spatial'][:,0], composite_1.obsm['spatial'][:,1], cmap='PuBu',c=pre_health1_pca_all_to01[:, np.where(health1.var_names == 'ADIRF')[0][0]],s=20)
plt.axis("off")

plt.subplot(3,6,5)
plt.scatter(fat2.obsm['spatial'][:,0], fat2.obsm['spatial'][:,1], cmap='PuBu',c=fat2.X[:, np.where(health1.var_names == 'ADIRF')[0][0]],s=1000)
plt.axis("off")

plt.subplot(3,6,6)
plt.scatter(composite_1.obsm['spatial'][:,0], composite_1.obsm['spatial'][:,1], cmap='PuBu',c=pre_fat2_pca_all_to01[:, np.where(fat2.var_names == 'ADIRF')[0][0]],s=20)
plt.axis("off")

plt.subplot(3,6, 7)
plt.scatter(composite_1.obsm['spatial'][:,0],
            composite_1.obsm['spatial'][:,1],
            c = data_slice1.X[:,np.where(data_slice1.var_names=='SERPINA1')[0][0]],
            vmin=data_slice1.X[:,np.where(data_slice1.var_names=='SERPINA1')[0][0]].min(),
            vmax=np.percentile(data_slice1.X[:,np.where(data_slice1.var_names=='SERPINA1')[0][0]], 99.5),
            cmap='PuBu',
            s=20)
plt.axis("off")

plt.subplot(3,6,8)
plt.scatter(composite_1.obsm['spatial'][:,0], composite_1.obsm['spatial'][:,1], cmap='PuBu',c=composite_1.X[:, np.where(composite_1.var_names == 'SERPINA1')[0][0]],s=20)
plt.axis("off")

plt.subplot(3,6,9)
plt.scatter(health1.obsm['spatial'][:,0], health1.obsm['spatial'][:,1], cmap='PuBu',c=health1.X[:, np.where(health1.var_names == 'SERPINA1')[0][0]],s=1000)
plt.axis("off")

plt.subplot(3,6,10)
plt.scatter(composite_1.obsm['spatial'][:,0], composite_1.obsm['spatial'][:,1], cmap='PuBu',c=pre_health1_pca_all_to01[:, np.where(health1.var_names == 'SERPINA1')[0][0]],s=20)
plt.axis("off")

plt.subplot(3,6,11)
plt.scatter(fat2.obsm['spatial'][:,0], fat2.obsm['spatial'][:,1], cmap='PuBu',c=fat2.X[:, np.where(health1.var_names == 'SERPINA1')[0][0]],s=1000)
plt.axis("off")

plt.subplot(3,6,12)
plt.scatter(composite_1.obsm['spatial'][:,0], composite_1.obsm['spatial'][:,1], cmap='PuBu',c=pre_fat2_pca_all_to01[:, np.where(fat2.var_names == 'SERPINA1')[0][0]],s=20)
plt.axis("off")

plt.subplot(3,6, 13)
plt.scatter(composite_1.obsm['spatial'][:,0],
            composite_1.obsm['spatial'][:,1],
            c = data_slice1.X[:,np.where(data_slice1.var_names=='FOS')[0][0]],
            vmin=data_slice1.X[:,np.where(data_slice1.var_names=='FOS')[0][0]].min(),
            vmax=np.percentile(data_slice1.X[:,np.where(data_slice1.var_names=='FOS')[0][0]], 99.5),
            cmap='PuBu',
            s=20)
plt.axis("off")

plt.subplot(3,6,14)
plt.scatter(composite_1.obsm['spatial'][:,0], composite_1.obsm['spatial'][:,1], cmap='PuBu',c=composite_1.X[:, np.where(composite_1.var_names == 'FOS')[0][0]],s=20)
plt.axis("off")

plt.subplot(3,6,15)
plt.scatter(health1.obsm['spatial'][:,0], health1.obsm['spatial'][:,1], cmap='PuBu',c=health1.X[:, np.where(health1.var_names == 'FOS')[0][0]],s=1000)
plt.axis("off")

plt.subplot(3,6,16)
plt.scatter(composite_1.obsm['spatial'][:,0], composite_1.obsm['spatial'][:,1], cmap='PuBu',c=pre_health1_pca_all_to01[:, np.where(health1.var_names == 'FOS')[0][0]],s=20)
plt.axis("off")

plt.subplot(3,6,17)
plt.scatter(fat2.obsm['spatial'][:,0], fat2.obsm['spatial'][:,1], cmap='PuBu',c=fat2.X[:, np.where(health1.var_names == 'FOS')[0][0]],s=1000)
plt.axis("off")

plt.subplot(3,6,18)
plt.scatter(composite_1.obsm['spatial'][:,0], composite_1.obsm['spatial'][:,1], cmap='PuBu',c=pre_fat2_pca_all_to01[:, np.where(fat2.var_names == 'FOS')[0][0]],s=20)
plt.axis("off")

plt.tight_layout()