---
title: "02 Working with larger VCF: T2T by chromosome"
author: "Vincent J. Carey, stvjc at channing.harvard.edu"
date: "`r format(Sys.time(), '%B %d, %Y')`"
vignette: >
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteIndexEntry{02 Working with larger VCF: T2T by chromosome}
  %\VignetteEncoding{UTF-8}
output:
  BiocStyle::html_document:
    highlight: pygments
    number_sections: yes
    theme: united
    toc: yes
---

# Overview

In this document we illustrate some issues with
large data volumes.

Here is how we acquire the genotypes for the thousand
genomes samples based on the T2T reference.  We obtained
bgzipped vcf via
```
AnVIL::gsutil_cp("gs://fc-47de7dae-e8e6-429c-b760-b4ba49136eee/1KGP/joint_genotyping/joint_vcfs/raw/chr22.genotyped.vcf.gz", ".")
```
Then we used hail in python to deal with the conversion to MatrixTable.
We could have done this in R, but we had to learn how to manipulate
the 'reference sequence'.  We have a conjectural 'reference sequence'
json document in the json folder of the BiocHail package, used here
as `t2tAnVIL.json`.
```
>>> import hail as h
>>> rg = h.get_reference('GRCh38')
Initializing Hail with default parameters...
Setting default log level to "WARN".
To adjust logging level use sc.setLogLevel(newLevel). For SparkR, use setLogLevel(newLevel).
Running on Apache Spark version 3.1.3
SparkUI available at http://756809c79837:4040
Welcome to
     __  __     <>__
    / /_/ /__  __/ /
   / __  / _ `/ / /
  /_/ /_/\_,_/_/_/   version 0.2.105-acd89e80c345
LOGGING: writing to /home/rstudio/hail-20221213-1558-0.2.105-acd89e80c345.log
>>> nn = rg.read("t2tAnVIL.json")
>>> h.import_vcf("chr22.genotyped.vcf.gz", force_bgz=True, reference_genome=nn).write("t2t22.mt")
```
This operation seems to take a long time even with 64 cores.  We
could not get exact timing owing to some connectivity problems.

Here we'll work with the genotype data for T2T chr17.  We assume
that the MatrixTable is located at the path given by the environment
variable `HAIL_T2T_CHR17`.  This MatrixTable is available
in the Open Storage Network at osn:/bir190004-bucket01/Bioc1KGt2t/t17.zip.
This is a 45 GB file.  It should be unzipped at the location pointed to
by `HAIL_T2T_CHR17`.  Instructions for using rclone to acquire the zip file
are given in the appendix.

```{r doini,message=FALSE}
library(BiocHail)
```

(If the above command indicates that BiocHail is not available,
see the Installation section near the end of this document.)

```{r do17, eval=TRUE}
hl <- hail_init()
# NB the following two commands are now encapsulated in the rg_update function
nn <- hl$get_reference('GRCh38')
nn <- nn$read(system.file("json/t2tAnVIL.json", package="BiocHail"))
# updates the hail reference genome
bigloc = Sys.getenv("HAIL_T2T_CHR17")
if (nchar(bigloc)>0) {
  mt17 <- hl$read_matrix_table(Sys.getenv("HAIL_T2T_CHR17"))
  mt17$count()
}
```

# Population stratification assessment via PCA

The following command would compute PCs based on all SNP.
```
pcastuff = hl$hwe_normalized_pca(mt17$GT)
```
This seems impractical.  We have found that with a 64 core machine
at terra.bio, PCA on samples of 1-5% of the 3.8 million loci on T2T chr17
takes 40 minutes.  Hail's code seems good at utilizing all the cores.

We saved the PC scores for PCA based on 38k randomly sampled loci,
and 191k randomly sampled loci.  Here's a simple view of the latter.

```{r lk191,eval=TRUE}
utils::data(pcs_191k)
graphics::pairs(pcs_191k[,1:5], pch=".")
```

# Exercises

1. Comment on the gain in information about geographic origin
achieved by using a 5% sample of loci instead of a 1% sample.

2. Find the geographic origins of donors of the 1000 genomes genotypes
and bind them to `mt17` using the methods given in vignette `01_gwas_tut`.
Use these codes to color the points in the PCA plot.

3. Produce an artificial "phenotype" for these donors via `rnorm(3202,0,1)`,
bind it to the genotype data, and perform a naive GWAS.  What are the
loci most strongly associated with the artificial phenotype?

4. Produce a new artificial phenotype which has some association with
geographic origin of donor, but no association with genotype.
Produce a new naive GWAS, and a third using some of the PCA scores
as covariates.  What are the effects of this covariate adjustment on
reasoning about genetic association with the artificial phenotype?


# Appendix: using rclone with docker to get the chr17 data

It can be painful to install and configure rclone.  We use
a docker container.  Let RC_DATADIR be an environment variable
evaluating to an available folder.  

Also, place the text file with contents
```
[osn]
type = s3
provider = AWS
endpoint = https://mghp.osn.xsede.org
acl = public
no_check_bucket = true
```
in a file `rclone.conf` in a folder pointed to by the environment
variable `RC_CONFDIR`.


Then the following

```
docker run -v $RC_DATADIR:/data -v $RC_CONFDIR:/config/rclone -t rclone/rclone:latest ls osn:/bir190004-bucket01/Bioc1KGt2t
```

will list the files with 1KG samples genotyped against the T2T reference.



Use the rclone `copyto` command to obtain a local copy of the zip file `t17.zip`
in the folder pointed to by `$RC_DATADIR`:

```
docker run -v $RC_DATADIR:/data -v $RC_CONFDIR:/config/rclone -t rclone/rclone:latest copyto osn:/bir190004-bucket01/Bioc1KGt2t/t17.zip ./t17.zip
```
This file should be unzipped in a folder to which the environment variable `HAIL_T2T_CHR17` will
point.

# Installing `BiocHail`

`BiocHail` should be installed as follows:

```{r, eval = FALSE}
if (!require("BiocManager"))
    install.packages("BiocManager")
BiocManager::install("BiocHail")
```

As of 1.0.0, a JDK for Java version `<=` 11.0 is necessary
to use the version of Hail that is installed with the package.
This package should be usable on MacOS with suitable java
support.  If Java version `>=` 8.x is used, warnings from
Apache Spark may be observed.  To the best of our knowledge
the conditions to which the warnings pertain do not affect program performance.


# SessionInfo {-}

```{r sessionInfo}
sessionInfo()
```