Tutorial 2 - Moran’s I
Demo data for the tutorial can be downloaded from Zenodo.
import numpy as np
import seaborn as sns
from matplotlib import pyplot as plt
import SplIsoFind
import warnings
warnings.filterwarnings("ignore", category=FutureWarning)
Construct isoform matrix with relative expression
The create isoform matrix saves a sparse matrix, and corresponding index and column labels to the output folder.
NOTE: this is the same step as we do before plotting the results, you only have to do this once
fn_allinfo = 'data/allinfo_ds.filtered.labeled.gz'
fn_CIDmap = 'data/sample1_barcodeToPos.CellID_ds.tsv.gz'
fn_adata = 'data/sample1_cellbin_adjusted.h5ad'
output = 'data/isoform_matrix'
SplIsoFind.pp.create_isoform_matrix(fn_allinfo,
fn_CIDmap,
fn_adata,
output)
Potentially interesting isoforms (total): 145
Potentially interesting isoforms (novel): 44
x_sparse, labels, isoforms = SplIsoFind.pp.load_sparse(input_dir = output)
Calculate Moran’s I
For all cells (‘All’) and for the four major cell types.
The running time mainly depends on the number of cells, number of tested isoforms, and number of permutations. Doing 1e6 permutations will take ~10min for this demo dataset.
The moran’s I function returns 2 dataframes:
pval: the uncorrected p-values
qval: Benjamini-Yekutieli corrected values
NaN indicates that a specific isoform-celltype combination did not fulfill the testing criteria (e.g. not enough cells).
mI_scores, pval, qval = SplIsoFind.sv.moransI_sparse(x_sparse,
labels,
isoforms,
celltypes=['All','ExciteNeuron','InhibNeuron','Astro','Oligo'],
output_dir='data/moransI',
n_jobs=8)
Number of significant isoforms per celltype
(qval <= 0.05).sum()
All 20
ExciteNeuron 7
InhibNeuron 2
Astro 0
Oligo 0
dtype: int64
Cell-type constrained permutations
For the significant isoforms using all cells, we run Moran’s I with cell-type constrained permutations. This way we can check whether significance is caused by one or more cell types causing splicing changes (p-value will still be significant) or because of different cell-type proportions between the regions (p-value won’t be significant anymore).
This function is significantly slower, so we only run it on the significant isoforms. This takes ~20 min on the demo dataset with n_jobs=8.
iso_sign = qval.index[qval['All'] <= 0.05]
iso_sign[:1]
Index(['ENSMUST00000051477.13'], dtype='str', name='Transcript ID')
res = SplIsoFind.sv.moransI_ctperm_sparse(x_sparse,
labels,
isoforms,
var_totest=iso_sign,
output_dir='data/moransI',
n_jobs=8)
g = sns.jointplot(
x=-np.log10(res['p-value (original)']),
y=-np.log10(res['p-value (new)']),
marginal_kws={"binwidth": 0.1, "binrange": [0.05, 4.05]},
height=2.5, s=5,
ratio=4
)
# Add x = y line
lims = [
np.min([g.ax_joint.get_xlim(), g.ax_joint.get_ylim()]), # min of both axes
np.max([g.ax_joint.get_xlim(), g.ax_joint.get_ylim()]) # max of both axes
]
g.ax_joint.plot(lims, lims, '--', color='gray', linewidth=0.5) # diagonal line
g.ax_joint.set_xlim(lims)
g.ax_joint.set_ylim(lims)
# Set axis labels
g.set_axis_labels('Normal permutations \n-log10(p-value)',
'Cell-type constrained \n-log10(p-value)')
plt.show()
For most isoforms, the p-value of the normal and cell-type constrained permutations are very similar. For some, however, we see a big difference. Let’s plot an example of both.
res.sort_values('p-value (new)')
| morans I | p-value (original) | p-value (new) | Num cells | Imbalance | |
|---|---|---|---|---|---|
| variable | |||||
| ENSMUST00000211044.2 | 0.134694 | 0.000100 | 0.000100 | 9937 | 0.721143 |
| ENSMUST00000211283.2 | 0.158175 | 0.000100 | 0.000100 | 9937 | 0.322029 |
| ENSMUST00000028727.11 | 0.046870 | 0.000100 | 0.000100 | 9762 | 0.070580 |
| ENSMUST00000110098.4 | 0.046870 | 0.000100 | 0.000100 | 9762 | 0.907191 |
| ENSMUST00000188368.7 | 0.033126 | 0.000100 | 0.000100 | 3459 | 0.908066 |
| ENSMUST00000234496.2 | 0.127653 | 0.000100 | 0.000100 | 2061 | 0.658903 |
| transcript5764.chr18.nnic | 0.090287 | 0.000100 | 0.000100 | 2061 | 0.737021 |
| ENSMUST00000107745.8 | 0.046292 | 0.000100 | 0.000100 | 769 | 0.585176 |
| ENSMUST00000075316.10 | 0.033701 | 0.000100 | 0.000100 | 769 | 0.568270 |
| ENSMUST00000121334.8 | 0.210834 | 0.000100 | 0.000100 | 989 | 0.407482 |
| ENSMUST00000120878.9 | 0.206663 | 0.000100 | 0.000100 | 989 | 0.608696 |
| ENSMUST00000025239.9 | 0.042852 | 0.000100 | 0.000200 | 2061 | 0.669093 |
| ENSMUST00000191011.7 | 0.020798 | 0.000100 | 0.000300 | 3459 | 0.226077 |
| ENSMUST00000165617.3 | 0.024483 | 0.000200 | 0.004700 | 1258 | 0.330684 |
| transcript32816.chr10.nnic | 0.027521 | 0.006699 | 0.006299 | 314 | 0.907643 |
| ENSMUST00000025931.14 | 0.023167 | 0.000100 | 0.009499 | 1258 | 0.669316 |
| ENSMUST00000051477.13 | 0.040551 | 0.000100 | 0.010399 | 1129 | 0.267493 |
| transcript44056.chr4.nnic | 0.039395 | 0.000100 | 0.011899 | 1129 | 0.735164 |
| ENSMUST00000030187.14 | 0.027611 | 0.001600 | 0.220878 | 544 | 0.606618 |
| transcript11577.chr4.nnic | 0.027611 | 0.001500 | 0.234177 | 544 | 0.389706 |
from PIL import Image
Image.MAX_IMAGE_PIXELS = 553190400
im_reg = Image.open('data/Sample1_ssDNA_regist.tif')
imarray = np.array(im_reg)
xlim1=5000
xlim2=17000
ylim1=4450
ylim2=18000
For example, within Snap25 has the lowest p-value possible using both normal and cell-type constrained permutations and we see spatial variation within excitatory and inhibitory neurons
isoform = 'ENSMUST00000028727.11'
res.loc[isoform]
morans I 0.04687
p-value (original) 0.00010
p-value (new) 0.00010
Num cells 9762.00000
Imbalance 0.07058
Name: ENSMUST00000028727.11, dtype: float64
SplIsoFind.pl.spatial_hexplot_sparse(x_sparse,
labels,
isoforms,
imarray=imarray,
varName=isoform,
celltype='',
hexsize=350,
plot_lim=[xlim1, xlim2, ylim1, ylim2])
plt.show()
for ct in ['ExciteNeuron', 'InhibNeuron', 'Oligo', 'Astro']:
SplIsoFind.pl.spatial_hexplot_sparse(x_sparse,
labels,
isoforms,
imarray=imarray,
varName=isoform,
celltype=ct,
hexsize=350,
plot_lim=[xlim1, xlim2, ylim1, ylim2])
plt.show()
A counter example is an isoform of Tln1. Oligodendrocytes mostly use this isoform (mean PI ≈ 0.9), whereas excitatory neurons mostly do not (≈ 0.2). The spatial pattern using all cells therefore mainly reflects where the oligodendrocytes are (e.g. white matter and midbrain) rather than variation within a cell type, and the p-value is no longer significant with cell-type constrained permutations.
isoform = 'ENSMUST00000030187.14'
res.loc[isoform]
morans I 0.027611
p-value (original) 0.001600
p-value (new) 0.220878
Num cells 544.000000
Imbalance 0.606618
Name: ENSMUST00000030187.14, dtype: float64
SplIsoFind.pl.spatial_hexplot_sparse(x_sparse,
labels,
isoforms,
imarray=imarray,
varName=isoform,
celltype='',
hexsize=350,
plot_lim=[xlim1, xlim2, ylim1, ylim2])
plt.show()
for ct in ['ExciteNeuron', 'InhibNeuron', 'Oligo', 'Astro']:
SplIsoFind.pl.spatial_hexplot_sparse(x_sparse,
labels,
isoforms,
imarray=imarray,
varName=isoform,
celltype=ct,
hexsize=350,
plot_lim=[xlim1, xlim2, ylim1, ylim2])
plt.show()
Less than 10 cells to plot, skipping.