saveRDS(snakemake, snakemake@log$snakemake)

suppressPackageStartupMessages({
  library(OUTRIDER)
  library(SummarizedExperiment)
  library(GenomicAlignments)
  library(ggplot2)
  library(ggthemes)
  library(cowplot)
  library(data.table)
  library(tidyr)
})

ods <- readRDS(snakemake@input$ods)
ods@colData$isExternal <- F

has_external <- any(as.logical(colData(ods)$isExternal))
cnts_mtx_local <- counts(ods, normalized = F)[,!as.logical(ods@colData$isExternal),drop=FALSE]
cnts_mtx <- counts(ods, normalized = F)

Number of samples:

Local: 541
External: 0

Count Quality Control

Compares total mapped vs counted reads.
The Mapped vs Counted Reads plot does not include external counts.
Consider removing samples with too low or too high size factors.

bam_coverage <- fread(snakemake@input$bam_cov)
bam_coverage[, sampleID := as.character(sampleID)]
setnames(bam_coverage, 'record_count', 'total_count')
coverage_dt <- merge(bam_coverage,
                     data.table(sampleID = colnames(ods),
                                read_count = colSums(cnts_mtx),
                                isExternal = ods@colData$isExternal),
                     by = "sampleID", sort = FALSE)
# read counts
coverage_dt[, count_rank := rank(read_count)]

# size factors 
ods <- estimateSizeFactors(ods)
coverage_dt[, size_factors := round(sizeFactors(ods), 3)]
coverage_dt[, sf_rank := rank(size_factors)]

# Show this plot only if there are external samples, as the next plot contains the same info
if(has_external == T){
  p_depth <- ggplot(coverage_dt, aes(x = count_rank, y = read_count, col = isExternal)) +
    geom_point(size = 3,show.legend = has_external) +
    theme_cowplot() + background_grid() +
    labs(title = "Obtained Read Counts", x="Sample Rank", y = "Counted Reads") +
    ylim(c(0,NA)) +
    scale_color_brewer(palette="Dark2")
  p_depth
}


p_comp <- ggplot(coverage_dt[isExternal == F], aes(x = total_count, y = read_count)) +
  geom_point(size = 3, show.legend = has_external, color = "#1B9E77") +
  theme_cowplot() + background_grid() +
  labs(title = "Total mapped vs. Counted Reads", x = "Mapped Reads", y = "Counted Reads") +
  xlim(c(0,NA)) + ylim(c(0,NA))
p_comp

# ggExtra::ggMarginal(p_comp, type = "histogram") # could be a nice add-on

p_sf <- ggplot(coverage_dt, aes(sf_rank, size_factors, col = isExternal)) +
  geom_point(size = 3,show.legend = has_external) +
  ylim(c(0,NA)) +
  theme_cowplot() + background_grid() +
  labs(title = 'Size Factors', x = 'Sample Rank', y = 'Size Factors') +
  scale_color_brewer(palette="Dark2")

p_sf

setnames(coverage_dt, old = c('total_count', 'read_count', 'size_factors'),
         new = c('Reads Mapped', 'Reads Counted', 'Size Factors'))
DT::datatable(coverage_dt[, .(sampleID, `Reads Mapped`, `Reads Counted`, `Size Factors`)][order(`Reads Mapped`)],
              caption = 'Reads summary statistics')

Filtering

local: A pre-filtered summary of counts using only the local (from BAM) counts. Omitted if no external counts
all: A pre-filtered summary of counts using only the merged local (from BAM) and external counts
passed FPKM: Passes the user defined FPKM cutoff in at least 5% of genes
min 1 read: minimum of 1 read expressed in 5% of genes
min 10 reads: minimum of 10 reads expressed in 5% of genes

quant <- .95

if(has_external){
    filter_mtx <- list(
      local = cnts_mtx_local,
      all = cnts_mtx,
      `passed FPKM` = cnts_mtx[rowData(ods)$passedFilter,],
      `min 1 read` = cnts_mtx[rowQuantiles(cnts_mtx, probs = quant) > 1, ],
      `min 10 reads` = cnts_mtx[rowQuantiles(cnts_mtx, probs = quant) > 10, ]
    )
    filter_dt <- lapply(names(filter_mtx), function(filter_name) {
      mtx <- filter_mtx[[filter_name]]
      data.table(gene_ID = rownames(mtx), median_counts = rowMeans(mtx), filter = filter_name)
    }) %>% rbindlist
    filter_dt[, filter := factor(filter, levels = c('local', 'all', 'passed FPKM', 'min 1 read', 'min 10 reads'))]
} else{
    filter_mtx <- list(
      all = cnts_mtx,
      `passed FPKM` = cnts_mtx[rowData(ods)$passedFilter,],
      `min 1 read` = cnts_mtx[rowQuantiles(cnts_mtx, probs = quant) > 1, ],
      `min 10 reads` = cnts_mtx[rowQuantiles(cnts_mtx, probs = quant) > 10, ]
    )
    filter_dt <- lapply(names(filter_mtx), function(filter_name) {
      mtx <- filter_mtx[[filter_name]]
      data.table(gene_ID = rownames(mtx), median_counts = rowMeans(mtx), filter = filter_name)
    }) %>% rbindlist
    filter_dt[, filter := factor(filter, levels = c('all', 'passed FPKM', 'min 1 read', 'min 10 reads'))]
}

