GSVA calculation for NanoString data

gene_x 0 like s 880 view s

Tags: R, sequencing, pipeline

  1. library("rmarkdown")
  2. library("tidyverse")
  3. library(rmarkdown)
  4. library("GeomxTools")
  5. library("GeoMxWorkflows")
  6. library("NanoStringNCTools")
  7. setwd("/home/jhuang/DATA/Data_Susanne_Carotis_spatialRNA_PUBLISHING/run_2023")
  8. render("run.Rmd", c("html_document"))
  12. # ------------------------ Violin Presentation 1 --------------------
  13. #identical(exprs(target_m666Data), assay(target_m666Data, "exprs")): contains the gene expression values that have been both normalized and log-transformed, which are often used for downstream analysis and visualization in many genomic studies.
  15. #assay(target_m666Data, "log_q"): These are likely log2-transformed values of the quantile-normalized data. Quantile normalization ensures that the distributions of gene expression values across samples are the same, and taking the log2 of these normalized values is a common step in the analysis of high-throughput gene expression data. It's a way to make the data more suitable for downstream statistical analysis and visualization.
  17. #assay(target_m666Data, "q_norm"): This typically represents the quantile-normalized gene expression values. Quantile normalization adjusts the expression values in each sample so that they have the same distribution. This normalization method helps remove systematic biases and makes the data more comparable across samples.
  19. library(readxl)
  21. # Path to the Excel file
  22. file_path <- "Signatures.xls"
  24. #example of a signature:
  25. #geneSymbol geneEntrezID    ENSEMBL GeneSet
  26. #CD160  11126   ENSG00000117281 Anergic or act. T cells
  27. #CD244  51744   ENSG00000122223 Anergic or act. T cells
  28. #CTLA4  1493    ENSG00000163599 Anergic or act. T cells
  29. #HAVCR2 84868   ENSG00000135077 Anergic or act. T cells
  30. #ICOS   29851   ENSG00000163600 Anergic or act. T cells
  31. #KLRG1  10219   ENSG00000139187 Anergic or act. T cells
  32. #LAG3   3902    ENSG00000089692 Anergic or act. T cells
  33. #PDCD1  5133    ENSG00000188389 Anergic or act. T cells
  34. #PDCD1  5133    ENSG00000276977 Anergic or act. T cells
  36. # Get the names of the sheets
  37. sheet_names <- excel_sheets(file_path)
  39. # Initialize an empty list to hold gene sets
  40. geneSets <- list()
  42. # Loop over each sheet, extract the ENSEMBL IDs, and add to the list
  43. for (sheet in sheet_names) {
  44.   # Read the sheet
  45.   data <- read_excel(file_path, sheet = sheet)
  47.   # Process the GeneSet names (replacing spaces with underscores, for example)
  48.   gene_set_name <- gsub(" ", "_", unique(data$GeneSet)[1])
  50.   # Add ENSEMBL IDs to the list
  51.   geneSets[[gene_set_name]] <- unique(as.character(data$geneSymbol))
  52. }
  54. # Print the result to check
  55. print(geneSets)
  57. library(gridExtra)
  58. library(ggplot2)
  59. library(GSVA)
  61. # 0. for the following calculation, the GSVA input requires a gene expression matrix 'exprs' where rows are genes and columns are samples.
  62. exprs <- exprs(target_m666Data)
  63. #exprs <- exprs(filtered_or_neg_target_m666Data)
  65. # 1. Compute GSVA scores:
  66. gsva_scores <- gsva(exprs, geneSets, method="gsva")
  68. # 2. Convert to data.frame for ggplot:
  69. gsva_df <- as.data.frame(t(gsva_scores))
  71. # 3. Add conditions to gsva_df:
  72. identical(rownames(pData(target_m666Data)), rownames(gsva_df))
  73. gsva_df$Condition <- pData(target_m666Data)$Grp
  74. #identical(rownames(gsva_df_filtered), rownames(pData(target_m666Data)) )
  75. gsva_df$SampleID <- pData(target_m666Data)$SampleID
  77. # 4. Filter the gsva_df to retain only the desired conditions:
  78. #group 1 vs. group 3 in the nanostring data
  79. gsva_df_filtered <- gsva_df[gsva_df$Condition %in% c("1", "3"), ]
  81. # 5. Define a function to plot violin plots:
  82. # Update the condition levels in gsva_df_filtered to ensure the desired order on x-axis:
  83. gsva_df_filtered$Condition <- gsub("1", "Group1", gsva_df_filtered$Condition)
  84. gsva_df_filtered$Condition <- gsub("3", "Group3", gsva_df_filtered$Condition)
  85. gsva_df_filtered$Condition <- factor(gsva_df_filtered$Condition, levels = c("Group1", "Group3"))
  86. plot_violin <- function(data, gene_name) {
  87.   # Calculate the t-test p-value for the two conditions
  88.   condition1_data <- data[data$Condition == "Group1", gene_name]
  89.   condition2_data <- data[data$Condition == "Group3", gene_name]
  90.   p_value <- t.test(condition1_data, condition2_data)$p.value
  92.   # Convert p-value to annotation
  93.   p_annotation <- ifelse(p_value < 0.01, "**", ifelse(p_value < 0.05, "*", ""))
  94.   rounded_p_value <- paste0("p = ", round(p_value, 2))
  96.   plot_title <- gsub("_", " ", gene_name)
  97.   p <- ggplot(data, aes(x=Condition, y=!!sym(gene_name))) +
  98.     geom_violin(linewidth=1.2) + 
  99.     labs(title=plot_title, y="GSVA Score") +
  100.     ylim(-1, 1) +
  101.     theme_light() +
  102.     theme(
  103.       axis.title.x = element_text(size=12),
  104.       axis.title.y = element_text(size=12),
  105.       axis.text.x  = element_text(size=10),
  106.       axis.text.y  = element_text(size=10),
  107.       plot.title   = element_text(size=12, hjust=0.5)
  108.     )
  110.   # Add p-value annotation to the plot
  111.   p <- p + annotate("text", x=1.5, y=0.9, label=paste0(p_annotation, " ", rounded_p_value), size=5, hjust=0.5)
  113.   return(p)
  114. }
  116. # 6. Generate the list of plots:
  117. #genes <- colnames(gsva_df_filtered)[!colnames(gsva_df_filtered) %in% "Condition"]
  118. #plots_list <- lapply(genes, function(gene) plot_violin(gsva_df_filtered, gene))
  120. # 6. Generate the list of plots in a predefined order:
  121. desired_order <- c("Platelets","Granulocytes","LDG","pDC","Anti-inflammation",  "Pro-inflam._IL-1","Dendritic_cells","MHC_II","Alt._complement","TNF",  "NLRP3_inflammasome","Unfolded_protein","B_cells","Monocyte_cell_surface","Inflammasome",  "Monocyte_secreted","IL-1_cytokines","SNOR_low_UP","CD40_activated","Lectin_complement",  "Classical_complement","Cell_cycle","Plasma_cells","IG_chains","Erythrocytes",  "IL-6R_complex","IFN","TCRB","TCRA","Cyt._act._T_cells",  "TCRG","T_cells","CD8T-NK-NKT","Anergic_or_act._T_cells","T_activated",  "NK_cells","TCRD","T_regs","SNOR_low_DOWN","Monocytes",  "Myeloid_cells","Neutrophils")
  122. genes <- colnames(gsva_df_filtered)[!colnames(gsva_df_filtered) %in% "Condition"]
  123. genes <- genes[match(desired_order, genes)]
  124. genes <- genes[!is.na(genes)]
  125. plots_list <- lapply(genes, function(gene) plot_violin(gsva_df_filtered, gene))
  127. # 7. Pad the list of plots:
  128. remaining_plots <- 6 - (length(plots_list) %% 6)
  129. if (remaining_plots != 6) {
  130.   plots_list <- c(plots_list, rep(list(NULL), remaining_plots))
  131. }
  133. # 7. Pad the list of plots:
  134. #plots_needed_for_full_grid <- ceiling(length(plots_list) / 6) * 6
  135. #remaining_plots <- plots_needed_for_full_grid - length(plots_list)
  136. #plots_list <- c(plots_list, rep(list(NULL), remaining_plots))
  138. # 8. Create the plots and arrange them in a grid:
  139. library(gridExtra)
  140. plots_matrix <- matrix(plots_list, ncol=6, byrow=T)
  141. #do.call("grid.arrange", c(plots_matrix, list(ncol=6)))
  143. # 9. Save the plots to a PNG:
  144. png("All_Violin_Plots_NanoString_3_vs_1.png", width=1000, height=1000)
  145. do.call("grid.arrange", c(plots_matrix, list(ncol=6)))
  146. dev.off()
  149. # ------------------------ Violin Presentation 2 --------------------
  151. # 2. Convert to data.frame for ggplot:
  152. gsva_df <- as.data.frame(t(gsva_scores))
  154. # 3. Add conditions to gsva_df:
  155. identical(rownames(pData(target_m666Data)), rownames(gsva_df))
  156. gsva_df$Condition <- pData(target_m666Data)$Group
  157. #identical(rownames(gsva_df_filtered), rownames(pData(target_m666Data)) )
  158. gsva_df$SampleID <- pData(target_m666Data)$SampleID
  160. # 4. Filter the gsva_df to retain only the desired conditions:
  161. #group 1 vs. group 3 in the nanostring data
  162. gsva_df_filtered <- gsva_df[gsva_df$Condition %in% c("1a", "3"), ]
  164. # 5. Define a function to plot violin plots:
  165. # Update the condition levels in gsva_df_filtered to ensure the desired order on x-axis:
  166. gsva_df_filtered$Condition <- gsub("1a", "Group1a", gsva_df_filtered$Condition)
  167. gsva_df_filtered$Condition <- gsub("3", "Group3", gsva_df_filtered$Condition)
  168. gsva_df_filtered$Condition <- factor(gsva_df_filtered$Condition, levels = c("Group1a", "Group3"))
  169. plot_violin <- function(data, gene_name) {
  170.   # Calculate the t-test p-value for the two conditions
  171.   condition1_data <- data[data$Condition == "Group1a", gene_name]
  172.   condition2_data <- data[data$Condition == "Group3", gene_name]
  173.   p_value <- t.test(condition1_data, condition2_data)$p.value
  175.   # Convert p-value to annotation
  176.   p_annotation <- ifelse(p_value < 0.01, "**", ifelse(p_value < 0.05, "*", ""))
  177.   rounded_p_value <- paste0("p = ", round(p_value, 2))
  179.   plot_title <- gsub("_", " ", gene_name)
  180.   p <- ggplot(data, aes(x=Condition, y=!!sym(gene_name))) +
  181.     geom_violin(linewidth=1.2) + 
  182.     labs(title=plot_title, y="GSVA Score") +
  183.     ylim(-1, 1) +
  184.     theme_light() +
  185.     theme(
  186.       axis.title.x = element_text(size=12),
  187.       axis.title.y = element_text(size=12),
  188.       axis.text.x  = element_text(size=10),
  189.       axis.text.y  = element_text(size=10),
  190.       plot.title   = element_text(size=12, hjust=0.5)
  191.     )
  193.   # Add p-value annotation to the plot
  194.   p <- p + annotate("text", x=1.5, y=0.9, label=paste0(p_annotation, " ", rounded_p_value), size=5, hjust=0.5)
  196.   return(p)
  197. }
  199. # 6. Generate the list of plots:
  200. #genes <- colnames(gsva_df_filtered)[!colnames(gsva_df_filtered) %in% "Condition"]
  201. #plots_list <- lapply(genes, function(gene) plot_violin(gsva_df_filtered, gene))
  203. # 6. Generate the list of plots in a predefined order:
  204. desired_order <- c("Platelets","Granulocytes","LDG","pDC","Anti-inflammation",  "Pro-inflam._IL-1","Dendritic_cells","MHC_II","Alt._complement","TNF",  "NLRP3_inflammasome","Unfolded_protein","B_cells","Monocyte_cell_surface","Inflammasome",  "Monocyte_secreted","IL-1_cytokines","SNOR_low_UP","CD40_activated","Lectin_complement",  "Classical_complement","Cell_cycle","Plasma_cells","IG_chains","Erythrocytes",  "IL-6R_complex","IFN","TCRB","TCRA","Cyt._act._T_cells",  "TCRG","T_cells","CD8T-NK-NKT","Anergic_or_act._T_cells","T_activated",  "NK_cells","TCRD","T_regs","SNOR_low_DOWN","Monocytes",  "Myeloid_cells","Neutrophils")
  205. genes <- colnames(gsva_df_filtered)[!colnames(gsva_df_filtered) %in% "Condition"]
  206. genes <- genes[match(desired_order, genes)]
  207. genes <- genes[!is.na(genes)]
  208. plots_list <- lapply(genes, function(gene) plot_violin(gsva_df_filtered, gene))
  210. # 7. Pad the list of plots:
  211. remaining_plots <- 6 - (length(plots_list) %% 6)
  212. if (remaining_plots != 6) {
  213.   plots_list <- c(plots_list, rep(list(NULL), remaining_plots))
  214. }
  216. # 7. Pad the list of plots:
  217. #plots_needed_for_full_grid <- ceiling(length(plots_list) / 6) * 6
  218. #remaining_plots <- plots_needed_for_full_grid - length(plots_list)
  219. #plots_list <- c(plots_list, rep(list(NULL), remaining_plots))
  221. # 8. Create the plots and arrange them in a grid:
  222. library(gridExtra)
  223. plots_matrix <- matrix(plots_list, ncol=6, byrow=T)
  224. #do.call("grid.arrange", c(plots_matrix, list(ncol=6)))
  226. # 9. Save the plots to a PNG:
  227. png("All_Violin_Plots_NanoString_3_vs_1a.png", width=1000, height=1000)
  228. do.call("grid.arrange", c(plots_matrix, list(ncol=6)))
  229. dev.off()
  234. # ------------------------ Heatmap Presentation --------------------
  236. # 3. Add conditions to gsva_df:
  237. identical(rownames(pData(target_m666Data)), rownames(gsva_df))
  238. gsva_df$Condition <- pData(target_m666Data)$Grp
  239. gsva_df$SampleID <- pData(target_m666Data)$SampleID
  241. # 4. Filter the gsva_df to retain only the desired conditions:
  242. gsva_df_filtered <- gsva_df[gsva_df$Condition %in% c("1", "3"), ]
  244. # 5. Define a function to plot violin plots:
  245. # Update the condition levels in gsva_df_filtered to ensure the desired order on x-axis:
  246. gsva_df_filtered$Condition <- gsub("1", "Group1", gsva_df_filtered$Condition)
  247. gsva_df_filtered$Condition <- gsub("3", "Group3", gsva_df_filtered$Condition)
  248. gsva_df_filtered$Condition <- factor(gsva_df_filtered$Condition, levels = c("Group1", "Group3"))
  251. # Filter the data for "Group1" samples
  252. group1_data <- gsva_df_filtered[gsva_df_filtered$Condition == "Group1", genes]
  254. # Set the median values for each column in the data.frame
  255. # Calculate median values for each column
  256. group1_medians <- apply(group1_data, 2, median)
  260. # 3. Add conditions to gsva_df:
  261. identical(rownames(pData(target_m666Data)), rownames(gsva_df))
  262. gsva_df$Condition <- pData(target_m666Data)$Group
  263. #identical(rownames(gsva_df_filtered), rownames(pData(target_m666Data)) )
  264. gsva_df$SampleID <- pData(target_m666Data)$SampleID
  266. # 4. Filter the gsva_df to retain only the desired conditions:
  267. #group 1 vs. group 3 in the nanostring data
  268. gsva_df_filtered <- gsva_df[gsva_df$Condition %in% c("1a", "1b", "3"), ]
  270. # 5. Define a function to plot violin plots:
  271. # Update the condition levels in gsva_df_filtered to ensure the desired order on x-axis:
  272. gsva_df_filtered$Condition <- gsub("1a", "Group1a", gsva_df_filtered$Condition)
  273. gsva_df_filtered$Condition <- gsub("1b", "Group1b", gsva_df_filtered$Condition)
  274. gsva_df_filtered$Condition <- gsub("3", "Group3", gsva_df_filtered$Condition)
  275. gsva_df_filtered$Condition <- factor(gsva_df_filtered$Condition, levels = c("Group1a", "Group1b", "Group3"))
  277. # Filter the data for "Group1a" samples
  278. group1a_data <- gsva_df_filtered[gsva_df_filtered$Condition == "Group1a", genes]
  279. group1a_medians <- apply(group1a_data, 2, median)
  281. # Filter the data for "Group1b" samples
  282. group1b_data <- gsva_df_filtered[gsva_df_filtered$Condition == "Group1b", genes]
  283. group1b_medians <- apply(group1b_data, 2, median)
  285. # Filter the data for "Group3" samples
  286. group3_data <- gsva_df_filtered[gsva_df_filtered$Condition == "Group3", genes]
  287. group3_medians <- apply(group3_data, 2, median)
  289. # Create a new data frame with the middle values for both groups
  290. heatmap_data_nanostring <- data.frame(Group1 = group1_medians, Group1a = group1a_medians, Group1b = group1b_medians, Group3 = group3_medians)
  292. # Merge two heatmap_data_* to a matrix
  293. heatmap_data <- merge(heatmap_data_nanostring, heatmap_data_rnaseq, by="row.names", all=TRUE)
  294. rownames(heatmap_data) <- heatmap_data$Row.names
  295. heatmap_data$Row.names <- NULL
  297. # Transpose the data matrix to switch rows and columns
  298. #heatmap_data <- t(heatmap_data)
  300. # Convert the data frame to a numeric matrix
  301. heatmap_matrix <- as.matrix(heatmap_data)
  303. # Transpose the data matrix to switch rows and columns
  304. #heatmap_data <- t(heatmap_data)
  306. # Convert the data frame to a numeric matrix
  307. heatmap_matrix <- as.matrix(heatmap_data)
  309. # Specify heatmap colors
  310. color_scheme <- colorRampPalette(c("blue", "white", "red"))(50)
  312. #TODO: adjust the color, so that it is not so dark and light!!!!, add RNAseq data for mock and infection!!
  313. # Create the heatmap without hierarchical clustering
  314. library(gplots)
  315. png("Heatmap.png", width=1000, height=1000)
  316. heatmap.2(
  317.   heatmap_matrix,
  318.   Colv = FALSE,  # Turn off column dendrogram
  319.   Rowv = FALSE,  # Turn off row dendrogram
  320.   trace = "none",  # Turn off row and column dendrograms
  321.   scale = "row",   # Scale the rows
  322.   col = color_scheme,
  323.   main = "GSVA Heatmap",  # Title of the heatmap
  324.   key = TRUE,            # Show color key
  325.   key.title = "GSVA Score",  # Color key title
  326.   key.xlab = NULL,       # Remove x-axis label from color key
  327.   density.info = "none",  # Remove density plot
  328.   cexRow = 1.2,  # Increase font size of row names
  329.   cexCol = 1.2,  # Increase font size of column names
  330.   srtCol = 40, keysize = 0.72, margins = c(5, 14)
  331. )
  332. dev.off()
  334. # Write the table to Excel-file
  335. library(writexl)
  336. # Specify the file path where you want to save the Excel file
  337. excel_file <- "heatmap_data.xlsx"
  339. # Save the data frame to the Excel file
  340. write_xlsx(heatmap_data, excel_file)
  342. # Verify that the Excel file has been saved
  343. if (file.exists(excel_file)) {
  344.   cat("Data frame has been saved to", excel_file, "\n")
  345. } else {
  346.   cat("Error: Failed to save the data frame to", excel_file, "\n")
  347. }

like unlike

点赞本文的读者

还没有人对此文章表态

本文有评论

没有评论

看文章，发评论，不要沉默

最受欢迎文章

  1. Why Do Significant Gene Lists Change After Adding Additional Conditions in Differential Gene Expression Analysis?
  2. Motif Discovery in Biological Sequences: A Comparison of MEME and HOMER
  3. Calling peaks using findPeaks of HOMER
  4. Updating Human Gene Identifiers using Ensembl BioMart: A Step-by-Step Guide
  5. Setup conda environments
  6. Should the inputs for GSVA be normalized or raw?
  7. pheatmap vs heatmap.2
  8. File format for single channel analysis of Agilent microarray data with Limma?
  9. Kraken2 Installation and Usage Guide
  10. 洛那法尼(lonafarnib)：从抗癌症到抗病毒的多功能药物

最新文章

最多评论文章

推荐相似文章

piRNA 2024 Ute from raw counts

Tn-seq analysis pipeline (improved)

Step-by-Step Analysis of Transposon Insertion Sites

QIIME + Phyloseq + MicrobiotaProcess (v1)

© 2023 XGenes.com Impressum