GWAS for fitness in Drosophila

Get the annotations for each DGRP variant

The following function temporarily loads the >1GB annotation file provided on the DGRP website at http://dgrp2.gnets.ncsu.edu/data/website/dgrp.fb557.annot.txt. We then extract the following variables for each variant, and save them in a SQLite database for memory-efficient searching inside R:

• The Flybase ID(s), if the variant is within or close to a gene
• The site class of the variant (e.g. intron, 5’-UTR…)
• The distance-to-gene in nucleotides (for UPSTREAM and DOWNSTREAM variants only)

get_variant_annotations <- function(){
  
  # Load up the big annotation file, get pertinent info. It's stored in some sort of text string format
  annot <- read.table("data/input/dgrp.fb557.annot.txt", header = FALSE, stringsAsFactors = FALSE)
  
  get.info <- function(rows){
    lapply(rows, function(row){
      site.class.field <- strsplit(annot$V3[row], split = "]")[[1]][1]
      num.genes <- str_count(site.class.field, ";") + 1
      output <- cbind(rep(annot$V1[row], num.genes), 
                      do.call("rbind", lapply(strsplit(site.class.field, split = ";")[[1]], 
                                              function(x) strsplit(x, split = "[|]")[[1]])))
      if(ncol(output) == 5) return(output[,c(1,2,4,5)]) # only return SNPs that have some annotation. Don't get the gene symbol
      else return(NULL)
    }) %>% do.call("rbind", .)
  }
  
  variant.details <- future_lapply(chunker(1:nrow(annot), max_chunk_size = 10000), get.info) %>% 
    do.call("rbind", .) %>% as.data.frame()
  
  names(variant.details) <- c("id", "FBID", "site.class", "distance.to.gene")
  variant.details$FBID <- unlist(str_extract_all(variant.details$FBID, "FBgn[:digit:]+")) # clean up text strings for Flybase ID
  variant.details %>%
    dplyr::filter(site.class != "FBgn0003638") %>% # NB this is a bug in the DGRP's annotation file
    mutate(chr = str_remove_all(substr(id, 1, 2), "_")) # get chromosome now for faster sorting later
}

Get the annotations for each Drosophila gene

The following function gets the annotations for the all the genes covered by DGRP variants, from the org.Dm.eg.db database object from Bioconductor. I don’t like the interface to those objects (it messes up dplyr), so I save the info as a second table in the SQLite database.

get_gene_annotations <- function(FBIDs_to_get){
  con <- dbconn(org.Dm.eg.db)
   tbl(con, "genes") %>%
    left_join(tbl(con, "flybase"), by = "_id") %>%
    left_join(tbl(con, "gene_info"), by = "_id") %>% 
    dplyr::select(flybase_id, gene_name, symbol, gene_id) %>%
    dplyr::rename(FBID = flybase_id, gene_symbol = symbol, entrez_id = gene_id) %>%
    dplyr::filter(FBID %in% FBIDs_to_get) %>%
    collect(n = Inf)
}

Create the SQLite database and add the two tables of annotations

if(file.exists("data/derived/annotations.sqlite3")) unlink("data/derived/annotations.sqlite3")

db <- src_sqlite("data/derived/annotations.sqlite3", create = TRUE)

db %>% copy_to(get_variant_annotations(), 
               "variants", temporary = FALSE, 
               indexes = list("id", "FBID", "chr", "site.class")) 

db %>% copy_to(get_gene_annotations(db %>% tbl("variants") %>% pull(FBID)), 
               "genes", temporary = FALSE)

Session information

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Sierra 10.12.6

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

locale:
[1] en_AU.UTF-8/en_AU.UTF-8/en_AU.UTF-8/C/en_AU.UTF-8/en_AU.UTF-8

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] bindrcpp_0.2.2       org.Dm.eg.db_3.6.0   AnnotationDbi_1.42.1
 [4] IRanges_2.14.10      S4Vectors_0.18.3     Biobase_2.40.0      
 [7] BiocGenerics_0.26.0  future.apply_1.0.1   future_1.9.0        
[10] stringr_1.3.1        dplyr_0.7.6         

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.18      dbplyr_1.2.2      compiler_3.5.1   
 [4] pillar_1.3.0      git2r_0.23.0      workflowr_1.1.1  
 [7] bindr_0.1.1       R.methodsS3_1.7.1 R.utils_2.7.0    
[10] tools_3.5.1       bit_1.1-14        digest_0.6.15    
[13] memoise_1.1.0     RSQLite_2.1.1     evaluate_0.11    
[16] tibble_1.4.2      pkgconfig_2.0.1   rlang_0.2.1      
[19] DBI_1.0.0         yaml_2.2.0        knitr_1.20       
[22] globals_0.12.2    bit64_0.9-7       rprojroot_1.3-2  
[25] tidyselect_0.2.4  glue_1.3.0        listenv_0.7.0    
[28] R6_2.2.2          rmarkdown_1.10    blob_1.1.1       
[31] purrr_0.2.5       magrittr_1.5      whisker_0.3-2    
[34] backports_1.1.2   codetools_0.2-15  htmltools_0.3.6  
[37] assertthat_0.2.0  stringi_1.2.4     crayon_1.3.4     
[40] R.oo_1.22.0