bg.go()

Description

Searchs the GO terms of the protein contained in a given set

Usage

bg.go(ids)

Arguments

ids either a vector containing the UniProt IDs of the background set or the path to the txt file containing the list of IDs acting as background.

Value

Returns a dataframe with two columns (Uniprot ID, GO terms) and as many rows as different proteins there are in the input set.

References

Rhee, Wood, Dolinkski & Draghici, Nature Reviews Genetics 2008, 9:509–515.

Details

The Gene Ontolory project (GO) provides a controlled vocabulary to describe gene and gene product attributes. Thus, a GO annotation is a statement describing some biological aspect. GO annotations are created by associating a gene or gene product with a GO term. Together, these statements comprise a “snapshot” of current biological knowledge. Hence, GO annotations capture statements about three disjoint categories: cellular component, molecular function and biological process. In this way, GO terms, using a controlled and hierarchical vocabulary, help to describe how a gene functions at the molecular level, where in the cell it functions, and what biological processes it helps to carry out.

The use of this vocabulary (GO annotations) has diverse applications, perhaps the most popular of them is the functional profiling. The goal of functional profiling is to determine which processes might be different in particular sets of genes (or gene products). That is, GO annotations are used to determine which biological processes, functions, and/or locations are significantly over- or under-represented in a given group of genes (or gene products). This approach is also known as GO term enrichment analysis.

Whenever we have a set of genes that are differentially expressed between different conditions (for example, cancerous versus healthy), we can apply GO enrichment analyses. For the sake of concretion let’s build some toy data that may be useful to illustrate the theory behind the GO term enrichment analysis. So, suppose we have analyzed the expression of 20 genes (N = 20) in cancerous and healthy cells. Let’s continue imagining: among these 20 genes, 5 of them are differentially expressed between the two conditions. This set of differentially expressed genes constitutes our target set (T), the remaining 15 genes belong to the complementary set ( $T^c$ ). Please, note that the union of these two sets gives the background set.

At this point, we must be aware that two types of questions can be addressed when performing functional profiling:

Hypothesis-generating query.
Hypothesis-driven query.

In the former (hypothesis-generating query) we assess which GO terms are significant in our target set, while in the later (hypothesis-driven query) we evaluate whether, for instance, apoptosis (or any other preselected process, function or component) is significantly enriched or depleted in our target set. Let’s start with the latter.

Thus, we start computing how many of the 20 genes forming the background set (B), are related to apoptosis and how many of the 5 differentially expressed genes are linked to apoptosis. Suppose, that after performing that computations we obtain the following result:

We can observe that the proportion of genes related to apoptosis within the target set is (3/5 = 0.6) much higher than that in the background set (5/20 = 0.25), but …, what can we conclude? The main problem here is that any enrichment value can occur just by chance. That is, we need an appropriate statistical model. Thus, we define our random variable X as the number of genes in the target set that are related to apoptosis. Now, as a first approach, we can assume that X follows a binomial distribution with parameters p = 7/20 (probability of randomly pick an apoptotic gene from the background set) and n = 5 (cardinal of the target set). That is, $X \sim Bin(p = 0.35, n = 5)$ . Therefore, the probability of getting just by random three apoptotic genes in the target set would be: $P[X = 3] = \binom{5}{3} 0.35^3 0.65^2 = 0.181$ . If we want to estimate the p-value, using the binomial model, we have to estimate the following probability:

1- pbinom(2, size = 5, p = 0.35)

## [1] 0.2351694

$p-value = P[X \geq 3] = 0.235$

A valid critic to the model we have chosen (the binomial one) is that the binomial model assumes sampling with replacement, that is, that the probability of picking from the background set a second gene (or protein) related to apoptosis once we have drawn a first apoptotic gene remain constant. However, such an assumption may not be acceptable, especially when the cardinal of the background set is not very high. For instance, in our toy example, the probability of picking an apoptotic gene is 7/20, but the probability of obtaining a second apoptotic gene is 6/19, and for a third one would be 5/18. Under these circumstances, a better suited model is the hypergeometric distribution $(X \sim HG(N, K, n))$ that describes the number of success (number of apoptotic genes in the target set) in a sequence of $n$ draws (cardinal of the target set) from a finite population (the background set of cardinal N) without replacement. Now, using this HG model we are going to compute again the probability of getting 3 apoptotic genes in our target set just by chance, to this end we will use the probability mass function of the hypergeometric distribution. $P[X = 3] = \frac{\binom{7}{3} \binom{13}{2}}{\binom{20}{5}}$

dhyper(3, 7, 13, 5)

## [1] 0.1760836

