Prepare analysis

Load libraries

library(tidyverse)
library(rhdf5)
#library(DropletUtils) # install but not load

Set filepaths

knitr::opts_knit$set(root.dir = rprojroot::find_rstudio_root_file(),
                     fig.width=15,
                     digit=5,
                     scipen=8)
options(digits=5, 
        scipen=8,
        future.globals.maxSize = +Inf)
project_dir <- rprojroot::find_rstudio_root_file()
if(is.null(project_dir)){
  project_dir <- getwd()
  warning(sprintf("No rstudio project root file  found. 
                  Setting project directory to current workflow.Rmd file location: %s. 
                  Override if needed.",
                  project_dir))
 
}
message(sprintf("Project directory: %s",
                project_dir))
Project directory: /home/rfarouni/Documents/index_hopping
code_dir <- file.path(project_dir, "code")
source(file.path(code_dir, "1_create_joined_counts_table.R"))
datasets_names <- c("hiseq4000_nonplexed", "hiseq4000_plexed")
names(datasets_names) <- datasets_names
datasets_names
  hiseq4000_nonplexed      hiseq4000_plexed 
"hiseq4000_nonplexed"    "hiseq4000_plexed" 
validation_output_dir <- file.path(project_dir, "data", "hiseq4000_validation")
dir.create(validation_output_dir)
'/home/rfarouni/Documents/index_hopping/data/hiseq4000_validation' already exists

Define functions

get_read_counts <- function(dataset_name) {
  data_dir <- file.path(project_dir, "data", dataset_name)
  output_dir <- file.path(data_dir, "output")
  input_dir <- file.path(data_dir, "input")
  read_counts_filepath <- file.path(
    output_dir,
    paste0(
      dataset_name,
      "_read_counts.rds"
    )
  )
  read_counts <- 
    create_joined_counts(
      input_dir,
      read_counts_filepath
  )
  return(read_counts)
}

Explore and prepare data

Load data

read_counts <- map(datasets_names, get_read_counts)
data_fj <-
  full_join(read_counts$hiseq4000_nonplexed %>%
    select(-outcome),
  read_counts$hiseq4000_plexed,
  by = c("cell", "gene", "umi"),
  suffix = c("_nonplexed", "_plexed")
  )

Anonymize

create keys

old_key <- unique(data_fj$gene)
new_key <-
  sample.int(
    n = length(old_key),
    size = length(old_key)
  )
names(new_key) <- old_key
new_key <-
  enframe(new_key,
    name = "gene",
    value = "gene_new"
  )

recode

data_fj <-
  left_join(data_fj, new_key, by = "gene") %>%
  mutate(
    gene = NULL,
    gene = gene_new
  ) %>%
  select(cell, gene, umi, everything()) %>%
  select(-gene_new) %>%
  mutate_if(is.double, as.integer) %>%
  set_names(c("cell", "gene", "umi", "s1_nonplexed", "s2_nonplexed", "s1_plexed", "s2_plexed", "outcome"))
data_fj
ABCDEFGHIJ0123456789
cell
<chr>
gene
<int>
umi
<int>
s1_nonplexed
<int>
s2_nonplexed
<int>
s1_plexed
<int>
s2_plexed
<int>
outcome
<chr>
AAACCTGAGAAACCTA999667523810NANANA
AAACCTGAGAAACCTA1310881055910NANANA
AAACCTGAGAAACCTA1528348648710NANANA
AAACCTGAGAAACCTA1412115900010101,0
AAACCTGAGAAACCTA1412111550820NANANA
AAACCTGAGAAACCTA810884239230202,0
AAACCTGAGAAACCTA2605100066010NANANA
AAACCTGAGAAACCTA26051850010NANANA
AAACCTGAGAAACCTA260513854610101,0
AAACCTGAGAAACCTA1114860601940NANANA
Next
123456
...
100
Previous
1-10 of 16,547,728 rows

Save data

write_tsv(data_fj,
          file.path(validation_output_dir, "hiseq4000_joined_datatable_plexed_nonplexed.txt"))

Retain molecules that are observed in both datasets

inner join

data <-  
  data_fj %>%
  drop_na()
data
ABCDEFGHIJ0123456789
 
 
cell
<chr>
gene
<int>
umi
<int>
s1_nonplexed
<int>
s2_nonplexed
<int>
s1_plexed
<int>
s2_plexed
<int>
outcome
<chr>
4AAACCTGAGAAACCTA1412115900010101,0
6AAACCTGAGAAACCTA810884239230202,0
9AAACCTGAGAAACCTA260513854610101,0
11AAACCTGAGAAACCTA272585581250202,0
12AAACCTGAGAAACCTA221157712420202,0
13AAACCTGAGAAACCTA533925023310303,0
17AAACCTGAGAAACCTA1961930254430303,0
19AAACCTGAGAAACCTA2088464417930202,0
20AAACCTGAGAAACCTA2614104337820202,0
21AAACCTGAGAAACCTA417075174932401851185,1
Next
123456
...
100
Previous
1-10 of 9,252,147 rows

Create goundtruth labels

data <-
  data %>% 
  mutate(label= case_when(
      s1_nonplexed != 0 & s2_nonplexed != 0 & s1_plexed!=0 & s2_plexed!=0~ "r,r",
      s1_nonplexed == 0 & s1_plexed == 0 ~ "0,r",
      s2_nonplexed == 0 & s2_plexed == 0 ~ "r,0",
      s1_nonplexed == 0 & s1_plexed != 0 & s2_plexed == 0 ~ "f,0",
      s1_nonplexed == 0 & s1_plexed != 0 & s2_plexed != 0 ~ "f,r",
      s2_nonplexed == 0 & s2_plexed != 0 & s1_plexed == 0 ~ "0,f",
      s2_nonplexed == 0 & s2_plexed != 0 & s1_plexed != 0 ~ "r,f",
      s1_nonplexed != 0 & s2_nonplexed != 0 & s1_plexed==0 & s2_plexed!=0~ "0,r",
      s1_nonplexed != 0 & s2_nonplexed != 0 & s1_plexed!=0 & s2_plexed==0~ "r,0",
      TRUE                      ~ "NA"
    )) 
