STM1001: Computer Lab 7B

Science/Health Science/Data Science Module

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

Welcome to the seventh computer lab for the Science, Health Science and Data Science modules. In this computer lab we will look at methods for adjusting \(p\)-values obtained when conducting multiple inferences.

This computer lab introduces some biology terminology. Therefore, it is recommended that you read over sections 2.1 to 2.3 of the Introduction to Bioinformatics in R supplement. The material in this supplement provides all the background information you will need.

By the end of this lab, you should understand how to adjust \(p\)-values to counteract the effects of conducting multiple tests simultaneously. Let’s get started!

1 Introduction

By now, you should be familiar with the concept of a hypothesis test. You should also feel comfortable conducting a statistical test that involves computing a \(p\)-value to compare a null hypothesis against an alternative hypothesis (a simple example would be an independent samples \(t\)-test to assess the difference in the mean value of a variable between two groups).

In many situations however, we may be interested in performing a hypothesis test on multiple features simultaneously. For example, in biology, it is now common for thousands of features to be tested against a null hypothesis (where the expectation is that a number of features will be statistically significant).

Unfortunately, when we conduct multiple tests, the probability of obtaining a false positive aka Type I error will become inflated, meaning that the original significance level of \(\alpha\) designated for a single test does not hold.

1.1 Example

As a simple example, suppose that we set \(\alpha = 0.05\) for a single test. Therefore, the probability of not observing a Type I error for this test is \(1-\alpha = 0.95\). So far, so good.

However, if we use this level of significance, and carry out 10 tests, then actually the probability of not observing a Type I error in these 10 tests is \((1 - \alpha)^{10} \approx 0.6\). In other words, by carrying out multiple testing, our probability of Type I error has increased from \(5\%\) to \(40\%\)! As the number of tests increases, this error rate will continue to rise.

2 The p.adjust function

Fortunately, there are several methods we can use in R to adjust \(p\)-values obtained from performing multiple tests. We have already come across one such adjustment in the Topic 7 Readings: the Tukey adjustment, which was used to carry out post-hoc tests after carrying out a One-way ANOVA.

In this computer lab, we will learn about other types of adjustments which may be more appropriate in different settings. Most of these adjustments can be applied using the R function p.adjust. Let’s take a look.

2.1

Suppose that we have compared gene expression levels between two groups of individuals (cases and controls) for 10 key genes of interest. As a result, we have obtained a set of 10 \(p\)-values (one for each gene). These \(p\)-values are presented below:

\(0.0003, 0.0085, 0.001, 0.0001, 0.045, 0.62, 0.009, 0.18, 0.92, 0.02.\)

Store these values in a vector called base_p_values.

2.2

If we assessed these \(p\)-values based on an \(\alpha\) value of \(0.05\), how many genes would we find to have significant differential expression between the cases and controls?

2.3

A simple way to adjust our initial \(p\)-values is to multiply each \(p\)-value by the number of \(p\)-values we have (with any resultant values over 1 scaled back to 1). This is known as the Bonferroni correction method.

Run the code below to conduct this adjustment:

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

2.4

Assessing the adjusted \(p\)-values, how many of the genes remain significantly differentially expressed, at the \(\alpha = 0.05\) level of significance?

2.5

As you will have noticed, the Bonferroni correction is quite conservative, and in practice can often lead to rejecting true positives (especially when dealing with larger numbers of tests). This is because the adjustment method focuses on controlling the probability of making one or more Type I errors when conducting multiple tests. This probability is known as the family-wise error rate (FWER).

2.6

The p.adjust function contains options for several other methods which focus on controlling the FWER. Take a look at the help file for the p.adjust function, and try using the other three FWER-control methods.

Compare and contrast your results, using the R code below as a guide:

# We can compute the number of significant results obtained by a method via:
sum(p.adjust(base_p_values, method = "holm") < .05)
# Note here we use the `holm` method

3 False Discovery Rate

A more flexible and more powerful alternative to controlling the FWER is to control the False Discovery Rate (FDR) when conducting multiple tests. When dealing with Big Data, using the FDR is generally preferred over more conservative adjustments such as Bonferonni correction. FDR control is also typically favoured over FWER control in biological studies involving Big Data¹.

The FDR can be defined as the expected proportion of errors (i.e. false discoveries) among the rejected hypotheses². Controlling the FDR helps to reduce false positives, but also helps minimise false negatives.

3.1

We can apply FDR correction to \(p\)-values using the R function p.adjust - just specify the method fdr.

Try this now, on our base_p_values vector of \(p\)-values.

3.2

Compare your FDR results from 3.1 to your Bonferroni correction results in 2.3. What do you notice?

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

To conclude this computer lab, let’s take a look at a real-world application of \(p\)-value adjustment.

For this example, we will look at \(p\)-values computed as part of an Oxford University expression profiling study relating to SARS-CoV-2 (Coronavirus). Information on the study can be found here³.

We won’t go into all the details of the study. Suffice to say that as part of the study, \(p\)-values were computed when comparing gene expression differences over time, for human cells infected with the VIC or B.1.1.7 (Alpha) strain of SARS-CoV-2.

We have collected a sample of the \(p\)-values from this data set, and these are available on the LMS in the file sars_data.csv. This data has been sourced from publicly available data on the NCBI Gene Expression Omnibus website⁴.

4.1

Download the file sars_data.csv now, and store it in your working directory.

4.2

Read the sars_data.csv file into R, and store it in the object sars_data.

Use the head and summary R commands to take a quick look at this data.

Hint: Recall that you can use the function read.csv to load a .csv file into R.

4.3

The \(p\)-values samples we have selected are computed for gene expression comparisons made between cases and controls at 2 hours, 8 hours, and 24 hours post-infection. For example, the column pvalue.Mock_vs_VIC_2h contains \(p\)-values for the differences between the cases and controls at 2 hours post-infection.

To begin, we will focus on the pvalue.Mock_vs_VIC_2h column.

Run the R code below to determine the number of genes that are considered significantly differentially expressed at the \(\alpha = 0.01\) level of significance (at the 2 hour post-infection mark).

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

4.4

Use the p.adjust function to perform FDR correction on the pvalue.Mock_vs_VIC_2h \(p\)-values, and store your results in a new column in the sars_data object. You can use the R code below as a starting point - just fill in the missing ...’s.

sars_data$pvalue_VIC_2h_fdr <- ...

4.5

Determine the number of genes at the 2 hour post-infection mark that are, according to your FDR results, significantly differentially expressed at the \(\alpha = 0.01\) level of significance. Compare this with your initial results.

Hint: Use the R code in 4.3 as a guide.

4.6

Repeat the process conducted in 4.3 to 4.5 for the other two time points, and record your findings. Why do you think it is important in this example to be able to determine gene expression differences at different time points?

Great job, that’s everything for today.

This was a shorter lab, so you might like to use any remaining time to work on your assessments, or take a look at the reading material for the next Topic.

If you have any questions about Assignment 2, please ask for help.

References

Benjamini, Y., and Y. Hochberg. 1995. “Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing.” J. R. Statist. Soc. B 57 (1): 289–300.

Storey, J. D., and R. Tibshirani. 2003. “Statistical Significance for Genomewide Studies.” PNAS 100 (16): 9440–45.

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 Mathematics and Statistics 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.

---

1. See e.g. Storey and Tibshirani (2003).↩︎

2. See Benjamini and Hochberg (1995).↩︎

3. Please note that at the time of writing, this research has not been published as an academic research paper.↩︎

4. Series GSE184932.↩︎