STM1001: Computer Lab 7B Solutions

Science/Health Science/Data Science Module

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

Example R code solutions for the Week 7 Science/Health Science/Data Science Computer Lab, 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

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.2

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

2.3

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.4

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

2.5

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.

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 removed one gene from contention.

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

4.1

No answer required

4.2

Example R code is provided below.

sars_data <- read.csv(file = "SARS_data_subset.csv", header = T)

4.3

# First, we select only the values in the 
# `pvalue.Alpha_8h_vs_VIC_8h` column that are less than 0.05
sars_data_8h_base <- sars_data[sars_data$pvalue.Alpha_8h_vs_VIC_8h < 0.05,]
# Then we assess the number of rows remaining - each row represents a gene
nrow(sars_data_8h_base)

## [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.4

sars_data$fdr_adjusted_pvalue_8h <- p.adjust(sars_data$pvalue.Alpha_8h_vs_VIC_8h, method = "fdr")

4.5

Example R code is provided below.

sars_sig_fdr_pvalue_genes_8h <- sars_data[sars_data$fdr_adjusted_pvalue_8h  < 0.05,]
nrow(sars_sig_fdr_pvalue_genes_8h)

## [1] 48

We can see that following FDR correction we now have 48 significant genes, a reduction of 655 genes.

4.6 Volcano Plots

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(sars_data$pvalue.Alpha_8h_vs_VIC_8h), 
     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(sars_data$pvalue.Alpha_8h_vs_VIC_8h <= 0.05) + 1)

plot(sars_data$log2FoldChange.Alpha_8h_vs_VIC_8h, 
     -log10(sars_data$fdr_adjusted_pvalue_8h), 
     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(sars_data$fdr_adjusted_pvalue_8h <= 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. The copyright for the material in these notes resides with the author 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.