> Modules non standards > Autres modules non standards > scanpy
scanpy
scanpy.datasets : contient des datasets de tests, par exemple adata = scanpy.datasets.pbmc3k() pour récupérer un dataset de PBMC.
adata.var_names_make_unique() : rend les noms de gènes uniques en rajoutant un -1 à la 2eme occurrence, puis -2 à la 3ème ...
Preprocessing
En général (sauf exception) :
- Par défaut, l'objet est modifié en place.
- si on passe inplace = False, pas de modification en place.
- si on passe copy = True, ne modifie pas l'objet, mais retourne une copie modifiée.
QC :
- adata.var['mt'] = adata.var_names.str.startswith('MT-') : rajout de l'info si le gène est mitochondrial.
- dfCell, dfGene = scanpy.pp.calculate_qc_metrics(adata, qc_vars = ['mt']) : ne modifie pas adata, mais renvoie 2 dataframes : un dataframe cellules pour les QC par cellule, et un dataframe par gène :
- dataframe par cellule : n_genes_by_counts (nombre de gènes pour chaque cellule), total_counts (total du nombre de reads), total_counts_mt (nombre de reads mitochondriaux si on a passé 'mt' à qc_vars).
- dataframe par gène avec n_cells_by_counts (nombre de cellules exprimant le gène), mean_counts (nombre moyen de reads de ce gène par cellule), pct_dropout_by_counts (pourcentage de drop-out du gène)
- scanpy.pp.calculate_qc_metrics(adata, qc_vars = ['mt'], inplace = True) : les informations sont rajoutées à adata.obs et adata.var (le défaut est False ici).
Filtrage des cellules et des gènes :
- scanpy.pp.filter_cells(adata, min_genes = 100, inplace = True) : garde seulement les cellules qui ont au moins 100 gènes 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. Autres paramètres :
- scanpy.pp.filter_cells(adata, max_genes = 5000) : le nombre de gènes max
- scanpy.pp.filter_cells(adata, min_counts = 1000) : le nombre de read min (rajoute n_counts).
- scanpy.pp.filter_cells(adata, max_counts = 50000) : le nombre de read max (rajoute n_counts).
- on ne peut pas donner plusieurs de ces paramètres en même temps, mais on peut faire des appels successifs !
- 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), mais on ne peut donner qu'un seul paramètre à la fois (faire des appels successifs).
Normalisation des comptes :
- adata.layers['counts'] = adata.X.copy() : permet de conserver une copie des comptes avant la normalisation.
- scanpy.pp.normalize_total(adata) : normalise chaque cellule à la médiane des comptes et remplace X par le résultat
- 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.normalize_total(adata, key_added = 'normFactor') : conserve dans la colonne normFactor de adata.obs le facteur de normalisation.
- scanpy.pp.log1p(adata) : transforme en log après normalisation en prenant le log(x + 1) (log népérien).
Identification des gènes fortement variables :
- scanpy.pp.highly_variable_genes(adata, n_top_genes = 2000, batch_key = 'sample') : calcule les 2000 gènes les plus variables,mais en considérant chaque batch séparemment, repéré par la clef sample dans adata.obs (pour éviter l'influence trop grande de certain batchs).
- modifie en place l'objet en rajoutant des colonnes à adata.var :
- highly_variable : booléen qui indique les gènes qui sont variables.
- means : la moyenne du gène.
- dispersions : la dispersion.
- dispersions_norm : la dispersion normalisée.
- scanpy.pl.highly_variable_genes(adata, show = False) : montre un plot avec les gènes fortement variables.
Etapes de processing
Réduction de dimension :
- scanpy.tl.pca(adata) (ou scanpy.pp.pca(adata), c'est la même fonction) : fait une PCA, par défaut sur les gènes variables si ceux-ci ont bien été calculé avant
- paramètres :
- n_comps = 50 : nombre de dimensions à conserver (défaut de 50).
- use_highly_variable = False : utilise tous les gènes plutôt que ceux les plus variables.
- Résultats :
- 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.pl.pca_variance_ratio(adata, n_pcs = 50, log = True, show = False) : montre la variance pour chaque dimension.
- scanpy.pl.pca(adata, color = ['sample', 'sample'], dimensions = [(0, 1), (2, 3)], ncols = 1, size = 40, show = False) : fait 2 graphes coloré par la clef sample de adata.obs sur les paires de dimensions (0, 1) et (2, 3)
Calcul des plus proches voisins :
- 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.
Clustering :
- nécessite d'avoir calculé les plus proches voisins.
- scanpy.tl.leiden(adata, resolution = 1) : fait un clustering avec la résolution donnée (le nombre de clusters augmente avec la résolution)
- key_added : le nom de la clef à rajouter à adata.obs, par défaut, 'leiden'.
- par défaut, une colonne leiden est rajoutée à adata.obs avec le numéro du cluster (qui est de type catégorie)
- on peut aussi faire un clustering de louvain : scanpy.tl.louvain(adata, resolution = 1) : fait un clustering (rajoute le numéro du cluster à adata.obs dans la colonne louvain)
Calcul de la UMAP :
- nécessite d'avoir calculé les plus proches voisins.
- scanpy.tl.umap(adata) : fait une projection UMAP.
- paramètres :
- min_dist : la distance minimum entre les points (défaut = 0.5)
- spread : l'échelle des points. Détermine avec min_dist si les points sont plus ou moins étalés (défaut = 1)
- rajoute les champs suivants :
- adata.obsm['X_umap'] : la matrice nombre de cellules x 2 avec les coordonnées UMAP des observations.
- on peut aussi faire une projection tSNE : scanpy.tl.tsne(adata). Rajoute une clef adata.obsm['X_tsne']
Force-directed graph drawing :
- alternative à UMAP/tSNE.
- nécessite d'avoir calculé les plus proches voisins.
- scanpy.tl.draw_graph(adata, layout = 'fa', show = False) : calcule le mapping en utilisant l'algorithme ForceAtlas2 (défaut). Ce défaut nécessite d'avoir le package python fa2-modified installé.
- autres possibilités pour layout : 'fr', 'grid_fr', 'lgl', 'drl', 'rt' et d'autres.
- une clef est rajoutée à adata.obsm : 'X_draw_graph_fa' par défaut, qui est une matrice nombre de cellules x 2 ('X_draw_graph_fr' si layout est 'fr', etc ...)
- random_state : permet d'avoir un résultat reproductible (vaut 0 par défaut).
On peut alors tracer la figure avec scanpy.pl.draw_graph :
- scanpy.pl.draw_graph(adata, color = 'sample', layout = 'fa') : trace le graphe correspondant avec des couleurs par sample, mais on peut aussi donner un nom de gène.
- size = 50 : taille des points (défaut est un calcul automatique.
- layout = 'fa' : le layout à utiliser (le défaut est d'utiliser le dernier calculé).
- sort_order = True : True si les points les plus forts doivent être tracés en dernier (le défaut).
- groups = ['sample1', 'sample3'] : pour restreindre à certaines catégories (les autres points sont en na_color).
- legend_loc = 'right margin'' : position de la légende ('on data', 'best', 'right margin', 'upper right', ...)
- legend_fontsize = 'medium : 'small', 'x-small', 'xx-small', 'large', ...
- color_map = 'YlOrRd' : pour les variables continues.
- palette' : pour les variables de catégorie.
- na_color = 'lightgrey' : la couleur des points masqués.
- vmin, vmax, vcenter : pour la colorbar
- ncols, wspace, hspace : quand on donne plusieurs valeurs à color.
Diffusion map :
- scanpy.tl.diffmap(adata, n_comps = 15, random_state = 0) : calcule une diffusion map à 15 dimensions (le défaut).
- adata2.obsm['X_diffmap'] la matrice nombre de cellules x 15 (si 15 dimensions).
Expression différentielle des gènes avec rank_genes_groups :
- il faut faire l'analyse sur les données transformées en log.
- scanpy.tl.rank_genes_groups(adata, groupby = 'cluster') : effectue la comparaison de chaque groupe déterminé par la valeur de 'cluster' dans les metadata contre le reste.
- attention, le champ dans groupby doit être de type string (entier pas possible, donc le convertir en string si besoin).
- si le champ raw a été défini, il est utilisé, à moins d'indiquer use_raw = False
- attention scanpy.tl.rank_genes_groups par défaut modifie l'objet adata sur place en rajoutant la clef "rank_genes_groups" dans le champ adata.uns (si elle existe déjà, elle est écrasée). Faire copy = True si on ne le veut pas, et alors, une copie de l'objet est renvoyée avec ce champ ".uns" rempli.
- adata.uns['rank_genes_groups'] est un dictionnaire avec les champs suivants :
- params : un dictionnaire avec les paramètres de la recherche.
- names : une recarray avec des records indexés par groupe et comme valeurs les noms des gènes.
- scores : recarray comme names, mais avec les z-score permettant le calcul de la p-value
- pvals : recarray comme names, mais avec les p-values.
- pvals_adj : recarray comme names, mais avec les p-values ajustées.
- logfoldchanges : recarray comme names, mais avec les log2 fold change
- pts : si demandé, dataframe avec les gènes en ligne et les groupes en colonne indiquant le pourcentage des cellules exprimant le gène dans le groupe considéré (pas dépendant d'une comparaison particulière).
- pts_rest : si demandé et si reference vaut 'rest' (le défaut), dataframe avec les gènes en ligne et les groupes en colonne indiquant le pourcentage des cellules exprimant le gène dans le reste du groupe considéré.
- groups = ['0', '1', '2'] : limite l'analyse à certains groupes, ici les numéros de clusters si groupby indique le champ cluster (le défaut est 'all')
- reference = '2' : indique que toutes les comparaisons sont à faire par rapport au groupe '2' (défaut est 'rest'). Si groups est fixé, on peut ou non mettre ce groupe référence dans groups, ça ne change rien.
- method = 't-test_overestim_var' : la méthode de comparaison utilisée. C'est le défaut. Les autres valeurs sont 't-test', 'wilcoxon' et 'logreg'.
- corr_method = 'benjamini-hochberg' : la méthode de correction de la p-value. C'est le défaut (autre valeur : 'bonferroni').
- pts = True : sort aussi les pourcentages de cellules qui expriment le gène (défaut est False)
- n_genes = 100: nombre de gènes à sortir (le défaut est None, c'est à dire tous).
- key_added = 'myComp' : la clef à rajouter à adata.uns. Le défaut est 'rank_genes_groups'.
- layer = 'myLayer' : utilise la layer indiquée (le défaut est None).
- use_raw = True : utilise adata.raw si existe (le défaut est None).
Pour récupérer le résultat d'une analyse d'expression différentielle des gènes sous forme d'un dataframe :
- df = scanpy.get.rank_genes_groups_df(adata, group = None) : retourne les résultats de scanpy.tl.rank_genes_groups_df sous forme d'un dataframe
- paramètres :
- si group est une liste, ne retourne que les comparaisons pour les groupes indiqué (défaut est None, tous les groupes)
- on peut donner key = 'myComp' : va chercher les résultats qui ont été générés avec scanpy.tl.rank_genes_groups en donnant key_added = 'myComp'.
- Attention : si un seul groupe dans la comparaison, le dataframe résultant n'aura pas de colonne group !
Calcul de la phase du cycle cellulaire :
- scanpy.tl.score_genes_cell_cycle(adata, s_genes = sGenes, g2m_genes = g2mGenes)
- doit se faire après normalisation et transformation log, sur tous les gènes mesurés.
- sGenes et g2mGenes (gènes spécifiques des phases S et G2M) doivent être inclus dans les gènes mesurés.
- rajoute à adata.obs les colonnes S_score, G2M_score et phase
Calcul d'un score pour une liste de gènes :
- scanpy.tl.score_genes(adata, geneList, score_name = 'myScore') : calcule un score sur les gènes indiqués pour chaque cellule et le met dans la colonne 'myScore'.
pseudobulk :
- adata2 = scanpy.get.aggregate(adata[:, [gene1, gene2]], by = 'sample', func = ['mean', 'var']) : renvoie un anndata de dimension nbr de sample x 2 (car ici 2 gènes)
- fonctions possibles : count_nonzero, mean, sum, var.
- l'objet résultant adata2 a autant de layers que de fonctions données, et on accède à la matrice avec adata2.layers['mean'] pour la moyenne.
- adata2.X est vide.
Récupération de données sous forme d'un 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).
Graphes
Plots :
- les plots sont sauvés par défaut dans un directory figures : pour éviter ça : scanpy.settings.figdir = ''. scanpy rajoute en plus un prefix au nom du fichier, par exemple umap ...
- une façon d'éviter ce comportement 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.umap :
- scanpy.pl.umap (adata, color = 'clusterNbr') : représentation d'une donnée de type catégorie. color indique le nom de la catégorie à utiliser (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 = 'YlOrRd') : 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.embedding(adata, basis="umap", frameon=False) : donne la même chose que scanpy.pl.umap
scanpy.pl.scatter :
- scanpy.pl.scatter(adata, x = 'total_counts', y = 'pct_counts_mito')
- scanpy.pl.scatter(adata, basis = 'umap', color = 'cellType')
- 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.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.
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.
Divers :
- scanpy.read_10x_h5('myFile.h5') : lit un fichier 10x au format h5.
- 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 (il faut avoir installé loompy).
- 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
