2018_Meetup.Rmd

---
title: "CNV Analysis Workflow"
author: Sehyun Oh
date: "November 29, 2018"
output:
  BiocStyle::html_document:
    number_sections: yes
    toc: true
vignette: >
  %\VignetteIndexEntry{Analysis and Resources of Multi-Assay Experiments}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, eval = FALSE)
```

# Prepare input files
```{r echo=FALSE}
system.file("extdata", package="PureCN")
```

```{bash echo=FALSE}
export PURECN="/path/to/PureCN/extdata"
```

## Interval file

**Download GRCh38 reference file (without ALT contigs)**
```{bash}
wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.15_GRCh38/seqs_for_alignment_pipelines.ucsc_ids/GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz
```

```{bash}
Rscript $PURECN/IntervalFile.R \
    --infile baits_hg38.bed \
    --fasta hg38.fasta \
    --outfile $OUT/baits_hg38_intervals.txt \
    --offtarget --genome hg38 \
    --mappability GRCh38_no_alt_76.bw \
    --force
```

*** 
## VCF and stats files

**Overview of MuTect outputs**   
(https://gatkforums.broadinstitute.org/gatk/discussion/4231/what-is-the-output-of-mutect-and-how-should-i-interpret-it)

* Call-stats file
The main output which people typically work with is the "call-stats" file. It is an exhaustive report of all the metrics and statistics available about the calls made by MuTect and the filters that are applied internally by default. 

* VCF file of candidate mutations
Upon request, MuTect can output a summary VCF file containing the mutation candidates annotated with `KEEP` or `REJECT` in the `FILTER` field.
```{bash}
java -jar mutect.jar \
    --analysis_type MuTect \
    -R hg38.fasta \
    --dbsnp $DBSNP_VCF \
    --cosmic $COSMIC_VCF \
# -I:normal $BAM_NORMAL \
    -I:tumor $BAM_TUMOR  \
    -o $OUT/${SAMPLEID}_mutect_stats.txt \
    -vcf $OUT/${SAMPLEID}_mutect.vcf
```

*** 
## Pool of normals

**About PON from GATK Forum**   
(https://gatkforums.broadinstitute.org/gatk/discussion/11053/panel-of-normals-pon)

A Panel of Normal or PON is a type of resource used in somatic variant analysis. Depending on the type of variant you're looking for, the PON will be generated differently. What all PONs have in common is that (1) they are made from normal samples (in this context, "normal" means derived from healthy tissue that is believed to not have any somatic alterations) and (2) their main purpose is to capture recurrent technical artifacts in order to improve the results of the variant calling analysis.

As a result, the most important selection criteria for choosing normals to include in any PON are the technical properties of how the data was generated. It's very important to use normals that are as technically similar as possible to the tumor (same exome or genome preparation methods, sequencing technology and so on). Additionally, the samples should come from subjects that were young and healthy to minimize the chance of using as normal a sample from someone who has an undiagnosed tumor. Normals are typically derived from blood samples.

There is no definitive rule for how many samples should be used to make a PON (even a small PON is better than no PON) but in practice we recommend aiming for a minimum of 40.

**`--artifact-detection-mode`**  
-   used when running the caller on a normal (as if it were a  tumor) to detect artifacts      
-   include variant calls that are clearly germline   
-   used for Panel-Of-Normals creation   

```{bash}
java -jar mutect.jar \
   --analysis_type MuTect \
   -R hg38.fasta \
   --artifact_detection_mode \
   --dbsnp $DBSNP_VCF \
   --cosmic $COSMIC_VCF \
   -dt None \
   -I:tumor $BAM_NORMAL  \
   -o $OUT/${SAMPLEID}_pon_stats.txt \
   -vcf $OUT/${SAMPLEID}_pon.vcf

java -jar GenomeAnalysisTK.jar \
   --analysis_type SelectVariants \
   -R hg38.fasta \
   --exclude_sample_expressions none \
   -V $OUT/${SAMPLEID}_bwa_mutect_artifact_detection_mode.vcf \
   -o $OUT/${SAMPLEID}_bwa_mutect_artifact_detection_mode_no_none.vcf 

java -jar GenomeAnalysisTK.jar \
   -T CombineVariants \
   --minimumN 5 \
   --genotypemergeoption UNSORTED \
   --variant $OUT/${SAMPLEID}_bwa_mutect_artifact_detection_mode_no_none.vcf \
   -o $OUT/normals.merged.min5.vcf

