(20231027) Circadian_Rhythms_DE
# load libraries
library(DESeq2)## Loading required package: S4Vectors## Loading required package: stats4## Loading required package: BiocGenerics##
## Attaching package: 'BiocGenerics'## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs## The following objects are masked from 'package:base':
##
## anyDuplicated, aperm, append, as.data.frame, basename, cbind,
## colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
## get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
## match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
## Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
## table, tapply, union, unique, unsplit, which.max, which.min##
## Attaching package: 'S4Vectors'## The following object is masked from 'package:utils':
##
## findMatches## The following objects are masked from 'package:base':
##
## expand.grid, I, unname## Loading required package: IRanges## Loading required package: GenomicRanges## Loading required package: GenomeInfoDb## Loading required package: SummarizedExperiment## Loading required package: MatrixGenerics## Loading required package: matrixStats##
## Attaching package: 'MatrixGenerics'## The following objects are masked from 'package:matrixStats':
##
## colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
## colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
## colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
## colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
## colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
## colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
## colWeightedMeans, colWeightedMedians, colWeightedSds,
## colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
## rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
## rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
## rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
## rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
## rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
## rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
## rowWeightedSds, rowWeightedVars## Loading required package: Biobase## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.##
## Attaching package: 'Biobase'## The following object is masked from 'package:MatrixGenerics':
##
## rowMedians## The following objects are masked from 'package:matrixStats':
##
## anyMissing, rowMedianslibrary(apeglm)
library(EnhancedVolcano)## Loading required package: ggplot2## Loading required package: ggrepellibrary(org.Mm.eg.db)## Loading required package: AnnotationDbi## library(magrittr)##
## Attaching package: 'magrittr'## The following object is masked from 'package:GenomicRanges':
##
## subtractlibrary(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.3 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ lubridate 1.9.3 ✔ tibble 3.2.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.0## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ lubridate::%within%() masks IRanges::%within%()
## ✖ dplyr::collapse() masks IRanges::collapse()
## ✖ dplyr::combine() masks Biobase::combine(), BiocGenerics::combine()
## ✖ dplyr::count() masks matrixStats::count()
## ✖ dplyr::desc() masks IRanges::desc()
## ✖ tidyr::expand() masks S4Vectors::expand()
## ✖ tidyr::extract() masks magrittr::extract()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::first() masks S4Vectors::first()
## ✖ dplyr::lag() masks stats::lag()
## ✖ ggplot2::Position() masks BiocGenerics::Position(), base::Position()
## ✖ purrr::reduce() masks GenomicRanges::reduce(), IRanges::reduce()
## ✖ dplyr::rename() masks S4Vectors::rename()
## ✖ lubridate::second() masks S4Vectors::second()
## ✖ lubridate::second<-() masks S4Vectors::second<-()
## ✖ dplyr::select() masks AnnotationDbi::select()
## ✖ purrr::set_names() masks magrittr::set_names()
## ✖ dplyr::slice() masks IRanges::slice()
## ✖ magrittr::subtract() masks GenomicRanges::subtract()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(GenomicFeatures)
# BiocManager::install("DESeq2")
# BiocManager::install("apeglm")
# BiocManager::install("EnhancedVolcano")Load the count matrix, and perform some data wrangling
counts <- readRDS("count_matrix_20230925.rds")
# Check dimensions of the matrix and use head()
dim(counts)## [1] 56945 73colnames(counts)## [1] "gene" "Bolz_072_Plate1_Strip9_H09_72__S132"
## [3] "Bolz_071_Plate1_Strip9_G09_71__S131" "Bolz_070_Plate1_Strip9_F09_70__S130"
## [5] "Bolz_069_Plate1_Strip9_E09_69__S129" "Bolz_068_Plate1_Strip9_D09_68__S128"
## [7] "Bolz_067_Plate1_Strip9_C09_67__S127" "Bolz_066_Plate1_Strip9_B09_66__S126"
## [9] "Bolz_065_Plate1_Strip9_A09_65__S125" "Bolz_064_Plate1_Strip9_H08_64__S124"
## [11] "Bolz_063_Plate1_Strip8_G08_63__S123" "Bolz_062_Plate1_Strip8_F08_62__S122"
## [13] "Bolz_061_Plate1_Strip8_E08_61__S121" "Bolz_060_Plate1_Strip8_D08_60__S120"
## [15] "Bolz_059_Plate1_Strip8_C08_59__S119" "Bolz_058_Plate1_Strip8_B08_58__S118"
## [17] "Bolz_057_Plate1_Strip8_A08_57__S117" "Bolz_056_Plate1_Strip7_H07_56__S116"
## [19] "Bolz_055_Plate1_Strip7_G07_55__S115" "Bolz_054_Plate1_Strip7_F07_54__S114"
## [21] "Bolz_053_Plate1_Strip7_E07_53__S113" "Bolz_052_Plate1_Strip7_D07_52__S112"
## [23] "Bolz_051_Plate1_Strip7_C07_51__S111" "Bolz_050_Plate1_Strip7_B07_50__S110"
## [25] "Bolz_049_Plate1_Strip7_A07_49__S109" "Bolz_048_Plate1_Strip6_H06_48__S108"
## [27] "Bolz_047_Plate1_Strip6_G06_47__S107" "Bolz_046_Plate1_Strip6_F06_46__S106"
## [29] "Bolz_045_Plate1_Strip6_E06_45__S105" "Bolz_044_Plate1_Strip6_D06_44__S104"
## [31] "Bolz_043_Plate1_Strip6_C06_43__S103" "Bolz_042_Plate1_Strip6_B06_42__S102"
## [33] "Bolz_041_Plate1_Strip6_A06_41__S101" "Bolz_040_Plate1_Strip5_H05_40__S100"
## [35] "Bolz_039_Plate1_Strip5_G05_39__S99" "Bolz_038_Plate1_Strip5_F05_38__S98"
## [37] "Bolz_037_Plate1_Strip5_E05_37__S97" "Bolz_036_Plate1_Strip5_D05_36__S96"
## [39] "Bolz_035_Plate1_Strip5_C05_35__S95" "Bolz_034_Plate1_Strip5_B05_34__S94"
## [41] "Bolz_033_Plate1_Strip5_A05_33__S93" "Bolz_032_Plate1_Strip4_H04_32__S92"
## [43] "Bolz_031_Plate1_Strip4_G04_31__S91" "Bolz_030_Plate1_Strip4_F04_30__S90"
## [45] "Bolz_029_Plate1_Strip4_E04_29__S89" "Bolz_028_Plate1_Strip4_D04_28__S88"
## [47] "Bolz_027_Plate1_Strip4_C04_27__S87" "Bolz_026_Plate1_Strip4_B04_26__S86"
## [49] "Bolz_025_Plate1_Strip4_A04_25__S85" "Bolz_024_Plate1_Strip3_H03_24__S84"
## [51] "Bolz_023_Plate1_Strip3_G03_23__S83" "Bolz_022_Plate1_Strip3_F03_22__S82"
## [53] "Bolz_021_Plate1_Strip3_E03_21__S81" "Bolz_020_Plate1_Strip3_D03_20__S80"
## [55] "Bolz_019_Plate1_Strip3_C03_19__S79" "Bolz_018_Plate1_Strip3_B03_18__S78"
## [57] "Bolz_017_Plate1_Strip3_A03_17__S77" "Bolz_016_Plate1_Strip2_H02_16__S76"
## [59] "Bolz_015_Plate1_Strip2_G02_15__S75" "Bolz_014_Plate1_Strip2_F02_14__S74"
## [61] "Bolz_013_Plate1_Strip2_E02_13__S73" "Bolz_012_Plate1_Strip2_D02_12__S72"
## [63] "Bolz_011_Plate1_Strip2_C02_11__S71" "Bolz_010_Plate1_Strip2_B02_10__S70"
## [65] "Bolz_009_Plate1_Strip2_A02_9__S69" "Bolz_008_Plate1_Strip1_H01_8__S68"
## [67] "Bolz_007_Plate1_Strip1_G01_7__S67" "Bolz_006_Plate1_Strip1_F01_6__S66"
## [69] "Bolz_005_Plate1_Strip1_E01_5__S65" "Bolz_004_Plate1_Strip1_D01_4__S64"
## [71] "Bolz_003_Plate1_Strip1_C01_3__S63" "Bolz_002_Plate1_Strip1_B01_2__S62"
## [73] "Bolz_001_Plate1_Strip1_A01_1__S61"# 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")## 'select()' returned 1:many mapping between keys and columns# how many gene symbols?
length(gene_symbols)## [1] 56941gene_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))## Mode TRUE
## logical 33148# We'll select out the first gene ENSEMBL ID among the duplicates.
counts <- counts[!duplicated(gene_symbols$gene),]# 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")## Import genomic features from the file as a GRanges object ... OK
## Prepare the 'metadata' data frame ... OK
## Make the TxDb object ...## Warning in .get_cds_IDX(mcols0$type, mcols0$phase): The "phase" metadata column contains non-NA values for features of type
## stop_codon. This information was ignored.## OK# 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)## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [31] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [46] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [61] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE# Add biological sex information in metadata
coldata <- coldata %>%
mutate(sex = ifelse(coldata$Y_exp > 100, "male", "female"))# 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)## Joining with `by = join_by(ENSEMBL)`counts <- counts %>%
column_to_rownames("gene") %>%
dplyr::select(-ENSEMBL)
colnames(counts)## [1] "Bolz_072_Plate1_Strip9_H09_72__S132" "Bolz_071_Plate1_Strip9_G09_71__S131"
## [3] "Bolz_070_Plate1_Strip9_F09_70__S130" "Bolz_069_Plate1_Strip9_E09_69__S129"
## [5] "Bolz_068_Plate1_Strip9_D09_68__S128" "Bolz_067_Plate1_Strip9_C09_67__S127"
## [7] "Bolz_066_Plate1_Strip9_B09_66__S126" "Bolz_065_Plate1_Strip9_A09_65__S125"
## [9] "Bolz_064_Plate1_Strip9_H08_64__S124" "Bolz_063_Plate1_Strip8_G08_63__S123"
## [11] "Bolz_062_Plate1_Strip8_F08_62__S122" "Bolz_061_Plate1_Strip8_E08_61__S121"
## [13] "Bolz_060_Plate1_Strip8_D08_60__S120" "Bolz_059_Plate1_Strip8_C08_59__S119"
## [15] "Bolz_058_Plate1_Strip8_B08_58__S118" "Bolz_057_Plate1_Strip8_A08_57__S117"
## [17] "Bolz_056_Plate1_Strip7_H07_56__S116" "Bolz_055_Plate1_Strip7_G07_55__S115"
## [19] "Bolz_054_Plate1_Strip7_F07_54__S114" "Bolz_053_Plate1_Strip7_E07_53__S113"
## [21] "Bolz_052_Plate1_Strip7_D07_52__S112" "Bolz_051_Plate1_Strip7_C07_51__S111"
## [23] "Bolz_050_Plate1_Strip7_B07_50__S110" "Bolz_049_Plate1_Strip7_A07_49__S109"
## [25] "Bolz_048_Plate1_Strip6_H06_48__S108" "Bolz_047_Plate1_Strip6_G06_47__S107"
## [27] "Bolz_046_Plate1_Strip6_F06_46__S106" "Bolz_045_Plate1_Strip6_E06_45__S105"
## [29] "Bolz_044_Plate1_Strip6_D06_44__S104" "Bolz_043_Plate1_Strip6_C06_43__S103"
## [31] "Bolz_042_Plate1_Strip6_B06_42__S102" "Bolz_041_Plate1_Strip6_A06_41__S101"
## [33] "Bolz_040_Plate1_Strip5_H05_40__S100" "Bolz_039_Plate1_Strip5_G05_39__S99"
## [35] "Bolz_038_Plate1_Strip5_F05_38__S98" "Bolz_037_Plate1_Strip5_E05_37__S97"
## [37] "Bolz_036_Plate1_Strip5_D05_36__S96" "Bolz_035_Plate1_Strip5_C05_35__S95"
## [39] "Bolz_034_Plate1_Strip5_B05_34__S94" "Bolz_033_Plate1_Strip5_A05_33__S93"
## [41] "Bolz_032_Plate1_Strip4_H04_32__S92" "Bolz_031_Plate1_Strip4_G04_31__S91"
## [43] "Bolz_030_Plate1_Strip4_F04_30__S90" "Bolz_029_Plate1_Strip4_E04_29__S89"
## [45] "Bolz_028_Plate1_Strip4_D04_28__S88" "Bolz_027_Plate1_Strip4_C04_27__S87"
## [47] "Bolz_026_Plate1_Strip4_B04_26__S86" "Bolz_025_Plate1_Strip4_A04_25__S85"
## [49] "Bolz_024_Plate1_Strip3_H03_24__S84" "Bolz_023_Plate1_Strip3_G03_23__S83"
## [51] "Bolz_022_Plate1_Strip3_F03_22__S82" "Bolz_021_Plate1_Strip3_E03_21__S81"
## [53] "Bolz_020_Plate1_Strip3_D03_20__S80" "Bolz_019_Plate1_Strip3_C03_19__S79"
## [55] "Bolz_018_Plate1_Strip3_B03_18__S78" "Bolz_017_Plate1_Strip3_A03_17__S77"
## [57] "Bolz_016_Plate1_Strip2_H02_16__S76" "Bolz_015_Plate1_Strip2_G02_15__S75"
## [59] "Bolz_014_Plate1_Strip2_F02_14__S74" "Bolz_013_Plate1_Strip2_E02_13__S73"
## [61] "Bolz_012_Plate1_Strip2_D02_12__S72" "Bolz_011_Plate1_Strip2_C02_11__S71"
## [63] "Bolz_010_Plate1_Strip2_B02_10__S70" "Bolz_009_Plate1_Strip2_A02_9__S69"
## [65] "Bolz_008_Plate1_Strip1_H01_8__S68" "Bolz_007_Plate1_Strip1_G01_7__S67"
## [67] "Bolz_006_Plate1_Strip1_F01_6__S66" "Bolz_005_Plate1_Strip1_E01_5__S65"
## [69] "Bolz_004_Plate1_Strip1_D01_4__S64" "Bolz_003_Plate1_Strip1_C01_3__S63"
## [71] "Bolz_002_Plate1_Strip1_B01_2__S62" "Bolz_001_Plate1_Strip1_A01_1__S61"Let’s also modify the column names and sample information, before running DESeq2
# Let's also add columns addressing their circadian status (time of the day)
# Load the .csv file including the timepoint data
tpdata <- read_csv("Sample_Info(Sept.2023).csv", col_names = F)## Rows: 72 Columns: 4
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): X2, X4
## dbl (2): X1, X3
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.# 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
head(counts)head(coldata)# At this point, we can generate the DESeqDataSetFromMatrix
dds <- DESeqDataSetFromMatrix(countData = counts,
colData = coldata,
design= ~ sex)## converting counts to integer mode## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors# Verify colnames of counts equals rownames of colData
summary(rownames(dds@colData) == colnames(dds@assays@data$counts))## Mode TRUE
## logical 72I think we can write the prefiltering step here too, where we follow the original vignette!
# We will filter out genes that have less than 10 reads.
# We also specify the minimum number of samples to be 36, which is the number of M/F samples
smallestGroupSize <- 36
keep <- rowSums(counts(dds) >= 10) >= smallestGroupSize
dds <- dds[keep,]
# Check the dimensions of the count matrix
dim(dds)## [1] 14794 72# At this point, I run DESeq. Compare before and after running DESeq
dds <- DESeq(dds)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testing## -- replacing outliers and refitting for 19 genes
## -- DESeq argument 'minReplicatesForReplace' = 7
## -- original counts are preserved in counts(dds)## estimating dispersions## fitting model and testing# Save the dds R object as a RDS file
# saveRDS(dds, "DESeq2_object_20231106.rds")
not_normalized_counts <- counts(dds, normalized = FALSE)
normalized_counts <- counts(dds, normalized = TRUE)
head(normalized_counts)## female_ZT23_072 male_ZT23_071 female_ZT19_070 male_ZT19_069
## Pcmtd1 2965.943602 2777.74362 3511.11763 2698.23369
## Uxs1 728.559839 712.79383 683.66255 764.57185
## Cdh19 57.579729 75.96382 125.50754 94.35142
## Dsel 585.198064 462.11323 586.43840 486.94008
## Asic4 8.225676 17.72489 10.60627 11.92949
## Rrs1 493.540536 591.25172 480.81762 593.22099
## female_ZT15_068 male_ZT15_067 female_ZT11_066 male_ZT11_065
## Pcmtd1 3287.073967 2948.79726 2912.5425 2716.55764
## Uxs1 726.251507 718.47139 701.5077 887.70572
## Cdh19 120.907576 88.53202 134.9760 74.22289
## Dsel 603.731830 543.67740 526.1308 618.52405
## Asic4 4.836303 28.37565 24.7915 11.87566
## Rrs1 572.295860 505.08652 513.2759 513.62237
## female_ZT7_064 male_ZT7_063 female_ZT3_062 male_ZT3_061 female_ZT23_060
## Pcmtd1 2899.305467 2984.50678 3348.95022 2991.33085 3256.40517
## Uxs1 675.152952 920.20769 849.16801 827.28846 680.37284
## Cdh19 98.333088 57.57009 131.23506 58.82097 71.55645
## Dsel 587.971039 687.18588 774.77809 617.62017 591.22054
## Asic4 7.096202 10.96573 12.63225 13.28215 39.88393
## Rrs1 489.637951 598.54622 514.41335 592.00458 480.95321
## male_ZT23_059 female_ZT19_058 male_ZT19_057 female_ZT15_056
## Pcmtd1 3099.867326 3368.02353 3411.652814 4172.546688
## Uxs1 781.659118 719.34883 763.212884 689.117496
## Cdh19 49.255232 218.85426 111.491922 80.662597
## Dsel 525.746065 670.01693 565.568114 701.067510
## Asic4 9.636893 17.04193 5.067815 9.958345
## Rrs1 584.638190 458.33822 508.808590 460.075554
## male_ZT15_055 female_ZT11_054 male_ZT11_053 female_ZT7_052 male_ZT7_051
## Pcmtd1 2951.41610 3893.43054 2964.561341 3401.5691 3281.02774
## Uxs1 700.46601 767.56202 683.369080 722.0522 652.75183
## Cdh19 67.74580 128.43264 63.023188 193.0681 79.43541
## Dsel 550.91396 720.03183 521.485988 618.2642 591.73623
## Asic4 16.61689 29.32714 7.414493 13.3920 10.36114
## Rrs1 483.16816 566.31717 511.599998 445.2841 556.04786
## female_ZT3_050 male_ZT3_049 female_ZT23_048 male_ZT23_047
## Pcmtd1 3325.56279 3078.940456 3832.7266 3124.081831
## Uxs1 691.02603 764.341575 787.4233 785.609297
## Cdh19 64.78369 47.257678 170.9833 66.079287
## Dsel 509.37765 626.677904 600.2415 555.555484
## Asic4 24.13510 4.109363 14.3986 8.565833
## Rrs1 532.24248 652.361425 452.6559 578.805604
## female_ZT19_046 male_ZT19_045 female_ZT15_044 male_ZT15_043
## Pcmtd1 3879.28549 3119.21749 3931.33773 3099.95279
## Uxs1 734.85571 777.23076 762.75336 758.29440
## Cdh19 104.97939 288.24452 176.84218 135.81392
## Dsel 612.04964 609.43127 592.71281 465.16269
## Asic4 15.84595 19.55945 15.54657 10.18604
## Rrs1 455.57093 504.42791 491.66013 437.99990
## female_ZT11_042 male_ZT11_041 female_ZT7_040 male_ZT7_039 female_ZT3_038
## Pcmtd1 3096.96218 3319.355619 3299.856461 4021.565867 3392.9820
## Uxs1 818.84506 704.916919 844.916488 863.255207 764.4523
## Cdh19 199.07701 119.964764 42.565062 119.112027 52.8542
## Dsel 631.97561 673.190701 563.987077 670.465638 721.9111
## Asic4 18.78085 4.957222 6.384759 7.367754 19.3369
## Rrs1 517.41242 504.645164 508.652496 475.220150 542.7224
## male_ZT3_037 female_ZT23_036 male_ZT23_035 female_ZT19_034 male_ZT19_033
## Pcmtd1 3597.804900 3665.62826 3865.082049 3712.78906 3201.232914
## Uxs1 799.991639 718.16390 789.843511 810.35196 866.562392
## Cdh19 81.984041 87.77559 106.457169 112.28406 83.290948
## Dsel 576.477254 694.22511 624.148079 621.79948 694.091236
## Asic4 6.903919 10.97195 1.717051 10.59284 4.206614
## Rrs1 641.201497 530.64333 513.398282 517.98969 551.066375
## female_ZT15_032 male_ZT15_031 female_ZT11_030 male_ZT11_029
## Pcmtd1 3223.67973 2856.60151 3626.24187 2720.407970
## Uxs1 751.27230 681.42832 665.31049 705.428225
## Cdh19 119.90645 99.85286 37.16149 70.419280
## Dsel 651.52711 532.99839 547.83224 492.934959
## Asic4 23.34462 16.19236 20.37888 2.470852
## Rrs1 548.59856 491.16814 471.11175 491.699533
## female_ZT7_028 male_ZT7_027 female_ZT3_026 male_ZT3_025 female_ZT23_024
## Pcmtd1 3697.10758 3301.130524 3381.7203 3207.25405 3682.954672
## Uxs1 694.40492 858.698983 793.2836 829.17543 872.156395
## Cdh19 45.49549 68.047844 114.9686 73.54882 109.716903
## Dsel 566.29919 732.324416 669.2817 822.17078 686.195544
## Asic4 15.56425 6.480747 17.2453 12.25814 8.368238
## Rrs1 475.30820 563.014899 451.6625 495.57898 535.567254
## male_ZT23_023 female_ZT19_022 male_ZT19_021 female_ZT15_020
## Pcmtd1 3559.85476 4439.13137 3679.231652 3061.27958
## Uxs1 685.23491 723.53903 776.087927 698.25105
## Cdh19 98.03216 113.65077 135.174474 94.32514
## Dsel 659.48909 632.39337 647.905236 600.25090
## Asic4 15.84358 15.75357 9.322377 15.92502
## Rrs1 515.90663 502.98903 521.276275 607.60091
## male_ZT15_019 female_ZT11_018 male_ZT11_017 female_ZT7_016 male_ZT7_015
## Pcmtd1 3530.66548 3547.356586 2961.540379 2647.51004 2977.32591
## Uxs1 779.02333 817.913824 783.579295 732.89802 826.32787
## Cdh19 121.72240 82.710387 135.057922 284.92251 148.81734
## Dsel 704.55786 811.480794 546.315378 510.17098 625.62027
## Asic4 22.91245 5.514026 8.517166 31.09774 22.51841
## Rrs1 476.14937 487.072277 585.250995 545.47112 493.44698
## female_ZT3_014 male_ZT3_013 female_ZT23_012 male_ZT23_011
## Pcmtd1 3487.48997 3328.791968 3633.75146 4095.21714
## Uxs1 714.62452 837.055179 768.75549 764.17928
## Cdh19 167.18748 93.096079 111.54491 171.89580
## Dsel 668.74990 669.482237 641.13203 686.69257
## Asic4 15.29154 4.047656 20.09818 16.03173
## Rrs1 515.83453 416.098996 567.77367 424.84093
## female_ZT19_010 male_ZT19_009 female_ZT15_008 male_ZT15_007
## Pcmtd1 3439.99527 3440.485090 3235.74617 3208.96173
## Uxs1 709.26808 851.027058 830.11851 770.22005
## Cdh19 247.83828 112.002436 81.51034 233.66226
## Dsel 696.65086 595.431755 464.39446 715.69885
## Asic4 16.22214 9.572858 10.72505 13.84665
## Rrs1 492.97286 552.353895 670.31533 516.65322
## female_ZT11_006 male_ZT11_005 female_ZT7_004 male_ZT7_003 female_ZT3_002
## Pcmtd1 3361.396049 3437.84252 3284.85245 2406.86214 2930.63765
## Uxs1 817.514670 676.61747 648.16446 690.59921 787.26251
## Cdh19 193.526703 132.97684 129.01852 79.72596 109.20619
## Dsel 614.944663 593.50694 674.78733 525.76040 626.71311
## Asic4 6.330313 21.51096 20.47913 14.00591 22.00423
## Rrs1 496.477382 510.39642 513.00221 493.43906 460.45892
## male_ZT3_001
## Pcmtd1 3368.71992
## Uxs1 808.49278
## Cdh19 98.95615
## Dsel 645.32041
## Asic4 10.52725
## Rrs1 488.46439male_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)## male_ZT3_001 male_ZT7_003 male_ZT11_005 male_ZT15_007 male_ZT19_009
## Pcmtd1 3368.71992 2406.86214 3437.84252 3208.96173 3440.485090
## Uxs1 808.49278 690.59921 676.61747 770.22005 851.027058
## Cdh19 98.95615 79.72596 132.97684 233.66226 112.002436
## Dsel 645.32041 525.76040 593.50694 715.69885 595.431755
## Asic4 10.52725 14.00591 21.51096 13.84665 9.572858
## Rrs1 488.46439 493.43906 510.39642 516.65322 552.353895
## male_ZT23_011 male_ZT3_013 male_ZT7_015 male_ZT11_017 male_ZT15_019
## Pcmtd1 4095.21714 3328.791968 2977.32591 2961.540379 3530.66548
## Uxs1 764.17928 837.055179 826.32787 783.579295 779.02333
## Cdh19 171.89580 93.096079 148.81734 135.057922 121.72240
## Dsel 686.69257 669.482237 625.62027 546.315378 704.55786
## Asic4 16.03173 4.047656 22.51841 8.517166 22.91245
## Rrs1 424.84093 416.098996 493.44698 585.250995 476.14937
## male_ZT19_021 male_ZT23_023 male_ZT3_025 male_ZT7_027 male_ZT11_029
## Pcmtd1 3679.231652 3559.85476 3207.25405 3301.130524 2720.407970
## Uxs1 776.087927 685.23491 829.17543 858.698983 705.428225
## Cdh19 135.174474 98.03216 73.54882 68.047844 70.419280
## Dsel 647.905236 659.48909 822.17078 732.324416 492.934959
## Asic4 9.322377 15.84358 12.25814 6.480747 2.470852
## Rrs1 521.276275 515.90663 495.57898 563.014899 491.699533
## male_ZT15_031 male_ZT19_033 male_ZT23_035 male_ZT3_037 male_ZT7_039
## Pcmtd1 2856.60151 3201.232914 3865.082049 3597.804900 4021.565867
## Uxs1 681.42832 866.562392 789.843511 799.991639 863.255207
## Cdh19 99.85286 83.290948 106.457169 81.984041 119.112027
## Dsel 532.99839 694.091236 624.148079 576.477254 670.465638
## Asic4 16.19236 4.206614 1.717051 6.903919 7.367754
## Rrs1 491.16814 551.066375 513.398282 641.201497 475.220150
## male_ZT11_041 male_ZT15_043 male_ZT19_045 male_ZT23_047 male_ZT3_049
## Pcmtd1 3319.355619 3099.95279 3119.21749 3124.081831 3078.940456
## Uxs1 704.916919 758.29440 777.23076 785.609297 764.341575
## Cdh19 119.964764 135.81392 288.24452 66.079287 47.257678
## Dsel 673.190701 465.16269 609.43127 555.555484 626.677904
## Asic4 4.957222 10.18604 19.55945 8.565833 4.109363
## Rrs1 504.645164 437.99990 504.42791 578.805604 652.361425
## male_ZT7_051 male_ZT11_053 male_ZT15_055 male_ZT19_057 male_ZT23_059
## Pcmtd1 3281.02774 2964.561341 2951.41610 3411.652814 3099.867326
## Uxs1 652.75183 683.369080 700.46601 763.212884 781.659118
## Cdh19 79.43541 63.023188 67.74580 111.491922 49.255232
## Dsel 591.73623 521.485988 550.91396 565.568114 525.746065
## Asic4 10.36114 7.414493 16.61689 5.067815 9.636893
## Rrs1 556.04786 511.599998 483.16816 508.808590 584.638190
## male_ZT3_061 male_ZT7_063 male_ZT11_065 male_ZT15_067 male_ZT19_069
## Pcmtd1 2991.33085 2984.50678 2716.55764 2948.79726 2698.23369
## Uxs1 827.28846 920.20769 887.70572 718.47139 764.57185
## Cdh19 58.82097 57.57009 74.22289 88.53202 94.35142
## Dsel 617.62017 687.18588 618.52405 543.67740 486.94008
## Asic4 13.28215 10.96573 11.87566 28.37565 11.92949
## Rrs1 592.00458 598.54622 513.62237 505.08652 593.22099
## male_ZT23_071
## Pcmtd1 2777.74362
## Uxs1 712.79383
## Cdh19 75.96382
## Dsel 462.11323
## Asic4 17.72489
## Rrs1 591.25172female_count <- counts(female_dds, normalized = TRUE)
female_count <- female_count[, ncol(female_dds@assays@data$counts):1]
head(female_count)## female_ZT3_002 female_ZT7_004 female_ZT11_006 female_ZT15_008
## Pcmtd1 2930.63765 3284.85245 3361.396049 3235.74617
## Uxs1 787.26251 648.16446 817.514670 830.11851
## Cdh19 109.20619 129.01852 193.526703 81.51034
## Dsel 626.71311 674.78733 614.944663 464.39446
## Asic4 22.00423 20.47913 6.330313 10.72505
## Rrs1 460.45892 513.00221 496.477382 670.31533
## female_ZT19_010 female_ZT23_012 female_ZT3_014 female_ZT7_016
## Pcmtd1 3439.99527 3633.75146 3487.48997 2647.51004
## Uxs1 709.26808 768.75549 714.62452 732.89802
## Cdh19 247.83828 111.54491 167.18748 284.92251
## Dsel 696.65086 641.13203 668.74990 510.17098
## Asic4 16.22214 20.09818 15.29154 31.09774
## Rrs1 492.97286 567.77367 515.83453 545.47112
## female_ZT11_018 female_ZT15_020 female_ZT19_022 female_ZT23_024
## Pcmtd1 3547.356586 3061.27958 4439.13137 3682.954672
## Uxs1 817.913824 698.25105 723.53903 872.156395
## Cdh19 82.710387 94.32514 113.65077 109.716903
## Dsel 811.480794 600.25090 632.39337 686.195544
## Asic4 5.514026 15.92502 15.75357 8.368238
## Rrs1 487.072277 607.60091 502.98903 535.567254
## female_ZT3_026 female_ZT7_028 female_ZT11_030 female_ZT15_032
## Pcmtd1 3381.7203 3697.10758 3626.24187 3223.67973
## Uxs1 793.2836 694.40492 665.31049 751.27230
## Cdh19 114.9686 45.49549 37.16149 119.90645
## Dsel 669.2817 566.29919 547.83224 651.52711
## Asic4 17.2453 15.56425 20.37888 23.34462
## Rrs1 451.6625 475.30820 471.11175 548.59856
## female_ZT19_034 female_ZT23_036 female_ZT3_038 female_ZT7_040
## Pcmtd1 3712.78906 3665.62826 3392.9820 3299.856461
## Uxs1 810.35196 718.16390 764.4523 844.916488
## Cdh19 112.28406 87.77559 52.8542 42.565062
## Dsel 621.79948 694.22511 721.9111 563.987077
## Asic4 10.59284 10.97195 19.3369 6.384759
## Rrs1 517.98969 530.64333 542.7224 508.652496
## female_ZT11_042 female_ZT15_044 female_ZT19_046 female_ZT23_048
## Pcmtd1 3096.96218 3931.33773 3879.28549 3832.7266
## Uxs1 818.84506 762.75336 734.85571 787.4233
## Cdh19 199.07701 176.84218 104.97939 170.9833
## Dsel 631.97561 592.71281 612.04964 600.2415
## Asic4 18.78085 15.54657 15.84595 14.3986
## Rrs1 517.41242 491.66013 455.57093 452.6559
## female_ZT3_050 female_ZT7_052 female_ZT11_054 female_ZT15_056
## Pcmtd1 3325.56279 3401.5691 3893.43054 4172.546688
## Uxs1 691.02603 722.0522 767.56202 689.117496
## Cdh19 64.78369 193.0681 128.43264 80.662597
## Dsel 509.37765 618.2642 720.03183 701.067510
## Asic4 24.13510 13.3920 29.32714 9.958345
## Rrs1 532.24248 445.2841 566.31717 460.075554
## female_ZT19_058 female_ZT23_060 female_ZT3_062 female_ZT7_064
## Pcmtd1 3368.02353 3256.40517 3348.95022 2899.305467
## Uxs1 719.34883 680.37284 849.16801 675.152952
## Cdh19 218.85426 71.55645 131.23506 98.333088
## Dsel 670.01693 591.22054 774.77809 587.971039
## Asic4 17.04193 39.88393 12.63225 7.096202
## Rrs1 458.33822 480.95321 514.41335 489.637951
## female_ZT11_066 female_ZT15_068 female_ZT19_070 female_ZT23_072
## Pcmtd1 2912.5425 3287.073967 3511.11763 2965.943602
## Uxs1 701.5077 726.251507 683.66255 728.559839
## Cdh19 134.9760 120.907576 125.50754 57.579729
## Dsel 526.1308 603.731830 586.43840 585.198064
## Asic4 24.7915 4.836303 10.60627 8.225676
## Rrs1 513.2759 572.295860 480.81762 493.540536# 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)## male_ZT3_001 male_ZT3_013 male_ZT3_025 male_ZT3_037 male_ZT3_049
## Pcmtd1 3368.71992 3328.791968 3207.25405 3597.804900 3078.940456
## Uxs1 808.49278 837.055179 829.17543 799.991639 764.341575
## Cdh19 98.95615 93.096079 73.54882 81.984041 47.257678
## Dsel 645.32041 669.482237 822.17078 576.477254 626.677904
## Asic4 10.52725 4.047656 12.25814 6.903919 4.109363
## Rrs1 488.46439 416.098996 495.57898 641.201497 652.361425
## male_ZT3_061 male_ZT7_003 male_ZT7_015 male_ZT7_027 male_ZT7_039
## Pcmtd1 2991.33085 2406.86214 2977.32591 3301.130524 4021.565867
## Uxs1 827.28846 690.59921 826.32787 858.698983 863.255207
## Cdh19 58.82097 79.72596 148.81734 68.047844 119.112027
## Dsel 617.62017 525.76040 625.62027 732.324416 670.465638
## Asic4 13.28215 14.00591 22.51841 6.480747 7.367754
## Rrs1 592.00458 493.43906 493.44698 563.014899 475.220150
## male_ZT7_051 male_ZT7_063 male_ZT11_005 male_ZT11_017 male_ZT11_029
## Pcmtd1 3281.02774 2984.50678 3437.84252 2961.540379 2720.407970
## Uxs1 652.75183 920.20769 676.61747 783.579295 705.428225
## Cdh19 79.43541 57.57009 132.97684 135.057922 70.419280
## Dsel 591.73623 687.18588 593.50694 546.315378 492.934959
## Asic4 10.36114 10.96573 21.51096 8.517166 2.470852
## Rrs1 556.04786 598.54622 510.39642 585.250995 491.699533
## male_ZT11_041 male_ZT11_053 male_ZT11_065 male_ZT15_007 male_ZT15_019
## Pcmtd1 3319.355619 2964.561341 2716.55764 3208.96173 3530.66548
## Uxs1 704.916919 683.369080 887.70572 770.22005 779.02333
## Cdh19 119.964764 63.023188 74.22289 233.66226 121.72240
## Dsel 673.190701 521.485988 618.52405 715.69885 704.55786
## Asic4 4.957222 7.414493 11.87566 13.84665 22.91245
## Rrs1 504.645164 511.599998 513.62237 516.65322 476.14937
## male_ZT15_031 male_ZT15_043 male_ZT15_055 male_ZT15_067 male_ZT19_009
## Pcmtd1 2856.60151 3099.95279 2951.41610 2948.79726 3440.485090
## Uxs1 681.42832 758.29440 700.46601 718.47139 851.027058
## Cdh19 99.85286 135.81392 67.74580 88.53202 112.002436
## Dsel 532.99839 465.16269 550.91396 543.67740 595.431755
## Asic4 16.19236 10.18604 16.61689 28.37565 9.572858
## Rrs1 491.16814 437.99990 483.16816 505.08652 552.353895
## male_ZT19_021 male_ZT19_033 male_ZT19_045 male_ZT19_057 male_ZT19_069
## Pcmtd1 3679.231652 3201.232914 3119.21749 3411.652814 2698.23369
## Uxs1 776.087927 866.562392 777.23076 763.212884 764.57185
## Cdh19 135.174474 83.290948 288.24452 111.491922 94.35142
## Dsel 647.905236 694.091236 609.43127 565.568114 486.94008
## Asic4 9.322377 4.206614 19.55945 5.067815 11.92949
## Rrs1 521.276275 551.066375 504.42791 508.808590 593.22099
## male_ZT23_011 male_ZT23_023 male_ZT23_035 male_ZT23_047 male_ZT23_059
## Pcmtd1 4095.21714 3559.85476 3865.082049 3124.081831 3099.867326
## Uxs1 764.17928 685.23491 789.843511 785.609297 781.659118
## Cdh19 171.89580 98.03216 106.457169 66.079287 49.255232
## Dsel 686.69257 659.48909 624.148079 555.555484 525.746065
## Asic4 16.03173 15.84358 1.717051 8.565833 9.636893
## Rrs1 424.84093 515.90663 513.398282 578.805604 584.638190
## male_ZT23_071
## Pcmtd1 2777.74362
## Uxs1 712.79383
## Cdh19 75.96382
## Dsel 462.11323
## Asic4 17.72489
## Rrs1 591.25172head(female_count_modified)## female_ZT3_002 female_ZT3_014 female_ZT3_026 female_ZT3_038
## Pcmtd1 2930.63765 3487.48997 3381.7203 3392.9820
## Uxs1 787.26251 714.62452 793.2836 764.4523
## Cdh19 109.20619 167.18748 114.9686 52.8542
## Dsel 626.71311 668.74990 669.2817 721.9111
## Asic4 22.00423 15.29154 17.2453 19.3369
## Rrs1 460.45892 515.83453 451.6625 542.7224
## female_ZT3_050 female_ZT3_062 female_ZT7_004 female_ZT7_016
## Pcmtd1 3325.56279 3348.95022 3284.85245 2647.51004
## Uxs1 691.02603 849.16801 648.16446 732.89802
## Cdh19 64.78369 131.23506 129.01852 284.92251
## Dsel 509.37765 774.77809 674.78733 510.17098
## Asic4 24.13510 12.63225 20.47913 31.09774
## Rrs1 532.24248 514.41335 513.00221 545.47112
## female_ZT7_028 female_ZT7_040 female_ZT7_052 female_ZT7_064
## Pcmtd1 3697.10758 3299.856461 3401.5691 2899.305467
## Uxs1 694.40492 844.916488 722.0522 675.152952
## Cdh19 45.49549 42.565062 193.0681 98.333088
## Dsel 566.29919 563.987077 618.2642 587.971039
## Asic4 15.56425 6.384759 13.3920 7.096202
## Rrs1 475.30820 508.652496 445.2841 489.637951
## female_ZT11_006 female_ZT11_018 female_ZT11_030 female_ZT11_042
## Pcmtd1 3361.396049 3547.356586 3626.24187 3096.96218
## Uxs1 817.514670 817.913824 665.31049 818.84506
## Cdh19 193.526703 82.710387 37.16149 199.07701
## Dsel 614.944663 811.480794 547.83224 631.97561
## Asic4 6.330313 5.514026 20.37888 18.78085
## Rrs1 496.477382 487.072277 471.11175 517.41242
## female_ZT11_054 female_ZT11_066 female_ZT15_008 female_ZT15_020
## Pcmtd1 3893.43054 2912.5425 3235.74617 3061.27958
## Uxs1 767.56202 701.5077 830.11851 698.25105
## Cdh19 128.43264 134.9760 81.51034 94.32514
## Dsel 720.03183 526.1308 464.39446 600.25090
## Asic4 29.32714 24.7915 10.72505 15.92502
## Rrs1 566.31717 513.2759 670.31533 607.60091
## female_ZT15_032 female_ZT15_044 female_ZT15_056 female_ZT15_068
## Pcmtd1 3223.67973 3931.33773 4172.546688 3287.073967
## Uxs1 751.27230 762.75336 689.117496 726.251507
## Cdh19 119.90645 176.84218 80.662597 120.907576
## Dsel 651.52711 592.71281 701.067510 603.731830
## Asic4 23.34462 15.54657 9.958345 4.836303
## Rrs1 548.59856 491.66013 460.075554 572.295860
## female_ZT19_010 female_ZT19_022 female_ZT19_034 female_ZT19_046
## Pcmtd1 3439.99527 4439.13137 3712.78906 3879.28549
## Uxs1 709.26808 723.53903 810.35196 734.85571
## Cdh19 247.83828 113.65077 112.28406 104.97939
## Dsel 696.65086 632.39337 621.79948 612.04964
## Asic4 16.22214 15.75357 10.59284 15.84595
## Rrs1 492.97286 502.98903 517.98969 455.57093
## female_ZT19_058 female_ZT19_070 female_ZT23_012 female_ZT23_024
## Pcmtd1 3368.02353 3511.11763 3633.75146 3682.954672
## Uxs1 719.34883 683.66255 768.75549 872.156395
## Cdh19 218.85426 125.50754 111.54491 109.716903
## Dsel 670.01693 586.43840 641.13203 686.195544
## Asic4 17.04193 10.60627 20.09818 8.368238
## Rrs1 458.33822 480.81762 567.77367 535.567254
## female_ZT23_036 female_ZT23_048 female_ZT23_060 female_ZT23_072
## Pcmtd1 3665.62826 3832.7266 3256.40517 2965.943602
## Uxs1 718.16390 787.4233 680.37284 728.559839
## Cdh19 87.77559 170.9833 71.55645 57.579729
## Dsel 694.22511 600.2415 591.22054 585.198064
## Asic4 10.97195 14.3986 39.88393 8.225676
## Rrs1 530.64333 452.6559 480.95321 493.540536#
# 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)# Let's try using the MetaCycle package at this point!
# MetaCycle::meta2d(infile = "DESeq2_normalized_male_counts_20231106.csv",
# filestyle = "csv",
# timepoints = rep(seq(3, 23, by=4), each=6),
# outdir = "metaout_20231106")
# MetaCycle::meta2d(infile = "DESeq2_normalized_female_counts_20231106.csv",
# filestyle = "csv",
# timepoints = rep(seq(3, 23, by=4), each=6),
# outdir = "metaout_20231106")
male_metaCycle <- read.csv("~/Circadian_Rhythms/metaout_20231106/meta2d_DESeq2_normalized_male_counts_20231106.csv")
female_metaCycle <- read.csv("~/Circadian_Rhythms/metaout_20231106/meta2d_DESeq2_normalized_female_counts_20231106.csv")
male_cycling <- male_metaCycle %>%
filter(meta2d_pvalue < 0.05) %>% pull(CycID)
female_cycling <- female_metaCycle %>%
filter(meta2d_pvalue < 0.05) %>% pull(CycID)
length(male_cycling)## [1] 1287length(female_cycling)## [1] 1721# Find the common genes that cycle between males and females
common_cycling <- intersect(male_cycling, female_cycling)
"Per2" %in% common_cycling## [1] TRUE"Bmal1" %in% common_cycling## [1] TRUE"Clock" %in% common_cycling## [1] TRUE"Npas2" %in% common_cycling## [1] TRUE"Cry1" %in% common_cycling## [1] TRUEI think it’s a good idea to check for adjusted p-values here too.
male_cycling_BH <- male_metaCycle %>%
filter(meta2d_BH.Q < 0.05) %>% pull(CycID)
female_cycling_BH <- female_metaCycle %>%
filter(meta2d_BH.Q < 0.05) %>% pull(CycID)
length(male_cycling_BH)## [1] 537length(female_cycling_BH)## [1] 889# Find the common genes that cycle between males and females
common_cycling_BH <- intersect(male_cycling_BH, female_cycling_BH)
length(common_cycling_BH)## [1] 323"Per2" %in% common_cycling_BH## [1] TRUE"Bmal1" %in% common_cycling_BH## [1] TRUE"Clock" %in% common_cycling_BH## [1] TRUE"Npas2" %in% common_cycling_BH## [1] TRUE"Cry1" %in% common_cycling_BH## [1] TRUEmale_phase <- male_metaCycle %>%
filter(CycID %in% common_cycling_BH) %>%
pull(meta2d_phase)
female_phase <- female_metaCycle %>%
filter(CycID %in% common_cycling_BH) %>%
pull(meta2d_phase)
test1 <- female_metaCycle %>%
filter(CycID %in% common_cycling_BH)
test1[(abs(male_phase - female_phase) > 6) & (abs(male_phase - female_phase) < 18),] %>%
pull(CycID) %>% cat()## Ttc9ctest2 <- male_metaCycle %>%
filter(CycID %in% common_cycling)
test2[(abs(male_phase - female_phase) > 6) & (abs(male_phase - female_phase) < 18),] %>%
pull(CycID) %>% cat()## Tubgcp2 Cxxc1# Check for the results from DESeq2
resultsNames(dds)## [1] "Intercept" "sex_male_vs_female"res <- results(dds, name = "sex_male_vs_female")
resLFC <- lfcShrink(dds, coef = "sex_male_vs_female", type = "apeglm")## using 'apeglm' for LFC shrinkage. If used in published research, please cite:
## Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
## sequence count data: removing the noise and preserving large differences.
## Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895# 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.
# 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"
)## Warning: One or more p-values is 0. Converting to 10^-1 * current lowest
## non-zero p-value...