library(Seurat)


########Merged Seurat v3######

Vastus56M_steven.data <- Read10X(data.dir = "~/Matrices/56-vastus/hg19/")
vastus20M.data <- Read10X(data.dir = "~/Matrices/427-Vastus20M/filtered_feature_bc_matrix/")
vastus55M.data <- Read10X(data.dir = "~/Matrices/429-Vastus55M/filtered_feature_bc_matrix/")
vastus57F.data <- Read10X(data.dir = "~/Matrices/431-Vastus57F/filtered_feature_bc_matrix/")
vastus52M.data <- Read10X(data.dir = "~/Matrices/430-Vastus52M/filtered_feature_bc_matrix/")
vastus66M.data <- Read10X(data.dir = "~/Matrices/416-Vastus_rerun/filtered_feature_bc_matrix/")
vastus73M.data <- Read10X(data.dir = "~/Matrices/428-Vastus73M/filtered_feature_bc_matrix/")
vastus83M.data <- Read10X(data.dir = "~/Matrices/423-Vastus83M/filtered_feature_bc_matrix/")


# Set up control object
Vastus56M_steven <- CreateSeuratObject(counts = Vastus56M_steven.data, project = "Vastus56M_steven", min.cells = 5)
Vastus56M_steven$sample <- "Vastus56M_steven"
vastus20M <- CreateSeuratObject(counts = vastus20M.data, project = "vastus20M", min.cells = 5)
vastus20M$sample <- "vastus20M"
vastus55M <- CreateSeuratObject(counts = vastus55M.data, project = "vastus55M", min.cells = 5)
vastus55M$sample <- "vastus55M"
vastus57F <- CreateSeuratObject(counts = vastus57F.data, project = "vastus57F", min.cells = 5)
vastus57F$sample <- "vastus57F"
vastus52M <- CreateSeuratObject(counts = vastus52M.data, project = "vastus52M", min.cells = 5)
vastus52M$sample <- "vastus52M"
vastus66M <- CreateSeuratObject(counts = vastus66M.data, project = "vastus66M", min.cells = 5)
vastus66M$sample <- "vastus66M"
vastus73M <- CreateSeuratObject(counts = vastus73M.data, project = "vastus73M", min.cells = 5)
vastus73M$sample <- "vastus73M"
vastus83M <- CreateSeuratObject(counts = vastus83M.data, project = "vastus83M", min.cells = 5)
vastus83M$sample <- "vastus83M"

# The [[ operator can add columns to object metadata. This is a great place to stash QC stats
Vastus56M_steven[["percent.mt"]] <- PercentageFeatureSet(object = Vastus56M_steven, pattern = "^MT-")
vastus20M[["percent.mt"]] <- PercentageFeatureSet(object = vastus20M, pattern = "^MT-")
vastus55M[["percent.mt"]] <- PercentageFeatureSet(object = vastus55M, pattern = "^MT-")
vastus57F[["percent.mt"]] <- PercentageFeatureSet(object = vastus57F, pattern = "^MT-")
vastus52M[["percent.mt"]] <- PercentageFeatureSet(object = vastus52M, pattern = "^MT-")
vastus66M[["percent.mt"]] <- PercentageFeatureSet(object = vastus66M, pattern = "^MT-")
vastus73M[["percent.mt"]] <- PercentageFeatureSet(object = vastus73M, pattern = "^MT-")
vastus83M[["percent.mt"]] <- PercentageFeatureSet(object = vastus83M, pattern = "^MT-")


muscleMT <- merge(x = vastus20M, y = list(vastus52M, vastus55M, vastus57F, vastus66M, vastus73M, vastus83M, Vastus56M_steven))
VlnPlot(muscleMT, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, split.by = "sample")
plot1 <- FeatureScatter(muscleMT, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(muscleMT, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
CombinePlots(plots = list(plot1, plot2))

muscles <- subset(muscleMT, subset =  nFeature_RNA > 500 & nFeature_RNA < 6000 & percent.mt < 10)
muscles.list <- SplitObject(muscles, split.by = "sample")

for (i in 1:length(muscles.list)) {
  muscles.list[[i]] <- NormalizeData(muscles.list[[i]], verbose = FALSE)
  muscles.list[[i]] <- FindVariableFeatures(muscles.list[[i]], selection.method = "vst", 
                                            nfeatures = 2000, verbose = FALSE)
}

muscles.anchors <- FindIntegrationAnchors(object.list = muscles.list, dims = 1:30)
muscles.integrated <- IntegrateData(anchorset = muscles.anchors, dims = 1:30)

DefaultAssay(object = muscles.integrated) <- "integrated"

#Scale the data
all.genes2 <- rownames(muscles.integrated)
muscle <- ScaleData(muscles.integrated, features = all.genes2)

muscle -> muscle1

muscle1 <- ScaleData(muscles.integrated, vars.to.regress = c("nCount_RNA", "percent.mt"), verbose = FALSE)
muscle1 <- RunPCA(muscle1, npcs = 30, verbose = FALSE)
muscle1 <- RunUMAP(muscle1, reduction = "pca", dims = 1:30)


muscle1 <- FindNeighbors(object = muscle1, reduction = "pca", dims = 1:30)
muscle1 <- FindClusters(object = muscle1, resolution = 0.5)

DimPlot(object = muscle1, reduction = "umap", legend = TRUE, label = TRUE)

saveRDS(muscle1, file = "Vasti8.rds")


# Differential Expression
DefaultAssay(muscle1) <- "RNA"
muscle.markers <- FindAllMarkers(object = muscle1, only.pos = TRUE, min.pct = 0.25, logfc.threshold = log(0.25),test.use = "MAST" )
head(muscle.markers)

write.csv(muscle.markers, file=paste("cluster_markers_vasti8MAST.csv", sep=""))



