Variant-Level-QC.Rmd

<!-- R Markdown Documentation, DO NOT EDIT THE PLAIN MARKDOWN VERSION OF THIS FILE -->

<!-- Copyright 2015 Google Inc. All rights reserved. -->

<!-- Licensed under the Apache License, Version 2.0 (the "License"); -->
<!-- you may not use this file except in compliance with the License. -->
<!-- You may obtain a copy of the License at -->

<!--     http://www.apache.org/licenses/LICENSE-2.0 -->

<!-- Unless required by applicable law or agreed to in writing, software -->
<!-- distributed under the License is distributed on an "AS IS" BASIS, -->
<!-- WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -->
<!-- See the License for the specific language governing permissions and -->
<!-- limitations under the License. -->

*This codelab was made in collaboration with [Google Genomics](https://github.com/googlegenomics).
New [Standard SQL](https://cloud.google.com/bigquery/docs/reference/standard-sql/)
versions of many of these queries can be found
[here](https://github.com/googlegenomics/codelabs/tree/master/R/PlatinumGenomes-QC/sql).*

# Part 3: Variant-Level QC

```{r echo=FALSE, eval=FALSE}
######################[ CHANGE ME ]##################################
# This codelab assumes that the current working directory is where the Rmd file resides.
setwd("/Users/gmcinnes/GitHub/mvp_aaa_codelabs")

# Set the Google Cloud Platform project id under which these queries will run.
project <- "gbsc-gcp-project-mvp"
#####################################################################
```

```{r echo=FALSE, eval=TRUE, message=FALSE, warning=FALSE}
# Set up for BigQuery access.
source("./rHelpers/setup.R")

require(scales)
```

In Part 3 of the codelab, we perform some quality control analyses that could help to identify any problematic variants which should be excluded from further analysis.  The appropriate cut off thresholds will depend upon the input dataset and/or other factors.

* [Setup](#setup)
* [Ti/Tv by Genomic Window](#titv-by-genomic-window)
* [Ti/Tv by Alternate Allele Counts](#titv-by-alternate-allele-counts)
* [Ti/Tv by Depth](#titv-by-depth)
* [Missingness Rate](#missingness-rate)
* [Hardy-Weinberg Equilibrium](#hardy-weinberg-equilibrium)
* [Heterozygous Haplotype](#heterozygous-haplotype)
* [Blacklisted Variants](#blacklisted-variants)


## Setup

```{r}
queryReplacements <- list("_THE_TABLE_"="va_aaa_pilot_data.all_genomes_gvcfs",
                          "_THE_EXPANDED_TABLE_"="va_aaa_pilot_data.all_genomes_expanded_vcfs_java2",
                          "_BLACKLISTED_TABLE_"="resources.blacklisted_positions",
                          "_PATIENT_INFO_"="va_aaa_pilot_data.patient_info")
sampleData <- read.csv("./data/patient_info.csv")
sampleInfo <- select(sampleData, call_call_set_name=Catalog.ID, gender=Gender)
```

ggplot2 themes to use throughout this page.
```{r}
plot_theme = theme_minimal(base_size = 14, base_family = "Helvetica") + 
  theme(axis.line = element_line(colour = "black"),
        panel.grid = element_blank())

boxPlotTheme = theme_minimal(base_size=14, base_family = "Helvetica") +
  theme(panel.grid = element_blank())
```


## Ti/Tv by Genomic Window

We want to check whether either variants are occuring more frequently than random chance at regions across the genome.  One way we can check this is to check whether the ratio of transitions to transversions is within an expected range.  In this analysis we group variants together within 100kb windows accross the entire genome and calculate the ti/tv ratio for each region.  Variants within windows that have a ti/tv ratio outside the expected range will be flagged.

```{r message=FALSE, warning=FALSE, comment=NA}
result <- DisplayAndDispatchQuery("./sql/ti-tv-ratio.sql",
                                  project=project,
                                  replacements=c("#_WHERE_"="WHERE reference_name = 'chr1'",
                                                 "_WINDOW_SIZE_"="100000",
                                                 queryReplacements))
```
Number of rows returned by this query: **`r if(is.null(result)) { "None" } else { nrow(result) }`**.

Displaying the first few rows of the dataframe of results:
```{r echo=FALSE, message=FALSE, warning=FALSE, comment=NA, results="asis"}
if(is.null(result)) { 
  cat("**None**") 
} else {
  print(xtable(head(result)), type="html", include.rownames=F)
}
```

Visualizing the results:
```{r titvByWindow, fig.align="center", fig.width=10, message=FALSE, comment=NA}
ggplot(result, aes(x=window_start, y=titv)) +
  geom_point() +
  stat_smooth() +
  scale_x_continuous(labels=comma) +
  xlab("Genomic Position") +
  ylab("Ti/Tv") +
  ggtitle("Ti/Tv by 100,000 base pair windows on Chromosome 1") +
  plot_theme
```

## Ti/Tv by Alternate Allele Counts

Check whether the ratio of transitions vs. transversions in SNPs appears to be resonable across the range of rare variants to common variants.  This query may help to identify problems with rare or common variants.  

```{r message=FALSE, warning=FALSE, comment=NA}
result <- DisplayAndDispatchQuery("./sql/ti-tv-by-alternate-allele-count.sql",
                                  project=project,
                                  replacements=c(queryReplacements))
```
Number of rows returned by this query: **`r if(is.null(result)) { "None" } else { nrow(result) }`**.

Displaying the first few rows of the dataframe of results:
```{r echo=FALSE, message=FALSE, warning=FALSE, comment=NA, results="asis"}
if(is.null(result)) { 
  cat("**None**") 
} else {
  print(xtable(head(result)), type="html", include.rownames=F)
}
```

Visualizing the results:
```{r titvByAlt, fig.align="center", fig.width=10, message=FALSE, comment=NA}
ggplot(result, aes(x=alternate_allele_count, y=titv)) +
  geom_point() +
  stat_smooth() +
  scale_x_continuous() +
  xlab("Total Number of Sample Alleles with the Variant") +
  ylab("Ti/Tv") +
  ggtitle("Ti/Tv by Alternate Allele Count") +
  plot_theme
```

## Ti/Tv by Depth

We next want to determine whether mutations at a specific depth seem to be occuring at various coverage depths are occuring more often than would be expected by random chance.  We do this by grouping variants from each genome together by coverage depth, e.g. if variant A and variant B both occur in the same sample with 30x coverage they will be grouped for this analysis.  Variants at coverage depths that have ti/tv ratios outside the expected range will be flagged.



```{r message=FALSE, warning=FALSE, comment=NA}
query <- "./sql/ti-tv-by-depth.sql"
result <- DisplayAndDispatchQuery(query,
                                  project=project,
                                  replacements=c(tableReplacement))
```

```{r echo=FALSE, message=FALSE, warning=FALSE, comment=NA, results="asis"}
print(xtable(head(result)), type="html", include.rownames=F)
```

```{r titv-by-depth, fig.align="center", fig.width=10, message=FALSE, comment=NA, warning=FALSE}
ggplot(result, aes(x=average_depth, y=titv_ratio, color=call_call_set_name)) + 
  geom_point() +
  ggtitle("Ti/Tv Ratio By Depth") +
  xlab("Coverage Depth") + 
  ylab("Ti/Tv") + 
  plot_theme +
  theme(legend.position="none")
```

## Missingness Rate

We next want to check how frequenctly or infrequenctly a variant is called across all samples.  If a variant has a very low call rate this may mean it is difficult to sequence and even high confidence calls at that position may be suspect.  

```{r message=FALSE, warning=FALSE, comment=NA}
sortAndLimit <- "ORDER BY missingness_rate DESC, reference_name, start, reference_bases, alternate_bases LIMIT 1000"
cutoff <- list("_CUTOFF_"="0.9")
result <- DisplayAndDispatchQuery("./sql/variant-level-missingness-fail.sql",
                                  project=project,
                                  replacements=c("#_ORDER_BY_"=sortAndLimit,
                                                 queryReplacements,
                                                 cutoff))
```
Number of rows returned by this query: **`r if(is.null(result)) { "None" } else { nrow(result) }`**.

Displaying the first few rows of the dataframe of results:
```{r echo=FALSE, message=FALSE, warning=FALSE, comment=NA, results="asis"}
if(is.null(result)) { 
  cat("**None**") 
} else {
  print(xtable(head(result)), type="html", include.rownames=F)
}
```

## Hardy-Weinberg Equilibrium

For each variant, compute the expected versus observed relationship between allele frequencies and genotype frequencies per the Hardy-Weinberg Equilibrium.

```{r message=FALSE, warning=FALSE, comment=NA}
sortAndLimit <- "ORDER BY ChiSq DESC, reference_name, start, alternate_bases LIMIT 1000"
result <- DisplayAndDispatchQuery("./sql/hardy-weinberg.sql",
                                  project=project,
                                  replacements=c("#_ORDER_BY_"=sortAndLimit,
                                                 queryReplacements))
```
Number of rows returned by this query: **`r if(is.null(result)) { "None" } else { nrow(result) }`**.

Displaying the first few rows of the dataframe of results:
```{r echo=FALSE, message=FALSE, warning=FALSE, comment=NA, results="asis"}
if(is.null(result)) { 
  cat("**None**") 
} else {
  print(xtable(head(result)), type="html", include.rownames=F)
}
```

## Heterozygous Haplotype
For each variant within the X and Y chromosome, identify heterozygous variants in male genomes.  Males should not have any heterozygous calls on the X-chromosome outside of the pseudo-autosomal regions.  If any are found they will be flagged. 

First we use our sample information to determine which genomes are male.  
```{r message=FALSE, warning=FALSE, comment=NA}
maleSampleIds <- paste("'", filter(sampleInfo, gender == "Male")$call_call_set_name, "'", sep="", collapse=",")
```

```{r message=FALSE, warning=FALSE, comment=NA}
sortAndLimit <- "ORDER BY reference_name, start, alternate_bases, call.call_set_name LIMIT 1000"
result <- DisplayAndDispatchQuery("./sql/sex-chromosome-heterozygous-haplotypes.sql",
                                  project=project,
                                  replacements=c("_MALE_SAMPLE_IDS_"=maleSampleIds,
                                                 "#_ORDER_BY_"=sortAndLimit,
                                                 queryReplacements))
```
Number of rows returned by this query: **`r if(is.null(result)) { "None" } else { nrow(result) }`**.

Displaying the first few rows of the dataframe of results:
```{r echo=FALSE, message=FALSE, warning=FALSE, comment=NA, results="asis"}
if(is.null(result)) { 
  cat("**None**") 
} else {
  print(xtable(head(result)), type="html", include.rownames=F)
}
```

## Blacklisted Variants
Identify all variants within our cohort that have been blacklisted.  For more information on what variants are blacklisted and why see [here](https://sites.google.com/site/anshulkundaje/projects/blacklists).  All variants lccuring within blacklisted regions will be flagged.

```{r message=FALSE, warning=FALSE, comment=NA}
result <- DisplayAndDispatchQuery("./sql/blacklisted-variants.sql",
                                  project=project,
                                  replacements=c(queryReplacements))
```
Number of rows returned by this query: **`r if(is.null(result)) { "None" } else { nrow(result) }`**.

Displaying the first few rows of the dataframe of results:
```{r echo=FALSE, message=FALSE, warning=FALSE, comment=NA, results="asis"}
print(xtable(head(result)), type="html", include.rownames=F)
```

--------------------------------------------------------
_Next_: [ Part 4: QC Implementation](./QC-Implementation.md)