cageminer: mining candidate genes by integrating GWAS and gene coexpression networks

class: center, middle, inverse, title-slide

# <p>cageminer: mining candidate genes by integrating GWAS and gene coexpression networks</p>
### Fabrício Almeida-Silva <span class="citation">@almeidasilvaf</span>
### UENF, Brazil
### August 6, 2021

---

## The problem

<br />

.center[GWAS can identify SNPs associated with traits, but not the underlying genes.

Hundreds of putative candidate genes. What then?
]

---

## *cageminer* to the rescue!

---

## Required input data

.panelset.sideways[
.panel[.panel-name[SNP positions]

A GRanges or GRangesList object.

```r
set.seed(123)
library(cageminer)
data(snp_pos)
snp_pos
```

```
## GRanges object with 116 ranges and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##     2    Chr02 149068682      *
##     3    Chr03   5274098      *
##     4    Chr05  27703815      *
##     5    Chr05  27761792      *
##     6    Chr05  27807397      *
##   ...      ...       ...    ...
##   114    Chr12 230514706      *
##   115    Chr12 230579178      *
##   116    Chr12 230812962      *
##   117    Chr12 230887290      *
##   118    Chr12 231022812      *
##   -------
##   seqinfo: 8 sequences from an unspecified genome; no seqlengths
```
]
.panel[.panel-name[Gene coordinates]

A GRanges object.

```r
data(gene_ranges)
head(gene_ranges)
```

```
## GRanges object with 6 ranges and 6 metadata columns:
##       seqnames        ranges strand |   source     type     score     phase
##          <Rle>     <IRanges>  <Rle> | <factor> <factor> <numeric> <integer>
##   [1]    Chr01   63209-63880      - |  PGA1.55     gene        NA      <NA>
##   [2]    Chr01 112298-112938      - |  PGA1.55     gene        NA      <NA>
##   [3]    Chr01 117979-118392      + |  PGA1.55     gene        NA      <NA>
##   [4]    Chr01 119464-119712      + |  PGA1.55     gene        NA      <NA>
##   [5]    Chr01 119892-120101      + |  PGA1.55     gene        NA      <NA>
##   [6]    Chr01 132987-134420      - |  PGA1.55     gene        NA      <NA>
##                ID          Parent
##       <character> <CharacterList>
##   [1]  CA01g00010                
##   [2]  CA01g00020                
##   [3]  CA01g00030                
##   [4]  CA01g00040                
##   [5]  CA01g00050                
##   [6]  CA01g00060                
##   -------
##   seqinfo: 12 sequences from an unspecified genome; no seqlengths
```
]
.panel[.panel-name[Guide genes]

Character vector or a 2-column data frame with genes in column 1 and gene functional class (e.g., pathway, GO) in column 2.

```r
data(guides)
head(guides)
```

```
##         Gene                               Description
## 1 CA10g07770                      response to stimulus
## 2 CA10g07770                        response to stress
## 3 CA10g07770             cellular response to stimulus
## 4 CA10g07770               cellular response to stress
## 6 CA10g07770 regulation of cellular response to stress
## 8 CA10g07770        regulation of response to stimulus
```
]
.panel[.panel-name[Expression data]

SummarizedExperiment object or matrix/data frame with genes in rows and samples in columns.

```r
data(pepper_se)
pepper_se
```

```
## class: SummarizedExperiment 
## dim: 3892 45 
## metadata(0):
## assays(1): ''
## rownames(3892): CA02g23440 CA02g05510 ... CA03g35110 CA02g12750
## rowData names(0):
## colnames(45): PL1 PL2 ... TMV-P0-3D TMV-P0-Up
## colData names(1): Condition
```
]

.panel[.panel-name[Coexpression network]

List object returned by `BioNERO::exp2gcn()`.

```r
# Infer network
gcn <- BioNERO::exp2gcn(pepper_se, SFTpower = 12)
```
]
]

---

## Exploratory analysis

<br />
🤔 &nbsp; Are SNPs evenly distributed across chromosomes?

.pull-left[

```r
plot_snp_distribution(snp_pos)
```
]
.pull-right[
<img src="slides_files/figure-html/unnamed-chunk-9-1.png" style="display: block; margin: auto;" />
]

---

## Exploratory analysis

<br />
🤔 &nbsp; Are SNPs evenly distributed across chromosomes?
 (Handling multiple traits)

