Comparison.Rmd

# Comparison of Google Genomics tools vs Academic Tools

The purpose of this codelab is to compare BigQuery queries to academic command line tools for genomic analysis. 

* [Setup](#setup)  
* [Singleton Rate](#singleton-rate)
* [Homozygosity Rate and Inbreeding Coefficient](#homozygosity-rate-and-inbreeding-coefficient)
* [Hardy-Weinberg Equilibrium](#hardy-weingberg-equilibrium)
* [Ti/Tv by Alternate Allele Counts](#titv-by-alternate-allele-counts)
* [Ti/Tv by Depth](#titv-by-depth)
* [Sample-Level Missingness](#sample-level-missingness)
* [Variant-Level Missingness](#variant-level-missingness)
* [Sex Inference](#sex-inference)
* [Heterozygous Haplotype](#heterozygous-haplotype)
* [Ethnicity Inference](#ethnicity-inference)
* [Genome Similarity](#genome-similarity)
* [Genotype Concordance](#genotype-concordance)

## Setup
```{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/qc")

# 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")
```


```{r}
tableReplacement <- list("_THE_TABLE_"="gbsc-gcp-project-mvp:va_aaa_pilot_data.5_genome_test_gvcfs_2",
                          "_THE_EXPANDED_TABLE_"="gbsc-gcp-project-mvp:va_aaa_pilot_data.5_genome_test_vcfs_2",
                         "_REF_TABLE_"="gbsc-gcp-project-mvp:qc_tables.5_genomes_ref_calls_brca1",
                         "_VARIANT_TABLE_"="gbsc-gcp-project-mvp:qc_tables.5_genomes_variants_brca1",
                         "_GENOTYPING_TABLE_"="gbsc-gcp-project-mvp:va_aaa_pilot_data.5_genome_test_genotyping_vcfs")
sampleData <- read.csv("./data/patient_info.csv")
sampleInfo <- select(sampleData, call_call_set_name=Catalog.ID, gender=Gender)

constraints <- list("#_WHERE_"="WHERE 
                    reference_name = 'chr17'
                    AND start BETWEEN 41196311
                    AND 41277499")

googleRepository = "https://raw.githubusercontent.com/googlegenomics/codelabs/master/R/PlatinumGenomes-QC/sql/"
```
 

[Top](#top)

## Singleton Rate
[Singletone Rate in Google Genomics](https://github.com/googlegenomics/codelabs/blob/master/R/PlatinumGenomes-QC/Sample-Level-QC.md#singleton-rate)

For each genome, count the number of variants shared by no other member of the cohort.  Too many private calls for a particular individual may indicate a problem.

*For our 5 genome sample set the singleton rate will be very high due to the small number of genomes.*

```{r message=FALSE, warning=FALSE, comment=NA}
query <- paste(googleRepository, "private-variants-brca1.sql", sep = "")
result <- DisplayAndDispatchQuery(query,
                                  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)
```

Compare to [output](./data/brca1/singletons.singletons) from vcftools, which has 85 some of which are for 0/0 genotypes from reference matching blocks (see the [vcftools command line](./data/brca1/singletons.e9636222) used to create this file).

Command line execution
```
vcftools --gzvcf brca1.merged.sample_set.vcf.gz --singletons --out singletons
```

```{r}
expectedResult <- read.table("./data/brca1/singletons.singletons", header=TRUE)
# Convert to zero-based coordinates
expectedResult <- mutate(expectedResult, POS = POS - 1)
# Clean colnames to match
colnames(expectedResult) <- gsub('\\.+', '_', colnames(expectedResult))
```

How many singletons do the two results have in common?
```{r, warning=FALSE}
nrow(inner_join(result, expectedResult))
```

How many singletons were only identified by BigQuery?
```{r, warning=FALSE}
onlyBQ <- anti_join(result, expectedResult)
nrow(onlyBQ)
```

```{r results='asis', echo=FALSE, warning=FALSE}
if(nrow(onlyBQ > 0)) {
  print(xtable(onlyBQ), type="html", include.rownames=F)
}
```

Which singletons were only identified by vcftools?
```{r results="asis", warning=FALSE}
onlyVcftools <- anti_join(expectedResult, result)
print(xtable(onlyVcftools), type="html", include.rownames=F)
```

Retrieving the gVCF data for the singletons identified only by vcftools:
```{r message=FALSE, warning=FALSE, comment=NA}
having <- paste("start = ", onlyVcftools$POS,
                sep="", collapse=" OR ")
query <- paste(googleRepository, "examine-data.sql", sep="")
result <- DisplayAndDispatchQuery(query,
                                  project=project,
                                  replacements=c(tableReplacement,
                                                 "_HAVING_"=having))
```

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

It appears that they correspond either to:

* A reference-matching block, so not actually a singleton and just perhaps violating an assumption in the vcftools code.  Most of the start sites queried return nothing because they are withing reference matching blocks.
* Or a non-singleon variant, perhaps due to a problem in converting the gVCF data to all-positions VCF via gvcftools?

[Top](#top)

## Homozygosity Rate and Inbreeding Coefficient
[Homozygosity Rate and Inbreeding Coefficient in Google Genomics](https://github.com/googlegenomics/codelabs/blob/master/R/PlatinumGenomes-QC/Sample-Level-QC.md#homozygosity-rate-and-inbreeding-coefficient)

For each genome, compare the expected and observed rates of homozygosity.

```{r message=FALSE, warning=FALSE, comment=NA}
query <- paste(googleRepository, "homozygous-variants.sql", sep="") 
result <- DisplayAndDispatchQuery(query,
                                  project=project,
                                  replacements=c(tableReplacement, constraints))
```
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 compare to what was found using vcftools.

Command line execution.

```
vcftools --gzvcf brca1.merged.sample_set.vcf.gz --het --out het_calls
```

```{r}
expectedResult <- read.table("./data/brca1/het_calls.het", header=TRUE)
```

```{r results="asis", echo=FALSE}
# Clean colnames to match
colnames(expectedResult) <- gsub('\\.+$', '', colnames(expectedResult))
colnames(expectedResult) <- gsub('\\.+', '_', colnames(expectedResult))
n = names(result)
n[1] = "INDV"
names(result) = n
joinedResult <- inner_join(expectedResult, result, by=c("INDV"))
print(xtable(joinedResult[,order(colnames(joinedResult))]), type="html", include.rownames=F)
```

The logic in the query looks similar to vcftools [output_het method](http://sourceforge.net/p/vcftools/code/HEAD/tree/trunk/cpp/variant_file_output.cpp#l165) but there is clearly a difference.  

TODO: investigate the difference further.

[Top](#top)

Hardy-Weinberg Equilibrium
-----------------------------------
[Hardy-Weinberg Equiibrium in Google Genomics](https://github.com/googlegenomics/codelabs/blob/master/R/PlatinumGenomes-QC/Variant-Level-QC.md#hardy-weinberg-equilibrium)

```{r message=FALSE, warning=FALSE, comment=NA}
query <- paste(googleRepository, "hardy-weinberg.sql", sep="")
result <- DisplayAndDispatchQuery(query,
                                  project=project,
                                  replacements=c(tableReplacement, constraints))
```
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)
```

Compare to [brca1.hwe](./data/brca1/hardy.hwe) (see the [vcftools command line](./data/hwe/brca1.log) used to create this file).

Command line
```
vcftools --gzvcf brca1.merged.sample_set.vcf.gz --hardy --out hardy
```

```{r}
require(dplyr)
df <- read.table("./data/brca1/hardy.hwe", header=TRUE)
obsSplitCol <- "OBS.HOM1.HET.HOM2."
obsTemp <- read.table(text=as.character(df[, obsSplitCol]), sep = "/")
names(obsTemp) <- c("OBS_HOM1", "OBS_HET", "OBS_HOM2")
eSplitCol <- "E.HOM1.HET.HOM2."
eTemp <- read.table(text=as.character(df[, eSplitCol]), sep = "/")
names(eTemp) <- c("E_HOM1", "E_HET", "E_HOM2")
expectedResult <- cbind(cbind(df[setdiff(names(df), c(obsSplitCol,eSplitCol))], obsTemp), eTemp)
# Convert to zero-based coordinates
expectedResult <- mutate(expectedResult, POS = POS - 1)
names(expectedResult)[1] = 'reference_name'
names(expectedResult)[2] = 'start'
```

How many results do the two results have in common?
```{r}
nrow(inner_join(result, expectedResult, by=c("reference_name", "start", "OBS_HOM1", "OBS_HET", "OBS_HOM2")))
```

How many result were only identified by BigQuery?
```{r}
onlyBQ <- anti_join(result, expectedResult, , by=c("reference_name", "start", "OBS_HOM1", "OBS_HET", "OBS_HOM2"))
nrow(onlyBQ)
```

Let's take a look at the first few.
```{r results="asis", echo=FALSE}
if(nrow(onlyBQ) > 0) {
  print(xtable(head(arrange(onlyBQ, reference_name, start))), type="html", include.rownames=F)
}
```

How many results were only identified by vcftools?
```{r, warning=FALSE}
onlyVcftools <- anti_join(expectedResult, result, , by=c("reference_name", "start", "OBS_HOM1", "OBS_HET", "OBS_HOM2"))
nrow(onlyVcftools)
```

Let's take a look at the first few.
```{r results='asis', echo=FALSE}
print(xtable(head(arrange(onlyVcftools, reference_name, start))), type="html", include.rownames=F)
```

Retrieving the gVCF data for the results identified only by vcftools:
```{r message=FALSE, warning=FALSE, comment=NA}
having <- list("_HAVING_" = paste("start <= ", onlyVcftools$start,
                "AND",
                "end >= ", onlyVcftools$start+1,
                collapse=" OR "))
query <- paste(googleRepository, "examine-data.sql", sep="")
result <- DisplayAndDispatchQuery(query,
                                  project=project,
                                  replacements=c(tableReplacement,
                                                 having))
```

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

[Top](#top)

Ti/Tv by Alternate Allele Counts
-----------------------------------
[Ti/Tv by Alternate Allele Counts in Google Genomics](https://github.com/googlegenomics/codelabs/blob/master/R/PlatinumGenomes-QC/Variant-Level-QC.md#titv-by-alternate-allele-counts)

```{r message=FALSE, warning=FALSE, comment=NA}
query <- paste(googleRepository, "ti-tv-by-alternate-allele-count.sql", sep="")
result <- DisplayAndDispatchQuery(query,
                                  project=project,
                                  replacements=c(tableReplacement, constraints))
```

```{r echo=FALSE, message=FALSE, warning=FALSE, comment=NA, results="asis"}
print(xtable(result), type="html", include.rownames=F)
```
    
Compare to [output](./data/brca1/TsTv-by-count.TsTv.count) from vcftools.

Command line
```
vcftools --gzvcf brca1.merged.sample_set.vcf.gz --TsTv-by-count --out TsTv-by-coun
```

```{r}
expectedResult <- read.table("./data/brca1/TsTv-by-count.TsTv.count", header=TRUE)
```

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

The results for BigQuery and vcftools are identical except for the case where there are 10 alternate alleles.  With only 5 genomes, 10 sites with an alternate is the maximum.  vcftools may ignore sites where all calls are the alternate allele.

[Top](#top)

Ti/Tv by Depth
-----------------------------------
```{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() +
  #geom_smooth(se=FALSE,    # Don't add shaded confidence region
  #              fullrange=T) +
  ggtitle("Ti/Tv Ratio By Depth") +
  xlab("Coverage Depth") + 
  ylab("Ti/Tv")
```


TODO: Find academic tool to compare against.

[Top](#top)

Sample-Level Missingness
-----------------------------------
[Sample-Level Missingness in Google Genomics](https://github.com/googlegenomics/codelabs/blob/master/R/PlatinumGenomes-QC/Sample-Level-QC.md#missingness-rate)

See [this document](./Sample_Missingness_Preparation.md) for details on how tables were created for this section.

```{r message=FALSE, warning=FALSE, comment=NA}
query <- "./sql/sample-missingness-missing-calls.sql"
result <- DisplayAndDispatchQuery(query,
                                  project=project,
                                  replacements=c(tableReplacement, constraints))
```

```{r echo=FALSE, message=FALSE, warning=FALSE, comment=NA, results="asis"}
print(xtable(result), type="html", include.rownames=F)
```
    
Compare to [output](./data/brca1/missingness.imiss) from vcftools.

Command line
```
vcftools --gzvcf brca1.merged.sample_set.vcf.gz --missing-indv --out missingness
```

```{r}
expectedResult <- read.table("./data/brca1/missingness.imiss", header=TRUE)
```

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

These results from BigQuery and vcftools are expected to be different.  Sample-level missingness in BigQuery is calculated by joining all positions in the genome against the reference genome to identify all sites that were not called by sequencing.  Vcftools only looks at sites that were called by at least one sample in the group.

[Top](#top)

Variant-Level Missingness
-----------------------------------
[Variant-Level Missingness in Google Genomics](https://github.com/googlegenomics/codelabs/blob/master/R/PlatinumGenomes-QC/Variant-Level-QC.md#missingness-rate)
```{r message=FALSE, warning=FALSE, comment=NA}
query <- "./sql/variant-level-missingness-missing-calls.sql"
result <- DisplayAndDispatchQuery(query,
                                  project=project,
                                  replacements=c(tableReplacement, constraints))
```

```{r echo=FALSE, message=FALSE, warning=FALSE, comment=NA, results="asis"}
print(xtable(head(result)), type="html", include.rownames=F)
```
    
    
Compare to [output](./data/brca1/missing_sites.lmiss) from vcftools.

Command line
```
vcftools --gzvcf brca1.merged.sample_set.vcf.gz --missing_site --out missing_sites
```

```{r}
expectedResult <- read.table("./data/brca1/missing_sites.lmiss", header=TRUE)
```

```{r echo=FALSE, message=FALSE, warning=FALSE, comment=NA, results="asis"}
```{r}
expectedResult <- mutate(expectedResult, POS = POS - 1)
n = c("reference_name", "start", "expected_count", "filtered", "expected_missing_count", "expected_rate")
names(expectedResult) = n
```

```{r results='asis'}
print(xtable(head(expectedResult)), type="html", include.rownames=F)
```

Let's join them and look at the differences
```{r}
joinedResult = data.table(inner_join(result, expectedResult, by=c("reference_name", "start")))
matches = joinedResult[missingness_rate == expected_rate]
differences = joinedResult[missingness_rate != expected_rate]
```

How many do they have in common?
```{r}
nrow(matches)
```

How many are different?
```{r}
nrow(differences)
```

Here we can see which rows where different
```{r results='asis', echo=FALSE}
print(xtable(differences), type="html", include.rownames=F)
```

The differences appear to be mostly form positions with insertions or deletions.

TODO: Look at the differences more closely.

[Top](#top)

Sex Inference
-----------------------------------
[Sex Inference in Google Genomics](https://github.com/googlegenomics/codelabs/blob/master/R/PlatinumGenomes-QC/Sample-Level-QC.md#sex-inference)

```{r message=FALSE, warning=FALSE, comment=NA}
query <- paste(googleRepository, "gender-check.sql", sep="")
result <- DisplayAndDispatchQuery(query,
                                  project=project,
                                  replacements=c(tableReplacement, constraints))
```

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)) +
  theme(axis.text.x=if(nrow(result) <= 20)
    {element_text(angle = 90, hjust = 1)} else {element_blank()}) +
  xlab("Sample") +
  ylab("Heterozygosity Rate ") +
  ggtitle("Heterozygosity Rate on the X Chromosome")

```

TODO: Compare to Plink

Command line
```
plink --vcf ~/scratch/chrX.vcf.gz --make-bed --split-x 'hg19' --out ~/scratch/chrX-split
plink --bfile ~/scratch/chrX-split --check-sex
```

[Top](#top)

Heterozygous Haplotype
-----------------------------------
[Heterozygous Haplotype in Google Genomics](https://github.com/googlegenomics/codelabs/blob/master/R/PlatinumGenomes-QC/Variant-Level-QC.md#heterozygous-haplotype)
```{r message=FALSE, warning=FALSE, comment=NA}
query <- paste(googleRepository, "sex-chromosome-heterozygous-haplotypes.sql", sep="")
male_sample_ids = paste("'", paste(sampleInfo[sampleInfo$gender == 'Male',]$call_call_set_name, collapse=",'", "'", sep=''))
male_sample_ids = gsub(" ", "", male_sample_ids, fixed = TRUE)
male_sub = list("_MALE_SAMPLE_IDS_" = male_sample_ids)
result <- DisplayAndDispatchQuery(query,
                                  project=project,
                                  replacements=c(tableReplacement, constraints, male_sub))
```

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

TODO Compare to command line

[Top](#top)

Ethnicity Inference
-----------------------------------
[Ethnicity Inference in Google Genomics](https://github.com/googlegenomics/codelabs/blob/master/R/PlatinumGenomes-QC/Sample-Level-QC.md#ethnicity-inference)

Status: Waiting for PCA release from Elmer.  See [here](https://github.com/elmer-garduno/spark-examples/tree/multiple_dataset_pca) for status.

[Top](#top)

Genome Similarity
-----------------------------------
[Genome Similarity in Google Genomics](https://github.com/googlegenomics/codelabs/blob/master/R/PlatinumGenomes-QC/Sample-Level-QC.md#genome-similarity)

Status: Installed apache maven and built Google Genomics dataflow package.  Ran into errors when trying to run dataflow job.

```
java -cp target/google-genomics-dataflow-v1beta2-0.2-SNAPSHOT.jar   com/google/cloud/genomics/dataflow/pipelines/IdentityByState   --project=my-project-id   --output=gs://my-bucket/localtest.txt   --genomicsSecretsFile=client_secrets.json
Error: Could not find or load main class com.google.cloud.genomics.dataflow.pipelines.IdentityByState
```

Command line
```
plink --genome --vcf chr22-merged.vcf.gz
```

[Top](#top)

Genotype Concordance
-----------------------------------

Calculate the concordance between sequencing data and genotyping data.

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

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