diff --git a/analyses/focal-cn-file-preparation/02-add-ploidy-consensus.Rmd b/analyses/focal-cn-file-preparation/02-add-ploidy-consensus.Rmd index 489f83cef3..7dc1733e25 100644 --- a/analyses/focal-cn-file-preparation/02-add-ploidy-consensus.Rmd +++ b/analyses/focal-cn-file-preparation/02-add-ploidy-consensus.Rmd @@ -34,6 +34,10 @@ library(tidyverse) ```{r} scratch_dir <- file.path("..", "..", "scratch") +output_dir <- file.path(scratch_dir, "cytoband_status", "segments") +if(!dir.exists(output_dir)) { + dir.create(output_dir, recursive = TRUE) +} ``` ```{r} @@ -146,3 +150,69 @@ add_status_df %>% output_file <- file.path("..", "..", "scratch", "consensus_seg_with_status.tsv") write_tsv(add_status_df, output_file) ``` + +### Prepare separate bed files for losses/gains for bedtools coverage function + +```{r} +bed_status_df <- add_status_df %>% + select(chrom, loc.start, loc.end, everything()) %>% + group_by(Kids_First_Biospecimen_ID) %>% + arrange(chrom, loc.start, loc.end) + +losses_bed_status_df <- bed_status_df %>% + filter(status == "loss") + +gains_bed_status_df <- bed_status_df %>% + filter(status == "gain") + +# make lists of data frames by sample +bed_status_list <- bed_status_df %>% + group_split() +names(bed_status_list) <- bed_status_df %>% + group_keys() %>% + pull(Kids_First_Biospecimen_ID) + +bed_loss_list <- losses_bed_status_df %>% + group_split() +names(bed_loss_list) <- losses_bed_status_df %>% + group_keys() %>% + pull(Kids_First_Biospecimen_ID) + +bed_gain_list <- gains_bed_status_df %>% + group_split() +names(bed_gain_list) <- gains_bed_status_df %>% + group_keys() %>% + pull(Kids_First_Biospecimen_ID) +``` + +### Write to file + +```{r} +temp <- purrr::imap(bed_status_list, + ~ write_tsv(.x, + file.path( + output_dir, paste0("consensus_callable.", .y, ".bed") + ), + col_names = FALSE)) + +temp_loss <- purrr::imap(bed_loss_list, + ~ write_tsv(.x, + file.path( + output_dir, paste0("consensus_loss.", .y, ".bed") + ), + col_names = FALSE)) + +temp_gain <- purrr::imap(bed_gain_list, + ~ write_tsv(.x, + file.path( + output_dir, paste0("consensus_gain.", .y, ".bed") + ), + col_names = FALSE)) +``` + +## Session Info + +```{r} +sessionInfo() +``` + diff --git a/analyses/focal-cn-file-preparation/02-add-ploidy-consensus.nb.html b/analyses/focal-cn-file-preparation/02-add-ploidy-consensus.nb.html index 57e9efd613..723380dc69 100644 --- a/analyses/focal-cn-file-preparation/02-add-ploidy-consensus.nb.html +++ b/analyses/focal-cn-file-preparation/02-add-ploidy-consensus.nb.html @@ -2955,8 +2955,8 @@

Libraries and functions

library(tidyverse)
- -
── Attaching packages ────────────────────────────────────── tidyverse 1.2.1 ──
+ +
── Attaching packages ────────────────────────────────────────────── tidyverse 1.2.1 ──
✔ ggplot2 3.2.0     ✔ purrr   0.3.2
@@ -2964,8 +2964,8 @@ 

Libraries and functions

✔ tidyr 0.8.3 ✔ stringr 1.4.0 ✔ readr 1.3.1 ✔ forcats 0.4.0
- -
── Conflicts ───────────────────────────────────────── tidyverse_conflicts() ──
+
+
── Conflicts ───────────────────────────────────────────────── tidyverse_conflicts() ──
 ✖ dplyr::filter() masks stats::filter()
 ✖ dplyr::lag()    masks stats::lag()
@@ -2976,8 +2976,12 @@

Libraries and functions

Files and directories

- -
scratch_dir <- file.path("..", "..", "scratch")
+ +
scratch_dir <- file.path("..", "..", "scratch")
+output_dir <- file.path(scratch_dir, "cytoband_status", "segments")
+if(!dir.exists(output_dir)) {
+  dir.create(output_dir, recursive = TRUE)
+}
@@ -3113,7 +3117,7 @@

Does seg.mean agree with status?

geom_jitter()
-

+

@@ -3142,10 +3146,127 @@

Write to file

write_tsv(add_status_df, output_file) + +
+

Prepare separate bed files for losses/gains for bedtools coverage function

+ + + +
bed_status_df <- add_status_df %>%
+  select(chrom, loc.start, loc.end, everything()) %>%
+  group_by(Kids_First_Biospecimen_ID) %>%
+  arrange(chrom, loc.start, loc.end)
+
+losses_bed_status_df <- bed_status_df %>%
+  filter(status == "loss")
+
+gains_bed_status_df <- bed_status_df %>%
+  filter(status == "gain")
+
+# make lists of data frames by sample
+bed_status_list <- bed_status_df %>%
+  group_split()
+names(bed_status_list) <- bed_status_df %>%
+  group_keys() %>%
+  pull(Kids_First_Biospecimen_ID)
+
+bed_loss_list <- losses_bed_status_df %>%
+  group_split()
+names(bed_loss_list) <- losses_bed_status_df %>%
+  group_keys() %>%
+  pull(Kids_First_Biospecimen_ID)
+
+bed_gain_list <- gains_bed_status_df %>%
+  group_split()
+names(bed_gain_list) <- gains_bed_status_df %>%
+  group_keys() %>%
+  pull(Kids_First_Biospecimen_ID)
+ + + +
+
+

Write to file

+ + + +
temp <- purrr::imap(bed_status_list,
+                    ~ write_tsv(.x,
+                                file.path(
+                                  output_dir, paste0("consensus_callable.", .y, ".bed")
+                                ),
+                                col_names = FALSE))
+
+temp_loss <- purrr::imap(bed_loss_list,
+                         ~ write_tsv(.x,
+                                     file.path(
+                                       output_dir, paste0("consensus_loss.", .y, ".bed")
+                                     ),
+                                     col_names = FALSE))
+
+temp_gain <- purrr::imap(bed_gain_list,
+                         ~ write_tsv(.x,
+                                     file.path(
+                                       output_dir, paste0("consensus_gain.", .y, ".bed")
+                                     ),
+                                     col_names = FALSE))
+ + + +
+ +
+

Session Info

