Sample-Level-QC.Rmd

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

<!-- Copyright 2015 Stanford University 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 2: Sample-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/src/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}
# Load packages
require(scales)

# Set up for BigQuery access.
source("./rHelpers/setup.R")
```

In Parts 2 & 3 of the codelab we describe the details of each quality control step implemented on our population of genomes.  In Part 2 we perform some quality control analyses that could help to identify any problematic genomes that should be removed from the cohort before proceeding with further analysis.  The appropriate cut off thresholds will depend upon the input dataset and/or other factors.

To skip the details and jump to the implementation of the quality control methods skip to [Part 4: QC Implementation](./QC-Implementation.md).
* [Setup](#setup)
* [Missingness Rate](#missingness-rate)
* [Singleton Rate](#singleton-rate)
* [Inbreeding Coefficient](#inbreeding-coefficient)
* [Heterozygosity Rate](#heterozygosity-rate)
* [Sex Inference](#sex-inference)
* [Genotyping Concordance](#genotyping-concordance)
* [Batch Effect](#batch-effect)
* [Ethnicity Inference](#ethnicity-inference)
* [Genome Similarity](#genome-similarity)

## Setup

```{r}
tableReplacement <- list("_THE_TABLE_"="va_aaa_pilot_data.genome_calls_seq_qc",
                          "_THE_EXPANDED_TABLE_"="va_aaa_pilot_data.multi_sample_variants_seq_qc",
                          "_GENOTYPING_TABLE_"="va_aaa_pilot_data.genotyping_data")

ibs <- read.table("./data/all-genomes-ibs.tsv",
                  col.names=c("sample1", "sample2", "ibsScore", "similar", "observed"))

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())
```



## Missingness Rate

Missingess is defined as the proportion of sites found in the reference genome that are not called in a given genome. We calculate the missingness rate of each genome in our cohort in order to identify samples that are potentially low quality.  If a sample has a high missingness rate it may be indicative of issues with sample preparation or sequencing.  Genomes with a missingness rate greater than 0.1 are removed from the cohort.

```{r message=FALSE, warning=FALSE, comment=NA}
result <- DisplayAndDispatchQuery("./sql/missingness-sample-level.sql",
                                  project=project,
                                  replacements=tableReplacement)
```
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)
```

And visualizing the results:
```{r sampleMissingness, fig.align="center", fig.width=10, message=FALSE, comment=NA}
ggplot(result) +
  geom_point(aes(x=sample_id, y=missingness)) +
  xlab("Sample") +
  ylab("Missingness Rate") +
  ggtitle("Genome-Specific Missingness") +
  plot_theme +
  theme(axis.ticks = element_blank(),
        axis.text.x = element_blank())
```

## Singleton Rate

Singleton rate is defined as the number of variants that are unique to a genome.  If a variant is found in only one genome in the cohort it is considered a singleton.  Genomes with singleton rates more than 3 standard deviations away from the mean are removed from the cohort.  

```{r message=FALSE, warning=FALSE, comment=NA}
result <- DisplayAndDispatchQuery("./sql/private-variants.sql",
                                  project=project,
                                  replacements=tableReplacement)
```
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)
```

And visualizing the results:
```{r singletons, fig.align="center", fig.width=10, message=FALSE, comment=NA}
ggplot(result) +
  geom_point(aes(x=call_call_set_name, y=private_variant_count)) +
  xlab("Sample") +
  ylab("Number of Singletons") +
  ggtitle("Count of Singletons Per Genome") +
  scale_y_continuous(labels=comma) +
  plot_theme +
  theme(axis.ticks = element_blank(),
        axis.text.x = element_blank())
```

## Inbreeding Coefficient

The inbreeding coefficient (F) is a measure of expected homozygosity rates vs observed homozygosity rates for individual genomes.  Here, we calculate the inbreeding coefficient using the method-of-moments estimator.  Genomes with an inbreeding coefficient more than 3 standard deviations away from the mean are removed from the cohort.  


```{r message=FALSE, warning=FALSE, comment=NA}
result <- DisplayAndDispatchQuery("./sql/homozygous-variants.sql",
                                  project=project,
                                  replacements=tableReplacement)