data
ABCDEFGHIJ0123456789
cell
<chr>
gene
<int>
umi
<int>
s1_nonplexed
<int>
s2_nonplexed
<int>
s1_plexed
<int>
s2_plexed
<int>
outcome
<chr>
label
<chr>
AAACCTGAGAAACCTA1412115900010101,0r,0
AAACCTGAGAAACCTA810884239230202,0r,0
AAACCTGAGAAACCTA260513854610101,0r,0
AAACCTGAGAAACCTA272585581250202,0r,0
AAACCTGAGAAACCTA221157712420202,0r,0
AAACCTGAGAAACCTA533925023310303,0r,0
AAACCTGAGAAACCTA1961930254430303,0r,0
AAACCTGAGAAACCTA2088464417930202,0r,0
AAACCTGAGAAACCTA2614104337820202,0r,0
AAACCTGAGAAACCTA417075174932401851185,1r,f
Next
123456
...
100
Previous
1-10 of 9,252,147 rows

Tally labels

label_tally <-
  data %>%
  group_by(label) %>%
  tally() %>%
  add_row(label = "TOTAL", n=sum(.$n))
label_tally
ABCDEFGHIJ0123456789
label
<chr>
n
<int>
0,f3484
0,r6637616
f,06511
f,r215147
r,02129080
r,f260257
r,r52
TOTAL9252147
8 rows

There are only 52 (real, real) colliding CUGs out of 9252147 so assumption I is validated. The collision rate is 0.00001.

Filter out colliding CUGs

data <-
  data %>% 
  filter(label!="r,r")

Save data

write_tsv(data, file.path(validation_output_dir,"hiseq4000_inner_joined_with_labels.txt"))

Summary statistics

number of reads (in millions)

data_fj %>%  
  ungroup() %>%
  summarize_at(vars(matches("^s")), list(~ sum(./10^6,  na.rm = TRUE)))
ABCDEFGHIJ0123456789
s1_nonplexed
<dbl>
s2_nonplexed
<dbl>
s1_plexed
<dbl>
s2_plexed
<dbl>
184.77184.1792.777.909
1 row
data  %>% 
  summarize_at(vars(matches("^s")), list(~ sum(./10^6)))
ABCDEFGHIJ0123456789
s1_nonplexed
<dbl>
s2_nonplexed
<dbl>
s1_plexed
<dbl>
s2_plexed
<dbl>
182.01176.6791.88476.294
1 row

number of molecules (in millions)

data_fj %>%  
  mutate_at(vars(matches("^s")), list(~ as.integer(.!=0)))  %>%
  ungroup() %>%
  summarize_at(vars(matches("^s")), list(~ sum(./10^6, na.rm = TRUE)))
ABCDEFGHIJ0123456789
s1_nonplexed
<dbl>
s2_nonplexed
<dbl>
s1_plexed
<dbl>
s2_plexed
<dbl>
3.906610.8373.22098.3142
1 row
 data  %>% 
  mutate_at(vars(matches("^s")), list(~ as.integer(.!=0)))  %>%
  ungroup()  %>% 
  summarize_at(vars(matches("^s")), list(~ sum(./10^6, na.rm = TRUE)))
ABCDEFGHIJ0123456789
s1_nonplexed
<dbl>
s2_nonplexed
<dbl>
s1_plexed
<dbl>
s2_plexed
<dbl>
2.39296.85932.6117.1165
1 row

Examine concordance between samples

Read counts with same CUG (cell-umi-gene) label

p1 <- ggplot(data, aes(x = s1_nonplexed,
                  y = s1_plexed)) 
p1 + geom_hex(bins = 500)+
  geom_abline(slope=.5,intercept=0)

p2 <- ggplot(data, 
            aes(x = s2_nonplexed,
                  y = s2_plexed)) 
p2 + 
  geom_hex(bins = 500)+
  geom_abline(slope=.5,intercept=0)

 #   geom_point()

Gene expression read counts with same cell-gene label

data_reads_cell_gene <- 
    data %>%  
  group_by(cell, gene) %>%
  summarize_at(vars(matches("^s")),
               list(~ sum(.)))
 ggplot(data_reads_cell_gene,
            aes(x = s1_nonplexed,
                  y = s1_plexed))  + 
  geom_hex(bins = 400) +
  geom_abline(slope=1,intercept=0)+
  scale_x_log10() +
   scale_y_log10() 

ggplot(data_reads_cell_gene,
            aes(x = s2_nonplexed,
                  y = s2_plexed))  + 
  geom_hex(bins = 400) +
    geom_abline(slope=1,intercept=0)+
  scale_x_log10() +
   scale_y_log10() 

NA
ggplot(data_reads_cell_gene,
            aes(x = s1_nonplexed,
                  y = s2_nonplexed))  + 
  geom_hex(bins = 400) +
  geom_abline(slope=1,intercept=0)+
  scale_x_log10() +
   scale_y_log10() 

NA
ggplot(data_reads_cell_gene,
            aes(x = s1_plexed,
                  y = s2_plexed))  + 
  geom_hex(bins = 400) +
  geom_abline(slope=1,intercept=0)+
  scale_x_log10() +
   scale_y_log10() 

NA
ggplot(data_reads_cell_gene,
            aes(x = s1_plexed,
                  y = s2_nonplexed))  + 
  geom_hex(bins = 400) +
  geom_abline(slope=1,intercept=0)+
  scale_x_log10() +
   scale_y_log10() 

NA

## Gene expression UMI counts with same cell-gene label

data_molec_cell_gene <- 
    data %>%  
  mutate_at(vars(matches("^s")), 
            list(~ as.integer(.!=0)))  %>%  
  group_by(cell, gene) %>%
  summarize_at(vars(matches("^s")),
               list(~ sum(.)))
 ggplot(data_molec_cell_gene,
            aes(x = s1_nonplexed,
                  y = s1_plexed))  + 
  geom_hex(bins = 400) +
  geom_abline(slope=1,intercept=0)+
  scale_x_log10() +
   scale_y_log10() 

ggplot(data_molec_cell_gene,
            aes(x = s2_nonplexed,
                  y = s2_plexed))  + 
  geom_hex(bins = 400) +
  geom_abline(slope=1,intercept=0)+
  scale_x_log10() +
   scale_y_log10() 

