Pde10a for Joel March 2025

Human primary myoblast differentiation data

HSMM-Pde10a

HSMM-Sik1

HSMM-Fos

HSMM-Fosb

HSMM-Il6

HSMM-Junb

HSMM-Nr4a2

HSMM-Ppargc1a

HSMM-Ckm

HSMM-Myh4

Figure 1: Primary human myoblast differentiation data from Trapnell et al., 2014. Note: uses first-generation single-cell RNA-seq technology (Fluidigm) so very few cells.

Mouse notexin injury data

FACS-sorted

DeMicheliFACS-Fos

DeMicheliFACS-Il6

DeMicheliFACS-Junb

DeMicheliFACS-Nr4a2

DeMicheliFACS-Pde10a

DeMicheliFACS-Ppargc1a

DeMicheliFACS-Sik1

Figure 2: FACS-sorted (on calcein) cells from notexin-injured mouse skeletal muscle at multiple time points post-injury from De Micheli et al., 2020. Note there are no replicates for this dataset. White squares indicate no detection of this cell type/timepoint combination.

Not FACS-sorted

DeMicheli-Fos

DeMicheli-Il6

DeMicheli-Junb

DeMicheli-Nr4a2

DeMicheli-Pde10a

DeMicheli-Ppargc1a

DeMicheli-Sik1

DeMicheli-bulk

Figure 3: Non-FACS-sorted cells from the same study as Figure 2. Either 2 or 3 replicates for each time point. The final subplot here averages expression across all cell types for each timepoint. White squares indicate no detection of this cell type/timepoint combination.

DeMicheli-Fos-dots

DeMicheli-Il6-dots

DeMicheli-Junb-dots

DeMicheli-Nr4a2-dots

DeMicheli-Pde10a-dots

DeMicheli-Ppargc1a-dots

DeMicheli-Sik1-dots

Figure 4: Alternative way of viewing Figure 3 (same data).

Mouse cardiotoxin injury data

Shihuan-Fos

Shihuan-Il6

Shihuan-Junb

Shihuan-Nr4a2

Shihuan-Pde10a

Shihuan-Ppargc1a

Shihuan-Sik1

Figure 5: From Oprescu et al., 2020. Only a single replicate for each time point (uninjured, 0.5 dpi, 2 dpi, 3.5 dpi, 5 dpi, 10 dpi, 21 dpi).

Code/reproducibility

HSMM data

library(HSMMSingleCell)
library(SingleCellExperiment)
library(scater)

data(HSMM_sample_sheet); data(HSMM_expr_matrix); data(HSMM_gene_annotation)
sce <- SingleCellExperiment(list(logcounts=HSMM_expr_matrix))
colData(sce) <- S4Vectors::DataFrame(HSMM_sample_sheet)
rowData(sce) <- S4Vectors::DataFrame(HSMM_gene_annotation)
rownames(sce) <- rowData(sce)$gene_short_name

plotExpression(sce, features="PDE10A", x="Hours", color_by="Media")
plotExpression(sce, features="SIK1", x="Hours", color_by="Media")
plotExpression(sce, features="FOS", x="Hours", color_by="Media")
plotExpression(sce, features="FOSB", x="Hours", color_by="Media")
plotExpression(sce, features="IL6", x="Hours", color_by="Media")
plotExpression(sce, features="JUNB", x="Hours", color_by="Media")
plotExpression(sce, features="NR4A2", x="Hours", color_by="Media")
plotExpression(sce, features="PPARGC1A", x="Hours", color_by="Media")
plotExpression(sce, features="CKM", x="Hours", color_by="Media")
plotExpression(sce, features="MYH7", x="Hours", color_by="Media")

De Micheli data

