#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"
cav1high <- subset(x = subset, subset = CAV1 > 2)
cav1_label <- sample(x = "cav1high", size = ncol(x = cav1high), replace = TRUE)
cav1high$groups <- cav1_label
my.data = FetchData(object = cav1high, vars = c("orig.ident","CAV1", "groups"))
head(my.data,10)

cav1low <- subset(x = subset, subset = CAV1 < 2)
cav1_label <- sample(x = "cav1low", size = ncol(x = cav1low), replace = TRUE)
cav1low$groups <- cav1_label
my.data = FetchData(object = cav1low, vars = c("orig.ident", "CAV1", "groups"))
head(my.data,10)


muscle.merged <- merge(x = cav1high, y = cav1low)
my.data = FetchData(object = muscle.merged, vars = c("orig.ident","CAV1", "groups"))
tail(my.data,10)
head(my.data,10)

DefaultAssay(muscle.merged) <- "RNA"

#Plots
cc.markers <- c("CDKN1A", "CDKN1B", "CDKN1C", "TP53", "RBL2", "CCND1", "CCND3", "MKI67", "ANLN", "BIRC5", "CCNA2", "CCNB1", "CCNE2", "CGOL1")
quiescence.markers <- c("KLF4", "MYC", "NDRG2", "SPRY1", "HEY1", "EGR1", "HES1", "CD82", "DAG1", "CHRNA1", "MEF2C", "SOX8")
DotPlot(object = muscle.merged, features = rev(x = quiescence.markers), cols = c("blue", "red"), 
        dot.scale = 8, group.by = "groups") + RotatedAxis()


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

# Find differentially expressed features between CAVhigh and CAV1low cells
DefaultAssay(muscle.merged) <- "RNA"
CAV1.MAST.markers <- FindMarkers(muscle.merged, ident.1 = "cav1high", ident.2 = "cav1low", test.use = "MAST")
