1. download data
## 2022-05-02 17:31:13 URL:https://cf.10xgenomics.com/samples/cell-arc/1.0.0/pbmc_granulocyte_sorted_10k/pbmc_granulocyte_sorted_10k_filtered_feature_bc_matrix.h5 [162282142/162282142] -> "pbmc_granulocyte_sorted_10k_filtered_feature_bc_matrix.h5.1" [1]
## 2022-05-02 17:31:15 URL:https://ucd327629ac5bfc194e3d8e5b76e.dl.dropboxusercontent.com/cd/0/inline/Bkjcu2D3KndG-jTuX0l0g8tj_zmnxZiRsamJLqt2-weTzd26q5E8z9iL4DgvPRdBLbfWAfl5294eK_nn_7yKIn5vA0NF8kYdzD1Z0ifNnjgA-vU1E-A03mFeATA4QHJrAe8XMg5HLSbwJdSO9FwOf0AFktkDsdc-f_Lutq98Ny5wHw/file [400178/400178] -> "PBMC-Multiom_annotation.tsv.1" [1]
3. Create Seurat Object
mtxs <- Read10X_h5("pbmc_granulocyte_sorted_10k_filtered_feature_bc_matrix.h5")
## Genome matrix has multiple modalities, returning a list of matrices for this genome
RNA <- mtxs[["Gene Expression"]]
ATAC <- mtxs[["Peaks"]]
peaks <- rownames(ATAC)
ATAC <- ATAC[which(startsWith(peaks, "chr")), ]
meta <- read.csv("PBMC-Multiom_annotation.tsv", sep="\t")
RNA <- RNA[, meta$barcode]
ATAC <- ATAC[, meta$barcode]
rownames(meta) <- meta$barcode
object <- CreateSeuratObject(counts=RNA, meta.data=meta, assay="RNA")
object[["Peaks"]] <- CreateAssayObject(counts=ATAC)
object$celltype <- object$annotation
rm(mtxs)
rm(RNA)
rm(ATAC)
gc()
## used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells 13004860 694.6 20968081 1119.9 NA 16555972 884.2
## Vcells 185197949 1413.0 540369956 4122.7 16384 570321102 4351.3
subset if need to save memory
4. RNA & ATAC dimension reductions
## RNA pre-processing and PCA dimension reduction
DefaultAssay(object) <- "RNA"
object <- NormalizeData(object, normalization.method = "LogNormalize", scale.factor = 10000, verbose=F)
object <- FindVariableFeatures(object, nfeatures=3000, verbose=F)
object <- ScaleData(object, verbose=F)
object <- RunPCA(object, npcs=50, reduction.name="RNA_PCA", verbose=F)
## RNA pre-processing and LSI dimension reduction
DefaultAssay(object) <- "Peaks"
object <- RunTFIDF(object, verbose=F)
## Warning in RunTFIDF.default(object = GetAssayData(object = object, slot =
## "counts"), : Some features contain 0 total counts
5. RUN MOJITOO
object <- mojitoo(
object=object,
reduction.list = list("RNA_PCA", "lsi"),
dims.list = list(1:50, 2:50), ## exclude 1st dimension of LSI
reduction.name='MOJITOO',
assay="RNA"
)
## processing RNA_PCA
## adding lsi
## 1 round cc 44
## Warning: Keys should be one or more alphanumeric characters followed by an
## underscore, setting key from MOJITOO to MOJITOO_
## Warning: All keys should be one or more alphanumeric characters followed by an
## underscore '_', setting key to MOJITOO_
DefaultAssay(object) <- "RNA"
embedd <- Embeddings(object[["MOJITOO"]])
object <- RunUMAP(object, reduction="MOJITOO", reduction.name="MOJITOO_UMAP", dims=1:ncol(embedd), verbose=F)
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
6. UMAP of True labels
object <- readRDS("seu.Rds") ## 1/3 cells
DimPlot(object, reduction = "MOJITOO_UMAP", group.by = "celltype", cols=ggsci::pal_igv()(30))
7. MOJITOO CCs
GeneCCDimPlot(object,
CCsToPlot = 1:4,
RNA.assay="RNA",
umap = "MOJITOO_UMAP",
MOJITOO.reduction="MOJITOO")
## adding MOJITOO assay...
## Warning: Keys should be one or more alphanumeric characters followed by an
## underscore, setting key from mojitoo.assay_ to mojitooassay_
## Centering and scaling data matrix
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
Heatmap
GeneCCHeatmap(object,
CCsToPlot = 1:2,
RNA.assay="RNA",
colorbar.group = "celltype",
MOJITOO.reduction="MOJITOO",
filter.mito = T,
filter.ribo = T,
topN = 10
)
## ccsybygenes...
## Centering and scaling data matrix
## 'magick' package is suggested to install to give better rasterization.
##
## Set `ht_opt$message = FALSE` to turn off this message.
## 'magick' package is suggested to install to give better rasterization.
##
## Set `ht_opt$message = FALSE` to turn off this message.
## $`1`
##
## $`2`
Download track dependent data
## 2022-05-02 17:33:13 URL:https://uc9676b5a9d2110aa87cb8ec746a.dl.dropboxusercontent.com/cd/0/inline2/BkjXeSMERgrBWvXdvXtygx2d4o5jS4XxsbKH5bgxKIXTVZWb2OsZfkhrsVRJ1WxT0IDwVWS0x3F-06IS3NVUyCImXs0jbYpW-8t692TicmYLFPUFsA5hrmjcwkCidaFbyiEx_1a0YFh7uMPxzewPOgGIkaci-tCjECLFU0iAjXmDuDQuudc8A7nfb_4LLnKB0Uq1ATmd_ZOrmKzcrcwGa8t_Is1NgYLO9AiZM3ho-GzSg834uPeD4ygW4P9zBFHI_JmgNmVqfOMtCb60n8EcflVfh3CIUkZZ5OdwLJ-BNLhJW5ohPPIZyJ4fIqSsmpcK9LY2EyW8R644ANSgSEjE9tKcbqGjTKJCPt2esWrWp5-gJiQFaJOWsq4ZPN2hCPmlYlwACC3_6jXt7L-fKYmRPhvtoFrOAlAfDhQSJp6xRKZR-g/file [24022019/24022019] -> "genes.gtf.zip" [1]
## Archive: genes.gtf.zip
## inflating: genes.gtf
## 2022-05-02 17:35:48 URL:https://uc594a3dd2b81eb581c3ec6dd20d.dl.dropboxusercontent.com/cd/0/inline2/BkjssFb0wBrcZkBUuNyCd7uyxfsd3Bnm-CtS-yGJ2-bdZnTYa4FJdbVMHIz4CLL7n_5hZs1ffCwGgLKW0aSYm4ra9huhRFZslFICCgXiFQm4VI3rjvED9hT0TFK2_IXtLkWqEpEp7KanEFpHVjQ4fyjiiMci-a19ulpk71h6S5bRg6LM--uq7LwHFAe4K-ePFvatSMwX6mZnp8UdjUijmPysYNkeh3DP6G8Iml7e9pvxN-Yp5FclK0zs-ohCDa8O8-8F3RHn-SkFTeLJ7qE5SE2o3a3UVy5pbdUG6-qoIujdVoZNdVwDmMo2QXa5PoPDiX-IucXG5cmaE1hmGqke8CzbMJZYYY-Jf2FLhCN2bfiM8pBhvre7ZtnUSIxHObwsBRUgQg5ol2HZUaWbHAASPU5078srMDbJQD3uMIDkziqiCQ/file [717556067/717556067] -> "bigwig.zip" [1]
## Archive: bigwig.zip
## creating: bigwig/
## inflating: bigwig/CD14+_Monocytes.bw
## inflating: bigwig/CD4_Memory.bw
## inflating: bigwig/Platelets.bw
## inflating: bigwig/CD4_Naive.bw
## inflating: bigwig/CD16+_Monocytes.bw
## inflating: bigwig/pre-B_cell.bw
## inflating: bigwig/CD8_Naive.bw
## inflating: bigwig/pDC.bw
## inflating: bigwig/Dendritic_cell.bw
## inflating: bigwig/NK_cell.bw
## inflating: bigwig/Double_negative_T_cell.bw
## inflating: bigwig/CD8_effector.bw
## inflating: bigwig/B_cell_progenitor.bw
Tracks
data_track_bws <- list(
"B cell progenitor" = "bigwig/B_cell_progenitor.bw",
"CD14+ Monocytes" = "bigwig/CD14+_Monocytes.bw",
"CD16+ Monocytes" = "bigwig/CD16+_Monocytes.bw",
"CD4 Memory" = "bigwig/CD4_Memory.bw",
"CD4 Naive" = "bigwig/CD4_Naive.bw",
"CD8 effector" = "bigwig/CD8_effector.bw",
"CD8 Naive" = "bigwig/CD8_Naive.bw",
"Dendritic cell" = "bigwig/Dendritic_cell.bw",
"Double negative T cell" = "bigwig/Double_negative_T_cell.bw",
"NK cell" = "bigwig/NK_cell.bw",
"pDC" = "bigwig/pDC.bw",
"Platelets" = "bigwig/Platelets.bw",
"pre-B cell" = "bigwig/pre-B_cell.bw"
)
suppressPackageStartupMessages(library(rtracklayer))
gene_model <- readGFF("genes.gtf")
gene_model$chromosome <- gene_model$seqid
gene_model$gene <- gene_model$gene_id
gene_model$transcript <- gene_model$transcript_id
gene_model$symbol <- gene_model$gene_name
gene_model$exon <- gene_model$exon_id
gene_model$width <- gene_model$end - gene_model$start + 1
gene_model$feature <- gene_model$transcript_type
gene_model <- subset(gene_model, !is.na(transcript) & !is.na(exon))
gtree <- ATACTrack(object,
CC = 1,
group.by="celltype",
bigwig.file.list=data_track_bws,
MOJITOO.reduction="MOJITOO",
Peak.assay="Peaks",
gene.model=gene_model,
cols =ggsci::pal_igv()(51),
ylim.datatrack=c(0,16),
fontsize.geneAxis=5,
fontsize.geneRegion=10,
fontsize.datatrack=8,
show.legend=T,
genome="hg38"
)
## No peaks input, will get top2 posi&nega peaks from CC1
## loading data track bigwig files
## column
## 1
## 2
## 3
## 4
SessionInfo
## R version 4.2.0 (2022-04-22)
## Platform: aarch64-apple-darwin21.3.0 (64-bit)
## Running under: macOS Monterey 12.3.1
##
## Matrix products: default
## BLAS: /opt/homebrew/Cellar/openblas/0.3.20/lib/libopenblasp-r0.3.20.dylib
## LAPACK: /opt/homebrew/Cellar/r/4.2.0/lib/R/lib/libRlapack.dylib
##
## locale:
## [1] C/UTF-8/C/C/C/C
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] rtracklayer_1.56.0 GenomicRanges_1.48.0 GenomeInfoDb_1.32.1
## [4] IRanges_2.30.0 S4Vectors_0.34.0 BiocGenerics_0.42.0
## [7] ggsci_2.9 MOJITOO_1.0 Signac_1.6.0
## [10] sp_1.4-7 SeuratObject_4.1.0 Seurat_4.1.1
##
## loaded via a namespace (and not attached):
## [1] rappdirs_0.3.3 SnowballC_0.7.0
## [3] scattermore_0.8 ragg_1.2.2
## [5] tidyr_1.2.0 ggplot2_3.3.5
## [7] bit64_4.0.5 knitr_1.39
## [9] irlba_2.3.5 DelayedArray_0.22.0
## [11] data.table_1.14.2 rpart_4.1.16
## [13] AnnotationFilter_1.20.0 KEGGREST_1.36.0
## [15] RCurl_1.98-1.6 doParallel_1.0.17
## [17] generics_0.1.2 GenomicFeatures_1.48.0
## [19] cowplot_1.1.1 RSQLite_2.2.13
## [21] RANN_2.6.1 future_1.25.0
## [23] bit_4.0.4 spatstat.data_2.2-0
## [25] xml2_1.3.3 httpuv_1.6.5
## [27] SummarizedExperiment_1.26.0 assertthat_0.2.1
## [29] xfun_0.30 hms_1.1.1
## [31] jquerylib_0.1.4 evaluate_0.15
## [33] promises_1.2.0.1 fansi_1.0.3
## [35] restfulr_0.0.13 progress_1.2.2
## [37] dbplyr_2.1.1 igraph_1.3.1
## [39] DBI_1.1.2 htmlwidgets_1.5.4
## [41] sparsesvd_0.2 spatstat.geom_2.4-0
## [43] purrr_0.3.4 ellipsis_0.3.2
## [45] RSpectra_0.16-1 ks_1.13.5
## [47] backports_1.4.1 dplyr_1.0.9
## [49] biomaRt_2.52.0 deldir_1.0-6
## [51] MatrixGenerics_1.8.0 vctrs_0.4.1
## [53] Biobase_2.56.0 here_1.0.1
## [55] ensembldb_2.20.0 ROCR_1.0-11
## [57] abind_1.4-5 cachem_1.0.6
## [59] ggforce_0.3.3 Gviz_1.40.0
## [61] BSgenome_1.64.0 progressr_0.10.0
## [63] checkmate_2.1.0 sctransform_0.3.3
## [65] GenomicAlignments_1.32.0 prettyunits_1.1.1
## [67] mclust_5.4.9 goftest_1.2-3
## [69] cluster_2.1.3 lazyeval_0.2.2
## [71] crayon_1.5.1 hdf5r_1.3.5
## [73] labeling_0.4.2 pkgconfig_2.0.3
## [75] slam_0.1-50 tweenr_1.0.2
## [77] ProtGenerics_1.28.0 nlme_3.1-157
## [79] nnet_7.3-17 rlang_1.0.2
## [81] globals_0.14.0 lifecycle_1.0.1
## [83] miniUI_0.1.1.1 filelock_1.0.2
## [85] BiocFileCache_2.4.0 dichromat_2.0-0.1
## [87] rprojroot_2.0.3 polyclip_1.10-0
## [89] matrixStats_0.62.0 lmtest_0.9-40
## [91] Matrix_1.4-1 ggseqlogo_0.1
## [93] Rhdf5lib_1.18.0 zoo_1.8-10
## [95] base64enc_0.1-3 ggridges_0.5.3
## [97] GlobalOptions_0.1.2 png_0.1-7
## [99] viridisLite_0.4.0 rjson_0.2.21
## [101] bitops_1.0-7 KernSmooth_2.23-20
## [103] rhdf5filters_1.8.0 Biostrings_2.64.0
## [105] blob_1.2.3 shape_1.4.6
## [107] stringr_1.4.0 parallelly_1.31.1
## [109] spatstat.random_2.2-0 jpeg_0.1-9
## [111] scales_1.2.0 memoise_2.0.1
## [113] magrittr_2.0.3 plyr_1.8.7
## [115] ica_1.0-2 zlibbioc_1.42.0
## [117] hdrcde_3.4 compiler_4.2.0
## [119] BiocIO_1.6.0 RColorBrewer_1.1-3
## [121] clue_0.3-60 fitdistrplus_1.1-8
## [123] Rsamtools_2.12.0 cli_3.3.0
## [125] XVector_0.36.0 listenv_0.8.0
## [127] patchwork_1.1.1 pbapply_1.5-0
## [129] htmlTable_2.4.0 Formula_1.2-4
## [131] MASS_7.3-57 mgcv_1.8-40
## [133] tidyselect_1.1.2 stringi_1.7.6
## [135] textshaping_0.3.6 highr_0.9
## [137] yaml_2.3.5 latticeExtra_0.6-29
## [139] ggrepel_0.9.1 grid_4.2.0
## [141] VariantAnnotation_1.42.0 sass_0.4.1
## [143] fastmatch_1.1-3 tools_4.2.0
## [145] future.apply_1.9.0 parallel_4.2.0
## [147] rstudioapi_0.13 circlize_0.4.14
## [149] foreign_0.8-82 foreach_1.5.2
## [151] lsa_0.73.2 gridExtra_2.3
## [153] farver_2.1.0 Rtsne_0.16
## [155] digest_0.6.29 rgeos_0.5-9
## [157] pracma_2.3.8 shiny_1.7.1
## [159] qlcMatrix_0.9.7 Rcpp_1.0.8.3
## [161] later_1.3.0 fda_6.0.3
## [163] RcppAnnoy_0.0.19 httr_1.4.2
## [165] AnnotationDbi_1.58.0 biovizBase_1.44.0
## [167] ComplexHeatmap_2.12.0 colorspace_2.0-3
## [169] rainbow_3.6 XML_3.99-0.9
## [171] fs_1.5.2 tensor_1.5
## [173] reticulate_1.24 splines_4.2.0
## [175] uwot_0.1.11 RcppRoll_0.3.0
## [177] spatstat.utils_2.3-0 pkgdown_2.0.3
## [179] ArchR_1.0.1 plotly_4.10.0
## [181] systemfonts_1.0.4 xtable_1.8-4
## [183] fds_1.8 jsonlite_1.8.0
## [185] corpcor_1.6.10 R6_2.5.1
## [187] Hmisc_4.7-0 ramify_0.3.3
## [189] pillar_1.7.0 htmltools_0.5.2
## [191] mime_0.12 glue_1.6.2
## [193] fastmap_1.1.0 BiocParallel_1.30.0
## [195] deSolve_1.32 codetools_0.2-18
## [197] pcaPP_2.0-1 mvtnorm_1.1-3
## [199] utf8_1.2.2 lattice_0.20-45
## [201] bslib_0.3.1 spatstat.sparse_2.1-1
## [203] tibble_3.1.6 curl_4.3.2
## [205] leiden_0.3.10 survival_3.3-1
## [207] rmarkdown_2.14 docopt_0.7.1
## [209] desc_1.4.1 munsell_0.5.0
## [211] GetoptLong_1.0.5 rhdf5_2.40.0
## [213] GenomeInfoDbData_1.2.8 iterators_1.0.14
## [215] reshape2_1.4.4 gtable_0.3.0
## [217] spatstat.core_2.4-2