binwidth <- .2
p_hist <- ggplot(filter_dt, aes(x = median_counts, fill = filter)) +
  geom_histogram(binwidth = binwidth) +
  scale_x_log10() +
  facet_wrap(.~filter) +
  labs(x = "Mean counts per gene", y = "Frequency", title = 'Mean Count Distribution') +
  guides(col = guide_legend(title = NULL)) +
  scale_fill_brewer(palette = "Paired") +
  theme_cowplot() +
  theme(legend.position = "none")

p_dens <- ggplot(filter_dt, aes(x = median_counts, col = filter)) +
  geom_density(aes(y=binwidth * ..count..), size = 1.2) +
  scale_x_log10() +
  labs(x = "Mean counts per gene", y = "Frequency") +
  guides(col = guide_legend(title = NULL)) +
  scale_color_brewer(palette = "Paired") +
  theme_cowplot() +
  theme(legend.position = "top",
        legend.justification="center",
        legend.background = element_rect(color = NA))
plot_grid(p_hist, p_dens)

Expressed Genes

exp_genes_cols <- c(Rank = "expressedGenesRank",`Expressed\ngenes` = "expressedGenes", 
                    `Union of\nexpressed genes` = "unionExpressedGenes", 
                    `Intersection of\nexpressed genes` = "intersectionExpressedGenes", 
                    `Genes passed\nfiltering` = "passedFilterGenes", `Is External` = "isExternal")

expressed_genes <- as.data.table(colData(ods)[,exp_genes_cols], keep.rownames = TRUE)
colnames(expressed_genes) <- c('RNA_ID', names(exp_genes_cols))
plotExpressedGenes(ods) + 
  theme_cowplot() +
  background_grid(major = "y") +
  geom_point(data = melt(expressed_genes, id.vars = c("RNA_ID", "Rank", "Is External")),
             aes(Rank, value, col = variable, shape = `Is External`), 
             show.legend = has_external)

if(has_external){
    DT::datatable(expressed_genes[order(Rank)],rownames = F)
} else{
    DT::datatable(expressed_genes[order(Rank),-"Is External"],rownames = F)
}

Considerations: The verifying of the samples sex is performed by comparing the expression levels of the genes XIST and UTY. In general, females should have high XIST and low UTY expression, and viceversa for males. For it to work, the sample annotation must have a column called ‘SEX’, with values male/female. If some other values exist, e.g., unknown, they will be matched. The prediction is performed via a linear discriminant analysis on the log2 counts.

