--- title: "(20231027) Circadian_Rhythms_DE" output: pdf_document --- ```{r} # load libraries library(DESeq2) library(apeglm) library(EnhancedVolcano) library(org.Mm.eg.db) library(magrittr) library(GenomicFeatures) library(tidyverse) # BiocManager::install("DESeq2") # BiocManager::install("apeglm") # BiocManager::install("EnhancedVolcano") ``` Load the count matrix, and perform some data wrangling ```{r Prepare count matrix} counts <- readRDS("~/Circadian_Rhythms/count_matrix_20230925.rds") # Check dimensions of the matrix and use head() dim(counts) colnames(counts) # Taking a glance of the matrix counts[1:10, 1:6] # we should remove the first four rows and make the gene column as the row names counts <- data.frame(counts) rownames(counts) <- counts$gene counts %<>% dplyr::select(-gene) counts <- counts[5:nrow(counts),] # Let's add the gene symbol at this point, and filter out those that don't have gene symbols from the counts matrix ensembl_ids <- rownames(counts) # This function checks for the gene symbols that match the particular ENSEMBL ID. gene_symbols <- mapIds(org.Mm.eg.db, keys = ensembl_ids, column = "SYMBOL", keytype = "ENSEMBL") # how many gene symbols? length(gene_symbols) gene_symbols %<>% as_tibble(rownames = "ENSEMBL") %>% dplyr::rename(gene = value) gene_symbols %>% group_by(gene) %>% dplyr::count() %>% filter(n != 1) # Filter out the genes where it doesn't have a gene symbol gene_symbols %<>% dplyr::filter(!is.na(gene)) # Save the genes that are NAs (invalid since we already filtered them out) # gene_symbols_NA <- gene_symbols %>% dplyr::filter(is.na(gene)) #check for duplicate ENSEMBL IDs gene_symbols %>% group_by(ENSEMBL) %>% dplyr::count() %>% filter(n != 1) #check for duplicate gene symbols gene_symbols %>% group_by(gene) %>% dplyr::count() %>% filter(n != 1) # Filter the counts matrix (these have a registered gene symbol) counts <- counts[gene_symbols$ENSEMBL,] # Check order for gene symbols summary(gene_symbols$ENSEMBL == rownames(counts)) # We'll select out the first gene ENSEMBL ID among the duplicates. counts <- counts[!duplicated(gene_symbols$gene),] ``` ```{r colData Wrangling} # Load the mouse gtf file, containing information on each transcript gtf_file <- "/home/youkim/Genome_indexing/Mus_musculus.GRCm39.110.gtf" txdb <- makeTxDbFromGFF(gtf_file, format = "gtf") # Select out the Y chromosome chromosome_name <- "Y" #Extract the Y chromosome Ensembl IDs genes_on_Y <- genes(txdb, filter = list(tx_chrom = "Y")) ensembl_gene_ids <- genes_on_Y@ranges@NAMES #Check for all Y chromosome genes Y_chr_expression <- counts[(rownames(counts) %in% ensembl_gene_ids),] %>% colSums() %>% data.frame(Y_exp = .) # We should also set the colData for DESeq2 coldata <- data.frame(sample = colnames(counts), Y_chr_expression) #check if the order of the rows are identical rownames(coldata) == rownames(Y_chr_expression) # Add biological sex information in metadata coldata <- coldata %>% mutate(sex = ifelse(coldata$Y_exp > 100, "male", "female")) ``` ```{r} # Now we can change the rownames of the count matrix to gene symbols! counts <- counts %>% data.frame() %>% rownames_to_column() %>% rename(ENSEMBL = rowname) %>% inner_join(gene_symbols) counts <- counts %>% column_to_rownames("gene") %>% dplyr::select(-ENSEMBL) colnames(counts) %>% head() ``` Let's also modify the column names and sample information, before running DESeq2 ```{r} # Let's also add columns addressing their circadian status (time of the day) # Load the .csv file including the timepoint data tpdata <- read_csv("~/Circadian_Rhythms/Sample_Info(Sept.2023).csv", col_names = F) # Modify the second column extracted_part <- str_extract(tpdata$X2, "(?<=_).*?(?=_)") %>% rev() # Attach the timepoint data to our bulk RNA-seq metadata. coldata <- coldata %>% mutate(time_point = extracted_part) # We want to modify our colnames of the count matrix, so that it can be compatible for NiteCap sample_numbers <- str_extract(coldata$sample, "(?<=_)[^_]+(?=_)") coldata <- coldata %>% mutate(sample_numbers = sample_numbers) new_colnames <- paste0(coldata$sex, "_", coldata$time_point, "_", coldata$sample_numbers) # Let's add the new column names colnames(counts) <- new_colnames rownames(coldata) <- new_colnames # Take an overview of the counts and coldata dataframes colnames(counts) %>% grepl("female") head(coldata) male_count <- counts[, !grepl("female", colnames(counts))] male_contrast <- bind_cols(male_count[,grepl("ZT7", colnames(male_count))], male_count[,grepl("ZT19", colnames(male_count))]) female_count <- counts[, grepl("female", colnames(counts))] female_contrast <- bind_cols(female_count[,grepl("ZT7", colnames(female_count))], female_count[,grepl("ZT19", colnames(female_count))]) # Let's also arrange the columns by time frames # tibble_colnames_1 <- tibble(name = colnames(male_count)) %>% # filter(grepl("male_", name)) %>% # mutate(timepoint = as.numeric(gsub(".*_ZT(.*)_.*", "\\1", name))) %>% # arrange(timepoint) # male_count_modified <- male_count[,tibble_colnames_1$name] # # tibble_colnames_2 <- tibble(name = colnames(female_count)) %>% # filter(grepl("female_", name)) %>% # mutate(timepoint = as.numeric(gsub(".*_ZT(.*)_.*", "\\1", name))) %>% # arrange(timepoint) # female_count_modified <- female_count[,tibble_colnames_2$name] # average_per_gene_M <- sapply(seq(1, ncol(male_count_modified), by = 6), function(start_col) { # end_col <- min(start_col + 5, ncol(male_count_modified)) # Determine the end column for the group # avg_exprs <- rowSums(male_count_modified[, start_col:end_col])/6 # return(avg_exprs) # }) %>% data.frame() # # colnames(average_per_gene_M) <- c("M_ZT3", "M_ZT7", "M_ZT11", "M_ZT15", "M_ZT19", "M_ZT23") # # average_per_gene_F <- sapply(seq(1, ncol(female_count_modified), by = 6), function(start_col) { # end_col <- min(start_col + 5, ncol(female_count_modified)) # Determine the end column for the group # avg_exprs <- rowSums(female_count_modified[, start_col:end_col])/6 # return(avg_exprs) # }) %>% data.frame() # # colnames(average_per_gene_F) <- c("F_ZT3", "F_ZT7", "F_ZT11", "F_ZT15", "F_ZT19", "F_ZT23") # counts_revised <- bind_cols(average_per_gene_M, average_per_gene_F) ``` Contrary to what we did for the timepoint-based analysis, we only need two count matrices. ```{r} # Use a for loop to generate the 6 matrices write_csv(male_contrast, "male_contrast.csv") write_csv(female_contrast, "female_contrast_(ZT7vsZT19).csv") ``` Let's revise the coldata ```{r} # Load relevant library library(stringr) # Regular expression pattern to match text between underscores pattern <- "_([^_]+)_" male_coldata_revised <- data.frame(sample = colnames(male_contrast)) %>% mutate(timepoint = str_match(colnames(male_contrast) ,pattern)[,2]) rownames(male_coldata_revised) <- male_coldata_revised$sample # write_csv(male_coldata_revised, "male_coldata_revised.csv") female_coldata_revised <- data.frame(sample = colnames(female_contrast)) %>% mutate(timepoint = str_match(colnames(female_contrast) ,pattern)[,2]) rownames(female_coldata_revised) <- female_coldata_revised$sample # write_csv(female_coldata_revised, "female_coldata_revised.csv") # At this point, we can generate the DESeqDataSetFromMatrix smallestGroupSize <- 6 male_contrast <- male_contrast %>% mutate_at(colnames(male_contrast), as.integer) male_dds_revised <- DESeqDataSetFromMatrix(countData = male_contrast, colData = male_coldata_revised, design= ~ timepoint) female_contrast <- female_contrast %>% mutate_at(colnames(female_contrast), as.integer) female_dds_revised <- DESeqDataSetFromMatrix(countData = female_contrast, colData = female_coldata_revised, design= ~ timepoint) summary(rownames(male_dds_revised@colData) == colnames(male_dds_revised@assays@data$counts)) summary(rownames(female_dds_revised@colData) == colnames(female_dds_revised@assays@data$counts)) # We will filter out genes that have less than 10 reads. keep <- rowSums(counts(male_dds_revised) >= 10) >= smallestGroupSize male_dds <- male_dds_revised[keep,] keep <- rowSums(counts(female_dds_revised) >= 10) >= smallestGroupSize female_dds <- female_dds_revised[keep,] # Check the dimensions of the count matrix dim(male_dds) %>% print() dim(female_dds) %>% print() # At this point, I run DESeq. Compare before and after running DESeq male_dds <- DESeq(male_dds) female_dds <- DESeq(female_dds) # Save the dds R object as a RDS file # saveRDS(male_dds, "male_dds_20240302.rds") saveRDS(female_dds, "female_dds_20240304.rds") ``` I think we can write the prefiltering step here too, where we follow the original vignette! ```{r} not_normalized_counts <- counts(female_dds, normalized = FALSE) normalized_counts <- counts(female_dds, normalized = TRUE) head(normalized_counts) ``` ```{r} # Let's recall our dds object male_dds <- readRDS("male_dds_20240302.rds") female_dds <- readRDS("female_dds_20240304.rds") resultsNames(male_dds) resultsNames(female_dds) male_res <- results(male_dds, name = "timepoint_ZT7_vs_ZT19") female_res <- results(female_dds, name = "timepoint_ZT7_vs_ZT19") pdf(file = "Volcanoplot_20240304.pdf", height = 8, width = 8) EnhancedVolcano(female_res, lab = rownames(female_res), x = "log2FoldChange", pCutoff = 0.05, FCcutoff = 1.0, y= "pvalue", drawConnectors = TRUE, widthConnectors = 0.75, max.overlaps = 100, title = "Female ZT7 vs ZT19 Volcano Plot") EnhancedVolcano(male_res, lab = rownames(male_res), x = "log2FoldChange", pCutoff = 0.05, FCcutoff = 1.0, y= "padj", drawConnectors = TRUE, widthConnectors = 0.75, max.overlaps = 100, title = "Male ZT7 vs ZT19 Volcano Plot (adjusted p-value)") EnhancedVolcano(female_res, lab = rownames(female_res), x = "log2FoldChange", pCutoff = 0.05, FCcutoff = 1.0, y= "padj", drawConnectors = TRUE, widthConnectors = 0.50, max.overlaps = 100, labSize = 4.0, title = "Female ZT7 vs ZT19 Volcano Plot (adjusted p-value)") EnhancedVolcano(female_res, lab = rownames(female_res), x = "log2FoldChange", pCutoff = 0.05, FCcutoff = 1.5, y= "padj", drawConnectors = TRUE, widthConnectors = 0.75, max.overlaps = 100, title = "Female ZT7 vs ZT19 Volcano Plot (adjusted p-value)") dev.off() male_contrast["Bmal1",] female_contrast["Cry1",] female_res[abs(female_res$log2FoldChange) > 1,] %>% data.frame() %>% dplyr::filter(is.na(.$padj) == FALSE & .$padj < 0.05) %>% dplyr::filter(log2FoldChange < -1) %>% rownames() %>% homologene::mouse2human() %>% pull(humanGene) %>% cat() ``` ```{r} male_samples <- dds@colData[dds@colData$sex == "male",] female_samples <- dds@colData[dds@colData$sex == "female",] male_dds <- dds[, (colnames(dds) %in% rownames(male_samples))] female_dds <- dds[, (colnames(dds) %in% rownames(female_samples))] male_count <- counts(male_dds, normalized = TRUE) male_count <- male_count[, ncol(male_dds@assays@data$counts):1] head(male_count) female_count <- counts(female_dds, normalized = TRUE) female_count <- female_count[, ncol(female_dds@assays@data$counts):1] head(female_count) # Let's also arrange the columns by time frames tibble_colnames_1 <- tibble(name = colnames(male_count)) %>% filter(grepl("male_", name)) %>% mutate(timepoint = as.numeric(gsub(".*_ZT(.*)_.*", "\\1", name))) %>% arrange(timepoint) male_count_modified <- male_count[,tibble_colnames_1$name] tibble_colnames_2 <- tibble(name = colnames(female_count)) %>% filter(grepl("female_", name)) %>% mutate(timepoint = as.numeric(gsub(".*_ZT(.*)_.*", "\\1", name))) %>% arrange(timepoint) female_count_modified <- female_count[,tibble_colnames_2$name] # Let's use head() to see if our columns look correct for males and females head(male_count_modified) head(female_count_modified) # # We save the count matrices as a new .csv file male_count_modified <- male_count_modified %>% data.frame() %>% rownames_to_column() %>% dplyr::rename(gene = rowname) female_count_modified <- female_count_modified %>% data.frame() %>% rownames_to_column() %>% dplyr::rename(gene = rowname) # write.csv(male_count_modified, "DESeq2_normalized_male_counts_20231106.csv", row.names = FALSE) # write.csv(female_count_modified, "DESeq2_normalized_female_counts_20231106.csv", row.names = FALSE) ``` ```{r EnhancedVolcano} # Check for the results from DESeq2 resultsNames(dds) res <- results(dds, name = "sex_male_vs_female") resLFC <- lfcShrink(dds, coef = "sex_male_vs_female", type = "apeglm") # Let's take a quick look at a visualization of the results plotMA(res, ylim = c(-2,2)) # MAplot for LFC results plotMA(resLFC, ylim = c(-2,2)) ``` I think at this point we can implement the EnhancedVolcano package, which is great when we use DESeq2. ```{r} # We will make two versions of the EnhancedVolcano plots # General results res_filtered <- res[!(rownames(res) %in% c("Xist", "Ddx3y", "Kdm5d", "Eif2s3y","Uty")),] pdf(file = "Volcanoplot_20231107.pdf", height = 12, width = 6) EnhancedVolcano(res_filtered, lab = rownames(res_filtered), x = "log2FoldChange", pCutoff = 10e-10, FCcutoff = 1.0, y= "pvalue" ) dev.off() EnhancedVolcano(resLFC, lab = rownames(resLFC), x = "log2FoldChange", y= "pvalue" ) ```