+ + + +
sessionInfo()
+ + +
R version 3.6.0 (2019-04-26)
+Platform: x86_64-pc-linux-gnu (64-bit)
+Running under: Debian GNU/Linux 9 (stretch)
+
+Matrix products: default
+BLAS/LAPACK: /usr/lib/libopenblasp-r0.2.19.so
+
+locale:
+ [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
+ [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
+ [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=C             
+ [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
+ [9] LC_ADDRESS=C               LC_TELEPHONE=C            
+[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
+
+attached base packages:
+[1] stats     graphics  grDevices utils     datasets  methods   base     
+
+other attached packages:
+[1] forcats_0.4.0   stringr_1.4.0   dplyr_0.8.3     purrr_0.3.2    
+[5] readr_1.3.1     tidyr_0.8.3     tibble_2.1.3    ggplot2_3.2.0  
+[9] tidyverse_1.2.1
+
+loaded via a namespace (and not attached):
+ [1] Rcpp_1.0.1       cellranger_1.1.0 pillar_1.4.2     compiler_3.6.0  
+ [5] base64enc_0.1-3  tools_3.6.0      digest_0.6.20    lubridate_1.7.4 
+ [9] jsonlite_1.6     evaluate_0.14    nlme_3.1-140     gtable_0.3.0    
+[13] lattice_0.20-38  pkgconfig_2.0.2  rlang_0.4.0      cli_1.1.0       
+[17] rstudioapi_0.10  yaml_2.2.0       haven_2.1.1      xfun_0.8        
+[21] withr_2.1.2      xml2_1.2.0       httr_1.4.0       knitr_1.23      
+[25] generics_0.0.2   hms_0.4.2        grid_3.6.0       tidyselect_0.2.5
+[29] glue_1.3.1       R6_2.4.0         readxl_1.3.1     rmarkdown_1.13  
+[33] modelr_0.1.4     magrittr_1.5     backports_1.1.4  scales_1.0.0    
+[37] htmltools_0.3.6  rvest_0.3.4      assertthat_0.2.1 colorspace_1.4-1
+[41] labeling_0.3     stringi_1.4.3    lazyeval_0.2.2   munsell_0.5.0   
+[45] broom_0.5.2      crayon_1.3.4    
+ + + +
-
LS0tCnRpdGxlOiAiQWRkIHBsb2lkeSBjb2x1bW4sIHN0YXR1cyB0byBjb25zZW5zdXMgU0VHIGZpbGUiCm91dHB1dDogCiAgaHRtbF9ub3RlYm9vazoKICAgIHRvYzogVFJVRQogICAgdG9jX2Zsb2F0OiBUUlVFCmF1dGhvcjogQ2hhbnRlIEJldGhlbGwgYW5kIEphY2x5biBUYXJvbmkgZm9yIEFMU0YgQ0NETApkYXRlOiAyMDIwCi0tLQoKVGhlIGBwYnRhLWhpc3RvbG9naWVzLnRzdmAgZmlsZSBjb250YWlucyBhIGB0dW1vcl9wbG9pZHlgIGNvbHVtbiwgd2hpY2ggaXMgdHVtb3IgcGxvaWR5IGFzIGluZmVycmVkIGJ5IENvbnRyb2xGcmVlQy4KVGhlIGNvcHkgbnVtYmVyIGluZm9ybWF0aW9uIHNob3VsZCBiZSBpbnRlcnByZXRlZCBpbiB0aGUgbGlnaHQgb2YgdGhpcyBpbmZvcm1hdGlvbiAoc2VlOiBbY3VycmVudCB2ZXJzaW9uIG9mIERhdGEgRm9ybWF0cyBmaWxlXShodHRwczovL2dpdGh1Yi5jb20vQWxleHNMZW1vbmFkZS9PcGVuUEJUQS1hbmFseXNpcy9ibG9iLzBlNjQyZWYwNjAyYWQwNGI4ZDBmYmQ4MGJlNjI2ZDQzNzRmZWU5MjgvZG9jL2RhdGEtZm9ybWF0cy5tZCNhLW5vdGUtb24tcGxvaWR5KSkuCgpUaGlzIG5vdGVib29rIGFkZHMgcGxvaWR5IGluZm9ybWF0aW9uIHRvIHRoZSBjb25zZW5zdXMgU0VHIGZpbGUgYW5kIGFkZHMgYSBzdGF0dXMgY29sdW1uIHRoYXQgZGVmaW5lcyBnYWluIGFuZCBsb3NzIGJyb2FkbHkuCioqTm90ZSoqIHRoYXQgdGhlIGNvbnNlbnN1cyBTRUcgZmlsZSBkb2VzIG5vdCBoYXZlIGNvcHkgbnVtYmVyIGluZm9ybWF0aW9uIGZvciBzZXggY2hyb21vc29tZXMuCgojIyBVc2FnZQoKVGhpcyBub3RlYm9vayBpcyBpbnRlbmRlZCB0byBiZSBydW4gZnJvbSB0aGUgY29tbWFuZCBsaW5lIHdpdGggdGhlIGZvbGxvd2luZyAoYXNzdW1lcyB5b3UgYXJlIGluIHRoZSByb290IGRpcmVjdG9yeSBvZiB0aGUgcmVwb3NpdG9yeSk6CgpgYGAKUnNjcmlwdCAtZSAicm1hcmtkb3duOjpyZW5kZXIoJ2FuYWx5c2VzL2ZvY2FsLWNuLWZpbGUtcHJlcGFyYXRpb24vMDItYWRkLXBsb2lkeS1jb25zZW5zdXMuUm1kJywgY2xlYW4gPSBUUlVFKSIKYGBgCgojIyBTZXQgdXAKCiMjIyBMaWJyYXJpZXMgYW5kIGZ1bmN0aW9ucwoKYGBge3J9CmxpYnJhcnkodGlkeXZlcnNlKQpgYGAKCiMjIyBGaWxlcyBhbmQgZGlyZWN0b3JpZXMKCmBgYHtyfQpzY3JhdGNoX2RpciA8LSBmaWxlLnBhdGgoIi4uIiwgIi4uIiwgInNjcmF0Y2giKQpgYGAKCmBgYHtyfQojIFRPRE86IHRoZSBjb25zZW5zdXMgU0VHIGZpbGUgaXMgbm90IGN1cnJlbnRseSBpbiB0aGUgZGF0YSBkb3dubG9hZCAtLSB3aGVuIGl0CiMgZ2V0cyBpbmNsdWRlZCB3ZSB3aWxsIGhhdmUgdG8gY2hhbmdlIHRoZSBmaWxlIHBhdGggaGVyZQpjb25zZW5zdXNfc2VnX2ZpbGUgPC0gZmlsZS5wYXRoKCIuLiIsICJjb3B5X251bWJlcl9jb25zZW5zdXNfY2FsbCIsICJyZXN1bHRzIiwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgInBidGEtY252LWNvbnNlbnN1cy5zZWcuZ3oiKQpoaXN0b2xvZ2llc19maWxlIDwtIGZpbGUucGF0aCgiLi4iLCAiLi4iLCAiZGF0YSIsICJwYnRhLWhpc3RvbG9naWVzLnRzdiIpCgpjb25zZW5zdXNfc2VnX2RmIDwtIHJlYWRfdHN2KGNvbnNlbnN1c19zZWdfZmlsZSkKaGlzdG9sb2dpZXNfZGYgPC0gcmVhZF90c3YoaGlzdG9sb2dpZXNfZmlsZSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgY29sX3R5cGVzID0gY29scygKICAgICAgICAgICAgICAgICAgICAgICAgICAgICBtb2xlY3VsYXJfc3VidHlwZSA9IGNvbF9jaGFyYWN0ZXIoKQogICAgICAgICAgICAgICAgICAgICAgICAgICApKQpgYGAKCiMjIEFkZCBwbG9pZHkgYW5kIHN0YXR1cwoKQXMgbm90ZWQgYWJvdmUsIHRoZSBzZXggY2hyb21vc29tZXMgZG8gbm90IGhhdmUgY29weSBudW1iZXIgaW5mb3JtYXRpb24gaW5jbHVkZWQgaW4gdGhlIGNvbnNlbnN1cyBTRUcgZmlsZSwgYnV0IG90aGVyIGNvcHkgbnVtYmVyIHZhbHVlcyBjYW4gYmUgbWlzc2luZy4KCkZyb20gdGhlIFtjb25zZW5zdXMgU0VHIGZpbGUgbWV0aG9kc10oaHR0cHM6Ly9naXRodWIuY29tL0FsZXhzTGVtb25hZGUvT3BlblBCVEEtYW5hbHlzaXMvdHJlZS8wZTY0MmVmMDYwMmFkMDRiOGQwZmJkODBiZTYyNmQ0Mzc0ZmVlOTI4L2FuYWx5c2VzL2NvcHlfbnVtYmVyX2NvbnNlbnN1c19jYWxsI2NvbnNlbnN1cy1jbnYtY3JlYXRpb24pOiAKCj4gVGhlIGBjb3B5Lm51bWAgY29sdW1uIGlzIHRoZSB3ZWlnaHRlZCBtZWRpYW4gb2YgQ05Wa2l0IHNlZ21lbnQgdmFsdWVzIHdoZXJlIHRoZXkgZXhpc3QsIG9yIENvbnRyb2wtRlJFRUMgdmFsdWVzIGluIHRoZSBhYnNlbmNlIG9mIENOVmtpdCBkYXRhLiBCZWNhdXNlIHNvbWUgc29mdHdhcmUgKG5vdGFibHkgR0lTVElDKSByZXF1aXJlcyBhbGwgc2FtcGxlcyB0byBoYXZlIHRoZSBzYW1lIHJlZ2lvbnMgY2FsbGVkLCB0aGUgY29weSBudW1iZXIgdmFyaWFudHMgZnJvbSBgY252X2NvbnNlbnN1cy50c3ZgIGFyZSBzdXBwbGVtZW50ZW50ZWQgd2l0aCAibmV1dHJhbCIgc2VnbWVudHMgd2hlcmUgbm8gY2FsbCB3YXMgbWFkZS4gVGhlc2UgaW5jbHVkZSBhbGwgbm9uLXZhcmlhbnQgcmVnaW9ucyBwcmVzZW50IGluIGByZWYvY252X2NhbGxhYmxlLmJlZGAgVGhlIG5ldXRyYWwgcmVnaW9ucyBhcmUgYXNzaWduZWQgY29weS5udW0gMiwgZXhjZXB0IG9uIGNoclggYW5kIGNoclksIHdoZXJlIHRoZSBjb3B5IG51bWJlciBpcyBsZWZ0IE5BLgoKRm9yIHNlZ21lbnRzIHdpdGggbWlzc2luZyBgY29weS5udW1gLCB3aGF0IGNocm9tb3NvbWVzIGFyZSB0aGV5IG9uPwoKYGBge3J9CmNvbnNlbnN1c19zZWdfZGYgJT4lCiAgZmlsdGVyKGlzLm5hKGNvcHkubnVtKSkgJT4lCiAgZ3JvdXBfYnkoY2hyb20pICU+JQogIHRhbGx5KCkKYGBgCgpUaGUgbWFqb3JpdHkgYXJlIG9uIHRoZSBzZXggY2hyb21vc29tZXMsIHdoaWNoIGlzIGV4cGVjdGVkLgpXZSBjYW4gcmVtb3ZlIHRoZXNlLCBhcyB3ZSBjYW4gbm90IGFkZCBzdGF0dXMgdXNpbmcgdGhlIHBsb2lkeSBpbmZvcm1hdGlvbiB3aXRob3V0IHRoZSBgY29weS5udW1gIGluZm9ybWF0aW9uLgpXZSBuZWVkIHRvIGFkZCBpbiB0aGUgYHR1bW9yX3Bsb2lkeWAgY29sdW1uIGZyb20gdGhlIGhpc3RvbG9naWVzIGZpbGUuCgpgYGB7cn0KYWRkX3Bsb2lkeV9kZiA8LSBjb25zZW5zdXNfc2VnX2RmICU+JQogIGZpbHRlcighaXMubmEoY29weS5udW0pKSAlPiUKICBpbm5lcl9qb2luKHNlbGVjdChoaXN0b2xvZ2llc19kZiwgCiAgICAgICAgICAgICAgICAgICAgS2lkc19GaXJzdF9CaW9zcGVjaW1lbl9JRCwgCiAgICAgICAgICAgICAgICAgICAgdHVtb3JfcGxvaWR5LAogICAgICAgICAgICAgICAgICAgIGdlcm1saW5lX3NleF9lc3RpbWF0ZSksIAogICAgICAgICAgICAgYnkgPSBjKCJJRCIgPSAiS2lkc19GaXJzdF9CaW9zcGVjaW1lbl9JRCIpKSAlPiUKICBzZWxlY3QoLXR1bW9yX3Bsb2lkeSwgLWdlcm1saW5lX3NleF9lc3RpbWF0ZSwgZXZlcnl0aGluZygpKQpgYGAKCiMjIyBIYW5kbGUgYXV0b3NvbWVzIGZpcnN0CgpgYGB7cn0KYWRkX2F1dG9zb21lc19kZiA8LSBhZGRfcGxvaWR5X2RmICU+JQogICMgeCBhbmQgeSBjaHJvbW9zb21lcyBtdXN0IGJlIGhhbmRsZWQgZGlmZmVyZW50bHkKICBmaWx0ZXIoIShjaHJvbSAlaW4lIGMoImNoclgiLCAiY2hyWSIpKSkgJT4lCiAgbXV0YXRlKHN0YXR1cyA9IGNhc2Vfd2hlbigKICAgICMgd2hlbiB0aGUgY29weSBudW1iZXIgaXMgbGVzcyB0aGFuIGluZmVycmVkIHBsb2lkeSwgbWFyayB0aGlzIGFzIGEgbG9zcwogICAgY29weS5udW0gPCB0dW1vcl9wbG9pZHkgfiAibG9zcyIsCiAgICAjIGlmIGNvcHkgbnVtYmVyIGlzIGhpZ2hlciB0aGFuIHBsb2lkeSwgbWFyayBhcyBhIGdhaW4KICAgIGNvcHkubnVtID4gdHVtb3JfcGxvaWR5IH4gImdhaW4iLAogICAgY29weS5udW0gPT0gdHVtb3JfcGxvaWR5IH4gIm5ldXRyYWwiCiAgKSkKYGBgCgojIyMgSGFuZGxlIHNleCBjaHJvbW9zb21lcwoKYGBge3J9CiMgdGhpcyBsb2dpYyBpcyBjb25zaXN0ZW50IHdpdGggd2hhdCBpcyBvYnNlcnZlZCBpbiB0aGUgY29udHJvbGZyZWVjIGZpbGUKIyBzcGVjaWZpY2FsbHksIGluIHNhbXBsZXMgd2hlcmUgZ2VybWxpbmUgc2V4IGVzdGltYXRlID0gRmVtYWxlLCBYIGNocm9tb3NvbWVzCiMgYXBwZWFyIHRvIGJlIHRyZWF0ZWQgdGhlIHNhbWUgd2F5IGFzIGF1dG9zb21lcwphZGRfeHlfZGYgPC0gYWRkX3Bsb2lkeV9kZiAlPiUKICBmaWx0ZXIoY2hyb20gJWluJSBjKCJjaHJYIiwgImNoclkiKSkgJT4lICAKICBtdXRhdGUoc3RhdHVzID0gY2FzZV93aGVuKAogICAgZ2VybWxpbmVfc2V4X2VzdGltYXRlID09ICJNYWxlIiAmIGNvcHkubnVtIDwgKDAuNSAqIHR1bW9yX3Bsb2lkeSkgfiAibG9zcyIsCiAgICBnZXJtbGluZV9zZXhfZXN0aW1hdGUgPT0gIk1hbGUiICYgY29weS5udW0gPiAoMC41ICogdHVtb3JfcGxvaWR5KSB+ICJnYWluIiwKICAgIGdlcm1saW5lX3NleF9lc3RpbWF0ZSA9PSAiTWFsZSIgJiBjb3B5Lm51bSA9PSAoMC41ICogdHVtb3JfcGxvaWR5KSB+ICJuZXV0cmFsIiwKICAgICMgdGhlcmUgYXJlIHNvbWUgaW5zdGFuY2VzIHdoZXJlIHRoZXJlIGFyZSBjaHJvbW9zb21lIFkgc2VnbWVudHMgYXJlIGluCiAgICAjIHNhbXBsZXMgd2hlcmUgdGhlIGdlcm1saW5lX3NleF9lc3RpbWF0ZSBpcyBGZW1hbGUKICAgIGNocm9tID09ICJjaHJZIiAmIGdlcm1saW5lX3NleF9lc3RpbWF0ZSA9PSAiRmVtYWxlIiAmIGNvcHkubnVtID4gMCB+ICJ1bmtub3duIiwKICAgIGNocm9tID09ICJjaHJZIiAmIGdlcm1saW5lX3NleF9lc3RpbWF0ZSA9PSAiRmVtYWxlIiAmIGNvcHkubnVtID09IDAgfiAibmV1dHJhbCIsCiAgICAjIG1pcnJvcmluZyBjb250cm9sZnJlZWMsIFggdHJlYXRlZCBzYW1lIGFzIGF1dG9zb21lcwogICAgY2hyb20gPT0gImNoclgiICYgZ2VybWxpbmVfc2V4X2VzdGltYXRlID09ICJGZW1hbGUiICYgY29weS5udW0gPCB0dW1vcl9wbG9pZHkgfiAibG9zcyIsCiAgICBjaHJvbSA9PSAiY2hyWCIgJiBnZXJtbGluZV9zZXhfZXN0aW1hdGUgPT0gIkZlbWFsZSIgJiBjb3B5Lm51bSA+IHR1bW9yX3Bsb2lkeSB+ICJnYWluIiwKICAgIGNocm9tID09ICJjaHJYIiAmIGdlcm1saW5lX3NleF9lc3RpbWF0ZSA9PSAiRmVtYWxlIiAmIGNvcHkubnVtID09IHR1bW9yX3Bsb2lkeSB+ICJuZXV0cmFsIgogICkpCgphZGRfc3RhdHVzX2RmIDwtIGJpbmRfcm93cyhhZGRfYXV0b3NvbWVzX2RmLCBhZGRfeHlfZGYpICU+JQogIHNlbGVjdCgtZ2VybWxpbmVfc2V4X2VzdGltYXRlKSAlPiUKICByZW5hbWUoS2lkc19GaXJzdF9CaW9zcGVjaW1lbl9JRCA9IElEKQpgYGAKCiMjIyBEb2VzIGBzZWcubWVhbmAgYWdyZWUgd2l0aCBgc3RhdHVzYD8KCmBgYHtyfQphZGRfc3RhdHVzX2RmICU+JQogIGZpbHRlcighaXMubmEoc2VnLm1lYW4pKSAlPiUKICBnZ3Bsb3QoYWVzKHggPSBzdGF0dXMsIHkgPSBzZWcubWVhbiwgZ3JvdXAgPSBzdGF0dXMpKSArCiAgZ2VvbV9qaXR0ZXIoKQpgYGAKCmBgYHtyfQphZGRfc3RhdHVzX2RmICU+JQogIGZpbHRlcighaXMubmEoc2VnLm1lYW4pKSAlPiUKICBncm91cF9ieShzdGF0dXMpICU+JQogIHN1bW1hcml6ZShtZWFuID0gbWVhbihzZWcubWVhbiksIHNkID0gc2Qoc2VnLm1lYW4pKQpgYGAKCiMjIyBXcml0ZSB0byBmaWxlCgpgYGB7cn0Kb3V0cHV0X2ZpbGUgPC0gZmlsZS5wYXRoKCIuLiIsICIuLiIsICJzY3JhdGNoIiwgImNvbnNlbnN1c19zZWdfd2l0aF9zdGF0dXMudHN2IikKd3JpdGVfdHN2KGFkZF9zdGF0dXNfZGYsIG91dHB1dF9maWxlKQpgYGAK
+
LS0tCnRpdGxlOiAiQWRkIHBsb2lkeSBjb2x1bW4sIHN0YXR1cyB0byBjb25zZW5zdXMgU0VHIGZpbGUiCm91dHB1dDogCiAgaHRtbF9ub3RlYm9vazoKICAgIHRvYzogVFJVRQogICAgdG9jX2Zsb2F0OiBUUlVFCmF1dGhvcjogQ2hhbnRlIEJldGhlbGwgYW5kIEphY2x5biBUYXJvbmkgZm9yIEFMU0YgQ0NETApkYXRlOiAyMDIwCi0tLQoKVGhlIGBwYnRhLWhpc3RvbG9naWVzLnRzdmAgZmlsZSBjb250YWlucyBhIGB0dW1vcl9wbG9pZHlgIGNvbHVtbiwgd2hpY2ggaXMgdHVtb3IgcGxvaWR5IGFzIGluZmVycmVkIGJ5IENvbnRyb2xGcmVlQy4KVGhlIGNvcHkgbnVtYmVyIGluZm9ybWF0aW9uIHNob3VsZCBiZSBpbnRlcnByZXRlZCBpbiB0aGUgbGlnaHQgb2YgdGhpcyBpbmZvcm1hdGlvbiAoc2VlOiBbY3VycmVudCB2ZXJzaW9uIG9mIERhdGEgRm9ybWF0cyBmaWxlXShodHRwczovL2dpdGh1Yi5jb20vQWxleHNMZW1vbmFkZS9PcGVuUEJUQS1hbmFseXNpcy9ibG9iLzBlNjQyZWYwNjAyYWQwNGI4ZDBmYmQ4MGJlNjI2ZDQzNzRmZWU5MjgvZG9jL2RhdGEtZm9ybWF0cy5tZCNhLW5vdGUtb24tcGxvaWR5KSkuCgpUaGlzIG5vdGVib29rIGFkZHMgcGxvaWR5IGluZm9ybWF0aW9uIHRvIHRoZSBjb25zZW5zdXMgU0VHIGZpbGUgYW5kIGFkZHMgYSBzdGF0dXMgY29sdW1uIHRoYXQgZGVmaW5lcyBnYWluIGFuZCBsb3NzIGJyb2FkbHkuCioqTm90ZSoqIHRoYXQgdGhlIGNvbnNlbnN1cyBTRUcgZmlsZSBkb2VzIG5vdCBoYXZlIGNvcHkgbnVtYmVyIGluZm9ybWF0aW9uIGZvciBzZXggY2hyb21vc29tZXMuCgojIyBVc2FnZQoKVGhpcyBub3RlYm9vayBpcyBpbnRlbmRlZCB0byBiZSBydW4gZnJvbSB0aGUgY29tbWFuZCBsaW5lIHdpdGggdGhlIGZvbGxvd2luZyAoYXNzdW1lcyB5b3UgYXJlIGluIHRoZSByb290IGRpcmVjdG9yeSBvZiB0aGUgcmVwb3NpdG9yeSk6CgpgYGAKUnNjcmlwdCAtZSAicm1hcmtkb3duOjpyZW5kZXIoJ2FuYWx5c2VzL2ZvY2FsLWNuLWZpbGUtcHJlcGFyYXRpb24vMDItYWRkLXBsb2lkeS1jb25zZW5zdXMuUm1kJywgY2xlYW4gPSBUUlVFKSIKYGBgCgojIyBTZXQgdXAKCiMjIyBMaWJyYXJpZXMgYW5kIGZ1bmN0aW9ucwoKYGBge3J9CmxpYnJhcnkodGlkeXZlcnNlKQpgYGAKCiMjIyBGaWxlcyBhbmQgZGlyZWN0b3JpZXMKCmBgYHtyfQpzY3JhdGNoX2RpciA8LSBmaWxlLnBhdGgoIi4uIiwgIi4uIiwgInNjcmF0Y2giKQpvdXRwdXRfZGlyIDwtIGZpbGUucGF0aChzY3JhdGNoX2RpciwgImN5dG9iYW5kX3N0YXR1cyIsICJzZWdtZW50cyIpCmlmKCFkaXIuZXhpc3RzKG91dHB1dF9kaXIpKSB7CiAgZGlyLmNyZWF0ZShvdXRwdXRfZGlyLCByZWN1cnNpdmUgPSBUUlVFKQp9CmBgYAoKYGBge3J9CiMgVE9ETzogdGhlIGNvbnNlbnN1cyBTRUcgZmlsZSBpcyBub3QgY3VycmVudGx5IGluIHRoZSBkYXRhIGRvd25sb2FkIC0tIHdoZW4gaXQKIyBnZXRzIGluY2x1ZGVkIHdlIHdpbGwgaGF2ZSB0byBjaGFuZ2UgdGhlIGZpbGUgcGF0aCBoZXJlCmNvbnNlbnN1c19zZWdfZmlsZSA8LSBmaWxlLnBhdGgoIi4uIiwgImNvcHlfbnVtYmVyX2NvbnNlbnN1c19jYWxsIiwgInJlc3VsdHMiLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAicGJ0YS1jbnYtY29uc2Vuc3VzLnNlZy5neiIpCmhpc3RvbG9naWVzX2ZpbGUgPC0gZmlsZS5wYXRoKCIuLiIsICIuLiIsICJkYXRhIiwgInBidGEtaGlzdG9sb2dpZXMudHN2IikKCmNvbnNlbnN1c19zZWdfZGYgPC0gcmVhZF90c3YoY29uc2Vuc3VzX3NlZ19maWxlKQpoaXN0b2xvZ2llc19kZiA8LSByZWFkX3RzdihoaXN0b2xvZ2llc19maWxlLAogICAgICAgICAgICAgICAgICAgICAgICAgICBjb2xfdHlwZXMgPSBjb2xzKAogICAgICAgICAgICAgICAgICAgICAgICAgICAgIG1vbGVjdWxhcl9zdWJ0eXBlID0gY29sX2NoYXJhY3RlcigpCiAgICAgICAgICAgICAgICAgICAgICAgICAgICkpCmBgYAoKIyMgQWRkIHBsb2lkeSBhbmQgc3RhdHVzCgpBcyBub3RlZCBhYm92ZSwgdGhlIHNleCBjaHJvbW9zb21lcyBkbyBub3QgaGF2ZSBjb3B5IG51bWJlciBpbmZvcm1hdGlvbiBpbmNsdWRlZCBpbiB0aGUgY29uc2Vuc3VzIFNFRyBmaWxlLCBidXQgb3RoZXIgY29weSBudW1iZXIgdmFsdWVzIGNhbiBiZSBtaXNzaW5nLgoKRnJvbSB0aGUgW2NvbnNlbnN1cyBTRUcgZmlsZSBtZXRob2RzXShodHRwczovL2dpdGh1Yi5jb20vQWxleHNMZW1vbmFkZS9PcGVuUEJUQS1hbmFseXNpcy90cmVlLzBlNjQyZWYwNjAyYWQwNGI4ZDBmYmQ4MGJlNjI2ZDQzNzRmZWU5MjgvYW5hbHlzZXMvY29weV9udW1iZXJfY29uc2Vuc3VzX2NhbGwjY29uc2Vuc3VzLWNudi1jcmVhdGlvbik6IAoKPiBUaGUgYGNvcHkubnVtYCBjb2x1bW4gaXMgdGhlIHdlaWdodGVkIG1lZGlhbiBvZiBDTlZraXQgc2VnbWVudCB2YWx1ZXMgd2hlcmUgdGhleSBleGlzdCwgb3IgQ29udHJvbC1GUkVFQyB2YWx1ZXMgaW4gdGhlIGFic2VuY2Ugb2YgQ05Wa2l0IGRhdGEuIEJlY2F1c2Ugc29tZSBzb2Z0d2FyZSAobm90YWJseSBHSVNUSUMpIHJlcXVpcmVzIGFsbCBzYW1wbGVzIHRvIGhhdmUgdGhlIHNhbWUgcmVnaW9ucyBjYWxsZWQsIHRoZSBjb3B5IG51bWJlciB2YXJpYW50cyBmcm9tIGBjbnZfY29uc2Vuc3VzLnRzdmAgYXJlIHN1cHBsZW1lbnRlbnRlZCB3aXRoICJuZXV0cmFsIiBzZWdtZW50cyB3aGVyZSBubyBjYWxsIHdhcyBtYWRlLiBUaGVzZSBpbmNsdWRlIGFsbCBub24tdmFyaWFudCByZWdpb25zIHByZXNlbnQgaW4gYHJlZi9jbnZfY2FsbGFibGUuYmVkYCBUaGUgbmV1dHJhbCByZWdpb25zIGFyZSBhc3NpZ25lZCBjb3B5Lm51bSAyLCBleGNlcHQgb24gY2hyWCBhbmQgY2hyWSwgd2hlcmUgdGhlIGNvcHkgbnVtYmVyIGlzIGxlZnQgTkEuCgpGb3Igc2VnbWVudHMgd2l0aCBtaXNzaW5nIGBjb3B5Lm51bWAsIHdoYXQgY2hyb21vc29tZXMgYXJlIHRoZXkgb24/CgpgYGB7cn0KY29uc2Vuc3VzX3NlZ19kZiAlPiUKICBmaWx0ZXIoaXMubmEoY29weS5udW0pKSAlPiUKICBncm91cF9ieShjaHJvbSkgJT4lCiAgdGFsbHkoKQpgYGAKClRoZSBtYWpvcml0eSBhcmUgb24gdGhlIHNleCBjaHJvbW9zb21lcywgd2hpY2ggaXMgZXhwZWN0ZWQuCldlIGNhbiByZW1vdmUgdGhlc2UsIGFzIHdlIGNhbiBub3QgYWRkIHN0YXR1cyB1c2luZyB0aGUgcGxvaWR5IGluZm9ybWF0aW9uIHdpdGhvdXQgdGhlIGBjb3B5Lm51bWAgaW5mb3JtYXRpb24uCldlIG5lZWQgdG8gYWRkIGluIHRoZSBgdHVtb3JfcGxvaWR5YCBjb2x1bW4gZnJvbSB0aGUgaGlzdG9sb2dpZXMgZmlsZS4KCmBgYHtyfQphZGRfcGxvaWR5X2RmIDwtIGNvbnNlbnN1c19zZWdfZGYgJT4lCiAgZmlsdGVyKCFpcy5uYShjb3B5Lm51bSkpICU+JQogIGlubmVyX2pvaW4oc2VsZWN0KGhpc3RvbG9naWVzX2RmLCAKICAgICAgICAgICAgICAgICAgICBLaWRzX0ZpcnN0X0Jpb3NwZWNpbWVuX0lELCAKICAgICAgICAgICAgICAgICAgICB0dW1vcl9wbG9pZHksCiAgICAgICAgICAgICAgICAgICAgZ2VybWxpbmVfc2V4X2VzdGltYXRlKSwgCiAgICAgICAgICAgICBieSA9IGMoIklEIiA9ICJLaWRzX0ZpcnN0X0Jpb3NwZWNpbWVuX0lEIikpICU+JQogIHNlbGVjdCgtdHVtb3JfcGxvaWR5LCAtZ2VybWxpbmVfc2V4X2VzdGltYXRlLCBldmVyeXRoaW5nKCkpCmBgYAoKIyMjIEhhbmRsZSBhdXRvc29tZXMgZmlyc3QKCmBgYHtyfQphZGRfYXV0b3NvbWVzX2RmIDwtIGFkZF9wbG9pZHlfZGYgJT4lCiAgIyB4IGFuZCB5IGNocm9tb3NvbWVzIG11c3QgYmUgaGFuZGxlZCBkaWZmZXJlbnRseQogIGZpbHRlcighKGNocm9tICVpbiUgYygiY2hyWCIsICJjaHJZIikpKSAlPiUKICBtdXRhdGUoc3RhdHVzID0gY2FzZV93aGVuKAogICAgIyB3aGVuIHRoZSBjb3B5IG51bWJlciBpcyBsZXNzIHRoYW4gaW5mZXJyZWQgcGxvaWR5LCBtYXJrIHRoaXMgYXMgYSBsb3NzCiAgICBjb3B5Lm51bSA8IHR1bW9yX3Bsb2lkeSB+ICJsb3NzIiwKICAgICMgaWYgY29weSBudW1iZXIgaXMgaGlnaGVyIHRoYW4gcGxvaWR5LCBtYXJrIGFzIGEgZ2FpbgogICAgY29weS5udW0gPiB0dW1vcl9wbG9pZHkgfiAiZ2FpbiIsCiAgICBjb3B5Lm51bSA9PSB0dW1vcl9wbG9pZHkgfiAibmV1dHJhbCIKICApKQpgYGAKCiMjIyBIYW5kbGUgc2V4IGNocm9tb3NvbWVzCgpgYGB7cn0KIyB0aGlzIGxvZ2ljIGlzIGNvbnNpc3RlbnQgd2l0aCB3aGF0IGlzIG9ic2VydmVkIGluIHRoZSBjb250cm9sZnJlZWMgZmlsZQojIHNwZWNpZmljYWxseSwgaW4gc2FtcGxlcyB3aGVyZSBnZXJtbGluZSBzZXggZXN0aW1hdGUgPSBGZW1hbGUsIFggY2hyb21vc29tZXMKIyBhcHBlYXIgdG8gYmUgdHJlYXRlZCB0aGUgc2FtZSB3YXkgYXMgYXV0b3NvbWVzCmFkZF94eV9kZiA8LSBhZGRfcGxvaWR5X2RmICU+JQogIGZpbHRlcihjaHJvbSAlaW4lIGMoImNoclgiLCAiY2hyWSIpKSAlPiUgIAogIG11dGF0ZShzdGF0dXMgPSBjYXNlX3doZW4oCiAgICBnZXJtbGluZV9zZXhfZXN0aW1hdGUgPT0gIk1hbGUiICYgY29weS5udW0gPCAoMC41ICogdHVtb3JfcGxvaWR5KSB+ICJsb3NzIiwKICAgIGdlcm1saW5lX3NleF9lc3RpbWF0ZSA9PSAiTWFsZSIgJiBjb3B5Lm51bSA+ICgwLjUgKiB0dW1vcl9wbG9pZHkpIH4gImdhaW4iLAogICAgZ2VybWxpbmVfc2V4X2VzdGltYXRlID09ICJNYWxlIiAmIGNvcHkubnVtID09ICgwLjUgKiB0dW1vcl9wbG9pZHkpIH4gIm5ldXRyYWwiLAogICAgIyB0aGVyZSBhcmUgc29tZSBpbnN0YW5jZXMgd2hlcmUgdGhlcmUgYXJlIGNocm9tb3NvbWUgWSBzZWdtZW50cyBhcmUgaW4KICAgICMgc2FtcGxlcyB3aGVyZSB0aGUgZ2VybWxpbmVfc2V4X2VzdGltYXRlIGlzIEZlbWFsZQogICAgY2hyb20gPT0gImNoclkiICYgZ2VybWxpbmVfc2V4X2VzdGltYXRlID09ICJGZW1hbGUiICYgY29weS5udW0gPiAwIH4gInVua25vd24iLAogICAgY2hyb20gPT0gImNoclkiICYgZ2VybWxpbmVfc2V4X2VzdGltYXRlID09ICJGZW1hbGUiICYgY29weS5udW0gPT0gMCB+ICJuZXV0cmFsIiwKICAgICMgbWlycm9yaW5nIGNvbnRyb2xmcmVlYywgWCB0cmVhdGVkIHNhbWUgYXMgYXV0b3NvbWVzCiAgICBjaHJvbSA9PSAiY2hyWCIgJiBnZXJtbGluZV9zZXhfZXN0aW1hdGUgPT0gIkZlbWFsZSIgJiBjb3B5Lm51bSA8IHR1bW9yX3Bsb2lkeSB+ICJsb3NzIiwKICAgIGNocm9tID09ICJjaHJYIiAmIGdlcm1saW5lX3NleF9lc3RpbWF0ZSA9PSAiRmVtYWxlIiAmIGNvcHkubnVtID4gdHVtb3JfcGxvaWR5IH4gImdhaW4iLAogICAgY2hyb20gPT0gImNoclgiICYgZ2VybWxpbmVfc2V4X2VzdGltYXRlID09ICJGZW1hbGUiICYgY29weS5udW0gPT0gdHVtb3JfcGxvaWR5IH4gIm5ldXRyYWwiCiAgKSkKCmFkZF9zdGF0dXNfZGYgPC0gYmluZF9yb3dzKGFkZF9hdXRvc29tZXNfZGYsIGFkZF94eV9kZikgJT4lCiAgc2VsZWN0KC1nZXJtbGluZV9zZXhfZXN0aW1hdGUpICU+JQogIHJlbmFtZShLaWRzX0ZpcnN0X0Jpb3NwZWNpbWVuX0lEID0gSUQpCmBgYAoKIyMjIERvZXMgYHNlZy5tZWFuYCBhZ3JlZSB3aXRoIGBzdGF0dXNgPwoKYGBge3J9CmFkZF9zdGF0dXNfZGYgJT4lCiAgZmlsdGVyKCFpcy5uYShzZWcubWVhbikpICU+JQogIGdncGxvdChhZXMoeCA9IHN0YXR1cywgeSA9IHNlZy5tZWFuLCBncm91cCA9IHN0YXR1cykpICsKICBnZW9tX2ppdHRlcigpCmBgYAoKYGBge3J9CmFkZF9zdGF0dXNfZGYgJT4lCiAgZmlsdGVyKCFpcy5uYShzZWcubWVhbikpICU+JQogIGdyb3VwX2J5KHN0YXR1cykgJT4lCiAgc3VtbWFyaXplKG1lYW4gPSBtZWFuKHNlZy5tZWFuKSwgc2QgPSBzZChzZWcubWVhbikpCmBgYAoKIyMjIFdyaXRlIHRvIGZpbGUKCmBgYHtyfQpvdXRwdXRfZmlsZSA8LSBmaWxlLnBhdGgoIi4uIiwgIi4uIiwgInNjcmF0Y2giLCAiY29uc2Vuc3VzX3NlZ193aXRoX3N0YXR1cy50c3YiKQp3cml0ZV90c3YoYWRkX3N0YXR1c19kZiwgb3V0cHV0X2ZpbGUpCmBgYAoKIyMjIFByZXBhcmUgc2VwYXJhdGUgYmVkIGZpbGVzIGZvciBsb3NzZXMvZ2FpbnMgZm9yIGJlZHRvb2xzIGNvdmVyYWdlIGZ1bmN0aW9uCgpgYGB7cn0KYmVkX3N0YXR1c19kZiA8LSBhZGRfc3RhdHVzX2RmICU+JQogIHNlbGVjdChjaHJvbSwgbG9jLnN0YXJ0LCBsb2MuZW5kLCBldmVyeXRoaW5nKCkpICU+JQogIGdyb3VwX2J5KEtpZHNfRmlyc3RfQmlvc3BlY2ltZW5fSUQpICU+JQogIGFycmFuZ2UoY2hyb20sIGxvYy5zdGFydCwgbG9jLmVuZCkKCmxvc3Nlc19iZWRfc3RhdHVzX2RmIDwtIGJlZF9zdGF0dXNfZGYgJT4lCiAgZmlsdGVyKHN0YXR1cyA9PSAibG9zcyIpCgpnYWluc19iZWRfc3RhdHVzX2RmIDwtIGJlZF9zdGF0dXNfZGYgJT4lCiAgZmlsdGVyKHN0YXR1cyA9PSAiZ2FpbiIpCgojIG1ha2UgbGlzdHMgb2YgZGF0YSBmcmFtZXMgYnkgc2FtcGxlCmJlZF9zdGF0dXNfbGlzdCA8LSBiZWRfc3RhdHVzX2RmICU+JQogIGdyb3VwX3NwbGl0KCkKbmFtZXMoYmVkX3N0YXR1c19saXN0KSA8LSBiZWRfc3RhdHVzX2RmICU+JQogIGdyb3VwX2tleXMoKSAlPiUKICBwdWxsKEtpZHNfRmlyc3RfQmlvc3BlY2ltZW5fSUQpCgpiZWRfbG9zc19saXN0IDwtIGxvc3Nlc19iZWRfc3RhdHVzX2RmICU+JQogIGdyb3VwX3NwbGl0KCkKbmFtZXMoYmVkX2xvc3NfbGlzdCkgPC0gbG9zc2VzX2JlZF9zdGF0dXNfZGYgJT4lCiAgZ3JvdXBfa2V5cygpICU+JQogIHB1bGwoS2lkc19GaXJzdF9CaW9zcGVjaW1lbl9JRCkKCmJlZF9nYWluX2xpc3QgPC0gZ2FpbnNfYmVkX3N0YXR1c19kZiAlPiUKICBncm91cF9zcGxpdCgpCm5hbWVzKGJlZF9nYWluX2xpc3QpIDwtIGdhaW5zX2JlZF9zdGF0dXNfZGYgJT4lCiAgZ3JvdXBfa2V5cygpICU+JQogIHB1bGwoS2lkc19GaXJzdF9CaW9zcGVjaW1lbl9JRCkKYGBgCgojIyMgV3JpdGUgdG8gZmlsZQoKYGBge3J9CnRlbXAgPC0gcHVycnI6OmltYXAoYmVkX3N0YXR1c19saXN0LAogICAgICAgICAgICAgICAgICAgIH4gd3JpdGVfdHN2KC54LAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGZpbGUucGF0aCgKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIG91dHB1dF9kaXIsIHBhc3RlMCgiY29uc2Vuc3VzX2NhbGxhYmxlLiIsIC55LCAiLmJlZCIpCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgKSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBjb2xfbmFtZXMgPSBGQUxTRSkpCgp0ZW1wX2xvc3MgPC0gcHVycnI6OmltYXAoYmVkX2xvc3NfbGlzdCwKICAgICAgICAgICAgICAgICAgICAgICAgIH4gd3JpdGVfdHN2KC54LAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgZmlsZS5wYXRoKAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBvdXRwdXRfZGlyLCBwYXN0ZTAoImNvbnNlbnN1c19sb3NzLiIsIC55LCAiLmJlZCIpCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICApLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgY29sX25hbWVzID0gRkFMU0UpKQoKdGVtcF9nYWluIDwtIHB1cnJyOjppbWFwKGJlZF9nYWluX2xpc3QsCiAgICAgICAgICAgICAgICAgICAgICAgICB+IHdyaXRlX3RzdigueCwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGZpbGUucGF0aCgKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgb3V0cHV0X2RpciwgcGFzdGUwKCJjb25zZW5zdXNfZ2Fpbi4iLCAueSwgIi5iZWQiKQogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgKSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGNvbF9uYW1lcyA9IEZBTFNFKSkKYGBgCgojIyBTZXNzaW9uIEluZm8KCmBgYHtyfQpzZXNzaW9uSW5mbygpCmBgYAoK
diff --git a/analyses/focal-cn-file-preparation/03-add-cytoband-status-consensus.Rmd b/analyses/focal-cn-file-preparation/03-add-cytoband-status-consensus.Rmd new file mode 100644 index 0000000000..1f57b471ac --- /dev/null +++ b/analyses/focal-cn-file-preparation/03-add-cytoband-status-consensus.Rmd @@ -0,0 +1,153 @@ +--- +title: "Add dominant status column to consensus SEG file with cytoband field" +output: + html_notebook: + toc: TRUE + toc_float: TRUE +author: Chante Bethell for ALSF CCDL +date: 2020 +--- + +This notebook adds dominant status information per cytoband to the consensus SEG files prepared in `run-prepare-cn.sh` using the `bedtools coverage` function (these files are stored in the project's scratch directory as bed files). + +## Usage + +This notebook is intended to be run from the command line with the following (assumes you are in the root directory of the repository): + +``` +Rscript -e "rmarkdown::render('analyses/focal-cn-file-preparation/03-add-cytoband-status-consensus.Rmd', clean = TRUE)" +``` + +## Set up + +### Libraries and functions + +```{r} +library(tidyverse) +``` + +### Files and directories + +```{r} +scratch_dir <- file.path("..", "..", "scratch") +``` + +### Read in files + +The files we are reading in here are the bed files prepared in the shell script using `bedtools coverage`. + +```{r message = FALSE} +all_callable_cytoband_status_df <- + read_tsv(file.path(scratch_dir, "intersect_with_cytoband_callable.bed"), + col_names = FALSE) + +loss_cytoband_status_df <- + read_tsv(file.path(scratch_dir, "intersect_with_cytoband_loss.bed"), + col_names = FALSE) + +gain_cytoband_status_df <- + read_tsv(file.path(scratch_dir, "intersect_with_cytoband_gain.bed"), + col_names = FALSE) +``` + +## Wrangle and merge consensus cytoband status data + +### Filter each of the cytoband status data.frames + +```{r} +all_callable_cytoband_status_df <- + all_callable_cytoband_status_df %>% + select( + chr = X1, + cytoband = X4, + band_length = X8, + callable_fraction = X9, + Kids_First_Biospecimen_ID = X10 + ) %>% + filter(!is.na(cytoband)) + +gain_cytoband_status_df <- gain_cytoband_status_df %>% + select( + chr = X1, + cytoband = X4, + gain_fraction = X9, + Kids_First_Biospecimen_ID = X10 + ) %>% + filter(!is.na(cytoband)) + +loss_cytoband_status_df <- loss_cytoband_status_df %>% + select( + chr = X1, + cytoband = X4, + loss_fraction = X9, + Kids_First_Biospecimen_ID = X10 + ) %>% + filter(!is.na(cytoband)) +``` + +### Join all data.frames together + +```{r} +final_df <- all_callable_cytoband_status_df %>% + left_join(gain_cytoband_status_df, + by = c("chr", "cytoband", "Kids_First_Biospecimen_ID")) %>% + inner_join(loss_cytoband_status_df, + by = c("chr", "cytoband", "Kids_First_Biospecimen_ID")) +``` + +### Add `dominant_status` field to final data.frame + +```{r} +# Create a dominant status column +final_df <- final_df %>% + replace_na(list( + gain_fraction = 0, + loss_fraction = 0 + )) %>% + mutate( + dominant_status = case_when( + callable_fraction < 0.5 ~ "uncallable", + gain_fraction / callable_fraction > 0.5 ~ "gain", + loss_fraction / callable_fraction > 0.5 ~ "loss", + (gain_fraction + loss_fraction) / callable_fraction > 0.5 ~ "unstable", + TRUE ~ "neutral" + ) + ) +``` + +### Add chromosome arm column + +```{r} +# Add a column that tells us the position of the p or q and then use this to +# split the cytoband column +final_df <- final_df %>% + mutate( + cytoband_with_arm = paste0(gsub("chr", "", chr), cytoband), + chromosome_arm = gsub("(p|q).*", "\\1", cytoband_with_arm) + ) %>% + select( + Kids_First_Biospecimen_ID, + chr, + cytoband, + dominant_status, + band_length, + everything(), + -cytoband_with_arm, + ) +``` + +### Display and save final table + +```{r} +# Display final table with `uncallable` value added to `dominant_status` column +final_df + +# Write to file +write_tsv(final_df, file.path("results", "consensus_seg_with_ucsc_cytoband_status.tsv.gz")) +``` + +## Session Info + +```{r} +sessionInfo() +``` diff --git a/analyses/focal-cn-file-preparation/03-add-cytoband-status-consensus.nb.html b/analyses/focal-cn-file-preparation/03-add-cytoband-status-consensus.nb.html new file mode 100644 index 0000000000..c02907ea45 --- /dev/null +++ b/analyses/focal-cn-file-preparation/03-add-cytoband-status-consensus.nb.html @@ -0,0 +1,3220 @@ + + + + + + + + + + + + + + + +Add dominant status column to consensus SEG file with cytoband field + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ + + +
+
+
+
+
+ +
+ + + + + + + + +

This notebook adds dominant status information per cytoband to the consensus SEG files prepared in run-prepare-cn.sh using the bedtools coverage function (these files are stored in the project’s scratch directory as bed files).

+
+

Usage

+

This notebook is intended to be run from the command line with the following (assumes you are in the root directory of the repository):

+
Rscript -e "rmarkdown::render('analyses/focal-cn-file-preparation/03-add-cytoband-status-consensus.Rmd', clean = TRUE)"
+
+
+

Set up

+
+

Libraries and functions

+ + + +
library(tidyverse)
+ + +
── Attaching packages ────────────────────────────────────────────── tidyverse 1.2.1 ──
+ + +
✔ ggplot2 3.2.0     ✔ purrr   0.3.2
+✔ tibble  2.1.3     ✔ dplyr   0.8.3
+✔ tidyr   0.8.3     ✔ stringr 1.4.0
+✔ readr   1.3.1     ✔ forcats 0.4.0
+ + +
── Conflicts ───────────────────────────────────────────────── tidyverse_conflicts() ──
+✖ dplyr::filter() masks stats::filter()
+✖ dplyr::lag()    masks stats::lag()
+ + + +
+
+

Files and directories

+ + + +
scratch_dir <- file.path("..", "..", "scratch")
+ + + +
+
+

Read in files

+

The files we are reading in here are the bed files prepared in the shell script using bedtools coverage.

+ + + +
all_callable_cytoband_status_df <-
+  read_tsv(file.path(scratch_dir, "intersect_with_cytoband_callable.bed"),
+                  col_names = FALSE)
+
+loss_cytoband_status_df <-
+  read_tsv(file.path(scratch_dir, "intersect_with_cytoband_loss.bed"),
+                  col_names = FALSE)
+
+gain_cytoband_status_df <-
+  read_tsv(file.path(scratch_dir, "intersect_with_cytoband_gain.bed"),
+                  col_names = FALSE)
+ + + +
+
+
+

Wrangle and merge consensus cytoband status data

+
+

Filter each of the cytoband status data.frames

+ + + +
all_callable_cytoband_status_df <-
+  all_callable_cytoband_status_df %>%
+  select(
+    chr = X1,
+    cytoband = X4,
+    band_length = X8,
+    callable_fraction = X9,
+    Kids_First_Biospecimen_ID = X10
+  ) %>%
+  filter(!is.na(cytoband))
+
+gain_cytoband_status_df <- gain_cytoband_status_df %>%
+  select(
+    chr = X1,
+    cytoband = X4,
+    gain_fraction = X9,
+    Kids_First_Biospecimen_ID = X10
+  ) %>%
+  filter(!is.na(cytoband))
+
+loss_cytoband_status_df <- loss_cytoband_status_df %>%
+  select(
+    chr = X1,
+    cytoband = X4,
+    loss_fraction = X9,
+    Kids_First_Biospecimen_ID = X10
+  ) %>%
+  filter(!is.na(cytoband))
+ + + +
+
+

Join all data.frames together

+ + + +
final_df <- all_callable_cytoband_status_df %>%
+  left_join(gain_cytoband_status_df,
+            by = c("chr", "cytoband", "Kids_First_Biospecimen_ID")) %>%
+  inner_join(loss_cytoband_status_df,
+            by = c("chr", "cytoband", "Kids_First_Biospecimen_ID"))
+ + + +
+
+

Add dominant_status field to final data.frame

+ + + +
# Create a dominant status column
+final_df <- final_df %>%
+  replace_na(list(
+    gain_fraction = 0,
+    loss_fraction = 0
+  )) %>%
+  mutate(
+    dominant_status = case_when(
+      callable_fraction < 0.5 ~ "uncallable",
+      gain_fraction / callable_fraction > 0.5 ~ "gain",
+      loss_fraction / callable_fraction > 0.5 ~ "loss",
+      (gain_fraction + loss_fraction) / callable_fraction > 0.5 ~ "unstable",
+      TRUE ~ "neutral"
+    )
+  )
+ + + +
+
+

Add chromosome arm column

+ + + +
# Add a column that tells us the position of the p or q and then use this to
+# split the cytoband column
+final_df <- final_df %>%
+  mutate(
+    cytoband_with_arm = paste0(gsub("chr", "", chr), cytoband),
+    chromosome_arm = gsub("(p|q).*", "\\1", cytoband_with_arm)
+  ) %>%
+  select(
+    Kids_First_Biospecimen_ID,
+    chr,
+    cytoband,
+    dominant_status,
+    band_length,
+    everything(),
+    -cytoband_with_arm,
+  )
+ + + +
+
+

Display and save final table

+ + + +
# Display final table with `uncallable` value added to `dominant_status` column
+final_df
+ +
+ +
+ +
# Write to file
+write_tsv(final_df, file.path("results", "consensus_seg_with_ucsc_cytoband_status.tsv.gz"))
+ + + +
+
+
+

Session Info

+ + + +
sessionInfo()
+ + +
R version 3.6.0 (2019-04-26)
+Platform: x86_64-pc-linux-gnu (64-bit)
+Running under: Debian GNU/Linux 9 (stretch)
+
+Matrix products: default
+BLAS/LAPACK: /usr/lib/libopenblasp-r0.2.19.so
+
+locale:
+ [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
+ [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
+ [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=C             
+ [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
+ [9] LC_ADDRESS=C               LC_TELEPHONE=C            
+[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
+
+attached base packages:
+[1] stats     graphics  grDevices utils     datasets  methods   base     
+
+other attached packages:
+[1] forcats_0.4.0   stringr_1.4.0   dplyr_0.8.3     purrr_0.3.2    
+[5] readr_1.3.1     tidyr_0.8.3     tibble_2.1.3    ggplot2_3.2.0  
+[9] tidyverse_1.2.1
+
+loaded via a namespace (and not attached):
+ [1] Rcpp_1.0.1       cellranger_1.1.0 pillar_1.4.2     compiler_3.6.0  
+ [5] base64enc_0.1-3  tools_3.6.0      digest_0.6.20    lubridate_1.7.4 
+ [9] jsonlite_1.6     evaluate_0.14    nlme_3.1-140     gtable_0.3.0    
+[13] lattice_0.20-38  pkgconfig_2.0.2  rlang_0.4.0      cli_1.1.0       
+[17] rstudioapi_0.10  yaml_2.2.0       haven_2.1.1      xfun_0.8        
+[21] withr_2.1.2      xml2_1.2.0       httr_1.4.0       knitr_1.23      
+[25] generics_0.0.2   hms_0.4.2        grid_3.6.0       tidyselect_0.2.5
+[29] glue_1.3.1       R6_2.4.0         readxl_1.3.1     rmarkdown_1.13  
+[33] modelr_0.1.4     magrittr_1.5     backports_1.1.4  scales_1.0.0    
+[37] htmltools_0.3.6  rvest_0.3.4      assertthat_0.2.1 colorspace_1.4-1
+[41] stringi_1.4.3    lazyeval_0.2.2   munsell_0.5.0    broom_0.5.2     
+[45] crayon_1.3.4    
+ + +
+ +
LS0tCnRpdGxlOiAiQWRkIGRvbWluYW50IHN0YXR1cyBjb2x1bW4gdG8gY29uc2Vuc3VzIFNFRyBmaWxlIHdpdGggY3l0b2JhbmQgZmllbGQiCm91dHB1dDogCiAgaHRtbF9ub3RlYm9vazoKICAgIHRvYzogVFJVRQogICAgdG9jX2Zsb2F0OiBUUlVFCmF1dGhvcjogQ2hhbnRlIEJldGhlbGwgZm9yIEFMU0YgQ0NETApkYXRlOiAyMDIwCi0tLQoKVGhpcyBub3RlYm9vayBhZGRzIGRvbWluYW50IHN0YXR1cyBpbmZvcm1hdGlvbiBwZXIgY3l0b2JhbmQgdG8gdGhlIGNvbnNlbnN1cyBTRUcgZmlsZXMgcHJlcGFyZWQgaW4gYHJ1bi1wcmVwYXJlLWNuLnNoYCB1c2luZyB0aGUgYGJlZHRvb2xzIGNvdmVyYWdlYCBmdW5jdGlvbiAodGhlc2UgZmlsZXMgYXJlIHN0b3JlZCBpbiB0aGUgcHJvamVjdCdzIHNjcmF0Y2ggZGlyZWN0b3J5IGFzIGJlZCBmaWxlcykuCgojIyBVc2FnZQoKVGhpcyBub3RlYm9vayBpcyBpbnRlbmRlZCB0byBiZSBydW4gZnJvbSB0aGUgY29tbWFuZCBsaW5lIHdpdGggdGhlIGZvbGxvd2luZyAoYXNzdW1lcyB5b3UgYXJlIGluIHRoZSByb290IGRpcmVjdG9yeSBvZiB0aGUgcmVwb3NpdG9yeSk6CgpgYGAKUnNjcmlwdCAtZSAicm1hcmtkb3duOjpyZW5kZXIoJ2FuYWx5c2VzL2ZvY2FsLWNuLWZpbGUtcHJlcGFyYXRpb24vMDMtYWRkLWN5dG9iYW5kLXN0YXR1cy1jb25zZW5zdXMuUm1kJywgY2xlYW4gPSBUUlVFKSIKYGBgCgojIyBTZXQgdXAKCiMjIyBMaWJyYXJpZXMgYW5kIGZ1bmN0aW9ucwoKYGBge3J9CmxpYnJhcnkodGlkeXZlcnNlKQpgYGAKCiMjIyBGaWxlcyBhbmQgZGlyZWN0b3JpZXMKCmBgYHtyfQpzY3JhdGNoX2RpciA8LSBmaWxlLnBhdGgoIi4uIiwgIi4uIiwgInNjcmF0Y2giKQpgYGAKCiMjIyBSZWFkIGluIGZpbGVzIAoKVGhlIGZpbGVzIHdlIGFyZSByZWFkaW5nIGluIGhlcmUgYXJlIHRoZSBiZWQgZmlsZXMgcHJlcGFyZWQgaW4gdGhlIHNoZWxsIHNjcmlwdCB1c2luZyBgYmVkdG9vbHMgY292ZXJhZ2VgLgoKYGBge3IgbWVzc2FnZSA9IEZBTFNFfQphbGxfY2FsbGFibGVfY3l0b2JhbmRfc3RhdHVzX2RmIDwtCiAgcmVhZF90c3YoZmlsZS5wYXRoKHNjcmF0Y2hfZGlyLCAiaW50ZXJzZWN0X3dpdGhfY3l0b2JhbmRfY2FsbGFibGUuYmVkIiksCiAgICAgICAgICAgICAgICAgIGNvbF9uYW1lcyA9IEZBTFNFKQoKbG9zc19jeXRvYmFuZF9zdGF0dXNfZGYgPC0KICByZWFkX3RzdihmaWxlLnBhdGgoc2NyYXRjaF9kaXIsICJpbnRlcnNlY3Rfd2l0aF9jeXRvYmFuZF9sb3NzLmJlZCIpLAogICAgICAgICAgICAgICAgICBjb2xfbmFtZXMgPSBGQUxTRSkKCmdhaW5fY3l0b2JhbmRfc3RhdHVzX2RmIDwtCiAgcmVhZF90c3YoZmlsZS5wYXRoKHNjcmF0Y2hfZGlyLCAiaW50ZXJzZWN0X3dpdGhfY3l0b2JhbmRfZ2Fpbi5iZWQiKSwKICAgICAgICAgICAgICAgICAgY29sX25hbWVzID0gRkFMU0UpCmBgYAoKIyMgV3JhbmdsZSBhbmQgbWVyZ2UgY29uc2Vuc3VzIGN5dG9iYW5kIHN0YXR1cyBkYXRhCgojIyMgRmlsdGVyIGVhY2ggb2YgdGhlIGN5dG9iYW5kIHN0YXR1cyBkYXRhLmZyYW1lcwoKYGBge3J9CmFsbF9jYWxsYWJsZV9jeXRvYmFuZF9zdGF0dXNfZGYgPC0KICBhbGxfY2FsbGFibGVfY3l0b2JhbmRfc3RhdHVzX2RmICU+JQogIHNlbGVjdCgKICAgIGNociA9IFgxLAogICAgY3l0b2JhbmQgPSBYNCwKICAgIGJhbmRfbGVuZ3RoID0gWDgsCiAgICBjYWxsYWJsZV9mcmFjdGlvbiA9IFg5LAogICAgS2lkc19GaXJzdF9CaW9zcGVjaW1lbl9JRCA9IFgxMAogICkgJT4lCiAgZmlsdGVyKCFpcy5uYShjeXRvYmFuZCkpCgpnYWluX2N5dG9iYW5kX3N0YXR1c19kZiA8LSBnYWluX2N5dG9iYW5kX3N0YXR1c19kZiAlPiUKICBzZWxlY3QoCiAgICBjaHIgPSBYMSwKICAgIGN5dG9iYW5kID0gWDQsCiAgICBnYWluX2ZyYWN0aW9uID0gWDksCiAgICBLaWRzX0ZpcnN0X0Jpb3NwZWNpbWVuX0lEID0gWDEwCiAgKSAlPiUKICBmaWx0ZXIoIWlzLm5hKGN5dG9iYW5kKSkKCmxvc3NfY3l0b2JhbmRfc3RhdHVzX2RmIDwtIGxvc3NfY3l0b2JhbmRfc3RhdHVzX2RmICU+JQogIHNlbGVjdCgKICAgIGNociA9IFgxLAogICAgY3l0b2JhbmQgPSBYNCwKICAgIGxvc3NfZnJhY3Rpb24gPSBYOSwKICAgIEtpZHNfRmlyc3RfQmlvc3BlY2ltZW5fSUQgPSBYMTAKICApICU+JQogIGZpbHRlcighaXMubmEoY3l0b2JhbmQpKQpgYGAKCiMjIyBKb2luIGFsbCBkYXRhLmZyYW1lcyB0b2dldGhlcgoKYGBge3J9CmZpbmFsX2RmIDwtIGFsbF9jYWxsYWJsZV9jeXRvYmFuZF9zdGF0dXNfZGYgJT4lCiAgbGVmdF9qb2luKGdhaW5fY3l0b2JhbmRfc3RhdHVzX2RmLAogICAgICAgICAgICBieSA9IGMoImNociIsICJjeXRvYmFuZCIsICJLaWRzX0ZpcnN0X0Jpb3NwZWNpbWVuX0lEIikpICU+JQogIGlubmVyX2pvaW4obG9zc19jeXRvYmFuZF9zdGF0dXNfZGYsCiAgICAgICAgICAgIGJ5ID0gYygiY2hyIiwgImN5dG9iYW5kIiwgIktpZHNfRmlyc3RfQmlvc3BlY2ltZW5fSUQiKSkKYGBgCgojIyMgQWRkIGBkb21pbmFudF9zdGF0dXNgIGZpZWxkIHRvIGZpbmFsIGRhdGEuZnJhbWUKCmBgYHtyfQojIENyZWF0ZSBhIGRvbWluYW50IHN0YXR1cyBjb2x1bW4KZmluYWxfZGYgPC0gZmluYWxfZGYgJT4lCiAgcmVwbGFjZV9uYShsaXN0KAogICAgZ2Fpbl9mcmFjdGlvbiA9IDAsCiAgICBsb3NzX2ZyYWN0aW9uID0gMAogICkpICU+JQogIG11dGF0ZSgKICAgIGRvbWluYW50X3N0YXR1cyA9IGNhc2Vfd2hlbigKICAgICAgY2FsbGFibGVfZnJhY3Rpb24gPCAwLjUgfiAidW5jYWxsYWJsZSIsCiAgICAgIGdhaW5fZnJhY3Rpb24gLyBjYWxsYWJsZV9mcmFjdGlvbiA+IDAuNSB+ICJnYWluIiwKICAgICAgbG9zc19mcmFjdGlvbiAvIGNhbGxhYmxlX2ZyYWN0aW9uID4gMC41IH4gImxvc3MiLAogICAgICAoZ2Fpbl9mcmFjdGlvbiArIGxvc3NfZnJhY3Rpb24pIC8gY2FsbGFibGVfZnJhY3Rpb24gPiAwLjUgfiAidW5zdGFibGUiLAogICAgICBUUlVFIH4gIm5ldXRyYWwiCiAgICApCiAgKQpgYGAKCiMjIyBBZGQgY2hyb21vc29tZSBhcm0gY29sdW1uCgpgYGB7cn0KIyBBZGQgYSBjb2x1bW4gdGhhdCB0ZWxscyB1cyB0aGUgcG9zaXRpb24gb2YgdGhlIHAgb3IgcSBhbmQgdGhlbiB1c2UgdGhpcyB0bwojIHNwbGl0IHRoZSBjeXRvYmFuZCBjb2x1bW4KZmluYWxfZGYgPC0gZmluYWxfZGYgJT4lCiAgbXV0YXRlKAogICAgY3l0b2JhbmRfd2l0aF9hcm0gPSBwYXN0ZTAoZ3N1YigiY2hyIiwgIiIsIGNociksIGN5dG9iYW5kKSwKICAgIGNocm9tb3NvbWVfYXJtID0gZ3N1YigiKHB8cSkuKiIsICJcXDEiLCBjeXRvYmFuZF93aXRoX2FybSkKICApICU+JQogIHNlbGVjdCgKICAgIEtpZHNfRmlyc3RfQmlvc3BlY2ltZW5fSUQsCiAgICBjaHIsCiAgICBjeXRvYmFuZCwKICAgIGRvbWluYW50X3N0YXR1cywKICAgIGJhbmRfbGVuZ3RoLAogICAgZXZlcnl0aGluZygpLAogICAgLWN5dG9iYW5kX3dpdGhfYXJtLAogICkKYGBgCgojIyMgRGlzcGxheSBhbmQgc2F2ZSBmaW5hbCB0YWJsZQoKYGBge3J9CiMgRGlzcGxheSBmaW5hbCB0YWJsZSB3aXRoIGB1bmNhbGxhYmxlYCB2YWx1ZSBhZGRlZCB0byBgZG9taW5hbnRfc3RhdHVzYCBjb2x1bW4KZmluYWxfZGYKCiMgV3JpdGUgdG8gZmlsZQp3cml0ZV90c3YoZmluYWxfZGYsIGZpbGUucGF0aCgicmVzdWx0cyIsICJjb25zZW5zdXNfc2VnX3dpdGhfdWNzY19jeXRvYmFuZF9zdGF0dXMudHN2Lmd6IikpCmBgYAoKIyMgU2Vzc2lvbiBJbmZvCgpgYGB7cn0Kc2Vzc2lvbkluZm8oKQpgYGAK
+ + +
+
+ +
+ + + + + + + + diff --git a/analyses/focal-cn-file-preparation/03-prepare-cn-file.R b/analyses/focal-cn-file-preparation/04-prepare-cn-file.R similarity index 99% rename from analyses/focal-cn-file-preparation/03-prepare-cn-file.R rename to analyses/focal-cn-file-preparation/04-prepare-cn-file.R index 78c5942d94..d3e6429c5e 100644 --- a/analyses/focal-cn-file-preparation/03-prepare-cn-file.R +++ b/analyses/focal-cn-file-preparation/04-prepare-cn-file.R @@ -11,7 +11,7 @@ # This script is intended to be run via the command line. # This example assumes it is being run from the root of the repository. # -# Rscript --vanilla analyses/oncoprint-landscape/03-prepare-cn-file.R \ +# Rscript --vanilla analyses/focal-cn-file-preparation/03-prepare-cn-file.R \ # --cnv_file data/pbta-cnv-controlfreec.tsv.gz \ # --gtf_file data/gencode.v27.primary_assembly.annotation.gtf.gz \ # --metadata data/pbta-histologies.tsv \ diff --git a/analyses/focal-cn-file-preparation/README.md b/analyses/focal-cn-file-preparation/README.md index e07991bcf9..33d288e43e 100644 --- a/analyses/focal-cn-file-preparation/README.md +++ b/analyses/focal-cn-file-preparation/README.md @@ -1,6 +1,6 @@ ## Focal copy number file preparation -**Module authors:** Chante Bethell ([@cbethell](https://github.com/cbethell)) and Jaclyn Taroni ([@jaclyn-taroni](https://github.com/jaclyn-taroni)) +**Module authors:** Chante Bethell ([@cbethell](https://github.com/cbethell)), Joshua Shapiro ([@jashapiro](https://github.com/jashapiro)), and Jaclyn Taroni ([@jaclyn-taroni](https://github.com/jaclyn-taroni)) The copy number data from OpenPBTA are provided as ranges or segments. The purpose of this module is to map from those ranges to gene identifiers for consumption by downstream analyses (e.g., OncoPrint plotting). @@ -12,6 +12,8 @@ To run this analysis _only on consensus SEG file_, use the following (from the r ``` bash analyses/focal-cn-file-preparation/run-prepare-cn.sh ``` +**Note**: The `run-bedtools.snakemake` script is implemented in `run-prepare-cn.sh` to run the bedtools coverage steps between the UCSC cytoband file and the samples in the copy number files produced in `02-add-ploidy-consensus.Rmd`. +This script currently takes a while to run, and therefore slows down the processing speed of the main shell script `run-prepare-cn.sh`. Running the following from the root directory of the repository will run the steps for the original copy number call files (CNVkit and ControlFreeC) in addition to the consensus SEG file: @@ -28,9 +30,14 @@ RUN_ORIGINAL=1 bash analyses/focal-cn-file-preparation/run-prepare-cn.sh * `02-add-ploidy-consensus.Rmd` - This is very similar to the CNVkit file prep (`01-add-ploidy-cnvkit.Rmd`). However, there are instances in the consensus SEG file where `copy.num` is `NA` which are removed. -See the notebook for more information. +See the notebook for more information. This notebook also prepares lists of copy number bed files by sample for use in the implementation of bedtools coverage in `run-bedtools.snakemake`. -* `03-prepare-cn-file.R` - This script performs the ranges to annotation mapping using the GENCODE v27 GTF included via the data download step; it takes the ControlFreeC file or a SEG (e.g., CNVkit, consensus SEG) file prepared with `01-add-ploidy-cnvkit.Rmd` and `02-add-ploidy-cnvkit.Rmd` as input. +* `03-add-cytoband-status-consensus.Rmd` - This notebook reads in the bedtools coverage output files and defines the dominant copy number status for each chromosome arm. The output of this notebook is a table with the following columns: + + | `Kids_First_Biospecimen_ID` | chr | cytoband | dominant_status | band_length | callable_fraction | gain_fraction | loss_fraction | chromosome_arm | + |----------------|--------|-------------|--------|---------|----------|-------------|---------|---------------| + +* `04-prepare-cn-file.R` - This script performs the ranges to annotation mapping using the GENCODE v27 GTF included via the data download step; it takes the ControlFreeC file or a SEG (e.g., CNVkit, consensus SEG) file prepared with `01-add-ploidy-cnvkit.Rmd` and `02-add-ploidy-cnvkit.Rmd` as input. **The mapping is limited to _exons_.** Mapping to cytobands is performed with the [`org.Hs.eg.db`](https://doi.org/doi:10.18129/B9.bioc.org.Hs.eg.db) package. A table with the following columns is returned: @@ -66,9 +73,43 @@ focal-cn-file-preparation ├── 01-add-ploidy-cnvkit.nb.html ├── 02-add-ploidy-consensus.Rmd ├── 02-add-ploidy-consensus.nb.html -├── 03-prepare-cn-file.R +├── 03-add-cytoband-status-consensus.Rmd +├── 03-add-cytoband-status-consensus.nb.html +├── 04-prepare-cn-file.R ├── README.md +├── annotation_files +│   └── txdb_from_gencode.v27.gtf.db ├── display-plots.md +├── gistic-results +│   └── pbta-cnv-cnvkit-gistic +│   ├── D.cap1.5.mat +│   ├── all_data_by_genes.txt +│   ├── all_lesions.conf_90.txt +│   ├── all_thresholded.by_genes.txt +│   ├── amp_genes.conf_90.txt +│   ├── amp_qplot.pdf +│   ├── amp_qplot.png +│   ├── broad_data_by_genes.txt +│   ├── broad_gistic_plot.pdf +│   ├── broad_significance_results.txt +│   ├── broad_values_by_arm.txt +│   ├── del_genes.conf_90.txt +│   ├── del_qplot.pdf +│   ├── del_qplot.png +│   ├── focal_dat.0.98.mat +│   ├── focal_data_by_genes.txt +│   ├── freqarms_vs_ngenes.pdf +│   ├── gistic_inputs.mat +│   ├── peak_regs.mat +│   ├── perm_ads.mat +│   ├── raw_copy_number.pdf +│   ├── raw_copy_number.png +│   ├── regions_track.conf_90.bed +│   ├── sample_cutoffs.txt +│   ├── sample_seg_counts.txt +│   ├── scores.0.98.mat +│   ├── scores.gistic +│   └── wide_peak_regs.mat ├── plots │   ├── cnvkit_annotated_cn_autosomes_polya_loss_cor_plot.png │   ├── cnvkit_annotated_cn_autosomes_polya_stacked_plot.png @@ -111,9 +152,11 @@ focal-cn-file-preparation │   ├── cnvkit_annotated_cn_x_and_y.tsv.gz │   ├── consensus_seg_annotated_cn_autosomes.tsv.gz │   ├── consensus_seg_annotated_cn_x_and_y.tsv.gz +│   ├── consensus_seg_with_ucsc_cytoband.tsv.gz │   ├── controlfreec_annotated_cn_autosomes.tsv.gz │   └── controlfreec_annotated_cn_x_and_y.tsv.gz ├── rna-expression-validation.R +├── run-bedtools.snakemake ├── run-prepare-cn.sh └── util └── rna-expression-functions.R diff --git a/analyses/focal-cn-file-preparation/results/consensus_seg_with_ucsc_cytoband_status.tsv.gz b/analyses/focal-cn-file-preparation/results/consensus_seg_with_ucsc_cytoband_status.tsv.gz new file mode 100644 index 0000000000..2a0175c872 Binary files /dev/null and b/analyses/focal-cn-file-preparation/results/consensus_seg_with_ucsc_cytoband_status.tsv.gz differ diff --git a/analyses/focal-cn-file-preparation/run-bedtools.snakemake b/analyses/focal-cn-file-preparation/run-bedtools.snakemake new file mode 100644 index 0000000000..05910be9d3 --- /dev/null +++ b/analyses/focal-cn-file-preparation/run-bedtools.snakemake @@ -0,0 +1,75 @@ + + +# Set up +scratch_dir = "../../scratch" +cytoband_dir = scratch_dir + "/cytoband_status" + +# Get the lists of samples to analyze +samples, = glob_wildcards(cytoband_dir + "/segments/consensus_callable.{sample}.bed") +gain_samples, = glob_wildcards(cytoband_dir + "/segments/consensus_gain.{sample}.bed") +loss_samples, = glob_wildcards(cytoband_dir + "/segments/consensus_loss.{sample}.bed") + +wildcard_constraints: + set = "callable|gain|loss", + +# Set output files +rule targets: + input: + scratch_dir + "/intersect_with_cytoband_callable.bed", + scratch_dir + "/intersect_with_cytoband_gain.bed", + scratch_dir + "/intersect_with_cytoband_loss.bed" + +##### +# Join each set of coverage stats separately +rule cat_coverage: + input: + expand(cytoband_dir + "/coverage/consensus_callable.{sample}.coverage.bed", sample = samples) + output: + scratch_dir + "/intersect_with_cytoband_callable.bed" + shell: + "cat {input} > {output}" + +rule cat_gains: + input: + expand(cytoband_dir + "/coverage/consensus_gain.{sample}.coverage.bed", sample = gain_samples) + output: + scratch_dir + "/intersect_with_cytoband_gain.bed" + shell: + "cat {input} > {output}" + +rule cat_losses: + input: + expand(cytoband_dir + "/coverage/consensus_loss.{sample}.coverage.bed", sample = loss_samples) + output: + scratch_dir + "/intersect_with_cytoband_loss.bed" + shell: + "cat {input} > {output}" + +#### +# Calculate coverage for each sample +rule bed_coverage: + input: + bed = cytoband_dir + "/segments/consensus_{set}.{sample}.bed", + bands = scratch_dir + "/ucsc_cytoband.bed" + output: + cytoband_dir + "/coverage/consensus_{set}.{sample}.coverage.bed" + shell: + "bedtools coverage " + " -a {input.bands}" + " -b {input.bed}" + " -sorted" + " | sed 's/$/\t{wildcards.sample}/' " # add a column to label the sample. + " > {output}" + + + +rule get_cytobands: + output: + scratch_dir + "/ucsc_cytoband.bed" + params: + url = "http://hgdownload.cse.ucsc.edu/goldenpath/hg38/database/cytoBand.txt.gz" + shell: + "wget -O - {params.url}" + " | gunzip -c" + " | grep '^chr[0-9XY]\+\s' " # filter to only canonical autosomes + " > {output}" diff --git a/analyses/focal-cn-file-preparation/run-prepare-cn.sh b/analyses/focal-cn-file-preparation/run-prepare-cn.sh index d9598393c6..9939d9f2b0 100644 --- a/analyses/focal-cn-file-preparation/run-prepare-cn.sh +++ b/analyses/focal-cn-file-preparation/run-prepare-cn.sh @@ -18,6 +18,7 @@ cd "$script_directory" || exit scratch_dir=../../scratch data_dir=../../data +results_dir=../../analyses/focal-cn-file-preparation/results histologies_file=${data_dir}/pbta-histologies.tsv gtf_file=${data_dir}/gencode.v27.primary_assembly.annotation.gtf.gz goi_file=../../analyses/oncoprint-landscape/driver-lists/brain-goi-list-long.txt @@ -26,66 +27,75 @@ independent_specimens_file=${data_dir}/independent-specimens.wgswxs.primary.tsv # Prep the consensus SEG file data Rscript --vanilla -e "rmarkdown::render('02-add-ploidy-consensus.Rmd', clean = TRUE)" +# Run snakemake script implementing `bedtools coverage` for each sample bed file in +# `scratch/cytoband_status` -- these files are generated in +# `02-add-ploidy-consensus.Rmd` +# currently runs 10 jobs in parallel, which should be fine for most implementations +snakemake -j 10 --snakefile run-bedtools.snakemake + +# Determine the dominant status for each chromosome arm and compare with GISTIC's arm status calls +Rscript --vanilla -e "rmarkdown::render('03-add-cytoband-status-consensus.Rmd', clean = TRUE)" + # Run annotation step for consensus file -Rscript --vanilla 03-prepare-cn-file.R \ +Rscript --vanilla 04-prepare-cn-file.R \ --cnv_file ${scratch_dir}/consensus_seg_with_status.tsv \ --gtf_file $gtf_file \ --metadata $histologies_file \ --filename_lead "consensus_seg_annotated_cn" \ --seg -libraryStrategies=("polya" "stranded") -chromosomesType=("autosomes" "x_and_y") -for strategy in ${libraryStrategies[@]}; do - - for chromosome_type in ${chromosomesType[@]}; do - - Rscript --vanilla rna-expression-validation.R \ - --annotated_cnv_file results/consensus_seg_annotated_cn_${chromosome_type}.tsv.gz \ - --expression_file ${data_dir}/pbta-gene-expression-rsem-fpkm-collapsed.${strategy}.rds \ - --independent_specimens_file $independent_specimens_file \ - --metadata $histologies_file \ - --goi_list $goi_file \ - --filename_lead "consensus_seg_annotated_cn"_${chromosome_type}_${strategy} - done -done - -# if we want to process the CNV data from the original callers -# (e.g., CNVkit, ControlFreeC) -if [ "$RUN_ORIGINAL" -gt "0" ]; then - - # Prep the CNVkit data - Rscript --vanilla -e "rmarkdown::render('01-add-ploidy-cnvkit.Rmd', clean = TRUE)" - - # Run annotation step for CNVkit - Rscript --vanilla 03-prepare-cn-file.R \ - --cnv_file ${scratch_dir}/cnvkit_with_status.tsv \ - --gtf_file $gtf_file \ - --metadata $histologies_file \ - --filename_lead "cnvkit_annotated_cn" \ - --seg - - # Run annotation step for ControlFreeC - Rscript --vanilla 03-prepare-cn-file.R \ - --cnv_file ${data_dir}/pbta-cnv-controlfreec.tsv.gz \ - --gtf_file $gtf_file \ - --metadata $histologies_file \ - --filename_lead "controlfreec_annotated_cn" \ - --controlfreec - - filenameLead=("cnvkit_annotated_cn" "controlfreec_annotated_cn") - for filename in ${filenameLead[@]}; do - for strategy in ${libraryStrategies[@]}; do - for chromosome_type in ${chromosomesType[@]}; do - Rscript --vanilla rna-expression-validation.R \ - --annotated_cnv_file results/${filename}_${chromosome_type}.tsv.gz \ - --expression_file ${data_dir}/pbta-gene-expression-rsem-fpkm-collapsed.${strategy}.rds \ - --independent_specimens_file $independent_specimens_file \ - --metadata $histologies_file \ - --goi_list $goi_file \ - --filename_lead ${filename}_${chromosome_type}_${strategy} - done - done - done - -fi +# libraryStrategies=("polya" "stranded") +# chromosomesType=("autosomes" "x_and_y") +# for strategy in ${libraryStrategies[@]}; do +# +# for chromosome_type in ${chromosomesType[@]}; do +# +# Rscript --vanilla rna-expression-validation.R \ +# --annotated_cnv_file results/consensus_seg_annotated_cn_${chromosome_type}.tsv.gz \ +# --expression_file ${data_dir}/pbta-gene-expression-rsem-fpkm-collapsed.${strategy}.rds \ +# --independent_specimens_file $independent_specimens_file \ +# --metadata $histologies_file \ +# --goi_list $goi_file \ +# --filename_lead "consensus_seg_annotated_cn"_${chromosome_type}_${strategy} +# done +# done +# +# # if we want to process the CNV data from the original callers +# # (e.g., CNVkit, ControlFreeC) +# if [ "$RUN_ORIGINAL" -gt "0" ]; then +# +# # Prep the CNVkit data +# Rscript --vanilla -e "rmarkdown::render('01-add-ploidy-cnvkit.Rmd', clean = TRUE)" +# +# # Run annotation step for CNVkit +# Rscript --vanilla 04-prepare-cn-file.R \ +# --cnv_file ${scratch_dir}/cnvkit_with_status.tsv \ +# --gtf_file $gtf_file \ +# --metadata $histologies_file \ +# --filename_lead "cnvkit_annotated_cn" \ +# --seg +# +# # Run annotation step for ControlFreeC +# Rscript --vanilla 04-prepare-cn-file.R \ +# --cnv_file ${data_dir}/pbta-cnv-controlfreec.tsv.gz \ +# --gtf_file $gtf_file \ +# --metadata $histologies_file \ +# --filename_lead "controlfreec_annotated_cn" \ +# --controlfreec +# +# filenameLead=("cnvkit_annotated_cn" "controlfreec_annotated_cn") +# for filename in ${filenameLead[@]}; do +# for strategy in ${libraryStrategies[@]}; do +# for chromosome_type in ${chromosomesType[@]}; do +# Rscript --vanilla rna-expression-validation.R \ +# --annotated_cnv_file results/${filename}_${chromosome_type}.tsv.gz \ +# --expression_file ${data_dir}/pbta-gene-expression-rsem-fpkm-collapsed.${strategy}.rds \ +# --independent_specimens_file $independent_specimens_file \ +# --metadata $histologies_file \ +# --goi_list $goi_file \ +# --filename_lead ${filename}_${chromosome_type}_${strategy} +# done +# done +# done +# +# fi