{"metadata":{"kernelspec":{"name":"ir","display_name":"R","language":"R"},"language_info":{"name":"R","codemirror_mode":"r","pygments_lexer":"r","mimetype":"text/x-r-source","file_extension":".r","version":"4.0.5"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"

PCA coloured by CD45-related genes and proteins that may be involved in post-translational/transcriptional regulation.

","metadata":{}},{"cell_type":"markdown","source":"

What we do here:

\n
\n \n
\n\nUsefool tools: \nAnimate ggplots with gganimate: CHEAT SHEET \nBecoming an AnimatoR","metadata":{}},{"cell_type":"markdown","source":"
\n\n

Data

\n\n
\n\nRDS files with sparse matrices of raw data and normalised counts data for Open Problems - Multimodal Single-Cell Integration (CITEseq2022). \n \nThe dataset for this competition comprises single-cell multiomics data (n = 70988 cells) collected from mobilized peripheral CD34+ hematopoietic stem and progenitor cells (HSPCs) isolated from healthy human donors.\nIn total, the CITEseq2022 train dataset we use here, contains 140 surface proteins and 22050 RNA form 3 days, 3 donors, 7 cell types. The gene names are given by {EnsemblID}_{GeneName} where EnsemblID refers to the Ensembl Gene ID and GeneName to the gene name (e.g. \"ENSG00000159840_ZYX\").\n \n
\n

Output

\n \n
\n
    \n \n
  1. 3D plots for CD45 isoforms, CD36, their RNAs, and some interesting RNAs
    2. \n
  2. GIFs with CD45 isoforms, CD36, their RNAs, and some interesting RNAs - change of expression by day
    3. \n
\n
\n
","metadata":{}},{"cell_type":"code","source":"suppressPackageStartupMessages({\n\n library(tictoc)\n library(tidyverse)\n library(Matrix)\n #devtools::install_github(\"satijalab/seurat\", ref = \"develop\")\n library(Seurat)\n \n # visualisation packages\n library(ggrepel)\n library(ggplot2)\n library(patchwork)\n library(ggpubr) #ggscatter \n library(pheatmap) #heatmap\n library(ggExtra) #ggMarginal\n library(gridExtra) # side by side plots\n library(plotly) #3D plot\n library(IRdisplay) #display 3Dplots\n library(gganimate) # animated plots\n library(gifski)# animated plots to gif\n \n install.packages(\"BiocManager\")\n BiocManager::install(\"glmGamPoi\")\n})\n\nfig <- function(width, heigth){\n options(repr.plot.width = width, repr.plot.height = heigth)\n}\n\n#fix colors range to have same colors on each heatmap\nbreaksList = seq(-1, 1, by = 0.01)\nmyColors <- c(colorRampPalette(c(\"darkblue\",\"white\"))(100), colorRampPalette(c(\"white\",\"red\"))(100))\n\n#We will save 3D plots in \"3D_plots\" folder and animations in \"GIFs\" folder\n\npath_3d <- \"3D_plots/\"\npath_gif <- \"GIFs/\"\n\ndir.create(file.path(path_3d), showWarnings = FALSE)\ndir.create(file.path(path_gif), showWarnings = FALSE)","metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","_kg_hide-output":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T04:49:49.432942Z","iopub.execute_input":"2023-04-16T04:49:49.434977Z","iopub.status.idle":"2023-04-16T04:58:43.602363Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Load data and create the Seurat object**","metadata":{}},{"cell_type":"code","source":"# randomly select x% of cells to save memory footprint\nset.seed(1)\nselected_prop <- 0.4\nselected_cells <- sample(1:70988,\n floor(selected_prop * 70988),\n replace = FALSE)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T04:58:43.606084Z","iopub.execute_input":"2023-04-16T04:58:43.636790Z","iopub.status.idle":"2023-04-16T04:58:43.656121Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Load raw RNA data\npath <- \"/kaggle/input/sparse-raw-counts-data-open-problems-multimodal/citeseq/sp_train_cite_inputs_raw.rds\"\nmat_RNA <- readRDS(path)\nmat_RNA <- mat_RNA[selected_cells, ]\nmat_RNA <- t(mat_RNA)\n\ncat(\"The dataset:\", nrow(mat_RNA), \"rows with RNA names and\", ncol(mat_RNA), \"columns with cell IDs\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T04:58:43.659780Z","iopub.execute_input":"2023-04-16T04:58:43.661526Z","iopub.status.idle":"2023-04-16T04:59:45.291419Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Additional cell-level metadata are stored in a data.frame where the rows are cell names and the columns are metadata fields. Row names in the metadata have to match the column names of the counts matrix.","metadata":{}},{"cell_type":"code","source":"# Load and filter metadata\nmetadata <- read.csv('../input/open-problems-multimodal/metadata.csv',row.names = 1)\n#summary(mutate_all(metadata, as.factor))\n\n# Filter citeseq, donor, day\nmetadata <- filter(metadata, technology == \"citeseq\" &\n donor %in% c(\"13176\", \"31800\", \"32606\") &\n day %in% c(\"2\",\"3\",\"4\"))\nmetadata <- filter(metadata, row.names(metadata) %in% colnames(mat_RNA)) #use xx% of cells for explaratory analysis\ncat(nrow(metadata), \"rows with cell IDs\")\nsummary(mutate_all(metadata, as.factor))\n#head(metadata)","metadata":{"_kg_hide-output":false,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T04:59:45.295072Z","iopub.execute_input":"2023-04-16T04:59:45.296759Z","iopub.status.idle":"2023-04-16T04:59:46.233500Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Cell types in this dataset: \nMasP = Mast Cell Progenitor \nMkP = Megakaryocyte Progenitor \nNeuP = Neutrophil Progenitor \nMoP = Monocyte Progenitor \nEryP = Erythrocyte Progenitor \nHSC = Hematoploetic Stem Cell \nBP = B-Cell Progenitor ","metadata":{}},{"cell_type":"markdown","source":"Now we create the Seurat object with metadata. As in [this notebook](https://www.kaggle.com/code/reminho/basic-scrnaseq-tutorial) we filtered for genes that are found in at least 100 cells (min.cells = 100) and a single cell should express at least 100 of our target genes (min.features = 100).","metadata":{}},{"cell_type":"code","source":"# Initialize the Seurat object with the raw (non-normalized data)\ntic()\nseur_RNA <- CreateSeuratObject(counts = mat_RNA, project = \"MmSCel\",\n min.cells = 100, min.features = 100, \n meta.data = metadata)\n\n# calculate mitochondrial percentage and asign this to the \"meta.data\" dataframe\nseur_RNA[[\"percent.mt\"]] <- PercentageFeatureSet(seur_RNA, pattern = \"MT-\")\nseur_RNA[[\"percent.rb\"]] <- PercentageFeatureSet(seur_RNA, pattern = \"-RP\")\n\n\ntoc()\nseur_RNA\ncat(\"meta data:\\n\")\nstr(seur_RNA@meta.data)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T04:59:46.235798Z","iopub.execute_input":"2023-04-16T04:59:46.237159Z","iopub.status.idle":"2023-04-16T05:00:16.519584Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rm(mat_RNA, metadata,path,selected_prop) ; gc()\nfor (thing in ls()) { message(thing); print(object.size(get(thing)), units='auto') }","metadata":{"_kg_hide-output":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:00:16.521800Z","iopub.execute_input":"2023-04-16T05:00:16.523112Z","iopub.status.idle":"2023-04-16T05:00:17.390928Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Normalizing the data**","metadata":{}},{"cell_type":"markdown","source":"As discussed [here](http://www.kaggle.com/datasets/alexandervc/research-project-01-around-multimodal-singlecell/discussion/378460), Pearson residuals from negative binomial regression, recommend e.g. by Hafemeister and Satija (Genome Biol 20:296, 2019, [vignette](http://satijalab.org/seurat/articles/sctransform_vignette.html)) outperform other methods in application for common downstream analytical tasks such as variable gene selection, dimensional reduction, and differential expression. We also compared different normalization techniques in [this notebook](http://www.kaggle.com/code/antoninadolgorukova/mmscel-raw-data-normalization-with-r).\n\nFor this we use SCTransform function from [sctransform packadge](http://github.com/saketkc/sctransform) (the last update - v2 regularization, vignette is [here](http://satijalab.org/seurat/articles/sctransform_v2_vignette.html#identify-differential-expressed-genes-across-conditions-1), accompanying article - [here](http://www.biorxiv.org/content/10.1101/2021.07.07.451498v1). Transformed data will be available in the SCT assay, which is set as the default after running sctransform.","metadata":{}},{"cell_type":"code","source":"tic() # ~3 min or 3.2 hours if conserve.memory = TRUE on 40% of cells\nseur_RNA <- SCTransform(seur_RNA, vst.flavor = \"v2\",\n method=\"glmGamPoi\", verbose = FALSE,\n #conserve.memory = TRUE # also note that return.only.var.genes = TRUE is the default\n ) \ntoc()\ngc()","metadata":{"_kg_hide-output":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:00:17.393151Z","iopub.execute_input":"2023-04-16T05:00:17.394419Z","iopub.status.idle":"2023-04-16T05:04:41.311889Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Dimensionality Reduction: PCA","metadata":{}},{"cell_type":"markdown","source":"Next we perform PCA. By default, only the previously determined variable features are used as input, but can be defined using features argument if you wish to choose a different subset.","metadata":{}},{"cell_type":"code","source":"seur_RNA <- RunPCA(seur_RNA, npcs = 50, verbose = FALSE)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:04:41.315746Z","iopub.execute_input":"2023-04-16T05:04:41.317468Z","iopub.status.idle":"2023-04-16T05:04:56.356807Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Examine and visualize PCA results a few different ways\ncat(\"5 most important positive and negative expression features in the first 2 principle components\\n\")\nprint(seur_RNA[[\"pca\"]], dims = 1:2, nfeatures = 10)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:04:56.366734Z","iopub.execute_input":"2023-04-16T05:04:56.373908Z","iopub.status.idle":"2023-04-16T05:04:56.409805Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Plot the first two PCs**","metadata":{}},{"cell_type":"code","source":"fig(25,8)\n#Standart plots with seurat DimPlot function\nplt1 <- DimPlot(seur_RNA, reduction = \"pca\", group.by = \"cell_type\") + theme_classic(base_size = 24)\nplt2 <- DimPlot(seur_RNA, reduction = \"pca\", group.by = \"day\") + theme_classic(base_size = 24)\nplt3 <- DimPlot(seur_RNA, reduction = \"pca\", group.by = \"donor\") + theme_classic(base_size = 24)\nplt1 + plt2 + plt3","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:04:56.413450Z","iopub.execute_input":"2023-04-16T05:04:56.415144Z","iopub.status.idle":"2023-04-16T05:05:00.471072Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Extract the data for plotting from the Seurat object**","metadata":{}},{"cell_type":"code","source":"pca_color <- seur_RNA@reductions$pca@cell.embeddings %>%\n as.data.frame %>% \n select(PC_1, PC_2) %>%\n mutate(day = as.factor(seur_RNA$day),\n donor = as.factor(seur_RNA$donor),\n cell_type = seur_RNA$cell_type,\n percent_mt = seur_RNA$percent.mt,\n percent_rb = seur_RNA$percent.rb,\n nCount_RNA = seur_RNA$nCount_RNA,\n nFeature_RNA = seur_RNA$nFeature_RNA)\n\npca_color$cell_type <- factor(pca_color$cell_type, levels = c(\n \"HSC\",\"MoP\", \"BP\", \"EryP\", \"MkP\", \"MasP\", \"NeuP\" \n ))\n#add working variables for merging and adding sample sizes in plots \npca_color$Row.names <- rownames(pca_color)\nss_per_day <- pca_color %>% group_by(day) %>% summarise(ss_per_day = n())\nss_per_donor <- pca_color %>% group_by(donor) %>% summarise(ss_per_donor = n())\npca_color <- pca_color %>% left_join(ss_per_day, by = \"day\") %>% mutate(ss_per_day = paste0(day, \", n = \", ss_per_day))\npca_color$ss_per_day <- factor(pca_color$ss_per_day, levels = unique(pca_color$ss_per_day))\npca_color <- pca_color %>% left_join(ss_per_donor, by = \"donor\") %>% mutate(ss_per_donor = paste0(donor, \", n = \", ss_per_donor))\nrm(ss_per_day, ss_per_donor)\n\nhead(pca_color)\n\n#write.csv(pca_color, \"PCA.csv\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:05:00.474981Z","iopub.execute_input":"2023-04-16T05:05:00.477797Z","iopub.status.idle":"2023-04-16T05:05:00.781444Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Animated plots in GIF format - in the output**","metadata":{}},{"cell_type":"code","source":"p1 <- ggplot(pca_color, aes(x = PC_1, y = PC_2, col = cell_type)) +\n geom_point(size = 0.5, alpha = 0.3) +\n theme_classic(base_size = 18) +\n #scale_color_manual(values=c(\"#2D735F\", \"#BF9039\", \"#BF5841\")) +\n theme(aspect.ratio = 1, legend.position=\"none\") +\n # Animating the plot\n labs(title = 'Cell type: {current_frame}') +\n transition_manual(cell_type, cumulative = TRUE)\nanimate(p1, fps = 10)\n\nfile_name <- paste0(path_gif,'Cell_type_PCA.gif')# \nanim_save(file.path(normalizePath(dirname(file_name)),basename(file_name)))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:05:00.785666Z","iopub.execute_input":"2023-04-16T05:05:00.787700Z","iopub.status.idle":"2023-04-16T05:05:04.487688Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"p1 <- ggplot(pca_color, aes(x = PC_1, y = PC_2, col = day)) +\n geom_point(size = 0.5, alpha = 0.3) +\n theme_classic(base_size = 18) +\n scale_color_manual(values=c(\"#2D735F\", \"#BF9039\", \"#BF5841\")) +\n theme(aspect.ratio = 1, legend.position=\"none\") +\n # Animating the plot\n labs(title = 'Day: {current_frame}') +\n transition_manual(ss_per_day, cumulative = TRUE)\nanimate(p1, fps = 30)\n\nfile_name <- paste0(path_gif,'DAY_PCA.gif') #path_gif, \nanim_save(file.path(normalizePath(dirname(file_name)),basename(file_name)))","metadata":{"_kg_hide-output":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:05:04.490351Z","iopub.execute_input":"2023-04-16T05:05:04.491817Z","iopub.status.idle":"2023-04-16T05:05:07.024146Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"p1 <- ggplot(pca_color, aes(x = PC_1, y = PC_2, col = donor)) +\n geom_point(size = 0.5, alpha = 0.3) +\n theme_classic(base_size = 18) +\n scale_color_manual(values=c(\"#2D735F\", \"#BF9039\", \"#BF5841\")) +\n theme(aspect.ratio = 1, legend.position=\"none\") +\n # Animating the plot\n labs(title = 'Donor: {current_frame}') +\n transition_manual(ss_per_donor, cumulative = TRUE)\nanimate(p1, fps = 30)\n\nfile_name <- paste0(path_gif, 'DONOR_PCA.gif') #path_gif,\nanim_save(file.path(normalizePath(dirname(file_name)),basename(file_name)))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:05:07.028188Z","iopub.execute_input":"2023-04-16T05:05:07.029828Z","iopub.status.idle":"2023-04-16T05:05:08.766559Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Ribosomal and mitochondrial genes**","metadata":{}},{"cell_type":"code","source":"fig(22,12)\n\nplt1 <- ggplot(pca_color, aes(x = PC_1, y = PC_2, col = percent_mt)) +\n geom_point(size=0.8) +\n theme_classic(base_size = 24) +\n scale_color_gradient(low='darkblue', high='yellow') +\n theme(aspect.ratio = 1)\n\nplt2 <- ggplot(pca_color, aes(x = PC_1, y = PC_2, col = percent_rb)) +\n geom_point(size=0.8) +\n theme_classic(base_size = 24) +\n scale_color_gradient(low='darkblue', high='yellow') +\n theme(aspect.ratio = 1)\n\nplt3 <- ggplot(pca_color, aes(x = PC_1, y = PC_2, col = nCount_RNA)) +\n geom_point(size=0.8) +\n theme_classic(base_size = 24) +\n scale_color_gradient(low='darkblue', high='yellow') +\n theme(aspect.ratio = 1)\n\nplt4 <- ggplot(pca_color, aes(x = PC_1, y = PC_2, col = nFeature_RNA)) +\n geom_point(size=0.8) +\n theme_classic(base_size = 24) +\n scale_color_gradient(low='darkblue', high='yellow') +\n theme(aspect.ratio = 1)\n\n\nplt1 + plt2 + plt3 + plt4","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:05:08.769887Z","iopub.execute_input":"2023-04-16T05:05:08.771644Z","iopub.status.idle":"2023-04-16T05:05:14.899775Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Load normalized protein and RNA data and sample the same cells**","metadata":{}},{"cell_type":"code","source":"# Load protein normalized data\n#path <- \"/kaggle/input/sparse-raw-counts-data-open-problems-multimodal/citeseq/sp_train_cite_targets_raw.rds\"\npath <- \"/kaggle/input/sparse-measurement-data-open-problems-multimodal/sp_train_cite_targets.rds\"\nmat_prot <- readRDS(path)\n\n#dgCMatrix to matrix\nmat_prot <- as.matrix(mat_prot)\nmat_prot <- mat_prot[selected_cells, ]\n\n# Load RNA normalized data\npath <- \"/kaggle/input/sparse-measurement-data-open-problems-multimodal/sp_train_cite_inputs.rds\"\nmat_RNA <- readRDS(path)\nmat_RNA <- mat_RNA[selected_cells, ]\n\ngc()","metadata":{"_kg_hide-output":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:05:14.903455Z","iopub.execute_input":"2023-04-16T05:05:14.906097Z","iopub.status.idle":"2023-04-16T05:06:30.617894Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Functions for subsetting expression data and plot PC1 vs PC2 colored by the proteins and RNAs of interest**","metadata":{}},{"cell_type":"markdown","source":"We will plot RNA and proteins (mostly, from [this notebook](http://www.kaggle.com/code/alexandervc/cite-seq-cd53-vs-cd45-correlation), and two enrichment analyses of top correlated genes unique for [CD45](http://www.kaggle.com/code/antoninadolgorukova/mmscel-cd45-iso-go-enrichment#4.-Translation,-transcription,-splicing) iso or [CD36](http://www.kaggle.com/code/antoninadolgorukova/mmscel-cd36-go-enrichment#--Terms-related-to-translation,-transcription,-splicing) and their RNA) that can be somehow related to CD45 and CD36 regulation","metadata":{}},{"cell_type":"code","source":"#The function subsets RNAs and proteins of interest, specified in the prots and rnas arguments,\n#from the main data matrices. By default, it uses exact matching, but with exact = FALSE\n#can find additional RNA and proteins that have a specified pattern in the name.\n#The function also adds a column with the mean expression of the RNA and proteins of interest.\nmake_dat_for_PCplot <- function(prots = c(), rnas = c(), pca_color, exact = TRUE) {\n \n #Proteins\n if(!is.null(prots)) {\n \n if(exact == TRUE) {\n \n prots <- paste0(prots, \"$\") #gsub(\"[[:>:]]\", \"$\", prots, perl=TRUE) \n } \n \n f_prots <- c(\n colnames(mat_prot)[grepl(paste0(prots, collapse = \"|\"),\n colnames(mat_prot))])\n my_prot <- mat_prot[, f_prots] %>% as.data.frame\n my_prot$Row.names <- rownames(my_prot)\n \n #add mean of the genes\n my_prot$Mean_prots <- apply(my_prot[, -ncol(my_prot)], 1, mean, na.rm = TRUE)\n \n pca_color_long <- merge(pca_color, my_prot, by = \"Row.names\")\n \n prots <- gsub(\"\\\\$\", \"\", prots)\n \n #print results\n cat(\"\\nProteins found:\", length(f_prots), \"\\n\")\n print(f_prots)\n cat(\"\\nProteins not found:\", length(setdiff(prots, f_prots)), \"\\n\")\n print(setdiff(prots, f_prots))\n \n } else {\n pca_color_long <- pca_color\n f_prots <- NULL\n } \n \n # RNAs\n if(!is.null(rnas)) {\n \n if(exact == TRUE) {\n \n rnas <- paste0(\"_\", rnas, \"$\") \n }\n \n f_rnas <- c(\n colnames(mat_RNA)[grepl(paste0(rnas, collapse = \"|\"),\n colnames(mat_RNA))])\n my_RNA <- mat_RNA[, f_rnas] %>% as.data.frame\n my_RNA$Row.names <- rownames(my_RNA)\n \n #add mean of the genes\n my_RNA$Mean_RNAs <- apply(my_RNA[, -ncol(my_RNA)], 1, mean, na.rm = TRUE)\n \n pca_color_long <- merge(pca_color_long, my_RNA, by = \"Row.names\")\n \n rnas <- gsub(\"\\\\$\", \"\", rnas)\n \n cat(\"\\nRNA found:\", length(f_rnas), \"\\n\")\n print(f_rnas)\n cat(\"\\nRNA not found:\", length(setdiff(gsub(\".*_\", \"\", rnas),\n gsub(\".*_\", \"\", f_rnas))), \"\\n\")\n print(setdiff(gsub(\".*_\", \"\", rnas),\n gsub(\".*_\", \"\", f_rnas)))\n \n \n } else {\n f_rnas <- NULL\n }\n \n pca_color_long <- pivot_longer(pca_color_long, cols = (ncol(pca_color)+1):ncol(pca_color_long),\n names_to = \"Name\", values_to = \"Expression\")\n \n # calculate median expression of each RNA for binarisation\n pca_color_long$`Expression\\nto median` <- NA\n for(nm in unique(pca_color_long$Name) ) {\n\n tmp <- pca_color_long[pca_color_long$Name == nm, ]\n med <- median(tmp$Expression)\n pca_color_long[pca_color_long$Name == nm, ]$`Expression\\nto median` <-\n ifelse(pca_color_long[pca_color_long$Name == nm, ]$Expression > med, \"higher\", \"lower\")\n }\n\n # exclude \"outliers\" - extrim values\n\n pca_color_long$`Expression\\n(corrected)` <- NA\n for(nm in unique(pca_color_long$Name) ) {\n\n tmp <- pca_color_long[pca_color_long$Name == nm, ]\n q5 <- quantile(tmp$Expression, 0.05)\n q95 <- quantile(tmp$Expression, 0.95)\n pca_color_long[pca_color_long$Name == nm, ]$`Expression\\n(corrected)` <-\n ifelse(pca_color_long[pca_color_long$Name == nm, ]$Expression > q95, q95,\n ifelse(pca_color_long[pca_color_long$Name == nm, ]$Expression < q5, q5,\n pca_color_long[pca_color_long$Name == nm, ]$Expression) )\n }\n\n return(pca_color_long)\n}\n\nprint_plots <- function(sbst, da = \"3\", do = \"32606\") {\n \n for (p in unique(sbst$Name) ) {\n \n plt1 <- ggplot(sbst[sbst$Name == p , ], #& sbst$cell_type != \"MkP\"\n aes(x = PC_1, y = PC_2, col = `Expression\\n(corrected)`) ) +\n geom_point(size = 0.7, alpha = 0.5) +\n theme_classic(base_size = 16) +\n scale_color_gradient(low='darkblue', high='yellow') +\n theme(aspect.ratio = 1) +\n facet_wrap( ~ Name) +\n ggtitle(paste0(\"All days, all donors\",\n \"\\nn = \",nrow(sbst[sbst$Name == p , ]), \" cells\"))\n\n plt2 <- ggplot(sbst[sbst$Name == p, ], # & sbst$cell_type != \"MkP\"\n aes(x = PC_1, y = PC_2, col = `Expression\\nto median`) ) +\n geom_point(size = 0.7, alpha = 0.5) +\n theme_classic(base_size = 16) +\n scale_colour_manual(values = c('yellow', 'darkblue')) +\n theme(aspect.ratio = 1) +\n facet_wrap( ~ Name) +\n ggtitle(paste0(\"All days, all donors\",\n \"\\nn = \",nrow(sbst[sbst$Name == p , ]), \" cells\"))\n \n dd_fixed <- sbst %>% filter(Name == p, day %in% da, donor %in% do) #, cell_type != \"MkP\"\n \n plt3 <- ggplot(dd_fixed,\n aes(x = PC_1, y = PC_2, col = `Expression\\n(corrected)`) ) +\n geom_point(size = 0.7, alpha = 0.5) +\n theme_classic(base_size = 16) +\n scale_color_gradient(low='darkblue', high='yellow') +\n theme(aspect.ratio = 1) +\n ggtitle(paste0(\"Day \", paste(unique(dd_fixed$day), collapse = \", \"),\n \", donor \", paste(unique(dd_fixed$donor), collapse = \", \"),\n \"\\nn = \",nrow(dd_fixed), \" cells\")) +\n facet_wrap( ~ Name)\n \n plt4 <- ggplot(dd_fixed,\n aes(x = PC_1, y = PC_2, col = `Expression\\nto median`) ) +\n geom_point(size = 0.7, alpha = 0.5) +\n theme_classic(base_size = 16) +\n scale_colour_manual(values = c('yellow', 'darkblue')) +\n theme(aspect.ratio = 1) +\n ggtitle(paste0(\"Day \", paste(unique(dd_fixed$day), collapse = \", \"),\n \", donor \", paste(unique(dd_fixed$donor), collapse = \", \"),\n \"\\nn = \",nrow(dd_fixed), \" cells\")) + \n facet_wrap( ~ Name)\n\n print(plt1 + plt2 + plt3 + plt4 + plot_layout(ncol = 4))\n print(paste0(\"-----------------------------\", p, \"----------------------------------\"))\n } \n}","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:06:30.620533Z","iopub.execute_input":"2023-04-16T05:06:30.622006Z","iopub.status.idle":"2023-04-16T05:06:30.639927Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# CD36 and CD45 isoforms vs PC_2/PC_1","metadata":{}},{"cell_type":"code","source":"rnas <- c('PTPRC','CD36')\nprots <- c('CD45', 'CD45RO', 'CD45RA', 'CD36')\nsbst <- make_dat_for_PCplot(prots, rnas, pca_color)\n\nhead(sbst, n = n_distinct(sbst$Name))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25,20)\nprot_of_interest <- \"CD36\"\nplt <- ggscatter(sbst[sbst$Name == prot_of_interest, ], x = \"Expression\", y = \"PC_2\", color = \"PC_1\",\n shape = 20, size = 3, #, alpha = 0.3\n add = \"reg.line\", conf.int = TRUE,\n cor.coef = TRUE, cor.coef.size = 6,\n cor.coeff.args = list(method = \"spearman\", label.sep = \"\\n\"),\n add.params = list(fill = \"lightgray\"),\n xlab = prot_of_interest,\n title = paste0(\"Protein expresion vs PC2, all cells (n = \", nrow(sbst[sbst$Name == prot_of_interest,]), \")\"),\n ggtheme = theme_bw(base_size = 20)) + theme(aspect.ratio = 1)\n\nplt1 <- ggMarginal(plt, type=\"histogram\")\n\nprot_of_interest <- \"CD45\"\nplt <- ggscatter(sbst[sbst$Name == prot_of_interest, ], x = \"Expression\", y = \"PC_2\", color = \"PC_1\",\n shape = 20, size = 3, #, alpha = 0.3\n add = \"reg.line\", conf.int = TRUE,\n cor.coef = TRUE, cor.coef.size = 6,\n cor.coeff.args = list(method = \"spearman\", label.sep = \"\\n\"),\n add.params = list(fill = \"lightgray\"),\n xlab = paste(prot_of_interest, \"RNA\"),\n title = paste0(\"Protein expresion vs PC2, all cells (n = \", nrow(sbst[sbst$Name == prot_of_interest,]), \")\"),\n ggtheme = theme_bw(base_size = 20)) + theme(aspect.ratio = 1)\n\nplt2 <- ggMarginal(plt, type=\"histogram\")\n\nprot_of_interest <- \"CD45RO\"\nplt <- ggscatter(sbst[sbst$Name == prot_of_interest, ], x = \"Expression\", y = \"PC_2\", color = \"PC_1\",\n shape = 20, size = 3, #, alpha = 0.3\n add = \"reg.line\", conf.int = TRUE,\n cor.coef = TRUE, cor.coef.size = 6,\n cor.coeff.args = list(method = \"spearman\", label.sep = \"\\n\"),\n add.params = list(fill = \"lightgray\"),\n xlab = paste(prot_of_interest, \"RNA\"),\n title = paste0(\"Protein expresion vs PC2, all cells (n = \", nrow(sbst[sbst$Name == prot_of_interest,]), \")\"),\n ggtheme = theme_bw(base_size = 20)) + theme(aspect.ratio = 1)\n\nplt3 <- ggMarginal(plt, type=\"histogram\")\n\nprot_of_interest <- \"CD45RA\"\nplt <- ggscatter(sbst[sbst$Name == prot_of_interest, ], x = \"Expression\", y = \"PC_2\", color = \"PC_1\",\n shape = 20, size = 3, #, alpha = 0.3\n add = \"reg.line\", conf.int = TRUE,\n cor.coef = TRUE, cor.coef.size = 6,\n cor.coeff.args = list(method = \"spearman\", label.sep = \"\\n\"),\n add.params = list(fill = \"lightgray\"),\n xlab = paste(prot_of_interest, \"RNA\"),\n title = paste0(\"Protein expresion vs PC2, all cells (n = \", nrow(sbst[sbst$Name == prot_of_interest,]), \")\"),\n ggtheme = theme_bw(base_size = 20)) + theme(aspect.ratio = 1)\n\nplt4 <- ggMarginal(plt, type=\"histogram\")\n\ngrid.arrange(plt1, plt2, plt3, plt4, ncol = 2)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**3D plots**","metadata":{}},{"cell_type":"code","source":"fig(15,20)\nprot_of_interest <- \"CD36\"\nplt <- plot_ly(sample_n(sbst[sbst$Name == prot_of_interest,], 5000), x = ~PC_1, y = ~PC_2, z = ~Expression, color = ~cell_type,\n size = 0.5) %>% \n add_markers() %>%\n layout(title = paste(prot_of_interest, \"expression vs first two principal components (sample = 5000)\"))\n#plt\nfile_name = paste0(path_3d, gsub(\".*_\", \"\", prot_of_interest), \"_PCA3d.HTML\")\nhtmlwidgets::saveWidget(partial_bundle(plt), file.path(normalizePath(dirname(file_name)),basename(file_name)))\ndisplay_html(paste0(''))","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(15,20) \nprot_of_interest <- \"ENSG00000135218_CD36\"\nplt <- plot_ly(sample_n(sbst[sbst$Name == prot_of_interest,], 5000), x = ~PC_1, y = ~PC_2, z = ~Expression, color = ~cell_type,\n size = 0.5) %>% \n add_markers() %>%\n layout(title = paste(prot_of_interest, \"expression vs first two principal components (sample = 5000)\"))\n#plt\nfile_name = paste0(path_3d, gsub(\".*_\", \"\", prot_of_interest), \"RNA_PCA3d.HTML\")\nhtmlwidgets::saveWidget(partial_bundle(plt), file.path(normalizePath(dirname(file_name)),basename(file_name)))\ndisplay_html(paste0(''))","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(15,20) \nprot_of_interest <- \"CD45RA\"\nplt <- plot_ly(sample_n(sbst[sbst$Name == prot_of_interest,], 5000), x = ~PC_1, y = ~PC_2, z = ~Expression, color = ~cell_type,\n size = 0.5) %>% \n add_markers() %>%\n layout(title = paste(prot_of_interest, \"expression vs first two principal components (sample = 5000)\"))\n#plt\nfile_name = paste0(path_3d, gsub(\".*_\", \"\", prot_of_interest), \"_PCA3d.HTML\")\nhtmlwidgets::saveWidget(partial_bundle(plt), file.path(normalizePath(dirname(file_name)),basename(file_name)))\ndisplay_html(paste0(''))","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(15,20) \nprot_of_interest <- \"CD45RO\"\nplt <- plot_ly(sample_n(sbst[sbst$Name == prot_of_interest,], 5000), x = ~PC_1, y = ~PC_2, z = ~Expression, color = ~cell_type,\n size = 0.5) %>% \n add_markers() %>%\n layout(title = paste(prot_of_interest, \"expression vs first two principal components (sample = 5000)\"))\n#plt\nfile_name = paste0(path_3d, gsub(\".*_\", \"\", prot_of_interest), \"_PCA3d.HTML\")\nhtmlwidgets::saveWidget(partial_bundle(plt), file.path(normalizePath(dirname(file_name)),basename(file_name)))\ndisplay_html(paste0(''))","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(15,20) \nprot_of_interest <- \"CD45\"\nplt <- plot_ly(sample_n(sbst[sbst$Name == prot_of_interest,], 5000), x = ~PC_1, y = ~PC_2, z = ~Expression, color = ~cell_type,\n size = 0.5) %>% \n add_markers() %>%\n layout(title = paste(prot_of_interest, \"expression vs first two principal components (sample = 5000)\"))\n#plt\nfile_name = paste0(path_3d, gsub(\".*_\", \"\", prot_of_interest), \"_PCA3d.HTML\")\nhtmlwidgets::saveWidget(partial_bundle(plt), file.path(normalizePath(dirname(file_name)),basename(file_name)))\ndisplay_html(paste0(''))","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(15,20) \nprot_of_interest <- \"ENSG00000081237_PTPRC\"\nplt <- plot_ly(sample_n(sbst[sbst$Name == prot_of_interest,], 5000), x = ~PC_1, y = ~PC_2, z = ~Expression, color = ~cell_type,\n size = 0.5) %>% \n add_markers() %>%\n layout(title = paste(prot_of_interest, \"expression vs first two principal components (sample = 5000)\"))\n#plt\nfile_name = paste0(path_3d, gsub(\".*_\", \"\", prot_of_interest), \"_PCA3d.HTML\")\nhtmlwidgets::saveWidget(partial_bundle(plt), file.path(normalizePath(dirname(file_name)),basename(file_name)))\ndisplay_html(paste0(''))","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**The 3D plots are stored in the output in HTML format**","metadata":{}},{"cell_type":"code","source":"# # EryP cells\n# eryp <- filter(sbst, grepl(\"CD36|CD45\", Name) & cell_type == \"EryP\")\n# eryp <- eryp %>% select(\"Row.names\", \"PC_1\", \"PC_2\", \"day\", \"donor\",\"cell_type\", \"Name\",\"Expression\") %>%\n# pivot_wider(names_from = \"Name\", values_from = \"Expression\")\n# cat(\"The head the subset with EryP cells, the protein of interest and its RNA\")\n# head(eryp)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# CD45 isoforms, CD36 and their RNA","metadata":{}},{"cell_type":"code","source":"fig(25,6)\nsbst$Name <- \n factor(sbst$Name,\n levels = c('CD45RA','CD45RO','CD45',\n 'ENSG00000081237_PTPRC',\n 'CD36','ENSG00000135218_CD36',\n 'Mean_RNAs', 'Mean_prots'))\nsbst <- sbst[order(sbst$Name), ]\nprint_plots(sbst[!sbst$Name %in% c('Mean_RNAs', 'Mean_prots'),])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Examine change by day with animated plot - in the output (GIF)**","metadata":{}},{"cell_type":"code","source":"prot_of_interest <- \"ENSG00000081237_PTPRC\"\n\np1 <- ggplot(sbst[sbst$Name == prot_of_interest, ],\n aes(x = PC_1, y = PC_2, color = `Expression\\nto median`)) +\n geom_point(size = 0.5) +\n theme_classic(base_size = 18) +\n scale_color_manual(values=c('yellow', 'darkblue')) +\n #scale_color_gradient(low='darkblue', high='yellow') +\n theme(aspect.ratio = 1) + #, legend.position=\"none\"\n # Animating the plot\n labs(title = paste0(prot_of_interest),\n subtitle = 'Day: {current_frame}') +\n transition_manual(ss_per_day, cumulative = FALSE)\nanimate(p1, fps = 30)\n\nfile_name <- paste0(path_gif, gsub(\".*_\",\"\", prot_of_interest), '_by_Day_PCA.gif')\nanim_save(file.path(normalizePath(dirname(file_name)),basename(file_name)))","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"prot_of_interest <- \"CD45RA\"\n\np1 <- ggplot(sbst[sbst$Name == prot_of_interest, ],\n aes(x = PC_1, y = PC_2, color = `Expression\\nto median`)) +\n geom_point(size = 0.5) +\n theme_classic(base_size = 18) +\n scale_color_manual(values=c('yellow', 'darkblue')) +\n #scale_color_gradient(low='darkblue', high='yellow') +\n theme(aspect.ratio = 1) + #, legend.position=\"none\"\n # Animating the plot\n labs(title = paste0(prot_of_interest),\n subtitle = 'Day: {current_frame}') +\n transition_manual(ss_per_day, cumulative = FALSE)\nanimate(p1, fps = 30)\n\nfile_name <- paste0(path_gif, gsub(\".*_\",\"\", prot_of_interest), '_by_Day_PCA.gif')\nanim_save(file.path(normalizePath(dirname(file_name)),basename(file_name)))","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"prot_of_interest <- \"CD45RO\"\n\np1 <- ggplot(sbst[sbst$Name == prot_of_interest, ],\n aes(x = PC_1, y = PC_2, color = `Expression\\nto median`)) +\n geom_point(size = 0.5) +\n theme_classic(base_size = 18) +\n scale_color_manual(values=c('yellow', 'darkblue')) +\n #scale_color_gradient(low='darkblue', high='yellow') +\n theme(aspect.ratio = 1) + #, legend.position=\"none\"\n # Animating the plot\n labs(title = paste0(prot_of_interest),\n subtitle = 'Day: {current_frame}') +\n transition_manual(ss_per_day, cumulative = FALSE)\nanimate(p1, fps = 30)\n\nfile_name <- paste0(path_gif, gsub(\".*_\",\"\", prot_of_interest), '_by_Day_PCA.gif')\nanim_save(file.path(normalizePath(dirname(file_name)),basename(file_name)))","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"prot_of_interest <- \"CD36\"\n\np1 <- ggplot(sbst[sbst$Name == prot_of_interest, ],\n aes(x = PC_1, y = PC_2, color = `Expression\\nto median`)) +\n geom_point(size = 0.5) +\n theme_classic(base_size = 18) +\n scale_color_manual(values=c('yellow', 'darkblue')) +\n #scale_color_gradient(low='darkblue', high='yellow') +\n theme(aspect.ratio = 1) + #, legend.position=\"none\"\n # Animating the plot\n labs(title = paste0(prot_of_interest),\n subtitle = 'Day: {current_frame}') +\n transition_manual(ss_per_day, cumulative = FALSE)\nanimate(p1, fps = 30)\n\nfile_name <- paste0(path_gif, gsub(\".*_\",\"\", prot_of_interest), '_by_Day_PCA.gif')\nanim_save(file.path(normalizePath(dirname(file_name)),basename(file_name)))","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"prot_of_interest <- \"ENSG00000135218_CD36\"\n\np1 <- ggplot(sbst[sbst$Name == prot_of_interest, ],\n aes(x = PC_1, y = PC_2, color = `Expression\\nto median`)) +\n geom_point(size = 0.5) +\n theme_classic(base_size = 18) +\n scale_color_manual(values=c('yellow', 'darkblue')) +\n #scale_color_gradient(low='darkblue', high='yellow') +\n theme(aspect.ratio = 1) + #, legend.position=\"none\"\n # Animating the plot\n labs(title = paste0(prot_of_interest),\n subtitle = 'Day: {current_frame}') +\n transition_manual(ss_per_day, cumulative = FALSE)\nanimate(p1, fps = 30)\n\nfile_name <- paste0(path_gif, gsub(\".*_\",\"\", prot_of_interest), 'RNA_by_Day_PCA.gif')\nanim_save(file.path(normalizePath(dirname(file_name)),basename(file_name)))","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Ribosomal genes","metadata":{}},{"cell_type":"markdown","source":"**Calculate correlation matrix**","metadata":{}},{"cell_type":"code","source":"f_rnas <- c(colnames(mat_RNA)[grepl(\"_RP\",\n colnames(mat_RNA),\n ignore.case = TRUE)])\ncat(\"Found\", f_rnas %>% length, \"genes starting with _RP\")\nrib_genes <- mat_RNA[, f_rnas] %>% as.data.frame\n\n#exclude genes with 0 expression in all cells\ncat(\"\\nZero expression in all cells:\")\nnames(colSums(rib_genes)[colSums(rib_genes) == 0])\n\nrib_genes <- select(rib_genes, -c(names(colSums(rib_genes)[colSums(rib_genes) == 0])) )\n\n# correlation matrix\nrib_cor_mat <- cor(rib_genes, method = \"spearman\")\n\ncat(\"\\nCorrelation matix for all found ribosomal genes with non-zero expression\")\nhead(rib_cor_mat, 5)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:45:50.249730Z","iopub.execute_input":"2023-04-16T05:45:50.251408Z","iopub.status.idle":"2023-04-16T05:45:58.654787Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Plot a clustered heatmap**","metadata":{}},{"cell_type":"code","source":"fig(20, 20)\nall_rib_heat <- pheatmap(rib_cor_mat, fontsize = 14, fontsize_row = 3, fontsize_col = 3,\n color = myColors,\n breaks = breaksList,\n angle_col = 90,\n main = paste0(nrow(rib_cor_mat), \" ribosomal genes with non-zero expression\"))\n\n#Re-order original data (genes) to match ordering in heatmap (top-to-bottom)\n#rownames(rib_cor_mat[rib_heat$tree_row[[\"order\"]],])","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:45:58.657340Z","iopub.execute_input":"2023-04-16T05:45:58.658830Z","iopub.status.idle":"2023-04-16T05:46:00.396813Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Locate and remove the biggest cluster with around 0 correlations**","metadata":{}},{"cell_type":"code","source":"cat(\"The number of genes per cluster (at the level shown by the red line)\")\nsort(cutree(all_rib_heat$tree_row, h = 2.4)) %>% table\n\nfig(30, 10)\nplot(all_rib_heat$tree_row)\nabline(h = 2.4, col = \"red\", lty = 3, lwd = 1)\n\n# mate a subset of genes without the biggest empty cluster\nfilt_rib <- cutree(all_rib_heat$tree_row, h = 2.4)\nfilt_rib <- filt_rib[filt_rib != 1] #%>% length\n\nrib_to_plot <- rib_cor_mat[names(filt_rib), names(filt_rib)]","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:46:00.400162Z","iopub.execute_input":"2023-04-16T05:46:00.402076Z","iopub.status.idle":"2023-04-16T05:46:00.961958Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Plot the remaining clusters**","metadata":{}},{"cell_type":"code","source":"fig(25, 25)\n\ncut_rib_heat <- pheatmap(rib_to_plot, cutree_cols = 4, cutree_rows = 4,\n color = myColors,\n breaks = breaksList,\n angle_col = 90, fontsize = 14,\n #filename = \"rib105_heatmap.png\",\n main = paste0(nrow(rib_to_plot), \" ribosomal genes left after removal of the empty cluster\"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:46:00.964524Z","iopub.execute_input":"2023-04-16T05:46:00.966093Z","iopub.status.idle":"2023-04-16T05:46:02.602288Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Make subsets by cluster and select mean RNAs to color the PCA plot**","metadata":{}},{"cell_type":"code","source":"cat(\"The number of genes per cluster (at the level shown by the red line)\")\nsort(cutree(cut_rib_heat$tree_row, h = 3.6)) %>% table\n\nrib_to_plot <- rib_to_plot %>% as.data.frame %>%\n mutate(cluster = cutree(cut_rib_heat$tree_row, h = 3.6) )\n\nfig(30, 10)\nplot(cut_rib_heat$tree_row)\nabline(h = 3.6, col = \"red\", lty = 3, lwd = 1)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:46:02.606972Z","iopub.execute_input":"2023-04-16T05:46:02.610689Z","iopub.status.idle":"2023-04-16T05:46:02.978553Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sbst1 <- rib_to_plot[rib_to_plot$cluster == 1, ]\ncat(\"Cluster 1:\", nrow(sbst1), \"genes\")\nrnas <- c(gsub(\".*_\", \"\", rownames(sbst1)))\nprots <- NULL\nsbst1 <- make_dat_for_PCplot(prots, rnas, pca_color, exact = TRUE)\nsbst1[sbst1$Name == \"Mean_RNAs\", ]$Name <- \"Mean_Rib_clust1\"\nsbst <- sbst1[sbst1$Name == \"Mean_Rib_clust1\", ]\n\nsbst1 <- rib_to_plot[rib_to_plot$cluster == 2, ]\ncat(\"\\nCluster 2:\", nrow(sbst1), \"genes\")\nrnas <- c(gsub(\".*_\", \"\", rownames(sbst1)))\nprots <- NULL\nsbst1 <- make_dat_for_PCplot(prots, rnas, pca_color, exact = TRUE)\nsbst1[sbst1$Name == \"Mean_RNAs\", ]$Name <- \"Mean_Rib_clust2\"\nsbst <- rbind(sbst, sbst1[sbst1$Name == \"Mean_Rib_clust2\", ])\n\nsbst1 <- rib_to_plot[rib_to_plot$cluster == 3, ]\ncat(\"\\nCluster 3:\", nrow(sbst1), \"genes\")\nrnas <- c(gsub(\".*_\", \"\", rownames(sbst1)))\nprots <- NULL\nsbst1 <- make_dat_for_PCplot(prots, rnas, pca_color, exact = TRUE)\nsbst1[sbst1$Name == \"Mean_RNAs\", ]$Name <- \"Mean_Rib_clust3\"\nsbst <- rbind(sbst, sbst1[sbst1$Name == \"Mean_Rib_clust3\", ])\n\nsbst1 <- rib_to_plot[rib_to_plot$cluster == 4, ]\ncat(\"\\nCluster 4:\", nrow(sbst1), \"genes\")\nrnas <- c(gsub(\".*_\", \"\", rownames(sbst1)))\nprots <- NULL\nsbst1 <- make_dat_for_PCplot(prots, rnas, pca_color, exact = TRUE)\nsbst1[sbst1$Name == \"Mean_RNAs\", ]$Name <- \"Mean_Rib_clust4\"\nsbst <- rbind(sbst, sbst1[sbst1$Name == \"Mean_Rib_clust4\", ])\n \nhead(sbst, n = 5)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:46:02.982070Z","iopub.execute_input":"2023-04-16T05:46:02.984389Z","iopub.status.idle":"2023-04-16T05:46:50.689167Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(22,6)\nprint_plots(sbst)","metadata":{"execution":{"iopub.status.busy":"2023-04-16T05:46:50.691502Z","iopub.execute_input":"2023-04-16T05:46:50.692874Z","iopub.status.idle":"2023-04-16T05:47:06.698348Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Mitochondrial genes","metadata":{}},{"cell_type":"code","source":"f_rnas <- c(colnames(mat_RNA)[grepl(\"MT-\",\n colnames(mat_RNA),\n ignore.case = TRUE)])\ncat(\"Found\", f_rnas %>% length, \"genes starting with MT-\")\ngenes <- mat_RNA[, f_rnas] %>% as.data.frame\n\n#exclude genes with 0 expression in all cells - no such genes here\n#cat(\"\\nZero expression in all cells:\")\n#names(colSums(genes)[colSums(genes) == 0])\n#genes <- select(genes, -c(names(colSums(genes)[colSums(genes) == 0])) )\n\n# correlation matrix\ncor_mat <- cor(genes, method = \"spearman\")\n\ncat(\"\\nCorrelation matix for all found mitochondrial genes\")\nhead(cor_mat, 5)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:55:24.072397Z","iopub.execute_input":"2023-04-16T05:55:24.074228Z","iopub.status.idle":"2023-04-16T05:55:24.684003Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25, 10)\nsum_expr_RNA <- apply(genes, 2, median) %>% sort(decreasing = TRUE)\nsum_expr_RNA <- data.frame(\"Median_expression\" = sum_expr_RNA,\n \"Name\" = names(sum_expr_RNA))\nsum_expr_RNA$Name <- factor(sum_expr_RNA$Name, levels = unique(sum_expr_RNA$Name ))\n#hline <- sum_expr_RNA[sum_expr_RNA$Name == \"PTPRC\",]$Median_expression\n\nggplot(sum_expr_RNA, aes(x = Name, y = Median_expression)) +\n geom_bar(stat = 'identity', fill=\"#BF9039\", col=\"grey\") + \n# geom_hline(yintercept=sum_expr_RNA23[sum_expr_RNA23$Name == \"PTPRC\",]$Median_expression,\n# linetype=\"dashed\")+\n theme_bw(base_size = 22) +\n xlab(\"\") +\n ggtitle(paste0(\"CITEseq 2022\")) +\n theme(axis.text.x = element_text(angle = 45, hjust=1))\n\ncat(\"Zero median expression:\\n\")\ncat(paste0('\"', sum_expr_RNA[sum_expr_RNA$Median_expression == 0,]$Name, '\"', collapse = \", \"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:47:07.369160Z","iopub.execute_input":"2023-04-16T05:47:07.370959Z","iopub.status.idle":"2023-04-16T05:47:07.985797Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"not_plot <- c(\"ENSG00000210127_MT-TA\", \"ENSG00000210140_MT-TC\", \"ENSG00000210154_MT-TD\",\n \"ENSG00000210194_MT-TE\", \"ENSG00000210049_MT-TF\", \"ENSG00000210164_MT-TG\",\n \"ENSG00000210176_MT-TH\", \"ENSG00000210100_MT-TI\", \"ENSG00000210156_MT-TK\",\n \"ENSG00000209082_MT-TL1\", \"ENSG00000210191_MT-TL2\", \"ENSG00000210112_MT-TM\",\n \"ENSG00000210135_MT-TN\", \"ENSG00000210196_MT-TP\", \"ENSG00000210107_MT-TQ\",\n \"ENSG00000210174_MT-TR\", \"ENSG00000210151_MT-TS1\", \"ENSG00000210184_MT-TS2\",\n \"ENSG00000210195_MT-TT\", \"ENSG00000210077_MT-TV\", \"ENSG00000210117_MT-TW\", \"ENSG00000210144_MT-TY\")\n\ncor_mat <- cor_mat[!rownames(cor_mat) %in% not_plot,\n !colnames(cor_mat) %in% not_plot]\n\nfig(10, 10)\nall_mit_heat <- pheatmap(cor_mat, fontsize = 14,\n color = myColors,\n breaks = breaksList,\n angle_col = 90,\n main = paste0(nrow(cor_mat), \" mitochondrial genes with non-zero expression\"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:47:07.988256Z","iopub.execute_input":"2023-04-16T05:47:07.989664Z","iopub.status.idle":"2023-04-16T05:47:08.217623Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rnas <- colnames(cor_mat)[!colnames(cor_mat) %in% not_plot]\nrnas <- gsub(\".*_\", \"\", rnas)\nprots <- NULL\n\nsbst <- make_dat_for_PCplot(prots, rnas, pca_color, exact = TRUE)\n\nhead(sbst, n = 5) #n_distinct(sbst$Name)","metadata":{"_kg_hide-input":false,"execution":{"iopub.status.busy":"2023-04-16T05:47:08.219985Z","iopub.execute_input":"2023-04-16T05:47:08.221363Z","iopub.status.idle":"2023-04-16T05:47:12.650077Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(22,6)\nprint_plots(sbst)","metadata":{"execution":{"iopub.status.busy":"2023-04-16T05:47:12.652488Z","iopub.execute_input":"2023-04-16T05:47:12.654298Z","iopub.status.idle":"2023-04-16T05:48:18.653173Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Mitochondrial translation","metadata":{}},{"cell_type":"code","source":"f_rnas <- c(colnames(mat_RNA)[grepl(\"MRPL\",\n colnames(mat_RNA),\n ignore.case = TRUE)])\ncat(\"Found\", f_rnas %>% length, \"genes starting with MRPL\")\ngenes <- mat_RNA[, f_rnas] %>% as.data.frame\n\n#exclude genes with 0 expression in all cells - no such genes here\n#cat(\"\\nZero expression in all cells:\")\n#names(colSums(genes)[colSums(genes) == 0])\n#genes <- select(genes, -c(names(colSums(genes)[colSums(genes) == 0])) )\n\n# correlation matrix\ncor_mat <- cor(genes, method = \"spearman\")\n\ncat(\"\\nCorrelation matix for all found genes\")\nhead(cor_mat, 5)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:48:18.655757Z","iopub.execute_input":"2023-04-16T05:48:18.657177Z","iopub.status.idle":"2023-04-16T05:48:19.473134Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25, 10)\nsum_expr_RNA <- apply(genes, 2, median) %>% sort(decreasing = TRUE)\nsum_expr_RNA <- data.frame(\"Median_expression\" = sum_expr_RNA,\n \"Name\" = names(sum_expr_RNA))\nsum_expr_RNA$Name <- factor(sum_expr_RNA$Name, levels = unique(sum_expr_RNA$Name ))\n#hline <- sum_expr_RNA[sum_expr_RNA$Name == \"PTPRC\",]$Median_expression\n\nggplot(sum_expr_RNA, aes(x = Name, y = Median_expression)) +\n geom_bar(stat = 'identity', fill=\"#BF9039\", col=\"grey\") + \n# geom_hline(yintercept=sum_expr_RNA23[sum_expr_RNA23$Name == \"PTPRC\",]$Median_expression,\n# linetype=\"dashed\")+\n theme_bw(base_size = 22) +\n xlab(\"\") +\n ggtitle(paste0(\"CITEseq 2022\")) +\n theme(axis.text.x = element_text(angle = 45, hjust=1))\n\ncat(\"Zero median expression:\\n\")\ncat(paste0('\"', sum_expr_RNA[sum_expr_RNA$Median_expression == 0,]$Name, '\"', collapse = \", \"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:48:19.475645Z","iopub.execute_input":"2023-04-16T05:48:19.477059Z","iopub.status.idle":"2023-04-16T05:48:20.091501Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"not_plot <- c(\"ENSG00000185414_MRPL30\", \"ENSG00000226253_MRPL35P3\", \"ENSG00000266946_MRPL37P1\",\n \"ENSG00000204316_MRPL38\", \"ENSG00000256037_MRPL40P1\", \"ENSG00000242810_MRPL42P6\",\n \"ENSG00000228782_MRPL45P2\", \"ENSG00000259494_MRPL46\", \"ENSG00000149792_MRPL49\", \"ENSG00000204822_MRPL53\")\n\ncor_mat <- cor_mat[!rownames(cor_mat) %in% not_plot,\n !colnames(cor_mat) %in% not_plot]\nfig(15, 15)\nall_mit_heat <- pheatmap(cor_mat, fontsize = 14,\n color = myColors,\n breaks = breaksList,\n angle_col = 90,\n main = paste0(nrow(cor_mat), \" mitochondrial genes with non-zero expression\"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:48:20.094265Z","iopub.execute_input":"2023-04-16T05:48:20.095803Z","iopub.status.idle":"2023-04-16T05:48:20.531732Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Here we plot RNA, top negatively correlated with CD45 RNA and/or CD36 protein and forming the enriched terms related to mitochondrial translation (found [here](http://www.kaggle.com/code/antoninadolgorukova/mmscel-cd45-iso-go-enrichment#4.-Translation,-transcription,-splicing) and [here](http://www.kaggle.com/code/antoninadolgorukova/mmscel-cd36-go-enrichment#--Terms-related-to-translation,-transcription,-splicing))","metadata":{}},{"cell_type":"code","source":"rnas <- c('MRPL41', 'MRPL52', 'MRPL27',\n 'AURKAIP1', 'TUFM', 'SEC61B','SEC61G', 'HSPA5'\n )\nprots <- NULL\nsbst <- make_dat_for_PCplot(prots, rnas, pca_color)\n\nhead(sbst, n = n_distinct(sbst$Name) ) ","metadata":{"execution":{"iopub.status.busy":"2023-04-16T05:48:20.537640Z","iopub.execute_input":"2023-04-16T05:48:20.539355Z","iopub.status.idle":"2023-04-16T05:48:22.979337Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(22,6)\nprint_plots(sbst)","metadata":{"execution":{"iopub.status.busy":"2023-04-16T05:48:22.983587Z","iopub.execute_input":"2023-04-16T05:48:22.985293Z","iopub.status.idle":"2023-04-16T05:49:00.237056Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Examine change by day with animated plot - in the output (GIF)**","metadata":{}},{"cell_type":"code","source":"prot_of_interest <- \"ENSG00000172590_MRPL52\"\n\np1 <- ggplot(sbst[sbst$Name == prot_of_interest, ],\n aes(x = PC_1, y = PC_2, color = `Expression\\nto median`)) +\n geom_point(size = 0.5) +\n theme_classic(base_size = 18) +\n scale_color_manual(values=c('yellow', 'darkblue')) +\n #scale_color_gradient(low='darkblue', high='yellow') +\n theme(aspect.ratio = 1) + #, legend.position=\"none\"\n # Animating the plot\n labs(title = paste0(prot_of_interest),\n subtitle = 'Day: {current_frame}') +\n transition_manual(ss_per_day, cumulative = FALSE)\nanimate(p1, fps = 30)\n\nfile_name <- paste0(path_gif, gsub(\".*_\",\"\", prot_of_interest), '_by_Day_PCA.gif')\nanim_save(file.path(normalizePath(dirname(file_name)),basename(file_name)))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:49:00.239583Z","iopub.execute_input":"2023-04-16T05:49:00.241063Z","iopub.status.idle":"2023-04-16T05:49:01.586084Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"p1 <- ggplot(sbst[sbst$Name == 'Mean_RNAs', ],\n aes(x = PC_1, y = PC_2, color = `Expression\\nto median`)) +\n geom_point(size = 0.5) +\n theme_classic(base_size = 18) +\n scale_color_manual(values=c('yellow', 'darkblue')) +\n #scale_color_gradient(low='darkblue', high='yellow') +\n theme(aspect.ratio = 1) + #, legend.position=\"none\"\n # Animating the plot\n labs(title = 'Mean expression\\nof mitochondrial translation-related RNAs',\n subtitle = 'Day: {current_frame}') +\n transition_manual(ss_per_day, cumulative = FALSE)\nanimate(p1, fps = 30)\n\nfile_name <- paste0(path_gif, 'Mean_Mitoch_by_Day_PCA.gif')\nanim_save(file.path(normalizePath(dirname(file_name)),basename(file_name)))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:49:01.589392Z","iopub.execute_input":"2023-04-16T05:49:01.592331Z","iopub.status.idle":"2023-04-16T05:49:02.935868Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Examine individual RNAs in the 3D plot - also in the output (HTML format)**","metadata":{}},{"cell_type":"code","source":"fig(15,20) \nprot_of_interest <- \"ENSG00000106803_SEC61B\" #ENSG00000106803_SEC61B #ENSG00000178952_TUFM\nplt <- plot_ly(sample_n(sbst[sbst$Name == prot_of_interest,], 5000), x = ~PC_1, y = ~PC_2, z = ~Expression, color = ~cell_type,\n size = 0.5) %>% \n add_markers() %>%\n layout(title = paste(prot_of_interest, \"expression vs first two principal components (sample = 5000)\"))\n#plt\nfile_name = paste0(path_3d, gsub(\".*_\", \"\", prot_of_interest), \"_PCA3d.HTML\")\nhtmlwidgets::saveWidget(partial_bundle(plt), file.path(normalizePath(dirname(file_name)),basename(file_name)))\ndisplay_html(paste0(''))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:49:02.939790Z","iopub.execute_input":"2023-04-16T05:49:02.941659Z","iopub.status.idle":"2023-04-16T05:49:03.964897Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Genes, regulating splicing","metadata":{}},{"cell_type":"markdown","source":"Here we plot RNA, top negatively correlated with CD45 RNA and forming the enriched terms related to splicing, transcription, translation ([found here](http://www.kaggle.com/code/antoninadolgorukova/mmscel-cd45-iso-go-enrichment#4.-Translation,-transcription,-splicing)). Also, SNRPE, FOXP1, SNRPD2 are the components of slicesiosome. \nExpression and correlation with other proteins and RNAs is examined [here](http://www.kaggle.com/code/antoninadolgorukova/cd45-reg-related-rna-in-2-citeseq-datasets#Splicing-genes-from-GO-terms-(CD45-isoforms))","metadata":{}},{"cell_type":"code","source":"rnas <- c('SNRPD2', 'POLR2E', 'SNRNP25', 'TXNL4A', 'SNRNP40',\n 'SNRPE', 'POLR2I', 'SNRPD3', 'POLR2F', 'SNRPG',\n 'SF3B5', 'POLR2L', 'SNU13', 'SNRPB', 'SRSF7',\n 'SRSF2', 'SNRPD1', 'SNRPF', 'YBX1', 'SRRM1',\n 'SNRNP70', 'PPIH', 'LSM3', 'HNRNPF', 'LSM7',\n 'SNRPA1', 'ALYREF', 'STRAP', 'SRSF3', 'LSM5',\n 'SNRPC', 'DDX39A', 'UBL5', 'SRPK1', 'HSPA8', 'C1QBP'\n )\nprots <- NULL\nsbst <- make_dat_for_PCplot(prots, rnas, pca_color)\n\nhead(sbst, n = 5) #n_distinct(sbst$Name) ","metadata":{"execution":{"iopub.status.busy":"2023-04-16T05:49:03.968497Z","iopub.execute_input":"2023-04-16T05:49:03.970151Z","iopub.status.idle":"2023-04-16T05:49:18.522144Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(22,6)\n\nsbst[sbst$Name == \"Mean_RNAs\", ]$Name <- \"Mean_GO_splicing_genes\"\n\nprint_plots(sbst[sbst$Name %in% c(\n\"ENSG00000108561_C1QBP\",\n\"ENSG00000143977_SNRPG\",\n\"ENSG00000125743_SNRPD2\",\n\"ENSG00000130332_LSM7\",\n\"ENSG00000100142_POLR2F\",\n\"ENSG00000169976_SF3B5\",\n\"ENSG00000065978_YBX1\",\n\"Mean_GO_splicing_genes\"), ])","metadata":{"execution":{"iopub.status.busy":"2023-04-16T05:49:18.524614Z","iopub.execute_input":"2023-04-16T05:49:18.526051Z","iopub.status.idle":"2023-04-16T05:49:50.231262Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Examine change by day with animated plot - in the output (GIF)**","metadata":{}},{"cell_type":"code","source":"prot_of_interest <- \"ENSG00000108561_C1QBP\"\n\np1 <- ggplot(sbst[sbst$Name == prot_of_interest, ],\n aes(x = PC_1, y = PC_2, color = `Expression\\nto median`)) +\n geom_point(size = 0.5) +\n theme_classic(base_size = 18) +\n scale_color_manual(values=c('yellow', 'darkblue')) +\n #scale_color_gradient(low='darkblue', high='yellow') +\n theme(aspect.ratio = 1) + #, legend.position=\"none\"\n # Animating the plot\n labs(title = paste0(prot_of_interest),\n subtitle = 'Day: {current_frame}') +\n transition_manual(ss_per_day, cumulative = FALSE)\nanimate(p1, fps = 30)\n\nfile_name <- paste0(path_gif, gsub(\".*_\",\"\", prot_of_interest), '_by_Day_PCA.gif')\nanim_save(file.path(normalizePath(dirname(file_name)),basename(file_name)))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:49:50.233894Z","iopub.execute_input":"2023-04-16T05:49:50.235482Z","iopub.status.idle":"2023-04-16T05:49:51.589942Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Examine individual RNAs in the 3D plot - also in the output (HTML format)**","metadata":{}},{"cell_type":"code","source":"fig(15,20) \nprot_of_interest <- \"ENSG00000108561_C1QBP\" #\nplt <- plot_ly(sample_n(sbst[sbst$Name == prot_of_interest,], 5000), x = ~PC_1, y = ~PC_2, z = ~Expression, color = ~cell_type,\n size = 0.5) %>% \n add_markers() %>%\n layout(title = paste(prot_of_interest, \"expression vs first two principal components (sample = 5000)\"))\n#plt\nfile_name = paste0(path_3d, gsub(\".*_\", \"\", prot_of_interest), \"_PCA3d.HTML\")\nhtmlwidgets::saveWidget(partial_bundle(plt), file.path(normalizePath(dirname(file_name)),basename(file_name)))\ndisplay_html(paste0(''))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:49:51.593807Z","iopub.execute_input":"2023-04-16T05:49:51.595610Z","iopub.status.idle":"2023-04-16T05:49:52.591023Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# HNRNPLL and related","metadata":{}},{"cell_type":"markdown","source":"We use the pattern to find all RNA in the data set with HNRNP in their name","metadata":{}},{"cell_type":"code","source":"f_rnas <- c(colnames(mat_RNA)[grepl(\"HNRNP\",\n colnames(mat_RNA),\n ignore.case = TRUE)])\ncat(\"Found\", f_rnas %>% length, \"genes starting with HNRNP\")\ngenes <- mat_RNA[, f_rnas] %>% as.data.frame\n\n# correlation matrix\ncor_mat <- cor(genes, method = \"spearman\")\n\ncat(\"\\nCorrelation matix for all found HNRNP genes with non-zero expression\")\nhead(cor_mat, 5)","metadata":{"execution":{"iopub.status.busy":"2023-04-16T05:49:52.594831Z","iopub.execute_input":"2023-04-16T05:49:52.596537Z","iopub.status.idle":"2023-04-16T05:49:53.269308Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25, 10)\nsum_expr_RNA <- apply(genes, 2, median) %>% sort(decreasing = TRUE)\nsum_expr_RNA <- data.frame(\"Median_expression\" = sum_expr_RNA,\n \"Name\" = names(sum_expr_RNA))\nsum_expr_RNA$Name <- factor(sum_expr_RNA$Name, levels = unique(sum_expr_RNA$Name ))\n#hline <- sum_expr_RNA[sum_expr_RNA$Name == \"PTPRC\",]$Median_expression\n\nggplot(sum_expr_RNA, aes(x = Name, y = Median_expression)) +\n geom_bar(stat = 'identity', fill=\"#BF9039\", col=\"grey\") + \n# geom_hline(yintercept=sum_expr_RNA23[sum_expr_RNA23$Name == \"PTPRC\",]$Median_expression,\n# linetype=\"dashed\")+\n theme_bw(base_size = 22) +\n xlab(\"\") +\n ggtitle(paste0(\"CITEseq 2022\")) +\n theme(axis.text.x = element_text(angle = 45, hjust=1))\n\ncat(\"Zero median expression:\\n\")\ncat(paste0('\"', sum_expr_RNA[sum_expr_RNA$Median_expression == 0,]$Name, '\"', collapse = \", \"))","metadata":{"execution":{"iopub.status.busy":"2023-04-16T05:49:53.271995Z","iopub.execute_input":"2023-04-16T05:49:53.273526Z","iopub.status.idle":"2023-04-16T05:49:53.913709Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"not_plot <- c(\"ENSG00000139675_HNRNPA1L2\", \"ENSG00000214223_HNRNPA1P10\", \"ENSG00000227638_HNRNPA1P14\",\n \"ENSG00000262333_HNRNPA1P16\", \"ENSG00000213049_HNRNPA1P34\", \"ENSG00000212961_HNRNPA1P40\",\n \"ENSG00000259512_HNRNPA1P5\", \"ENSG00000257195_HNRNPA1P50\", \"ENSG00000236539_HNRNPA1P54\",\n \"ENSG00000230946_HNRNPA1P68\", \"ENSG00000255141_HNRNPA1P76\", \"ENSG00000260689_HNRNPA3P11\",\n \"ENSG00000227688_HNRNPA3P2\", \"ENSG00000213300_HNRNPA3P6\", \"ENSG00000258900_HNRNPCP1\",\n \"ENSG00000204253_HNRNPCP2\", \"ENSG00000263179_HNRNPCP4\", \"ENSG00000213305_HNRNPCP6\",\n \"ENSG00000220305_HNRNPH1P1\", \"ENSG00000227347_HNRNPKP2\", \"ENSG00000243547_HNRNPKP4\",\n \"ENSG00000259335_HNRNPMP1\", \"ENSG00000223984_HNRNPRP1\")\n\ncor_mat <- cor_mat[!rownames(cor_mat) %in% not_plot,\n !colnames(cor_mat) %in% not_plot]\n \nfig(10, 10)\nall_heat <- pheatmap(cor_mat, fontsize = 14,\n color = myColors,\n breaks = breaksList,\n angle_col = 90,\n main = paste0(nrow(cor_mat), \" HNRNP** genes with non-zero expression\"))","metadata":{"execution":{"iopub.status.busy":"2023-04-16T05:49:53.916618Z","iopub.execute_input":"2023-04-16T05:49:53.918139Z","iopub.status.idle":"2023-04-16T05:49:54.150452Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We select genes with median expression above 0.","metadata":{}},{"cell_type":"code","source":"rnas <- c(\"HNRNPA0\", \"HNRNPA1\", \"HNRNPA1P48\", \"HNRNPA2B1\",\n \"HNRNPA3\", \"HNRNPAB\", \"HNRNPC\", \"HNRNPD\", \"HNRNPDL\",\n \"HNRNPF\", \"HNRNPH1\", \"HNRNPH2\", \"HNRNPH3\", \"HNRNPK\",\n \"HNRNPL\", \"HNRNPLL\", \"HNRNPM\", \"HNRNPR\", \"HNRNPU\",\n \"HNRNPUL1\", \"HNRNPUL2\")\nprots <- NULL\nsbst <- make_dat_for_PCplot(prots, rnas, pca_color)\n\nhead(sbst, n = 5) #n_distinct(sbst$Name)","metadata":{"execution":{"iopub.status.busy":"2023-04-16T05:49:54.153095Z","iopub.execute_input":"2023-04-16T05:49:54.154577Z","iopub.status.idle":"2023-04-16T05:50:01.179774Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(22,6)\nsbst[sbst$Name == \"Mean_RNAs\", ]$Name <- \"Mean_HNRNPs\"\nprint_plots(sbst[sbst$Name %in% c(\n \"ENSG00000135486_HNRNPA1\",\n\"ENSG00000122566_HNRNPA2B1\",\n\"ENSG00000197451_HNRNPAB\",\n\"ENSG00000169045_HNRNPH1\",\n\"ENSG00000143889_HNRNPLL\",\n\"ENSG00000169813_HNRNPF\",\n\"Mean_HNRNPs\"), ])","metadata":{"execution":{"iopub.status.busy":"2023-04-16T05:50:01.182158Z","iopub.execute_input":"2023-04-16T05:50:01.183576Z","iopub.status.idle":"2023-04-16T05:50:29.355513Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Examine change by day with animated plot - in the output (GIF)**","metadata":{}},{"cell_type":"code","source":"prot_of_interest <- \"ENSG00000135486_HNRNPA1\"\n\np1 <- ggplot(sbst[sbst$Name == prot_of_interest, ],\n aes(x = PC_1, y = PC_2, color = `Expression\\nto median`)) +\n geom_point(size = 0.5) +\n theme_classic(base_size = 18) +\n scale_color_manual(values=c('yellow', 'darkblue')) +\n #scale_color_gradient(low='darkblue', high='yellow') +\n theme(aspect.ratio = 1) + #, legend.position=\"none\"\n # Animating the plot\n labs(title = paste0(prot_of_interest),\n subtitle = 'Day: {current_frame}') +\n transition_manual(ss_per_day, cumulative = FALSE)\nanimate(p1, fps = 30)\n\nfile_name <- paste0(path_gif, gsub(\".*_\",\"\", prot_of_interest), '_by_Day_PCA.gif')\nanim_save(file.path(normalizePath(dirname(file_name)),basename(file_name)))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:50:29.358811Z","iopub.execute_input":"2023-04-16T05:50:29.360655Z","iopub.status.idle":"2023-04-16T05:50:30.775421Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"prot_of_interest <- \"ENSG00000143889_HNRNPLL\"\n\np1 <- ggplot(sbst[sbst$Name == prot_of_interest, ],\n aes(x = PC_1, y = PC_2, color = `Expression\\nto median`)) +\n geom_point(size = 0.5) +\n theme_classic(base_size = 18) +\n scale_color_manual(values=c('yellow', 'darkblue')) +\n #scale_color_gradient(low='darkblue', high='yellow') +\n theme(aspect.ratio = 1) + #, legend.position=\"none\"\n # Animating the plot\n labs(title = paste0(prot_of_interest),\n subtitle = 'Day: {current_frame}') +\n transition_manual(ss_per_day, cumulative = FALSE)\nanimate(p1, fps = 30)\n\nfile_name <- paste0(path_gif, gsub(\".*_\",\"\", prot_of_interest), '_by_Day_PCA.gif')\nanim_save(file.path(normalizePath(dirname(file_name)),basename(file_name)))","metadata":{"execution":{"iopub.status.busy":"2023-04-16T05:50:30.779428Z","iopub.execute_input":"2023-04-16T05:50:30.781324Z","iopub.status.idle":"2023-04-16T05:50:32.157035Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Examine individual RNAs in the 3D plot - also in the output (HTML format)**","metadata":{}},{"cell_type":"code","source":"fig(15,20)\nprot_of_interest <- \"ENSG00000143889_HNRNPLL\" \nplt <- plot_ly(sample_n(sbst[sbst$Name == prot_of_interest,], 5000), x = ~PC_1, y = ~PC_2, z = ~Expression, color = ~cell_type,\n size = 0.5) %>% \n add_markers() %>%\n layout(title = paste(prot_of_interest, \"expression vs first two principal components (sample = 5000)\"))\n#plt\nfile_name = paste0(path_3d, gsub(\".*_\", \"\", prot_of_interest), \"_PCA3d.HTML\")\nhtmlwidgets::saveWidget(partial_bundle(plt), file.path(normalizePath(dirname(file_name)),basename(file_name)))\ndisplay_html(paste0(''))","metadata":{"execution":{"iopub.status.busy":"2023-04-16T05:50:32.161239Z","iopub.execute_input":"2023-04-16T05:50:32.163059Z","iopub.status.idle":"2023-04-16T05:50:33.068368Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Next we plot only well-expressed genes that forme a correlated clucter (see [this notebook](http://www.kaggle.com/code/antoninadolgorukova/cd45-reg-related-rna-in-2-citeseq-datasets#Alternative-Splicing-of-CD45))","metadata":{}},{"cell_type":"code","source":"to_plot <- c(\n\n\"HNRNPA0\", \"HNRNPA2B1\", \"HNRNPA3\", \"HNRNPAB\",\n \"HNRNPC\", \"HNRNPD\", \"HNRNPDL\", \"HNRNPF\", \n \"HNRNPH2\", \"HNRNPH3\", \"HNRNPK\", \"HNRNPM\", \"HNRNPR\", \"HNRNPU\"\n\n)\nrnas <- to_plot\nprots <- NULL\nsbst <- make_dat_for_PCplot(prots, rnas, pca_color, exact = TRUE)\n\nhead(sbst, n = 5) #n_distinct(sbst$Name)","metadata":{"execution":{"iopub.status.busy":"2023-04-16T05:50:33.072199Z","iopub.execute_input":"2023-04-16T05:50:33.073869Z","iopub.status.idle":"2023-04-16T05:50:37.116300Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(22,6)\nsbst[sbst$Name == \"Mean_RNAs\", ]$Name <- \"Mean_14HNRNPs_cluster\"\nprint_plots(sbst[sbst$Name == \"Mean_14HNRNPs_cluster\", ])","metadata":{"execution":{"iopub.status.busy":"2023-04-16T05:50:37.119261Z","iopub.execute_input":"2023-04-16T05:50:37.120796Z","iopub.status.idle":"2023-04-16T05:50:41.251978Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Examine change by day with animated plot - in the output (GIF)**","metadata":{}},{"cell_type":"code","source":"prot_of_interest <- \"Mean_14HNRNPs_cluster\"\n\np1 <- ggplot(sbst[sbst$Name == prot_of_interest, ],\n aes(x = PC_1, y = PC_2, color = `Expression\\nto median`)) +\n geom_point(size = 0.5) +\n theme_classic(base_size = 18) +\n scale_color_manual(values=c('yellow', 'darkblue')) +\n #scale_color_gradient(low='darkblue', high='yellow') +\n theme(aspect.ratio = 1) + #, legend.position=\"none\"\n # Animating the plot\n labs(title = paste0('Mean expression of HNRNP** RNAs'),\n subtitle = 'Day: {current_frame}') +\n transition_manual(ss_per_day, cumulative = FALSE)\nanimate(p1, fps = 30)\n\nfile_name <- paste0(path_gif, 'Mean_14HNRNPs_cluster_by_Day_PCA.gif')\nanim_save(file.path(normalizePath(dirname(file_name)),basename(file_name)))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:50:41.254706Z","iopub.execute_input":"2023-04-16T05:50:41.256322Z","iopub.status.idle":"2023-04-16T05:50:42.643254Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Next we separately plot genes, related to HNRNPLL, using the same approach","metadata":{}},{"cell_type":"code","source":"f_rnas <- c(colnames(mat_RNA)[grepl(paste0(\n c('GALM', 'WDR83OS', 'ALK', 'ARL','HAA', 'TEM'), collapse = \"|\"),\n colnames(mat_RNA),\n ignore.case = TRUE)])\ncat(\"Found\", f_rnas %>% length, \"genes\")\ngenes <- mat_RNA[, f_rnas] %>% as.data.frame\n\n# correlation matrix\ncor_mat <- cor(genes, method = \"spearman\")\n\ncat(\"\\nCorrelation matix for all found genes with non-zero expression\")\nhead(cor_mat, 5)","metadata":{"execution":{"iopub.status.busy":"2023-04-16T05:50:42.647294Z","iopub.execute_input":"2023-04-16T05:50:42.649121Z","iopub.status.idle":"2023-04-16T05:50:43.352796Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25, 10)\nsum_expr_RNA <- apply(genes, 2, median) %>% sort(decreasing = TRUE)\nsum_expr_RNA <- data.frame(\"Median_expression\" = sum_expr_RNA,\n \"Name\" = names(sum_expr_RNA))\nsum_expr_RNA$Name <- factor(sum_expr_RNA$Name, levels = unique(sum_expr_RNA$Name ))\n#hline <- sum_expr_RNA[sum_expr_RNA$Name == \"PTPRC\",]$Median_expression\n\nggplot(sum_expr_RNA, aes(x = Name, y = Median_expression)) +\n geom_bar(stat = 'identity', fill=\"#BF9039\", col=\"grey\") + \n# geom_hline(yintercept=sum_expr_RNA23[sum_expr_RNA23$Name == \"PTPRC\",]$Median_expression,\n# linetype=\"dashed\")+\n theme_bw(base_size = 22) +\n xlab(\"\") +\n ggtitle(paste0(\"CITEseq 2022\")) +\n theme(axis.text.x = element_text(angle = 45, hjust=1))\n\ncat(\"Zero median expression:\\n\")\ncat(paste0('\"', sum_expr_RNA[sum_expr_RNA$Median_expression == 0,]$Name, '\"', collapse = \", \"))","metadata":{"execution":{"iopub.status.busy":"2023-04-16T05:50:43.356791Z","iopub.execute_input":"2023-04-16T05:50:43.358597Z","iopub.status.idle":"2023-04-16T05:50:44.026855Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"not_plot <- c(\"ENSG00000189292_ALKAL2\", \"ENSG00000100601_ALKBH1\", \"ENSG00000189046_ALKBH2\",\n \"ENSG00000166199_ALKBH3\", \"ENSG00000244926_ALKBH3-AS1\", \"ENSG00000160993_ALKBH4\",\n \"ENSG00000239382_ALKBH6\", \"ENSG00000137760_ALKBH8\", \"ENSG00000175414_ARL10\",\n \"ENSG00000152213_ARL11\", \"ENSG00000169379_ARL13B\", \"ENSG00000185305_ARL15\",\n \"ENSG00000185829_ARL17A\", \"ENSG00000228696_ARL17B\", \"ENSG00000122644_ARL4A\",\n \"ENSG00000218996_ARL4AP5\", \"ENSG00000188042_ARL4C\", \"ENSG00000175906_ARL4D\",\n \"ENSG00000162980_ARL5A\", \"ENSG00000165997_ARL5B\", \"ENSG00000113966_ARL6\",\n \"ENSG00000143862_ARL8A\", \"ENSG00000196503_ARL9\", \"ENSG00000156958_GALK2\",\n \"ENSG00000143891_GALM\", \"ENSG00000162882_HAAO\", \"ENSG00000105583_WDR83OS\")\n\ncor_mat <- cor_mat[!rownames(cor_mat) %in% not_plot,\n !colnames(cor_mat) %in% not_plot]\n\nfig(10, 10)\nall_heat <- pheatmap(cor_mat, fontsize = 14,\n color = myColors,\n breaks = breaksList,\n angle_col = 90,\n main = paste0(nrow(cor_mat), \" mitochondrial genes with non-zero expression\"))","metadata":{"execution":{"iopub.status.busy":"2023-04-16T05:50:44.029677Z","iopub.execute_input":"2023-04-16T05:50:44.031239Z","iopub.status.idle":"2023-04-16T05:50:44.256724Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rnas <- colnames(cor_mat)[!colnames(cor_mat) %in% not_plot]\nrnas <- gsub(\".*_\", \"\", rnas)\nprots <- NULL\nsbst <- make_dat_for_PCplot(prots, rnas, pca_color, exact = TRUE)\n\nhead(sbst, n = 5) #n_distinct(sbst$Name)","metadata":{"execution":{"iopub.status.busy":"2023-04-16T05:50:44.259548Z","iopub.execute_input":"2023-04-16T05:50:44.261154Z","iopub.status.idle":"2023-04-16T05:50:48.690588Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(22,6)\nprint_plots(sbst[sbst$Name %in% c(\n \"ENSG00000125652_ALKBH7\",\n \"ENSG00000213465_ARL2\",\n \"ENSG00000134108_ARL8B\"\n), ]) #mean makes no sense here, too different genes","metadata":{"execution":{"iopub.status.busy":"2023-04-16T05:50:48.694501Z","iopub.execute_input":"2023-04-16T05:50:48.696231Z","iopub.status.idle":"2023-04-16T05:51:00.761552Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Examine change by day with animated plot - in the output (GIF)**","metadata":{}},{"cell_type":"code","source":"prot_of_interest <- \"ENSG00000213465_ARL2\"\n\np1 <- ggplot(sbst[sbst$Name == prot_of_interest, ],\n aes(x = PC_1, y = PC_2, color = `Expression\\nto median`)) +\n geom_point(size = 0.5) +\n theme_classic(base_size = 18) +\n scale_color_manual(values=c('yellow', 'darkblue')) +\n #scale_color_gradient(low='darkblue', high='yellow') +\n theme(aspect.ratio = 1) + #, legend.position=\"none\"\n # Animating the plot\n labs(title = paste0(prot_of_interest),\n subtitle = 'Day: {current_frame}') +\n transition_manual(ss_per_day, cumulative = FALSE)\nanimate(p1, fps = 30)\n\nfile_name <- paste0(path_gif, gsub(\".*_\",\"\", prot_of_interest), '_by_Day_PCA.gif')\nanim_save(file.path(normalizePath(dirname(file_name)),basename(file_name)))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:51:00.764350Z","iopub.execute_input":"2023-04-16T05:51:00.766005Z","iopub.status.idle":"2023-04-16T05:51:02.193440Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# STAT family","metadata":{}},{"cell_type":"code","source":"rnas <- c('STAT')\nprots <- NULL\nsbst <- make_dat_for_PCplot(prots, rnas, pca_color, exact = FALSE)\n\nhead(sbst, n = n_distinct(sbst$Name) )","metadata":{"_kg_hide-input":false,"execution":{"iopub.status.busy":"2023-04-16T05:51:02.200852Z","iopub.execute_input":"2023-04-16T05:51:02.203117Z","iopub.status.idle":"2023-04-16T05:51:04.256114Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(22,6)\nprint_plots(sbst)","metadata":{"execution":{"iopub.status.busy":"2023-04-16T05:51:04.259727Z","iopub.execute_input":"2023-04-16T05:51:04.261464Z","iopub.status.idle":"2023-04-16T05:51:36.657094Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Examine change by day with animated plot - in the output (GIF)**","metadata":{}},{"cell_type":"code","source":"prot_of_interest <- \"Mean_RNAs\"\n\np1 <- ggplot(sbst[sbst$Name == prot_of_interest, ],\n aes(x = PC_1, y = PC_2, color = `Expression\\nto median`)) +\n geom_point(size = 0.5) +\n theme_classic(base_size = 18) +\n scale_color_manual(values=c('yellow', 'darkblue')) +\n #scale_color_gradient(low='darkblue', high='yellow') +\n theme(aspect.ratio = 1) + #, legend.position=\"none\"\n # Animating the plot\n labs(title = paste0('Mean expression of STAT** RNAs'),\n subtitle = 'Day: {current_frame}') +\n transition_manual(ss_per_day, cumulative = FALSE)\nanimate(p1, fps = 30)\n\nfile_name <- paste0(path_gif, 'Mean_STAT_by_Day_PCA.gif')\nanim_save(file.path(normalizePath(dirname(file_name)),basename(file_name)))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:51:36.659626Z","iopub.execute_input":"2023-04-16T05:51:36.661115Z","iopub.status.idle":"2023-04-16T05:51:38.011819Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# CD53 and related","metadata":{}},{"cell_type":"markdown","source":"Here we plot genes extracted from [the paper](http://doi.org/10.1016/j.celrep.2022.111006) on CD53-CD45 relation.","metadata":{}},{"cell_type":"code","source":"rnas <- c('CD53', 'SP1','MUC1','STAR', 'LCK'\n #'APC', 'BCR', # evenly distributed\n #'SDS','LCK','NHS','ZAP70' #empty \n )\nprots <- NULL\nsbst <- make_dat_for_PCplot(prots, rnas, pca_color)\n\nhead(sbst, n = n_distinct(sbst$Name) )","metadata":{"execution":{"iopub.status.busy":"2023-04-16T05:51:38.016570Z","iopub.execute_input":"2023-04-16T05:51:38.018534Z","iopub.status.idle":"2023-04-16T05:51:39.630072Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(22,6)\nprint_plots(sbst[!sbst$Name == \"Mean_RNAs\", ]) #mean makes no sense here, too different genes","metadata":{"_kg_hide-input":false,"execution":{"iopub.status.busy":"2023-04-16T05:51:39.633918Z","iopub.execute_input":"2023-04-16T05:51:39.635771Z","iopub.status.idle":"2023-04-16T05:52:01.664050Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Examine change by day with animated plot - in the output (GIF)**","metadata":{}},{"cell_type":"code","source":"prot_of_interest <- \"ENSG00000143119_CD53\"\n\np1 <- ggplot(sbst[sbst$Name == prot_of_interest, ],\n aes(x = PC_1, y = PC_2, color = `Expression\\nto median`)) +\n geom_point(size = 0.5) +\n theme_classic(base_size = 18) +\n scale_color_manual(values=c('yellow', 'darkblue')) +\n #scale_color_gradient(low='darkblue', high='yellow') +\n theme(aspect.ratio = 1) + #, legend.position=\"none\"\n # Animating the plot\n labs(title = paste0(prot_of_interest),\n subtitle = 'Day: {current_frame}') +\n transition_manual(ss_per_day, cumulative = FALSE)\nanimate(p1, fps = 30)\n\nfile_name <- paste0(path_gif, gsub(\".*_\",\"\", prot_of_interest), '_by_Day_PCA.gif')\nanim_save(file.path(normalizePath(dirname(file_name)),basename(file_name)))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:52:01.672688Z","iopub.execute_input":"2023-04-16T05:52:01.676623Z","iopub.status.idle":"2023-04-16T05:52:03.835898Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Examine individual RNAs in the 3D plot - also in the output (HTML format)**","metadata":{}},{"cell_type":"code","source":"fig(15,20) \nprot_of_interest <- \"ENSG00000143119_CD53\"\nplt <- plot_ly(sample_n(sbst[sbst$Name == prot_of_interest,], 5000), x = ~PC_1, y = ~PC_2, z = ~Expression, color = ~cell_type,\n size = 0.5) %>% \n add_markers() %>%\n layout(title = paste(prot_of_interest, \"expression vs first two principal components (sample = 5000)\"))\n#plt\nfile_name = paste0(path_3d, gsub(\".*_\", \"\", prot_of_interest), \"_PCA3d.HTML\")\nhtmlwidgets::saveWidget(partial_bundle(plt), file.path(normalizePath(dirname(file_name)),basename(file_name)))\ndisplay_html(paste0(''))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:52:03.840181Z","iopub.execute_input":"2023-04-16T05:52:03.842135Z","iopub.status.idle":"2023-04-16T05:52:05.503728Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# MALAT1, NEAT1","metadata":{}},{"cell_type":"code","source":"rnas <- c('MALAT1','NEAT1'\n )\nprots <- NULL\nsbst <- make_dat_for_PCplot(prots, rnas, pca_color)\n\nhead(sbst, n = n_distinct(sbst$Name) )","metadata":{"execution":{"iopub.status.busy":"2023-04-16T05:52:05.512282Z","iopub.execute_input":"2023-04-16T05:52:05.517909Z","iopub.status.idle":"2023-04-16T05:52:07.160702Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(22,6)\nprint_plots(sbst, da = \"3\", do = unique(sbst$donor))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Examine change by day with animated plot - in the output (GIF)**","metadata":{}},{"cell_type":"code","source":"prot_of_interest <- \"ENSG00000245532_NEAT1\"\n\np1 <- ggplot(sbst[sbst$Name == prot_of_interest, ],\n aes(x = PC_1, y = PC_2, color = `Expression\\nto median`)) +\n geom_point(size = 0.5) +\n theme_classic(base_size = 18) +\n scale_color_manual(values=c('yellow', 'darkblue')) +\n #scale_color_gradient(low='darkblue', high='yellow') +\n theme(aspect.ratio = 1) + #, legend.position=\"none\"\n # Animating the plot\n labs(title = paste0(prot_of_interest),\n subtitle = 'Day: {current_frame}') +\n transition_manual(ss_per_day, cumulative = FALSE)\nanimate(p1, fps = 30)\n\nfile_name <- paste0(path_gif, gsub(\".*_\",\"\", prot_of_interest), '_by_Day_PCA.gif')\nanim_save(file.path(normalizePath(dirname(file_name)),basename(file_name)))","metadata":{"_kg_hide-output":false,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:52:28.111804Z","iopub.execute_input":"2023-04-16T05:52:28.113464Z","iopub.status.idle":"2023-04-16T05:52:29.490262Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Examine individual RNAs in the 3D plot - also in the output (HTML format)**","metadata":{}},{"cell_type":"code","source":"fig(15,20) \nprot_of_interest <- \"ENSG00000245532_NEAT1\"\nplt <- plot_ly(sample_n(sbst[sbst$Name == prot_of_interest,], 5000), x = ~PC_1, y = ~PC_2, z = ~Expression, color = ~cell_type,\n size = 0.5) %>% \n add_markers() %>%\n layout(title = paste(prot_of_interest, \"expression vs first two principal components (sample = 5000)\"))\n#plt\nfile_name = paste0(path_3d, gsub(\".*_\", \"\", prot_of_interest), \"_PCA3d.HTML\")\nhtmlwidgets::saveWidget(partial_bundle(plt), file.path(normalizePath(dirname(file_name)),basename(file_name)))\ndisplay_html(paste0(''))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:52:29.494524Z","iopub.execute_input":"2023-04-16T05:52:29.496446Z","iopub.status.idle":"2023-04-16T05:52:30.521314Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Proteins","metadata":{}},{"cell_type":"code","source":"rnas <- NULL\nprots <- c('CD3E','CD3D','CD3G',\n #'CD3','CD4', # evenly distributed\n 'CD28', 'CD37','CD38', 'CD69', 'LCK', 'CD44','CD48')\nsbst <- make_dat_for_PCplot(prots, rnas, pca_color)\n\nhead(sbst, n = n_distinct(sbst$Name) ) ","metadata":{"execution":{"iopub.status.busy":"2023-04-16T05:52:30.523885Z","iopub.execute_input":"2023-04-16T05:52:30.525370Z","iopub.status.idle":"2023-04-16T05:52:31.999672Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(22,6)\nprint_plots(sbst)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Examine change by day with animated plot - in the output (GIF)**","metadata":{}},{"cell_type":"code","source":"prot_of_interest <- \"CD38\"\n\np1 <- ggplot(sbst[sbst$Name == prot_of_interest, ],\n aes(x = PC_1, y = PC_2, color = `Expression\\nto median`)) +\n geom_point(size = 0.5) +\n theme_classic(base_size = 18) +\n scale_color_manual(values=c('yellow', 'darkblue')) +\n #scale_color_gradient(low='darkblue', high='yellow') +\n theme(aspect.ratio = 1) + #, legend.position=\"none\"\n # Animating the plot\n labs(title = paste0(prot_of_interest),\n subtitle = 'Day: {current_frame}') +\n transition_manual(ss_per_day, cumulative = FALSE)\nanimate(p1, fps = 30)\n\nfile_name <- paste0(path_gif, gsub(\".*_\",\"\", prot_of_interest), '_by_Day_PCA.gif')\nanim_save(file.path(normalizePath(dirname(file_name)),basename(file_name)))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:52:56.563781Z","iopub.execute_input":"2023-04-16T05:52:56.565721Z","iopub.status.idle":"2023-04-16T05:52:57.918893Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Examine individual RNAs in the 3D plot - also in the output (HTML format)**","metadata":{}},{"cell_type":"code","source":"fig(15,20) \nprot_of_interest <- \"CD38\"\nplt <- plot_ly(sample_n(sbst[sbst$Name == prot_of_interest,], 5000), x = ~PC_1, y = ~PC_2, z = ~Expression, color = ~cell_type,\n size = 0.5) %>% \n add_markers() %>%\n layout(title = paste(prot_of_interest, \"expression vs first two principal components (sample = 5000)\"))\n#plt\nfile_name = paste0(path_3d, gsub(\".*_\", \"\", prot_of_interest), \"_PCA3d.HTML\")\nhtmlwidgets::saveWidget(partial_bundle(plt), file.path(normalizePath(dirname(file_name)),basename(file_name)))\ndisplay_html(paste0(''))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:52:57.924109Z","iopub.execute_input":"2023-04-16T05:52:57.926059Z","iopub.status.idle":"2023-04-16T05:52:59.020299Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Ubiquitin–proteasome system","metadata":{}},{"cell_type":"markdown","source":"Here we examine 39 gene with non-zero expression identified in [this notebook](http://www.kaggle.com/code/antoninadolgorukova/cd45-reg-related-rna-in-2-citeseq-datasets).","metadata":{}},{"cell_type":"code","source":"rnas <- c(\"CUEDC2\", \"OTUB1\", \"OTUD6B-AS1\",\n \"RPS27A\", \"RPS27AP16\", \"UBA1\", \"UBA2\", \"UBA3\",\n \"UBA5\", \"UBA52\", \"UBA6\", \"UBAC1\", \"UBAC2\", \"UBALD2\",\n \"UBAP2\", \"UBAP2L\", \"UBB\", \"UBC\", \"UCHL3\", \"UCHL5\",\n \"USP1\", \"USP10\", \"USP11\", \"USP14\", \"USP15\", \"USP16\",\n \"USP22\", \"USP3\", \"USP33\", \"USP34\", \"USP39\", \"USP4\",\n \"USP47\", \"USP48\", \"USP5\", \"USP7\", \"USP8\", \"USP9X\"\n )\nprots <- NULL\nsbst <- make_dat_for_PCplot(prots, rnas, pca_color)\n\nhead(sbst, n = 5 ) ","metadata":{"execution":{"iopub.status.busy":"2023-04-16T05:52:59.025727Z","iopub.execute_input":"2023-04-16T05:52:59.028128Z","iopub.status.idle":"2023-04-16T05:53:16.464857Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(22,6)\nprint_plots(sbst[!sbst$Name %in% c(\n\"ENSG00000167770_OTUB1\",\n\"ENSG00000253738_OTUD6B-AS1\",\n\"ENSG00000224631_RPS27AP16\",\n\"ENSG00000130985_UBA1\",\n \"ENSG00000144744_UBA3\",\n\"ENSG00000081307_UBA5\",\n\"ENSG00000134882_UBAC2\",\n\"ENSG00000118939_UCHL3\",\n\"ENSG00000143569_UBAP2L\",\n\"ENSG00000116750_UCHL5\",\n\"ENSG00000162607_USP1\",\n\"ENSG00000103194_USP10\",\n\"ENSG00000102226_USP11\",\n\"ENSG00000101557_USP14\",\n\"ENSG00000135655_USP15\",\n\"ENSG00000156256_USP16\",\n\"ENSG00000124422_USP22\",\n\"ENSG00000077254_USP33\",\n\"ENSG00000168883_USP39\",\n\"ENSG00000114316_USP4\",\n\"ENSG00000111667_USP5\",\n\"ENSG00000187555_USP7\",\n\"ENSG00000138592_USP8\",\n\"ENSG00000124486_USP9X\"\n), ])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Examine change by day with animated plot - in the output (GIF)**","metadata":{}},{"cell_type":"code","source":"prot_of_interest <- \"ENSG00000143947_RPS27A\" #ENSG00000221983_UBA52\n\np1 <- ggplot(sbst[sbst$Name == prot_of_interest, ],\n aes(x = PC_1, y = PC_2, color = `Expression\\n(corrected)`)) +\n geom_point(size = 0.5) +\n theme_classic(base_size = 18) +\n #scale_color_manual(values=c('yellow', 'darkblue')) +\n scale_color_gradient(low='darkblue', high='yellow') +\n theme(aspect.ratio = 1) + #, legend.position=\"none\"\n # Animating the plot\n labs(title = paste0(prot_of_interest),\n subtitle = 'Day: {current_frame}') +\n transition_manual(ss_per_day, cumulative = FALSE)\nanimate(p1, fps = 30)\n\nfile_name <- paste0(path_gif, gsub(\".*_\",\"\", prot_of_interest), '_cor_by_Day_PCA.gif')\nanim_save(file.path(normalizePath(dirname(file_name)),basename(file_name)))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:54:19.795079Z","iopub.execute_input":"2023-04-16T05:54:19.796705Z","iopub.status.idle":"2023-04-16T05:54:21.224318Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"prot_of_interest <- \"ENSG00000221983_UBA52\" #\n\np1 <- ggplot(sbst[sbst$Name == prot_of_interest, ],\n aes(x = PC_1, y = PC_2, color = `Expression\\n(corrected)`)) +\n geom_point(size = 0.5) +\n theme_classic(base_size = 18) +\n #scale_color_manual(values=c('yellow', 'darkblue')) +\n scale_color_gradient(low='darkblue', high='yellow') +\n theme(aspect.ratio = 1) + #, legend.position=\"none\"\n # Animating the plot\n labs(title = paste0(prot_of_interest),\n subtitle = 'Day: {current_frame}') +\n transition_manual(ss_per_day, cumulative = FALSE)\nanimate(p1, fps = 30)\n\nfile_name <- paste0(path_gif, gsub(\".*_\",\"\", prot_of_interest), '_cor_by_Day_PCA.gif')\nanim_save(file.path(normalizePath(dirname(file_name)),basename(file_name)))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:54:21.228120Z","iopub.execute_input":"2023-04-16T05:54:21.230192Z","iopub.status.idle":"2023-04-16T05:54:22.763633Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Examine individual RNAs in the 3D plot - also in the output (HTML format)**","metadata":{}},{"cell_type":"code","source":"fig(15,20) \nprot_of_interest <- \"ENSG00000143947_RPS27A\"\nplt <- plot_ly(sample_n(sbst[sbst$Name == prot_of_interest,], 5000), x = ~PC_1, y = ~PC_2, z = ~Expression, color = ~cell_type,\n size = 0.5) %>% \n add_markers() %>%\n layout(title = paste(prot_of_interest, \"expression vs first two principal components (sample = 5000)\"))\n#plt\nfile_name = paste0(path_3d, gsub(\".*_\", \"\", prot_of_interest), \"_PCA3d.HTML\")\nhtmlwidgets::saveWidget(partial_bundle(plt), file.path(normalizePath(dirname(file_name)),basename(file_name)))\ndisplay_html(paste0(''))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:54:22.768166Z","iopub.execute_input":"2023-04-16T05:54:22.770345Z","iopub.status.idle":"2023-04-16T05:54:23.769773Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(15,20) \nprot_of_interest <- \"ENSG00000221983_UBA52\"\nplt <- plot_ly(sample_n(sbst[sbst$Name == prot_of_interest,], 5000), x = ~PC_1, y = ~PC_2, z = ~Expression, color = ~cell_type,\n size = 0.5) %>% \n add_markers() %>%\n layout(title = paste(prot_of_interest, \"expression vs first two principal components (sample = 5000)\"))\n#plt\nfile_name = paste0(path_3d, gsub(\".*_\", \"\", prot_of_interest), \"_PCA3d.HTML\")\nhtmlwidgets::saveWidget(partial_bundle(plt), file.path(normalizePath(dirname(file_name)),basename(file_name)))\ndisplay_html(paste0(''))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:54:23.772651Z","iopub.execute_input":"2023-04-16T05:54:23.774300Z","iopub.status.idle":"2023-04-16T05:54:26.108457Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# RNA editing genes","metadata":{}},{"cell_type":"code","source":"rnas <- c('ADAR', 'APOBEC', 'ADAT', 'RED1', 'DNMT2',\n 'ADARB1', 'ADARB3', 'ADARB4',\n #'ADARB2', #emply\n 'ADARB5', 'ADARB6', 'ADARB7', 'ADARB8',\n 'ADARB9', 'ADARB10', 'ADARB11', 'ADARB12',\n 'ADARB13', 'ADARB14', 'ADARB15', 'ADARB16',\n 'ADARB17', 'ADARB18', 'ADARB19', 'ADARB20',\n 'ADARB21', 'ADARB22', 'ADARB23', 'ADARB24',\n 'ADARB25', 'ADARB26', 'ADARB27', 'ADARB28',\n 'ADARB29', 'ADARB30', 'ADARB31', 'ADARB32',\n 'ADARB33', 'ADARB34', 'ADARB35', 'ADARB36',\n 'ADARB37', 'ADARB38', 'ADARB39', 'ADARB40',\n 'ADARB41', 'ADARB42', 'ADARB43', 'ADARB44',\n 'ADARB45', 'ADARB46', 'ADARB47', 'ADARB48',\n 'ADARB49', 'ADARB50', 'ADARB51')\nprots <- NULL\n\nsbst <- make_dat_for_PCplot(prots, rnas, pca_color)\n\nhead(sbst, n = n_distinct(sbst$Name) ) ","metadata":{"execution":{"iopub.status.busy":"2023-04-16T05:54:26.111143Z","iopub.execute_input":"2023-04-16T05:54:26.112698Z","iopub.status.idle":"2023-04-16T05:54:27.439169Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(22,6)\nsbst[sbst$Name == \"Mean_RNAs\", ]$Name <- \"Mean_RNA_editing\"\nprint_plots(sbst)","metadata":{"_kg_hide-input":false,"execution":{"iopub.status.busy":"2023-04-16T05:54:27.441645Z","iopub.execute_input":"2023-04-16T05:54:27.443206Z","iopub.status.idle":"2023-04-16T05:54:39.124782Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Cell cycle","metadata":{}},{"cell_type":"code","source":"S_phase_genes_Tirosh = c('MCM5', 'PCNA', 'TYMS', 'FEN1', 'MCM2', 'MCM4', 'RRM1',\n 'UNG', 'GINS2', 'MCM6', 'CDCA7', 'DTL', 'PRIM1', 'UHRF1',\n 'MLF1IP', 'HELLS', 'RFC2', 'RPA2', 'NASP', 'RAD51AP1',\n 'GMNN', 'WDR76', 'SLBP', 'CCNE2', 'UBR7', 'POLD3',\n 'MSH2', 'ATAD2', 'RAD51', 'RRM2', 'CDC45', 'CDC6', 'EXO1',\n 'TIPIN', 'DSCC1', 'BLM', 'CASP8AP2', 'USP1', 'CLSPN',\n 'POLA1', 'CHAF1B', 'BRIP1', 'E2F8') \n \nG2_M_genes_Tirosh = c('HMGB2', 'CDK1', 'NUSAP1', 'UBE2C', 'BIRC5', 'TPX2', 'TOP2A',\n 'NDC80', 'CKS2', 'NUF2', 'CKS1B', 'MKI67', 'TMPO', 'CENPF',\n 'TACC3', 'FAM64A', 'SMC4', 'CCNB2', 'CKAP2L', 'CKAP2', 'AURKB',\n 'BUB1', 'KIF11', 'ANP32E', 'TUBB4B', 'GTSE1', 'KIF20B', 'HJURP',\n 'CDCA3', 'HN1', 'CDC20', 'TTK', 'CDC25C', 'KIF2C', 'RANGAP1',\n 'NCAPD2', 'DLGAP5', 'CDCA2', 'CDCA8', 'ECT2', 'KIF23', 'HMMR',\n 'AURKA', 'PSRC1', 'ANLN', 'LBR', 'CKAP5', 'CENPE', 'CTCF',\n 'NEK2', 'G2E3', 'GAS2L3', 'CBX5', 'CENPA')","metadata":{"_kg_hide-input":false,"execution":{"iopub.status.busy":"2023-04-16T05:54:39.127509Z","iopub.execute_input":"2023-04-16T05:54:39.129041Z","iopub.status.idle":"2023-04-16T05:54:39.145066Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"f_rnas <- c(colnames(mat_RNA)[grepl(paste(paste0(\"_\", c(S_phase_genes_Tirosh, G2_M_genes_Tirosh), \"$\"), collapse = \"|\"),\n colnames(mat_RNA),\n ignore.case = TRUE)])\ncat(\"Found\", f_rnas %>% length, \"genes of cell cycle\")\ncc_genes <- mat_RNA[, f_rnas] %>% as.data.frame\n\n#exclude genes with 0 expression in all cells\n#cat(\"\\nZero expression in all cells:\")\n#names(colSums(cc_genes)[colSums(cc_genes) == 0])\n#cc_genes <- select(cc_genes, -c(names(colSums(cc_genes)[colSums(cc_genes) == 0])) )\n\n# correlation matrix\ncc_cor_mat <- cor(cc_genes, method = \"spearman\")\n\ncat(\"\\nCorrelation matix for all found genes with non-zero expression\")\nhead(cc_cor_mat, 5)\n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:54:39.147630Z","iopub.execute_input":"2023-04-16T05:54:39.149036Z","iopub.status.idle":"2023-04-16T05:54:40.904366Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(30, 10)\nsum_expr_RNA <- apply(cc_genes, 2, median) %>% sort(decreasing = TRUE)\nsum_expr_RNA <- data.frame(\"Median_expression\" = sum_expr_RNA,\n \"Name\" = names(sum_expr_RNA))\nsum_expr_RNA$Name <- factor(sum_expr_RNA$Name, levels = unique(sum_expr_RNA$Name ))\n#hline <- sum_expr_RNA[sum_expr_RNA$Name == \"PTPRC\",]$Median_expression\n\nggplot(sum_expr_RNA, aes(x = Name, y = Median_expression)) +\n geom_bar(stat = 'identity', fill=\"#BF9039\", col=\"grey\") + \n# geom_hline(yintercept=sum_expr_RNA23[sum_expr_RNA23$Name == \"PTPRC\",]$Median_expression,\n# linetype=\"dashed\")+\n theme_bw(base_size = 22) +\n xlab(\"\") +\n ggtitle(paste0(\"CITEseq 2022\")) +\n theme(axis.text.x = element_text(angle = 45, hjust=1))\n\ncat(\"Zero median expression:\\n\")\ncat(paste0('\"', sum_expr_RNA[sum_expr_RNA$Median_expression == 0,]$Name, '\"', collapse = \", \"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:54:40.907060Z","iopub.execute_input":"2023-04-16T05:54:40.908614Z","iopub.status.idle":"2023-04-16T05:54:41.789542Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"not_plot <- c(\"ENSG00000011426_ANLN\", \"ENSG00000087586_AURKA\", \"ENSG00000136492_BRIP1\",\n \"ENSG00000175305_CCNE2\", \"ENSG00000158402_CDC25C\", \"ENSG00000184661_CDCA2\",\n \"ENSG00000111665_CDCA3\", \"ENSG00000115163_CENPA\", \"ENSG00000159259_CHAF1B\",\n \"ENSG00000169607_CKAP2L\", \"ENSG00000129173_E2F8\", \"ENSG00000114346_ECT2\",\n \"ENSG00000174371_EXO1\", \"ENSG00000092140_G2E3\", \"ENSG00000139354_GAS2L3\",\n \"ENSG00000075218_GTSE1\", \"ENSG00000123485_HJURP\", \"ENSG00000137807_KIF23\",\n \"ENSG00000142945_KIF2C\", \"ENSG00000080986_NDC80\", \"ENSG00000117650_NEK2\",\n \"ENSG00000101868_POLA1\", \"ENSG00000134222_PSRC1\", \"ENSG00000051180_RAD51\",\n \"ENSG00000112742_TTK\", \"ENSG00000012963_UBR7\")\n\ncc_cor_mat <- cc_cor_mat[!rownames(cc_cor_mat) %in% not_plot,\n !colnames(cc_cor_mat) %in% not_plot]\n# Add annotation\nannotation <- data.frame(genes = rownames(cc_cor_mat)) %>%\n mutate(group = ifelse(gsub(\".*_\", \"\", genes) %in% S_phase_genes_Tirosh, \"S_phase\", \"G2_M_phase\"))\nrownames(annotation) <- rownames(cc_cor_mat) \nannotation <- select(annotation, -c(genes))\n\nfig(25, 18)\nall_cc_heat <- pheatmap(cc_cor_mat, fontsize = 14, annotation = annotation,\n color = myColors,\n breaks = breaksList,\n angle_col = 90,\n main = paste0(nrow(cc_cor_mat), \" mitochondrial genes with non-zero expression\"))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:54:41.792021Z","iopub.execute_input":"2023-04-16T05:54:41.793416Z","iopub.status.idle":"2023-04-16T05:54:42.450637Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rnas <- c(G2_M_genes_Tirosh[!G2_M_genes_Tirosh %in% gsub(\".*_\", \"\", not_plot)])\n\nprots <- NULL\nsbst <- make_dat_for_PCplot(prots, rnas, pca_color)\n\nhead(sbst, n = 5 ) ","metadata":{"_kg_hide-input":false,"execution":{"iopub.status.busy":"2023-04-16T05:54:42.453197Z","iopub.execute_input":"2023-04-16T05:54:42.454668Z","iopub.status.idle":"2023-04-16T05:54:56.119326Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(22,6)\nsbst[sbst$Name == \"Mean_RNAs\", ]$Name <- \"Mean_G2_M_genes_Tirosh\"\nprint_plots(sbst[sbst$Name %in% c(\n \"ENSG00000148773_MKI67\",\n \"Mean_G2_M_genes_Tirosh\"\n), ])","metadata":{"_kg_hide-input":false,"execution":{"iopub.status.busy":"2023-04-16T05:54:56.123128Z","iopub.execute_input":"2023-04-16T05:54:56.124874Z","iopub.status.idle":"2023-04-16T05:55:04.010584Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"prot_of_interest <- \"Mean_G2_M_genes_Tirosh\"\n\np1 <- ggplot(sbst[sbst$Name == prot_of_interest, ],\n aes(x = PC_1, y = PC_2, color = `Expression\\nto median`)) +\n geom_point(size = 0.5) +\n theme_classic(base_size = 18) +\n scale_color_manual(values=c('yellow', 'darkblue')) +\n #scale_color_gradient(low='darkblue', high='yellow') +\n theme(aspect.ratio = 1) + #, legend.position=\"none\"\n # Animating the plot\n labs(title = paste0(prot_of_interest),\n subtitle = 'Day: {current_frame}') + \n transition_manual(ss_per_day, cumulative = FALSE)\nanimate(p1, fps = 30)\n#path_gif, \nfile_name <- paste0(path_gif, gsub(\".*_\",\"\", prot_of_interest), '_by_Day_PCA.gif')\nanim_save(file.path(normalizePath(dirname(file_name)),basename(file_name)))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:55:04.013056Z","iopub.execute_input":"2023-04-16T05:55:04.014545Z","iopub.status.idle":"2023-04-16T05:55:05.404102Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rnas <- c(S_phase_genes_Tirosh[!S_phase_genes_Tirosh %in% gsub(\".*_\", \"\", not_plot)])\n\nprots <- NULL\nsbst <- make_dat_for_PCplot(prots, rnas, pca_color)\n\nhead(sbst, n = 5 ) ","metadata":{"_kg_hide-input":false,"execution":{"iopub.status.busy":"2023-04-16T05:55:05.408703Z","iopub.execute_input":"2023-04-16T05:55:05.410562Z","iopub.status.idle":"2023-04-16T05:55:18.603884Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(22,6)\nsbst[sbst$Name == \"Mean_RNAs\", ]$Name <- \"S_phase_genes_Tirosh\"\nprint_plots(sbst[sbst$Name %in% c(\n \"S_phase_genes_Tirosh\"\n), ])","metadata":{"_kg_hide-input":false,"execution":{"iopub.status.busy":"2023-04-16T05:55:18.606499Z","iopub.execute_input":"2023-04-16T05:55:18.607969Z","iopub.status.idle":"2023-04-16T05:55:22.618726Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"prot_of_interest <- \"S_phase_genes_Tirosh\"\n\np1 <- ggplot(sbst[sbst$Name == prot_of_interest, ],\n aes(x = PC_1, y = PC_2, color = `Expression\\nto median`)) +\n geom_point(size = 0.5) +\n theme_classic(base_size = 18) +\n scale_color_manual(values=c('yellow', 'darkblue')) +\n #scale_color_gradient(low='darkblue', high='yellow') +\n theme(aspect.ratio = 1) + #, legend.position=\"none\"\n # Animating the plot\n labs(title = paste0(prot_of_interest),\n subtitle = 'Day: {current_frame}') + \n transition_manual(ss_per_day, cumulative = FALSE)\nanimate(p1, fps = 30)\n#path_gif, \nfile_name <- paste0(path_gif, gsub(\".*_\",\"\", prot_of_interest), '_by_Day_PCA.gif')\nanim_save(file.path(normalizePath(dirname(file_name)),basename(file_name)))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-16T05:55:22.621308Z","iopub.execute_input":"2023-04-16T05:55:22.622817Z","iopub.status.idle":"2023-04-16T05:55:23.992092Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Your RNAs and proteins","metadata":{}},{"cell_type":"markdown","source":"'name1', 'name2' - RNA names in SYMBOL format ('ZNRF1', not 'ENSG00000186187_ZNRF1') \nUse **exact = FALSE** argument (TRUE is default) in the **make_dat_for_PCplot** function to find all RNA containing the pattern. With the default exact = TRUE you search only exact matches","metadata":{}},{"cell_type":"code","source":"# rnas <- c('name_pattern1', 'name_pattern2') # OR NULL\n# prots <- c('name1','name2','name3') # OR NULL\n\n# sbst <- make_dat_for_PCplot(prots, rnas, pca_color, exact = TRUE)\n\n# head(sbst, n = n_distinct(sbst$Name) ) ","metadata":{"execution":{"iopub.status.busy":"2023-04-16T05:55:23.996086Z","iopub.execute_input":"2023-04-16T05:55:23.997862Z","iopub.status.idle":"2023-04-16T05:55:24.010863Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# fig(22,6)\n# print_plots(sbst)","metadata":{"execution":{"iopub.status.busy":"2023-04-16T05:55:24.014528Z","iopub.execute_input":"2023-04-16T05:55:24.016186Z","iopub.status.idle":"2023-04-16T05:55:24.028657Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# prot_of_interest <- \"full_name1\"\n\n# p1 <- ggplot(sbst[sbst$Name == prot_of_interest, ],\n# aes(x = PC_1, y = PC_2, color = `Expression\\nto median`)) +\n# geom_point(size = 0.5) +\n# theme_classic(base_size = 18) +\n# scale_color_manual(values=c('yellow', 'darkblue')) +\n# #scale_color_gradient(low='darkblue', high='yellow') +\n# theme(aspect.ratio = 1) + #, legend.position=\"none\"\n# # Animating the plot\n# labs(title = paste0(prot_of_interest),\n# subtitle = 'Day: {current_frame}') +\n# transition_manual(ss_per_day, cumulative = FALSE)\n# animate(p1, fps = 30)\n\n# file_name <- paste0(path_gif, gsub(\".*_\",\"\", prot_of_interest), '_by_Day_PCA.gif')\n# anim_save(file.path(normalizePath(dirname(file_name)),basename(file_name)))","metadata":{"execution":{"iopub.status.busy":"2023-04-16T05:55:24.032334Z","iopub.execute_input":"2023-04-16T05:55:24.034017Z","iopub.status.idle":"2023-04-16T05:55:24.048304Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# fig(15,20) \n# prot_of_interest <- \"full_name1\"\n# plt <- plot_ly(sample_n(sbst[sbst$Name == prot_of_interest,], 5000), x = ~PC_1, y = ~PC_2, z = ~Expression, color = ~cell_type,\n# size = 0.5) %>% \n# add_markers() %>%\n# layout(title = paste(prot_of_interest, \"expression vs first two principal components (sample = 5000)\"))\n# #plt\n# file_name = paste0(path_3d, gsub(\".*_\", \"\", prot_of_interest), \"_PCA3d.HTML\")\n# htmlwidgets::saveWidget(partial_bundle(plt), file.path(normalizePath(dirname(file_name)),basename(file_name)))\n# display_html(paste0(''))","metadata":{"_kg_hide-input":false,"execution":{"iopub.status.busy":"2023-04-16T05:55:24.051995Z","iopub.execute_input":"2023-04-16T05:55:24.054301Z","iopub.status.idle":"2023-04-16T05:55:24.068166Z"},"trusted":true},"execution_count":null,"outputs":[]}]}