# download De Micheli data
wget https://ftp.ncbi.nlm.nih.gov/geo/series/GSE143nnn/GSE143435/suppl/GSE143435%5FDeMicheli%5FD0%5FFACSatlas%5Fmetadata%2Etxt%2Egz
wget https://ftp.ncbi.nlm.nih.gov/geo/series/GSE143nnn/GSE143435/suppl/GSE143435%5FDeMicheli%5FD0%5FFACSatlas%5Fnormalizeddata%2Etxt%2Egz
wget https://ftp.ncbi.nlm.nih.gov/geo/series/GSE143nnn/GSE143435/suppl/GSE143435%5FDeMicheli%5FD2%5FFACSatlas%5Fmetadata%2Etxt%2Egz
wget https://ftp.ncbi.nlm.nih.gov/geo/series/GSE143nnn/GSE143435/suppl/GSE143435%5FDeMicheli%5FD2%5FFACSatlas%5Fnormalizeddata%2Etxt%2Egz
wget https://ftp.ncbi.nlm.nih.gov/geo/series/GSE143nnn/GSE143435/suppl/GSE143435%5FDeMicheli%5FD5%5FFACSatlas%5Fmetadata%2Etxt%2Egz
wget https://ftp.ncbi.nlm.nih.gov/geo/series/GSE143nnn/GSE143435/suppl/GSE143435%5FDeMicheli%5FD5%5FFACSatlas%5Fnormalizeddata%2Etxt%2Egz
wget https://ftp.ncbi.nlm.nih.gov/geo/series/GSE143nnn/GSE143435/suppl/GSE143435%5FDeMicheli%5FD7%5FFACSatlas%5Fmetadata%2Etxt%2Egz
wget https://ftp.ncbi.nlm.nih.gov/geo/series/GSE143nnn/GSE143435/suppl/GSE143435%5FDeMicheli%5FD7%5FFACSatlas%5Fnormalizeddata%2Etxt%2Egz

import scanpy as sc
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

# hard-coded but seems to work for all
def read_dataset(dset):
    df = pd.read_table(f"GSE143435_DeMicheli_{dset}_FACSatlas_normalizeddata.txt.gz")
    df = df.set_index('Unnamed: 0')
    adata = sc.AnnData(X=df.to_numpy().T)
    adata.obs_names = df.columns
    obs = pd.read_table(f"GSE143435_DeMicheli_{dset}_FACSatlas_metadata.txt.gz")
    obs = obs.set_index("Unnamed: 0")
    adata.obs = obs
    adata.var_names = df.index
    return adata

d0 = read_dataset("D0")
d2 = read_dataset("D2")
d5 = read_dataset("D5")
d7 = read_dataset("D7")

merged = sc.concat([d0, d2, d5, d7])
gene_list = ["Pde10a", "Fos", "Nr4a2", "Sik1", "Junb", "Ppargc1a", "Il6"]

adata_agg = merged[:, gene_list]
adata_agg2 = sc.get.aggregate(adata_agg, by=["sampleID", "cell_annotation"], func=["mean", "count_nonzero"])
adata_agg2.X = adata_agg2.layers["mean"]
df = adata_agg2.to_df()
df["sampleID"] = adata_agg2.obs["sampleID"]
df["cell_annotation"] = adata_agg2.obs["cell_annotation"]

for g in gene_list:
    pivoted = df.pivot(index="sampleID", columns="cell_annotation", values=g)
    fig, ax = plt.subplots()
    sns.heatmap(pivoted, linewidths=.5, cmap="viridis", square=True, linecolor="black", ax=ax)
    fig.savefig(f"heatmap_{g}.png")

adata_agg2 = sc.get.aggregate(adata_agg, by=["sampleID"], func=["mean"])
adata_agg2.X = adata_agg2.layers["mean"]
df = adata_agg2.to_df()

fig, ax = plt.subplots()
sns.heatmap(df, linewidths=.5, cmap="viridis", square=True, linecolor="black")
plt.savefig("heatmap_bulk_timepoints.png")

# get non-FACS-sorted De Micheli data
wget https://ftp.ncbi.nlm.nih.gov/geo/series/GSE143nnn/GSE143437/suppl/GSE143437%5FDeMicheli%5FMuSCatlas%5Fnormalizeddata%2Etxt%2Egz
wget https://ftp.ncbi.nlm.nih.gov/geo/series/GSE143nnn/GSE143437/suppl/GSE143437%5FDeMicheli%5FMuSCatlas%5Fmetadata%2Etxt%2Egz

import scanpy as sc
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

df = pd.read_table("GSE143437_DeMicheli_MuSCatlas_normalizeddata.txt.gz")
df.set_index('Unnamed: 0', inplace=True)

meta = pd.read_table("GSE143437_DeMicheli_MuSCatlas_metadata.txt.gz")
meta.set_index('Unnamed: 0', inplace=True)