Again, to compute the p-value we have to resort to the cumulative distribution function instead to the probability mass function. That is, $P[X \geq 3] = \sum_{i=3}^5 \frac{\binom{7}{i} \binom{13}{5-i}}{\binom{20}{5}}$

1- phyper(2, 7, 13, 5)

## [1] 0.2067853

In general, if we have summarized the results in a contingency table as the one shown in the figure above, the question we want to answer is: knowing that we have $a + b$ genes related to the process (function or component) of interest in the background set, which is formed by $N$ genes, and assuming the null hypothesis (any observed enrichment is due to chance), what is the probability that the these genes of interest would be so unevenly (or even more) distributed as in our observation? In other words, what is the probability of observing data as extreme or more extreme than our observation? To answer this question we use, as we have already noted, the hypergeometric distribution, although in this context we may refer to the statistic model as the Fisher’s exact test. In any case, the reasoning is as follows: In a population of $N$ genes, we have $a + b$ belonging to the process of interest. Now, if we randomly draw (that is the null hypothesis) $a + c$ genes from the background set, what is the probability that $a$ of them would belong to the process of interest?

$P[X = a] = \frac{\binom{a+b}{a} \binom{c+d}{c}}{\binom{N}{a+c}} = \frac{(a + b)! (c + d)! (a + c)! (b + d)!}{a! b! c! d! N!}$

To obtain the p-value, the probability we have to compute is $P[X \geq a]$ , which can be expressed as $1 - F_X(a -1)$ , where $F_X$ is the cumulative distribution function of the random variable X. Remember that $F_X(x) = P[X \leq x]$ . In R:

#  p-value = 1 - phyper(a-1, a+b, c+d, a+c)

If instead of an enrichment test we would like to perform a depletion test, using the same contingency table we should compute: $P[X \leq a] = F_X(a)$ :

#  p-value = phyper(a, a+b, c+d, a+c)

Hitherto, we have addressed the so-called hypothesis-driven approach. Now, we will present the hypothesis-generating one. In this approach we wonder what GO terms are significantly enriched in our test set. To this end, instead of performing a single test as the one we have described to test for enrichment of apoptotic genes, we carry out multiple tests for different categories of genes. Whenever using this second approach, it is important to correct for the fact that many test performed in parallel will greatly increase the number of false positives in the target set. For instance, in the background set shown in the figure above, suppose that the 20 genes belong to 5 disjoint GO categories. If we pick randomly 5 genes from this background set, the probability of the 5 being from the GO category ‘red’ is

dhyper(5, 7, 13, 5)

this probability is as low as 0.001. However, if we compute the probability of the 5 genes belonging to the same category, we find that this probability has been multiplied by the number of tests we have to performed (as many as colors are being tested).

Thus, it is obvious that we have to make a correction for multiple tests. The simplest correction would multiply the p-values of all terms with the number of parallel test performed. This correction is known as the Bonferroni correction. Bonferroni correction is very stringent, that is, the Bonferroni correction can be too conservative in the sense that while it reduces the number of false positives, it also reduces the number of true discoveries. The False Discovery Rate (FDR) approach is a more recent development. This approach also determines adjusted p-values for each test. However, it controls the number of false discoveries in those tests that result in a discovery (i.e. a significant result). Because of this, it is less conservative that the Bonferroni approach and has greater ability (i.e. power) to find truly significant results. The FDR threshold is often called the q-value. Another way to look at the difference between p-values and q-values is that a p-value of 0.05 implies that 5% of all tests will result in false positives, while q-value of 0.05 implies that 5% of significant tests will result in false positives.

How do we compute the q-value? To illustrate its calculation, suppose we have carried out 55 tests for as many GO categories, and the results obtained, ranked from lowest to highest p-value, were:

We are going to add a new column to this table to shown the “adjusted p-value”, which is computed as the p-value times the number of test divided by the rank:

We have emphasized with red the adjusted p-value lower than our significance level $alpha = 0.05$ with the highest rank number. Next, we assign this value to all the rows with a rank number lower than that:

Q-values can be compute from p-values using the qvalue() function from the qvalue package:

# if (!requireNamespace("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")
# 
# BiocManager::install("qvalue")
# 
# library(qvalue)
# p <- c(0.001, 0.002, 0.003, 0.0035, 0.005) # p-values vector
# qvalue(p = p) # computes the corresponding q-values

After this rather long introduction to GO term enrichment analysis, let’s point out that the package ptm offers several functions that will assist you in the process of GO analysis:

search.go
term.go
get.go
bg.go (current doc)
hdfisher.go
gorilla
net.go

