--- title: "(20231027) Circadian_Rhythms_DE" output: pdf_document --- ```{r} # load libraries library(DESeq2) library(apeglm) library(EnhancedVolcano) library(org.Mm.eg.db) # 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) ``` 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))] female_count <- counts[, grepl("female", colnames(counts))] # 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) ``` At this point, we need to generate 6 count matrices, based on different timepoints. ```{r} # Use a for loop to generate the 6 matrices for(i in 1:6){ a <- grepl(paste0("ZT", 4*i-1),colnames(male_count_modified)) df <- bind_cols(male_count_modified[,a],female_count_modified[,a]) saveRDS(df, paste0("ZT", 4*i-1,".csv")) } ``` Let's revise the coldata ```{r} # Let's use the for loop again for(i in 1:6){ a <- readRDS(paste0("ZT",4*i-1,".csv")) coldata_revised <- data.frame(sample = colnames(a)) %>% mutate(sex = substr(colnames(b),1,1)) rownames(coldata_revised) <- coldata_revised$sample saveRDS(coldata_revised, paste0("coldata_ZT",4*i-1,".csv")) } # At this point, we can generate the DESeqDataSetFromMatrix smallestGroupSize <- 6 for(i in 1:6){ a <- readRDS(paste0("ZT",4*i-1,".csv")) a <- a %>% mutate_at(colnames(a), as.integer) coldata_revised <- readRDS(paste0("coldata_ZT",4*i-1,".csv")) dds_revised <- DESeqDataSetFromMatrix(countData = a, colData = coldata_revised, design= ~ sex) summary(rownames(dds_revised@colData) == colnames(dds_revised@assays@data$counts)) %>% print() # We will filter out genes that have less than 10 reads. keep <- rowSums(counts(dds_revised) >= 10) >= smallestGroupSize dds <- dds_revised[keep,] # Check the dimensions of the count matrix dim(dds) %>% print() # At this point, I run DESeq. Compare before and after running DESeq dds <- DESeq(dds) # Save the dds R object as a RDS file saveRDS(dds, paste0("DDS_ZT",4*i-1,"_20240227.rds")) } ``` I think we can write the prefiltering step here too, where we follow the original vignette! ```{r} not_normalized_counts <- counts(dds, normalized = FALSE) normalized_counts <- counts(dds, normalized = TRUE) head(normalized_counts) ``` ```{r} # To remove Y chromosome genes, 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 # This function checks for the gene symbols that match the particular ENSEMBL ID. gene_symbols <- mapIds(org.Mm.eg.db, keys = ensembl_gene_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) dim(gene_symbols) # Let's recall our dds object for(i in 1:6){ dds <- readRDS(paste0("DDS_ZT",4*i-1,"_20240227.rds")) # check results resultsNames(dds) %>% print() dim(dds) %>% print() res <- results(dds, name = "sex_m_vs_f") # Now remove the genes from our results that are located in the Y chromosome. Y_genes <- intersect(gene_symbols$gene, rownames(res)) res <- res[!(rownames(res) %in% Y_genes),] res_filtered <- res[!(rownames(res) %in% c("Xist")),] print("Now generating volcano plot") pdf(file = paste0("Volcanoplot_ZT",4*i-1,"_20240227.pdf"), height = 8, width = 8) print(EnhancedVolcano(res, lab = rownames(res), x = "log2FoldChange", pCutoff = 0.05, FCcutoff = 2.0, y= "pvalue", drawConnectors = TRUE, widthConnectors = 0.75, max.overlaps = 100, title = paste0("ZT", 4*i-1, " Volcano Plot"))) print(EnhancedVolcano(res_filtered, lab = rownames(res_filtered), x = "log2FoldChange", y= "pvalue", pCutoff = 0.05, FCcutoff = 1.0, drawConnectors = TRUE, widthConnectors = 0.75, max.overlaps = 100, title = paste0("ZT", 4*i-1, " Volcano Plot"))) print(EnhancedVolcano(res, lab = rownames(res), x = "log2FoldChange", pCutoff = 0.05, FCcutoff = 1.0, y= "padj", drawConnectors = TRUE, widthConnectors = 0.75, max.overlaps = 100, title = paste0("ZT", 4*i-1, " Volcano Plot"))) print(EnhancedVolcano(res_filtered, lab = rownames(res_filtered), x = "log2FoldChange", y= "padj", pCutoff = 0.05, FCcutoff = 1.0, drawConnectors = TRUE, widthConnectors = 0.75, max.overlaps = 100, title = paste0("ZT", 4*i-1, " Volcano Plot"))) dev.off() } interim <- res_filtered[!is.na(res_filtered$padj < 0.05),] interim[(interim$padj < 0.05),] interim[(interim$pvalue < 0.05),] ``` ```{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" ) ``` We notice that three sex-specific genes are too strong in the volcano plot ```{r} dds ```