# Get sex column and proceed if it exists
sex_idx <- which('SEX' == toupper(colnames(colData(ods))))
if(isEmpty(sex_idx)){
  print('Sex column not found in sample annotation')
} else{
  
  # Verify if both XIST and UTY were counted
  xist_id <- 'XIST'
  uty_id <- 'UTY'
  
  if(grepl('ENSG0', rownames(ods)[1])){
    xist_id <- 'ENSG00000229807'
    uty_id <- 'ENSG00000183878'
  }
  xist_idx <- grep(xist_id, rownames(ods))
  uty_idx <- grep(uty_id, rownames(ods))
  
  if(isEmpty(xist_idx) | isEmpty(uty_idx)){
    print('Either XIST or UTY is not expressed')
  }else{
    
    sex_counts <- counts(ods)[c(xist_idx, uty_idx), ]
    sex_dt <- data.table(sampleID = colnames(ods), 
                         XIST = counts(ods)[xist_idx,], 
                         UTY = counts(ods)[uty_idx,])
    sex_dt <- merge(sex_dt, as.data.table(colData(ods))[,c(1, sex_idx), with = F], sort = F)
    colnames(sex_dt) <- toupper(colnames(sex_dt))
    sex_dt[, SEX := tolower(SEX)]
    sex_dt[is.na(SEX), SEX := '']
    
    # Train only in male/female in case there are other values
    train_dt <- sex_dt[SEX %in% c('f', 'm', 'female', 'male')]
    
    library("MASS")
    lda <- lda(SEX ~ log2(XIST+1) + log2(UTY+1), data = train_dt)
    
    sex_dt[, predicted_sex := predict(lda, sex_dt)$class]
    sex_dt[, match_sex := SEX == predicted_sex]
    table(sex_dt[, .(SEX, predicted_sex)])
    
    g <- ggplot(sex_dt, aes(XIST+1, UTY+1)) + 
      geom_point(aes(col = SEX, shape = predicted_sex, alpha = match_sex)) + 
      scale_x_log10(limits = c(1,NA)) + scale_y_log10(limits = c(1,NA)) +
      scale_alpha_manual(values = c(1, .1)) + 
      theme_cowplot() + background_grid(major = 'xy', minor = 'xy') + 
      annotation_logticks(sides = 'bl') + 
      labs(color = 'Sex', shape = 'Predicted sex', alpha = 'Matches sex')
    plot(g)
    
    DT::datatable(sex_dt[match_sex == F], caption = 'Sex mismatches')
    
    # Write table
    fwrite(sex_dt, gsub('ods_unfitted.Rds', 'xist_uty.tsv', snakemake@input$ods), 
           sep = '\t', quote = F)
  }
}
## [1] "Sex column not found in sample annotation"
IyctLS0KIycgdGl0bGU6ICJDb3VudHMgU3VtbWFyeTogYHIgcGFzdGUoc25ha2VtYWtlQHdpbGRjYXJkcyRkYXRhc2V0LCBzbmFrZW1ha2VAd2lsZGNhcmRzJGFubm90YXRpb24sIHNlcCA9ICctLScpYCIKIycgYXV0aG9yOiAKIycgd2I6CiMnICBsb2c6CiMnICAgLSBzbmFrZW1ha2U6ICdgc20gc3RyKHRtcF9kaXIgLyAiQUUiIC8gInthbm5vdGF0aW9ufSIgLyAie2RhdGFzZXR9IiAvICJjb3VudF9zdW1tYXJ5LlJkcyIpYCcKIycgIGlucHV0OiAKIycgICAgLSBvZHM6ICdgc20gY2ZnLmdldFByb2Nlc3NlZFJlc3VsdHNEaXIoKSArCiMnICAgICAgICAgICAgIi9hYmVycmFudF9leHByZXNzaW9uL3thbm5vdGF0aW9ufS9vdXRyaWRlci97ZGF0YXNldH0vb2RzX3VuZml0dGVkLlJkcyJgJwojJyAgICAtIGJhbV9jb3Y6ICdgc20gcnVsZXMuYWJlcnJhbnRFeHByZXNzaW9uX21lcmdlQmFtU3RhdHMub3V0cHV0YCcKIycgIG91dHB1dDoKIycgICAtIHdCaHRtbDogJ2BzbSBjb25maWdbImh0bWxPdXRwdXRQYXRoIl0gKwojJyAgICAgICAgICAgICAgIi9BYmVycmFudEV4cHJlc3Npb24vQ291bnRpbmcve2Fubm90YXRpb259L1N1bW1hcnlfe2RhdGFzZXR9Lmh0bWwiYCcKIycgIHR5cGU6IG5vaW5kZXgKIycgb3V0cHV0OgojJyAgaHRtbF9kb2N1bWVudDoKIycgICBjb2RlX2ZvbGRpbmc6IGhpZGUKIycgICBjb2RlX2Rvd25sb2FkOiBUUlVFCiMnLS0tCgpzYXZlUkRTKHNuYWtlbWFrZSwgc25ha2VtYWtlQGxvZyRzbmFrZW1ha2UpCgpzdXBwcmVzc1BhY2thZ2VTdGFydHVwTWVzc2FnZXMoewogIGxpYnJhcnkoT1VUUklERVIpCiAgbGlicmFyeShTdW1tYXJpemVkRXhwZXJpbWVudCkKICBsaWJyYXJ5KEdlbm9taWNBbGlnbm1lbnRzKQogIGxpYnJhcnkoZ2dwbG90MikKICBsaWJyYXJ5KGdndGhlbWVzKQogIGxpYnJhcnkoY293cGxvdCkKICBsaWJyYXJ5KGRhdGEudGFibGUpCiAgbGlicmFyeSh0aWR5cikKfSkKCm9kcyA8LSByZWFkUkRTKHNuYWtlbWFrZUBpbnB1dCRvZHMpCm9kc0Bjb2xEYXRhJGlzRXh0ZXJuYWwgPC0gRgoKaGFzX2V4dGVybmFsIDwtIGFueShhcy5sb2dpY2FsKGNvbERhdGEob2RzKSRpc0V4dGVybmFsKSkKY250c19tdHhfbG9jYWwgPC0gY291bnRzKG9kcywgbm9ybWFsaXplZCA9IEYpWywhYXMubG9naWNhbChvZHNAY29sRGF0YSRpc0V4dGVybmFsKSxkcm9wPUZBTFNFXQpjbnRzX210eCA8LSBjb3VudHMob2RzLCBub3JtYWxpemVkID0gRikKCiMnICMjIE51bWJlciBvZiBzYW1wbGVzOiAgCiMnIExvY2FsOiBgciBzdW0oIWFzLmxvZ2ljYWwob2RzQGNvbERhdGEkaXNFeHRlcm5hbCkpYCAgCiMnIEV4dGVybmFsOiBgciBzdW0oYXMubG9naWNhbChvZHNAY29sRGF0YSRpc0V4dGVybmFsKSlgICAKIycgCiMnICMgQ291bnQgUXVhbGl0eSBDb250cm9sCiMnIAojJyBDb21wYXJlcyB0b3RhbCBtYXBwZWQgdnMgY291bnRlZCByZWFkcy4gIAojJyBUaGUgYE1hcHBlZCB2cyBDb3VudGVkIFJlYWRzYCBwbG90IGRvZXMgbm90IGluY2x1ZGUgZXh0ZXJuYWwgY291bnRzLiAgCiMnIENvbnNpZGVyIHJlbW92aW5nIHNhbXBsZXMgd2l0aCB0b28gbG93IG9yIHRvbyBoaWdoIHNpemUgZmFjdG9ycy4KIycgIApiYW1fY292ZXJhZ2UgPC0gZnJlYWQoc25ha2VtYWtlQGlucHV0JGJhbV9jb3YpCmJhbV9jb3ZlcmFnZVssIHNhbXBsZUlEIDo9IGFzLmNoYXJhY3RlcihzYW1wbGVJRCldCnNldG5hbWVzKGJhbV9jb3ZlcmFnZSwgJ3JlY29yZF9jb3VudCcsICd0b3RhbF9jb3VudCcpCmNvdmVyYWdlX2R0IDwtIG1lcmdlKGJhbV9jb3ZlcmFnZSwKICAgICAgICAgICAgICAgICAgICAgZGF0YS50YWJsZShzYW1wbGVJRCA9IGNvbG5hbWVzKG9kcyksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgcmVhZF9jb3VudCA9IGNvbFN1bXMoY250c19tdHgpLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGlzRXh0ZXJuYWwgPSBvZHNAY29sRGF0YSRpc0V4dGVybmFsKSwKICAgICAgICAgICAgICAgICAgICAgYnkgPSAic2FtcGxlSUQiLCBzb3J0ID0gRkFMU0UpCiMgcmVhZCBjb3VudHMKY292ZXJhZ2VfZHRbLCBjb3VudF9yYW5rIDo9IHJhbmsocmVhZF9jb3VudCldCgojIHNpemUgZmFjdG9ycyAKb2RzIDwtIGVzdGltYXRlU2l6ZUZhY3RvcnMob2RzKQpjb3ZlcmFnZV9kdFssIHNpemVfZmFjdG9ycyA6PSByb3VuZChzaXplRmFjdG9ycyhvZHMpLCAzKV0KY292ZXJhZ2VfZHRbLCBzZl9yYW5rIDo9IHJhbmsoc2l6ZV9mYWN0b3JzKV0KCiMgU2hvdyB0aGlzIHBsb3Qgb25seSBpZiB0aGVyZSBhcmUgZXh0ZXJuYWwgc2FtcGxlcywgYXMgdGhlIG5leHQgcGxvdCBjb250YWlucyB0aGUgc2FtZSBpbmZvCmlmKGhhc19leHRlcm5hbCA9PSBUKXsKICBwX2RlcHRoIDwtIGdncGxvdChjb3ZlcmFnZV9kdCwgYWVzKHggPSBjb3VudF9yYW5rLCB5ID0gcmVhZF9jb3VudCwgY29sID0gaXNFeHRlcm5hbCkpICsKICAgIGdlb21fcG9pbnQoc2l6ZSA9IDMsc2hvdy5sZWdlbmQgPSBoYXNfZXh0ZXJuYWwpICsKICAgIHRoZW1lX2Nvd3Bsb3QoKSArIGJhY2tncm91bmRfZ3JpZCgpICsKICAgIGxhYnModGl0bGUgPSAiT2J0YWluZWQgUmVhZCBDb3VudHMiLCB4PSJTYW1wbGUgUmFuayIsIHkgPSAiQ291bnRlZCBSZWFkcyIpICsKICAgIHlsaW0oYygwLE5BKSkgKwogICAgc2NhbGVfY29sb3JfYnJld2VyKHBhbGV0dGU9IkRhcmsyIikKICBwX2RlcHRoCn0KCgpwX2NvbXAgPC0gZ2dwbG90KGNvdmVyYWdlX2R0W2lzRXh0ZXJuYWwgPT0gRl0sIGFlcyh4ID0gdG90YWxfY291bnQsIHkgPSByZWFkX2NvdW50KSkgKwogIGdlb21fcG9pbnQoc2l6ZSA9IDMsIHNob3cubGVnZW5kID0gaGFzX2V4dGVybmFsLCBjb2xvciA9ICIjMUI5RTc3IikgKwogIHRoZW1lX2Nvd3Bsb3QoKSArIGJhY2tncm91bmRfZ3JpZCgpICsKICBsYWJzKHRpdGxlID0gIlRvdGFsIG1hcHBlZCB2cy4gQ291bnRlZCBSZWFkcyIsIHggPSAiTWFwcGVkIFJlYWRzIiwgeSA9ICJDb3VudGVkIFJlYWRzIikgKwogIHhsaW0oYygwLE5BKSkgKyB5bGltKGMoMCxOQSkpCnBfY29tcAojIGdnRXh0cmE6OmdnTWFyZ2luYWwocF9jb21wLCB0eXBlID0gImhpc3RvZ3JhbSIpICMgY291bGQgYmUgYSBuaWNlIGFkZC1vbgoKcF9zZiA8LSBnZ3Bsb3QoY292ZXJhZ2VfZHQsIGFlcyhzZl9yYW5rLCBzaXplX2ZhY3RvcnMsIGNvbCA9IGlzRXh0ZXJuYWwpKSArCiAgZ2VvbV9wb2ludChzaXplID0gMyxzaG93LmxlZ2VuZCA9IGhhc19leHRlcm5hbCkgKwogIHlsaW0oYygwLE5BKSkgKwogIHRoZW1lX2Nvd3Bsb3QoKSArIGJhY2tncm91bmRfZ3JpZCgpICsKICBsYWJzKHRpdGxlID0gJ1NpemUgRmFjdG9ycycsIHggPSAnU2FtcGxlIFJhbmsnLCB5ID0gJ1NpemUgRmFjdG9ycycpICsKICBzY2FsZV9jb2xvcl9icmV3ZXIocGFsZXR0ZT0iRGFyazIiKQoKcF9zZgoKc2V0bmFtZXMoY292ZXJhZ2VfZHQsIG9sZCA9IGMoJ3RvdGFsX2NvdW50JywgJ3JlYWRfY291bnQnLCAnc2l6ZV9mYWN0b3JzJyksCiAgICAgICAgIG5ldyA9IGMoJ1JlYWRzIE1hcHBlZCcsICdSZWFkcyBDb3VudGVkJywgJ1NpemUgRmFjdG9ycycpKQpEVDo6ZGF0YXRhYmxlKGNvdmVyYWdlX2R0WywgLihzYW1wbGVJRCwgYFJlYWRzIE1hcHBlZGAsIGBSZWFkcyBDb3VudGVkYCwgYFNpemUgRmFjdG9yc2ApXVtvcmRlcihgUmVhZHMgTWFwcGVkYCldLAogICAgICAgICAgICAgIGNhcHRpb24gPSAnUmVhZHMgc3VtbWFyeSBzdGF0aXN0aWNzJykKCiMnICMgRmlsdGVyaW5nCiMnICoqbG9jYWwqKjogQSBwcmUtZmlsdGVyZWQgc3VtbWFyeSBvZiBjb3VudHMgdXNpbmcgb25seSB0aGUgbG9jYWwgKGZyb20gQkFNKSBjb3VudHMuIE9taXR0ZWQgaWYgbm8gZXh0ZXJuYWwgY291bnRzICAKIycgKiphbGwqKjogQSBwcmUtZmlsdGVyZWQgc3VtbWFyeSBvZiBjb3VudHMgdXNpbmcgb25seSB0aGUgbWVyZ2VkIGxvY2FsIChmcm9tIEJBTSkgYW5kIGV4dGVybmFsIGNvdW50cyAgCiMnICoqcGFzc2VkIEZQS00qKjogUGFzc2VzIHRoZSB1c2VyIGRlZmluZWQgRlBLTSBjdXRvZmYgaW4gYXQgbGVhc3QgNSUgb2YgZ2VuZXMgIAojJyAqKm1pbiAxIHJlYWQqKjogbWluaW11bSBvZiAxIHJlYWQgZXhwcmVzc2VkIGluIDUlIG9mIGdlbmVzICAKIycgKiptaW4gMTAgcmVhZHMqKjogbWluaW11bSBvZiAxMCByZWFkcyBleHByZXNzZWQgaW4gNSUgb2YgZ2VuZXMgIAoKcXVhbnQgPC0gLjk1CgppZihoYXNfZXh0ZXJuYWwpewogICAgZmlsdGVyX210eCA8LSBsaXN0KAogICAgICBsb2NhbCA9IGNudHNfbXR4X2xvY2FsLAogICAgICBhbGwgPSBjbnRzX210eCwKICAgICAgYHBhc3NlZCBGUEtNYCA9IGNudHNfbXR4W3Jvd0RhdGEob2RzKSRwYXNzZWRGaWx0ZXIsXSwKICAgICAgYG1pbiAxIHJlYWRgID0gY250c19tdHhbcm93UXVhbnRpbGVzKGNudHNfbXR4LCBwcm9icyA9IHF1YW50KSA+IDEsIF0sCiAgICAgIGBtaW4gMTAgcmVhZHNgID0gY250c19tdHhbcm93UXVhbnRpbGVzKGNudHNfbXR4LCBwcm9icyA9IHF1YW50KSA+IDEwLCBdCiAgICApCiAgICBmaWx0ZXJfZHQgPC0gbGFwcGx5KG5hbWVzKGZpbHRlcl9tdHgpLCBmdW5jdGlvbihmaWx0ZXJfbmFtZSkgewogICAgICBtdHggPC0gZmlsdGVyX210eFtbZmlsdGVyX25hbWVdXQogICAgICBkYXRhLnRhYmxlKGdlbmVfSUQgPSByb3duYW1lcyhtdHgpLCBtZWRpYW5fY291bnRzID0gcm93TWVhbnMobXR4KSwgZmlsdGVyID0gZmlsdGVyX25hbWUpCiAgICB9KSAlPiUgcmJpbmRsaXN0CiAgICBmaWx0ZXJfZHRbLCBmaWx0ZXIgOj0gZmFjdG9yKGZpbHRlciwgbGV2ZWxzID0gYygnbG9jYWwnLCAnYWxsJywgJ3Bhc3NlZCBGUEtNJywgJ21pbiAxIHJlYWQnLCAnbWluIDEwIHJlYWRzJykpXQp9IGVsc2V7CiAgICBmaWx0ZXJfbXR4IDwtIGxpc3QoCiAgICAgIGFsbCA9IGNudHNfbXR4LAogICAgICBgcGFzc2VkIEZQS01gID0gY250c19tdHhbcm93RGF0YShvZHMpJHBhc3NlZEZpbHRlcixdLAogICAgICBgbWluIDEgcmVhZGAgPSBjbnRzX210eFtyb3dRdWFudGlsZXMoY250c19tdHgsIHByb2JzID0gcXVhbnQpID4gMSwgXSwKICAgICAgYG1pbiAxMCByZWFkc2AgPSBjbnRzX210eFtyb3dRdWFudGlsZXMoY250c19tdHgsIHByb2JzID0gcXVhbnQpID4gMTAsIF0KICAgICkKICAgIGZpbHRlcl9kdCA8LSBsYXBwbHkobmFtZXMoZmlsdGVyX210eCksIGZ1bmN0aW9uKGZpbHRlcl9uYW1lKSB7CiAgICAgIG10eCA8LSBmaWx0ZXJfbXR4W1tmaWx0ZXJfbmFtZV1dCiAgICAgIGRhdGEudGFibGUoZ2VuZV9JRCA9IHJvd25hbWVzKG10eCksIG1lZGlhbl9jb3VudHMgPSByb3dNZWFucyhtdHgpLCBmaWx0ZXIgPSBmaWx0ZXJfbmFtZSkKICAgIH0pICU+JSByYmluZGxpc3QKICAgIGZpbHRlcl9kdFssIGZpbHRlciA6PSBmYWN0b3IoZmlsdGVyLCBsZXZlbHMgPSBjKCdhbGwnLCAncGFzc2VkIEZQS00nLCAnbWluIDEgcmVhZCcsICdtaW4gMTAgcmVhZHMnKSldCn0KCmJpbndpZHRoIDwtIC4yCnBfaGlzdCA8LSBnZ3Bsb3QoZmlsdGVyX2R0LCBhZXMoeCA9IG1lZGlhbl9jb3VudHMsIGZpbGwgPSBmaWx0ZXIpKSArCiAgZ2VvbV9oaXN0b2dyYW0oYmlud2lkdGggPSBiaW53aWR0aCkgKwogIHNjYWxlX3hfbG9nMTAoKSArCiAgZmFjZXRfd3JhcCgufmZpbHRlcikgKwogIGxhYnMoeCA9ICJNZWFuIGNvdW50cyBwZXIgZ2VuZSIsIHkgPSAiRnJlcXVlbmN5IiwgdGl0bGUgPSAnTWVhbiBDb3VudCBEaXN0cmlidXRpb24nKSArCiAgZ3VpZGVzKGNvbCA9IGd1aWRlX2xlZ2VuZCh0aXRsZSA9IE5VTEwpKSArCiAgc2NhbGVfZmlsbF9icmV3ZXIocGFsZXR0ZSA9ICJQYWlyZWQiKSArCiAgdGhlbWVfY293cGxvdCgpICsKICB0aGVtZShsZWdlbmQucG9zaXRpb24gPSAibm9uZSIpCgpwX2RlbnMgPC0gZ2dwbG90KGZpbHRlcl9kdCwgYWVzKHggPSBtZWRpYW5fY291bnRzLCBjb2wgPSBmaWx0ZXIpKSArCiAgZ2VvbV9kZW5zaXR5KGFlcyh5PWJpbndpZHRoICogLi5jb3VudC4uKSwgc2l6ZSA9IDEuMikgKwogIHNjYWxlX3hfbG9nMTAoKSArCiAgbGFicyh4ID0gIk1lYW4gY291bnRzIHBlciBnZW5lIiwgeSA9ICJGcmVxdWVuY3kiKSArCiAgZ3VpZGVzKGNvbCA9IGd1aWRlX2xlZ2VuZCh0aXRsZSA9IE5VTEwpKSArCiAgc2NhbGVfY29sb3JfYnJld2VyKHBhbGV0dGUgPSAiUGFpcmVkIikgKwogIHRoZW1lX2Nvd3Bsb3QoKSArCiAgdGhlbWUobGVnZW5kLnBvc2l0aW9uID0gInRvcCIsCiAgICAgICAgbGVnZW5kLmp1c3RpZmljYXRpb249ImNlbnRlciIsCiAgICAgICAgbGVnZW5kLmJhY2tncm91bmQgPSBlbGVtZW50X3JlY3QoY29sb3IgPSBOQSkpCgojKyBtZWFuQ291bnRzLCBmaWcuaGVpZ2h0PTYsIGZpZy53aWR0aD0xMgpwbG90X2dyaWQocF9oaXN0LCBwX2RlbnMpCgojJyAjIyMgRXhwcmVzc2VkIEdlbmVzCmV4cF9nZW5lc19jb2xzIDwtIGMoUmFuayA9ICJleHByZXNzZWRHZW5lc1JhbmsiLGBFeHByZXNzZWRcbmdlbmVzYCA9ICJleHByZXNzZWRHZW5lcyIsIAogICAgICAgICAgICAgICAgICAgIGBVbmlvbiBvZlxuZXhwcmVzc2VkIGdlbmVzYCA9ICJ1bmlvbkV4cHJlc3NlZEdlbmVzIiwgCiAgICAgICAgICAgICAgICAgICAgYEludGVyc2VjdGlvbiBvZlxuZXhwcmVzc2VkIGdlbmVzYCA9ICJpbnRlcnNlY3Rpb25FeHByZXNzZWRHZW5lcyIsIAogICAgICAgICAgICAgICAgICAgIGBHZW5lcyBwYXNzZWRcbmZpbHRlcmluZ2AgPSAicGFzc2VkRmlsdGVyR2VuZXMiLCBgSXMgRXh0ZXJuYWxgID0gImlzRXh0ZXJuYWwiKQoKZXhwcmVzc2VkX2dlbmVzIDwtIGFzLmRhdGEudGFibGUoY29sRGF0YShvZHMpWyxleHBfZ2VuZXNfY29sc10sIGtlZXAucm93bmFtZXMgPSBUUlVFKQpjb2xuYW1lcyhleHByZXNzZWRfZ2VuZXMpIDwtIGMoJ1JOQV9JRCcsIG5hbWVzKGV4cF9nZW5lc19jb2xzKSkKCiMrIGV4cHJlc3NlZEdlbmVzLCBmaWcuaGVpZ2h0PTYsIGZpZy53aWR0aD04CnBsb3RFeHByZXNzZWRHZW5lcyhvZHMpICsgCiAgdGhlbWVfY293cGxvdCgpICsKICBiYWNrZ3JvdW5kX2dyaWQobWFqb3IgPSAieSIpICsKICBnZW9tX3BvaW50KGRhdGEgPSBtZWx0KGV4cHJlc3NlZF9nZW5lcywgaWQudmFycyA9IGMoIlJOQV9JRCIsICJSYW5rIiwgIklzIEV4dGVybmFsIikpLAogICAgICAgICAgICAgYWVzKFJhbmssIHZhbHVlLCBjb2wgPSB2YXJpYWJsZSwgc2hhcGUgPSBgSXMgRXh0ZXJuYWxgKSwgCiAgICAgICAgICAgICBzaG93LmxlZ2VuZCA9IGhhc19leHRlcm5hbCkKCmlmKGhhc19leHRlcm5hbCl7CiAgICBEVDo6ZGF0YXRhYmxlKGV4cHJlc3NlZF9nZW5lc1tvcmRlcihSYW5rKV0scm93bmFtZXMgPSBGKQp9IGVsc2V7CiAgICBEVDo6ZGF0YXRhYmxlKGV4cHJlc3NlZF9nZW5lc1tvcmRlcihSYW5rKSwtIklzIEV4dGVybmFsIl0scm93bmFtZXMgPSBGKQp9CgojJyAqKkNvbnNpZGVyYXRpb25zOioqCiMnIFRoZSB2ZXJpZnlpbmcgb2YgdGhlIHNhbXBsZXMgc2V4IGlzIHBlcmZvcm1lZCBieSBjb21wYXJpbmcgdGhlIGV4cHJlc3Npb24gbGV2ZWxzIG9mIAojJyB0aGUgZ2VuZXMgWElTVCBhbmQgVVRZLiBJbiBnZW5lcmFsLCBmZW1hbGVzIHNob3VsZCBoYXZlIGhpZ2ggWElTVCBhbmQgbG93IFVUWSBleHByZXNzaW9uLAojJyBhbmQgdmljZXZlcnNhIGZvciBtYWxlcy4gRm9yIGl0IHRvIHdvcmssIHRoZSBzYW1wbGUgYW5ub3RhdGlvbiBtdXN0IGhhdmUgYSBjb2x1bW4gY2FsbGVkICdTRVgnLAojJyB3aXRoIHZhbHVlcyBtYWxlL2ZlbWFsZS4gSWYgc29tZSBvdGhlciB2YWx1ZXMgZXhpc3QsIGUuZy4sIHVua25vd24sIHRoZXkgd2lsbCBiZSBtYXRjaGVkLiAKIycgVGhlIHByZWRpY3Rpb24gaXMgcGVyZm9ybWVkIHZpYSBhIGxpbmVhciBkaXNjcmltaW5hbnQgYW5hbHlzaXMgb24gdGhlIGxvZzIgY291bnRzLgoKIyBHZXQgc2V4IGNvbHVtbiBhbmQgcHJvY2VlZCBpZiBpdCBleGlzdHMKc2V4X2lkeCA8LSB3aGljaCgnU0VYJyA9PSB0b3VwcGVyKGNvbG5hbWVzKGNvbERhdGEob2RzKSkpKQppZihpc0VtcHR5KHNleF9pZHgpKXsKICBwcmludCgnU2V4IGNvbHVtbiBub3QgZm91bmQgaW4gc2FtcGxlIGFubm90YXRpb24nKQp9IGVsc2V7CiAgCiAgIyBWZXJpZnkgaWYgYm90aCBYSVNUIGFuZCBVVFkgd2VyZSBjb3VudGVkCiAgeGlzdF9pZCA8LSAnWElTVCcKICB1dHlfaWQgPC0gJ1VUWScKICAKICBpZihncmVwbCgnRU5TRzAnLCByb3duYW1lcyhvZHMpWzFdKSl7CiAgICB4aXN0X2lkIDwtICdFTlNHMDAwMDAyMjk4MDcnCiAgICB1dHlfaWQgPC0gJ0VOU0cwMDAwMDE4Mzg3OCcKICB9CiAgeGlzdF9pZHggPC0gZ3JlcCh4aXN0X2lkLCByb3duYW1lcyhvZHMpKQogIHV0eV9pZHggPC0gZ3JlcCh1dHlfaWQsIHJvd25hbWVzKG9kcykpCiAgCiAgaWYoaXNFbXB0eSh4aXN0X2lkeCkgfCBpc0VtcHR5KHV0eV9pZHgpKXsKICAgIHByaW50KCdFaXRoZXIgWElTVCBvciBVVFkgaXMgbm90IGV4cHJlc3NlZCcpCiAgfWVsc2V7CiAgICAKICAgIHNleF9jb3VudHMgPC0gY291bnRzKG9kcylbYyh4aXN0X2lkeCwgdXR5X2lkeCksIF0KICAgIHNleF9kdCA8LSBkYXRhLnRhYmxlKHNhbXBsZUlEID0gY29sbmFtZXMob2RzKSwgCiAgICAgICAgICAgICAgICAgICAgICAgICBYSVNUID0gY291bnRzKG9kcylbeGlzdF9pZHgsXSwgCiAgICAgICAgICAgICAgICAgICAgICAgICBVVFkgPSBjb3VudHMob2RzKVt1dHlfaWR4LF0pCiAgICBzZXhfZHQgPC0gbWVyZ2Uoc2V4X2R0LCBhcy5kYXRhLnRhYmxlKGNvbERhdGEob2RzKSlbLGMoMSwgc2V4X2lkeCksIHdpdGggPSBGXSwgc29ydCA9IEYpCiAgICBjb2xuYW1lcyhzZXhfZHQpIDwtIHRvdXBwZXIoY29sbmFtZXMoc2V4X2R0KSkKICAgIHNleF9kdFssIFNFWCA6PSB0b2xvd2VyKFNFWCldCiAgICBzZXhfZHRbaXMubmEoU0VYKSwgU0VYIDo9ICcnXQogICAgCiAgICAjIFRyYWluIG9ubHkgaW4gbWFsZS9mZW1hbGUgaW4gY2FzZSB0aGVyZSBhcmUgb3RoZXIgdmFsdWVzCiAgICB0cmFpbl9kdCA8LSBzZXhfZHRbU0VYICVpbiUgYygnZicsICdtJywgJ2ZlbWFsZScsICdtYWxlJyldCiAgICAKICAgIGxpYnJhcnkoIk1BU1MiKQogICAgbGRhIDwtIGxkYShTRVggfiBsb2cyKFhJU1QrMSkgKyBsb2cyKFVUWSsxKSwgZGF0YSA9IHRyYWluX2R0KQogICAgCiAgICBzZXhfZHRbLCBwcmVkaWN0ZWRfc2V4IDo9IHByZWRpY3QobGRhLCBzZXhfZHQpJGNsYXNzXQogICAgc2V4X2R0WywgbWF0Y2hfc2V4IDo9IFNFWCA9PSBwcmVkaWN0ZWRfc2V4XQogICAgdGFibGUoc2V4X2R0WywgLihTRVgsIHByZWRpY3RlZF9zZXgpXSkKICAgIAogICAgZyA8LSBnZ3Bsb3Qoc2V4X2R0LCBhZXMoWElTVCsxLCBVVFkrMSkpICsgCiAgICAgIGdlb21fcG9pbnQoYWVzKGNvbCA9IFNFWCwgc2hhcGUgPSBwcmVkaWN0ZWRfc2V4LCBhbHBoYSA9IG1hdGNoX3NleCkpICsgCiAgICAgIHNjYWxlX3hfbG9nMTAobGltaXRzID0gYygxLE5BKSkgKyBzY2FsZV95X2xvZzEwKGxpbWl0cyA9IGMoMSxOQSkpICsKICAgICAgc2NhbGVfYWxwaGFfbWFudWFsKHZhbHVlcyA9IGMoMSwgLjEpKSArIAogICAgICB0aGVtZV9jb3dwbG90KCkgKyBiYWNrZ3JvdW5kX2dyaWQobWFqb3IgPSAneHknLCBtaW5vciA9ICd4eScpICsgCiAgICAgIGFubm90YXRpb25fbG9ndGlja3Moc2lkZXMgPSAnYmwnKSArIAogICAgICBsYWJzKGNvbG9yID0gJ1NleCcsIHNoYXBlID0gJ1ByZWRpY3RlZCBzZXgnLCBhbHBoYSA9ICdNYXRjaGVzIHNleCcpCiAgICBwbG90KGcpCiAgICAKICAgIERUOjpkYXRhdGFibGUoc2V4X2R0W21hdGNoX3NleCA9PSBGXSwgY2FwdGlvbiA9ICdTZXggbWlzbWF0Y2hlcycpCiAgICAKICAgICMgV3JpdGUgdGFibGUKICAgIGZ3cml0ZShzZXhfZHQsIGdzdWIoJ29kc191bmZpdHRlZC5SZHMnLCAneGlzdF91dHkudHN2Jywgc25ha2VtYWtlQGlucHV0JG9kcyksIAogICAgICAgICAgIHNlcCA9ICdcdCcsIHF1b3RlID0gRikKICB9Cn0K