{"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":"

Differentiation of possible protein and RNA regulation pathways. Part 2

","metadata":{}},{"cell_type":"markdown","source":"## Chapter: CD36","metadata":{}},{"cell_type":"markdown","source":"

What we do here:

\n
\n \n Reference tutorial\n
","metadata":{}},{"cell_type":"markdown","source":"
\n\n

Data

\n\n
\n
    \n
  1. RDS files with sparse matrices of normalised counts data for Open Problems - Multimodal Single-Cell Integration (CITEseq2022). The 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 (we use the train dataset: 3 days, 3 donors, 7 cell types).\n
  2. RDS files with top300 lists of RNA that correlate positively or negatively with proteins and/or RNA of 112 protein-RNA pairs (lists were prepared in the Part 1 notebook).\n\n
\n
","metadata":{}},{"cell_type":"code","source":"#function for figure size adjusment\nfig <- function(width, heigth) {\n options(repr.plot.width = width, repr.plot.height = heigth)\n}\n\noptions(ggrepel.max.overlaps = Inf)\n\nsuppressPackageStartupMessages({\n library(Matrix)\n library(stringr)\n library(ggplot2)\n library(dplyr)\n library(tidyr)\n\n library(patchwork) #arrange plots\n\n BiocManager::install('org.Hs.eg.db')\n library(org.Hs.eg.db)\n\n BiocManager::install('clusterProfiler')\n library(clusterProfiler)\n\n BiocManager::install(\"ggnewscale\")\n library(ggnewscale) \n library(enrichplot)\n\n BiocManager::install(\"ReactomePA\")\n library(\"ReactomePA\")\n \n BiocManager::install(\"GOSim\")\n library(GOSim)\n \n})","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Load data**","metadata":{}},{"cell_type":"code","source":"path <- \"/kaggle/input/sparse-measurement-data-open-problems-multimodal/sp_train_cite_inputs.rds\"\nmat_RNA <- readRDS(path)\n\ncat(\"\\nRNA data:\", ncol(mat_RNA),\n \"genes in columns and\", nrow(mat_RNA), \"cells in rows\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-22T17:19:59.850167Z","iopub.execute_input":"2023-03-22T17:19:59.852001Z","iopub.status.idle":"2023-03-22T17:20:36.908472Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"prot_of_interest <- \"CD36\"","metadata":{"_kg_hide-input":false,"execution":{"iopub.status.busy":"2023-03-22T17:20:36.912226Z","iopub.execute_input":"2023-03-22T17:20:36.913951Z","iopub.status.idle":"2023-03-22T17:20:36.928273Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Positive correlations**","metadata":{}},{"cell_type":"code","source":"n = 300\ndat <- read.csv(paste0(\"/kaggle/input/mmscel-difference-in-prot-rna-correlations-p1/top\", n, \"_pos_norib_mit_EIF_corr.csv\"))\n\npos_dat <- dat[grep(prot_of_interest, dat$Pair), ]","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-22T17:20:36.931689Z","iopub.execute_input":"2023-03-22T17:20:36.933249Z","iopub.status.idle":"2023-03-22T17:20:37.132001Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cat(\"Top\", n, \"lists of RNA that correlate positively with proteins and/or RNA of protein-RNA pairs were matched and separated by:\n\\nCommon RNA - present in both top\", n, \"lists (for the protein and the RNA of the pair)\n\\nUnique for protein - present only in the top\", n, \"list for the protein\n\\nUnique for RNA - present only in the top\", n, \"list for the RNA\")\n\npos_dat %>%\n pivot_wider(id_cols = !corr.coef, names_from = \"RNA_list\",\n values_from = \"RNA\", values_fn = list)\n\ncat(\"Correlation coefficients and sample sizes\")\npos_dat %>% pivot_wider(names_from = \"RNA_list\", id_cols = -RNA,\n values_from = \"corr.coef\",\n values_fn = function(x) paste(round(min(x), 3), \"-\", round(max(x), 3), \", n = \", length(x)) )","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-22T17:20:37.134726Z","iopub.execute_input":"2023-03-22T17:20:37.136327Z","iopub.status.idle":"2023-03-22T17:20:37.207163Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Negative correlations**","metadata":{}},{"cell_type":"code","source":"n = 300\ndat <- read.csv(paste0(\"/kaggle/input/mmscel-difference-in-prot-rna-correlations-p1/top\", n, \"_neg_norib_mit_EIF_corr.csv\"))\n\nneg_dat <- dat[grep(prot_of_interest, dat$Pair), ]","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-22T17:20:37.211197Z","iopub.execute_input":"2023-03-22T17:20:37.212987Z","iopub.status.idle":"2023-03-22T17:20:37.415413Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cat(\"Top\", n, \"lists of RNA that correlate negatively with proteins and/or RNA of protein-RNA pairs were matched and separated by:\n\\nCommon RNA - present in both top\", n, \"lists (for the protein and the RNA of the pair)\n\\nUnique for protein - present only in the top\", n, \"list for the protein\n\\nUnique for RNA - present only in the top\", n, \"list for the RNA\")\nneg_dat %>%\n pivot_wider(id_cols = !corr.coef, names_from = \"RNA_list\",\n values_from = \"RNA\", values_fn = list)\n\ncat(\"Correlation coefficients and sample sizes\")\nneg_dat %>% pivot_wider(names_from = \"RNA_list\", id_cols = -RNA,\n values_from = \"corr.coef\",\n values_fn = function(x) paste(round(min(x), 3), \"-\", round(max(x), 3), \", n = \", length(x)) )","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-22T17:20:37.419751Z","iopub.execute_input":"2023-03-22T17:20:37.421572Z","iopub.status.idle":"2023-03-22T17:20:37.500703Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#For convenience we remove full name of the RNA for the analysis.\npos_dat$Pair <- gsub(\" -.*\", \"\", pos_dat$Pair) #%>% head\nneg_dat$Pair <- gsub(\" -.*\", \"\", neg_dat$Pair) #%>% head","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-22T17:20:37.504759Z","iopub.execute_input":"2023-03-22T17:20:37.507145Z","iopub.status.idle":"2023-03-22T17:20:37.532360Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Make a list af all genes in the dataset with Entrez IDs**","metadata":{}},{"cell_type":"code","source":"geneUniverse <- gsub(\"_.*\",\"\", colnames(mat_RNA)) #ENSEMBL\n#geneUniverse <- gsub(\".*_\",\"\", colnames(mat_RNA)) #SYMBOL\n\n#get Entrez IDs of genes\ngeneUniverse <- AnnotationDbi::select(org.Hs.eg.db, keys = geneUniverse,\n columns = c('ENTREZID'), keytype = 'ENSEMBL') #SYMBOL\ncat(\"All RNA (genes collection): For\", sum(is.na(geneUniverse$ENTREZID)),\"out of\", nrow(geneUniverse), \"genes the Entrez ID was not found\")\n#geneUniverse[is.na(geneUniverse$ENTREZID), ]$ENTREZID\n\ngeneUniverse %>% head","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-22T17:20:37.535017Z","iopub.execute_input":"2023-03-22T17:20:37.536635Z","iopub.status.idle":"2023-03-22T17:20:37.882229Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Add missing Gene IDs manually**","metadata":{}},{"cell_type":"code","source":"gen_id <- read.csv(\"/kaggle/input/research-project-01-around-multimodal-singlecell/convers_TXNAME_GENEID_SYMBOL_ENSG.csv\")\ncat(\"The avalilable data of the data for conversion\")\ngen_id <- gen_id[gen_id$ENSG %in% geneUniverse[is.na(geneUniverse$ENTREZID), ]$ENSEMBL, ]\ngen_id %>% nrow","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-22T17:20:37.884965Z","iopub.execute_input":"2023-03-22T17:20:37.886548Z","iopub.status.idle":"2023-03-22T17:20:38.235358Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for(s in gen_id$ENSEMBL) {\n \n geneUniverse[geneUniverse$ENSEMBL == s, ]$ENTREZID <-\n unique(gen_id[gen_id$ENSG == s, ]$GENEID)\n}\ncat(\"the remaining missimg gene IDs\")\ngeneUniverse[is.na(geneUniverse$ENTREZID), ] %>% nrow","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-22T17:20:38.238053Z","iopub.execute_input":"2023-03-22T17:20:38.239627Z","iopub.status.idle":"2023-03-22T17:20:38.274731Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#

Common RNAs for CD45 and CD45 proteins and their RNA

","metadata":{}},{"cell_type":"markdown","source":"**Select the set of genes for enrichment analysis**","metadata":{}},{"cell_type":"code","source":"cat(\"------Positively correlated RNA------\\n\\n\")\n\ngene_set_pos <- pos_dat[pos_dat$Pair == prot_of_interest & pos_dat$RNA_list == \"Common for prot and RNA\", ] %>% dplyr::select(\"RNA\",\"corr.coef\")\ncat(\"Common RNA for\", prot_of_interest, \"and its RNA:\", nrow(gene_set_pos))\ngene_set_pos$RNA\n\n\ncat(\"------Negatively correlated RNA------\\n\\n\")\n\ngene_set_neg <- neg_dat[neg_dat$Pair == prot_of_interest & neg_dat$RNA_list == \"Common for prot and RNA\", ] %>% dplyr::select(\"RNA\",\"corr.coef\")\ncat(\"Common RNA for\", prot_of_interest, \"and its RNA:\", nrow(gene_set_neg))\ngene_set_neg$RNA","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-22T17:20:38.278630Z","iopub.execute_input":"2023-03-22T17:20:38.280232Z","iopub.status.idle":"2023-03-22T17:20:38.344802Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Get Entrez IDs of genes**","metadata":{}},{"cell_type":"code","source":"#-------------positive-------------\n#get Entrez IDs of genes\ngene_an <- AnnotationDbi::select(org.Hs.eg.db, keys = gene_set_pos$RNA,\n columns = c('ENTREZID'), keytype = 'SYMBOL')\n\ncat(\"For\", sum(is.na(gene_an$ENTREZID)), \"genes the Entrez ID was not found\")\n\n#check for duplicastes\n#gene_an[gene_an$SYMBOL == gene_an$SYMBOL[duplicated(gene_an$SYMBOL)], ]\n#colnames(mat_RNA)[grepl(gene_an$SYMBOL[duplicated(gene_an$SYMBOL)], colnames(mat_RNA))]\ngene_an <- dplyr::filter(gene_an, ENTREZID != \"100187828\")\n\n#join to the rest data\ngene_set_pos <- gene_set_pos %>% left_join(gene_an, by = c(\"RNA\" = \"SYMBOL\"))\ngene_set_pos[is.na(gene_set_pos$ENTREZID), ]\n\n#-------------negative-------------\n#get Entrez IDs of genes\ngene_an <- AnnotationDbi::select(org.Hs.eg.db, keys = gene_set_neg$RNA,\n columns = c('ENTREZID'), keytype = 'SYMBOL')\n\ncat(\"For\", sum(is.na(gene_an$ENTREZID)), \"genes the Entrez ID was not found\")\ngene_an[is.na(gene_an$ENTREZID), ]\n\n#check for duplicastes\n#gene_an[gene_an$SYMBOL == gene_an$SYMBOL[duplicated(gene_an$SYMBOL)], ]\n#colnames(mat_RNA)[grepl(gene_an$SYMBOL[duplicated(gene_an$SYMBOL)], colnames(mat_RNA))]\ngene_an <- dplyr::filter(gene_an, ENTREZID != \"100187828\")\n\n#join to the rest data\ngene_set_neg <- gene_set_neg %>% left_join(gene_an, by = c(\"RNA\" = \"SYMBOL\"))\ngene_set_neg[is.na(gene_set_neg$ENTREZID), ]","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-22T17:20:38.349106Z","iopub.execute_input":"2023-03-22T17:20:38.351060Z","iopub.status.idle":"2023-03-22T17:20:38.800150Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Add missing Gene IDs manually or omit\ngene_set_pos <- na.omit(gene_set_pos)\ngene_set_neg <- na.omit(gene_set_neg)\n\ngeneUniverse <- na.omit(geneUniverse)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-22T17:20:38.804113Z","iopub.execute_input":"2023-03-22T17:20:38.805992Z","iopub.status.idle":"2023-03-22T17:20:38.831271Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Make a list of genes to compare**","metadata":{}},{"cell_type":"code","source":"common <- list(\"CD36_pos\" = gene_set_pos$ENTREZID, \n \"CD36_neg\" = gene_set_neg$ENTREZID)\nstr(common)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-22T17:20:38.835571Z","iopub.execute_input":"2023-03-22T17:20:38.837559Z","iopub.status.idle":"2023-03-22T17:20:38.861870Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#

- GO enrichment over-representation test (biological processes)

","metadata":{}},{"cell_type":"code","source":"comGO <- compareCluster(geneCluster = common, OrgDb = org.Hs.eg.db, fun = enrichGO, ont = \"BP\")\ncomGO <- setReadable(comGO, OrgDb = org.Hs.eg.db, keyType=\"ENTREZID\")\n\ndf_comGO <- comGO %>% as.data.frame\nrownames(df_comGO) <- NULL\n\ncat(\"GO enrichment analysis:\\n\")\ncat(\"\\nCD36_pos:\", nrow(df_comGO[df_comGO$Cluster == \"CD36_pos\",]), \"enriched terms found\",\n \"\\nCD36_neg:\", nrow(df_comGO[df_comGO$Cluster == \"CD36_neg\",]), \"enriched terms found\",\n \"\\nTotal:\", nrow(df_comGO))\n\ndf_comGO$geneID <- str_wrap(df_comGO$geneID, width=30)\ndf_comGO %>% head(10) %>% select(1,3,4,5,7,9)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-22T17:20:38.865873Z","iopub.execute_input":"2023-03-22T17:20:38.867700Z","iopub.status.idle":"2023-03-22T17:21:08.616875Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"comGO <- simplify(comGO, cutoff = 0.7, by = \"p.adjust\", select_fun = min)\ncat(\"The clusterProfiler simplify method reduced GO terms to:\")\ncomGO %>% as.data.frame %>% nrow","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-22T17:21:08.620790Z","iopub.execute_input":"2023-03-22T17:21:08.622455Z","iopub.status.idle":"2023-03-22T17:21:14.726251Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Visualize enriched GO terms**","metadata":{}},{"cell_type":"code","source":"fig(25,20)\nn = 30\ndotplot(comGO, showCategory = n) + \n scale_y_discrete(labels=function(x) str_wrap(x, width=100)) + \n theme_bw(base_size = 18) + plot_annotation(\n title = paste0(\"enrichGO\\nRNA, common for the \", prot_of_interest, \" protein and its RNA\",\n \"\\nEnriched terms shown (per cluster): \",n, \" of \", comGO %>% as.data.frame %>% nrow)) & \ntheme(text = element_text(size = 20) ) ","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-22T17:21:14.728955Z","iopub.execute_input":"2023-03-22T17:21:14.730579Z","iopub.status.idle":"2023-03-22T17:21:15.757937Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25,15)\nn = 30 #comGO %>% as.data.frame %>% nrow\ncnetplot(comGO, showCategory = n)+ \n scale_color_gradient2(name='associated data', low='darkgreen', high='firebrick') + \n theme_void(base_size = 22) + plot_annotation(\n title = paste0(\"enrichGO\\nRNA, common for the \", prot_of_interest, \" protein and its RNA\",\n \"\\nEnriched terms shown: \",n, \" of \", comGO %>% as.data.frame %>% nrow)) & \ntheme(text = element_text(size = 22) ) \n\n#slise top terms in groups by their p.adjust value in the enrichment result\n# plt <- arrange(comGO, p.adjust) %>% \n# group_by(Cluster) %>% \n# slice(1:20)\n\n# plt %>% as.data.frame\n\n#color by R\n#vals <- setNames(gene_set_pos$corr.coef, gene_set_pos$RNA)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-22T17:21:15.762269Z","iopub.execute_input":"2023-03-22T17:21:15.764860Z","iopub.status.idle":"2023-03-22T17:21:23.515261Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#

- Reactome pathway over-representation analysis

","metadata":{}},{"cell_type":"code","source":"comPath <- compareCluster(geneCluster = common, fun = enrichPathway, readable = TRUE)\n\ndf_comPath <- comPath %>% as.data.frame\nrownames(df_comPath) <- NULL\n\ncat(\"GO enrichment analysis:\\n\")\ncat(\"\\nCD36_pos:\", nrow(df_comPath[df_comPath$Cluster == \"CD36_pos\",]), \"enriched terms found\",\n \"\\nCD36_neg:\", nrow(df_comPath[df_comPath$Cluster == \"CD36_neg\",]), \"enriched terms found\",\n \"\\nTotal:\", nrow(df_comPath)\n)\n\ndf_comPath$geneID <- str_wrap(df_comPath$geneID, width=30)\ndf_comPath %>% head(10) %>% select(1,3,4,5,7,9)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-22T17:21:23.518534Z","iopub.execute_input":"2023-03-22T17:21:23.520178Z","iopub.status.idle":"2023-03-22T17:21:42.019221Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Visualize enriched GO terms**","metadata":{}},{"cell_type":"code","source":"fig(25,20)\nn = nrow(df_comPath)\ndotplot(comPath, showCategory=n) + \n scale_y_discrete(labels=function(x) str_wrap(x, width=70)) + \n theme_bw(base_size = 18) + plot_annotation(\n title = paste0(\"enrichPathway\\nRNA, common for the \", prot_of_interest, \" protein and its RNA\",\n \"\\nEnriched terms shown: \",n, \" of \", nrow(df_comPath))) & \ntheme(text = element_text(size = 22) ) ","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-22T17:21:42.021912Z","iopub.execute_input":"2023-03-22T17:21:42.023452Z","iopub.status.idle":"2023-03-22T17:21:43.072593Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25,15)\nn = nrow(df_comPath)\ncnetplot(comPath, showCategory = n) + \n theme_void(base_size = 22) + plot_annotation(\n title = paste0(\"enrichPathway\\nRNA, common for the \", prot_of_interest, \" protein and its RNA\",\n \"\\nEnriched terms shown: \",n, \" of \", nrow(df_comPath))) & \ntheme(text = element_text(size = 22) ) ","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-22T17:21:43.075387Z","iopub.execute_input":"2023-03-22T17:21:43.077040Z","iopub.status.idle":"2023-03-22T17:21:50.979285Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#

Unique RNA for CD36 protein and its RNA

","metadata":{}},{"cell_type":"markdown","source":"We use RNA lists from the main table above with unique RNA for CD36 proteins and its RNA to see if they form different enrichment terms. This may help us to differentiate protein and RNA regulation.","metadata":{}},{"cell_type":"markdown","source":"**Select the set of genes for enrichment analysis**","metadata":{}},{"cell_type":"code","source":"cat(\"------Positively correlated RNA------\\n\\n\")\n\ngene_set_pos_prot <- pos_dat[pos_dat$Pair == prot_of_interest & pos_dat$RNA_list == \"Unique for protein\", ] %>% dplyr::select(\"RNA\",\"corr.coef\")\ncat(\"Unique RNA for\", prot_of_interest, \":\", nrow(gene_set_pos_prot))\ngene_set_pos_prot$RNA\n\ngene_set_pos_rna <- pos_dat[pos_dat$Pair == prot_of_interest & pos_dat$RNA_list == \"Unique for RNA\", ] %>% dplyr::select(\"RNA\",\"corr.coef\")\ncat(\"Unique RNA for\", prot_of_interest, \" RNA:\", nrow(gene_set_pos_rna))\ngene_set_pos_rna$RNA\n\ncat(\"------Negatively correlated RNA------\\n\\n\")\n\ngene_set_neg_prot <- neg_dat[neg_dat$Pair == prot_of_interest & neg_dat$RNA_list == \"Unique for protein\", ] %>% dplyr::select(\"RNA\",\"corr.coef\")\ncat(\"Common RNA for\", prot_of_interest, \"and its RNA:\", nrow(gene_set_neg_prot))\ngene_set_neg_prot$RNA\n\ngene_set_neg_rna <- neg_dat[neg_dat$Pair == prot_of_interest & neg_dat$RNA_list == \"Unique for RNA\", ] %>% dplyr::select(\"RNA\",\"corr.coef\")\ncat(\"Unique RNA for\", prot_of_interest, \" RNA:\", nrow(gene_set_neg_rna))\ngene_set_neg_rna$RNA","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-22T17:21:50.982124Z","iopub.execute_input":"2023-03-22T17:21:50.983734Z","iopub.status.idle":"2023-03-22T17:21:51.066496Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Get Entrez IDs of genes**","metadata":{}},{"cell_type":"code","source":"#-------------positive-------------\n#get Entrez IDs of genes\ngene_an <- AnnotationDbi::select(org.Hs.eg.db, keys = gene_set_pos_prot$RNA,\n columns = c('ENTREZID'), keytype = 'SYMBOL')\n\ncat(\"Subset CD\", prot_of_interest, \"protein: For\", sum(is.na(gene_an$ENTREZID)), \"genes the Entrez ID was not found\")\n\n#check for duplicastes\n#gene_an[gene_an$SYMBOL == gene_an$SYMBOL[duplicated(gene_an$SYMBOL)], ]\n#colnames(mat_RNA)[grepl(gene_an$SYMBOL[duplicated(gene_an$SYMBOL)], colnames(mat_RNA))]\n\n#join to the rest data\ngene_set_pos_prot <- gene_set_pos_prot %>% left_join(gene_an, by = c(\"RNA\" = \"SYMBOL\"))\ngene_set_pos_prot[is.na(gene_set_pos_prot$ENTREZID), ]\n\n\n#get Entrez IDs of genes\ngene_an <- AnnotationDbi::select(org.Hs.eg.db, keys = gene_set_pos_rna$RNA,\n columns = c('ENTREZID'), keytype = 'SYMBOL')\n\ncat(\"Subset CD\", prot_of_interest, \"protein: For\", sum(is.na(gene_an$ENTREZID)), \"genes the Entrez ID was not found\")\n\n#check for duplicastes\n#gene_an[gene_an$SYMBOL == gene_an$SYMBOL[duplicated(gene_an$SYMBOL)], ]\n#colnames(mat_RNA)[grepl(gene_an$SYMBOL[duplicated(gene_an$SYMBOL)], colnames(mat_RNA))]\n\n#join to the rest data\ngene_set_pos_rna <- gene_set_pos_rna %>% left_join(gene_an, by = c(\"RNA\" = \"SYMBOL\"))\ngene_set_pos_rna[is.na(gene_set_pos_rna$ENTREZID), ]\n\n#-------------negative-------------\n#get Entrez IDs of genes\ngene_an <- AnnotationDbi::select(org.Hs.eg.db, keys = gene_set_neg_prot$RNA,\n columns = c('ENTREZID'), keytype = 'SYMBOL')\n\ncat(\"Subset CD\", prot_of_interest, \"protein: For\", sum(is.na(gene_an$ENTREZID)), \"genes the Entrez ID was not found\")\n\n#check for duplicastes\n#gene_an[gene_an$SYMBOL == gene_an$SYMBOL[duplicated(gene_an$SYMBOL)], ]\n#colnames(mat_RNA)[grepl(gene_an$SYMBOL[duplicated(gene_an$SYMBOL)], colnames(mat_RNA))]\n\n#join to the rest data\ngene_set_neg_prot <- gene_set_neg_prot %>% left_join(gene_an, by = c(\"RNA\" = \"SYMBOL\"))\ngene_set_neg_prot[is.na(gene_set_neg_prot$ENTREZID), ]\n\n\n#get Entrez IDs of genes\ngene_an <- AnnotationDbi::select(org.Hs.eg.db, keys = gene_set_neg_rna$RNA,\n columns = c('ENTREZID'), keytype = 'SYMBOL')\n\ncat(\"Subset CD\", prot_of_interest, \"protein: For\", sum(is.na(gene_an$ENTREZID)), \"genes the Entrez ID was not found\")\n\n#check for duplicastes\n#gene_an[gene_an$SYMBOL == gene_an$SYMBOL[duplicated(gene_an$SYMBOL)], ]\n#colnames(mat_RNA)[grepl(gene_an$SYMBOL[duplicated(gene_an$SYMBOL)], colnames(mat_RNA))]\n\n#join to the rest data\ngene_set_neg_rna <- gene_set_neg_rna %>% left_join(gene_an, by = c(\"RNA\" = \"SYMBOL\"))\ngene_set_neg_rna[is.na(gene_set_neg_rna$ENTREZID), ]","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-22T17:21:51.069163Z","iopub.execute_input":"2023-03-22T17:21:51.070601Z","iopub.status.idle":"2023-03-22T17:21:51.810110Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Add missing Gene IDs manually or omit**","metadata":{}},{"cell_type":"code","source":"gene_set_pos_prot <- na.omit(gene_set_pos_prot)\ngene_set_pos_rna <- na.omit(gene_set_pos_rna)\ngene_set_neg_prot <- na.omit(gene_set_neg_prot)\ngene_set_neg_rna <- na.omit(gene_set_neg_rna)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-22T17:21:51.812805Z","iopub.execute_input":"2023-03-22T17:21:51.814307Z","iopub.status.idle":"2023-03-22T17:21:51.834427Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Make a list of genes to compare**","metadata":{}},{"cell_type":"code","source":"uniq <- list(\"CD36_prot_pos\" = gene_set_pos_prot$ENTREZID,\n \"CD36_prot_neg\" = gene_set_neg_prot$ENTREZID,\n \"CD36_RNA_pos\" = gene_set_pos_rna$ENTREZID,\n \"CD36_RNA_neg\" = gene_set_neg_rna$ENTREZID)\nstr(uniq)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-22T17:21:51.836908Z","iopub.execute_input":"2023-03-22T17:21:51.838345Z","iopub.status.idle":"2023-03-22T17:21:51.859315Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#

- GO enrichment analysis (biological processes)

","metadata":{}},{"cell_type":"code","source":"uniqGO <- compareCluster(geneCluster = uniq, OrgDb = org.Hs.eg.db, fun = enrichGO, ont = \"BP\")\nuniqGO <- setReadable(uniqGO, OrgDb = org.Hs.eg.db, keyType=\"ENTREZID\")\n\ndf_uniqGO <- uniqGO %>% as.data.frame\nrownames(df_uniqGO) <- NULL\n\ncat(\"GO enrichment analysis:\\n\")\ncat(\"\\nCD36_prot_pos:\", nrow(df_uniqGO[df_uniqGO$Cluster == \"CD36_prot_pos\",]), \"enriched terms found\",\n \"\\nCD36_prot_neg:\", nrow(df_uniqGO[df_uniqGO$Cluster == \"CD36_prot_neg\",]), \"enriched terms found\",\n \"\\nCD36_RNA_pos:\", nrow(df_uniqGO[df_uniqGO$Cluster == \"CD36_RNA_pos\",]), \"enriched terms found\",\n \"\\nCD36_RNA_neg:\", nrow(df_uniqGO[df_uniqGO$Cluster == \"CD36_RNA_neg\",]), \"enriched terms found\",\n \"\\nTotal:\", nrow(df_uniqGO))\n\ndf_uniqGO$geneID <- str_wrap(df_uniqGO$geneID, width=30)\ndf_uniqGO %>% head(10) %>% select(1,3,4,5,7,9)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-22T17:21:51.861846Z","iopub.execute_input":"2023-03-22T17:21:51.863297Z","iopub.status.idle":"2023-03-22T17:22:40.605520Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"uniqGO <- simplify(uniqGO, cutoff = 0.7, by = \"p.adjust\", select_fun = min)\ncat(\"The clusterProfiler simplify method reduced GO terms to:\")\nuniqGO %>% as.data.frame %>% nrow","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-22T17:22:40.607969Z","iopub.execute_input":"2023-03-22T17:22:40.609394Z","iopub.status.idle":"2023-03-22T17:22:40.716331Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Visualize enriched GO terms**","metadata":{}},{"cell_type":"code","source":"fig(25,12)\nn = uniqGO %>% as.data.frame %>% nrow\ndotplot(uniqGO, showCategory = n) + \n scale_y_discrete(labels=function(x) str_wrap(x, width=100)) + \n theme_bw(base_size = 18) + plot_annotation(\n title = paste0(\"enrichGO\\nRNA, unique for the \", prot_of_interest, \" protein and its RNA\",\n \"\\nEnriched terms shown (per cluster): \",n, \" of \", uniqGO %>% as.data.frame %>% nrow)) & \ntheme(text = element_text(size = 20) ) ","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-22T17:22:40.718845Z","iopub.execute_input":"2023-03-22T17:22:40.720311Z","iopub.status.idle":"2023-03-22T17:22:41.877708Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25,15)\nn = uniqGO %>% as.data.frame %>% nrow\ncnetplot(uniqGO, showCategory = n)+ \n scale_color_gradient2(name='associated data', low='darkgreen', high='firebrick') + \n theme_void(base_size = 22) + plot_annotation(\n title = paste0(\"enrichGO\\nRNA, unique for the \", prot_of_interest, \" protein and its RNA\",\n \"\\nEnriched terms shown: \",n, \" of \", uniqGO %>% as.data.frame %>% nrow)) & \ntheme(text = element_text(size = 22) ) ","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-22T17:22:41.880346Z","iopub.execute_input":"2023-03-22T17:22:41.881936Z","iopub.status.idle":"2023-03-22T17:22:47.736821Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#

- Reactome pathway over-representation analysis

","metadata":{}},{"cell_type":"code","source":"uniqPath <- compareCluster(geneCluster = uniq, fun = enrichPathway, readable = TRUE)\n\ndf_uniqPath <- uniqPath %>% as.data.frame\nrownames(df_uniqPath) <- NULL\n\ncat(\"GO enrichment analysis:\\n\")\ncat(\"\\nCD36_prot_pos:\", nrow(df_uniqPath[df_uniqPath$Cluster == \"CD36_prot_pos\",]), \"enriched terms found\",\n \"\\nCD36_prot_neg:\", nrow(df_uniqPath[df_uniqPath$Cluster == \"CD36_prot_neg\",]), \"enriched terms found\",\n \"\\nCD36_RNA_pos:\", nrow(df_uniqPath[df_uniqPath$Cluster == \"CD36_RNA_pos\",]), \"enriched terms found\",\n \"\\nCD36_RNA_neg:\", nrow(df_uniqPath[df_uniqPath$Cluster == \"CD36_RNA_neg\",]), \"enriched terms found\",\n \"\\nTotal:\", nrow(df_uniqPath)\n)\n\ndf_uniqPath$geneID <- str_wrap(df_uniqPath$geneID, width=30)\ndf_uniqPath %>% head(10) %>% select(1,3,4,5,7,9)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-22T17:22:47.739354Z","iopub.execute_input":"2023-03-22T17:22:47.740818Z","iopub.status.idle":"2023-03-22T17:23:21.393241Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25,8)\nn = nrow(df_uniqPath)\ndotplot(uniqPath, showCategory=n) + \n scale_y_discrete(labels=function(x) str_wrap(x, width=70)) + \n theme_bw(base_size = 18) + plot_annotation(\n title = paste0(\"enrichPathway\\nRNA, common for the \", prot_of_interest, \" protein and its RNA\",\n \"\\nEnriched terms shown: \",n, \" of \", nrow(df_uniqPath))) & \ntheme(text = element_text(size = 22) ) ","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-22T17:23:21.396217Z","iopub.execute_input":"2023-03-22T17:23:21.397886Z","iopub.status.idle":"2023-03-22T17:23:22.017905Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25,10)\nn = nrow(df_uniqPath)\ncnetplot(uniqPath, showCategory = n) + \n theme_void(base_size = 22) + plot_annotation(\n title = paste0(\"enrichPathway\\nRNA, common for the \", prot_of_interest, \" protein and its RNA\",\n \"\\nEnriched terms shown: \",n, \" of \", nrow(df_uniqPath))) & \ntheme(text = element_text(size = 22) ) ","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-22T17:23:22.020394Z","iopub.execute_input":"2023-03-22T17:23:22.021809Z","iopub.status.idle":"2023-03-22T17:23:23.404832Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#

- Terms related to translation, transcription, splicing

","metadata":{}},{"cell_type":"code","source":"dat <- rbind(df_comGO,\n df_comPath,\n df_uniqGO,\n df_uniqPath)\ndat <- dat[grepl(\"splicing|trancrition|translation\", dat$Description, ignore.case = TRUE), ]\ndat","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-22T17:23:23.407524Z","iopub.execute_input":"2023-03-22T17:23:23.409107Z","iopub.status.idle":"2023-03-22T17:23:23.450429Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The corresponding information is directly extracted from the \"GO\" library. The result depends on the currently set ontology (\"BP\",\"MF\",\"CC\"), i.e. only GO terms within the actual ontology are considered. The shown GO information refers to the actually installed GO library.","metadata":{}},{"cell_type":"code","source":"genes <- paste0(dat[grepl(\"trancrition|translation\", dat$Description, ignore.case = TRUE), ]$geneID)\ngenes <- gsub(\"\\\\\\n\", \"\", genes)\ngenes <- str_split(genes, \"/\")\ngenes <- genes %>% unlist %>% unique\ncat(\"Genes:\\n\", \"'\", paste(genes, collapse = \"', '\"),\"'\", sep = \"\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-22T17:23:23.453053Z","iopub.execute_input":"2023-03-22T17:23:23.454526Z","iopub.status.idle":"2023-03-22T17:23:23.478184Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#getGOInfo(\"92906\") #HNRNPLL\n# genes <- c('SNRPD3','POLR2F','SNRPG','YBX1',\n# 'SF3B5','POLR2L','SNU13','SNRPB','SRSF7',\n# 'SRSF2','SNRPD1','UBL5','HSPA8','C1QBP')\n# genes <- AnnotationDbi::select(org.Hs.eg.db, keys = genes,\n# columns = c('ENTREZID'), keytype = 'SYMBOL') #SYMBOL\n# getGOInfo(genes$ENTREZID)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-22T17:23:23.480602Z","iopub.execute_input":"2023-03-22T17:23:23.482026Z","iopub.status.idle":"2023-03-22T17:23:23.493922Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-03-21T20:39:29.966585Z","iopub.execute_input":"2023-03-21T20:39:29.968578Z","iopub.status.idle":"2023-03-21T20:39:29.982868Z"},"trusted":true},"execution_count":null,"outputs":[]}]}