inst/doc/multi_tissue_analysis.R

In RITAN: Rapid Integration of Term Annotation and Network resources

## ----bgee1, echo=TRUE, eval=FALSE, warning=FALSE, fig.width = 7, fig.height = 7, fig.align='center'----
#  # Load RITAN
#  library(RITANdata)
#  library(RITAN)
#  
#  # Install the Bgee package. GO.db is a dependency of a dependency and may need to be installed seperately.
#  for (pkg in c('GO.db','BgeeDB','biomaRt')){
#    if (! (pkg %in% rownames(installed.packages()) )){
#      if (!requireNamespace("BiocManager", quietly=TRUE))
#          install.packages("BiocManager")
#      BiocManager::install(pkg)
#    }
#    library(pkg, character.only = TRUE)
#  }
#  
#  for (pkg in c('tidyselect','venn','magrittr','ggplot2','igraph')){
#    if (! (pkg %in% rownames(installed.packages()) )){
#      install.packages(pkg)
#    }
#    library(pkg, character.only = TRUE)
#  }
#  
#  # Setup Bgee query & get data (this may take some time)
#  bgee <- Bgee$new(species = "Homo_sapiens", dataType = "rna_seq", release = "13.2")
#  data <- getData(bgee)
#  e    <- formatData(bgee, data[[1]], callType = "present", stats = "rpkm")
#  
#  # Explore the dataset with: str(sampleNames(e)), str(featureNames(e)), str(phenoData(e))
#  table(phenoData(e)@data$Anatomical.entity.name)
#  
#  ## -------------------- -
#  ## Get expression in two tissues
#  tmp <- exprs(e)[ , phenoData(e)@data$Anatomical.entity.name == "heart" ]
#  i   <- apply( tmp, 1, function(x){ any(is.na(x)) })
#  expr_heart <- tmp[ !i, ]
#  
#  tmp <- exprs(e)[ , phenoData(e)@data$Anatomical.entity.name == "skeletal muscle tissue" ]
#  i   <- apply( tmp, 1, function(x){ any(is.na(x)) })
#  expr_skele <- tmp[ !i, ]
#  
#  venn::venn( list(Heart = rownames(expr_heart),
#                   Skeletal = rownames(expr_skele) ),
#              cexil= 1, cexsn = 1, zcolor = "style" )
#  
#  ## -------------------- -
#  ensembl <- useMart("ensembl", dataset = "hsapiens_gene_ensembl", "http://Aug2017.archive.ensembl.org" ) # version 90
#  
#  map_heart <- getBM( attributes = c('ensembl_gene_id','ensembl_transcript_id','hgnc_symbol'),
#                      filters = 'ensembl_gene_id', values = rownames(expr_heart), mart = ensembl )
#  
#  map_skele <- getBM( attributes = c('ensembl_gene_id','ensembl_transcript_id','hgnc_symbol'),
#                      filters = 'ensembl_gene_id', values = rownames(expr_skele), mart = ensembl )
#  
#  ## -------------------- -
#  ## Functions associated with each tissue's top genes
#  ## Important: the p-values reported here are observational, not inferential.
#  
#  mh <- apply( expr_heart, 1, mean )
#  top_heart <- map_heart$hgnc_symbol[ map_heart$ensembl_gene_id %in% rownames( expr_heart )[ mh > quantile(mh, .975) ] ] %>% setdiff(.,'')
#  
#  ms <- apply( expr_skele, 1, mean )
#  top_skele <- map_skele$hgnc_symbol[ map_skele$ensembl_gene_id %in% rownames( expr_skele )[ ms > quantile(ms, .975) ] ] %>% setdiff(.,'')
#  
#  e <- term_enrichment_by_subset( list( Heart    = top_heart,
#                                        Skeletal = top_skele ),
#                                  resources = 'GO_slim_PIR', all_symbols = cached_coding_genes )
#  
#  plot( e[ apply(e[, c(3:4)], 1, max) >= 12, ], cap=40, label_size_y = 8, wrap_y_labels = FALSE )
#  
#  ## -------------------- -
#  ## Network Interactions Within Each Tissue
#  
#  net_h <- network_overlap( top_heart, resources = c('CCSB','dPPI','HumanNet') )
#  net_s <- network_overlap( top_skele, resources = c('CCSB','dPPI','HumanNet') )
#  
#  net2g <- function(x){
#    edges <- as.matrix( x[, c(1,3)] )
#    G <- igraph::make_undirected_graph( c(t(edges)) )
#    return(G)
#  }
#  
#  g_h <- net2g( net_h )
#  g_s <- net2g( net_s )
#  
#  g_dif <- igraph::difference( g_h, g_s )
#  g_int <- igraph::intersection( g_h, g_s )
#  
#  cat(sprintf('
#  Of the top expressed genes, %d are shared and %d differ.
#  ', length(V(g_int)), length(V(g_dif)) ))
#  
#  par(mar=rep(0,4))
#  plot(g_dif, vertex.size = 2, vertex.label = NA, vertex.frame.color = 'white', layout = layout_nicely )