
---
title: "Assignment 5"
---

This week's assignment will involve summarizing values for groups of
observations in the METABRIC data set and combining clinical and mutation data
to perform exploratory analysis involving both.

The mutation data come from targeted sequencing of 2433 primary tumours from the
METABRIC study using a panel of 173 of the most frequently mutated breast cancer
genes.

[Pereira *et al.*, Nature Communications 7:11479, 2016](https://www.ncbi.nlm.nih.gov/pubmed/27161491)

We'll start by loading the two tables we're going to need. The first is the
clinical data that you should by now be quite familiar with. This table does not
contain the additional mRNA expression columns we've used in previous
assignments.

```{r message = FALSE, warning = FALSE}
library(tidyverse)
metabric <- read_csv("metabric_clinical_data.csv")
# convert the Integrative_cluster variable into a categorical variable with
# levels in the correct order
metabric <- mutate(metabric, Integrative_cluster = factor(Integrative_cluster, levels = c("1", "2", "3", "4ER-", "4ER+", "5", "6", "7", "8", "9", "10")))
```

The second table contains the mutations that were detected in each patient
tumour sample.

```{r}
mutations <- read_csv("metabric_mutations.csv")
mutations
```

The METABRIC cohort contains 2509 patients but not all patient tumour samples
were sequenced and for some patient samples that were sequenced, no mutations
were detected. For one of the exercises we need to know which patients were
sequenced, so we'll also read the case list of samples that were sequenced and
add a column indicating this in the clinical data table.

```{r}
cases_sequenced <- read_tsv("cases_nat_comm_2016.txt", col_names = "Patient_ID")
metabric <- mutate(metabric, Sequenced = Patient_ID %in% cases_sequenced$Patient_ID)
count(metabric, Sequenced)
```

The clinical data, gene expression values and mutation data were all downloaded
from [cBioPortal](https://www.cbioportal.org/study/summary?id=brca_metabric).

---

**1. Compare the average survival time and other attributes between ER-negative and ER-positive patients**

Compute the average survival time for the ER-negative and ER-positive groups.
Note that such a comparison only makes sense for those patients that are
deceased so apply the appropriate filter first. Add a column for the number
of patients in each group.

```{r}
metabric %>%
  filter(Survival_status == "DECEASED") %>%
  group_by(ER_status) %>%
  summarize(`Average survival time` = mean(Survival_time), N = n())
```

Compute the average tumour size, number of positive lymph nodes and Nottingham
prognostic index within ER-negative and ER-positive patients.

```{r}
metabric %>%
  group_by(ER_status) %>%
  summarize_at(vars(Tumour_size, Lymph_nodes_examined_positive, Nottingham_prognostic_index), mean, na.rm = TRUE)
```

Compute the proportion of tumour samples that are ER positive within each of
PAM50 groups.

```{r}
metabric %>%
  group_by(PAM50) %>%
  summarize(Samples = n(), ER_positive_count = sum(ER_status == "Positive"), ER_positive_proportion = mean(ER_status == "Positive"))
```

---

**2. Find the top 10 most frequently mutated genes**

Count the numbers of mutations for each gene and display the top 10 most
frequently mutated genes.

```{r}
mutations %>%
  count(Gene) %>%
  arrange(desc(n)) %>%
  head(10)
```

Some genes contain multiple mutations within a single tumour sample so we will
be over-counting the numbers of samples containing a mutated gene. Modify your
calculation to prevent this. The count should be the number of patient samples
in which the gene is mutated.

```{r}
mutations %>%
  distinct(Patient_ID, Gene) %>%
  count(Gene) %>%
  arrange(desc(n)) %>%
  head(10)
```

_Hint: use the `distinct()` function_

We're usually only interested in non-silent mutations, i.e. those that may have
some functional effect. One of the findings of the METABRIC study is that 7
genes harbour coding mutations in at least 10% of the samples. A total of 2433
patient samples were sequenced.

Calculate the percentages of samples containing non-silent mutations for each
gene and filter those genes mutated in more than 10% of samples. Compare your
gene list and percentages with those given in the METABRIC 2016 publication
(second paragraph in the Results section).

```{r}
number_of_cases_sequenced <- nrow(cases_sequenced)
mutations %>%
  filter(Variant_Classification != "Silent") %>%
  distinct(Patient_ID, Gene) %>%
  count(Gene) %>%
  mutate(percentage = 100.0 * n / number_of_cases_sequenced) %>%
  filter(percentage >= 10) %>%
  arrange(desc(n))
```

---

**3. Find mutation hotspots in PIK3CA**

PIK3CA and TP53 are the most frequently mutated genes in breast cancers. But do
the mutations occur randomly within these genes or are they usually found in
specific locations?

Find the most frequently mutated codons in PIK3CA.

```{r}
mutations %>%
  filter(Gene == "PIK3CA") %>%
  count(Codon) %>%
  arrange(desc(n))
```

These frequently mutated loci are known as "hotspots".

What type of mutations are found at the most prevalent hotspot in PIK3CA, codon
1047? What changes are occurring at the DNA and protein level?

```{r}
mutations %>%
  filter(Gene == "PIK3CA", Codon == 1047) %>%
  count(Variant_Classification, HGVSc, HGVSp, HGVSp_Short)
```

_Hint: take a look at what columns are available in the `mutations` table and count the values_

---

**4. Compare the prevalance of TP53 mutations across the Integrative Cluster breast cancer subtypes**

Figure 5a from the METABRIC 2016 paper compares the percentage of samples that
have non-silent mutations for various genes within each of the Integrative
Clusters. In this exercise, we'll reproduce a version of one of these bar charts
for the TP53 gene.

First, we'll need to distinguish between those patients that have a non-silent
TP53 mutation and those that don't. Create a data frame containing the number
of non-silent TP53 mutations in each patient sample. It should have two columns,
`Patient_ID` and `TP53_mutation_count`.

```{r}
tp53_mutation_counts <- mutations %>%
  filter(Gene == "TP53", Variant_Classification != "Silent") %>%
  count(Patient_ID, name = "TP53_mutation_count")
```

Now join this table to the `metabric` data frame so that the latter has the
added `TP53_mutation_count` column.

```{r}
metabric <- left_join(metabric, tp53_mutation_counts, by = "Patient_ID")
```

Create a new column called `Has_TP53_mutation` containing just `TRUE` or `FALSE`
values, i.e. a logical variable.

```{r}
metabric <- metabric %>%
  mutate(Has_TP53_mutation = !is.na(TP53_mutation_count))
```

_Hint: use `is.na()` function_

Compute the proportion of patient samples within each Integrative Cluster with a
TP53 mutation. Exclude those samples that have not been sequenced and those that
have not been classified into one of the Integrative Cluster subtypes.

```{r}
tp53_percentages <- metabric %>%
  filter(Sequenced) %>%
  filter(!is.na(Integrative_cluster)) %>%
  group_by(Integrative_cluster) %>%
  summarize(Percentage_TP53_mutation = 100 * mean(Has_TP53_mutation))
tp53_percentages
```

Finally, plot these percentages as a bar chart using `geom_col()`. Look at the
help page for `geom_bar` and `geom_col` to see why we'd want to use the latter
in this case.

```{r}
ggplot(data = tp53_percentages) +
  geom_col(mapping = aes(x = Integrative_cluster, y = Percentage_TP53_mutation), show.legend = FALSE)
```

Customize your plot in the following ways:

* add a title
* change the axis labels so they don't contain underscores
* set the limits on the y axis so that the entire 0 - 100% range is displayed
* set breaks on the y axis to be at 20% intervals
* reduce the amount of space between the bottom of the bars and the tick marks for each integrative cluster
* apply the Viridis colour scheme for the bars (`scale_colour_viridis_d`)
* choose the black and white theme by appending `+ theme_bw()` to the plot
* remove the vertical grid lines by appending `+ theme(panel.grid.major.x = element_blank())`

The final two touches make stylistic changes to the theme using customizations
that we'll cover next week.

```{r}
ggplot(data = tp53_percentages) +
  geom_col(mapping = aes(x = Integrative_cluster, y = Percentage_TP53_mutation, fill = Integrative_cluster), show.legend = FALSE) +
  labs(
    title = "Prevalance of TP53 mutations across Integrative Clusters",
    x = "Integrative Cluster",
    y = "% samples"
  ) +
  scale_y_continuous(breaks = c(0, 20, 40, 60, 80), expand = expand_scale(mult = 0), limit = c(0, 100)) +
  scale_fill_viridis_d() +
  theme_bw() +
  theme(panel.grid.major.x = element_blank())
```

Compare the plot with the equivalent plot in figure 5a from the METABRIC
mutation paper. Do the bars look in roughly the correct proportions? Note that
only 9 of the integrative clusters are shown in figure 5a.

---

**5. Compute the tumour suppressor gene (TSG) score for each gene**

The tumour suppressor gene (TSG) score is the proportion of inactivating
mutations for each gene out of the total number of mutations for that gene.

Inactivating mutations are nonsense SNVs, frameshift substitutions and variants
that affected splice sites.

Compute the TSG score for each gene and find the top 10 tumour suppressor genes,
i.e. those with the highest TSG scores.

```{r}
inactivating_mutation_types <- c("Nonsense_Mutation", "Frame_Shift_Del", "Frame_Shift_Ins", "Splice_Region", "Splice_Site")

tsg_scores <- mutations %>%
  mutate(Inactivating = Variant_Classification %in% inactivating_mutation_types) %>%
  group_by(Gene) %>%
  summarize(TSG_score = mean(Inactivating))

tsg_scores %>%
  arrange(desc(TSG_score)) %>%
  head(10)
```

_Hint: use `%in%` to determine if a mutation is an inactivating mutation_

**Extension**

Download Supplementary Data 3 and create a scatter plot comparing the TSG
scores in the `allSamples_TSGScore` column with those you've calculated.

```{r}
library(readxl)
supplementary_table_3 <- read_excel("ncomms11479-s3.xlsx", skip = 1)
supplementary_table_3 %>%
  select(HGNC_symbol, allSamples_TSGScore) %>%
  left_join(tsg_scores, by = c("HGNC_symbol" = "Gene")) %>%
  ggplot() +
  geom_point(mapping = aes(x = allSamples_TSGScore, TSG_score))
```

