anndata et scanpy

anndata

package qui permet de stocker des matrices annotées

Définition d'un objet :

values = numpy.random.normal(1, size = (100, 200)).astype(numpy.float32); adata = anndata.AnnData(values)
les variables sont en colonnes, les observations en lignes (comme un dataframe pandas).
adata.n_vars : le nombre de variables.
adata.n_obs : le nombre d'observations.
adata.shape : paire nombre de lignes (observations), nombre de colonnes (variables).
adata.X : la matrice des observations. C'est une matrice sparse de type scipy.sparse.csr.csr_matrix.
adata.X.toarray() : récupère la matrix dense de type numpy array, mais attention à ne pas faire ça sur une très grosse matrice !
on peut donner des noms aux observations : adata.obs_names = ['cell' + str(i) for i in range(100)]
on peut donner des noms aux variables : adata.var_names = ['gene' + str(i) for i in range(200)]

Metadata :

on peut rajouter des metadata sur les variables (colonnes) et les observations (lignes)
attaché à l'objet, on a 2 dataframes de metadata : adata.obs pour les observations, adata.var pour les variables. Ils sont initialement vides.
adata.obs['cellType'] = ['typeA' if i < 50 else 'typeB' for i in range(100)] : rajout d'une colonne au dataframe adata.obs
si c'est une variable de type catégorie, faire plutôt : adata.obs['cellType'] = pandas.Categorical(['typeA' if i < 50 else 'typeB' for i in range(100)])
on peut aussi rajouter des metadata avec plusieurs dimensions :
- adata.obsm['umap'] = numpy.zeros((100, 2)) : la première dimension de la matrice doit être égale au nombre d'observations.
- adata.varm['geneProp'] = numpy.zeros((200, 5)) : la première dimension de la matrice doit être égale au nombre de variables.
adata.uns : on peut rajouter aussi des informations générales dans le champ uns (unstructured) : adata.uns['info'] = {'a': 1, 'b': 3} (les valeurs peuvent être de n'importe quel type.
var_keys(), obs_keys(), varm_keys(), obsm_keys, uns_keys() : méthodes qui retournent les clefs des différents champs.
on peut rajouter aussi des layers qui doivent avoir les mêmes dimensions que la matrice principale : adata.layers['logTransf'] = numpy.log10(adata.X). Ce sont des array numpy.
adata.layers.keys() : les noms des layers
adata.layers['logTransf'] : pour accéder à une layer avec le nom donné.
adata.layers['logTransf'][:, ['CD8', 'CD25']] : pour accéder à certaines valeurs.

Manipulation :

adata[0:5, 0:3] : subsetting avec les index.
adata[[0, 1], [0, 2]] : subsetting avec les index.
adata[['cell0', 'cell1'], ['gene0', 'gene2']] : subsetting grâce aux noms donnés.
adata[adata.obs.cellType == 'typeA'] : subsetting grâce aux metadata.
quand on fait un subsetting d'un objet AnnData, ça crée seulement une vue sur l'objet existant, pas de nouvelle création.
si on veut créer un objet indépendant, il faut faire une copie : adata2 = adata.copy()
adata.to_df() : renvoie un dataframe de la matrice.
adata.to_df(layer = 'logTransf') : renvoie un dataframe d'une autre layer que la principale.
adata.transpose() : transpose tout l'objet.
adata.chunk_X() : renvoie un sample aléatoire (mais reproductible) de adata.X (1000 lignes et colonnes par défaut).

Sauvegarde dans un fichier :

adata.write('myFile.h5ad', compression = 'gzip') (= adata.write_h5ad) : sauvegarde dans un fichier HDF5
adata = anndata.read('myFile.h5ad') (= anndata.read_h5ad) : pour relire le fichier.
on peut lire le fichier partiellement avec adata = anndata.read('myFile.h5ad', backed = 'r'). L'objet mémorise alors son fichier associé (avec adata.filename) et il faut le fermer quand on a fini : adata.file.close()
adata.write_csvs('myDir', skip_data = False, sep = '\t') : écrit dans myDir toutes les données au format csv avec une tabulation comme séparateur (par défaut, skip_data = False, et il n'écrit pas la matrice). Les noms des fichiers sont les noms des champs : X.csv, obs.csv, var.csv, uns.csv, ...

Concatenation d'objets selon les observations : adata = anndata.AnnData.concatenate(adata1, adata2, , adata3, batch_categories = ['a', 'b', 'c']) :

adata.obs a une colonne batch avec la valeur indiquée.
les noms des observations sont suffixés par le nom de leur batch (ici, en rajoutant -a, -b ou -c selon le batch).

scanpy

Récupération de dataframe :

scanpy.get.obs_df(adata, ['leiden', 'Phase', 'B2M']) : récupère un dataframe avec une ligne par cellule avec les valeurs indiquées. On peut donner comme colonnes les valeurs dans adata.obs.columns et dans adata.var_names. S'il y a des gènes, ce sont les valeurs mise à l'échelle qui sont récupérées.
scanpy.get.obs_df(adata, ['leiden', 'Phase', 'B2M'], use_raw = True) : récupère un dataframe avec une ligne par cellule avec les valeurs indiquées. Valeurs brutes pour les gènes (défaut : use_raw = False).
scanpy.get.var_df(adata, ['total_counts', 'pct_dropout_by_counts']) : récupère un dataframe avec une ligne par gène. On peut donner comme colonnes les valeurs dans adata.var.columns et éventullement dans adata.obs_names (barcodes des cellules).

Pré-processing :

scanpy.pp.filter_cells(adata, min_counts = 200, inplace = True) : garde seulement les cellules qui ont au moins 200 counts et modifie l'objet sur place (défaut) en rajoutant n_counts à adata.obs. Si inplace = False, ne modifie pas l'objet mais renvoie un tuple : array de booléens disant s'il faut garder la cellule, nombre de comptes pour chaque cellule.
scanpy.pp.filter_cells(adata, min_genes = 200, inplace = True) : même chose, mais en faisant les comptes de genes et en rajoutant n_genes à adata.obs (il y a aussi max_counts et max_genes).
scanpy.pp.filter_genes(adata, min_cells = 5, inplace = True) : garde seulement les gènes qui sont exprimés dans au moins 5 cellules et rajoute n_cells à adata.var. Idem avec max_cells, min_counts, max_counts (pour les deux derniers, rajoute n_counts à adata.var).
annData.var['mito'] = annData.var_names.str.startswith('MT-') : rajoute une annotation sur les gènes mitochrondriaux.
scanpy.pp.calculate_qc_metrics(annData, qc_vars = ['mito'], log1p = False, inplace = True) : si on a rajouté une annotation sur les gènes mitochrondriaux, rajoute :
- à adata.var : n_cells_by_counts, mean_counts, pct_dropout_by_counts, total_counts.
- à adata.obs : n_genes_by_counts, total_counts, pct_counts_in_top_50_genes, pct_counts_in_top_100_genes, pct_counts_in_top_200_genes, pct_counts_in_top_500_genes, total_counts_mito, pct_counts_mito.
scanpy.pp.normalize_total(adata, target_sum = 10000) : normalise les comptes pour que chaque cellule ait 10000 comptes au total (on peut le vérifier avec annData.X.sum(axis = 1)).
scanpy.pp.log1p(adata) : normalise les données en prenant le log(x + 1).
- par défaut, l'objet est modifié en place
- si on passe copy = True, ne modifie pas l'objet, mais retourne une copie modifiée.
- scanpy.pp.log1p(adata, layer = 'normalized') : normalise les données de la layer indiquée qui doit avoit été crée avant (permet de conserver les données d'origine dans l'objet)
scanpy.pp.highly_variable_genes(adata, min_mean = 0.0125, max_mean = 3, min_disp = 0.5) : calcule les gènes hautement variables :
- modifie en place l'objet en rajoutant des colonnes à adata.var.
- adata.var['highly_variable'] : booléen qui indique les gènes qui sont variables.
- on peut afficher un graphe qui montre les gènes sélectionnés : scanpy.pl.highly_variable_genes(adata) (scatterplot utilisant les champs means, dispersions, dispersions_norm et highly_variable de adata.var).
on peut alors travailler uniquement sur ces gènes en faisant : adata = adata[:, adata.var['highly_variable']
mais il est pratique de conserver les valeurs sur tous les gènes en faisant avant adata.raw = adata. On peut accéder à l'objet d'origine avec tous les gènes par : adata.raw.to_adata()
On peut enlever l'influence de certaines variables par : scanpy.pp.regress_out(adata, ['total_counts', 'pct_counts_mt'], copy = False)
- par défaut, update sur place la matrice (copy = False). Utiliser copy = True pour renvoyer une copie.
- adata.X devient alors une matrice dense (numpy array).
scanpy.pp.scale(adata) :
- calcule les données normalisées pour avoir une moyenne 0 et écart-type de 1 pour chaque variable.
- si copy = False (défaut), rajoute dans adata.var les clefs mean et std et modifie les données sur place, sinon renvoie une copie.

Processing :

scanpy.tl.pca(adata) (ou scanpy.pp.pca(adata), c'est la même fonction) : fait une PCA (par défaut, 50 dimensions) et rajoute les champs suivants :
- n_comps = 30 : 30 dimensions pour la PCA.
- adata.osbm['X_pca'] : une array numpy n_obs x nbr dimensions de la PCA avec les coordonnées des observations sur les axes.
- adata.varm['PCs'] : les loadings (n_genes x nbr dimensions de la PCA).
- adata.uns['pca']['variance'] : les variances le long de chaque axe.
- adata.uns['pca']['variance_ratio'] : le ratio de chaque variance de chaque axe par rapport à leur somme (la somme des ratios vaut 1).
scanpy.pp.neighbors(adata) : calcule le graphe des plus proches voisins. On peut donner :
- n_neighbors : le nombre de voisins (défaut = 15).
- n_pcs : le nombre d'axes principaux à utiliser.
la fonction rajoute les champs suivants :
- adata.obsp['distances'] : la matrice sparse des distances entre les observations.
- adata.obsp['connectivities'] : la matrice d'adjacence sparse du graphe avec les connectivités comme poids.
- adata.uns['neighbors'] : les paramètres d'appel de la fonction.
scanpy.tl.umap(adata) : fait une projection UMAP et rajoute le champs suivant :
- adata.obsm['X_umap'] : les coordonnées UMAP des observations.
scanpy.tl.tsne(adata) : fait une projection tSNE et rajoute le champs suivant :
- adata.obsm['X_tsne] : les coordonnées tSNE des observations.
scanpy.tl.leiden(adata) : fait un clustering de Leiden (amélioration de Louvain). Paramètres :
- resolution : règle le nombre de cluster (défaut vaut 1, le nombre de cluster augmente avec la resolution.
cela rajoute les champs suivants :
- adata.obs['leiden'] : numéro du cluster pour chaque observation (qui est de type catégorie).
- adata.uns['leiden'] : les paramètres utilisés.
scanpy.tl.louvain(adata) : fait un clustering de Louvain. Paramètres :
- resolution : règle le nombre de cluster (défaut vaut 1, le nombre de cluster augmente avec la resolution.
cela rajoute les champs suivants :
- adata.obs['louvain'] : numéro du cluster pour chaque observation.
- adata.uns['louvain'] : les paramètres utilisés.

PAGA : permet de calculer des trajectoires :

scanpy.tl.paga(adata, groups = 'leiden') : calcule le graphe PAGA. Les données sont rajoutées dans adata.uns['paga'].
scanpy.pl.paga(adata, color = ['leiden', 'MPO'], cmap = 'Reds') : trace les graphes PAGA, soit colorés par le numéro du cluster, soit coloré par l'intensité du gène MPO (2 graphes similaires).
scanpy.tl.umap(adata, init_pos = 'paga') : recalcule une umap, mais en prenant comme point de départ le graphe PAGA.
scanpy.pl.paga_compare(adata) : affiche à la fois la umap initialisée à partir du graphe PAGA et le graphe PAGA, pour comparer les deux.

Expression différentielle des gènes :

scanpy.tl.rank_genes_groups(adata, groupby = 'leiden', method = 'wilcoxon', pts = True)
scanpy.tl.rank_genes_groups(adata, groupby = ['0'], reference = '1', n_genes = 20, method = 'wilcoxon', pts = True) : fait la comparaison seulement sur le cluster 0 et par rapport au cluster 1 (et non toutes les autres cellules) et en ne retenant que les 20 premiers gènes.

le résultat est dans adata.uns['rank_genes_groups'] sous forme d'un dictionnaire dont les valeurs sont des recarray indexées par l'identifiant du groupe :

result = adata.uns['rank_genes_groups']
groups = result['names'].dtype.names
for group in groups:
    df = pandas.DataFrame({x:result[x][group] for x in ['names', 'pvals_adj', 'logfoldchanges', 'pts', 'pts_rest']})
    print('resultats pour la comparaison : ' + group)
    df.to_csv(sys.stdout, sep = '\t')

scanpy.get.rank_genes_groups_df(adata, 'HSC') : dataframe avec les gènes differentiellement exprimés. Les colonnes sont : names, scores, logfoldchanges, pvals, pvals_adj, pct_nz_group (% non zero dans le groupe), pct_nz_reference (% non zero dans la réference).
scanpy.get.rank_genes_groups_df(adata, ['HSC', 'GMP']) permet d'avoir un seul tableau pour tous les groupes, avec une colonne supplémentaire indiquant le groupe.
on peut voir le résultat des premiers gènes graphiquement par : scanpy.pl.rank_genes_groups(adata, ncols = 2).

Plots :

les plots sont sauvés par défaut dans un directory figures : pour éviter ça : scanpy.settings.figdir = ''. Malheureusement, scanpy rajoute en plus un prefix au nom du fichier, par exemple umap ...

une façon d'éviter ce comportement génant et de fixer la taille et la résolution juste pour la figure :

with pyplot.rc_context({"figure.figsize": (8, 8), "figure.dpi": (300)}):
    scanpy.pl.umap(adata, color = 'cellTypePredictionAz', show = False)
    pyplot.savefig('outfile.png', bbox_inches = 'tight')

scanpy.pl.scatter(adata, x = 'total_counts', y = 'pct_counts_mito')
scanpy.pl.violin(adata, ['n_genes_by_counts', 'total_counts', 'pct_counts_mito'], jitter = 0.4, multi_panel = True) : graphes de QC.
scanpy.pl.violin(adata, keys = ['CD4', 'CD8'], groupby = 'sampleName', log = True)
scanpy.pl.pca_variance_ratio(adata) : graphe avec le pourcentage de variance explique pour chaque axe (utilise en fait les données de adata.uns['pca']['variance_ratio']).
scanpy.pl.pca(adata, color = 'EEF1A1') : les points sur les deux premiers axes principaux, colorés par l'expression du gène EEF1A1 :
- par défaut, c'est adata.raw qui est utilisé.
- dimensions = [0, 2] : les numéros des axes principaux à utiliser (origine à 0).
- color_map = 'Reds' : la color map à utiliser.
scanpy.pl.umap(adata) :
- scanpy.pl.umap (adata, color = 'clusterNbr') : représentation d'une donnée de type catégorie (attention, si numéro de cluster, convertir éventuellement d'abord en catégorie avec adata.obs['clusterNbr'] = adata.obs['clusterNbr'].astype('category')
- scanpy.pl.umap(adata, color = 'clusterNbr', groups = [1, 4, 5]) : ne montre que certaines valeurs (certains clusters ici)
- scanpy.pl.umap(adata, color = ['HBB', 'MPO', 'THY1'], color_map = 'hot_r') : fait un graphe par gènes, en utilisant par défaut les données de raw.
- legend_loc = 'on data' : met la légende directement sur les données.
- palette : donne les couleurs, pour une variable de type catégorie.
- title = ['gene1', 'gene2', 'gene3'] : indique les titres à mettre pour chaque graphe.
- size = 2 : la taille des points.
- ncols = 2 : 2 graphes par ligne (donc 2 colonnes).
- vmin et vmax : les valeurs min et max pour l'échelle de couleurs.
scanpy.pl.highest_expr_genes(adata, n_top = 20) : fait un boxplot des 20 gènes les plus exprimés
scanpy.pl.dotplot(adata, ['B2M', 'TUBB', 'MX1', 'GATA1'], groupby = 'leiden') : dot plot par groupe des gènes indiqués :
- dendrogram = True : rajoute un dendrogramme sur les groupes.
- swap_axes = True : les gènes sont en ligne et les groupes en colonne (défaut = False).
scanpy.pl.stacked_violin(adata, ['B2M', 'TUBB', 'MX1', 'GATA1'], groupby = 'leiden', rotation = 90) : idem que dotplot, mais sous forme de violinplot tout petit.

Divers :

adata.write('myFile.h5') : écrit l'objet sous forme de fichier hdf5.
adata.write_loom('myFile.h5') : écrit un fichier au format loom.
adata.write_csvs('myDir', sep = '\t') : crée un directory dans lequel il écrit toutes les metadata : obs, var, obsm, varm, uns.
annData = scanpy.read_loom('myFile.loom') : pour lire un fichier loom.
si les metadata sur les cellules sont dans un dataframe df :
- (df.index != annData.obs.index).sum() : pour vérifier que les deux ont bien le même index constitué des barcodes des cellules.
- annData.obs = df : on peut alors attribuer les metadata à l'objet annData.

programmer en python, tutoriel python, graphes en python, Aymeric Duclert