Get started

Simon J. Larsen

Source: vignettes/higana.Rmd

Preparing the ontology

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

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

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
test <- results$test[order(results$test$p.value),]
print(head(test))

Visualizing results

library(DiagrammeR)

test$adj.p <- p.adjust(test$p.value, "bonferroni")
p.adj <- p.adjust(pvalues, "bonferroni")
sig.terms <- names(which(p.adj < 0.05))
plot_hierarchy(go, "tree.dot", sig.terms, pvalues=p.adj)
grViz("tree.dot")