01c_seuratObject_sct_integrated_dataLaura.R

#################################################
## Project: Create Seurat Object for Laura's data, use SCT and integration
## Date: 09.06.2020
## Author: Nathalie
##################################################

# setwd("/net/PE1/raid1/LAURA/SC_ANALYSIS/markerGeneDefinition")
# 
# library(scrattch.io, lib.loc="pkg_r")
# library(Seurat, lib.loc="pkg_r")
# library(dplyr)
# library(stringr)

setwd("~/Documents/PostmortemBrain/analysis/markerGeneDefinition")

library(scrattch.io)
library(Seurat)
library(dplyr)
library(stringr)

options(future.globals.maxSize = 3000 * 1024^2)

### DATA LAURA #################

# Read data as Seurat object
# <- Read10X_h5("/net/PE1/raid1/LAURA/SC_ANALYSIS/aggr/outs/filtered_feature_bc_matrix.h5")
data_aggr <- Read10X_h5("dataLaura/filtered_feature_bc_matrix.h5")
data <- CreateSeuratObject(counts = data_aggr,
                           project = "data_aggr",
                           min.cells = 3,
                           min.features = 100)

# add donor to metadata
cov_donor <- str_sub(dimnames(data)[[2]], -1, -1)
data <- AddMetaData(data, 
                    metadata = cov_donor,
                    col.name = "donor")

# split Seurat object into objects for each donor
data.list <- SplitObject(data, split.by = "donor")
rm(data)

# perform SC Transformation on each of the donor Seurat objects
for (i in 1:length(data.list)){
  data.list[[i]][["percent.mt"]] <- PercentageFeatureSet(data.list[[i]], pattern = "^MT-")
  data.list[[i]] <- subset(data.list[[i]], subset = nFeature_RNA > 500 & nFeature_RNA < 7500 & percent.mt <= 15)
  data.list[[i]] <- SCTransform(data.list[[i]], vars.to.regress = c("percent.mt"), verbose = FALSE)
}

# Select features for integration
data.features <- SelectIntegrationFeatures(object.list = data.list, nfeatures = 3000)
data.list <- PrepSCTIntegration(object.list = data.list, anchor.features = data.features, 
                                verbose = FALSE)

# Identify anchors and integrate datasets
data.anchors <- FindIntegrationAnchors(object.list = data.list, normalization.method = "SCT", 
                                       anchor.features = data.features, verbose = FALSE)
data.integrated <- IntegrateData(anchorset = data.anchors, normalization.method = "SCT", 
                                 verbose = FALSE)
rm(data.list)

# Perform PCA linear dimensional reduction
data.integrated <- RunPCA(data.integrated, verbose = FALSE)
DimPlot(data.integrated, reduction = "pca")

# Determine the "dimensionality" of the data
ElbowPlot(data.integrated)

# Run non-linear dimensional reduction (UMAP/tSNE)
data.integrated <- RunUMAP(data.integrated, dims = 1:10)
DimPlot(data.integrated, reduction = "umap")

# Cluster the cells
data.integrated <- FindNeighbors(data.integrated, dims = 1:8)
data.integrated <- FindClusters(data.integrated, resolution = 0.3)
DimPlot(data.integrated, reduction = "umap")
DimPlot(data.integrated, reduction = "umap", group.by = "donor")

saveRDS(data.integrated, file = "dataLaura/data_object_sct_integrated.rds")