NA
ggplot(data_molec_cell_gene,
            aes(x = s1_nonplexed,
                  y = s2_nonplexed))  + 
  geom_hex(bins = 400) +
  geom_abline(slope=1,intercept=0)+
  scale_x_log10() +
   scale_y_log10() 

NA

A visual demonstration of the effects of index hopping

ggplot(data_molec_cell_gene,
            aes(x = s1_plexed,
                  y = s2_plexed))  + 
  geom_hex(bins = 400) +
  geom_abline(slope=1,intercept=0)+
  scale_x_log10() +
   scale_y_log10() 

NA
ggplot(data_molec_cell_gene,
            aes(x = s1_plexed,
                  y = s2_nonplexed))  + 
  geom_hex(bins = 400) +
  geom_abline(slope=1,intercept=0)+
  scale_x_log10() +
   scale_y_log10() 

NA
ggplot(data_molec_cell_gene,
            aes(x = s2_plexed,
                  y = s1_nonplexed))  + 
  geom_hex(bins = 400) +
  geom_abline(slope=1,intercept=0)+
  scale_x_log10() +
   scale_y_log10() 

NA

## Gene expression UMI counts (nonhopping) with same cell-gene label

Here we retain molecules for which the reads didn’t hop

data_molec_cell_gene_nonhopping <- 
    data %>%   
  filter(label%in% c("r,0","0,r"))%>%
  mutate_at(vars(matches("^s")), 
            list(~ as.integer(.!=0)))  %>%  
  group_by(cell, gene) %>%
  summarize_at(vars(matches("^s")),
               list(~ sum(.)))
 ggplot(data_molec_cell_gene_nonhopping,
            aes(x = s1_nonplexed,
                  y = s1_plexed))  + 
  geom_hex(bins = 400) +
  geom_abline(slope=1,intercept=0)+
  scale_x_log10() +
   scale_y_log10() 

ggplot(data_molec_cell_gene_nonhopping,
            aes(x = s2_nonplexed,
                  y = s2_plexed))  + 
  geom_hex(bins = 400) +
  geom_abline(slope=1,intercept=0)+
  scale_x_log10() +
   scale_y_log10() 

NA
ggplot(data_molec_cell_gene_nonhopping,
            aes(x = s1_nonplexed,
                  y = s2_nonplexed))  + 
  geom_hex(bins = 400) +
  geom_abline(slope=1,intercept=0)+
  scale_x_log10() +
   scale_y_log10() 

NA
ggplot(data_molec_cell_gene_nonhopping,
            aes(x = s1_plexed,
                  y = s2_plexed))  + 
  geom_hex(bins = 400) +
  geom_abline(slope=1,intercept=0)+
  scale_x_log10() +
   scale_y_log10() 

NA
ggplot(data_molec_cell_gene_nonhopping,
            aes(x = s1_plexed,
                  y = s2_nonplexed))  + 
  geom_hex(bins = 400) +
  geom_abline(slope=1,intercept=0)+
  scale_x_log10() +
   scale_y_log10() 

NA
sessionInfo()
R version 3.6.0 (2019-04-26)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.2 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8   
 [6] LC_MESSAGES=en_US.UTF-8    LC_PAPER=en_US.UTF-8       LC_NAME=C                  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] hexbin_1.27.3   rhdf5_2.28.0    forcats_0.4.0   stringr_1.4.0   dplyr_0.8.1     purrr_0.3.2     readr_1.3.1     tidyr_0.8.3     tibble_2.1.2   
[10] ggplot2_3.1.1   tidyverse_1.2.1

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.1       cellranger_1.1.0 pillar_1.4.1     compiler_3.6.0   plyr_1.8.4       tools_3.6.0      digest_0.6.19    jsonlite_1.6    
 [9] lubridate_1.7.4  nlme_3.1-140     gtable_0.3.0     lattice_0.20-38  pkgconfig_2.0.2  rlang_0.3.4      cli_1.1.0        rstudioapi_0.10 
