01a_seuratObject_sct_allan.R

##################################################
## Project: Create Seurat object for Allan Brain Atlas
## Date: 06.04.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)

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

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

### HUMAN DATA #################

# Read tome object
tome <- "allan_human/transcrip.tome"
exons <- read_tome_dgCMatrix(tome, "data/t_exon")
introns <- read_tome_dgCMatrix(tome, "data/t_intron")
sample_name <- read_tome_sample_names(tome)
gene_name <- read_tome_gene_names(tome)

# Read sample annotations
anno <- read.table("allan_human/sample_annotations.csv",
                   header = TRUE,
                   sep = ",",
                   comment.char = "$",
                   row.names = 1)

# Create Seurat object
#norm_data <- logCPM(exons+introns)
counts <- exons+introns
colnames(counts) <- sample_name
rownames(counts) <- gene_name

data <- CreateSeuratObject(counts,
                           project = "allan_human",
                           assay = "RNA",
                           min.cells = 3,
                           min.features = 100)

data <- AddMetaData(
  object = data,
  metadata = anno
)
print(head(data@meta.data))

rm(counts)
rm(exons)
rm(introns)

# Visualize QC metrics as a violin plot
VlnPlot(data, features = c("nFeature_RNA", "nCount_RNA"), ncol = 2)
# Visualize relationship between nFeature und nCount
FeatureScatter(data, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")

# Remove cells with very few or too much features
data <- subset(data, subset = nFeature_RNA > 500 & nFeature_RNA < 7500)

# SC transformation (replaces NormalizeData, ScaleData and FindVariableFeatures)
data <- SCTransform(data, verbose = FALSE)

# Perform linear dimensional reduction
data <- RunPCA(data)
print(data[["pca"]], dims = 1:5, nfeatures = 5)
VizDimLoadings(data, dims = 1:2, reduction = "pca")
DimPlot(data, reduction = "pca")
DimHeatmap(data, dims = 1, cells = 500, balanced = TRUE)

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

# Cluster the cells
data <- FindNeighbors(data, dims = 1:8)
data <- FindClusters(data, resolution = 0.5)
head(Idents(data), 5)

# Run non-linear dimensional reduction (UMAP/tSNE)
data <- RunUMAP(data, dims = 1:8)
DimPlot(data, reduction = "umap") + NoLegend()
DimPlot(data, reduction = "umap", group.by = "subclass_label")

saveRDS(data, file = "allan_human/data_object_sct.rds")