Sample-QC.Rmd

<!-- R Markdown Documentation, DO NOT EDIT THE PLAIN MARKDOWN VERSION OF THIS FILE -->
  
  <!-- Copyright 2015 Stanford University All rights reserved. -->
  
  Performing QC on AAA genomes using Google Genomics
================================================
  
  The following example makes use of 5 genomes from the MVP.

Setting Up and Describing the Data
----------------------------------
  
  ```{r echo=FALSE, eval=FALSE}
# This codelab assumes that the current working directory is where the Rmd file resides
setwd("/Users/gmcinnes/GitHub/VA_AAA_BigQuery")
```


```{r comment=NA, message=FALSE, warning=FALSE}
# Setup for BigQuery access
require(bigrquery)
require(xtable)
require(RCurl)
require(dplyr)

project <- "gbsc-gcp-project-mvp"                   # put your projectID here

DisplayAndDispatchQuery <- function(queryUri, replacements=list()) {
  if(grepl("^https.*", queryUri)) {
    querySql <- getURL(queryUri, ssl.verifypeer=FALSE)
  } else {
    querySql <- readChar(queryUri, nchars=1e6)
  }
  for(replacement in names(replacements)) {
    querySql <- gsub(replacement, replacements[[replacement]], querySql, fixed=TRUE)
  }
  cat(querySql)
  query_exec(querySql, project)
}

table_replacement <- list("_THE_TABLE_"="gbsc-gcp-project-mvp:va_aaa_pilot_data.sample_gvcfs",
                          "_THE_EXPANDED_TABLE_"="gbsc-gcp-project-mvp:va_aaa_pilot_data.sample_vcfs")
```

Check Singletons
-----------------------------------

```{r message=FALSE, warning=FALSE, comment=NA}
limits = "WHERE 
          reference_name = 'chr22'
          AND call.QUAL >= 30"
result <- DisplayAndDispatchQuery("./sql/private-variants-brca1.sql",
                                  replacements=c(table_replacement, "_LIMITS_" = limits))
```
Number of rows returned by this query: `r 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)
```

Check Individual Heterozygosity
-----------------------------------

```{r message=FALSE, warning=FALSE, comment=NA}
limits = "WHERE 
          reference_name = 'chr22'
          AND call.QUAL >= 30"
result <- DisplayAndDispatchQuery("./sql/homozygous-variants.sql",
                                  replacements=c(table_replacement, "_LIMITS_"=limits))
```
Number of rows returned by this query: `r nrow(result)`.

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

Cohort Level QC
===============

Check Hardy-Weinberg Equilibrium
-----------------------------------
```{r message=FALSE, warning=FALSE, comment=NA}
limits = "WHERE 
          reference_name = 'chr22'
          AND call.QUAL >= 30"
result <- DisplayAndDispatchQuery("./sql/hardy-weinberg-brca1-expanded.sql",
                                  replacements=c(table_replacement, "_LIMITS_"=limits))
```
Number of rows returned by this query: `r nrow(result)`.

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

===============

Check Transition-Transversion Ratio
-----------------------------------
```{r message=FALSE, warning=FALSE, comment=NA}
limits = "AND reference_name = 'chr22'
          AND call.QUAL >= 30"
result <- DisplayAndDispatchQuery("./sql/ti-tv-ratio.sql",
                                  replacements=c(table_replacement,"_LIMITS_"=limits))
```
The result:
```{r echo=FALSE, message=FALSE, warning=FALSE, comment=NA, results="asis"}
print(xtable(result), type="html", include.rownames=F)
```