# Nov 09, 2023 and Nov 21, 2024 # R 4.2.2 # run DESeq2 on incubation data, ***North Sea samples*** # CH4 as 3 groups: 0.1%, 1%, 5% # temperature and season as confounding factors # run by Julia Engelmann # Clear the R workspace rm(list = ls()) set.seed(987654) # Install BiocManager if not already installed # Check if BiocManager is installed, and install it if not if (!requireNamespace("BiocManager", quietly = TRUE)) { install.packages("BiocManager") } # Install necessary Bioconductor and CRAN packages # DESeq2 is a Bioconductor package if (!requireNamespace("DESeq2", quietly = TRUE)) { BiocManager::install("DESeq2") } # pheatmap is a CRAN package if (!requireNamespace("pheatmap", quietly = TRUE)) { install.packages("pheatmap") } # compositions is a CRAN package if (!requireNamespace("compositions", quietly = TRUE)) { install.packages("compositions") } # Cairo is a CRAN package if (!requireNamespace("Cairo", quietly = TRUE)) { install.packages("Cairo") } # Load the installed libraries library(DESeq2) library(pheatmap) library(compositions) library(Cairo) # Load the ASV table asv_table <- read.csv(file='asv_table.csv', row.names=1) dim(asv_table) # Rows/ASVs: 27012 Columns/samples: 318 # read sample metadata (for 318 samples) metadata <- read.csv(file='filtered_metadata.csv', row.names=1, stringsAsFactors = FALSE) # select only incubation samples from the NS # replicates are indicated with _A, _B, _C, remove from label asv.ns <- asv_table[, (metadata$location=='north sea' & metadata$type=='incubation')] dim(asv.ns) # 67 samples, 27012 ASVs # reduce metadata to NS incubation samples metadata.ns <- metadata[metadata$location=='north sea' & metadata$type=='incubation',] # 67 samples table(metadata.ns$season) # summer has 31, autumn 36 # 12 samples per condition (5 different combinations of CH4 and temperature), 8 for T0 # save condition in vector cond.ns <- gsub("_[A,B,C]$", "", metadata.ns[,'description']) table(cond.ns) # aggregate tax on genus level # paste levels up to Genus level # taxonomy of full ASV table taxo <- read.csv(file='taxotable.csv', row.names = 1) taxo.gen <- unlist(apply(taxo, 1, function(x){ paste(x[c(2:6)], collapse='-') })) # aggregate ASVs per genus (EXCLUDE singletons first, likely artifacts!) asv1.ns <- asv.ns[rowSums(asv.ns) > 1,] nrow(asv1.ns) # 4655 ASVs remain # reduce tax string accordingly tax1.gen.ns <- taxo.gen[rownames(asv1.ns)] # aggregate on genus level ns.gen <-aggregate(asv1.ns, by = list(tax1.gen.ns), FUN = sum) dim(ns.gen) # 667 68 # move names to rownames rownames(ns.gen) <- ns.gen[,1] ns.gen <- ns.gen[,-1] # check singletons ns.gen <- ns.gen[rowSums(ns.gen) > 1,] # no 0 and singleton genera # If genus NA, automatically summarized at a higher level, this is ok. # rm feature with all NA (missing taxonomy) ns.gen['NA-NA-NA-NA-NA',] idx <- which(rownames(ns.gen)=='NA-NA-NA-NA-NA') ns.gen <- ns.gen[-idx,] dim(ns.gen) # 666 67 # set CH4 for T0 to 0.1% here (in situ conc. is very close to that), # as the 5% CH4 applied didn't have any time to cause an effect (samples sacrificed immediately) idx <- which(cond.ns=="T0") metadata.ns$ch4..nM[idx] <- 0.1 # change column name from ch4.hs to CH4, and temp for better readability colnames(metadata.ns) <- gsub('ch4..nM*', 'CH4', colnames(metadata.ns)) colnames(metadata.ns) <- gsub('temp...C.', 'temp', colnames(metadata.ns)) # season variables # reference is summer, the only other season is autumn autumn <- as.numeric(metadata.ns$season=='autumn') # covariates for DESeq2 covs <- data.frame('CH4'= metadata.ns$CH4, 'temp'=metadata.ns$temp, 'autumn'=autumn) covs <- apply(covs, 2, as.numeric) covs <- apply(covs, 2, as.factor) # estimate dispersion can't handle zeros so add a pseudocount dds <- DESeqDataSetFromMatrix(countData = ns.gen+1, colData = covs, design = ~ CH4 + temp + autumn) # run DESeq test ds = DESeq(dds, test="Wald", fitType="local", useT=TRUE) ##### moderated log2 fold changes #### resultsNames(ds) #[1] "Intercept" "CH4_1_vs_0.1" "CH4_5_vs_0.1" "temp_25_vs_15" #[5] "temp_30_vs_15" "autumn_1_vs_0" # for the different contrasts/coefficients # if lfcShrink is applied on ds, results is called first internally # "CH4_1_vs_0.1" resLFC.1perc <- lfcShrink(ds, coef=2, type="ashr") # "CH4_5_vs_0.1" resLFC.5perc <- lfcShrink(ds, coef=3, type="ashr") # "temp_25_vs_15" resLFC.temp25 <- lfcShrink(ds, coef=4, type="ashr") # "temp_30_vs_15" resLFC.temp30 <- lfcShrink(ds, coef=5, type="ashr") # "autumn" (aut vs summer) resLFC.spr <- lfcShrink(ds, coef=6, type="ashr") ######## # write all genera with padj < 0.05 to file # CH4 5% out <- resLFC.5perc[which(resLFC.5perc$padj < 0.05), c(1:3,5)] out <- apply(out, 2, signif, digits=4) write.table(out, file='DESEq2_p0.05_NS_5percCH4_20231109.txt', sep='\t', quote=FALSE) # CH4 1% out <- resLFC.1perc[which(resLFC.1perc$padj < 0.05), c(1:3,5)] out <- apply(out, 2, signif, digits=4) write.table(out, file='DESEq2_p0.05_NS_1percCH4_20231109.txt', sep='\t', quote=FALSE) # 25C out <- resLFC.temp25[which(resLFC.temp25$padj < 0.05), c(1:3,5)] out <- apply(out, 2, signif, digits=4) write.table(out, file='DESEq2_p0.05_NS_25vs15C_20231109.txt', sep='\t', quote=FALSE) # 30 C out <- resLFC.temp30[which(resLFC.temp30$padj < 0.05), c(1:3,5)] out <- apply(out, 2, signif, digits=4) write.table(out, file='DESEq2_p0.05_NS_30vs15C_20231109.txt', sep='\t', quote=FALSE) # autumn out <- resLFC.spr[which(resLFC.spr$padj < 0.05), c(1:3,5)] out <- apply(out, 2, signif, digits=4) write.table(out, file='DESEq2_p0.05_NS_autumn_20231109.txt', sep='\t', quote=FALSE) ################################ # select MOB # read MOB table mobsH <- read.csv(file='known_MOB.csv', header=TRUE) # rm empty cols at the end mobsH <- mobsH[, 1:7] # get tax of all genera used also for DESeq2 # taxo is on ASV level, collapse to unique tax strings tax.spec <- unique(taxo) dim(tax.spec) # [1] 1921 8 # are there leading or trailing spaces that impair mapping? no grep('^\ ', tax.spec[, 'Family']) # extract tax levels (genera and families) in Helge's MOB list # need to get *all* hits, not only the first, genus can have multiple species! gens <- tax.spec[,'Genus'] %in% mobsH[, 'Genus'] sum(gens) # 25 found tax.mob.gen <- tax.spec[gens,] # add families that are completely MOB # need to get *all* hits, not the first! fam <- unique(mobsH[,'Family']) # check note col from Helge by eye, ok fam <- c("Methylomonadaceae", "Methylococcaceae", "Methylohalobiaceae", "Methylacidiphilaceae") fams <- tax.spec[,'Family'] %in% fam tax.mob.fam <- tax.spec[fams,] dim(tax.mob.fam) # 22 fam %in% tax.spec[, 'Family'] # combine genus and fam level, then take unique tax strings on genus level tax.mob <- rbind(tax.mob.gen, tax.mob.fam) tax.mob <- unique(tax.mob) dim(tax.mob) # 27 8 # family-NA features: # Methylomonadaceae # Methylacidiphilaceae # make tax strings Phylum up to genus level to match with DESeq2 tables. tax.str.mob <- unlist(apply(tax.mob, 1, function(x){ # 25 genera, 2 families with NA genus paste(x[c(2:6)], collapse='-') })) # rm duplicates, bc there can be multiple species in one genus! tax.str.mob <- unique(tax.str.mob) length(tax.str.mob) # only the ones in the dataset tax.str.mob1 <- tax.str.mob[tax.str.mob %in% rownames(counts(ds))] ##### # read DESeq2 lists and write MOBs to separate file # NS lists with all genera p.adj <= 0.05 files <- c('DESEq2_p0.05_NS_5percCH4_20231109.txt', 'DESEq2_p0.05_NS_1percCH4_20231109.txt', 'DESEq2_p0.05_NS_25vs15C_20231109.txt', 'DESEq2_p0.05_NS_30vs15C_20231109.txt', 'DESEq2_p0.05_NS_autumn_20231109.txt') for(f in files){ inFile <- read.delim(f) idx <- match(tax.str.mob1, rownames(inFile)) idx <- idx[!is.na(idx)] mob.out <- inFile[idx, ] mob.out <- mob.out[order(mob.out$log2FoldChange, decreasing=TRUE),] # order by lfc mob.out <- mob.out[which(mob.out[, 'baseMean'] >= 10), ] # keep if baseMean >= 10 outFile <- gsub('20231109', '20241121', f) # update timestamp of new files outFile <- gsub('DESEq2_', 'DESeq2_MOB_', outFile) write.table(mob.out, file=outFile, sep='\t', quote=FALSE) } sessionInfo() #R version 4.4.0 (2024-04-24 ucrt) #Platform: x86_64-w64-mingw32/x64 #Running under: Windows 11 x64 (build 22621) #Matrix products: default #locale: # [1] LC_COLLATE=English_Europe.utf8 LC_CTYPE=English_Europe.utf8 #[3] LC_MONETARY=English_Europe.utf8 LC_NUMERIC=C #[5] LC_TIME=English_Europe.utf8 #time zone: Europe/Amsterdam #tzcode source: internal #attached base packages: # [1] stats4 stats graphics grDevices utils datasets methods base #other attached packages: # [1] Cairo_1.6-2 compositions_2.0-8 pheatmap_1.0.12 #[4] DESeq2_1.44.0 SummarizedExperiment_1.34.0 Biobase_2.64.0 #[7] MatrixGenerics_1.16.0 matrixStats_1.3.0 GenomicRanges_1.56.1 #[10] GenomeInfoDb_1.40.1 IRanges_2.38.0 S4Vectors_0.42.0 #[13] BiocGenerics_0.50.0 tidyr_1.3.1 readr_2.1.5 #[16] dplyr_1.1.4 ggplot2_3.5.1 vegan_2.6-6.1 #[19] lattice_0.22-6 permute_0.9-7 #loaded via a namespace (and not attached): # [1] tidyselect_1.2.1 farver_2.1.2 tensorA_0.36.2.1 #[4] lifecycle_1.0.4 cluster_2.1.6 invgamma_1.1 #[7] magrittr_2.0.3 compiler_4.4.0 rlang_1.1.4 #[10] tools_4.4.0 utf8_1.2.4 S4Arrays_1.4.1 #[13] labeling_0.4.3 bit_4.0.5 DelayedArray_0.30.1 #[16] plyr_1.8.9 RColorBrewer_1.1-3 pkgload_1.4.0 #[19] abind_1.4-5 BiocParallel_1.38.0 withr_3.0.0 #[22] purrr_1.0.2 grid_4.4.0 fansi_1.0.6 #[25] colorspace_2.1-0 scales_1.3.0 MASS_7.3-60.2 #[28] cli_3.6.2 crayon_1.5.2 generics_0.1.3 #[31] rstudioapi_0.16.0 robustbase_0.99-2 httr_1.4.7 #[34] tzdb_0.4.0 bayesm_3.1-6 zlibbioc_1.50.0 #[37] splines_4.4.0 parallel_4.4.0 XVector_0.44.0 #[40] vctrs_0.6.5 Matrix_1.7-0 jsonlite_1.8.8 #[43] gridGraphics_0.5-1 hms_1.1.3 mixsqp_0.3-54 #[46] bit64_4.0.5 irlba_2.3.5.1 locfit_1.5-9.9 #[49] colorBlindness_0.1.9 glue_1.7.0 DEoptimR_1.1-3 #[52] codetools_0.2-20 cowplot_1.1.3 gtable_0.3.5 #[55] UCSC.utils_1.0.0 munsell_0.5.1 tibble_3.2.1 #[58] pillar_1.9.0 GenomeInfoDbData_1.2.12 truncnorm_1.0-9 #[61] R6_2.5.1 vroom_1.6.5 SQUAREM_2021.1 #[64] ashr_2.2-63 Rcpp_1.0.12 SparseArray_1.4.8 #[67] nlme_3.1-164 mgcv_1.9-1 pkgconfig_2.0.3