---
output:
html_document: default
pdf_document: default
---
---
title: 'RGMQL Example R Notebook: Use case 3'
author: "Silvia Cascianelli"
date: "`r Sys.Date()`"
output:
pdf_document: default
html_notebook: default
html_document: default
BiocStyle::html_document:
chunk_output_type: inline
---
```{r setup, include = FALSE}
knitr::opts_chunk$set(cache = TRUE)
```
In this example we investigate the ENCODE Chip-Seq narrow and broad data, extracted from the GMQL remote repository to automate all the steps needed to identify transcription factor (TF) high accumulation DNA zones using RGQML together with TFHAZ, another R/Bioconductor package . The knowledge of DNA regions in which transcription factors bind, in particular the HOT (High Occupancy Target) regions occupied by many different factors, is crucial to understand cancer genesis and develop new targeted therapies.
Load the RGMQL package and initialize the remote GMQL context of scalable data management engine, specifying remote_processing = TRUE, and, possibly, an authenticated login:
```{r, initialization}
library(RGMQL)
remote_url <- "http://www.gmql.eu/gmql-rest"
init_gmql(url = remote_url, remote_processing = TRUE)
#, username = 'XXXX', password = 'XXXX')
```
Download and extract the list of datasets in the curated remote repository and focus on those concerning ENCODE:
```{r, available_datasets}
dataset_list <- show_datasets_list(remote_url)
list <- unlist(lapply(dataset_list[["datasets"]], function(x) x$name))
grep(pattern = 'ENCODE', x = list, value = TRUE)
```
Select ChIP-seq data from the ENCODE NARROW dataset AUG_2017, aligned to HG19:
```{r, read_NARROW}
Enc_Narrow <- read_gmql("public.HG19_ENCODE_NARROW_AUG_2017",
is_local = FALSE)
HM_TF_rep_narrow <- filter(Enc_Narrow, assay == "ChIP-seq" & assembly
== "hg19" & project == "ENCODE" & file_status == "released" & biosample_term_name == "H1-hESC" & output_type == "optimal idr thresholded peaks")
```
Select ChIP-seq data from the latest ENCODE BROAD dataset AUG_2017, aligned to HG19 and:
```{r, read_BROAD}
Enc_Broad <- read_gmql("public.HG19_ENCODE_BROAD_AUG_2017",
is_local = FALSE)
HM_TF_rep_broad <- filter(Enc_Broad, assay == "ChIP-seq" & assembly ==
"hg19" & project == "ENCODE" & file_status ==
"released" & biosample_term_name == "H1-hESC" & output_type == "optimal idr thresholded peaks")
```
Take the union of the two previously generated datasets:
```{r, union}
HM_TF_rep <- union(HM_TF_rep_broad, HM_TF_rep_narrow)
```
Filter out samples subjected to pharmacological treatment or with specific "audit" marker:
```{r, cleaning}
HM_TF_rep_good_0 <- filter(HM_TF_rep, !biosample_treatments == "*" & !
(audit_error == "extremely low read depth" | audit_error == "extremely low read
length") & !(audit_warning == "insufficient read depth") & !(audit_not_compliant ==
"insufficient read depth" | audit_not_compliant =="insufficient replicate concordance" | audit_not_compliant == "missing input
control" | audit_not_compliant == "severe bottlenecking" | audit_not_compliant
== "unreplicated experiment"))
```
Filter out samples related to HM histone modifications:
```{r, discard_HM}
TF_rep_good_0 <- filter(HM_TF_rep_good_0, !(experiment_target == "H2AFZhuman" | experiment_target == "H3F3A-human" | experiment_target == "H3K27ac-human" | experiment_target == "H3K27me3-human" | experiment_target == "H3K36me3-human" | experiment_target == "H3K4me1-human" | experiment_target == "H3K4me2-human" | experiment_target == "H3K4me3-human" | experiment_target == "H3K79me2-human" | experiment_target == "H3K9ac-human" | experiment_target == "H3K9me1-human" | experiment_target == "H3K9me2-human" | experiment_target == "H3K9me3-human" | experiment_target == "H4K20me1-human"))
```
Create new regions attribute: Length of each region and for each sample, compute the number of regions and the sum of each region length just created:
```{r, removal}
TF_rep_good_1 <- select(TF_rep_good_0, regions_update = list(length = right - left))
TF_rep_good <- extend(TF_rep_good_1, Region_number = COUNT(),
sum_length = SUM("length"))
```
(1) Processing to obtain a dataset on which finding the filtering threshold;
```{r, filtering_threshold_1}
TF_rep_good_merged <- aggregate(TF_rep_good, groupBy =
conds(default = c("biosample_term_name")))
TF_rep_good_ordered <- arrange(TF_rep_good_merged,
regions_ordering = list(ASC("length")))
collect(TF_rep_good_ordered, name = "TF_rep_good_ordered")
job <- execute()
```
(1.1)Monitor the job status:
```{r, filtering_threshold_1_job_monitoring}
trace_job(remote_url , job$id)
```
(2) Processing to obtain a dataset on which finding the filtering threshold; once the job status is 'SUCCESS' download the resulting dataset obtained remotely in the working directory of the local File System:
```{r, filtering_threshold_2_download_in_FS}
dataset_name <- job$datasets[[1]]$name
print(dataset_name)
grl_TF_rep_good_ordered <- download_as_GRangesList(remote_url, dataset_name)
download_dataset(remote_url, datasetName = dataset_name, path = './Results_use_case_3')
```
(3) Processing to obtain a dataset on which finding the filtering threshold; once the needed dataset is saved, perform in the R environment all the remaining operations required to retrieve the threshold value:
```{r, filtering_threshold_3}
name_sample <- names(grl_TF_rep_good_ordered)# ? campione singolo
g <- grl_TF_rep_good_ordered[[name_sample]] #estraggo GRanges
Region_number_tot <- length(g) # conto tutte le regioni
n_up <- Region_number_tot * 0.95 #95percentile del numero di regioni
n_up_1 <- n_up + 1
index <- which(g$order >= ceiling(n_up) & g$order <= floor(n_up_1))
region <- g[index]
```
(4) Processing to obtain a dataset on which finding the filtering threshold; the chosen threshold is the length of the 95percentile of the region number:
```{r, filtering_threshold_4}
threshold <- region$length
threshold <- as.numeric(threshold)
threshold
```
Going back to RGQML remote processing, take only the regions with region lengths greater than 100 and smaller than the threshold:
```{r, threshold_filtering}
TF_rep_good_filtered_0 <- filter(TF_rep_good, r_predicate = length >= 100 & length <= threshold)
```
Create new metadata for each sample, with number of filtered regions and the sum of their lengths:
```{r, attributes_after_filtering}
TF_rep_good_filtered <- extend(TF_rep_good_filtered_0,
Region_number_filtered = COUNT(),
sum_length_filtered = SUM("length"))
```
Combine multiple replicate samples of the same TF experiment:
```{r, combine_TF_exp}
TF_0 <- cover(TF_rep_good_filtered, 1, ANY(), groupBy =
conds("experiment_target"))
```
Add new region attribute as length of each region after sample combination:
```{r, regions_update}
TF_1 <- select(TF_0, regions_update = list(length = right - left))
```
Create new metadata for each sample, with number of combined regions, and min, max and sum of their lengths:
```{r, attributes_after_cover}
TF <- extend(TF_1, Region_number_cover = COUNT(), sum_length_cover =
SUM("length"), min_length_cover = MIN("length"), max_length_cover = MAX("length"))
```
Materialize TF dataset into repository and downlad it on mass memory but also in main memory as GRangesList
```{r, main_execution}
collect(TF, name= "TF_res")
res <- execute()
```
```{r main_job_monitoring}
#Monitor job status:
trace_job(remote_url, res$id)
```
```{r, download_in_FS}
# Download dataset as folder:
res_name <- res$datasets[[1]]$name
download_dataset(remote_url, res_name, path = './Results_use_case_3')
```
```{r, GRangesList}
# Download dataset as GRangesList:
samples <- download_as_GRangesList(remote_url, res_name)
```
Post-processing before the analysis with TFHAZ
```{r, Post-processing}
TF=vector()
len_rep <- sapply(samples, function(x) len <- length(x))
TF_rep <- mapply(function(x, l){
exp <- x$experiment_target
TF_temp <- rep(exp, l)
TF <- append(TF, TF_temp)
}, samples@metadata, len_rep) #metadata(samples)
TF <- unlist(TF_rep)
H1_hESC <- unlist(samples)
data <- data.frame(H1_hESC, TF)
H1_hESC <- as(data, "GRanges")
```
After loading the TFHAZ package, find the transcription factor HOT DNA zones focusing on one chromosome at at time, by executing the following instructions:
```{r, TFHAZ_analysis}
library(TFHAZ)
TF_acc_21_w_0 <- accumulation(H1_hESC, "TF", "chr21", 0)
d_zones <- high_accumulation_zones(TF_acc_21_w_0, method =
"overlaps", threshold = "std")
print(d_zones)
```
Log out from remote engine:
```{r, logout}
logout_gmql(remote_url)
```
