library(Seurat)SeuratObject = readRDS('SeuratObject.RDS')Use clusters 13 and 19 for astrocytes
result1 = NULL
for(c1 in c(13,19)){
cat('Do cluster', c1, '...\n')
r1 <- try(FindMarkers(SeuratObject, slot='data', ident.1 = "early-pathology", ident.2 = 'no-pathology', group.by = 'pathology.group', test.use = "MAST", subset.ident = c1, logfc.threshold = 0), silent = TRUE)
if(!inherits(r1, 'try-error')) result1 = rbind(result1, data.frame(Symbol = rownames(r1), Cluster = c1, Contrast = 'early-vs-nopathology', r1, stringsAsFactors = FALSE))
r1 <- try(FindMarkers(SeuratObject, slot='data', ident.1 = "late-pathology", ident.2 = 'no-pathology', group.by = 'pathology.group', test.use = "MAST", subset.ident = c1, logfc.threshold = 0), silent = TRUE)
if(!inherits(r1, 'try-error')) result1 = rbind(result1, data.frame(Symbol = rownames(r1), Cluster = c1, Contrast = 'late-vs-nopathology', r1, stringsAsFactors = FALSE))
}## Do cluster 13 ...##
## Done!## Combining coefficients and standard errors## Warning in melt(coefAndCI, as.is = TRUE): The melt generic in data.table has
## been passed a array and will attempt to redirect to the relevant reshape2
## method; please note that reshape2 is deprecated, and this redirection is now
## deprecated as well. To continue using melt methods from reshape2 while both
## libraries are attached, e.g. melt.list, you can prepend the namespace like
## reshape2::melt(coefAndCI). In the next version, this warning will become an
## error.## Calculating log-fold changes## Warning in melt(lfc): The melt generic in data.table has been passed a list
## and will attempt to redirect to the relevant reshape2 method; please note that
## reshape2 is deprecated, and this redirection is now deprecated as well. To
## continue using melt methods from reshape2 while both libraries are attached,
## e.g. melt.list, you can prepend the namespace like reshape2::melt(lfc). In the
## next version, this warning will become an error.## Calculating likelihood ratio tests## Refitting on reduced model...##
## Done!## Warning in melt(llrt): The melt generic in data.table has been passed a list
## and will attempt to redirect to the relevant reshape2 method; please note that
## reshape2 is deprecated, and this redirection is now deprecated as well. To
## continue using melt methods from reshape2 while both libraries are attached,
## e.g. melt.list, you can prepend the namespace like reshape2::melt(llrt). In the
## next version, this warning will become an error.##
## Done!## Combining coefficients and standard errors## Warning in melt(coefAndCI, as.is = TRUE): The melt generic in data.table has
## been passed a array and will attempt to redirect to the relevant reshape2
## method; please note that reshape2 is deprecated, and this redirection is now
## deprecated as well. To continue using melt methods from reshape2 while both
## libraries are attached, e.g. melt.list, you can prepend the namespace like
## reshape2::melt(coefAndCI). In the next version, this warning will become an
## error.## Calculating log-fold changes## Warning in melt(lfc): The melt generic in data.table has been passed a list
## and will attempt to redirect to the relevant reshape2 method; please note that
## reshape2 is deprecated, and this redirection is now deprecated as well. To
## continue using melt methods from reshape2 while both libraries are attached,
## e.g. melt.list, you can prepend the namespace like reshape2::melt(lfc). In the
## next version, this warning will become an error.## Calculating likelihood ratio tests## Refitting on reduced model...##
## Done!## Warning in melt(llrt): The melt generic in data.table has been passed a list
## and will attempt to redirect to the relevant reshape2 method; please note that
## reshape2 is deprecated, and this redirection is now deprecated as well. To
## continue using melt methods from reshape2 while both libraries are attached,
## e.g. melt.list, you can prepend the namespace like reshape2::melt(llrt). In the
## next version, this warning will become an error.## Do cluster 19 ...##
## Done!## Combining coefficients and standard errors## Warning in melt(coefAndCI, as.is = TRUE): The melt generic in data.table has
## been passed a array and will attempt to redirect to the relevant reshape2
## method; please note that reshape2 is deprecated, and this redirection is now
## deprecated as well. To continue using melt methods from reshape2 while both
## libraries are attached, e.g. melt.list, you can prepend the namespace like
## reshape2::melt(coefAndCI). In the next version, this warning will become an
## error.## Calculating log-fold changes## Warning in melt(lfc): The melt generic in data.table has been passed a list
## and will attempt to redirect to the relevant reshape2 method; please note that
## reshape2 is deprecated, and this redirection is now deprecated as well. To
## continue using melt methods from reshape2 while both libraries are attached,
## e.g. melt.list, you can prepend the namespace like reshape2::melt(lfc). In the
## next version, this warning will become an error.## Calculating likelihood ratio tests## Refitting on reduced model...##
## Done!## Warning in melt(llrt): The melt generic in data.table has been passed a list
## and will attempt to redirect to the relevant reshape2 method; please note that
## reshape2 is deprecated, and this redirection is now deprecated as well. To
## continue using melt methods from reshape2 while both libraries are attached,
## e.g. melt.list, you can prepend the namespace like reshape2::melt(llrt). In the
## next version, this warning will become an error.##
## Done!## Combining coefficients and standard errors## Warning in melt(coefAndCI, as.is = TRUE): The melt generic in data.table has
## been passed a array and will attempt to redirect to the relevant reshape2
## method; please note that reshape2 is deprecated, and this redirection is now
## deprecated as well. To continue using melt methods from reshape2 while both
## libraries are attached, e.g. melt.list, you can prepend the namespace like
## reshape2::melt(coefAndCI). In the next version, this warning will become an
## error.## Calculating log-fold changes## Warning in melt(lfc): The melt generic in data.table has been passed a list
## and will attempt to redirect to the relevant reshape2 method; please note that
## reshape2 is deprecated, and this redirection is now deprecated as well. To
## continue using melt methods from reshape2 while both libraries are attached,
## e.g. melt.list, you can prepend the namespace like reshape2::melt(lfc). In the
## next version, this warning will become an error.## Calculating likelihood ratio tests## Refitting on reduced model...##
## Done!## Warning in melt(llrt): The melt generic in data.table has been passed a list
## and will attempt to redirect to the relevant reshape2 method; please note that
## reshape2 is deprecated, and this redirection is now deprecated as well. To
## continue using melt methods from reshape2 while both libraries are attached,
## e.g. melt.list, you can prepend the namespace like reshape2::melt(llrt). In the
## next version, this warning will become an error.write.table(result1, file = 'cluster_astrocytes.disease_DEGs.MAST.tsv', row.names = FALSE, quote = FALSE, col.names = TRUE, sep = '\t')result2 = NULL
for(c1 in c(13,19)){
cat('Do cluster', c1, '...\n')
r1 <- try(FindMarkers(SeuratObject, slot='data', ident.1 = "early-pathology", ident.2 = 'no-pathology', group.by = 'pathology.group', test.use = "wilcox", subset.ident = c1, logfc.threshold = 0), silent = TRUE)
if(!inherits(r1, 'try-error')) result2 = rbind(result2, data.frame(Symbol = rownames(r1), Cluster = c1, Contrast = 'early-vs-nopathology', r1, stringsAsFactors = FALSE))
r1 <- try(FindMarkers(SeuratObject, slot='data', ident.1 = "late-pathology", ident.2 = 'no-pathology', group.by = 'pathology.group', test.use = "wilcox", subset.ident = c1, logfc.threshold = 0), silent = TRUE)
if(!inherits(r1, 'try-error')) result2 = rbind(result2, data.frame(Symbol = rownames(r1), Cluster = c1, Contrast = 'late-vs-nopathology', r1, stringsAsFactors = FALSE))
}## Do cluster 13 ...
## Do cluster 19 ...write.table(result2, file = 'cluster_astrocytes.disease_DEGs.wilcox.tsv', row.names = FALSE, quote = FALSE, col.names = TRUE, sep = '\t')library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionresult1 = result1 %>% mutate(Direction = ifelse(avg_logFC>0,"up","down"))
result2 = result2 %>% mutate(Direction = ifelse(avg_logFC>0,"up","down"))
result1_sum = result1 %>% group_by(Cluster, Direction) %>% summarize(MAST = n())## `summarise()` regrouping output by 'Cluster' (override with `.groups` argument)result2_sum = result2 %>% group_by(Cluster, Direction) %>% summarize(wilcox = n())## `summarise()` regrouping output by 'Cluster' (override with `.groups` argument)result_common = intersect(result1[,c("Symbol","Cluster","Direction")], result2[,c("Symbol","Cluster","Direction")]) %>% group_by(Cluster, Direction) %>% summarize(common = n())## `summarise()` regrouping output by 'Cluster' (override with `.groups` argument)result_combined = full_join(full_join(result1_sum, result2_sum), result_common)## Joining, by = c("Cluster", "Direction")## Joining, by = c("Cluster", "Direction")head(result_combined) ## Registered S3 method overwritten by 'cli':
## method from
## print.boxx spatstat## # A tibble: 4 x 5
## # Groups: Cluster [2]
## Cluster Direction MAST wilcox common
## <dbl> <chr> <int> <int> <int>
## 1 13 down 4529 4529 2878
## 2 13 up 3105 3105 2309
## 3 19 down 1397 1397 926
## 4 19 up 1461 1461 1124write.table(result_combined, file = 'cluster_astrocytes.disease_DEGs.MAST_wilcox_compare.tsv', row.names = FALSE, quote = FALSE, col.names = TRUE, sep = '\t')sessionInfo()## R version 3.6.3 (2020-02-29)
## Platform: x86_64-conda-linux-gnu (64-bit)
## Running under: CentOS Linux 7 (Core)
##
## Matrix products: default
## BLAS/LAPACK: /sc/arion/work/xup04/miniconda2_new/lib/libopenblasp-r0.3.12.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] dplyr_1.0.2 Seurat_3.2.3
##
## loaded via a namespace (and not attached):
## [1] Rtsne_0.15 colorspace_2.0-0
## [3] deldir_0.2-10 ellipsis_0.3.1
## [5] ggridges_0.5.2 XVector_0.26.0
## [7] GenomicRanges_1.38.0 spatstat.data_2.0-0
## [9] leiden_0.3.3 listenv_0.8.0
## [11] ggrepel_0.8.2 fansi_0.4.2
## [13] codetools_0.2-16 splines_3.6.3
## [15] knitr_1.30 polyclip_1.10-0
## [17] jsonlite_1.7.2 ica_1.0-2
## [19] cluster_2.1.1 png_0.1-7
## [21] uwot_0.1.10 shiny_1.6.0
## [23] sctransform_0.3.2 compiler_3.6.3
## [25] httr_1.4.2 Matrix_1.3-2
## [27] fastmap_1.1.0 lazyeval_0.2.2
## [29] cli_3.0.1 limma_3.42.2
## [31] later_1.1.0.1 prettyunits_1.1.1
## [33] htmltools_0.5.1.1 tools_3.6.3
## [35] rsvd_1.0.3 igraph_1.2.6
## [37] GenomeInfoDbData_1.2.2 gtable_0.3.0
## [39] glue_1.4.2 RANN_2.6.1
## [41] reshape2_1.4.4 rappdirs_0.3.3
## [43] Rcpp_1.0.6 spatstat_1.64-1
## [45] scattermore_0.7 Biobase_2.46.0
## [47] vctrs_0.3.6 nlme_3.1-152
## [49] lmtest_0.9-38 xfun_0.19
## [51] stringr_1.4.0 globals_0.14.0
## [53] mime_0.10 miniUI_0.1.1.1
## [55] lifecycle_1.0.0 irlba_2.3.3
## [57] goftest_1.2-2 future_1.21.0
## [59] zlibbioc_1.32.0 MASS_7.3-53
## [61] zoo_1.8-8 scales_1.1.1
## [63] MAST_1.12.0 hms_0.5.3
## [65] promises_1.2.0.1 spatstat.utils_2.0-0
## [67] parallel_3.6.3 SummarizedExperiment_1.16.1
## [69] RColorBrewer_1.1-2 SingleCellExperiment_1.8.0
## [71] reticulate_1.18 pbapply_1.4-3
## [73] gridExtra_2.3 ggplot2_3.3.2
## [75] rpart_4.1-15 stringi_1.5.3
## [77] S4Vectors_0.24.4 BiocGenerics_0.32.0
## [79] BiocParallel_1.20.1 GenomeInfoDb_1.22.1
## [81] bitops_1.0-6 rlang_0.4.10
## [83] pkgconfig_2.0.3 matrixStats_0.58.0
## [85] evaluate_0.14 lattice_0.20-41
## [87] ROCR_1.0-11 purrr_0.3.4
## [89] tensor_1.5 patchwork_1.0.1
## [91] htmlwidgets_1.5.3 cowplot_1.1.0
## [93] tidyselect_1.1.0 parallelly_1.23.0
## [95] RcppAnnoy_0.0.18 plyr_1.8.6
## [97] magrittr_2.0.1 R6_2.5.0
## [99] IRanges_2.20.2 generics_0.1.0
## [101] DelayedArray_0.12.3 pillar_1.5.0
## [103] mgcv_1.8-34 fitdistrplus_1.1-3
## [105] RCurl_1.98-1.2 survival_3.2-7
## [107] abind_1.4-5 tibble_3.1.0
## [109] future.apply_1.7.0 crayon_1.4.1
## [111] KernSmooth_2.23-18 utf8_1.1.4
## [113] plotly_4.9.2.1 rmarkdown_2.5
## [115] progress_1.2.2 grid_3.6.3
## [117] data.table_1.14.0 digest_0.6.27
## [119] xtable_1.8-4 tidyr_1.1.2
## [121] httpuv_1.5.5 stats4_3.6.3
## [123] munsell_0.5.0 viridisLite_0.3.0