Let’s focus on the function bg.go(). This function returns a dataframe with the GO terms associated to the proteins of interest. For instance, if we wish to use as background the set of all the proteins whose modification of any of their methionine residues has been reported to have an effect on the protein activity, we start getting the Uniprot ID of the proteins forming such a set:

sites <- meto.search(highthroughput.group = FALSE,
                 bodyguard.group = FALSE,
                 regulatory.group = TRUE)
ids <- unique(sites$prot_id)

Now, we are ready to use bg.go():

bs <- bg.go(ids)
kable(head(bs))

up_id	GO_id
A0A0E1M553	GO:0003700
Q6PHZ2	GO:0043194, GO:0005954, GO:0005737, GO:0005829, GO:0014704, GO:0031594, GO:0043005, GO:0043025, GO:0005634, GO:0048471, GO:0032991, GO:0016529, GO:0033017, GO:0030315, GO:0005524, GO:0005516, GO:0004683, GO:0042802, GO:0044325, GO:0050998, GO:0042803, GO:0004672, GO:0004674, GO:0019871, GO:0031432, GO:0006816, GO:0086003, GO:0060048, GO:0061049, GO:0030007, GO:0032469, GO:0000082, GO:1902306, GO:2000650, GO:0018105, GO:0018107, GO:0010666, GO:0010613, GO:2000573, GO:0070374, GO:0010971, GO:0035022, GO:0014911, GO:1904707, GO:0046777, GO:0006468, GO:1902514, GO:0098901, GO:0010649, GO:0060341, GO:0003254, GO:1903076, GO:1901897, GO:0010880, GO:0002028, GO:0055119, GO:0001666
P11035	GO:0005829, GO:0005634, GO:0005886, GO:0005773, GO:0071949, GO:0020037, GO:0030151, GO:0043546, GO:0009703, GO:0050464, GO:0008940, GO:0042128, GO:0006809, GO:0009635, GO:0009416, GO:0009610
P52901	GO:0005829, GO:0005967, GO:0005739, GO:0005634, GO:0004739, GO:0006086, GO:0006096, GO:0046686
P08510	GO:0016021, GO:0008076, GO:0005251, GO:0022843, GO:0005249, GO:0001508, GO:0048675, GO:0048150, GO:1903351, GO:0007619, GO:0009584, GO:0008345, GO:0007611, GO:0048047, GO:0045938, GO:0045838, GO:0071805, GO:0006813, GO:0007637, GO:0051260, GO:0045187, GO:0034765, GO:0060025, GO:0050909, GO:0030431
Q12791	GO:0016324, GO:0005901, GO:0016021, GO:0005886, GO:0045211, GO:0008076, GO:0003779, GO:0015269, GO:0042802, GO:0060072, GO:0046872, GO:0005249, GO:0030007, GO:0060073, GO:0045794, GO:0043065, GO:0006813, GO:0034765, GO:0042391, GO:0060087, GO:0051592, GO:0034465, GO:0001666, GO:0006970, GO:0060083

Thus the column bs$GO_id will form the background GO terms set. Of course, we can explore further any of the proteins forming the background set. For instance, let’s put the focus on the protein in the second row. We can get its name and a vector containing its GO terms:

id.features(bs$up_id[2], features = "protein_names")$Protein_names

## [1] "Calcium/calmodulin-dependent protein kinase type II subunit delta (CaM kinase II subunit delta) (CaMK-II subunit delta) (EC 2.7.11.17)"

strsplit(bs$GO_id[2], split = ",")[[1]]

##  [1] "GO:0043194"  " GO:0005954" " GO:0005737" " GO:0005829" " GO:0014704" " GO:0031594"
##  [7] " GO:0043005" " GO:0043025" " GO:0005634" " GO:0048471" " GO:0032991" " GO:0016529"
## [13] " GO:0033017" " GO:0030315" " GO:0005524" " GO:0005516" " GO:0004683" " GO:0042802"
## [19] " GO:0044325" " GO:0050998" " GO:0042803" " GO:0004672" " GO:0004674" " GO:0019871"
## [25] " GO:0031432" " GO:0006816" " GO:0086003" " GO:0060048" " GO:0061049" " GO:0030007"
## [31] " GO:0032469" " GO:0000082" " GO:1902306" " GO:2000650" " GO:0018105" " GO:0018107"
## [37] " GO:0010666" " GO:0010613" " GO:2000573" " GO:0070374" " GO:0010971" " GO:0035022"
## [43] " GO:0014911" " GO:1904707" " GO:0046777" " GO:0006468" " GO:1902514" " GO:0098901"
## [49] " GO:0010649" " GO:0060341" " GO:0003254" " GO:1903076" " GO:1901897" " GO:0010880"
## [55] " GO:0002028" " GO:0055119" " GO:0001666"