Ligand-receptor analysis - Julian

[giovp]

• download a data of choice from 10x genomics with sc.datasets_visium_sge() (you can browse them here)
• run the scanpy tutorial, similar to here: https://scanpy-tutorials.readthedocs.io/en/latest/spatial/basic-analysis.html
• compare marker genes in spatial dataset with marker genes in immune atlas (Scanpy function is sc.tl.rank_genes_group

Then we'll discuss next week about ligand-receptor interaction analysis with our squidpy package.

Have a look at analysis done here, potentially interesting for you: https://www.biorxiv.org/content/10.1101/2021.01.02.425073v1.full

@julianalbers

[giovp]

hi @julianalbers

so after haveing isntalled squidpy, you can run these commands for ligand receptor interaction analysis:

import squidpy as sq
sq.gr.ligrec(
    adata,
    cluster_key="cluster",
)

there are quite exhaustive docsstring, and several arguments to play around. Defaults ones are also good though. It might throw you an error re adata.raw issing. No worries, just make sure that just after loading hte anndata (so before normalization/clsutering) you assign raw counts to the slot like this adata.raw = adata

This will save results in adata.uns["<your cluster id>_ligrec"]. The object is quite complex, you can slice it as a list, for instance

adata.uns["<your cluster id>_ligrec"][0] # the means (as in cellphonedb)
adata.uns["<your cluster id>_ligrec"][1] # the pvalues (as in cellphonedb)
adata.uns["<your cluster id>_ligrec"][2] # rich metadata related to the LR pairs

we also have a convenient plotting function, with several arguments for filtering etc. For instance

# here, my cluster_key was "cluster" and I was looking at some brain clusters
sq.pl.ligrec(
    adata,
    cluster_key="cluster",
    source_groups="Hippocampus",
    target_groups=["Pyramidal_layer", "Pyramidal_layer_dentate_gyrus"],
    means_range=(3, np.inf),
    alpha=1e-4,
    swap_axes=True,
)

you can visualize this plot for a bunch of interesting clusters you find interesting. Also would be interesting to compute some statistiscs between LR pairs in clusters of visium vs. LR pairs in clusters of cancer immune cell atlas (filter for cancer of interest). For example, at a given pval thres, how many LR pairs are present in each cluster for thw two datasets ? it's like a confusion matrix with intersection counts. Or anything else you find interesting!

[julianalbers]

Hi Giovanni,

looking forward to the meeting tomorrow. I was just wondering if you know of any human lymph node scRNA-seq dataset? I checked the database by Valentin Svensson (https://docs.google.com/spreadsheets/d/1En7-UV0k0laDiIfjFkdn7dggyR7jIk3WH8QgXaMOZF0/edit#gid=0) but it was not successful.

For the human lymphnode a lot of interactions are well described, so I would love to try to compare ligrec results for these data with the spatial transcriptomics data.

Thanks, see you tomorrow. Julian

[giovp]

looked at bit but couldn't find any sorry :(

[julianalbers]

This is the publication, where they specifically analyzed the leading edge, where the tumor interacts with non-tumor tissue. Figure 3: https://www.sciencedirect.com/science/article/pii/S0092867420306723

[giovp]

for enighborhood enrichment

sq.gr.nhood_enrichemnt(adata)
sq.pl.nhood_enrichment(adata)

for deconvolution: https://github.com/theislab/AutoGeneS

[giovp]

lymphnode paper where they did diffmap on it: https://www.biorxiv.org/content/10.1101/2020.11.15.378125v1

[giovp]

De analysis method comparison, relevant for everyboyd @ngocphuongtran @julianalbers @PhilippStaedter

https://www.nature.com/articles/nmeth.4612

[giovp]

massive paper on trajectory inference benchmark: https://www.nature.com/articles/s41587-019-0071-9

[giovp]

so after installing cellrank, you can try this for gene trends

import cellrank
import scanpy

sc.tl.diffmap(adata)
scanpy.tl.dpt(adata)

cr.tl.terminal_states(
    adata,
    cluster_key="clusters",
    weight_connectivities=0.2,
    n_states=3,
    softmax_scale=4,
    show_progress_bar=False,
)
cr.tl.lineages(adata)

model = cr.ul.models.GAM(adata)

cr.pl.gene_trends(
    adata,
    model,
    ["Map2", "Dcx"],
    data_key="X",
    time_key="dpt_pseudotime",
    show_progress_bar=False,
)

in your case, youd' proably have to subset the anndata first to only select b cell populations.

Also, have a look at the lig rec example, quite inteersting settings I think. You can find a similar idea here (but without ligrec):

lymphnode paper where they did diffmap on it: biorxiv.org/content/10.1101/2020.11.15.378125v1

[giovp]

very cool paper: https://www.biorxiv.org/content/10.1101/2021.02.16.431475v1.full.pdf

[julianalbers]

yes, super nice.

I'll try to do something similar like this with the lymph node data