.pull-left[

```r
# Simulate multiple traits by resampling
snp_list <- GenomicRanges::GRangesList(
  Trait1 = sample(snp_pos, 20),
  Trait2 = sample(snp_pos, 20),
  Trait3 = sample(snp_pos, 20),
  Trait4 = sample(snp_pos, 20)
)

# Visualize SNP distribution
plot_snp_distribution(snp_list)
```
]
.pull-right[
<img src="slides_files/figure-html/unnamed-chunk-10-1.png" style="display: block; margin: auto;" />
]

---

## Exploratory analysis

<br />
🤔 &nbsp; Are SNPs physically close to each other?

.pull-left[

```r
data(chr_length)
plot_snp_circos(chr_length, gene_ranges, 
                snp_pos)
```
]
.pull-right[
<img src="slides_files/figure-html/unnamed-chunk-11-1.png" style="display: block; margin: auto;" />
]

---

## Exploratory analysis

<br />
🤔 &nbsp; Are SNPs physically close to each other? (handling multiple traits)

.pull-left[

```r
data(chr_length)
plot_snp_circos(chr_length, gene_ranges, 
*               snp_list)
```
]
.pull-right[
<img src="slides_files/figure-html/unnamed-chunk-12-1.png" style="display: block; margin: auto;" />
]

---

## Gene mining

.panelset.sideways[

.panel[.panel-name[Step 1]

```r
candidates1 <- mine_step1(gene_ranges, snp_pos)
length(candidates1)
```

```
## [1] 1265
```
]
.panel[.panel-name[Step 2]

```r
candidates2 <- mine_step2(pepper_se, gcn = gcn, guides = guides$Gene,
                          candidates = candidates1$ID)
str(candidates2)
```

```
## List of 2
##  $ candidates: chr [1:37] "CA10g08490" "CA03g01790" "CA10g12640" "CA12g21230" ...
##  $ enrichment:'data.frame':	1 obs. of  7 variables:
##   ..$ TermID: chr "guide"
##   ..$ genes : num 323
##   ..$ all   : num 1303
##   ..$ pval  : num 2.58e-05
##   ..$ padj  : num 5.15e-05
##   ..$ GeneID: chr "CA01g00320,CA01g00480,CA01g00740,CA01g01120,CA01g01470,CA01g03330,CA01g03580,CA01g05790,CA01g06190,CA01g07040,C"| __truncated__
##   ..$ Module: chr "cyan"
```
]
.panel[.panel-name[Step 3]

```r
candidates3 <- mine_step3(pepper_se, candidates = candidates2$candidates,
                          sample_group = "PRR_stress")
```

```r
candidates3
```

```
##           gene      trait        cor       pvalue
## 243 CA12g18400 PRR_stress  0.5963394 1.540534e-05
## 187 CA10g02780 PRR_stress  0.3201048 3.205995e-02
## 19  CA02g16620 PRR_stress  0.3113993 3.732128e-02
## 33  CA02g16900 PRR_stress -0.3204388 3.187100e-02
## 103 CA03g03310 PRR_stress -0.3983772 6.720204e-03
```
]
.panel[.panel-name[Wrapping up]

`mine_candidates()` is a wrapper that calls all 3 `mine_step*()` functions sequentially.

```r
candidates <- mine_candidates(gene_ranges = gene_ranges, 
                              marker_ranges = snp_pos, 
                              exp = pepper_se,
                              gcn = gcn, guides = guides$Gene,
                              sample_group = "PRR_stress")
```
]
]

---

## Scoring and ranking candidates

<br />
Useful to pick the top *n* genes for validation, for instance.

.font140[
`$$S_i = r_{pb} \kappa$$`
]

.center[
where:

`$\kappa$` = 2 if the gene encodes a transcription factor

`$\kappa$` = 2 if the gene is a hub

`$\kappa$` = 3 is the gene is a hub + encodes a transcription factor

]

---

## Scoring and ranking candidates

```r
data(tfs)
data(hubs)

# Score and rank candidates
scored <- score_genes(candidates3, hubs$Gene, tfs$Gene_ID)
scored
```

```
##           gene      trait        cor       pvalue      score
## 243 CA12g18400 PRR_stress  0.5963394 1.540534e-05  0.5963394
## 103 CA03g03310 PRR_stress -0.3983772 6.720204e-03 -0.3983772
## 33  CA02g16900 PRR_stress -0.3204388 3.187100e-02 -0.3204388
## 187 CA10g02780 PRR_stress  0.3201048 3.205995e-02  0.3201048
## 19  CA02g16620 PRR_stress  0.3113993 3.732128e-02  0.3113993
```

---

## Summary

From 1265 genes, `cageminer` found 5 high-confidence candidates (>99% reduction).

---

class: sydney-yellow, middle, center

## Here's where you can find me: