sv-report.Rmd

---
title: Structural Variation Report
author: UPD-UCSC
header-includes:
- \usepackage{booktabs}
- \usepackage{makecell}
urlcolor: teal
output: html_document
---

```{r include=FALSE}
knitr::opts_chunk$set(echo=FALSE, message=FALSE, warning=FALSE, fig.width=10)
```

```{r setup}
## Read input arguments
## 1: VCF file from Sniffles
## 2: VCF file from SVIM
## 3: gene list file
## 4: gene and LoF intolerance score file
## 5: simple repeat annotation
## 6: long-read SV catalog
## 7: DGV SV catalog
## 8: ClinGen Pathogenic CNVs
## 9: CTCF peaks
## 10: common gnomAD SV but with no homozygous alternate
## 11: regulatory regions for kidney
args = commandArgs(TRUE)
## args = c('na12878-chr11.sniffles.ann.freqGnomADcov10.vcf','na12878-chr11.svim.ann.freqGnomADcov10.vcf','gene_position_info.tsv', 'gnomad.v2.1.1.lof_metrics.by_gene.txt.bgz', 'simpleRepeat.txt.gz', 'hsvlr.vcf.gz', 'GRCh38_hg38_variants_2016-08-31.txt', 'iscaPathogenic.txt.gz', 'ENCFF010WHH.bed.gz', 'gnomad_v2_sv.sites.pass.lifted.vcf.gz', 'ENCFF166QIT.lifted.bed.gz')

## Load packages
library(sveval)
library(VariantAnnotation)
library(knitr)
library(kableExtra)
library(DT)
library(dplyr)

out.format='html'
if(is_latex_output()){
  out.format='pdf'
}

## Reciproal overlap
rOverlap <- function(qgr, sgr){
  findOverlaps(qgr, sgr) %>% as.data.frame %>%
    mutate(qw=width(qgr)[queryHits], sw=width(sgr)[subjectHits],
           qsw=width(pintersect(qgr[queryHits], sgr[subjectHits])),
           rol=ifelse(qw>sw, qsw/qw, qsw/sw))
}

## Parse the SnpEff field EFF
parseEff <- function(eff){
  effect = gsub('(.*)\\(.*\\)', '\\1', eff)
  target = gsub('.*\\((.*)\\)', '\\1', eff)
  target = strsplit(target, '\\|')[[1]]
  res = tibble(effect=effect, impact=target[[1]], gene=target[[6]],
               target.type=target[[8]])
  res = res %>% mutate(effect=ifelse(grepl('FUSION', effect), paste0(effect, '_', gene), effect))
  res2 = res %>% filter(grepl('+', gene)) %>% mutate(gene=gsub('.*\\+(.*)', '\\1', gene))
  res = res %>% mutate(gene=gsub('(.*)\\+.*', '\\1', gene)) %>% rbind(res2)
  return(res)
}
parseEff.nbpc <- function(eff){
  target = gsub('.*\\((.*)\\)', '\\1', eff)
  target = strsplit(target, '\\|')[[1]]
  res = tibble(gene=target[[6]], target.type=target[[8]])
  return(nrow(unique(subset(res, target.type=='CODING'))))
}
## Format links for genes (NCBI), genomic position (UCSC Genome Browser), or pli (gnomAD)
formatLink <- function(labels, urls, format='pdf'){
  if(format=='html'){
    return(paste0('[', labels, '](', urls,')'))
  }
  if(format=='rawhtml'){
    return(paste0('<a href="', urls, '" target="_blank">', labels, '</a>'))
  }
  if(format=='pdf'){
    return(paste0('{\\href{', urls, '}{', labels, '}}'))
  }
}
linkGenes <- function(genes, format='html'){
  urls = paste0('https://www.ncbi.nlm.nih.gov/gene?term=(', genes, '%5BGene%20Name%5D)%20AND%20Human%5BOrganism%5D')
  return(formatLink(genes, urls, format=format))
}
linkPli <- function(scores, genes, digits=3, format='html'){
  urls = paste0('https://gnomad.broadinstitute.org/gene/', genes)
  labels = round(scores, digits)
  return(formatLink(labels, urls, format=format))
}
linkPos <- function(chr, pos, type, size, flanks=500, format='html'){
  size = ifelse(type %in% c('DEL','DUP','INV'), size, 1)
  urls = paste0('https://genome.ucsc.edu/cgi-bin/hgTracks?db=hg38&position=', chr, '%3A',
                pos-flanks, '-', pos+size+flanks, '&highlight=hg38.', chr, '%3A',
                pos, '-', pos+size, '%23E18F99')
  labels = paste0(chr, ':', pos)
  return(formatLink(labels, urls, format=format))
}

## Make a data.frame for a particular gene from an effect data.frame and the two VCF objects
makeGeneTable <- function(gene.name, eff.df, vcf.sn, vcf.sv){
  eff.df = eff.df %>% filter(gene==gene.name)
  vars = tibble()
  if(any(eff.df$method=='Sniffles')){
    id.ii = eff.df %>% filter(method=='Sniffles') %>% .$id %>% unique
    vcf = vcf.sn[id.ii]
    vars = tibble(
      id=id.ii,
      chr=as.character(seqnames(rowRanges(vcf))),
      pos=start(rowRanges(vcf)),
      type=info(vcf)$SVTYPE,
      size=abs(info(vcf)$SVLEN),
      reads=info(vcf)$RE,
      freq=info(vcf)$AFMAX,
      simp.rep=info(vcf)$ol.simple.repeat,
      hsvlr=info(vcf)$HSVLR,
      comhet=info(vcf)$ComHET,
      dgv=info(vcf)$dgv.loss | info(vcf)$dgv.gain,
      clingen=info(vcf)$clingen,
      ctcf=info(vcf)$ctcf,
      cres=info(vcf)$cres,
      gt=geno(vcf)$GT[,1],
      method=info(vcf)$METHOD,
      nmeths=info(vcf)$METHODS) %>% rbind(vars)
  }
  if(any(eff.df$method=='SVIM')){
    id.ii = eff.df %>% filter(method=='SVIM') %>% .$id %>% unique
    vcf = vcf.sv[id.ii]
    vars = tibble(
      id=id.ii,
      chr=as.character(seqnames(rowRanges(vcf))),
      pos=start(rowRanges(vcf)),
      type=info(vcf)$SVTYPE,
      size=abs(info(vcf)$SVLEN),
      reads=info(vcf)$SUPPORT,
      freq=info(vcf)$AFMAX,
      simp.rep=info(vcf)$ol.simple.repeat,
      hsvlr=info(vcf)$HSVLR,
      comhet=info(vcf)$ComHET,
      dgv=info(vcf)$dgv.loss | info(vcf)$dgv.gain,
      clingen=info(vcf)$clingen,
      ctcf=info(vcf)$ctcf,
      cres=info(vcf)$cres,
      gt=geno(vcf)$GT[,1],
      method=info(vcf)$METHOD,
      nmeths=info(vcf)$METHODS) %>% rbind(vars)
  }
  vars = eff.df %>% mutate(effect=gsub('\\+', ' \\+ ', effect)) %>%
    group_by(method, id) %>% summarize(effect=paste(unique(effect), collapse='; ')) %>%
    merge(vars, .) %>% ungroup %>% select(-id)
  vars$pli = eff.df$pLI[1]
  vars$gene = gene.name
  vars %>% dplyr::select(gene, pli, everything())
}
## Format some fields of a data.frame to convert them into links to external resources
formatGeneTable <- function(df){
  if(nrow(df)==0) return(df)
  if('type' %in% colnames(df)){
    df$type = factor(df$type)
  }
  if('gt' %in% colnames(df)){
    df$gt = factor(df$gt)
  }
  if('pos' %in% colnames(df)){
    df$pos = linkPos(df$chr, df$pos, df$type, df$size, format='rawhtml')
    df$chr = NULL
  }
  if('gene' %in% colnames(df)){
    if('pli' %in% colnames(df)){
      df$pli = linkPli(df$pli, df$gene, format='rawhtml')
    }
    df$gene = linkGenes(df$gene, format='rawhtml')
  }
  df
}

## Read both VCFs
vcf.sn = readSVvcf(args[1], vcf.object=TRUE)
vcf.sv = readSVvcf(args[2], vcf.object=TRUE)

## Make GRanges object for the overlap
gr.sn = rowRanges(vcf.sn)
end(gr.sn) = info(vcf.sn)$END
gr.sn$size = info(vcf.sn)$SIZE
gr.sn$type = info(vcf.sn)$SVTYPE
gr.sn$GT = 'hom'
gr.sn$ID = 1:length(vcf.sn)
gr.sv = rowRanges(vcf.sv)
end(gr.sv) = info(vcf.sv)$END
gr.sv$size = info(vcf.sv)$SIZE
gr.sv$type = info(vcf.sv)$SVTYPE
gr.sv$GT = 'hom'
gr.sv$ID = 1:length(vcf.sv)

## Overlap and add a column with overlap
ol.o = svOverlap(gr.sn, gr.sv, max.ins.dist=100)
ol.o$query = subset(ol.o$query, cov.prop>.5)
ol.o$subject = subset(ol.o$subject, cov.prop>.5)
info(vcf.sn)$METHODS = 1
info(vcf.sn)$METHODS[ol.o$query$ID] = 2
info(vcf.sv)$METHODS = 1
info(vcf.sv)$METHODS[ol.o$subject$ID] = 2
info(vcf.sn)$METHOD = 'Sniffles'
info(vcf.sv)$METHOD = 'SVIM'

## Read gene list
genes = read.table(args[3], as.is=TRUE, header=TRUE) %>% dplyr::rename(gene=gene_name) %>%
  mutate(gene_list=factor(gene_list))
genes.list = genes %>% select(gene, gene_list) %>% unique
genes.gr = makeGRangesFromDataFrame(genes, keep.extra.columns=TRUE)

## Read pli score
pli.df = read.table(args[4], as.is=TRUE, header=TRUE, sep='\t')
pli.df = pli.df %>% select(gene, pLI)

## Simple repeats
sr = read.table(args[5])
sr = sr[,c(2,3,4,6,7,17)]
colnames(sr) = c('chrom', 'start', 'end', 'period', 'copyNum', 'sequence')
sr = makeGRangesFromDataFrame(sr, keep.extra.columns=TRUE)
info(vcf.sn)$ol.simple.repeat = overlapsAny(vcf.sn, sr, maxgap=10)
info(vcf.sv)$ol.simple.repeat = overlapsAny(vcf.sv, sr, maxgap=10)

## SVs from public long-read studies
lr.gr = readSVvcf(args[6], right.trim=FALSE)
ol.o = svOverlap(gr.sn, lr.gr, max.ins.dist=100)
info(vcf.sn)$HSVLR = FALSE
info(vcf.sn)$HSVLR[subset(ol.o$query, cov.prop>.5)$ID] = TRUE
ol.o = svOverlap(gr.sv, lr.gr, max.ins.dist=100)
info(vcf.sv)$HSVLR = FALSE
info(vcf.sv)$HSVLR[subset(ol.o$query, cov.prop>.5)$ID] = TRUE

## Common SVs with no homozygous individuals in the gnomAD database
gnomad = readSVvcf(args[10], right.trim=FALSE, vcf.object=TRUE)
af = sapply(info(gnomad)$AF, '[', 1)
gnomad = gnomad[which(info(gnomad)$FREQ_HOMALT==0 & af>.01)]
chet.gr = rowRanges(gnomad)
chet.gr$size = info(gnomad)$SIZE
chet.gr$type = info(gnomad)$SVTYPE
ol.o = svOverlap(gr.sn, chet.gr, max.ins.dist=100)
info(vcf.sn)$ComHET = FALSE
info(vcf.sn)$ComHET[subset(ol.o$query, cov.prop>.5)$ID] = TRUE
ol.o = svOverlap(gr.sv, chet.gr, max.ins.dist=100)
info(vcf.sv)$ComHET = FALSE
info(vcf.sv)$ComHET[subset(ol.o$query, cov.prop>.5)$ID] = TRUE

## DGV catalog
dgv = read.table(args[7], as.is=TRUE, header=TRUE, sep='\t')
dgv = dgv[,c('chr','start','end', 'variantsubtype')]
dgv$chr = paste0('chr', dgv$chr)
dgv = makeGRangesFromDataFrame(dgv, keep.extra.columns=TRUE)
dgv.loss = subset(dgv, variantsubtype %in% c('loss', 'deletion', 'gain+loss'))
info(vcf.sn)$dgv.loss = overlapsAny(vcf.sn, dgv.loss, maxgap=10)
info(vcf.sv)$dgv.loss = overlapsAny(vcf.sv, dgv.loss, maxgap=10)
dgv.gain = subset(dgv, variantsubtype %in% c('duplication', 'gain', 'insertion',
                                             'mobile element insertion',
                                             'novel sequence insertion',
                                             'tandem duplication', 'gain+loss'))
info(vcf.sn)$dgv.gain = overlapsAny(vcf.sn, dgv.gain, maxgap=10)
info(vcf.sv)$dgv.gain = overlapsAny(vcf.sv, dgv.gain, maxgap=10)

## ClinGenn pathogenic CNVs
clingen = read.table(args[8], as.is=TRUE, sep='\t')
clingen = GRanges(clingen$V2, IRanges(clingen$V3, clingen$V4))
ol.o = rOverlap(gr.sn, clingen)
info(vcf.sn)$clingen = FALSE
info(vcf.sn)$clingen[subset(ol.o, rol>.5)$queryHits] = TRUE
ol.o = rOverlap(gr.sv, clingen)
info(vcf.sv)$clingen = FALSE
info(vcf.sv)$clingen[subset(ol.o, rol>.5)$queryHits] = TRUE

## CTCF peaks
ctcf = read.table(args[9], as.is=TRUE)
ctcf = with(ctcf, GRanges(V1, IRanges(V2, V3), score=V5))
info(vcf.sn)$ctcf = overlapsAny(vcf.sn, ctcf)
info(vcf.sv)$ctcf = overlapsAny(vcf.sv, ctcf)

## Kidney regulatory regions
cres = read.table(args[11], as.is=TRUE)
cres = with(cres, GRanges(V1, IRanges(V2, V3), score=V5))
info(vcf.sn)$cres = overlapsAny(vcf.sn, cres)
info(vcf.sv)$cres = overlapsAny(vcf.sv, cres)
```