```
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)
```

And visualizing the results:
```{r inbreeding, fig.align="center", fig.width=10, message=FALSE, comment=NA}
limits <- c(min(result$O_HOM, result$E_HOM),
            max(result$O_HOM, result$E_HOM))

ggplot(result) +
  geom_point(aes(x=O_HOM, y=E_HOM, label=call_call_set_name), alpha=1/1.5) +
  geom_abline(color="darkslateblue") +
  scale_x_continuous(limits=limits, labels=comma) + 
  scale_y_continuous(limits=limits, labels=comma) +
  xlab("Observed Homozygous Variants") +
  ylab("Expected Homozygous Variants") +
  ggtitle("Homozygosity") +
  plot_theme
```


## Heterozygosity Rate 

Heterozygosity rate is simply the the number of heterozygous calls in a genome.  Genomes with a heterozygosity rate more than 3 standard deviations away from the mean are removed from the cohort.  

```{r message=FALSE, warning=FALSE, comment=NA}
result <- DisplayAndDispatchQuery("./sql/heterozygous-calls-count.sql",
                                  project=project,
                                  replacements=tableReplacement)
```
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)
```

And visualizing the results:
```{r heterozygosity, fig.align="center", fig.width=10, message=FALSE, comment=NA}
ggplot(result) +
  geom_point(aes(y=O_HET, x=call_call_set_name, label=call_call_set_name), alpha=1/1.5) +
  xlab("Sample") +
  ylab("Observed Heterozygous Call Counts") +
  ggtitle("Hetrozygosity Rates") +
  scale_y_continuous(labels=comma) +
  plot_theme +
  theme(axis.ticks = element_blank(),
        axis.text.x = element_blank())
```



## Sex Inference

Gender is inferred for each genome by calculating the heterozygosity rate on the X chromosome.  Genomes who's inferred sex is different from that of the reported sex are removed from the cohort.  Although it is possible for people to be genotypically male and phenotypically female, it is more likely that samples or phenotypic records were mislabeled.

```{r message=FALSE, warning=FALSE, comment=NA}
result <- DisplayAndDispatchQuery("./sql/gender-check.sql",
                                  project=project,
                                  replacements=tableReplacement)
```
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)
```

Let's join this with the sample information:
```{r message=FALSE, warning=FALSE, comment=NA}
joinedResult <- inner_join(result, sampleInfo)
```

And visualize the results:

```{r gender, fig.align="center", fig.width=10, message=FALSE, comment=NA}
ggplot(joinedResult) +
  geom_point(aes(x=call_call_set_name, y=perct_het_alt_in_snvs, color=gender), size=3) +
  xlab("Sample") +
  ylab("Heterozygosity Rate") +
  ggtitle("Heterozygosity Rate on the X Chromosome") +
  scale_colour_brewer(palette="Set1", name="Gender") +
  plot_theme +
  theme(axis.ticks = element_blank(),
    axis.text.x = element_blank())
```

## Genotyping Concordance

We next want to look at the concordance between SNPs called from the sequencing data and those called through the use genotyping.  This allows us to identify samples that may have been mixed up in the laboratory.  Samples with low concordance (>99%) should be removed from the cohort.

```{r message=FALSE, warning=FALSE, comment=NA}
concordanceResult <- DisplayAndDispatchQuery("./sql/genotyping-concordance.sql",
                                  project=project,
                                  replacements=tableReplacement)
```
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(concordanceResult)), type="html", include.rownames=F)
```

Get the sample preparation plate for each sample
```{r message=FALSE, warning=FALSE, comment=NA}
plate = substr(concordanceResult$sample_id, 1, 9)
concordanceResult = cbind(concordanceResult, plate)
```

Visualizing the results:
```{r concordance, fig.align="center", fig.width=10, message=FALSE, comment=NA}
ggplot(concordanceResult) +
  geom_point(aes(x=sample_id, y=concordance, color=plate), size=3) +
  xlab("Sample") +
  ylab("Concordance") +
  ggtitle("Concordance with Genotyping Data") +
  scale_colour_brewer(name="Sample Prep Plate", palette="Set1") +
  plot_theme +
  theme(axis.ticks = element_blank(),
        axis.text.x = element_blank())
