Create a heatmap with p-values

plotHeat(
  mat,
  pmat = NULL,
  type = "full",
  textUpp = "mat",
  textLow = "code",
  methUpp = "color",
  methLow = "color",
  diag = TRUE,
  title = "",
  mar = c(0, 0, 1, 0),
  labPos = "lt",
  labCol = "gray40",
  labCex = 1.1,
  textCol = "black",
  textCex = 1,
  textFont = 1,
  digits = 2L,
  legendPos = "r",
  colorPal = NULL,
  addWhite = TRUE,
  nCol = 51L,
  colorLim = NULL,
  revCol = FALSE,
  color = NULL,
  bg = "white",
  argsUpp = NULL,
  argsLow = NULL
)

# Load data sets from American Gut Project (from SpiecEasi package)
data("amgut2.filt.phy")

# Split data into two groups: with and without seasonal allergies
amgut_season_yes <- phyloseq::subset_samples(amgut2.filt.phy, 
                                      SEASONAL_ALLERGIES == "yes")
amgut_season_no <- phyloseq::subset_samples(amgut2.filt.phy, 
                                     SEASONAL_ALLERGIES == "no")

# Sample sizes
phyloseq::nsamples(amgut_season_yes)
#> [1] 121
phyloseq::nsamples(amgut_season_no)
#> [1] 163

# Make sample sizes equal to ensure comparability
n_yes <- phyloseq::nsamples(amgut_season_yes)

amgut_season_no <- phyloseq::subset_samples(amgut_season_no, X.SampleID %in% 
                                     get_variable(amgut_season_no, 
                                     "X.SampleID")[1:n_yes])
#> Error in h(simpleError(msg, call)): error in evaluating the argument 'table' in selecting a method for function '%in%': error in evaluating the argument 'object' in selecting a method for function 'sample_data': object 'amgut_season_no' not found

# Network construction
amgut_net <- netConstruct(data = amgut_season_yes,
                          data2 = amgut_season_no,
                          measure = "pearson",
                          filtTax = "highestVar",
                          filtTaxPar = list(highestVar = 50),
                          zeroMethod = "pseudoZO", 
                          normMethod = "clr",
                          sparsMethod = "thresh",
                          thresh = 0.4,
                          seed = 123456)
#> Checking input arguments ... 
#> Done.
#> Data filtering ...
#> 95 taxa removed in each data set.
#> 1 rows with zero sum removed in group 1.
#> 1 rows with zero sum removed in group 2.
#> 43 taxa and 120 samples remaining in group 1.
#> 43 taxa and 162 samples remaining in group 2.
#> 
#> Zero treatment in group 1:
#> Zero counts replaced by 1
#> 
#> Zero treatment in group 2:
#> Zero counts replaced by 1
#> 
#> Normalization in group 1:
#> Execute clr(){SpiecEasi} ... 
#> Done.
#> 
#> Normalization in group 2:
#> Execute clr(){SpiecEasi} ... 
#> Done.
#> 
#> Calculate 'pearson' associations ... 
#> Done.
#> 
#> Calculate associations in group 2 ... 
#> Done.
#> 
#> Sparsify associations via 'threshold' ... 
#> Done.
#> 
#> Sparsify associations in group 2 ... 
#> Done.

# Estimated and sparsified associations of group 1
plotHeat(amgut_net$assoEst1, textUpp = "none", labCex = 0.6)

plotHeat(amgut_net$assoMat1, textUpp = "none", labCex = 0.6)


# Compute graphlet correlation matrices and perform significance tests
adja1 <- amgut_net$adjaMat1
adja2 <- amgut_net$adjaMat2

gcm1 <- calcGCM(adja1)
gcm2 <- calcGCM(adja2)

gcmtest <- testGCM(obj1 = gcm1, obj2 = gcm2)
#> Perform Student's t-test for GCM1 ... 
#>  Adjust for multiple testing ... 
#> 
#>  Proportion of true null hypotheses: 0.22
#> Done.
#> 
#> Perform Student's t-test for GCM2 ... 
#>  Adjust for multiple testing ... 
#> 
#>  Proportion of true null hypotheses: 0.08
#> Done.
#> 
#> Test GCM1 and GCM2 for differences ... 
#>  Adjust for multiple testing ... 
#> 
#>  Proportion of true null hypotheses: 0.64
#> Done.

# Mixed heatmap of GCM1 and significance codes
plotHeat(mat = gcmtest$gcm1, 
         pmat = gcmtest$pAdjust1,
         type = "mixed", 
         textLow = "code")


# Mixed heatmap of GCM2 and p-values (diagonal disabled)
plotHeat(mat = gcmtest$gcm1, 
         pmat = gcmtest$pAdjust1,
         diag = FALSE,
         type = "mixed", 
         textLow = "pmat")


# Mixed heatmap of differences (GCM1 - GCM2) and significance codes
plotHeat(mat = gcmtest$diff, 
         pmat = gcmtest$pAdjustDiff,
         type = "mixed", 
         textLow = "code",
         title = "Differences between GCMs (GCM1 - GCM2)",
         mar = c(0, 0, 2, 0))


# Heatmap of differences (insignificant values are blank)
plotHeat(mat = gcmtest$diff, 
         pmat = gcmtest$pAdjustDiff,
         type = "full", 
         textUpp = "sigmat")


# Same as before but with higher significance level
plotHeat(mat = gcmtest$diff, 
         pmat = gcmtest$pAdjustDiff,
         type = "full", 
         textUpp = "sigmat",
         argsUpp = list(sig.level = 0.1))


# Heatmap of absolute differences
# (different position of labels and legend)
plotHeat(mat = gcmtest$absDiff, 
         type = "full", 
         labPos = "d",
         legendPos = "b")


# Mixed heatmap of absolute differences 
# (different methods, text options, and color palette)
plotHeat(mat = gcmtest$absDiff,
         type = "mixed",
         textLow = "mat",
         methUpp = "number",
         methLow = "circle",
         labCol = "black",
         textCol = "gray50",
         textCex = 1.3,
         textFont = 2,
         digits = 1L,
         colorLim = range(gcmtest$absDiff),
         colorPal = "Blues",
         nCol = 21L,
         bg = "darkorange",
         addWhite = FALSE)


# Mixed heatmap of differences 
# (different methods, text options, and color palette)
plotHeat(mat = gcmtest$diff,
         type = "mixed",
         textLow = "none",
         methUpp = "number",
         methLow = "pie",
         textCex = 1.3,
         textFont = 2,
         digits = 1L,
         colorLim = range(gcmtest$diff),
         colorPal = "PiYG",
         nCol = 21L,
         bg = "gray80")


# Heatmap of differences with given color vector
plotHeat(mat = gcmtest$diff, 
         nCol = 21L,
         color = grDevices::colorRampPalette(c("blue", "white", "orange"))(31))