STM1001: Computer Lab 8B Solutions

Science/Health/Data Science Streams

Topic 8B: Big Data II (\(p\)-value adjustments)

Example R code solutions for the Science/Health/Data Science Computer Lab 8, which uses data from Series GSE184932 on the NCBI Gene Expression Omnibus website are presented below.

This computer lab is designed to run alongside the content in the Foundational Biology for Analyses of Biological Data supplement. It might be helpful to have this material open as you look through these solutions.

1 Introduction

1.1 Type I error probability example

No answer required.

2 Adjusting \(p\)-values

2.1

No answer required.

2.2

Example R code is provided below.

base_p_values <- c(0.0003, 0.0085, 0.001, 0.0001, 0.045, 0.62, 0.009, 0.18, 0.92, 0.02)

2.3

There are 7 \(p\)-values below 0.05, so we would find 7 genes to be significantly differentially expressed.

2.4

p.adjust(base_p_values, method = "bonferroni")

##  [1] 0.003 0.085 0.010 0.001 0.450 1.000 0.090 1.000 1.000 0.200

2.5

Using the Bonferroni correction, we find that only 3 of the \(p\)-values are now below 0.05.

2.6

Example R code is provided below.

# Check all p-values
p.adjust(base_p_values, method = "holm")

##  [1] 0.0027 0.0595 0.0080 0.0010 0.1800 1.0000 0.0595 0.5400 1.0000 0.1000

p.adjust(base_p_values, method = "hochberg")

##  [1] 0.0027 0.0540 0.0080 0.0010 0.1800 0.9200 0.0540 0.5400 0.9200 0.1000

p.adjust(base_p_values, method = "hommel")

##  [1] 0.0027 0.0510 0.0080 0.0010 0.1800 0.9200 0.0540 0.5400 0.9200 0.1000

# Compute number of significant p-values
c(sum(p.adjust(base_p_values, method = "holm") < .05), 
  sum(p.adjust(base_p_values, method = "hochberg") < .05),
  sum(p.adjust(base_p_values, method = "hommel") < .05)
)

## [1] 3 3 3

All three methods result in 3 \(p\)-values below 0.05. The Hochberg and Hommel methods produce near-identical results.

2.7

Example R code is provided below.

# Compute number of significant p-values
c(sum(p.adjust(base_p_values, method = "holm") < .01), 
  sum(p.adjust(base_p_values, method = "hochberg") < .01),
  sum(p.adjust(base_p_values, method = "hommel") < .01)
)

## [1] 3 3 3

All three methods result in 3 \(p\)-values below 0.01. Note that the \(p\)-values themselves have not changed from 2.6, we are simply using a different threshold (level of significance) for determining if results are statistically significant or not.

3 False Discovery Rate

3.1

Example R code is provided below.

p.adjust(base_p_values, method = "fdr")

##  [1] 0.001500000 0.018000000 0.003333333 0.001000000 0.064285714 0.688888889
##  [7] 0.018000000 0.225000000 0.920000000 0.033333333

3.2

We see that there are 6 \(p\)-values below 0.05. The FDR method has been less stringent than the FWER methods, but has still reduced the number of genes considered to be statistically significant.

4 SARS-CoV-2 Study \(p\)-values

4.1

No answer required

4.2

For this question, you must have downloaded the SARS_data_subset.csv file from LMS, and either:

• stored it in your working directory (RStudio users) or
• loaded it into jamovi (jamovi users)

Example R code for the subsequent step is provided below.

# For RStudio, you would have run this code:
sars_data <- read.csv(file = "SARS_data_subset.csv", header = T)

# For jamovi, you would have run this code:
sars_data <- data

4.3

# Store base p-values
covid_base_p_values <- sars_data$pvalue.Alpha_8h_vs_VIC_8h

4.4

# Number of significant p-values
sum(covid_base_p_values < .05)

## [1] 703

So based on the initial \(p\)-values, there are 703 differentially expressed genes at the 8 hour post-infection time point, between the two strains.

4.5

Example R code is provided below.

covid_bonferroni_p_values <- p.adjust(covid_base_p_values, method = "bonferroni")

4.6

Example R code is provided below.

covid_fdr_p_values <- p.adjust(covid_base_p_values, method = "fdr")

4.7

sum(covid_bonferroni_p_values < 0.05)

## [1] 19

sum(covid_fdr_p_values < 0.05)

## [1] 48

Initially, we had 703 statistically significantly differentially expressed genes at the 8 hour post-infection time point, between the two strains.

We can see that following Bonferroni correction we have 19 statistically significantly differentially expressed genes, a reduction of 684 genes.

If we instead use FDR correction, we obtain 48 statistically significantly differentially expressed genes. This is a reduction of 655 genes compared to the initial un-adjusted \(p\)-values. FDR correction is more lenient than the Bonferroni correction approach, with 29 more genes considered statistically significantly differentially expressed.

4.8 Volcano Plots

No answer required.

4.9

The two volcano plots are shown below.

Note that the y-axis shows the negative of the log of the \(p\)-values, meaning that small \(p\)-values appear high on the y-axis, while large \(p\)-values are close to 0 on the y-axis.

Recall also that on the x-axis, a log-fold change of 0 signifies that a gene’s expression status has not changed between the two groups (here Alpha and VIC), while a gene with a log-fold change which is large in magnitude is likely to be biologically significant.

plot(sars_data$log2FoldChange.Alpha_8h_vs_VIC_8h, 
     -log10(covid_base_p_values), 
     main = "Volcano Plot for SARS-CoV-2 
     VIC and Alpha strains
     8h post-infection, using initial p-values",
     xlab = "log FC", ylab = "-log10(pval)",
     cex = 0.5, col = as.numeric(covid_base_p_values <= 0.05) + 1)

plot(sars_data$log2FoldChange.Alpha_8h_vs_VIC_8h, 
     -log10(covid_fdr_p_values), 
     main = "Volcano Plot for SARS-CoV-2 
     VIC and Alpha strains
     8h post-infection, using FDR-corrected p-values",
     xlab = "log FC", ylab = "-log10(pval)",
     cex = 0.5, col = as.numeric(covid_fdr_p_values <= 0.05) + 1)

Note here how many of the FDR-corrected \(p\)-values are no longer considered statistically significant. However, the genes with low FDR-corrected \(p\)-values (and thus high y-axis scores in the volcano plot) appear to have log-fold changes in gene expression which are non-zero and relatively large in magnitude. This suggests that the genes coloured red in the second volcano plot are both biologically and statistically significant - so the FDR-correction has been beneficial, and has helped filter out some genes which, while having low initial \(p\)-values, were not biologically significant in this context.

For this example, it is worth considering why it might be important to compare gene expression differences between the two coronavirus strains. For example, identifying genes which are statistically significantly differentially expressed between the two strains could help researchers develop treatments which were personalised to the different strains. Such treatments could be more effective than using a generic treatment which may have varying levels of efficacy for different strains. Can you think of any other reasons?

That’s everything covered, well done.

References

These notes have been prepared by Rupert Kuveke and Amanda Shaker. The copyright for the material in these notes resides with the authors named above, with the Department of Mathematical and Physical Sciences and with La Trobe University. Copyright in this work is vested in La Trobe University including all La Trobe University branding and naming. Unless otherwise stated, material within this work is licensed under a Creative Commons Attribution-Non Commercial-Non Derivatives License BY-NC-ND.