Load libraries
library(SingleCellExperiment)
library(scran)
library(dplyr)
library(batchelor)
library(scater)
library(BiocSingular)
Initial settings
# Metadata info
KO_sample <- c(2, 5, 6)
Input_samples <- KO_sample
cond = rep('KO', 3)
# Input parameters
cluster_num <- 15
cut_off_ptprc <- 1
cut_off_Epcam <- 4
Setting up the data
all.sce <- list(
'sce.MS002' = list(),
'sce.MS005' = list(),
'sce.MS006' = list()
)
# Dataset
for (n in 1:length(all.sce)) {
all.sce[[n]] <-
read.table(
paste(
'./aomdss/rawdata_AOMDSS/counts/',
dir(path = './aomdss/rawdata_AOMDSS/counts/', pattern = 'TranscriptCounts')[Input_samples[n]],
sep = "/"
),
header = TRUE,
sep = "\t",
stringsAsFactors = FALSE
)
row.names(all.sce[[n]]) <- all.sce[[n]][, 1]
all.sce[[n]] <- as.matrix(all.sce[[n]][, -1])
colnames(all.sce[[n]]) <-
paste(sub("_counts_matrix", "", names(all.sce)[n]), colnames(all.sce[[n]]), sep = "_")
# Extract ERCC spike from counts_matrix
MS00X_spike_counts <-
all.sce[[n]][grep('ERCC-', rownames(all.sce[[n]])),]
MS00X_spike_se <-
SummarizedExperiment(list(counts = MS00X_spike_counts))
all.sce[[n]] <-
all.sce[[n]][-grep('ERCC-', rownames(all.sce[[n]])),]
# Metadata
MS00X_cell_metadata <-
data.frame(sample = c(rep(
sub("_counts_matrix", "", names(all.sce)[n]), ncol(all.sce[[n]])
)), condition = c(rep(cond[n], (ncol(
all.sce[[n]]
)))))
# Single cell experiment object
all.sce[[n]] <-
SingleCellExperiment(assays = list(counts = all.sce[[n]]),
colData = MS00X_cell_metadata)
altExp(all.sce[[n]], "ERCC") <- MS00X_spike_se
}
Quality control and normalization
stats <-
high.mito <-
high.Ptprc <-
high.Epcam <-
high.ercc <-
low.lib <-
low.nexprs <- list()
for (n in names(all.sce)) {
current <- all.sce[[n]]
is.mito <- grep("_chrM", rownames(current))
is.Ptprc <- grep('Ptprc__*', rownames(current))
is.Epcam <- grep('Epcam__*', rownames(current))
stats[[n]] <-
perCellQCMetrics(current,
subsets = list(
Mito = is.mito,
Ptprc = is.Ptprc,
Epcam = is.Epcam
))
high.mito[[n]] <-
isOutlier(stats[[n]]$subsets_Mito_percent, type = "higher")
high.Ptprc[[n]] <- stats[[n]]$subsets_Ptprc_sum > cut_off_ptprc
high.Epcam[[n]] <- stats[[n]]$subsets_Epcam_sum > cut_off_Epcam
high.ercc[[n]] <-
isOutlier(stats[[n]]$altexps_ERCC_percent, type = "higher")
low.lib[[n]] <-
isOutlier(stats[[n]]$sum, log = TRUE, type = "lower") # Log transformation
low.nexprs[[n]] <-
isOutlier(stats[[n]]$detected, log = TRUE, type = "lower") # Log transformation
all.sce[[n]] <-
current[, !(high.mito[[n]] | high.Ptprc[[n]] |high.Epcam[[n]])]
all.sce[[n]] <-
addPerCellQC(all.sce[[n]],
subsets = list(
Mito = is.mito,
Ptprc = is.Ptprc,
Epcam = is.Epcam
))
}
# Normalization
all.sce <- lapply(all.sce, logNormCounts, pseudo_count = 1)
Integrating datasets
# Create common features
common_features <- Reduce(intersect, (lapply(all.sce, rownames)))
common_features_ercc <-
Reduce(intersect, (lapply(lapply(all.sce, altExp), rownames)))
# Subsetting SingleCellExperiment
all.sce <- lapply(all.sce, function(x)
x[common_features])
for (sce_names in names(all.sce)) {
altExp(all.sce[[sce_names]]) <- altExp(all.sce[[sce_names]])[common_features_ercc]
}
combined_sce = cbind(all.sce$sce.MS002,
all.sce$sce.MS005,
all.sce$sce.MS006)
outlier_cells_samp_005 <-
as.vector(unlist(
read.csv2(
'./aomdss/rawdata_AOMDSS/metadata/out_cells_samp_005.txt',
header = FALSE
)
))
combined_sce <-
combined_sce[, -which(colnames(combined_sce) %in% outlier_cells_samp_005)]
# Feature selection
dec.block <- modelGeneVar(combined_sce, block = combined_sce$sample)
chosen.hvgs.block <- dec.block$bio > 0
# Batch correction
set.seed(01001001)
sce.data <- correctExperiments(
combined_sce,
batch = combined_sce$sample,
subset.row = chosen.hvgs.block,
correct.all = TRUE,
PARAM = FastMnnParam(
d = 50,
k = 20,
BSPARAM = RandomParam(deferred = TRUE)
)
)
Clustering
graph_snn <-
buildSNNGraph(sce.data,
k = cluster_num,
use.dimred = 'corrected')
clusters.mnn <-
igraph::cluster_walktrap(graph_snn)$membership
colLabels(sce.data) <-
factor(clusters.mnn)
set.seed(0010101010)
sce.data <- runTSNE(sce.data, dimred = "corrected")
plotReducedDim(sce.data, "TSNE", colour_by = "label")
Marker gene
markers <-
findMarkers(sce.data, direction = "up", lfc = 1)
markers_all <-
findMarkers(sce.data, direction = "any", lfc = 1)
Save results
save(sce.data, file = './aomdss/out_AOMDSS/KO/sce.data.RData')
print(sessioninfo::session_info())
## - Session info ---------------------------------------------------------------
## setting value
## version R version 4.0.4 (2021-02-15)
## os Windows 10 x64
## system x86_64, mingw32
## ui RTerm
## language (EN)
## collate German_Germany.1252
## ctype German_Germany.1252
## tz Europe/Berlin
## date 2024-01-24
##
## - Packages -------------------------------------------------------------------
## package * version date lib source
## assertthat 0.2.1 2019-03-21 [1] CRAN (R 4.0.0)
## batchelor * 1.6.2 2020-11-26 [1] Bioconductor
## beachmat 2.6.4 2020-12-20 [1] Bioconductor
## beeswarm 0.3.1 2021-03-07 [1] CRAN (R 4.0.4)
## Biobase * 2.50.0 2020-10-27 [1] Bioconductor
## BiocGenerics * 0.36.0 2020-10-27 [1] Bioconductor
## BiocNeighbors 1.8.2 2020-12-07 [1] Bioconductor
## BiocParallel 1.24.1 2020-11-06 [1] Bioconductor
## BiocSingular * 1.6.0 2020-10-27 [1] Bioconductor
## bitops 1.0-6 2013-08-17 [1] CRAN (R 4.0.0)
## bluster 1.0.0 2020-10-27 [1] Bioconductor
## cli 3.2.0 2022-02-14 [1] CRAN (R 4.0.5)
## colorspace 2.0-0 2020-11-11 [1] CRAN (R 4.0.3)
## cowplot 1.1.1 2020-12-30 [1] CRAN (R 4.0.3)
## crayon 1.4.1 2021-02-08 [1] CRAN (R 4.0.3)
## DBI 1.1.1 2021-01-15 [1] CRAN (R 4.0.3)
## DelayedArray 0.16.1 2021-01-22 [1] Bioconductor
## DelayedMatrixStats 1.12.3 2021-02-09 [1] Bioconductor
## digest 0.6.27 2020-10-24 [1] CRAN (R 4.0.3)
## dplyr * 1.0.4 2021-02-02 [1] CRAN (R 4.0.3)
## dqrng 0.2.1 2019-05-17 [1] CRAN (R 4.0.2)
## edgeR 3.32.1 2021-01-14 [1] Bioconductor
## ellipsis 0.3.1 2020-05-15 [1] CRAN (R 4.0.0)
## evaluate 0.14 2019-05-28 [1] CRAN (R 4.0.0)
## fansi 0.4.2 2021-01-15 [1] CRAN (R 4.0.3)
## farver 2.1.0 2021-02-28 [1] CRAN (R 4.0.5)
## fastmap 1.1.0 2021-01-25 [1] CRAN (R 4.0.3)
## generics 0.1.2 2022-01-31 [1] CRAN (R 4.0.5)
## GenomeInfoDb * 1.26.2 2020-12-08 [1] Bioconductor
## GenomeInfoDbData 1.2.4 2021-02-18 [1] Bioconductor
## GenomicRanges * 1.42.0 2020-10-27 [1] Bioconductor
## ggbeeswarm 0.6.0 2017-08-07 [1] CRAN (R 4.0.2)
## ggplot2 * 3.3.3 2020-12-30 [1] CRAN (R 4.0.3)
## glue 1.4.2 2020-08-27 [1] CRAN (R 4.0.2)
## gridExtra 2.3 2017-09-09 [1] CRAN (R 4.0.2)
## gtable 0.3.0 2019-03-25 [1] CRAN (R 4.0.0)
## highr 0.8 2019-03-20 [1] CRAN (R 4.0.0)
## htmltools 0.5.7 2023-11-03 [1] CRAN (R 4.0.4)
## igraph 1.2.6 2020-10-06 [1] CRAN (R 4.0.3)
## IRanges * 2.24.1 2020-12-12 [1] Bioconductor
## irlba 2.3.3 2019-02-05 [1] CRAN (R 4.0.2)
## jquerylib 0.1.3 2020-12-17 [1] CRAN (R 4.0.3)
## knitr 1.31 2021-01-27 [1] CRAN (R 4.0.3)
## labeling 0.4.2 2020-10-20 [1] CRAN (R 4.0.3)
## lattice 0.20-41 2020-04-02 [1] CRAN (R 4.0.4)
## lifecycle 1.0.0 2021-02-15 [1] CRAN (R 4.0.4)
## limma 3.46.0 2020-10-27 [1] Bioconductor
## locfit 1.5-9.4 2020-03-25 [1] CRAN (R 4.0.2)
## magrittr 2.0.1 2020-11-17 [1] CRAN (R 4.0.3)
## Matrix 1.3-2 2021-01-06 [1] CRAN (R 4.0.4)
## MatrixGenerics * 1.2.1 2021-01-30 [1] Bioconductor
## matrixStats * 0.58.0 2021-01-29 [1] CRAN (R 4.0.3)
## munsell 0.5.0 2018-06-12 [1] CRAN (R 4.0.0)
## pillar 1.5.1 2021-03-05 [1] CRAN (R 4.0.5)
## pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.0.0)
## purrr 0.3.4 2020-04-17 [1] CRAN (R 4.0.0)
## R6 2.5.0 2020-10-28 [1] CRAN (R 4.0.3)
## Rcpp 1.0.6 2021-01-15 [1] CRAN (R 4.0.3)
## RCurl 1.98-1.2 2020-04-18 [1] CRAN (R 4.0.0)
## ResidualMatrix 1.0.0 2020-10-27 [1] Bioconductor
## rlang 1.1.3 2024-01-10 [1] CRAN (R 4.0.4)
## rmarkdown 2.25 2023-09-18 [1] CRAN (R 4.0.4)
## rstudioapi 0.13 2020-11-12 [1] CRAN (R 4.0.3)
## rsvd 1.0.3 2020-02-17 [1] CRAN (R 4.0.2)
## Rtsne 0.15 2018-11-10 [1] CRAN (R 4.0.2)
## S4Vectors * 0.28.1 2020-12-09 [1] Bioconductor
## scales 1.1.1 2020-05-11 [1] CRAN (R 4.0.5)
## scater * 1.18.5 2021-02-16 [1] Bioconductor
## scran * 1.18.5 2021-02-04 [1] Bioconductor
## scuttle 1.0.4 2020-12-17 [1] Bioconductor
## sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 4.0.0)
## SingleCellExperiment * 1.12.0 2020-10-28 [1] Bioconductor
## sparseMatrixStats 1.2.1 2021-02-02 [1] Bioconductor
## statmod 1.4.35 2020-10-19 [1] CRAN (R 4.0.3)
## stringi 1.5.3 2020-09-09 [1] CRAN (R 4.0.2)
## stringr 1.4.0 2019-02-10 [1] CRAN (R 4.0.0)
## SummarizedExperiment * 1.20.0 2020-10-28 [1] Bioconductor
## tibble 3.0.6 2021-01-29 [1] CRAN (R 4.0.3)
## tidyselect 1.1.0 2020-05-11 [1] CRAN (R 4.0.2)
## utf8 1.2.1 2021-03-12 [1] CRAN (R 4.0.5)
## vctrs 0.3.6 2020-12-17 [1] CRAN (R 4.0.3)
## vipor 0.4.5 2017-03-22 [1] CRAN (R 4.0.2)
## viridis 0.5.1 2018-03-29 [1] CRAN (R 4.0.2)
## viridisLite 0.3.0 2018-02-01 [1] CRAN (R 4.0.0)
## withr 2.5.0 2022-03-03 [1] CRAN (R 4.0.5)
## xfun 0.41 2023-11-01 [1] CRAN (R 4.0.4)
## XVector 0.30.0 2020-10-28 [1] Bioconductor
## yaml 2.2.1 2020-02-01 [1] CRAN (R 4.0.0)
## zlibbioc 1.36.0 2020-10-28 [1] Bioconductor
##
## [1] C:/Users/guptapa/R-4.0.4/library