```{r highimpact}
eff.df = lapply(list(vcf.sn, vcf.sv), function(vcf){
  high.ii = which(unlist(lapply(info(vcf)$EFF, function(effs) any(grepl('HIGH', effs)))))
  moderate.ii = which(unlist(lapply(info(vcf)$EFF, function(effs) any(grepl('MODERATE', effs)))))
  cons.ii = which(unlist(lapply(info(vcf)$EFF, function(effs) any(grepl('CONSERVED', effs)))))
  cds.ii = which(unlist(lapply(info(vcf)$EFF, function(effs) any(grepl('CDS', effs)))))
  reg.ii = which(unlist(lapply(info(vcf)$EFF, function(effs) any(grepl('REGULATION', effs)))))
  union.ii = unique(c(high.ii, moderate.ii, cons.ii, cds.ii, reg.ii))
  union.ii = unique(c(high.ii))
  eff.df = lapply(union.ii, function(ii){
    effs = lapply(info(vcf)$EFF[[ii]], parseEff)
    effs = do.call(rbind, effs)
    effs$id = ii
    effs$nmeths = info(vcf)$METHODS[ii]
    effs
  })
  eff.df = do.call(rbind, eff.df)
  eff.df$method = info(vcf)$METHOD[1]
  eff.df
})
eff.df = do.call(rbind, eff.df)
```

