# Dec 11, 2024 # R 4.2.2 # DESeq2 on incubation data, WS samples # CH4 as 3 groups: 0.1%, 1%, 5% # salinity, temp and season as confounding/additional factors # clear workspace rm(list = ls()) # Install BiocManager if not already installed if (!requireNamespace("BiocManager", quietly = TRUE)) { install.packages("BiocManager") } if (!requireNamespace("ashr", quietly = TRUE)) { install.packages("ashr") } # Install packages required_packages <- c("DESeq2") for (pkg in required_packages) { if (!requireNamespace(pkg, quietly = TRUE)) { BiocManager::install(pkg) } } # Load libraries library(DESeq2) ### # Load the ASV table asv_table <- read.csv(file='asv_table.csv', row.names=1) # load metadata metadata <- read.csv('filtered_metadata.csv', row.names=1, stringsAsFactors = FALSE) # reduce both to WS incubation samples idx <- which(metadata$type =='incubation' & metadata$location=='wadden sea') length(idx) # 98 all.equal(colnames(asv_table), rownames(metadata)) # TRUE, use idx on both # WS asv table asv.ws <- asv_table[,idx] # WS metadata metadata.ws <- metadata[idx,] # replicate samples have _A, _B, _C in label, remove and make condition vector cond.ws <- gsub("_[A,B,C]$", "", metadata.ws$description) # 12 samples per condition (7 different combinations), 14 for T0 # aggregate tax on genus level # paste levels up to Genus level # taxonomy of full ASV table taxo <- read.csv('taxotable.csv', row.names = 1) taxo.gen <- unlist(apply(taxo, 1, function(x){ paste(x[c(2:6)], collapse='-') })) ################### Wadden Sea incubations ###################### # aggregate asvs per genus (EXCLUDE singletons first, likely artifacts) asv1.ws <- asv.ws[rowSums(asv.ws) > 1,] nrow(asv1.ws) # 7700 ASVs ncol(asv1.ws) # 98 samples # reduce tax string accordingly tax1.gen.ws <- taxo.gen[rownames(asv1.ws)] # aggregate on genus level ws.gen <-aggregate(asv1.ws, by = list(tax1.gen.ws), FUN = sum) dim(ws.gen) # 890 98 # move names to rownames rownames(ws.gen) <- ws.gen[,1] ws.gen <- ws.gen[,-1] # double-check for singletons ws.gen <- ws.gen[rowSums(ws.gen) > 1,] # no 0 and singleton genera # If genus NA, automatically summarized at a higher level, this is ok. # rm feature with all NA ws.gen['NA-NA-NA-NA-NA',] idx <- which(rownames(ws.gen)=='NA-NA-NA-NA-NA') ws.gen <- ws.gen[-idx,] dim(ws.gen) # 889 98 # ok # CH4 for t0 were to 5%, but samples were sacrificed immediately so # set it to 0.1% as this resembles in situ conc. idx <- which(cond.ws=="T0") # change column name from ch4.hs to CH4, sal and temp for better readability colnames(metadata.ws) <- gsub('ch4..nM*', 'CH4', colnames(metadata.ws)) colnames(metadata.ws) <- gsub('temp...C.', 'temp', colnames(metadata.ws)) colnames(metadata.ws) <- gsub('sal..psu.', 'sal', colnames(metadata.ws)) metadata.ws$CH4[idx] <- 0.1 # season variables # reference is summer spring <- as.numeric(metadata.ws$season=='spring') autumn <- as.numeric(metadata.ws$season=='autumn') winter <- as.numeric(metadata.ws$season=='winter') table(metadata.ws$season) # summer has 26, the others 24 # scale salinity by dividing by standard deviation to give salinity similar weight than other explanatory variables sal.adj <- as.numeric(metadata.ws$sal) sal.adj <- sal.adj/sd(sal.adj) # collect covariates/explanatory variables in a data frame covs.sc <- data.frame('CH4'= metadata.ws$CH4, 'temp'=metadata.ws$temp, 'spring'=spring, 'autumn'=autumn, 'winter'=winter, 'sal'=sal.adj) covs.sc[, 1:5] <- apply(covs.sc[,1:5], 2, as.factor) covs.sc[, 'sal'] <- as.numeric(covs.sc[,'sal']) # only treat sal as numeric, the other coefs as factors # estimate dispersion can't handle zeros so add a pseudocount # do not use an intercept dds.sc <- DESeqDataSetFromMatrix(countData = ws.gen+1, colData = covs.sc, design = ~ CH4 + temp + sal + spring + autumn + winter) ds.sc = DESeq(dds.sc, test="Wald", fitType="local", useT=TRUE) # output tables for the different coefficients resultsNames(ds.sc) # [1] "Intercept" "CH4_1_vs_0.1" "CH4_5_vs_0.1" "temp_25_vs_15" "temp_30_vs_15" # "sal" "spring_1_vs_0" "autumn_1_vs_0" "winter_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.sc <- lfcShrink(ds.sc, coef=2, type="ashr") # "CH4_5_vs_0.1" resLFC.5perc.sc <- lfcShrink(ds.sc, coef=3, type="ashr") # "temp_25_vs_15" resLFC.temp25.sc <- lfcShrink(ds.sc, coef=4, type="ashr") # "temp_30_vs_15" resLFC.temp30.sc <- lfcShrink(ds.sc, coef=5, type="ashr") # salinity effects (per unit scaled) resLFC.sal.sc <- lfcShrink(ds.sc, coef=6, type="ashr") # "spring" resLFC.spr.sc <- lfcShrink(ds.sc, coef=7, type="ashr") # "autumn" resLFC.aut.sc <- lfcShrink(ds.sc, coef=8, type="ashr") # "winter" resLFC.win.sc <- lfcShrink(ds.sc, coef=9, type="ashr") # 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) # 18 are present in dataset # only the ones in the dataset tax.str.mob1 <- tax.str.mob[tax.str.mob %in% rownames(counts(ds.sc))] # write MOB lists to file resTodo <- ls(pattern='resLFC.*.sc$') for( i in resTodo){ out <- get(i)[tax.str.mob1,c(1:3,5)] # get the correct results object out <- apply(out, 2, signif, digits=4) write.table(out, file = paste('DESEq2_MOBs_WS_', i, '.txt', sep = ''), sep = '\t', quote = FALSE) } # sign. MOBS only for( i in resTodo){ out <- get(i)[tax.str.mob1,c(1:3,5)] # get the correct results object out <- apply(out, 2, signif, digits=4) out <- out[which(out[,'padj'] <=0.05), ] # p-val < 0.05 out <- out[order(out[, 'log2FoldChange'], decreasing=TRUE), ] # order by lfc out <- out[which(out[, 'baseMean'] >= 10), ] # baseMean >= 10 write.table(out, file=paste('DESEq2_SIGN_MOBs_WS_', i, '.txt', sep=''), 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