bgzip $OUT/normals.merged.min5.vcf
tabix $OUT/normals.merged.min5.vcf.gz
```

*** 
## GC-normalized coverage files
This process is applied to both tumor and (process-matched) normal bam files, separately. 
```{bash}
Rscript $PURECN/Coverage.R \
   --outdir $OUT/normal_cov \
   --bam $BAM_NORMAL \
   --intervals $OUT/baits_hg38_intervals.txt

Rscript $PURECN/Coverage.R \
   --outdir $OUT/tumor_cov \
   --bam $BAM_TUMOR \
   --intervals $OUT/baits_hg38_intervals.txt
```

*** 
## NormalDB
To build a normal database for coverage normalization, copy the paths to all GC-normalized normal coverage files in a single text file, line-by-line:

```{bash}
ls -a $OUT/normal_cov/*_coverage_loess.txt | cat > $OUT/normalDB/normalDB.list
```

Recommendations: 

* Do not mix normal data obtained with different capture kits!   
* Provide a normal panel VCF here to precompute mapping bias for faster runtimes.
```{bash}
Rscript $PURECN/NormalDB.R \
   --outdir $OUT/normalDB \
   --coveragefiles $OUT/normalDB/normalDB.list \
   --normal_panel $OUT/normals.merged.min5.vcf.gz \
   --genome hg38 --force
```

*** 
## Other files
### manifest
```{r message=FALSE}
library(TCGAutils)
library(jsonlite)
library(curl)
library(downloader)
library(GenomicDataCommons)
library(magrittr)
```
Check the available values for GDCQuery filters
```{r}
available_values('files','cases.project.project_id')
available_values('files','experimental_strategy')
available_values('files','data_format')
```

```{r}
# Prepare GDC manifest file
manifest <- GenomicDataCommons::files() %>%
  GenomicDataCommons::filter(~ cases.project.project_id == "TCGA-LUAD" &
           experimental_strategy == "WXS" &
           data_format == "BAM") %>%
  GenomicDataCommons::manifest()

# Translate study identifiers from UUID to barcode
manifest <- cbind(manifest, 
                 UUIDtoBarcode(manifest$id, id_type = "file_id", 
                               end_point = "center"))
names(manifest)[7] <- "barcode"

# Extract biospecimen data from the TCGA barcode
manifest <- cbind(manifest, TCGAutils::TCGAbiospec(manifest$barcode))
```

For a quick test run, just use the first 10 samples for the following annotation
```{r}
manifest_all <- manifest
manifest <- manifest[c(1:10),]
```

Target_capture_kit information
```{r}
res <- lapply(manifest$id, function(uuid) {
  con = curl::curl(paste0("https://api.gdc.cancer.gov/files/", uuid, "?pretty=true&fields=analysis.metadata.read_groups.target_capture_kit_target_region,analysis.metadata.read_groups.target_capture_kit_name,analysis.metadata.read_groups.target_capture_kit_vendor,analysis.metadata.read_groups.target_capture_kit_catalog_number"))
  x = jsonlite::fromJSON(con)
  return(x)
})

y <- lapply(res, function(x) unique(x$data$analysis$metadata$read_groups))
y <- do.call(rbind, y)
manifest <- cbind(manifest, y)
```

Cleaning target_region (BED file) information:   
Seperate the name of bedfiles from 'target_capture_kit_target_region'
```{r}
bedfiles <- tail(unlist(strsplit(as.character(manifest$target_capture_kit_target_region)[1], split="/", fixed=TRUE)), n = 1)
manifest$bedfiles <- bedfiles
```

*** 
### liftOver
`hg19ToHg38` chain file can be downloaded from here:   
http://hgdownload.soe.ucsc.edu/goldenPath/hg19/liftOver/hg19ToHg38.over.chain.gz
```{r}
data.dir <- "/path/to/data" 
BED <- "SeqCap_EZ_Exome_v3_capture"
  
# Change chromosome notation from '1, 2, 3, ...' to 'chr1, chr2, chr3, ...
bed <- file.path(data.dir, paste0(BED, ".bed"))

library(rtracklayer)
ch <- import.chain(file.path(data.dir, "hg19ToHg38.over.chain"))
bed_gr <- import(bed, format = "bed")
bed_hg38 <- liftOver(bed_gr, ch)
bed_hg38 <- unlist(bed_hg38)

export.bed(bed_hg38, file.path(data.dir, paste0(BED,"_hg38.bed")))
```

*** 
### mappability
**Download GEM library**
```{bash}
wget https://sourceforge.net/projects/gemlibrary/files/gem-library/Binary%20pre-release%203/GEM-binaries-Linux-x86_64-core_i3-20130406-045632.tbz2/download
tar xvf download
```

```{bash}
# To calculate mappability, set kmer size to length of mapped reads
THREADS=24
KMER=76

PREF="/path/to/GCA_000001405.15_GRCh38_no_alt_analysis_set"
REFERENCE="${PREF}.fna"

gem-indexer -T ${THREADS} -i ${REFERENCE} -o ${PREF}_index
gem-mappability -T ${THREADS} -I ${PREF}_index.gem -l ${KMER} -o ${PREF}_${KMER} -m 2 -e 2
gem-2-wig -I ${PREF}_index.gem -i ${PREF}_${KMER}.mappability -o ${PREF}_${KMER}

# Convert to bigWig format, for example using the UCSC wigToBigWig tool
cut -f1,2 ${REFERENCE}.fai > ${PREF}.sizes 

# I found the unexpected letter, "AC" in my .wig file --> remove it
cp ${PREF}_${KMER}.wig GRCh38_no_alt_${KMER}.wig
sed -e s/AC//g -i GRCh38_no_alt_${KMER}.wig
wigToBigWig GRCh38_no_alt_${KMER}.wig ${PREF}.sizes ${PREF}_${KMER}.bw
```

**Download GRCh38 mappability files**   
For GRCh38, you can download recommended 76-kmer or 100-kmer mappability files.   
https://s3.amazonaws.com/purecn/GCA_000001405.15_GRCh38_no_alt_analysis_set_76.bw
https://s3.amazonaws.com/purecn/GCA_000001405.15_GRCh38_no_alt_analysis_set_100.bw

*** 
### snpblacklist
Generate `hg38_simpleRepeats.bed` file
```{r}
downloadFromUCSC <- TRUE
if (downloadFromUCSC) {
  library(rtracklayer)
  mySession <- browserSession("UCSC")
  genome(mySession) <- "hg38"
  simpleRepeats <- track(ucscTableQuery(mySession,
                                        track="Simple Repeats", 
                                        table="simpleRepeat"))
  export(simpleRepeats, "hg38_simpleRepeats.bed")
}
snp.blacklist <- import("hg38_simpleRepeats.bed", format = "bed")
```

# PureCN

```{bash}
Rscript $PURECN/PureCN.R \
    --out $OUT/$SAMPLEID \
    --tumor $OUT/tumor_cov/${SAMPLEID}_coverage_loess.txt \
    --SAMPLEID ${SAMPLEID} \
    --vcf $OUT/${SAMPLEID}_mutect.vcf \
    --statsfile $OUT/${SAMPLEID}_mutect_stats.txt \
    --normaldb $OUT/normalDB/normalDB_hg38.rds \
    --normal_panel $OUT/normalDB/mapping_bias_hg38.rds \
    --intervals $OUT/baits_hg38_intervals.txt \
    --intervalweightfile $OUT/normalDB/interval_weights_hg38.txt \
    --snpblacklist hg38_simpleRepeats.bed \
    --genome hg38 \
    --force --postoptimize --seed 123
```

# Output
Restrict mutation burden calculation to coding sequences
```{bash}
# Generate a BED file with callable regions
java -jar GenomeAnalysisTK.jar \
   -T CallableLoci \
   -R hg38.fasta \
   -I:tumor $BAM_TUMOR  \
   --summary $OUT/${SAMPLEID}_table.txt \
   -o $OUT/${SAMPLEID}_callable_status.bed \
   --minDepth 30

# BED file with callable regions
grep CALLABLE $OUT/${SAMPLEID}_callable_status.bed > \
   $OUT/${SAMPLEID}_callable_status_filtered.bed

# Restrict mutation burden calculation to coding sequences
Rscript $PURECN/FilterCallableLoci.R \
   --genome hg38 \
   --infile $OUT/${SAMPLEID}_callable_status_filtered.bed \
   --outfile $OUT/${SAMPLEID}_callable_status_filtered_cds.bed
```

Dx.R extracts copy number and mutation metrics from PureCN.R output.
```{bash}
Rscript $PURECN/Dx.R \
   --out $OUT/$SAMPLEID/${SAMPLEID} \
   --rds $OUT/$SAMPLEID/${SAMPLEID}.rds \
   --callable $OUT/${SAMPLEID}_callable_status_filtered_cds.bed \
   --exclude hg38_simpleRepeats.bed \
   --signatures
```