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 :
  • 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) : la normalisation se fait à la médiane du compte par cellule avant normalisation.
• 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 retenir les 2000 gènes les plus variables avec : n_top_genes = 2000
• 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.
  • use_highly_variable = True : utilise les gènes les plus variables (c'est le défaut s'ils ont été calculés précédemment).
  • 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 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).
• 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 = 'rank_genes_groups' : la clef à rajouter à adata.uns (on peut ainsi utiliser autre chose que ce défaut).
• 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 sous forme d'un dataframe, on peut faire :

dfList = []
for group in res.uns['rank_genes_groups']['names'].dtype.names:
    df = pandas.DataFrame({field:res.uns['rank_genes_groups'][field][group] \
        for field in ['names', 'pvals_adj', 'logfoldchanges']})
    df.insert(0, 'group', group)
    dfList.append(df)
df = pandas.concat(dfList)
  

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.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.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.