In dzhang32/ODER: Optimising the Definition of Expressed Regions

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(tidyverse)
devtools::load_all("~/packages/ODER")
# library(ODER)

• For this example, we will use the GTEx v6 data made available through the recount2 project.
• To begin, we require a bigwig file and total AUC for each sample. Both of which can be obtained through the recount R package.
• For speed, this workflow only demonstrates this analysis on chromosome 21 & 22 using data from only 2 GTEx cerebellar samples. Though it would be advisable to use the bigwigs from as many samples as you have to gain a more accurate estimate of mean coverage.

Download data

gtex_metadata <- recount::all_metadata('gtex') 

gtex_metadata <- gtex_metadata %>% 
  as.data.frame() %>% 
  dplyr::filter(smtsd == "Brain - Cerebellar Hemisphere", 
                smafrze == "USE ME")

# for(i in seq_along(gtex_metadata$bigwig_file)){
# 
#   download.file(url = str_c("http://duffel.rail.bio/recount/SRP012682/bw/", gtex_metadata$bigwig_file[i]),
#                 destfile = str_c("~/data/", gtex_metadata$bigwig_file[i]),
#                 method = "wget", 
#                 quiet = T)
# 
# }

Quick-start

• If you are only interested in the output (a set of ERs with optimised definitions) then you can just run the wrapper function ODER::ODER().

opt_ERs <- 
  ODER(bw_paths = stringr::str_c("~/data/", gtex_metadata$bigwig_file)[1:2], 
       aucs = gtex_metadata$auc[1:2],
       chr_to_filter = stringr::str_c("chr", 21:22),
       target_auc = (40e6 * 100), # target 40 million reads with 100 bp length
       MCCs = c(5, 10), 
       MRGs = c(10, 20), 
       gtf = "/data/references/ensembl/gtf_gff3/v95/Homo_sapiens.GRCh38.95.gtf", 
       ucsc_chr = T, 
       ignore.strand = T) # bws are not strand-specific, set this to TRUE

opt_ERs

---

Step-by-step

• Below shows the intermediate functions and outputs that ODER::ODER() calls.

Calculate mean coverage across your samples

chrs_mean_cov <-
  calc_mean_cov(bw_paths = str_c("~/data/", gtex_metadata$bigwig_file)[1:2],
                aucs = gtex_metadata$auc[1:2],
                target_auc = (40e6 * 100), # target 40 million coverage with 100 bp length reads
                chr_to_filter = c("chr21", "chr22")) # defaults to chr1-22, chrX, chrY, chrM

chrs_mean_cov

Generate ERs across different MCCs and MRGs

ERs_MCCs_MRGs <- gen_ERs(chrs_mean_cov, 
                         MCCs = c(5, 10), 
                         MRGs = c(10, 20))

ERs_MCCs_MRGs$MCC_5$MRG_10

Get non-overlapping exons

exons_no_overlap_gr <- get_no_overlap_exons(gtf = "/data/references/ensembl/gtf_gff3/v95/Homo_sapiens.GRCh38.95.gtf", 
                                            ucsc_chr = T, ignore.strand = T)

head(exons_no_overlap_gr)

Calculate delta between ERs and exons

delta_df <- calc_ERs_delta(ERs_MCCs_MRGs, 
                           opt_gr = exons_no_overlap_gr)

delta_df

Select optimal ERs

opt_ERs <- get_opt_ERs(ERs_MCCs_MRGs, delta_df)

opt_ERs

---

dzhang32/ODER documentation built on April 23, 2021, 11:21 p.m.