/
higana.Rmd
62 lines (48 loc) · 1.35 KB
/
higana.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
---
title: "Get started"
author: "Simon J. Larsen"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{get-started}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
## Preparing the ontology
```r
library(magrittr)
library(higana)
# Load 'go.obo' from http://geneontology.org/page/download-ontology
go <- read_obo("go.obo", c("is_a","part_of"))
# Load 'goa_human.gaf' from http://geneontology.org/page/download-go-annotations
anno <- read_gaf("goa_human.gaf", exclude_evidence="IEA")
go <- go %>%
filter_obsolete() %>%
unite_roots() %>%
annotate(anno) %>%
filter_unannotated() %>%
propagate_annotations() %>%
collapse_redundant_terms(threshold=0.9)
```
## Computing principal components
```r
library(snpStats)
go <- readRDS("ontology.rds")
bed <- read_bed("geno.bed", hg19")
pcs <- compute_term_pcs("terms.pcs", go, bed, num_parts=10)
```
## Performing association tests
```r
covars <- read.table("covars.tsv")
# Test each term
results <- test_terms("terms.pcs", class ~ sex+age+PC1+PC2+PC3+PC4, covars, go, family=binomial("logit"))
# Print top 10 terms
print(head(results$test, 10))
```
## Visualizing results
```r
library(DiagrammeR)
adj.p <- p.adjust(results$test$p.value, "bonferroni")
sig.terms <- results$test$term[adj.p < 0.05]
plot_hierarchy("tree.dot", go, terms=sig.terms, test=results)
grViz("tree.dot")
```