# {.tabset}

## Methods

The structural variants found by [Sniffles](https://github.com/fritzsedlazeck/Sniffles) and [SVIM](https://github.com/eldariont/svim) were annotated by [SnpEff](http://snpeff.sourceforge.net/) and with the frequency of variants in the [gnoma-SV catalog](https://macarthurlab.org/2019/03/20/structural-variants-in-gnomad/).
To focus on high-confidence SV calls, we only look at calls made by both Sniffles and SVIM.

The variants' "effects" in this report come from SnpEff and follow the format [described in SnpEff documentation](http://snpeff.sourceforge.net/SnpEff_manual.html#eff).

The *pLI* score was computed by the [gnomAD project](https://gnomad.broadinstitute.org/).
It represents the probability that the gene is intolerant to loss-of-function variants.
A variant affecting a gene with a high pLI score (e.g. >0.9) is more likely to have a biological impact.
The variants in each section are ordered to show those affecting genes with the highest pLI first.

In the following, we also removed common variants, i.e. either:

- frequency higher than 1% in gnomAD-SV
- seen in the SV catalog from long-read studies

Clarifications about some column names:

- `pli` prob of  loss-of-function intolerance described above.
- `type` SV type
- `gt` genotype, usually `0/1` for het or `1/1` for hom.
- `freq` allele frequency of similar SV in gnomAD-SV (>10K genomes sequenced with Illumina whole-genome sequencing).
- `dgv` does the variant overlap any variant in DGV?
- `clingen` any similar ClinGen pathogenic variants (as found at the UCSC Genome Browser)? "Similar" defined as reciprocal overlap > 50%.
- `ctcf` does the variant overlap a CTCF binding site? From ENCODE track for kidney.
- `cres` does the variant overlap a regulatory region? From ENCODE track for kidney.
- `simp.rep` is the variant in or close to a simple repeat (see simple repeat track in the UCSC Genome Browser).
- `hsvlr` any similar variant in SV catalogs from public long-read sequencing studies?
- `comhet` any similar variant in gnomAD-SV catalogs that are common but no homozygous individual?
- `effect` predicted effect extracted from SnpEff annotation.
- `distance` distance to the gene boundaries.
- `nb.pc.genes` number of protein-coding genes overlapped by the variant.

## Genes of interest

There are `r length(unique(genes$gene))` genes in `r length(unique(genes$gene_list))` gene lists:

```{r genes}
genes %>% group_by(gene_list) %>% summarize(genes=n()) %>% kable
```

---

They contain the following genes:

```{r genes2}
genes %>% formatGeneTable %>% datatable(filter='top', escape=FALSE, options=list(pageLength=50))
```

## High impact {.tabset}

*HIGH* impact as defined by SnpEff.

### Genes of interest

```{r selhigh}
high.df = eff.df %>% filter(impact=='HIGH', gene %in% genes$gene, nmeths==2) %>%
  merge(pli.df, all.x=TRUE) %>% arrange(desc(pLI)) %>% unique
  
### Write table
t.df = lapply(unique(high.df$gene), function(gene.name){
  makeGeneTable(gene.name, high.df, vcf.sn, vcf.sv)
})
t.df = do.call(rbind, t.df)

if(!is.null(t.df)){
  t.df %>% filter(freq<.01, !hsvlr) %>% merge(genes.list) %>% formatGeneTable %>%
    select(gene_list, gene, pli, type, size, pos, gt, freq, dgv, clingen, simp.rep, effect) %>% 
    datatable(filter='top', escape=FALSE, options=list(pageLength=50))
}
```

### Other coding genes

```{r othershigh}
high.df = eff.df %>% filter(impact=='HIGH', target.type=='CODING', nmeths==2,
                            !(gene %in% genes$gene)) %>%
  merge(pli.df, all.x=TRUE) %>% arrange(desc(pLI)) %>% unique

### Write table
t.df = lapply(unique(high.df$gene), function(gene.name){
  makeGeneTable(gene.name, high.df, vcf.sn, vcf.sv)
})
t.df = do.call(rbind, t.df)

if(!is.null(t.df)){
  t.df %>% filter(freq<.01, !hsvlr) %>% formatGeneTable %>%
    select(gene, pli, type, size, pos, gt, freq, dgv, clingen, simp.rep, effect) %>% 
    datatable(filter='top', escape=FALSE, options=list(pageLength=50))
}
```

### Genes of interest - Common SVs but no homozygous controls

```{r selhighcomhet}
high.df = eff.df %>% filter(impact=='HIGH', gene %in% genes$gene, nmeths==2) %>%
  merge(pli.df, all.x=TRUE) %>% arrange(desc(pLI)) %>% unique
  
### Write table
t.df = lapply(unique(high.df$gene), function(gene.name){
  makeGeneTable(gene.name, high.df, vcf.sn, vcf.sv)
})
t.df = do.call(rbind, t.df)

if(!is.null(t.df)){
  t.df %>% filter(comhet) %>% merge(genes.list) %>% formatGeneTable %>%
    select(gene_list, gene, pli, type, size, pos, gt, freq, dgv, clingen, simp.rep, effect) %>% 
    datatable(filter='top', escape=FALSE, options=list(pageLength=50))
}
```

### Other coding genes - Common SVs but no homozygous controls

```{r othershighcomhet}
high.df = eff.df %>% filter(impact=='HIGH', target.type=='CODING', nmeths==2,
                            !(gene %in% genes$gene)) %>%
  merge(pli.df, all.x=TRUE) %>% arrange(desc(pLI)) %>% unique

### Write table
t.df = lapply(unique(high.df$gene), function(gene.name){
  makeGeneTable(gene.name, high.df, vcf.sn, vcf.sv)
})
t.df = do.call(rbind, t.df)

if(!is.null(t.df)){
  t.df %>% filter(comhet) %>% formatGeneTable %>%
    select(gene, pli, type, size, pos, gt, freq, dgv, clingen, simp.rep, effect) %>% 
    datatable(filter='top', escape=FALSE, options=list(pageLength=50))
}
```

## Moderate impact {.tabset}

### Genes of interest

*MODERATE* or *MODIFIER* impacts as defined by SnpEff.

```{r selmod}
eff.g.df = lapply(list(vcf.sn, vcf.sv), function(vcf){
  gene.ii = unlist(lapply(genes$gene, function(genen){
    which(unlist(lapply(info(vcf)$EFF, function(effs) any(grepl(genen, effs)))))
  }))
  eff.df = lapply(unique(gene.ii), function(ii){
    effs = lapply(info(vcf)$EFF[[ii]], parseEff)
    effs = do.call(rbind, effs)
    effs$id = ii
    effs$nmeths = info(vcf)$METHODS[ii]
    effs
  })
  eff.df = do.call(rbind, eff.df)
  eff.df$method = info(vcf)$METHOD[1]
  eff.df
})
eff.g.df = do.call(rbind, eff.g.df)

high.df = eff.g.df %>% filter(impact %in% c('MODERATE', 'MODIFIER'),
                              gene %in% genes$gene, nmeths==2) %>%
  merge(pli.df, all.x=TRUE) %>% arrange(desc(pLI)) %>% unique
  
### Write table
t.df = lapply(unique(high.df$gene), function(gene.name){
  makeGeneTable(gene.name, high.df, vcf.sn, vcf.sv)
})
t.df = do.call(rbind, t.df)

if(!is.null(t.df)){
  t.df %>% filter(freq<.01, !hsvlr) %>% merge(genes.list) %>% formatGeneTable %>%
    select(gene_list, gene, pli, type, size, pos, gt, freq, dgv, clingen, cres, ctcf, simp.rep, effect) %>% 
    datatable(filter='top', escape=FALSE, options=list(pageLength=50))
}
```

### Other coding genes

*MODERATE* impact or either *CONSERVED* or *CDS* or *REGULATION* effects (as defined by SnpEff).

```{r othersmod}
high.df = rbind(
  eff.df %>% filter(impact=='MODERATE', nmeths==2, !(gene %in% genes$gene)),
  eff.df %>% filter(grepl('CONSERVED', effect) | effect %in% c('CDS','REGULATION'),
                    nmeths==2, !(gene %in% genes$gene))) %>%
  merge(pli.df, all.x=TRUE) %>% arrange(desc(pLI)) %>% unique
  
### Write table
t.df = lapply(unique(high.df$gene), function(gene.name){
  makeGeneTable(gene.name, high.df, vcf.sn, vcf.sv)
})
t.df = do.call(rbind, t.df)

if(!is.null(t.df)){
  t.df %>% filter(freq<.01, !hsvlr) %>% formatGeneTable %>%
    select(gene, pli, type, size, pos, gt, freq, dgv, clingen, cres, ctcf, simp.rep, effect) %>% 
    datatable(filter='top', escape=FALSE, options=list(pageLength=50))
}
```

### Genes of interest - Common SVs but no homozygous controls

*MODERATE* or *MODIFIER* impacts as defined by SnpEff.

```{r selmodcomhet}
eff.g.df = lapply(list(vcf.sn, vcf.sv), function(vcf){
  gene.ii = unlist(lapply(genes$gene, function(genen){
    which(unlist(lapply(info(vcf)$EFF, function(effs) any(grepl(genen, effs)))))
  }))
  eff.df = lapply(unique(gene.ii), function(ii){
    effs = lapply(info(vcf)$EFF[[ii]], parseEff)
    effs = do.call(rbind, effs)
    effs$id = ii
    effs$nmeths = info(vcf)$METHODS[ii]
    effs
  })
  eff.df = do.call(rbind, eff.df)
  eff.df$method = info(vcf)$METHOD[1]
  eff.df
})
eff.g.df = do.call(rbind, eff.g.df)

high.df = eff.g.df %>% filter(impact %in% c('MODERATE', 'MODIFIER'),
                              gene %in% genes$gene, nmeths==2) %>%
  merge(pli.df, all.x=TRUE) %>% arrange(desc(pLI)) %>% unique
  
### Write table
t.df = lapply(unique(high.df$gene), function(gene.name){
  makeGeneTable(gene.name, high.df, vcf.sn, vcf.sv)
})
t.df = do.call(rbind, t.df)

if(!is.null(t.df)){
  t.df %>% filter(comhet) %>% merge(genes.list) %>% formatGeneTable %>%
    select(gene_list, gene, pli, type, size, pos, gt, freq, dgv, clingen, cres, ctcf, simp.rep, effect) %>% 
    datatable(filter='top', escape=FALSE, options=list(pageLength=50))
}
```

### Other coding genes - Common SVs but no homozygous controls

*MODERATE* impact or either *CONSERVED* or *CDS* or *REGULATION* effects (as defined by SnpEff).

```{r othersmodcomhet}
high.df = rbind(
  eff.df %>% filter(impact=='MODERATE', nmeths==2, !(gene %in% genes$gene)),
  eff.df %>% filter(grepl('CONSERVED', effect) | effect %in% c('CDS','REGULATION'),
                    nmeths==2, !(gene %in% genes$gene))) %>%
  merge(pli.df, all.x=TRUE) %>% arrange(desc(pLI)) %>% unique
  
### Write table
t.df = lapply(unique(high.df$gene), function(gene.name){
  makeGeneTable(gene.name, high.df, vcf.sn, vcf.sv)
})
t.df = do.call(rbind, t.df)

if(!is.null(t.df)){
  t.df %>% filter(comhet) %>% formatGeneTable %>%
    select(gene, pli, type, size, pos, gt, freq, dgv, clingen, cres, ctcf, simp.rep, effect) %>% 
    datatable(filter='top', escape=FALSE, options=list(pageLength=50))
}
```


## Rare SVs {.tabset}

### Around genes of interest

```{r raresel}
vcf.rare = subset(vcf.sn, METHODS==2 & AFMAX<.01)
d.df = distanceToNearest(vcf.rare, genes.gr) %>% as.data.frame
info(vcf.rare)$SELGENEDIST = NA
info(vcf.rare)$SELGENEDIST[d.df$queryHits] = d.df$distance
info(vcf.rare)$SELGENE = NA
info(vcf.rare)$SELGENE[d.df$queryHits] = genes.gr$gene[d.df$subjectHits]
## vcf.rare = subset(vcf.rare, SELGENEDIST<1e6)

svs.df = tibble(
  gene=info(vcf.rare)$SELGENE,
  chr=as.character(seqnames(rowRanges(vcf.rare))),
  pos=start(rowRanges(vcf.rare)),
  type=info(vcf.rare)$SVTYPE,
  size=abs(info(vcf.rare)$SVLEN),
  distance=info(vcf.rare)$SELGENEDIST,
  freq=info(vcf.rare)$AFMAX,
  simp.rep=info(vcf.rare)$ol.simple.repeat,
  hsvlr=info(vcf.rare)$HSVLR,
  comhet=info(vcf.rare)$ComHET,
  dgv=info(vcf.rare)$dgv.loss | info(vcf.rare)$dgv.gain,
  clingen=info(vcf.rare)$clingen,
  ctcf=info(vcf.rare)$ctcf,
  cres=info(vcf.rare)$cres,
  gt=geno(vcf.rare)$GT[,1],
  nb.pc.genes = unlist(lapply(info(vcf.rare)$EFF, parseEff.nbpc)))
if(info(vcf.rare)$METHOD[1]=='Sniffles'){
  svs.df$reads=info(vcf.rare)$RE
} else {
  svs.df$reads=info(vcf.rare)$SUPPORT
}

svs.sel = svs.df %>% filter(distance<1e5) %>% arrange(distance)
t.df = lapply(unique(svs.sel$gene), function(gene.name){
  pli = subset(pli.df, gene==gene.name)$pLI
  vars = svs.sel %>% filter(gene==gene.name) %>% select(-gene)
  vars$pli = ifelse(is.null(pli), NA, pli[1])
  vars$gene = gene.name
  vars %>% dplyr::select(gene, pli, everything())
})
t.df = do.call(rbind, t.df)

if(!is.null(t.df)){
  t.df %>% filter(freq<.01, !hsvlr) %>% merge(genes.list) %>% formatGeneTable %>%
    select(gene_list, gene, pli, distance, type, size, pos, gt, freq, dgv, clingen, cres, ctcf, simp.rep) %>% 
    datatable(filter='top', escape=FALSE, options=list(pageLength=50))
}
```

### Around genes of interest - Common SVs but no homozygous controls

```{r rareselcomhet}
if(!is.null(t.df)){
  t.df %>% filter(comhet) %>% merge(genes.list) %>% formatGeneTable %>%
    select(gene_list, gene, pli, distance, type, size, pos, gt, freq, dgv, clingen, cres, ctcf, simp.rep) %>% 
    datatable(filter='top', escape=FALSE, options=list(pageLength=50))
}
```

### Around genes of interest - Transposable element insertions

```{r meisel}
if(!is.null(t.df)){
  t.df %>% filter(freq<.01, !hsvlr, type=='INS', abs(300-size)<100 | abs(6000-size)<500) %>%
    merge(genes.list) %>% formatGeneTable %>%
    select(gene_list, gene, pli, distance, size, pos, gt, freq, dgv, clingen, cres, ctcf, simp.rep) %>% 
      datatable(filter='top', escape=FALSE, options=list(pageLength=50))
  t.df %>% filter(freq<.01, !hsvlr, type=='INS', abs(300-size)<100 | abs(6000-size)<500) %>%
    merge(genes.list) %>% write.table(file='sv-te-like-insertions.tsv', sep='\t', row.names=FALSE,
                                      quote=FALSE)
}
```

### Large

Large and rare variants tend to have a higher biological impact.
Interesting profiles to look for:

- large deletions spanning a CTCF binding region (could lead to TAD reorganization and ectopic gene expression).
- Overlap with pathogenic CNV in the ClinGen database. *clingen* column shows variants with 50% reciprocal overlap with a pathogenic CNV.
- Large SV affecting multiple protein-coding genes.

```{r rarelarge}
svs.df %>% filter(size>1e3, freq<.01, !hsvlr) %>%
  arrange(desc(abs(size))) %>% 
  select(-gene) %>% formatGeneTable %>%
  select(type, size, pos, gt, freq, dgv, clingen, cres, ctcf, simp.rep, nb.pc.genes) %>% 
  datatable(filter='top', escape=FALSE, options=list(pageLength=50))
```