[17] yaml_2.2.0       haven_2.1.0      xfun_0.7         withr_2.1.2      xml2_1.2.0       httr_1.4.0       knitr_1.23       generics_0.0.2  
[25] hms_0.4.2        rprojroot_1.3-2  grid_3.6.0       tidyselect_0.2.5 glue_1.3.1       R6_2.4.0         readxl_1.3.1     Rhdf5lib_1.6.0  
[33] modelr_0.1.4     magrittr_1.5     backports_1.1.4  scales_1.0.0     rvest_0.3.4      assertthat_0.2.1 colorspace_1.4-1 labeling_0.3    
[41] stringi_1.4.3    lazyeval_0.2.2   munsell_0.5.0    broom_0.5.2      crayon_1.3.4
LS0tCnRpdGxlOiAiUGhhbnRvbSBQdXJnZSIKc3VidGl0bGU6ICJWYWxpZGF0aW9uIEFuYWx5c2lzOiBQYXJ0IEkiCmF1dGhvcjogCi0gbmFtZTogUmljayBGYXJvdW5pCiAgYWZmaWxpYXRpb246CiAgLSAmY3J1ayBHw6lub21lIFF1w6liZWMgSW5ub3ZhdGlvbiBDZW50cmUsIE1jR2lsbCBVbml2ZXJzaXR5LCBNb250cmVhbCwgQ2FuYWRhCmRhdGU6ICdgciBmb3JtYXQoU3lzLkRhdGUoKSwgIiVZLSVCLSVkIilgJwpvdXRwdXQ6CiAgaHRtbF9ub3RlYm9vazoKICAgIGRmX3ByaW50OiBwYWdlZAogICAgY29kZV9mb2xkaW5nOiBzaG93CiAgICB0b2M6IG5vCiAgICB0b2NfZmxvYXQ6IAogICAgICBjb2xsYXBzZWQ6IGZhbHNlCiAgICAgIHNtb290aF9zY3JvbGw6IGZhbHNlCi0tLQoKIyBQcmVwYXJlIGFuYWx5c2lzCgojIyMgTG9hZCBsaWJyYXJpZXMKCmBgYHtyIG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0V9CmxpYnJhcnkodGlkeXZlcnNlKQpsaWJyYXJ5KHJoZGY1KQojbGlicmFyeShEcm9wbGV0VXRpbHMpICMgaW5zdGFsbCBidXQgZG8gbm90IGxvYWQKYGBgCgoKIyMjIFNldCBmaWxlcGF0aHMgCgoKYGBge3Igc2V0dXB9CmtuaXRyOjpvcHRzX2tuaXQkc2V0KHJvb3QuZGlyID0gcnByb2pyb290OjpmaW5kX3JzdHVkaW9fcm9vdF9maWxlKCksCiAgICAgICAgICAgICAgICAgICAgIGZpZy53aWR0aD0xNSwKICAgICAgICAgICAgICAgICAgICAgZGlnaXQ9NSwKICAgICAgICAgICAgICAgICAgICAgc2NpcGVuPTgpCm9wdGlvbnMoZGlnaXRzPTUsIAogICAgICAgIHNjaXBlbj04LAogICAgICAgIGZ1dHVyZS5nbG9iYWxzLm1heFNpemUgPSArSW5mKQpgYGAKCgoKYGBge3J9CnByb2plY3RfZGlyIDwtIHJwcm9qcm9vdDo6ZmluZF9yc3R1ZGlvX3Jvb3RfZmlsZSgpCgppZihpcy5udWxsKHByb2plY3RfZGlyKSl7CiAgcHJvamVjdF9kaXIgPC0gZ2V0d2QoKQogIHdhcm5pbmcoc3ByaW50ZigiTm8gcnN0dWRpbyBwcm9qZWN0IHJvb3QgZmlsZSAgZm91bmQuIAogICAgICAgICAgICAgICAgICBTZXR0aW5nIHByb2plY3QgZGlyZWN0b3J5IHRvIGN1cnJlbnQgd29ya2Zsb3cuUm1kIGZpbGUgbG9jYXRpb246ICVzLiAKICAgICAgICAgICAgICAgICAgT3ZlcnJpZGUgaWYgbmVlZGVkLiIsCiAgICAgICAgICAgICAgICAgIHByb2plY3RfZGlyKSkKIAp9Cm1lc3NhZ2Uoc3ByaW50ZigiUHJvamVjdCBkaXJlY3Rvcnk6ICVzIiwKICAgICAgICAgICAgICAgIHByb2plY3RfZGlyKSkKYGBgCgoKCmBgYHtyfQpjb2RlX2RpciA8LSBmaWxlLnBhdGgocHJvamVjdF9kaXIsICJjb2RlIikKc291cmNlKGZpbGUucGF0aChjb2RlX2RpciwgIjFfY3JlYXRlX2pvaW5lZF9jb3VudHNfdGFibGUuUiIpKQoKZGF0YXNldHNfbmFtZXMgPC0gYygiaGlzZXE0MDAwX25vbnBsZXhlZCIsICJoaXNlcTQwMDBfcGxleGVkIikKbmFtZXMoZGF0YXNldHNfbmFtZXMpIDwtIGRhdGFzZXRzX25hbWVzCmRhdGFzZXRzX25hbWVzCmBgYAoKCmBgYHtyfQp2YWxpZGF0aW9uX291dHB1dF9kaXIgPC0gZmlsZS5wYXRoKHByb2plY3RfZGlyLCAiZGF0YSIsICJoaXNlcTQwMDBfdmFsaWRhdGlvbiIpCmRpci5jcmVhdGUodmFsaWRhdGlvbl9vdXRwdXRfZGlyKQpgYGAKCiMjIyBEZWZpbmUgZnVuY3Rpb25zCgpgYGB7cn0KZ2V0X3JlYWRfY291bnRzIDwtIGZ1bmN0aW9uKGRhdGFzZXRfbmFtZSkgewogIGRhdGFfZGlyIDwtIGZpbGUucGF0aChwcm9qZWN0X2RpciwgImRhdGEiLCBkYXRhc2V0X25hbWUpCiAgb3V0cHV0X2RpciA8LSBmaWxlLnBhdGgoZGF0YV9kaXIsICJvdXRwdXQiKQogIGlucHV0X2RpciA8LSBmaWxlLnBhdGgoZGF0YV9kaXIsICJpbnB1dCIpCgogIHJlYWRfY291bnRzX2ZpbGVwYXRoIDwtIGZpbGUucGF0aCgKICAgIG91dHB1dF9kaXIsCiAgICBwYXN0ZTAoCiAgICAgIGRhdGFzZXRfbmFtZSwKICAgICAgIl9yZWFkX2NvdW50cy5yZHMiCiAgICApCiAgKQoKCiAgcmVhZF9jb3VudHMgPC0gCiAgICBjcmVhdGVfam9pbmVkX2NvdW50cygKICAgICAgaW5wdXRfZGlyLAogICAgICByZWFkX2NvdW50c19maWxlcGF0aAogICkKCiAgcmV0dXJuKHJlYWRfY291bnRzKQp9CmBgYAoKCiMgRXhwbG9yZSBhbmQgcHJlcGFyZSBkYXRhCgojIyBMb2FkIGRhdGEKCmBgYHtyfQpyZWFkX2NvdW50cyA8LSBtYXAoZGF0YXNldHNfbmFtZXMsIGdldF9yZWFkX2NvdW50cykKYGBgCgoKYGBge3J9CmRhdGFfZmogPC0KICBmdWxsX2pvaW4ocmVhZF9jb3VudHMkaGlzZXE0MDAwX25vbnBsZXhlZCAlPiUKICAgIHNlbGVjdCgtb3V0Y29tZSksCiAgcmVhZF9jb3VudHMkaGlzZXE0MDAwX3BsZXhlZCwKICBieSA9IGMoImNlbGwiLCAiZ2VuZSIsICJ1bWkiKSwKICBzdWZmaXggPSBjKCJfbm9ucGxleGVkIiwgIl9wbGV4ZWQiKQogICkKYGBgCgoKCiMjIEFub255bWl6ZQoKY3JlYXRlIGtleXMKCmBgYHtyfQpvbGRfa2V5IDwtIHVuaXF1ZShkYXRhX2ZqJGdlbmUpCgpuZXdfa2V5IDwtCiAgc2FtcGxlLmludCgKICAgIG4gPSBsZW5ndGgob2xkX2tleSksCiAgICBzaXplID0gbGVuZ3RoKG9sZF9rZXkpCiAgKQpuYW1lcyhuZXdfa2V5KSA8LSBvbGRfa2V5Cm5ld19rZXkgPC0KICBlbmZyYW1lKG5ld19rZXksCiAgICBuYW1lID0gImdlbmUiLAogICAgdmFsdWUgPSAiZ2VuZV9uZXciCiAgKQpgYGAKCnJlY29kZQoKYGBge3J9CmRhdGFfZmogPC0KICBsZWZ0X2pvaW4oZGF0YV9maiwgbmV3X2tleSwgYnkgPSAiZ2VuZSIpICU+JQogIG11dGF0ZSgKICAgIGdlbmUgPSBOVUxMLAogICAgZ2VuZSA9IGdlbmVfbmV3CiAgKSAlPiUKICBzZWxlY3QoY2VsbCwgZ2VuZSwgdW1pLCBldmVyeXRoaW5nKCkpICU+JQogIHNlbGVjdCgtZ2VuZV9uZXcpICU+JQogIG11dGF0ZV9pZihpcy5kb3VibGUsIGFzLmludGVnZXIpICU+JQogIHNldF9uYW1lcyhjKCJjZWxsIiwgImdlbmUiLCAidW1pIiwgInMxX25vbnBsZXhlZCIsICJzMl9ub25wbGV4ZWQiLCAiczFfcGxleGVkIiwgInMyX3BsZXhlZCIsICJvdXRjb21lIikpCgpkYXRhX2ZqCmBgYAoKCgojIyMgU2F2ZSBkYXRhCgpgYGB7cn0Kd3JpdGVfdHN2KGRhdGFfZmosCiAgICAgICAgICBmaWxlLnBhdGgodmFsaWRhdGlvbl9vdXRwdXRfZGlyLCAiaGlzZXE0MDAwX2pvaW5lZF9kYXRhdGFibGVfcGxleGVkX25vbnBsZXhlZC50eHQiKSkKYGBgCgoKIyMjIFJldGFpbiBtb2xlY3VsZXMgdGhhdCBhcmUgb2JzZXJ2ZWQgaW4gYm90aCBkYXRhc2V0cwoKaW5uZXIgam9pbgoKYGBge3J9CmRhdGEgPC0gIAogIGRhdGFfZmogJT4lCiAgZHJvcF9uYSgpCmRhdGEKYGBgCgojIyMgQ3JlYXRlIGdvdW5kdHJ1dGggbGFiZWxzCgpgYGB7cn0KZGF0YSA8LQogIGRhdGEgJT4lIAogIG11dGF0ZShsYWJlbD0gY2FzZV93aGVuKAogICAgICBzMV9ub25wbGV4ZWQgIT0gMCAmIHMyX25vbnBsZXhlZCAhPSAwICYgczFfcGxleGVkIT0wICYgczJfcGxleGVkIT0wfiAicixyIiwKICAgICAgczFfbm9ucGxleGVkID09IDAgJiBzMV9wbGV4ZWQgPT0gMCB+ICIwLHIiLAogICAgICBzMl9ub25wbGV4ZWQgPT0gMCAmIHMyX3BsZXhlZCA9PSAwIH4gInIsMCIsCiAgICAgIHMxX25vbnBsZXhlZCA9PSAwICYgczFfcGxleGVkICE9IDAgJiBzMl9wbGV4ZWQgPT0gMCB+ICJmLDAiLAogICAgICBzMV9ub25wbGV4ZWQgPT0gMCAmIHMxX3BsZXhlZCAhPSAwICYgczJfcGxleGVkICE9IDAgfiAiZixyIiwKICAgICAgczJfbm9ucGxleGVkID09IDAgJiBzMl9wbGV4ZWQgIT0gMCAmIHMxX3BsZXhlZCA9PSAwIH4gIjAsZiIsCiAgICAgIHMyX25vbnBsZXhlZCA9PSAwICYgczJfcGxleGVkICE9IDAgJiBzMV9wbGV4ZWQgIT0gMCB+ICJyLGYiLAogICAgICBzMV9ub25wbGV4ZWQgIT0gMCAmIHMyX25vbnBsZXhlZCAhPSAwICYgczFfcGxleGVkPT0wICYgczJfcGxleGVkIT0wfiAiMCxyIiwKICAgICAgczFfbm9ucGxleGVkICE9IDAgJiBzMl9ub25wbGV4ZWQgIT0gMCAmIHMxX3BsZXhlZCE9MCAmIHMyX3BsZXhlZD09MH4gInIsMCIsCiAgICAgIFRSVUUgICAgICAgICAgICAgICAgICAgICAgfiAiTkEiCiAgICApKSAKZGF0YQpgYGAKCiMjIyBUYWxseSBsYWJlbHMKCgpgYGB7cn0KbGFiZWxfdGFsbHkgPC0KICBkYXRhICU+JQogIGdyb3VwX2J5KGxhYmVsKSAlPiUKICB0YWxseSgpICU+JQogIGFkZF9yb3cobGFiZWwgPSAiVE9UQUwiLCBuPXN1bSguJG4pKQpsYWJlbF90YWxseQpgYGAKClRoZXJlIGFyZSBvbmx5IGByIGxhYmVsX3RhbGx5ICU+JSBmaWx0ZXIobGFiZWw9PSJyLHIiKSAlPiUgcHVsbChuKWAgKHJlYWwsIHJlYWwpIGNvbGxpZGluZyBDVUdzIG91dCBvZiBgciBsYWJlbF90YWxseSAlPiUgZmlsdGVyKGxhYmVsPT0iVE9UQUwiKSAlPiUgcHVsbChuKWAgc28gYXNzdW1wdGlvbiBJIGlzIHZhbGlkYXRlZC4gVGhlIGNvbGxpc2lvbiByYXRlIGlzIGByIChsYWJlbF90YWxseSAlPiUgZmlsdGVyKGxhYmVsPT0icixyIikgJT4lIHB1bGwobikpLyhsYWJlbF90YWxseSAlPiUgZmlsdGVyKGxhYmVsPT0iVE9UQUwiKSAlPiUgcHVsbChuKSlgLgoKCiMjIyBGaWx0ZXIgb3V0IGNvbGxpZGluZyBDVUdzCgpgYGB7cn0KZGF0YSA8LQogIGRhdGEgJT4lIAogIGZpbHRlcihsYWJlbCE9InIsciIpCmBgYAoKIyMjIFNhdmUgZGF0YQoKYGBge3J9CndyaXRlX3RzdihkYXRhLCBmaWxlLnBhdGgodmFsaWRhdGlvbl9vdXRwdXRfZGlyLCJoaXNlcTQwMDBfaW5uZXJfam9pbmVkX3dpdGhfbGFiZWxzLnR4dCIpKQpgYGAKCgojIyBTdW1tYXJ5IHN0YXRpc3RpY3MKCm51bWJlciBvZiByZWFkcyAoaW4gbWlsbGlvbnMpCgpgYGB7cn0KZGF0YV9maiAlPiUgIAogIHVuZ3JvdXAoKSAlPiUKICBzdW1tYXJpemVfYXQodmFycyhtYXRjaGVzKCJecyIpKSwgbGlzdCh+IHN1bSguLzEwXjYsICBuYS5ybSA9IFRSVUUpKSkKYGBgCgoKYGBge3J9CmRhdGEgICU+JSAKICBzdW1tYXJpemVfYXQodmFycyhtYXRjaGVzKCJecyIpKSwgbGlzdCh+IHN1bSguLzEwXjYpKSkKYGBgCgoKbnVtYmVyIG9mIG1vbGVjdWxlcyAoaW4gbWlsbGlvbnMpCgpgYGB7cn0KZGF0YV9maiAlPiUgIAogIG11dGF0ZV9hdCh2YXJzKG1hdGNoZXMoIl5zIikpLCBsaXN0KH4gYXMuaW50ZWdlciguIT0wKSkpICAlPiUKICB1bmdyb3VwKCkgJT4lCiAgc3VtbWFyaXplX2F0KHZhcnMobWF0Y2hlcygiXnMiKSksIGxpc3QofiBzdW0oLi8xMF42LCBuYS5ybSA9IFRSVUUpKSkKYGBgCgoKYGBge3J9CiBkYXRhICAlPiUgCiAgbXV0YXRlX2F0KHZhcnMobWF0Y2hlcygiXnMiKSksIGxpc3QofiBhcy5pbnRlZ2VyKC4hPTApKSkgICU+JQogIHVuZ3JvdXAoKSAgJT4lIAogIHN1bW1hcml6ZV9hdCh2YXJzKG1hdGNoZXMoIl5zIikpLCBsaXN0KH4gc3VtKC4vMTBeNiwgbmEucm0gPSBUUlVFKSkpCmBgYAoKCiMgRXhhbWluZSBjb25jb3JkYW5jZSBiZXR3ZWVuIHNhbXBsZXMKCgojIyBSZWFkIGNvdW50cyB3aXRoIHNhbWUgQ1VHIChjZWxsLXVtaS1nZW5lKSBsYWJlbAoKYGBge3IsIGZpZy5oZWlnaHQ9MTB9CnAxIDwtIGdncGxvdChkYXRhLCBhZXMoeCA9IHMxX25vbnBsZXhlZCwKICAgICAgICAgICAgICAgICAgeSA9IHMxX3BsZXhlZCkpIApwMSArIGdlb21faGV4KGJpbnMgPSA1MDApKwogIGdlb21fYWJsaW5lKHNsb3BlPS41LGludGVyY2VwdD0wKQpgYGAKCgoKCmBgYHtyLCBmaWcuaGVpZ2h0PTEwfQpwMiA8LSBnZ3Bsb3QoZGF0YSwgCiAgICAgICAgICAgIGFlcyh4ID0gczJfbm9ucGxleGVkLAogICAgICAgICAgICAgICAgICB5ID0gczJfcGxleGVkKSkgCnAyICsgCiAgZ2VvbV9oZXgoYmlucyA9IDUwMCkrCiAgZ2VvbV9hYmxpbmUoc2xvcGU9LjUsaW50ZXJjZXB0PTApCiAjICAgZ2VvbV9wb2ludCgpIApgYGAKCgoKCiMjIEdlbmUgZXhwcmVzc2lvbiByZWFkIGNvdW50cyB3aXRoIHNhbWUgY2VsbC1nZW5lIGxhYmVsCgoKYGBge3J9CmRhdGFfcmVhZHNfY2VsbF9nZW5lIDwtIAogICAgZGF0YSAlPiUgIAogIGdyb3VwX2J5KGNlbGwsIGdlbmUpICU+JQogIHN1bW1hcml6ZV9hdCh2YXJzKG1hdGNoZXMoIl5zIikpLAogICAgICAgICAgICAgICBsaXN0KH4gc3VtKC4pKSkKYGBgCgoKCgoKYGBge3IsIGZpZy5oZWlnaHQ9MTB9CiBnZ3Bsb3QoZGF0YV9yZWFkc19jZWxsX2dlbmUsCiAgICAgICAgICAgIGFlcyh4ID0gczFfbm9ucGxleGVkLAogICAgICAgICAgICAgICAgICB5ID0gczFfcGxleGVkKSkgICsgCiAgZ2VvbV9oZXgoYmlucyA9IDQwMCkgKwogIGdlb21fYWJsaW5lKHNsb3BlPTEsaW50ZXJjZXB0PTApKwogIHNjYWxlX3hfbG9nMTAoKSArCiAgIHNjYWxlX3lfbG9nMTAoKSAKYGBgCgoKCgpgYGB7ciwgZmlnLmhlaWdodD0xMH0KZ2dwbG90KGRhdGFfcmVhZHNfY2VsbF9nZW5lLAogICAgICAgICAgICBhZXMoeCA9IHMyX25vbnBsZXhlZCwKICAgICAgICAgICAgICAgICAgeSA9IHMyX3BsZXhlZCkpICArIAogIGdlb21faGV4KGJpbnMgPSA0MDApICsKICAgIGdlb21fYWJsaW5lKHNsb3BlPTEsaW50ZXJjZXB0PTApKwogIHNjYWxlX3hfbG9nMTAoKSArCiAgIHNjYWxlX3lfbG9nMTAoKSAKIApgYGAKCmBgYHtyLCBmaWcuaGVpZ2h0PTEwfQpnZ3Bsb3QoZGF0YV9yZWFkc19jZWxsX2dlbmUsCiAgICAgICAgICAgIGFlcyh4ID0gczFfbm9ucGxleGVkLAogICAgICAgICAgICAgICAgICB5ID0gczJfbm9ucGxleGVkKSkgICsgCiAgZ2VvbV9oZXgoYmlucyA9IDQwMCkgKwogIGdlb21fYWJsaW5lKHNsb3BlPTEsaW50ZXJjZXB0PTApKwogIHNjYWxlX3hfbG9nMTAoKSArCiAgIHNjYWxlX3lfbG9nMTAoKSAKIApgYGAKCgpgYGB7ciwgZmlnLmhlaWdodD0xMH0KZ2dwbG90KGRhdGFfcmVhZHNfY2VsbF9nZW5lLAogICAgICAgICAgICBhZXMoeCA9IHMxX3BsZXhlZCwKICAgICAgICAgICAgICAgICAgeSA9IHMyX3BsZXhlZCkpICArIAogIGdlb21faGV4KGJpbnMgPSA0MDApICsKICBnZW9tX2FibGluZShzbG9wZT0xLGludGVyY2VwdD0wKSsKICBzY2FsZV94X2xvZzEwKCkgKwogICBzY2FsZV95X2xvZzEwKCkgCiAKYGBgCgpgYGB7ciwgZmlnLmhlaWdodD0xMH0KZ2dwbG90KGRhdGFfcmVhZHNfY2VsbF9nZW5lLAogICAgICAgICAgICBhZXMoeCA9IHMxX3BsZXhlZCwKICAgICAgICAgICAgICAgICAgeSA9IHMyX25vbnBsZXhlZCkpICArIAogIGdlb21faGV4KGJpbnMgPSA0MDApICsKICBnZW9tX2FibGluZShzbG9wZT0xLGludGVyY2VwdD0wKSsKICBzY2FsZV94X2xvZzEwKCkgKwogICBzY2FsZV95X2xvZzEwKCkgCiAKYGBgCgogIyMgR2VuZSBleHByZXNzaW9uIFVNSSBjb3VudHMgd2l0aCBzYW1lIGNlbGwtZ2VuZSBsYWJlbAogCgpgYGB7cn0KZGF0YV9tb2xlY19jZWxsX2dlbmUgPC0gCiAgICBkYXRhICU+JSAgCiAgbXV0YXRlX2F0KHZhcnMobWF0Y2hlcygiXnMiKSksIAogICAgICAgICAgICBsaXN0KH4gYXMuaW50ZWdlciguIT0wKSkpICAlPiUgIAogIGdyb3VwX2J5KGNlbGwsIGdlbmUpICU+JQogIHN1bW1hcml6ZV9hdCh2YXJzKG1hdGNoZXMoIl5zIikpLAogICAgICAgICAgICAgICBsaXN0KH4gc3VtKC4pKSkKYGBgCgoKCmBgYHtyLCBmaWcuaGVpZ2h0PTEwfQogZ2dwbG90KGRhdGFfbW9sZWNfY2VsbF9nZW5lLAogICAgICAgICAgICBhZXMoeCA9IHMxX25vbnBsZXhlZCwKICAgICAgICAgICAgICAgICAgeSA9IHMxX3BsZXhlZCkpICArIAogIGdlb21faGV4KGJpbnMgPSA0MDApICsKICBnZW9tX2FibGluZShzbG9wZT0xLGludGVyY2VwdD0wKSsKICBzY2FsZV94X2xvZzEwKCkgKwogICBzY2FsZV95X2xvZzEwKCkgCmBgYAoKCgoKYGBge3IsIGZpZy5oZWlnaHQ9MTB9CmdncGxvdChkYXRhX21vbGVjX2NlbGxfZ2VuZSwKICAgICAgICAgICAgYWVzKHggPSBzMl9ub25wbGV4ZWQsCiAgICAgICAgICAgICAgICAgIHkgPSBzMl9wbGV4ZWQpKSAgKyAKICBnZW9tX2hleChiaW5zID0gNDAwKSArCiAgZ2VvbV9hYmxpbmUoc2xvcGU9MSxpbnRlcmNlcHQ9MCkrCiAgc2NhbGVfeF9sb2cxMCgpICsKICAgc2NhbGVfeV9sb2cxMCgpIAogCmBgYAoKYGBge3IsIGZpZy5oZWlnaHQ9MTB9CmdncGxvdChkYXRhX21vbGVjX2NlbGxfZ2VuZSwKICAgICAgICAgICAgYWVzKHggPSBzMV9ub25wbGV4ZWQsCiAgICAgICAgICAgICAgICAgIHkgPSBzMl9ub25wbGV4ZWQpKSAgKyAKICBnZW9tX2hleChiaW5zID0gNDAwKSArCiAgZ2VvbV9hYmxpbmUoc2xvcGU9MSxpbnRlcmNlcHQ9MCkrCiAgc2NhbGVfeF9sb2cxMCgpICsKICAgc2NhbGVfeV9sb2cxMCgpIAogCmBgYAoKIyMjIEEgdmlzdWFsIGRlbW9uc3RyYXRpb24gb2YgdGhlIGVmZmVjdHMgb2YgaW5kZXggaG9wcGluZwoKYGBge3IsIGZpZy5oZWlnaHQ9MTB9CmdncGxvdChkYXRhX21vbGVjX2NlbGxfZ2VuZSwKICAgICAgICAgICAgYWVzKHggPSBzMV9wbGV4ZWQsCiAgICAgICAgICAgICAgICAgIHkgPSBzMl9wbGV4ZWQpKSAgKyAKICBnZW9tX2hleChiaW5zID0gNDAwKSArCiAgZ2VvbV9hYmxpbmUoc2xvcGU9MSxpbnRlcmNlcHQ9MCkrCiAgc2NhbGVfeF9sb2cxMCgpICsKICAgc2NhbGVfeV9sb2cxMCgpIAogCmBgYAoKCgoKYGBge3IsIGZpZy5oZWlnaHQ9MTB9CmdncGxvdChkYXRhX21vbGVjX2NlbGxfZ2VuZSwKICAgICAgICAgICAgYWVzKHggPSBzMV9wbGV4ZWQsCiAgICAgICAgICAgICAgICAgIHkgPSBzMl9ub25wbGV4ZWQpKSAgKyAKICBnZW9tX2hleChiaW5zID0gNDAwKSArCiAgZ2VvbV9hYmxpbmUoc2xvcGU9MSxpbnRlcmNlcHQ9MCkrCiAgc2NhbGVfeF9sb2cxMCgpICsKICAgc2NhbGVfeV9sb2cxMCgpIAogCmBgYAoKCgoKCgpgYGB7ciwgZmlnLmhlaWdodD0xMH0KZ2dwbG90KGRhdGFfbW9sZWNfY2VsbF9nZW5lLAogICAgICAgICAgICBhZXMoeCA9IHMyX3BsZXhlZCwKICAgICAgICAgICAgICAgICAgeSA9IHMxX25vbnBsZXhlZCkpICArIAogIGdlb21faGV4KGJpbnMgPSA0MDApICsKICBnZW9tX2FibGluZShzbG9wZT0xLGludGVyY2VwdD0wKSsKICBzY2FsZV94X2xvZzEwKCkgKwogICBzY2FsZV95X2xvZzEwKCkgCiAKYGBgCiAKICMjIEdlbmUgZXhwcmVzc2lvbiBVTUkgY291bnRzIChub25ob3BwaW5nKSB3aXRoIHNhbWUgY2VsbC1nZW5lIGxhYmVsCgoKSGVyZSB3ZSByZXRhaW4gbW9sZWN1bGVzIGZvciB3aGljaCB0aGUgcmVhZHMgZGlkbid0IGhvcAoKYGBge3J9CmRhdGFfbW9sZWNfY2VsbF9nZW5lX25vbmhvcHBpbmcgPC0gCiAgICBkYXRhICU+JSAgIAogIGZpbHRlcihsYWJlbCVpbiUgYygiciwwIiwiMCxyIikpJT4lCiAgbXV0YXRlX2F0KHZhcnMobWF0Y2hlcygiXnMiKSksIAogICAgICAgICAgICBsaXN0KH4gYXMuaW50ZWdlciguIT0wKSkpICAlPiUgIAogIGdyb3VwX2J5KGNlbGwsIGdlbmUpICU+JQogIHN1bW1hcml6ZV9hdCh2YXJzKG1hdGNoZXMoIl5zIikpLAogICAgICAgICAgICAgICBsaXN0KH4gc3VtKC4pKSkKYGBgCgoKCmBgYHtyLCBmaWcuaGVpZ2h0PTEwfQogZ2dwbG90KGRhdGFfbW9sZWNfY2VsbF9nZW5lX25vbmhvcHBpbmcsCiAgICAgICAgICAgIGFlcyh4ID0gczFfbm9ucGxleGVkLAogICAgICAgICAgICAgICAgICB5ID0gczFfcGxleGVkKSkgICsgCiAgZ2VvbV9oZXgoYmlucyA9IDQwMCkgKwogIGdlb21fYWJsaW5lKHNsb3BlPTEsaW50ZXJjZXB0PTApKwogIHNjYWxlX3hfbG9nMTAoKSArCiAgIHNjYWxlX3lfbG9nMTAoKSAKYGBgCgoKCgpgYGB7ciwgZmlnLmhlaWdodD0xMH0KZ2dwbG90KGRhdGFfbW9sZWNfY2VsbF9nZW5lX25vbmhvcHBpbmcsCiAgICAgICAgICAgIGFlcyh4ID0gczJfbm9ucGxleGVkLAogICAgICAgICAgICAgICAgICB5ID0gczJfcGxleGVkKSkgICsgCiAgZ2VvbV9oZXgoYmlucyA9IDQwMCkgKwogIGdlb21fYWJsaW5lKHNsb3BlPTEsaW50ZXJjZXB0PTApKwogIHNjYWxlX3hfbG9nMTAoKSArCiAgIHNjYWxlX3lfbG9nMTAoKSAKIApgYGAKCmBgYHtyLCBmaWcuaGVpZ2h0PTEwfQpnZ3Bsb3QoZGF0YV9tb2xlY19jZWxsX2dlbmVfbm9uaG9wcGluZywKICAgICAgICAgICAgYWVzKHggPSBzMV9ub25wbGV4ZWQsCiAgICAgICAgICAgICAgICAgIHkgPSBzMl9ub25wbGV4ZWQpKSAgKyAKICBnZW9tX2hleChiaW5zID0gNDAwKSArCiAgZ2VvbV9hYmxpbmUoc2xvcGU9MSxpbnRlcmNlcHQ9MCkrCiAgc2NhbGVfeF9sb2cxMCgpICsKICAgc2NhbGVfeV9sb2cxMCgpIAogCmBgYAoKCmBgYHtyLCBmaWcuaGVpZ2h0PTEwfQpnZ3Bsb3QoZGF0YV9tb2xlY19jZWxsX2dlbmVfbm9uaG9wcGluZywKICAgICAgICAgICAgYWVzKHggPSBzMV9wbGV4ZWQsCiAgICAgICAgICAgICAgICAgIHkgPSBzMl9wbGV4ZWQpKSAgKyAKICBnZW9tX2hleChiaW5zID0gNDAwKSArCiAgZ2VvbV9hYmxpbmUoc2xvcGU9MSxpbnRlcmNlcHQ9MCkrCiAgc2NhbGVfeF9sb2cxMCgpICsKICAgc2NhbGVfeV9sb2cxMCgpIAogCmBgYAoKCgoKYGBge3IsIGZpZy5oZWlnaHQ9MTB9CmdncGxvdChkYXRhX21vbGVjX2NlbGxfZ2VuZV9ub25ob3BwaW5nLAogICAgICAgICAgICBhZXMoeCA9IHMxX3BsZXhlZCwKICAgICAgICAgICAgICAgICAgeSA9IHMyX25vbnBsZXhlZCkpICArIAogIGdlb21faGV4KGJpbnMgPSA0MDApICsKICBnZW9tX2FibGluZShzbG9wZT0xLGludGVyY2VwdD0wKSsKICBzY2FsZV94X2xvZzEwKCkgKwogICBzY2FsZV95X2xvZzEwKCkgCiAKYGBgCgoKCmBgYHtyfQpzZXNzaW9uSW5mbygpCmBgYAoKCg==