SC-7-25.R

source("https://raw.githubusercontent.com/dosorio/utilities/master/idConvert/hsa_SYMBOL2ENTREZ.R")
source("https://raw.githubusercontent.com/dosorio/utilities/master/enrichments/hsa_KEGG_ENTREZ.R")

setwd("../Downloads/")
# 1K Cells
download.file("http://cf.10xgenomics.com/samples/cell-exp/3.0.0/pbmc_1k_v3/pbmc_1k_v3_filtered_feature_bc_matrix.tar.gz", destfile = "PBMC1K.tar.gz")
# 10K Cells
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("PBMC1K.tar.gz", exdir = "1K")
untar("PBMC10K.tar.gz", exdir = "10K")

library(Seurat)
library(sva)

# Reading the cells
PBMC_1K <- Read10X("1K/filtered_feature_bc_matrix/")
PBMC_10K <- Read10X("10K/filtered_feature_bc_matrix/")

# Selecting 1000 random cells
PBMC_10K <- PBMC_10K[,sample(colnames(PBMC_10K), 1000)]

# Simulating a batch effect
PBMC_1K <- t(t(as.matrix(PBMC_1K))/apply(PBMC_1K, 2, sum)) * 1e6
PBMC_10K <- t(t(as.matrix(PBMC_10K))/apply(PBMC_10K, 2, sum)) * 0.5e6

# Renaming the cells
colnames(PBMC_1K) <- paste0("1K_",colnames(PBMC_1K))
colnames(PBMC_10K) <- paste0("10K_",colnames(PBMC_10K))

# QC
QC <- function(X){
  mtCounts <- apply(X[grepl("MT-", rownames(X)),],2,sum)
  allCounts <- apply(X,2,sum)
  mtRate <- mtCounts/allCounts
  X <- X[apply(X,1,sum) > 0,mtRate < 0.1]
  return(X)
}
PBMC_1K <- QC(PBMC_1K)
PBMC_10K <- QC(PBMC_10K)

# Intersection
sharedGenes <- intersect(rownames(PBMC_1K),rownames(PBMC_10K))
sharedGenes <- cbind(PBMC_1K[sharedGenes,], PBMC_10K[sharedGenes,])
par(mfrow=c(1,2), mar=c(3,3,1,1), mgp=c(1.5,0.5,0))
plot(sharedGenes["MALAT1",], main = "BEFORE CORRECTION")

# Batch correction
B <- unlist(lapply(strsplit(colnames(sharedGenes), "_"), function(X){X[1]}))
sharedGenes <- as.matrix(sharedGenes[apply(sharedGenes,1,var) > 0,])
sharedGenes <- ComBat(dat = sharedGenes, batch = B)
sharedGenes <- round(sharedGenes)
sharedGenes[sharedGenes < 0] <- 0
plot(sharedGenes["MALAT1",], main = "AFTER CORRECTION")

# PreProcessing
PP <- function(X){
  X <- CreateSeuratObject(X)
  X <- ScaleData(X)
  X <- FindVariableFeatures(X)
  return(X)
}
PBMC_1K <- PP(sharedGenes[,B == "1K"])
PBMC_10K <- PP(sharedGenes[,B == "10K"])

# List of HVG
HVG <- intersect(VariableFeatures(PBMC_1K), VariableFeatures(PBMC_10K))

# CCA
outCCA <- RunCCA(PBMC_1K, PBMC_10K)
DimPlot(outCCA, reduction = "cca")

# Clusters
PBMC_1K <- RunPCA(PBMC_1K, verbose = FALSE)
PBMC_1K <- RunTSNE(PBMC_1K)
PBMC_1K <- FindNeighbors(PBMC_1K)
PBMC_1K <- FindClusters(PBMC_1K, resolution = 0.1)
TSNEPlot(PBMC_1K)

# DE
DE <- FindMarkers(PBMC_1K, ident.1 = "0", ident.2 = "2", test.use = "MAST")
DE <- DE[DE$p_val_adj < 0.5,]
# DE_U <- DE[DE$avg_logFC > 0,]
# DE_D <- DE[DE$avg_logFC < 0,]

# # Enrichments
# enrich_U <- hsa_KEGG_ENTREZ(hsa_SYMBOL2ENTREZ(rownames(DE_U))[,2])
# enrich_D <- hsa_KEGG_ENTREZ(hsa_SYMBOL2ENTREZ(rownames(DE_D))[,2])