```


## Batch Effect

We next want to determine if there are any batch effects that were introduced during sequencing.  Since the genomes in our study were sequenced over time it is important to check if we can identify sequencing batches by variant data alone.  We can check for batch effects by running principal component analysis on all of the variants from our dataset.  If any of the batches are outliers in this part of the analysis we may want to remove them from the cohort.

Note: This analysis is run on a Google Compute Engine rather than in BigQuery.  For full instructions on how to run PCA see the [Google Genomics Cookbook](https://googlegenomics.readthedocs.org/en/latest/use_cases/compute_principal_coordinate_analysis/1-way-pca.html?highlight=pca).

Create cluster with Apache Spark installed
```
./bdutil -e extensions/spark/spark_env.sh deploy
```

Copy client_secrets to head node
```
gcloud compute copy-files client_secrets.json hadoop-m:~/
```

SSH into the head node
```
gcloud compute ssh hadoop-m
```

Install SBT
```
echo "deb http://dl.bintray.com/sbt/debian /" | sudo tee -a /etc/apt/sources.list.d/sbt.list
sudo apt-get update
sudo apt-get install sbt
```

Install git and clone spark-examples repository.
```
sudo apt-get install git
git clone https://github.com/googlegenomics/spark-examples.git
```

Compile the jar
```
cd spark-examples
sbt assembly
cp target/scala-2.10/googlegenomics-spark-examples-assembly-*.jar ~/
cd ~/
```

Running the job.  

Variant set definitions:
2442214609072563359: Full variant set for AAA dataset.

```
spark-submit \
  --class com.google.cloud.genomics.spark.examples.VariantsPcaDriver \
  --conf spark.shuffle.spill=true \
  --master spark://hadoop-m:7077 \
  /home/gmcinnes/spark-examples/target/scala-2.10/googlegenomics-spark-examples-assembly-1.0.jar \
  --client-secrets /home/gmcinnes/client_secrets.json \
  --variant-set-id 2442214609072563359 \
  --all-references \
  --num-reduce-partitions 15 \
  --output-path gs://gbsc-gcp-project-mvp-va_aaa_hadoop/spark/output/aaa-batch-effect-pca.tsv
```



Next, we need to download the files that were output from the job to our local machine.

```
gsutil cat gs://gbsc-gcp-project-mvp-va_aaa_hadoop/spark/output/aaa-batch-effect-pca.tsv* > aaa-batch-effect-pca.tsv
```

Now let's visualize the results.

```{r message=FALSE, warning=FALSE, comment=NA}
pcaFile = './data/aaa-batch-effect-pca.tsv'
pcaResult = read.table(pcaFile)
names(pcaResult) = c('sample_id','pc1', 'pc2', 'extra')
plate = substr(pcaResult$sample_id, 1, 9)
pcaResult = cbind(pcaResult, plate)
```

```{r pca-batch-effect-publication, fig.align="center", fig.width=10, message=FALSE, comment=NA}
ggplot(pcaResult, aes(pc1, pc2, color=plate)) + 
  geom_point(size=4, alpha=0.5) +
  ggtitle("Batch Effect PCA") +
  xlab("Principal component 1") + 
  ylab("Principal component 2") +
  scale_colour_brewer(name="Library Prep Plate", palette="Set1") +
  plot_theme +
  theme(axis.text.x=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks=element_blank(),
        legend.position = c(0.85,0.75))
