> Modules non standards > Autres modules non standards > anndata et scanpy
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.
Copyright python-simple.com
programmer en python, tutoriel python, graphes en python, Aymeric Duclert