adata = sc.AnnData(X=df.to_numpy().T)
adata.var_names = df.index
adata.obs = meta

gene_list = ["Pde10a", "Fos", "Nr4a2", "Sik1", "Junb", "Ppargc1a", "Il6"]
adata = adata[:, gene_list]

adata.write("de-micheli-bigdataset.h5ad", compression="gzip")

adata = sc.read_h5ad("de-micheli-bigdataset.h5ad")

adata_agg2 = sc.get.aggregate(adata, by=["sampleID", "cell_annotation"], func=["mean", "count_nonzero"])
adata_agg2.X = adata_agg2.layers["mean"]
df2 = adata_agg2.to_df()
df2["sampleID"] = adata_agg2.obs["sampleID"]
df2["cell_annotation"] = adata_agg2.obs["cell_annotation"]

for g in adata.var_names:
    pivoted = df2.pivot(index="sampleID", columns="cell_annotation", values=g)
    
    fig, ax = plt.subplots()
    sns.heatmap(pivoted, linewidths=.5, cmap="viridis", square=True, linecolor="black", ax=ax)
    fig.savefig(f"heatmap_bigdataset_{g}.png")

    pivoted["sample"] = pivoted.index.copy()
    pivoted["time"] = [0, 0, 0, 2, 2, 5, 5, 5, 7, 7]
    grid = sns.FacetGrid(data=pivoted.melt(id_vars=["time", "sample"]), col="cell_annotation", col_wrap=4, sharey=False, height=2)
    grid.map(sns.regplot, "time", "value")
    grid.set_titles(col_template="{col_name}")
    grid.figure.suptitle(f"{g}", y=1.00, fontsize="xx-large", fontweight="bold")
    grid.savefig(f"dotplot_grid_{g}.png")
    plt.close()

Shihuan data

wget https://ftp.ncbi.nlm.nih.gov/geo/series/GSE138nnn/GSE138826/suppl/GSE138826%5Fregen%5Fdata%2Erds%2Egz
gunzip GSE138826_regen_data.rds.gz

library(tidyverse)
library(SingleCellExperiment)
library(scater)

d <- readr::read_rds("GSE138826_regen_data.rds")
sce <- SingleCellExperiment(list(logcounts=d$SCT@data))
colData(sce) <- S4Vectors::DataFrame(d$SCT@varMetadata)
colData(sce)$timepoint <- factor(colData(sce)$timepoint, levels = c("non-inj", "0_5dpi", "2dpi", "3_5dpi", "5dpi", "10dpi", "21dpi"))

scater::plotExpression(sce, features = "Pde10a", x = "timepoint", color_by = "metacluster") + facet_wrap(~colData(sce)$metacluster, nrow = 4) + ggpubr::rotate_x_text(45) + ggtitle("Pde10a")
scater::plotExpression(sce, features = "Fos", x = "timepoint", color_by = "metacluster") + facet_wrap(~colData(sce)$metacluster, nrow = 4) + ggpubr::rotate_x_text(45) + ggtitle("Fos")
scater::plotExpression(sce, features = "Nr4a2", x = "timepoint", color_by = "metacluster") + facet_wrap(~colData(sce)$metacluster, nrow = 4) + ggpubr::rotate_x_text(45) + ggtitle("Nr4a2")
scater::plotExpression(sce, features = "Sik1", x = "timepoint", color_by = "metacluster") + facet_wrap(~colData(sce)$metacluster, nrow = 4) + ggpubr::rotate_x_text(45) + ggtitle("Sik1")
scater::plotExpression(sce, features = "Junb", x = "timepoint", color_by = "metacluster") + facet_wrap(~colData(sce)$metacluster, nrow = 4) + ggpubr::rotate_x_text(45) + ggtitle("Junb")
scater::plotExpression(sce, features = "Ppargc1a", x = "timepoint", color_by = "metacluster") + facet_wrap(~colData(sce)$metacluster, nrow = 4) + ggpubr::rotate_x_text(45) + ggtitle("Ppargc1a")
scater::plotExpression(sce, features = "Il6", x = "timepoint", color_by = "metacluster") + facet_wrap(~colData(sce)$metacluster, nrow = 4) + ggpubr::rotate_x_text(45) + ggtitle("Il6")