```


## Ethnicity Inference

Another application of principle component analysis is inferring ethnicity.  We can do so by performing PCA with our data merged with data from the 1,000 Genomes Project.  When we do this we expect to see that genomes from the 1,000 Genome Project and genomes from our dataset will cluster by ethnicity.  Since almost all of the genomes from our dataset are of European descent the European subpopulations from 1,000 Genomes and our genomes should cluster together.  Any genomes that cluster with an unexpected ethnic group should be removed from the cohort.

Full instructions for how to compute principal component analysis on the intersection of two datasets using Google Cloud can be found [here](http://googlegenomics.readthedocs.org/en/latest/use_cases/compute_principal_coordinate_analysis/2-way-pca.html).

#### Step by step instructions for running 2-way PCA

These instructions are identical to the Batch Effect step except that we run 2-way PCA instead of 1-way PCA.

Create cluster with Apache Spark installed
```
./bdutil -e extensions/spark/spark_env.sh deploy
```

Copy client_secrets to head node
```
gcloud compute copy-files client_secrets.json hadoop-m:~/
```

SSH into the head node
```
gcloud compute ssh hadoop-m
```

Install SBT
```
echo "deb http://dl.bintray.com/sbt/debian /" | sudo tee -a /etc/apt/sources.list.d/sbt.list
sudo apt-get update
sudo apt-get install sbt
```

Install git and clone spark-examples repository.
```
sudo apt-get install git
git clone https://github.com/googlegenomics/spark-examples.git
```

Compile the jar
```
cd spark-examples
sbt assembly
cp target/scala-2.10/googlegenomics-spark-examples-assembly-*.jar ~/
cd ~/
```

#### Run the job

Variant set definitions:
18444522614861607145: SNPs from the AAA dataset to use for ethnicity inference.  
12511568612726359419: 1000 genomes variant set.

```
spark-submit   \
  --class com.google.cloud.genomics.spark.examples.VariantsPcaDriver \
  --conf spark.shuffle.spill=true \
  --master spark://hadoop-m:7077 \
  /home/gmcinnes/spark-examples/target/scala-2.10/googlegenomics-spark-examples-assembly-1.0.jar \
  --client-secrets /home/gmcinnes/client_secrets.json \
  --variant-set-id 18444522614861607145 12511568612726359419 \
  --all-references \
  --num-reduce-partitions 15 \
  --output-path gs://gbsc-gcp-project-mvp-va_aaa_hadoop/spark/output/aaa-vs-1kg-pca.txt
```

Next, we need to download the files that were output from the job to our local machine.

```
gsutil cat gs://gbsc-gcp-project-mvp-va_aaa_hadoop/spark/output/aaa-vs-1kg-pca.txt* > aaa-vs-1kg-pca.tsv
```

Now let's visualize the results.

```{r message=FALSE, warning=FALSE, comment=NA}
pca = read.table('./data/aaa-vs-1kg-pca.tsv')
names(pca) <- c("sample_id","pc1","pc2")
```

```{r}
populations = read.csv('./data/1kg_info.csv')
pca = join(pca, populations, by='sample_id')
```

```{r ethnicity-inference, fig.align="center", fig.width=10, message=FALSE, comment=NA}
ggplot(pca) +
  geom_point(aes(pc1,pc2, color=Population)) +
  ggtitle("Ethnicity inference") +
  scale_colour_brewer(name="1KG super population", palette="Set1") +
  xlab("Principal component 1") +
  ylab("Principal component 2") +
  plot_theme + 
  theme(axis.ticks=element_blank(),
        axis.text=element_blank(),
        legend.position = c(0.85, 0.7))
