
#Load object

#subset satellite cells
subset = subset(muscle1, idents = c("0", "1", "2", "3", "4", "5", "6", "7", "8", "10", "12", "15"))

DimPlot(object = subset, reduction = "umap")

DefaultAssay(subset) <- "RNA"

##subset data
hey1spry1high <- subset(x = muscle, subset = HEY1 > 1 & SPRY1 > 1)
hey1spry1high_label <- sample(x = "hey1spry1high", size = ncol(x = hey1spry1high), replace = TRUE)
hey1spry1high$groups <- hey1spry1high_label
my.data = FetchData(object = hey1spry1high, vars = c("orig.ident","seurat_clusters", "SPRY1", "HEY1", "groups"))
head(my.data,10)


hey1spry1low <- subset(x = muscle, subset = HEY1 < 1 & SPRY1 < 1)
hey1spry1low_label <- sample(x = "hey1spry1low", size = ncol(x = hey1spry1low), replace = TRUE)
hey1spry1low$groups <- hey1spry1low_label
my.data = FetchData(object = hey1spry1low, vars = c("orig.ident","seurat_clusters", "SPRY1", "HEY1", "groups"))
head(my.data,10)

## Merge back
muscle.merged <- merge(x = hey1spry1high, y = hey1spry1low)
my.data = FetchData(object = muscle.merged, vars = c("orig.ident","seurat_clusters","SPRY1","HEY1", "groups"))
tail(my.data,10)
head(my.data,10)

#check differentially expressed genes
VlnPlot(muscle.merged, features = c("SPRY1", "HEY1", "CAV1","HES1"), pt.size = 0.0, group.by = "groups",ncol = 4)
VlnPlot(muscle.merged, features = c("PAX7","MYF5", "MYOD1", "MYOG"), pt.size = 0.0, group.by = "groups", ncol = 4)
VlnPlot(muscle.merged, features = c("MYF6", "FN1", "CHRNA1", "DAG1"), pt.size = 0.0, group.by = "groups", ncol = 4)
VlnPlot(muscle.merged, features = c("DLK1", "COL1A2", "CHRDL2", "CDKN1C"), pt.size = 0.0, group.by = "groups", ncol = 4)


#save
saveRDS(muscle.merged, file = "SC_subset_hey1spry1_cutoff_1.rds")


# Set identity classes to an existing column in meta data
Idents(object = muscle.merged) <- "groups"
levels(muscle.merged)


#Differential Expression using MAST package
DefaultAssay(muscle.merged) <- "RNA"
hey1spry1.MAST.markers <- FindMarkers(muscle.merged, ident.1 = "hey1spry1high", ident.2 = "hey1spry1low", test.use = "MAST")
head(n=10, hey1spry1.MAST.markers)

write.csv(hey1spry1.MAST.markers, file=paste("MAST_SCsubset_Hey1Spry1.csv", sep=""))


