SC-3-05.R

SC-3-05.R

# Install Packages
install.packages('scTenifoldNet')
BiocManager::install('monocle')
install.packages('reticulate')

library(reticulate)
reticulate::install_miniconda()
reticulate::conda_install(packages = 'scipy', envname = 'scipy')
reticulate::use_condaenv('scipy')

# Loading libraries
library(Seurat)
library(monocle)
library(MAST)
library(scTypeGSEA)
library(scTenifoldNet)
library(Matrix)


# Moving to downloads
setwd("../Downloads/")

# Downloading files
download.file("http://cf.10xgenomics.com/samples/cell-exp/3.0.0/pbmc_10k_v3/pbmc_10k_v3_filtered_feature_bc_matrix.tar.gz", destfile = "PBMC10K.tar.gz")
untar("PBMC10K.tar.gz", exdir = "10K")

# Reading the cells
PBMC <- Read10X("10K/filtered_feature_bc_matrix/")

# QC
source('https://raw.githubusercontent.com/dosorio/utilities/master/singleCell/scQC.R')
PBMC <- scQC(PBMC)

# Cell Type
A <- scTypeGSEA::assignCellType(PBMC, test.use = 'MAST', resolution = 0.1, logfc.threshold = 0.5)
B <- scTypeGSEA::labelCelltype(A$Seurat_obj, A$cluster_celltype)
PBMC <- B$obj
rm(A)
rm(B)

# UMAP 
PBMC <- RunUMAP(PBMC, dims = 1:20)
UMAPPlot(PBMC)

# Cell Types
PBMC <- Seurat::SplitObject(PBMC)
names(PBMC)

# T cells
BCells <- PBMC$B_cells

# Count Matrix
BCells <- GetAssayData(BCells, slot = 'counts')

# Filtering out genes in less than the 10% of the cells
dim(BCells)
BCells <- BCells[rowMeans(BCells != 0) > 0.1,]
dim(BCells)

# Pseudotime
fd <- data.frame('gene_short_name' = rownames(BCells))
rownames(fd) <- rownames(BCells)
fd <- new("AnnotatedDataFrame", data = fd)
cds <- newCellDataSet(as.matrix(BCells), featureData = fd, expressionFamily = negbinomial.size())
cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds)
cds <- reduceDimension(cds, reduction_method = "DDRTree", verbose = TRUE, maxIter = 10)
cds <- orderCells(cds)
plot_cell_trajectory(cds)

# Developmental process
colFun <- colorRampPalette(c('blue','red'))
colFun <- colFun(length(cds$Pseudotime))
plot(t(cds@reducedDimS)[order(cds$Pseudotime),], pch = 16, col = colFun, xlab = 'DDRTree 1', ylab = 'DDRTree 2')

# New Categories
cellCat <- ifelse(cds$Pseudotime > median(cds$Pseudotime), 'old', 'young')

# B Cells
BCells <- PBMC$B_cells
Idents(BCells) <- cellCat
BCells <- RunTSNE(BCells, features = 1:20)
TSNEPlot(BCells)
UMAPPlot(BCells)


# DE
DE <- FindMarkers(BCells, ident.1 = 'young', ident.2 = 'old', test.use = 'MAST')
head(DE)
# Go to: https://amp.pharm.mssm.edu/Enrichr/
writeLines(rownames(DE)[DE$p_val_adj < 0.05])

# DR
CountMatrix <- GetAssayData(BCells, slot = 'counts')
CountMatrix <- CountMatrix[order(rowMeans(CountMatrix), decreasing = TRUE) <= 1500,]

Y <- CountMatrix[,cellCat == 'young']
O <- CountMatrix[,cellCat == 'old']

DR <- scTenifoldNet(Y,O, qc_minLibSize = 0)
DR$diffRegulation