```

## Genome Similarity

The final step in the Sample QC portion of our workflow is a check for unexpected relatedness among our genomes.  To do this we perfrom an analysis called Identity By State (IBS).  This analysis checks the similarity between all positions in every genome against every other genome.  Any unexpected similarity among our genomes may be the result of unknown family relationships or the result of cross-sample contamination.  Any samples with unexpected relatedness should be removed from the cohort.

Full instructions for setting up and running this job can be found in the [Google Genomics Cookbook](http://googlegenomics.readthedocs.org/en/latest/use_cases/compute_identity_by_state/).

Note that this `n^2` analysis is a cluster compute job instead of a BigQuery query.


#### Run the job
```
java -cp target/google-genomics-dataflow*.jar \
  com.google.cloud.genomics.dataflow.pipelines.IdentityByState \
  --runner=BlockingDataflowPipelineRunner \
  --project=gbsc-gcp-project-mvp \
  --stagingLocation=gs://gbsc-gcp-project-mvp-va_aaa_jobs/dataflow/all-genomes-staging \
  --output=gs://gbsc-gcp-project-mvp-va_aaa_hadoop/dataflow/output/all-genomes-ibs.tsv \
  --numWorkers=40 \
  --genomicsSecretsFile=client_secrets.json \
  --datasetId=2442214609072563359 \
  --hasNonVariantSegments=true \
  —allReferences=true
```




```{r}
require(reshape2)
require(dplyr)
ibsDataFlowFilename = '/Users/gmcinnes/data/all-genomes-ibs-2.tsv'
ReadIBSFile <- function(ibsFilename, header=FALSE, rowNames=NULL) {
  ibsData <- read.table(file=ibsFilename, header=header,
                        row.names=rowNames, stringsAsFactors=FALSE)
  return (ibsData)
}
ibsDataflowData <- ReadIBSFile(ibsDataFlowFilename)

ColumnNames <- function(ibsData) { 
  if(3 == ncol(ibsData)) {
    colnames(ibsData) <- c("sample1", "sample2", "ibsScore")
  } else {
    colnames(ibsData) <- c("sample1", "sample2", "ibsScore", "similar", "observed")
  }
}
colnames(ibsDataflowData) <- ColumnNames(ibsDataflowData)

MakeIBSDataSymmetric <- function(ibsData) {
  ibsPairsMirrored <- data.frame(sample1=ibsData$sample2,
                                 sample2=ibsData$sample1,
                                 ibsScore=ibsData$ibsScore)
  ibsData <- rbind(ibsData[,1:3], ibsPairsMirrored)
}
ibsDataflowData <- MakeIBSDataSymmetric(ibsDataflowData)

ExcludeDiagonal <- function(ibsData) {
  ibsData <- filter(ibsData, ibsData$sample1 != ibsData$sample2)
  return (ibsData)
}
ibsDataflowDataSample <- ExcludeDiagonal(ibsDataflowData)

SampleIBSMatrix <- function(ibsData, sampleSize=50) {
  individuals <- unique(ibsData$sample1)
  sample <- sample(individuals, sampleSize)
  ibsData <- subset(ibsData, ibsData$sample1 %in% sample)
  ibsData <- subset(ibsData, ibsData$sample2 %in% sample)
  return (ibsData)
}
ibsDataflowDataSubset <- SampleIBSMatrix(ibsDataflowDataSample)
```

Let's plot a subset of the data to understand the plot
```{r ibs-1, fig.align="center", fig.width=12, message=FALSE, comment=NA}
DrawHeatMap <- function(ibsData) {
  p <- ggplot(data=ibsData, aes(x=sample1, y=sample2)) +
    theme(axis.ticks=element_blank(), axis.text=element_blank()) +
    geom_tile(aes(fill=ibsScore), colour="white") +
    scale_fill_gradient(low="white", high="steelblue", na.value="black",
                        guide=guide_colourbar(title= "IBS Score")) +
    labs(list(title="Identity By State (IBS) Heat Map",
              x="Sample", y="Sample"))
  p
}
DrawHeatMap(ibsDataflowDataSubset)
```

Let's take a look at the most similar genomes.
```{r echo=FALSE, message=FALSE, warning=FALSE, comment=NA, results="asis"}
related = ibsDataflowDataSample[ibsDataflowDataSample$ibsScore > 0.06,]
print(xtable(head(related)), type="html", include.rownames=F)
```


--------------------------------------------------------
_Next_: [Part 3: Variant-Level QC](./Variant-Level-QC.md)