In SimonLarsen/higana: Ontology-based analysis of genetic variants.

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
print(head(results$test, 10))

Visualizing results

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")