inst/shiny_admixtools/app.R
In uqrmaie1/admixtools: Inferring demographic history from genetic data
library(admixtools)
library(igraph)
library(htmlwidgets)
library(shiny)
library(shinythemes)
library(plotly)
library(tidyverse)
library(magrittr)
library(RColorBrewer)
library(shinyjs)
library(shinyBS)
library(shinyFiles)
library(shinyalert)
library(shinyWidgets)
library(shinydashboard)
#options(shiny.sanitize.errors = FALSE)
edgemul = 1000
dtstyle = 'font-size:70%'
pl = 15
do = list(pageLength = pl, lengthMenu = list(c(pl, -1), c(pl, 'All')),
buttons = c('copy', 'csv', 'excel'), scrollX = TRUE, dom = 'Blfrtip')
actionButton = function(...) actionBttn(..., style = 'unite', size = 'sm')
#cols = gg_color_hue(5, 97, 15)
cols = c("#FFF1F0", "#F7F8E4", "#E2FDF2", "#E8F9FF", "#FFF1FF")
mutfuns = c(`SPR leaves`='spr_leaves', `SPR all`='spr_all', `Swap leaves`='swap_leaves', `Move admixture edge`='move_admixedge_once', `Flip admix`='flipadmix_random', `Change root`='place_root_random', `Mutate n`='mutate_n')
box = function (..., title = NULL, footer = NULL, status = NULL, solidHeader = FALSE,
background = NULL, width = 6, height = NULL, collapsible = FALSE,
collapsed = FALSE)
{
boxClass <- "box"
if (solidHeader || !is.null(background)) {
boxClass <- paste(boxClass, "box-solid")
}
if (!is.null(status)) {
validateStatus(status)
boxClass <- paste0(boxClass, " box-", status)
}
if (collapsible && collapsed) {
boxClass <- paste(boxClass, "collapsed-box")
}
# if (!is.null(background)) {
# validateColor(background)
# boxClass <- paste0(boxClass, " bg-", background)
# }
style <- NULL
if (!is.null(height)) {
style <- paste0("background:", background, "; height: ", validateCssUnit(height))
}
titleTag <- NULL
if (!is.null(title)) {
titleTag <- h3(class = "box-title", title)
}
collapseTag <- NULL
if (collapsible) {
buttonStatus <- status %OR% "default"
collapseIcon <- if (collapsed)
"plus"
else "minus"
collapseTag <- div(class = "box-tools pull-right", tags$button(class = paste0("btn btn-box-tool"),
`data-widget` = "collapse", shiny::icon(collapseIcon)))
}
headerTag <- NULL
if (!is.null(titleTag) || !is.null(collapseTag)) {
headerTag <- div(class = "box-header", titleTag, collapseTag)
}
div(class = if (!is.null(width))
paste0("col-sm-", width), div(class = boxClass, style = if (!is.null(style))
style, headerTag, div(class = "box-body", style = paste('background:', background), ...), if (!is.null(footer))
div(class = "box-footer", footer)))
}
tt = c('data' = 'Load data here',
'dir' = 'Directory where precursors for f-statistics will be stored or have been stored.',
'popfilediv' = 'Select a PLINK or EIGENSTRAT sample metadata textfile.',
'graphfilediv' = 'A graph in ADMIXTOOLS format, or as a two column edge list.',
'genofile1' = 'Select a Packedancestrymap .geno file to extract data from. This step should only be necessary once. After that, the extracted data can be used for any further analyses.',
'adjpopdiv' = 'Assign individuals to populations. Population names can be changed, but should match population names in the graph file, if one is provided. Data will be read from disk and f2-statistics will be computed when switching to another sidebar item.',
'bookmarkdiv' = 'Click this to save the current session. This will write files to disk and generate a URL which can be used to restore the state of the progam.',
'qpgraph_fit' = 'Click this to see qpgraph fit statistics',
'qpgraph_similar_minus1' = 'Fit all graphs that result from removing one admixture edge from the current graph',
'qpgraph_similar_minusplus' = 'Fit all graphs that result from removing one admixture edge from the current graph and adding back one admixture edge',
'qpgraph_similar_plus1' = 'Fit all graphs that result from adding one admixture edge to the current graph',
'qpgraph_similar_decomposed' = 'Fit all trees that result from removing one edge from each admixture node',
'qpgraph_similar_treeneighbors' = 'Fit all trees that result from removing one edge from each admixture node and then reattaching each edge to each other edge',
'qpgraph_similar_flipadmix' = 'Fit all graphs that result from flipping the direction of each admixture edge (unless this would introduce a cycle)',
'qpgraph_update' = 'Set the current graph to the displayed graph',
'qpgraph_revert' = 'Display the original graph',
'addleaf' = 'Choose a population which should be added to the current graph',
'qpgraph_add' = 'Fit all graphs that result from adding this population to each edge',
'optim_run' = 'Run optimization. The current graph is modified in a number of different ways, and each modified graph is evaluated. Can be run iteratively across several generations.',
'optim_ngen' = 'Number of generations. In each generation the best fitting graphs are selected for the next generation.',
'optim_sec' = 'Minimum number of seconds to run',
'mutfuns' = 'SPR leaves: Subree prune and regraft leaf edges\n
SPR all: Subree prune and regraft any edges\\n
Swap leaves: Swap two random leaves\\\n
Move admixture edge: Shift one admixture edge\\\\n
Flip admix: Flip the direction of one admixture edge\\\\\n
Change root: Place the root of the tree on a random edge\\\\\\n
Mutate n: Combine two or more of the other mutation functions',
'econstraints_update_div' = 'Set limits to minimum and maximum edge lengths',
'aconstraints_update_div' = 'Set limits to admixture weights',
'qpgraphbin' = 'Choose the qpGraph executable on your computer',
'qpadmbin' = 'Choose the qpAdm executable on your computer',
'qpgraph_genofile' = 'Choose the .geno or .bed file for qpGraph to use',
'qpadm_genofile' = 'Choose the .geno or .bed file for qpAdm to use',
'usef3precomp' = 'Differences between the original and R version of qpGraph can be due to the estimation of f-statistics, or due to the fitting step. Click this to use the f-statistics generated in the original qpGraph to use in qpGraph in R.',
'init_resample' = 'Re-evaluate graph with different random intializations of starting weights. This is always done by default, but only the parameters for the starting weights that minimize the score are shown. Here, the table will show all estimated graph parameters (score, edge lengths, admixture weights) for all sets of random starting parameters.',
'ind_resample' = 'Re-evaluate the graph leaving out each sample at a time. Samples are only left out if they are not the only sample in their respective population. This can take long, as data has to be read and f2 computed separately for each sample that is being left out.',
'snp_resample' = 'Re-evaluate the graph leaving out SNP blocks. With jackknife, each SNP block is left out at a time. With bootstrap, SNP blocks are resampled with replacement to match the original number of SNP blocks (usually around 700). The number of random resamplings can be specified.',
'clear_edge' = 'Clear the current selection',
'delete_edge' = 'Delete the selected edge or node. Has to be an admixture edge or leaf node.',
'add_edge' = 'Add an admixture edge by connecting two selected non-admixture edges. Fails if it introduces cycles or attempts to connect admixture edges.',
'randgraph' = 'Generate a random admixture graph. Uses the current populations, a specified number of admixture nodes, and a specified outgroup.',
'lsqmode' = 'Least squares mode. By default, the computation of the graph likelihood score uses the inverse of the f3-statistic covariance matrix. When fitting many populations, this inverse can be unstable. Least squares mode sets the offdiagonal elements of the f3-statistic covariance matrix to zero, which makes the matrix inversion stable, but can introduce bias.',
'lambdascale' = 'Scale the estimated f2-statistics by a constant factor.',
'seed' = 'Set a random seed to always pick the same initial weights and to make the results reproducible.',
'multiprocess' = 'Enable parallel evaluation of graphs in multiple R processes. Faster, but may not work on all operating systems. See "A Future for R: A Comprehensive Overview" for more details.',
'collapse_edges' = 'Collapse all nodes which are separated by less then a certain threshold of estimated drift.',
'f2corr' = 'Untick this to compute f-statistics without bias correction factor.',
'addpop' = 'Add a population',
'removepop' = 'Remove a population',
'showseledge' = 'Selected edges or nodes are shown here',
'downloadgraph' = 'Download the fit of the current graph as a .tsv file.',
'numstart' = 'How often optimization should be initialized with random weights',
'qpgraph_diag' = 'Factor to add to diagonal elements of f3-statistic covariance matrix before inversion. Increase this value if optimization doesn\'t converge',
'qpadm_fit' = 'Return to the qpadm fit statistics',
'qpadm_diag' = 'Factor to add to diagonal elements of f4-statistic covariance matrix before inversion. Increase this value if optimization doesn\'t converge',
'qpadm_constrained' = 'Click this to constrain the admixture weights to be between 0 and 1',
'data_load' = 'Only a placeholder at the moment. Data is loaded when switching from the "Data" tab to another tab',
'b' = '',
'c' = '',
'f2' = 'These are the estimated f2 statistics',
'f3' = 'These are the estimated and fitted f3 statistics',
'opt' = 'This table contains information about the admixture weight optimization step for each iteration, such as the randomly drawn initial weights, as well as the final score ("value", the quantity that is being minimized) and convergence parameters from the R "optim" function. Run "Resample / Initial weights" to see how other estimated graph parameters depend on the initial weights.')
format_table = function(dat) {
if(is.null(dat)) return(tibble())
print('ft1')
print(names(dat))
popcols = which(str_sub(names(dat), 1, 3) == 'pop')
names = unique(c(as.matrix(dat[,popcols])))
print('ft2')
nam = shortest_unique_prefixes(names, 3)
pcols = which(names(dat) %in% c('p', 'p.value'))
print('ft3')
numcols = setdiff(which(sapply(dat, class) == 'numeric'), pcols)
print('ft4')
if(length(popcols) > 0) dat %<>% mutate_at(popcols, ~recode(., !!!purrr::set_names(nam, names)))
if('graph' %in% names(dat)) dat %<>% select(-graph)
if('edges' %in% names(dat)) dat %<>% select(-edges)
print('ft5')
out = dat %>%
mutate_at(numcols, ~round(., max(2, 2 - ceiling(log10(max(abs(na.omit(.)))))))) %>%
mutate_at(pcols, format.pval) %>%
rename_at(vars(pcols), ~'p.value')
out
}
aa = '.shiny-notification { height: 100px; width: 400px; position: fixed; top: 20%; left: calc(50% - 200px);;}'
ui = function(request) {
tagList(useShinyjs(),
useShinyalert(),
chooseSliderSkin('Modern', color = 'black'),
bsTooltip('tooltips', 'Click this to turn tooltips on or off. Needs to be reactivated after switching tabs.'),
dashboardPage(
dashboardHeader(title = 'ADMIXTOOLS 2', titleWidth = 300,
tags$li(class = 'dropdown', actionLink('tooltips', label = '', icon = icon('question'))),
tags$li(class = 'dropdown', id = 'bookmarkdiv', actionLink('._bookmark_', label = '', icon = icon('save')))),
dashboardSidebar(tags$head(tags$style(HTML(aa))), width = 300,
sidebarMenu(id = 'navbar',
menuItem('Data', tabName = 'data', expandedName = 'data', id = 'datax', icon = icon('database'), startExpanded = TRUE,
menuItem('Extract settings', tabName = 'extract_settings',
checkboxInput('fix_populations', 'Fix populations', TRUE),
numericInput('max_miss', 'max missing', value = 0, step = 0.01, min = 0, max = 1),
splitLayout(
numericInput('min_maf', 'min MAF', value = 0, step = 0.01, min = 0, max = 0.5),
numericInput('max_maf', 'max MAF', value = 0.5, step = 0.01, min = 0, max = 0.5)),
radioButtons('trans_extract', 'Mutations', choices = c('both', 'only transitions', 'only transversions')),
fileInput('keepsnps', NULL, buttonLabel = 'SNP list'),
numericInput('maxmem', 'max RAM in GB', value = 15, step = 1, min = 0)
)),
menuItem('f2', tabName = 'f2', expandedName = 'f2', icon = icon('dice-two'),
menuItem('Options', tabName = 'f2_options')),
menuItem('f3', tabName = 'f3', expandedName = 'f3', icon = icon('dice-three'),
menuItem('Options', tabName = 'f3_options')),
menuItem('f4', tabName = 'f4', expandedName = 'f4', icon = icon('dice-four'),
menuItem('Options', tabName = 'f4_options')),
menuItem('qpAdm', tabName = 'qpadm', expandedName = 'qpadm', icon = icon('balance-scale'),
actionLink('qpadm_fit', 'Fit'),
menuItem('Compare', href = 'qpadm_comparison', tabName = 'qpadm_comparison',
splitLayout(
shinyFilesButton('qpadmbin', 'qpAdm bin', 'Select qpAdm bin', FALSE),
shinyFilesButton('qpadm_genofile', 'Geno file', 'Select Packedancestrymap geno file', FALSE)),
p(HTML('<center>or</center>')),
fileInput('qpadm_out', NULL, buttonLabel = 'qpAdm output'),
hr(),
actionButton('run_qpadm', 'Run')),
menuItem('Resample', tabName = 'qpadm_resample'),
menuItem('Options', tabName = 'qpadm_options',
numericInput('qpadm_diag', 'diag', value = 0.0001, step = 0.0001),
checkboxInput('qpadm_constrained', 'Constrain weights'),
checkboxInput('qpadm_results', 'Show results', value = TRUE),
checkboxInput('qpadm_rotate', 'Rotate left', value = FALSE))),
menuItem('qpGraph', tabName = 'qpgraph', expandedName = 'qpgraph', id = 'qpgraph', icon = icon('project-diagram'),
checkboxInput('evaluate_graph', 'Evaluate graph', FALSE),
menuItem('Fit',
splitLayout(
shiny::actionButton('qpgraph_update', 'Update graph'),
shiny::actionButton('qpgraph_revert', 'Revert graph')),
#actionButton('options_update', 'Update'),
splitLayout(numericInput('qpgraph_diag', 'diag', value = 0.0001, step = 0.0001),
numericInput('numstart', '# init', value = 10, min = 1)),
splitLayout(numericInput('lambdascale', 'lambdascale', value = 1, step = 0.001),
numericInput('seed', 'Random seed', value = NULL)),
splitLayout(checkboxInput('lsqmode', 'lsqmode'),
checkboxInput('return_fstats', 'return_fstats'))),
menuItem('Modify', expandedName = 'qpgraph_modify', id = 'qpgraph_modify', tabName = 'qpgraph_modify',
hr(),
p('Selected'),
div(verbatimTextOutput('showseledge', placeholder = TRUE), tags$head(tags$style('#showseledge{max-height: 300px; background: white}'))),
splitLayout(
shiny::actionButton('clear_edge', 'Clear'),
shiny::actionButton('delete_edge', 'Delete'),
shiny::actionButton('add_edge', 'Add')),
uiOutput('qpgraph_add'),
hr(),
p('Randomize graph'),
uiOutput('nadmix'),
#splitLayout(checkboxInput('fixoutgroup', 'Fix outgroup', value = TRUE),
actionButton('randgraph', 'Randomize'),
hr()),
menuItem('Optimize', tabName = 'qpgraph_optim',
actionButton('optim_run', 'Run'),
numericInput('optim_ngen', '# generations', 10, 1, 1000),
numericInput('optim_sec', 'min seconds', 10, 1, 1000),
checkboxInput('optim_initrand', 'Start randomly'),
checkboxGroupInput('mutfuns', 'Graph modifications', choices = mutfuns, selected = mutfuns)),
menuItem('Constrain', tabName = 'qpgraph_constraints',
menuItem('Drift', tabName = 'qpgraph_econstraints',
uiOutput('econstraints')),
menuItem('Admixture', tabName = 'qpgraph_aconstraints',
uiOutput('aconstraints'))),
menuItem('Compare', tabName = 'qpgraph_comparison',
splitLayout(
shinyFilesButton('qpgraphbin', 'qpGraph bin', 'Select qpGraph bin', FALSE),
shinyFilesButton('qpgraph_genofile', 'Geno file', 'Select Packedancestrymap geno file', FALSE)),
checkboxInput('usef3precomp', 'Use qpGraph fstats in R'),
p(HTML('<center>or</center>')),
fileInput('qpgraph_out', NULL, buttonLabel = 'qpGraph output'),
hr(),
actionButton('run_qpgraph', 'Run')),
menuItem('Resample', tabName = 'qpgraph_resample',
menuItem('Initial weights', tabName = 'qpgraph_resample_init',
sliderInput('numstart2', '# init', 1, 100, 10),
actionButton('init_resample', 'Run')),
menuItem('Individuals', tabName = 'qpgraph_resample_individuals',
actionButton('ind_resample', 'Run')),
menuItem('SNPs', tabName = 'qpgraph_resample_snps',
checkboxInput('boot', 'Bootstrap'),
conditionalPanel('input.boot',
sliderInput('bootnum', '# boot', 0, 1000, 100)),
actionButton('snp_resample', 'Run'))),
menuItem('Load', fileInput('graphfile', NULL, placeholder = '', buttonLabel = 'Graph file')),
menuItem('Save', tabName = 'qpgraph_download', expandedName = 'qpgraph_download',
div(radioButtons('downloadgraph_format', 'Graph format',
choices = c(`ADMIXTOOLS`='admixtools', `Edge list`='edgelist', `DOT`='dot'),
selected = 'admixtools'),
downloadButton('downloadgraph', 'Save graph'))),
menuItem('Plot options',
checkboxGroupInput('plotopt', '',
choices = c(`Highlight unidentifiable` = 'highlight_unidentifiable',
`Reorder edges`='reorder_edges',
`Shift edges down`='shift_down',
`Collapse edges`='collapse_edges'),
selected = c('shift_down')),
sliderInput('collapse_threshold', 'Log10 collapse threshold', -6, 2, -3, step = 0.1))),
menuItem('Options', tabName = 'options', expandedName = 'options', id = 'options', icon = icon('cogs'),
checkboxInput('multiprocess', 'multiprocess', value = FALSE),
checkboxInput('f2corr', 'f2 bias correction', value = TRUE),
checkboxInput('afprod', 'allele frequency products', value = FALSE)))),
dashboardBody(uiOutput('dashboardbody'), tags$head(tags$style(HTML('.skin-black .sidebar .shiny-download-link a { color: #444; } .content {background-color: white } .main-header .logo {font-family: "GillSans", "Skia", "Avenir-Medium", "Herculanum", "Hiragino Maru Gothic Pro", "Maru Gothic Pro", "Hiragino", "Comic Sans MS"; font-weight: normal;};')))),
skin = 'black'
))
}
server = function(input, output, session) {
showLog()
dp = '~/mnt/o2/projects/admixprograms/indpairs_v42.1/'
dp = '~/Downloads/testinds/'
dp = '/Users/robert/Downloads/countdata_eas5/'
dp = '~/Downloads/countdata_iosifgraph/'
dp = ''
#dp = '~/Documents/countdata_afafam/'
#dp = '~/Documents/countdata_pavel2/'
global = reactiveValues(countdir = dp,
allinds = list.dirs(paste0(dp, '/pairs'),F,F),
poplist = list())
global$graph = NULL
home = normalizePath('~')
volumes = getVolumes()
#volumes = expr(c(home = normalizePath('~'), getVolumes()() %>% purrr::set_names()))
shinyDirChoose(input, 'dir', roots = volumes)
shinyFileChoose(input, 'genofile1', roots = volumes)
shinyFileChoose(input, 'qpgraph_genofile', roots = volumes)
shinyFileChoose(input, 'qpadm_genofile', roots = volumes)
shinyFileChoose(input, 'qpgraphbin', roots = volumes)
shinyFileChoose(input, 'qpadmbin', roots = volumes)
#global$show_qpadm_rotate = FALSE
#global$show_qpadm_results = TRUE
shinyjs::disable('qpgraph_update')
shinyjs::disable('qpgraph_revert')
output$dirout = renderText({global$countdir})
output$dirout = renderText({str_replace(global$countdir, '/Users/robert', '')})
output$show_popadjust = renderText({as.numeric(show_popadjust())})
output$show_extract = renderText({as.numeric(show_extract())})
output$show_indselect = renderText({as.numeric(show_indselect())})
output$genofile1out = renderText({parseFilePaths(eval(volumes), input$genofile1)$datapath %>% str_replace('\\.geno$', '') %>% str_replace('/Users/robert', '')})
observeEvent(input$dir, {
#if (!"path" %in% names(dir())) return()
#volumes = eval(volumes)
print('volumes')
print(volumes)
print('input$dir')
print(c(input$dir, class(input$dir)))
print('xxx')
print(parseDirPath(volumes, input$dir))
print(str(parseDirPath(volumes, input$dir)))
print(as.character(parseDirPath(volumes, input$dir)))
global$countdir = parseDirPath(volumes, input$dir)
global$iscountdata = 'indivs' %in% list.dirs(global$countdir,F,F)
global$isf2data = !'indivs' %in% list.dirs(global$countdir,F,F) &
'block_lengths_f2.rds' %in% list.files(global$countdir)
print('global$iscountdata')
print(global$iscountdata)
if(global$iscountdata) {
global$allinds = list.dirs(paste0(global$countdir, '/pairs'),F,F)
} else if(global$isf2data) {
global$allinds = list.dirs(global$countdir,F,F)
global$poplist = global$allinds %>% purrr::set_names(global$allinds)
print('global$allinds')
print(global$allinds)
} else {
}
})
observeEvent(input$textdirinput, {
nam = ifelse(!is.null(input$textdirinput) && dir.exists(input$textdirinput), 'Load', 'Create')
updateActionButton(session, "textdir", label = nam)
})
observeEvent(input$textdir, {
justcreated = FALSE
if(!dir.exists(input$textdirinput)) {
#shinyalert('Could not find directory!', "Please write the path of a directory into the box. This should be an empty directory if you haven't extracted f2-statistics/counts yet, or otherwise the directory with extracted data.")
dir.create(input$textdirinput)
justcreated = TRUE
#return()
}
global$countdir = normalizePath(input$textdirinput)
global$iscountdata = 'indivs' %in% list.dirs(global$countdir,F,F)
global$isf2data = !'indivs' %in% list.dirs(global$countdir,F,F) &
'block_lengths_f2.rds' %in% list.files(global$countdir)
print('global$iscountdata')
print(global$iscountdata)
if(global$iscountdata) {
global$allinds = list.dirs(paste0(global$countdir, '/pairs'),F,F)
} else if(global$isf2data) {
global$allinds = list.dirs(global$countdir,F,F)
global$poplist = global$allinds %>% purrr::set_names(global$allinds)
print('global$allinds')
print(global$allinds)
} else {
if(!justcreated) shinyalert("Directory does not seem to have f2 data!")
}
})
observeEvent(input$navbar, {
print('navbar detected')
print(input$navbar)
})
tolisten = reactive({
list(input$sidebarItemExpanded, input$qpadm_rotate, input$qpadm_results)
#!is.null(input$sidebarItemExpanded) && (is.null(global$exp) || input$sidebarItemExpanded != global$exp)
})
get_bod = eventReactive(tolisten(), {
print('expanded!')
exp = input$sidebarItemExpanded
print(exp)
#if(is.null(exp) || !is.null(global$exp) && exp == global$exp) return(global$bod)
if(is.null(exp)) exp = global$exp
global$exp = exp
if(exp == 'options') {
} else if(exp == 'data') {
print('navbar == data!')
global$databod = get_loaddata()
global$bod = global$databod
} else {
nam = map(paste0('pop', seq_len(length(global$poplist))), ~input[[.]])
if(is.null(nam[[1]])) {
shinyalert('Global poplist is NULL or has not been processed correctly!', global$poplist)
return()
}
global$poplist = map(nam, ~input[[.]]) %>% purrr::set_names(nam)
print('switch away from data')
if(length(unlist(global$poplist)) == 0) {
shinyalert('No samples selected!', '')
return()
}
if(exp == 'f2') {
print('switch to f2')
global$bod = fluidRow(
column(1, div(radioGroupButtons('fst', choices = c('f2', 'fst')))),
column(11, tabsetPanel(tabPanel('Pop pairs', plotlyOutput(outputId = 'f2heatmap', height = '800', width='auto')),
tabPanel('Ind pairs', plotlyOutput(outputId = 'f2heatmap_indiv', height = '800', width='auto')),
tabPanel('Pop table', dto('f2stats')))))
} else if(exp == 'f3') {
print('switch to f3')
global$bod = fluidRow(
column(12, tabsetPanel(tabPanel('Pop table', dto('f3stats')))))
} else if(exp == 'f4') {
print('switch to f4')
choices = names(global$poplist)
global$f4_left = div(map(1:4, ~selectizeInput(paste0('f4pops', .),
paste0('Population ', .),
choices = choices,
multiple = TRUE,
selected = choices[seq(.*(.-1)/2+1, .*(.-1)/2+.)])),
checkboxInput('f4_allowdup', 'Allow duplicates'),
sliderInput('f4_se', 'Standard errors', min = 0, max = 5, value = 0, step = 0.01))
#checkboxInput('f4_se', 'Show standard errors'))
popsinuse = unique(c(input$f4pops1, input$f4pops2, input$f4pops3, input$f4pops4))
# update
map(1:4, ~{
fp = paste0('f4pops', .)
sel = input[[fp]]
updateSelectizeInput(session, fp, selected = sel,
choices = union(sel, setdiff(names(global$poplist), popsinuse)))})
global$bod = fluidRow(
column(3, global$f4_left),
column(9, tabsetPanel(tabPanel('Plot', plotlyOutput(outputId = 'f4plot', height = '800', width='auto')),
tabPanel('Plot2', plotlyOutput(outputId = 'f4plot2', height = '800', width='auto')),
tabPanel('Table', dto('f4stats')))))
} else if(exp == 'qpadm') {
print('switch to qpAdm')
choices = names(global$poplist)
if(is.null(input$qpadmpops1) || is.null(input$qpadmpops2) || is.null(input$qpadmpops3)) {
global$qpadmpops1 = choices[2:3]
global$qpadmpops2 = choices[4:6]
global$qpadmpops3 = choices[1]
} else {
global$qpadmpops1 = input$qpadmpops1
global$qpadmpops2 = input$qpadmpops2
global$qpadmpops3 = input$qpadmpops3
}
#global$qpadm_leftpanel = qpadm_leftpanel()
global$bod = fluidRow(
column(3, qpadm_leftpanel()),
column(9, qpadm_rightpanel()))
} else if(exp == 'qpgraph') {
print('switch to qpGraph')
if(is.null(global$graph) || !all(get_leafnames(global$graph) %in% names(global$poplist))) {
numadm = 2
if(!is.null(global$graph)) {
shinyalert('Generating random graph because populations don\'t match',
setdiff(get_leafnames(global$graph), names(global$poplist)))
numadm = numadmix(isolate(global$graph))
}
op = if(is.null(input$outpop) || input$outpop == '< undefined >') NULL else input$outpop
global$graph = random_admixturegraph(na.omit(names(global$poplist)), numadm, TRUE, outpop = op)
}
#global$qpg_right = qpg_right_fit()
#global$qpg_right_show = 'fit'
#global$qpg_right = qpg_right()
global$bod = fluidRow(
column(8, plotlyOutput(outputId = 'graphplot', height = '800', width='auto')),
column(4, tabsetPanel(tabPanel('Fit', tabsetPanel(tabPanel('f2', dto('f2')),
tabPanel('f3', dto('f3')),
tabPanel('f4', dto('f4')),
tabPanel('opt', dto('opt')))),
tabPanel('Similar', tabsetPanel(tabPanel('-1', dto('minus1')),
tabPanel('+1', dto('plus1')),
tabPanel('flip', dto('flipadmix')))),
tabPanel('History', dto('history')),
tabPanel('Compare', plotlyOutput('graphcomparison')),
tabPanel('Resample', tabsetPanel(tabPanel('Starting values', dto('initresample')),
tabPanel('Individuals', dto('indresample')),
tabPanel('SNPs', dto('snpresample')))),
id = 'qpgraph_tabset')))
} else {
shinyalert('Error', 'implement tab')
}
}
global$bod
})
observeEvent(global$graph, {
print('observe global$graph')
updateSelectInput(session, 'addleaf',
choices = setdiff(names(global$poplist), get_leafnames(global$graph)))
})
get_graph = reactive({
print('get graph')
global$graph
})
# observeEvent(global$qpg_right, {
# print('global$qpg_right change detected!')
# global$qpgraphbod = fluidRow(
# column(8, plotlyOutput(outputId = 'graphplot', height = '800', width='auto')),
# column(4, global$qpg_right))
# })
# observeEvent(global$databod, {global$bod = global$databod})
# observeEvent(global$qpadmbod, {global$bod = global$qpadmbod})
observeEvent(global$qpgraphbod, {print('global$qpgraphbod change detected!'); global$bod = global$qpgraphbod})
observeEvent(global$bod, {print('global$bod change detected!'); b = global$bod; global$bod = b})
observeEvent(input$navbar, {
print('navbar selected')
print(input$navbar)
})
observeEvent(input$addpop, {
pops = names(global$poplist)
i = 0
repeat({i=i+1; newpop = paste0('pop', i); if(!newpop %in% pops) break})
global$poplist[newpop] = NA
})
observeEvent(input$removepop, {
global$poplist = head(global$poplist, -1)
})
observeEvent(input$popfile, {
inddat = read_table2(input$popfile$datapath, col_names = FALSE, col_types = cols())
if(ncol(inddat) == 3) colnames(inddat) = c('ind', 'sex', 'pop')
else if(ncol(inddat) == 6) colnames(inddat) = c('pop', 'ind', 'p1', 'p2', 'sex', 'pheno')
else shinyalert('Error', paste0('sample file has ',ncol(inddat),' columns!'))
global$poplist = inddat %>%
select(ind, pop) %$% split(ind, factor(pop, levels = unique(pop)))
if(length(global$allinds) == 0) { # f2data dir
global$allinds = unlist(global$poplist)
}
print('observe')
imap(names(global$poplist), ~{
bttn = paste0('group', .y)
print(bttn)
observeEvent(input[[bttn]], {
isolate({
if(input[[bttn]] == 0) return()
poplist = global$poplist
#poplistbak = global$poplistbak
countdir = global$countdir
inds = poplist[[.x]]
#if(is.null(global$poplistbak[[.x]])) {
isgrouped = length(inds) == 1 && admixtools:::is_group(countdir, inds)
if(!isgrouped) {
if(.x %in% global$allinds) {
shinyalert('Error', 'Group name cannot be identical to an existing sample or group name!')
return()
}
print('add group')
withProgress(message = paste0('grouping ', length(inds), ' samples...'), {
group_samples(global$countdir, inds, .x, overwrite = TRUE)
})
global$allinds = union(global$allinds, .x)
#global$poplistbak[[.x]] = inds
global$poplist[[.x]] = .x
updateSelectizeInput(session, .x, selected = .x, choices = global$allinds)
#updateButton(session, bttn, 'ungroup')
} else {
print('delete group')
withProgress(message = paste0('ungrouping ', length(inds), ' samples...'), {
global$poplist[[.x]] = delete_groups(global$countdir, .x)
})
#global$poplist[[.x]] = global$poplistbak[[.x]]
#global$poplistbak[[.x]] = NULL
global$allinds = setdiff(global$allinds, .x)
updateSelectizeInput(session, .x, selected = global$poplist[[.x]], choices = global$allinds)
#updateButton(session, bttn, 'group')
}
})
})
})
})
observeEvent(input$graphfile, {
gf = input$graphfile$datapath
print('read graphfile')
print(gf)
colnums = read_lines(gf) %>% str_split('\\s+') %>% map_dbl(length)
if(length(unique(colnums)) > 1) {
print('admixtools format')
global$graph = parse_qpgraph_graphfile(gf) %>% graph_from_edgelist()
} else {
print('two column format')
global$graph = read_table2(gf, col_types = cols()) %>% select(1:2) %>%
as.matrix %>% graph_from_edgelist()
}
print('read graphfile done')
})
observeEvent(input$qpgraph_out, {
global$qpgraphout = parse_qpgraph_output(input$qpgraph_out$datapath)
#global$qpg_right = plotlyOutput('graphcomparison')
})
observeEvent(input$randgraph, {
print('randg')
op = if(is.null(input$outpop) || input$outpop == '< undefined >') NULL else input$outpop
global$graph = random_admixturegraph(get_leafnames(global$graph), input$nadmix, TRUE, outpop = op)
})
observeEvent(input$qpgraph_update, {
print('qpgraph_update')
global$graph = global$selgraph
fit = get_fit()
global$edges = fit$edges
global$score = fit$score
global$worst_residual = fit$worst_residual
add_graph_to_history(global)
global$tempgraph = FALSE
shinyjs::disable('qpgraph_update')
shinyjs::disable('qpgraph_revert')
})
observeEvent(input$qpgraph_revert, {
print('qpgraph_revert')
global$selgraph = global$graph
fit = get_fit()
global$edges = fit$edges
global$score = fit$score
global$worst_residual = fit$worst_residual
add_graph_to_history(global)
global$tempgraph = FALSE
shinyjs::disable('qpgraph_update')
shinyjs::disable('qpgraph_revert')
print(global$score)
})
observeEvent(input$navbar, {
print('navbar:')
print(input$navbar)
})
observeEvent(input$optim_run, {
numgen = input$optim_ngen
stop_sec = input$optim_sec
diag = input$qpgraph_diag
lambdascale = as.numeric(input$lambdascale)
f2blocks = get_f2blocks()
pops = dimnames(f2blocks)[[1]]
f2blocks = f2blocks[pops, pops,]
outpop = get_outpop(global$graph)
nadmix = numadmix(global$graph)
g = NULL
if(!input$optim_initrand) g = global$graph
print('running opt...')
print(pops)
print(dim(f2blocks))
print(paste(numgen, nadmix))
mutfuns = sapply(input$mutfuns, get)
withProgress(message = paste('Running find_graphs...\nCheck R console to see progress'), {
# opt_results = find_graphs(f2blocks, outpop = outpop, numrep = numrep,
# numgraphs = numgraphs, numgen = numgen, numsel = numsel,
# numadmix = nadmix, initgraph = g, mutfuns = mutfuns, diag = diag)
opt_results = find_graphs(f2blocks, outpop = outpop, stop_gen = numgen, stop_sec = stop_sec,
numadmix = nadmix, initgraph = g, mutfuns = mutfuns, diag = diag)
})
print(opt_results)
opt_results %<>% filter(!is.na(score))
print('opt done')
if(nrow(opt_results) == 0) {
shinyalert('No valid graphs in output!', 'Try increasing the "diag" parameter.')
} else {
winner = opt_results %>% slice_min(score) %$% graph[[1]]
oldgraph = global$graph
oldscore = global$score
oldedges = global$edges
oldwr = global$worst_residual
withProgress(message = 'Fitting best graph...', {
global$graph = winner
fit = get_fit()
})
print('fit')
print(str(fit))
if(fit$score < oldscore) {
shinyalert('Better graph found!', paste('Old score: ', round(oldscore, 2),
'\nNew score: ', round(fit$score, 2)))
} else {
global$graph = oldgraph
global$score = oldscore
global$edges = oldedges
global$worst_residual = oldwr
add_graph_to_history(global)
fit = get_fit()
shinyalert('No better graph found!')
}
print('graph')
print(global$graph)
}
global$seed = NULL
})
qpgraphfun = reactive({
req(input$evaluate_graph)
print('qpgraphfun')
seed = isolate(as.numeric(input$seed))
if(is.na(seed)) seed = global$seed
if(input$usef3precomp && !is.null(global$precomp)) f3precomp = global$precomp
else f3precomp = NULL
#isolate({
numstart = as.numeric(input$numstart)
lambdascale = as.numeric(input$lambdascale)
if(is.null(lambdascale) || length(lambdascale) == 0) lambdascale = 1
diag = input$qpgraph_diag
lsqmode = input$lsqmode
return_fstats = input$return_fstats
print(paste('qpgraphfun:', numstart, seed, diag, lsqmode, lambdascale))
function(x, y, ...) {
args = list(...)
if(!'numstart' %in% names(args)) args[['numstart']] = numstart
args = c(list(x, y), args, diag = diag, lambdascale = lambdascale, lsqmode = lsqmode,
seed = seed, f3precomp = list(f3precomp), return_fstats = return_fstats)
do.call(quietly(qpgraph), args)$result
}
#})
})
get_fit = reactive({
print('get_fit')
if(!input$evaluate_graph) return(list(edges=global$edges, score=global$score, worst_residual=global$worst_residual))
f2blocks = get_f2blocks()
g = global$graph
input$qpgraph_revert
input$qpgraph_update
input$aconstraints_update
input$econstraints_update
global$qpgraph_ranges = NULL
global$qpgraph_scorerange = NULL
req(g, f2blocks)
replace_null = function(x, y) ifelse(is.null(x), y, x)
edges = g %>% igraph::as_edgelist() %>% as_tibble(.name_repair = ~c('from', 'to')) %>%
group_by(to) %>% mutate(type = ifelse(n() > 1, 'admix', 'normal')) %>% ungroup %>%
mutate(eid = paste(from, to, sep = ' -> '), eid2 = str_remove_all(eid, '\\.|>| '))
isolate({
if(isTRUE(global$useconstraints)) {
edges %<>% mutate(lower = map_dbl(eid2, ~replace_null(input[[.]][1], 0)),
upper = map_dbl(eid2, ~replace_null(input[[.]][2], 1e9))) %>%
mutate(upper = ifelse(type == 'normal', upper/edgemul, upper),
lower = ifelse(type == 'normal', lower/edgemul, lower)) %>% select(-type, -eid2)
#print(edges)
}
})
withProgress(message = 'Fitting graph...', {
fit = qpgraphfun()(f2blocks, edges)
})
global$edges = fit$edges
global$score = fit$score
global$worst_residual = fit$worst_residual
add_graph_to_history(global)
fit
})
get_history = reactive({
print('get_history')
isolate({
hist = global$history %>% purrr::transpose() %>% as_tibble %>%
transmute(graph, edges,
nadmix = as.numeric(nadmix),
score = as.numeric(score),
worst_residual = map_dbl(worst_residual, ~ifelse(is.null(.), NA, .))) %>%
mutate(n = 1:n(), .before = 1) %>% arrange(-n)
})
hist %>% mutate(x = list(get_fit())) %>% select(-x)
})
get_pdat = reactive({
print('get_pdat')
if(is.null(global$edges)) return()
print('get_pdat2')
init_plotdata()
})
get_f2blocks = reactive({
print('get f2')
poplist = global$poplist
req(poplist, global$countdir)
if(length(poplist) == 0) return()
pops = rep(names(poplist), sapply(poplist, length))
inds = NULL
if(global$iscountdata) inds = unlist(poplist)
withProgress(message = 'Reading data and computing f2...', {
f2blocks = f2_from_precomp(global$countdir, inds = inds, pops = pops,
apply_corr = input$f2corr, afprod = input$afprod)
})
print('get f2 done')
f2blocks
})
get_fstblocks = reactive({
print('get fst')
poplist = global$poplist
req(poplist, global$countdir)
if(length(poplist) == 0) return()
pops = rep(names(poplist), sapply(poplist, length))
inds = NULL
if(global$iscountdata) inds = unlist(poplist)
withProgress(message = 'Reading data and computing fst...', {
f2blocks = f2_from_precomp(global$countdir, inds = inds, pops = pops,
apply_corr = input$f2corr, afprod = input$afprod, fst = TRUE)
})
print('get fst done')
f2blocks
})
get_f2dat_indiv = reactive({
print('get f2 indiv')
poplist = global$poplist
req(poplist, global$countdir)
if(length(poplist) == 0) return()
withProgress(message = 'Reading data and computing f2...', {
f2dat = f2_from_precomp(global$countdir, inds = unlist(poplist), pops = unlist(poplist),
return_array = FALSE, afprod = input$afprod)
})
print('get f2 indiv done')
f2dat
})
init_plotdata = reactive({
print('init')
input$qpgraph_update
input$qpgraph_revert
input$evaluate_graph
graph = global$edges
if('igraph' %in% class(graph)) {
edges = graph %>% as_edgelist %>% as_tibble(.name_repair = ~c('from', 'to'))
} else {
edges = graph %>% as_tibble(.name_repair = make.unique)
if(nrow(edges) == 0) browser()
graph %<>% select(1:2) %>% as.matrix %>% graph_from_edgelist()
}
names(edges)[1:2] = c('from', 'to')
if('collapse_edges' %in% input$plotopt) {
withProgress(message = 'collapsing edges...', {
edges %<>% admixtools:::collapse_edges(10^input$collapse_threshold)
})
}
admixnodes = unique(edges[[2]][edges[[2]] %in% names(which(table(edges[[2]]) > 1))])
graph = edges %>% select(1:2) %>% as.matrix %>% graph_from_edgelist()
pos = data.frame(names(V(graph)), igraph::layout_as_tree(graph), stringsAsFactors = F) %>%
set_colnames(c('node', 'x', 'y'))
eg = edges %>% left_join(pos, by=c('from'='node')) %>%
left_join(pos %>% transmute(to=node, xend=x, yend=y), by='to') %>%
mutate(type = ifelse(to %in% admixnodes, 'admix', 'normal'))
if('reorder_edges' %in% input$plotopt) {
withProgress(message = 'reordering edges...', {
eg = admixtools:::fix_layout(eg %>% rename(name = from), graph) %>% rename(from = name)
})
}
if('shift_down' %in% input$plotopt) {
withProgress(message = 'reordering edges...', {
eg = admixtools:::fix_shiftdown(eg %>% rename(name = from), graph) %>% rename(from = name)
})
}
# if('simplify_graph' %in% input$plotopt) {
# withProgress(message = 'simplifying...', {
# isolate({
# global$graph = admixtools:::edges_to_igraph(eg) %>% admixtools:::simplify_graph()
# global$edges = edges = global$graph %>%
# igraph::as_edgelist() %>% as_tibble(.name_repair = ~c('from', 'to'))
# })
# eg = edges %>%
# left_join(pos, by=c('from'='node')) %>%
# left_join(pos %>% transmute(to=node, xend=x, yend=y), by='to') %>%
# mutate(type = ifelse(to %in% admixnodes, 'admix', 'normal'))
# })
# }
if('weight' %in% names(edges)) {
fa = function(x) paste0(round(x*100), '%')
fe = function(x) round(x*edgemul)
edges$label = ifelse(edges$type == 'admix', fa(edges$weight), fe(edges$weight))
if(!is.null(global$qpgraph_ranges)) {
edges %<>% left_join(global$qpgraph_ranges, by = c('from', 'to')) %>%
mutate(label = ifelse(type == 'admix',
paste(fa(mid), paste0('[', fa(lo), '-', fa(hi), ']'), sep = '\n'),
paste(fe(mid), paste0('[', fe(lo), '-', fe(hi), ']'), sep = '\n')))
}
}
if(!'label' %in% names(edges)) edges %<>% mutate(label='')
e2 = edges %>% transmute(from, to, label) %>% left_join(count(., to), by = c('from'='to'))
eg %<>% left_join(e2, by=c('from', 'to')) %>% mutate(indegree = replace_na(n, 0))
nodes = eg %>% filter(to %in% get_leafnames(graph)) %>% rename(x=xend, y=yend, xend=x, yend=y) %>%
transmute(name = to, x, y, xend, yend, to=NA, rownum = 1:n())
namedList(nodes, eg)
})
show_popadjust = reactive(length(global$poplist) > 0)
show_extract = reactive(isFALSE(global$iscountdata) && isFALSE(global$isf2data))
show_indselect = reactive(show_extract() || isTRUE(global$iscountdata))
observeEvent(input$minus1_cell_clicked, {
row = input$minus1_rows_selected
req(row)
sel = get_minus1() %>% slice(row)
global$selgraph = sel$graph[[1]]
global$edges = sel$edges[[1]]
global$score = sel$score[[1]]
global$worst_residual = sel$worst_residual[[1]]
add_graph_to_history(global)
global$tempgraph = TRUE
shinyjs::enable('qpgraph_update')
shinyjs::enable('qpgraph_revert')
print('newgraph set...')
print(numadmix(global$graph))
})
observeEvent(input$plus1_cell_clicked, {
row = input$plus1_rows_selected
req(row)
sel = get_plus1() %>% slice(row)
global$selgraph = sel$graph[[1]]
global$edges = sel$edges[[1]]
global$score = sel$score[[1]]
global$tempgraph = TRUE
shinyjs::enable('qpgraph_update')
shinyjs::enable('qpgraph_revert')
})
observeEvent(input$flipadmix_cell_clicked, {
row = input$flipadmix_rows_selected
req(row)
sel = get_flipadmix() %>% slice(row)
global$selgraph = sel$graph[[1]]
global$edges = sel$edges[[1]]
global$score = sel$score[[1]]
global$worst_residual = sel$worst_residual[[1]]
add_graph_to_history(global)
global$tempgraph = TRUE
shinyjs::enable('qpgraph_update')
shinyjs::enable('qpgraph_revert')
})
observeEvent(input$history_cell_clicked, {
row = input$history_rows_selected
req(row)
sel = get_history() %>% slice(row)
global$selgraph = sel$graph[[1]]
global$edges = sel$edges[[1]]
global$score = sel$score[[1]]
global$worst_residual = sel$worst_residual[[1]]
#add_graph_to_history(global)
global$tempgraph = TRUE
shinyjs::enable('qpgraph_update')
shinyjs::enable('qpgraph_revert')
})
observeEvent(input$addleaf_cell_clicked, {
row = input$addleaf_rows_selected
req(row)
sel = get_addleaf() %>% slice(row)
global$selgraph = sel$graph[[1]]
global$edges = sel$edges[[1]]
global$score = sel$score[[1]]
global$worst_residual = sel$worst_residual[[1]]
add_graph_to_history(global)
global$tempgraph = TRUE
shinyjs::enable('qpgraph_update')
shinyjs::enable('qpgraph_revert')
})
observeEvent(input$initresample_cell_clicked, {
row = input$initresample_rows_selected
req(row)
sel = get_initresample0() %>% slice(row)
global$edges = sel$edges[[1]]
global$score = sel$score[[1]]
global$worst_residual = sel$worst_residual[[1]]
add_graph_to_history(global)
global$tempgraph = TRUE
shinyjs::enable('qpgraph_update')
shinyjs::enable('qpgraph_revert')
})
observeEvent(input$indresample_cell_clicked, {
row = input$indresample_rows_selected
req(row)
sel = get_indresample0() %>% slice(row)
global$edges = sel$edges[[1]]
global$score = sel$score[[1]]
global$worst_residual = sel$worst_residual[[1]]
add_graph_to_history(global)
global$tempgraph = TRUE
shinyjs::enable('qpgraph_update')
shinyjs::enable('qpgraph_revert')
})
observeEvent(input$snpresample_cell_clicked, {
row = input$snpresample_rows_selected
req(row)
sel = get_snpresample0() %>% slice(row)
global$edges = sel$edges[[1]]
global$score = sel$score[[1]]
global$worst_residual = sel$worst_residual[[1]]
add_graph_to_history(global)
global$tempgraph = TRUE
shinyjs::enable('qpgraph_update')
shinyjs::enable('qpgraph_revert')
})
observeEvent(input$qpadm_rot_cell_clicked, {
row = input$qpadm_rot_rows_selected
req(row)
sel = get_qpadm_rotate() %>% slice(row)
global$qpadmpops_left_remain = sel$left[[1]]
global$qpadmpops_right_add = sel$right[[1]]
})
observeEvent(input$multiprocess, {
print('multiprocess')
if(input$multiprocess) future::plan('multiprocess')
else future::plan('sequential')
})
make_reactive_graphtabs = function(graphfun, inp) {
reactive({
print('graph neighbors')
f2blocks = get_f2blocks()
print('graph neighbors2')
#g = global$graph
g = get_graph()
withProgress(message = 'Finding graphs...', {
out = g %>% graphfun
})
print(paste('found', nrow(out), 'graphs'))
withProgress(message = paste('Evaluating', nrow(out), 'graphs...'), {
qpgfun = qpgraphfun()
out %<>% mutate(res = furrr::future_map(graph, ~qpgfun(f2blocks, .))) %>%
unnest_wider(res) %>%
select(starts_with('from'), starts_with('to'), starts_with('score'), starts_with('graph'), starts_with('edges'))
})
if(nrow(out) == 0) {
shinyalert('No graphs could be evaluated!')
return()
}
out %>% arrange(score)
})
}
get_minus1 = make_reactive_graphtabs(graph_minusone, input$qpgraph_similar_minus1)
get_plus1 = make_reactive_graphtabs(graph_plusone, input$qpgraph_similar_plus1)
get_flipadmix = make_reactive_graphtabs(admixtools::graph_flipadmix, input$qpgraph_similar_flipadmix)
get_addleaf = eventReactive(input$qpgraph_add_run, {
f2blocks = get_f2blocks()
g = global$graph
withProgress(message = 'Finding graphs...', {
out = g %>% graph_addleaf(input$addleaf)
})
withProgress(message = paste('Evaluating', nrow(out), 'graphs...'), {
qpgfun = qpgraphfun()
out %<>% mutate(res = furrr::future_map(graph, ~qpgfun(f2blocks, .))) %>%
unnest_wider(res) %>%
select(starts_with('from'), starts_with('to'), starts_with('score'), starts_with('graph'), starts_with('edges'))
})
if(nrow(out) > 0) out %<>% arrange(score)
out
})
qpgraph_scorerange = function(out, low = 0, high = 1) {
out %>%
select(score) %>%
filter(!is.na(score)) %>%
arrange(score) %>%
slice(ceiling((low+1e-10)*n()), round(n()/2), ceiling(n()*high)) %>%
add_column(type = c('lo', 'mid', 'hi'), .before = 1) %>%
deframe
}
qpgraph_ranges = function(out, low = 0, high = 1) {
out %>%
select(from, to, weight) %>%
filter(!is.na(weight)) %>%
group_by(from, to) %>%
arrange(weight) %>%
slice(ceiling((low+1e-10)*n()), round(n()/2), ceiling(n()*high)) %>%
mutate(type = c('lo', 'mid', 'hi')) %>%
ungroup %>%
select(from, to, type, weight) %>%
pivot_wider(names_from = type, values_from = weight)
}
add_graph_to_history = function(global) {
print('add_graph_to_history')
isolate({
new = list(graph = global$graph,
edges = global$edges,
nadmix = numadmix(global$graph),
score = global$score,
worst_residual = if(is.null(global$worst_residual)) NA else global$worst_residual)
global$history %<>% c(list(new))
})
}
get_initresample0 = reactive({
print('get initresample')
input$init_resample
f2blocks = get_f2blocks()
g = global$graph
numstart = as.numeric(isolate(input$numstart2))
withProgress(message = paste('Evaluating graph ', numstart, ' times...'), {
qpgfun = qpgraphfun()
out = tibble(id = seq_len(numstart)) %>%
mutate(res = furrr::future_map(id, ~qpgfun(f2blocks, g, numstart = 1))) %>%
unnest_wider(res) %>%
select(id, starts_with('score'), starts_with('opt'), starts_with('graph'), starts_with('edges'))
})
})
get_initresample = reactive({
out = get_initresample0()
if(nrow(out) > 0 && 'edges' %in% names(out)) {
global$qpgraph_scorerange = qpgraph_scorerange(out, low = 0.025, high = 0.975)
out %<>% unnest(edges)
global$qpgraph_ranges = qpgraph_ranges(out, low = 0.025, high = 0.975)
out %<>% mutate(edge = paste(from, to, sep = ' -> ')) %>%
select(id, opt, score, edge, weight) %>%
pivot_wider(names_from = edge, values_from = weight) %>%
unnest(opt)
}
print('initresample done')
out
})
get_indresample0 = reactive({
print('get indresample')
input$ind_resample
poplist = global$poplist
countdir = global$countdir
inds = unlist(poplist)
numsingle = sum(sapply(poplist, length) == 1)
pops = rep(names(poplist), sapply(poplist, length))
f2blocks = get_f2blocks()
g = global$graph
qpgfun = qpgraphfun()
fun = admixtools:::make_resample_inds_fun(~qpgfun(..., y = g))
withProgress(message = paste('Evaluating graph ', length(inds)-numsingle, ' times...'), {
out = fun(countdir, inds = inds, pops = pops, numstart = 1)
})
out
})
get_indresample = reactive({
out = get_indresample0()
if(nrow(out) > 0) {
print('head(out)')
print(head(out))
print(head(out$error))
global$qpgraph_scorerange = qpgraph_scorerange(out, low = 0.025, high = 0.975)
out %<>% unnest(edges)
global$qpgraph_ranges = qpgraph_ranges(out, low = 0.025, high = 0.975)
out %<>%
mutate(edge = paste(from, to, sep = ' -> ')) %>%
select(ind, pop, edge, score, weight) %>%
pivot_wider(names_from = edge, values_from = weight)
}
print('indresample done')
out
})
get_snpresample0 = reactive({
f2blocks = get_f2blocks()
input$snp_resample
g = global$graph
bootnum = ifelse(input$boot, input$bootnum, FALSE)
withProgress(message = paste('Evaluating graph ', ifelse(bootnum, bootnum, dim(f2blocks)[[3]]),' times...'), {
out = admixtools:::make_resample_snps_fun(~qpgraphfun()(..., y = g))(f2blocks, boot = bootnum, numstart = 1)
})
})
get_snpresample = reactive({
out = get_snpresample0()
if(nrow(out) > 0 && 'edges' %in% names(out)) {
print('head(out)')
print(head(out))
global$qpgraph_scorerange = qpgraph_scorerange(out, low = 0.025, high = 0.975)
out %<>% unnest(edges)
global$qpgraph_ranges = qpgraph_ranges(out, low = 0.025, high = 0.975)
out %<>%
mutate(edge = paste(from, to, sep = ' -> ')) %>%
select(id, edge, score, weight) %>%
pivot_wider(names_from = edge, values_from = weight)
}
print('snpres done')
out
})
get_qpadm_rightpanel = reactive({
print('get qpadm_rightpanel')
global$qpadm_rightpanel
})
output$dashboardbody = renderUI({
print('dashboardbody')
get_bod()
})
output$popselectors = renderUI({
poplist = global$poplist
indnames = unique(global$allinds)
print('popselectors')
print(indnames)
sellist = imap(names(poplist), ~{
#buttlab = ifelse(is.null(global$poplistbak[[.x]]), 'group', 'ungroup')
isgrouped = length(poplist[[.x]]) == 1 && admixtools:::is_group(global$countdir, poplist[[.x]])
buttlab = ifelse(isgrouped, 'ungroup', 'group')
fillCol(tagList(textInput(paste0('pop', .y), NULL, .x, placeholder = 'pop name'),
shiny::actionButton(paste0('group', .y), buttlab),
selectizeInput(.x, NULL, indnames,
poplist[[.x]],
multiple = TRUE, width = 'auto'),
tags$style(type="text/css", paste0('.selectize-control{white-space:pre-wrap}'))), height = '1000px')})
print('popselectors2')
do.call(splitLayout, sellist)
})
output$econstraints = renderUI({
global$useconstraints = FALSE
print('econstraints')
edgecols = isolate(get_edgecols())
edges = isolate(get_fit()$edges) %>%
mutate(nam = paste(from, to, sep = ' -> '), nam2 = str_remove_all(nam, '\\.|>| ')) %>%
left_join(edgecols, by = c('nam2'='edge')) %>%
filter(type == 'edge') %>% select(nam, nam2, weight, col)
maxlen = ceiling(max(edges$weight))*edgemul
print('econstraints 2')
div(div(id = 'econstraints_update_div', actionButton('econstraints_update', 'Update')),
sliderInput('any_econstraints', 'All edges', 0, maxlen, c(0, maxlen), dragRange = FALSE),
edges %>%
pmap(~{div(tags$style(HTML(paste0('[for=',..2,']+span>.irs-bar, [for=',..2,']+span>.irs>.irs-from, [for=',..2,']+span>.irs>.irs-to {background: ', ..4,'}'))),
sliderInput(..2, ..1, min = 0, max = maxlen, value = c(0, maxlen),
dragRange = FALSE, width = '100%'))}))
})
output$aconstraints = renderUI({
global$useconstraints = FALSE
print('aconstraints')
input$add_edge
input$delete_edge
edgecols = isolate(get_edgecols())
edges = isolate(get_fit()$edges) %>%
mutate(nam = paste(from, to, sep = ' -> '), nam2 = str_remove_all(nam, '\\.|>| ')) %>%
left_join(edgecols, by = c('nam2'='edge')) %>%
filter(type == 'admix') %>% select(nam, nam2, weight, col)
div(div(id = 'aconstraints_update_div', actionButton('aconstraints_update', 'Update')),
edges %>%
pmap(~{div(tags$style(HTML(paste0('[for=',..2,']+span>.irs-bar, [for=',..2,']+span>.irs>.irs-from, [for=',..2,']+span>.irs>.irs-to {background: ', ..4,'}'))),
sliderInput(..2, ..1, min = 0, max = 1, value = c(0, 1),
dragRange = FALSE, width = '100%'))})
)
})
output$nadmix = renderUI({
nad = 3
if(!is.null(global$graph)) nad = numadmix(global$graph)
div(sliderInput('nadmix', '# admix', 0, 10, nad),
selectizeInput('outpop', 'Outgroup', c('< undefined >', names(global$poplist)), selected = names(global$poplist)[1]))
})
output$qpgraph_add = renderUI({
print('selectInput addleaf')
choices = setdiff(names(global$poplist), get_leafnames(global$graph))
selectInput('addleaf', '', choices = choices)
})
output$showseledge = renderText({
input$clear_edge
print('showseledge1')
se = global$seledge
print('se')
print(global$seledge)
print(se)
if(!is.null(se) && (is.na(se) || se == 'NA')) {
print('is na')
se = NULL
}
print('se2')
print(se)
se
})
observeEvent(input$any_econstraints, {
print('update triggered by any_econstraints')
req(global$graph)
print('any_econstraints 2')
edges = global$graph %>% igraph::as_edgelist() %>% as_tibble(.name_repair = ~c('from', 'to')) %>%
group_by(to) %>% mutate(type = ifelse(n() > 1, 'admix', 'normal')) %>% ungroup %>%
mutate(eid = paste(from, to, sep = '->') %>% str_remove_all('\\.|>| ')) %$% eid
print('any_econstraints 3')
if(is.null(isolate(input[[edges[1]]]))) return()
edges %>% walk(~{updateSliderInput(session, ., value = c(input$any_econstraints[1],
input$any_econstraints[2]))})
print('any_econstraints 4')
})
observeEvent(input$delete_edge, {
print('del 1')
eg = get_seledge()
if(is.null(eg)) return()
gnew = global$graph
leaves = get_leafnames(gnew)
print('del 2')
print(eg)
eg %<>% str_split('\n') %>% `[[`(1)
if(str_detect(eg[1], ' -> ')) {
eg %<>% str_split(' -> ')
print(eg)
for(i in 1:length(eg)) {
#browser()
gnew = delete_admix(gnew, eg[[i]][1], eg[[i]][2])
# g <<- g
# gnew <<- gnew
if(!admixtools:::is_valid(gnew)) {
shinyalert('Error', 'Could not delete edge!')
return()
}
}
} else {
for(i in 1:length(eg)) {
gnew = delete_leaf(gnew, eg[i])
if(!admixtools:::is_valid(gnew)) {
browser()
shinyalert('Error', 'Could not delete node!')
return()
}
}
}
newleaves = setdiff(get_leafnames(gnew), leaves)
if(length(newleaves) > 0) {
# newleaves <<- newleaves
# oldgraph <<- global$graph
# newgraph <<- g
shinyalert('New leaves!', newleaves)
}
global$graph = gnew
global$seledge = NULL
global$edges = gnew %>% as_edgelist %>% as_tibble(.name_repair = ~c('from', 'to'))
global$score = NA
global$worst_residual = NA
shinyjs::disable('clear_edge')
shinyjs::disable('add_edge')
shinyjs::disable('delete_edge')
print('delete_edge done')
})
observeEvent(input$clear_edge, {
print('clear')
global$seledge = NULL
shinyjs::disable('clear_edge')
shinyjs::disable('add_edge')
shinyjs::disable('delete_edge')
})
observeEvent(input$add_edge, {
print('add')
eg = get_seledge() %>% str_split('\n') %>% `[[`(1) %>% str_split(' -> ')
g = global$graph
leaves = get_leafnames(g)
if(length(eg) == 1) {
print('input$addleaf')
print(input$addleaf)
if(is.null(input$addleaf)) {
shinyalert('Error', 'Can only add edge which connects exactly two other edges!')
return()
}
from = eg[[1]][1]
to = eg[[1]][2]
gnew = admixtools:::insert_leaf(g, input$addleaf, from, to)
if(!admixtools:::is_valid(gnew)) {
shinyalert('Error', 'Could not add population!')
return()
}
global$graph = gnew
global$seledge = NULL
global$edges = gnew %>% as_edgelist %>% as_tibble(.name_repair = ~c('from', 'to'))
global$score = NA
global$worst_residual = NA
shinyjs::disable('clear_edge')
shinyjs::disable('add_edge')
shinyjs::disable('delete_edge')
} else if(length(eg) == 2) {
#from = eg[[1]][2]
#to = eg[[2]][2]
#print(paste(from, to))
alert = function(x) shinyalert('Could not insert edge!', as.character(x))
tryCatch({
#gnew = admixtools:::insert_admix_old(g, from, to, allow_below_admix = TRUE)
gnew = admixtools:::insert_admix(g, eg[[1]][1], eg[[1]][2], eg[[2]][1], eg[[2]][2])
}, warning = alert, error = alert)
if(!exists('gnew') || is.null(gnew)) return()
if(!admixtools:::is_valid(gnew)) {
shinyalert('Error', 'Could not insert edge!')
return()
}
global$graph = gnew
global$seledge = NULL
global$edges = gnew %>% as_edgelist %>% as_tibble(.name_repair = ~c('from', 'to'))
global$score = NA
global$worst_residual = NA
shinyjs::disable('clear_edge')
shinyjs::disable('add_edge')
shinyjs::disable('delete_edge')
newleaves = setdiff(get_leafnames(gnew), leaves)
if(length(newleaves) > 0) shinyalert('New leaves!', newleaves)
} else {
shinyalert('Error', 'Too many items selected!')
return()
}
})
observeEvent(input$run_qpadm, {
print('observe run_qpadm')
global$qpadmout = get_qpadmout()
global$qpadm_rightpanel = plotOutput('qpadmcomparison')
})
observeEvent(input$run_qpgraph, {
print('observe run_qpgraph')
global$qpgraphout = get_qpgraphout()
#global$qpg_right = plotlyOutput('graphcomparison')
})
map(1:4,
~observeEvent(input[[paste0('f4pops', .)]], {
if(input$f4_allowdup) return()
popsinuse = unique(c(input$f4pops1, input$f4pops2, input$f4pops3, input$f4pops4))
map(setdiff(1:4, .), ~{
fp = paste0('f4pops', .)
sel = input[[fp]]
updateSelectizeInput(session, fp, selected = sel,
choices = union(sel, setdiff(names(global$poplist), popsinuse)))})
}))
observeEvent(input$f4_allowdup, {
map(1:4, ~{
fp = paste0('f4pops', .)
sel = input[[fp]]
updateSelectizeInput(session, fp, selected = sel,
choices = names(global$poplist))})
})
map(1:3, ~observeEvent(input[[paste0('qpadmpops', .)]], {
popsinuse = unique(c(input$qpadmpops1, input$qpadmpops2, input$qpadmpops3))
map(setdiff(1:3, .), ~{
fp = paste0('qpadmpops', .)
sel = input[[fp]]
global[[fp]] = sel
choices = union(sel, setdiff(names(global$poplist), popsinuse))
updateSelectizeInput(session, fp, selected = sel, choices = choices)})
}))
observeEvent(input$qpadm_randomize, {
nam = sample(names(global$poplist))
nl = floor(length(nam)/2)
choices = list(nam[1], nam[2:nl], nam[(nl+1):length(nam)])
map(1:3, ~{
global[[paste0('qpadmpops', .)]] = choices[[.]]
updateSelectizeInput(session, paste0('qpadmpops', .),
selected = choices[[.]], choices = nam)})
})
observeEvent(input$qpadm_rotate, {
#global$show_qpadm_rotate = input$qpadm_rotate
global$qpadmpops1 = input$qpadmpops1
global$qpadmpops2 = input$qpadmpops2
global$qpadmpops3 = input$qpadmpops3
})
observeEvent(input$qpadm_results, {
global$qpadmpops1 = input$qpadmpops1
global$qpadmpops2 = input$qpadmpops2
global$qpadmpops3 = input$qpadmpops3
#global$show_qpadm_results = input$qpadm_results
})
get_qpadm_rotate = reactive({
f2blocks = get_f2blocks()
p1 = input$qpadmpops1
p2 = input$qpadmpops2
p3 = input$qpadmpops3
req(f2blocks, p1, p2, p3)
withProgress(message = paste0('Generating qpadm models...'), {
lr = admixtools:::all_lr2(p1, length(p2))
})
withProgress(message = paste0('Evaluating ', length(lr[[1]]), ' qpadm models...'), {
out = admixtools:::qpadm_eval_rotate(f2blocks, p3, lr, p2, verbose = F) %>%
select(chisq, p, feasible, left, right)
})
out
})
get_qpadmout = reactive({
print('get_qpadmout')
get_qpadm()
volumes %<>% eval
pref = parseFilePaths(volumes, input$qpadm_genofile)$datapath
pref %<>% str_replace('\\.geno$', '') %>% normalizePath(mustWork = FALSE)
bin = parseFilePaths(volumes, input$qpadmbin)$datapath %>% normalizePath(mustWork = FALSE)
#if(length(bin) == 0) bin = '~/Downloads/admixtoools_hub_reldec19/bin/qpAdm'
#if(length(pref) == 0) pref = '~/Documents/v42.1_afafam'
target = input$qpadmpops1
left = input$qpadmpops2
right = input$qpadmpops3
withProgress(message = 'running original qpAdm...', {
dir = tempdir()
out = tryCatch(qpadm_wrapper(left, right, target, bin = bin, pref = pref, outdir = dir),
error = function(e) shinyalert('Error!', as.character(e)))
})
out
})
get_qpgraphout = reactive({
print('get_qpgraphout')
get_fit()
volumes %<>% eval
pref = parseFilePaths(eval(volumes), input$qpgraph_genofile)$datapath
pref %<>% str_replace('\\.geno$', '') %>% normalizePath(mustWork = FALSE)
bin = parseFilePaths(volumes, input$qpgraphbin)$datapath %>% normalizePath(mustWork = FALSE)
print('bin')
print(bin)
# if(length(bin) == 0) bin = '~/Downloads/admixtoools_hub_reldec19/bin/qpGraph'
# if(length(pref) == 0) pref = '~/Documents/v42.1_afafam'
g = global$graph
withProgress(message = 'running original qpGraph...', {
dir = tempdir() %>% normalizePath()
print('pref')
print(pref)
print('dir')
print(dir)
out = qpgraph_wrapper(g, bin, pref, outdir = dir, zthresh = 0)
out$f2 = NULL
print('parse fstats')
fstatsfile = paste0(dir, '/fstats.out')
if(file.exists(fstatsfile)) global$precomp = admixtools:::parse_fstats(fstatsfile)
})
out
})
observeEvent(input$econstraints_update, {global$useconstraints = TRUE})
observeEvent(input$aconstraints_update, {global$useconstraints = TRUE})
dto = function(x) div(DT::DTOutput(x), style = dtstyle)
#dtof4 = reactive({ if(input$return_fstats) dto('f4') else renderText("Activate 'return_fstats' under 'Fit'!") })
observeEvent(input$qpadm_fit, {global$qpadm_rightpanel = qpadm_rightpanel_fit()})
# observeEvent(input$qpgraph_fit, {
# print('fit!')
# global$qpgraph_ranges = NULL
# global$qpgraph_scorerange = NULL
# global$qpg_right = qpg_right_fit()
# })
qpadm_leftpanel = reactive({
print('reactive qpadm_leftpanel')
all = names(global$poplist)
choices1 = setdiff(all, c(global$qpadmpops2, global$qpadmpops3))
choices2 = setdiff(all, c(global$qpadmpops1, global$qpadmpops3))
choices3 = setdiff(all, c(global$qpadmpops1, global$qpadmpops2))
div(
selectizeInput('qpadmpops1', 'Left', choices = choices1, multiple = TRUE, selected = global$qpadmpops1),
selectizeInput('qpadmpops2', 'Right', choices = choices2, multiple = TRUE, selected = global$qpadmpops2),
selectizeInput('qpadmpops3', 'Target', choices = choices3, multiple = FALSE, selected = global$qpadmpops3),
splitLayout(actionButton('qpadm_randomize', 'Randomize')))
})
qpadm_rightpanel_fit = reactive({
tabsetPanel(id = 'qpadm_tabset',
tabPanel('Weights', dto('qpadm_weights')),
tabPanel('f4', dto('qpadm_f4')),
tabPanel('Rank drop', dto('qpadm_rankdrop')),
tabPanel('Pop drop', dto('qpadm_popdrop')))
})
qpadm_rightpanel_rotate = reactive(div(dto('qpadm_rot'), style='padding-right:50px;'))
qpadm_rightpanel = reactive({
#res = global$show_qpadm_results
#rot = global$show_qpadm_rotate
res = if(is.null(input$qpadm_results)) TRUE else input$qpadm_results
rot = if(is.null(input$qpadm_rotate)) FALSE else input$qpadm_rotate
if(rot && res) splitLayout(cellWidths = c('50%', '50%'), qpadm_rightpanel_rotate(), qpadm_rightpanel_fit())
else if(rot) qpadm_rightpanel_rotate()
else if(res) qpadm_rightpanel_fit()
else div()
})
get_loaddata = reactive({
print('load_data')
input$textdirinput
bh = '140px'
div(
fluidRow(
box(width=4, height=bh, background = cols[1], h4('Select data directory'),
div(splitLayout(textInput('textdirinput', NULL, placeholder = 'Data directory'),
actionButton('textdir', 'Create')))),
conditionalPanel('output.show_extract == "1"',
fluidRow(box(width=4, height=bh, background = cols[2], h4('Extract data'),
splitLayout(
div(
textInput('textgenofile1', NULL, placeholder = 'genotype file')
),
actionButton('extract_data', 'Extract data'))))),
conditionalPanel('output.show_indselect == "1"', (
box(width=4, height=bh, background = cols[3], h4('Select .ind or .fam file'),
div(fileInput('popfile', NULL, placeholder = '', buttonLabel = 'Ind file'), id = 'popfilediv')))),
fluidRow(column(12, conditionalPanel('input.extract_data > 0', verbatimTextOutput('console')))),
#box(width=6, textOutput('console')),
div(style = 'visibility: hidden', verbatimTextOutput('show_popadjust')),
div(style = 'visibility: hidden', verbatimTextOutput('show_extract')),
div(style = 'visibility: hidden', verbatimTextOutput('show_indselect')),
conditionalPanel('output.show_popadjust == "1"', fluidRow(column(12,
box(width=12, background = cols[5],
fluidRow(column(2, h4('Adjust populations', id = 'adjpopdiv')),
column(1, div(actionButton('removepop', '–'),
actionButton('addpop', '+')))),
fluidRow(column(11, uiOutput('popselectors')))))))))
})
observeEvent(input$extract_data, {
if(is.null(input$textgenofile1) || !file.exists(input$textgenofile1)) {
shinyalert('Error!', 'Please provide path and name for a .geno or .bed file!')
return()
}
pref = input$textgenofile1 %>% str_replace('\\.geno$|\\.bed$', '')
oldnam = names(global$poplist)
nam = map(paste0('pop', seq_len(length(global$poplist))), ~input[[.]])
if(is.null(nam[[1]])) return()
global$poplist = map(oldnam, ~input[[.]]) %>% purrr::set_names(nam)
poplist = global$poplist
inds = unlist(poplist)
pops = rep(names(poplist), sapply(poplist, length))
transitions = input$trans_extract %in% c('both', 'only transitions')
transversions = input$trans_extract %in% c('both', 'only transversions')
if(input$fix_populations) extractfun = function(...) extract_f2(pops = pops, ...)
else extractfun = extract_counts
print('inds')
print(inds)
print('global$countdir')
print(global$countdir)
withCallingHandlers({
shinyjs::html('console', '')
print(list(pref, global$countdir, inds = inds,
maxmiss = input$max_miss, minmaf = input$min_maf, maxmaf = input$max_maf,
transitions = transitions, transversions = transversions,
keepsnps = input$keepsnps, maxmem = input$maxmem*1e3))
extractfun(pref, global$countdir, inds = inds,
maxmiss = input$max_miss, minmaf = input$min_maf, maxmaf = input$max_maf,
transitions = transitions, transversions = transversions,
keepsnps = input$keepsnps, maxmem = input$maxmem*1e3)
},
message = function(m) {
shinyjs::html(id = 'console', html = m$message, add = TRUE)
})
global$iscountdata = 'indivs' %in% list.dirs(global$countdir,F,F)
global$isf2data = !'indivs' %in% list.dirs(global$countdir,F,F) &
'block_lengths_f2.rds' %in% list.files(global$countdir)
if(global$iscountdata) {
global$allinds = list.dirs(paste0(global$countdir, '/pairs'),F,F)
} else if(global$isf2data) {
global$allinds = list.dirs(global$countdir,F,F)
global$poplist = global$allinds %>% purrr::set_names(global$allinds)
}
})
observeEvent(event_data('plotly_click', source = 'src'), {
print('observe seledge')
ed = event_data('plotly_click', source = 'src')
isolate({
plt = plotly_graph()
newe = ed_to_name(ed, plt)
oldedge = global$seledge
if(!is.null(oldedge) && is.na(oldedge)) {
print('seledge is na. why?')
oldedge = NULL
}
global$seledge = newe
shinyjs::enable('clear_edge')
shinyjs::enable('add_edge')
shinyjs::enable('delete_edge')
})
if(!is.null(newe) && !is.null(oldedge) && newe != oldedge) global$seledge = paste(c(oldedge, newe), collapse = '\n')
})
get_seledge = reactive({
print('get_seledge')
print(global$seledge)
global$seledge
})
ed_to_name = function(ed, plt) {
dat = plt$x$data[[ed$curveNumber + 1]]
type = ifelse(length(dat) == 14, 'text', ifelse('hovertext' %in% names(dat), 'edge', 'node'))
pt = ed$pointNumber + 1
if (type == 'edge') {
name = dat$hovertext[pt]
} else if (type == 'node') {
col = ifelse(duplicated(dat$text)[pt], 5, 3)
name = dat$text[pt]
} else {
name = dat$text[pt]
}
print('ed_to_name22')
name
}
get_edgecols = reactive({
plt = plotly_graph()
edgecols = plt$x$data %>% map(~pluck(., 'line', 'color')) %>% discard(is.null)
txt = plt$x$data[seq_len(length(edgecols))] %>% map(~pluck(., 'text'))
tibble(edge = str_remove_all(unlist(txt), '\\.|>| '),
col = rep(as.vector(edgecols), map_dbl(txt, length))) %>%
unnest(col) %>% filter(!duplicated(.), !is.na(edge))
})
plotly_graph = reactive({
print('plotly_graph')
poplist = global$poplist
global$graph
pdat = get_pdat()
req(pdat, poplist)
nadmix = numadmix(graph_from_edgelist(as.matrix(global$edges)[,1:2]))
textsize = 2.5
print('plotly_graph 2')
eg = pdat$eg
nodelabs = imap(poplist, ~paste(.y, paste(.x, collapse = '\n'), sep = '\n')) %>%
unlist %>% enframe(name = 'name', value = 'text')
#nodelabs = poplist %>% group_by(pop) %>%
# summarize(text = paste(pop[1], paste(ind, collapse='\n'), sep = '\n'))
nodes = pdat$nodes %>% left_join(nodelabs, by = c('name'))
allnodes = eg %>% transmute(from, x, y, xend=0, yend=0, to=0) %>% filter(!duplicated(.))
segtext = "ifelse(indegree == 1, to, from)"
segtext = "paste(from, to, sep = ' -> ')"
print('plotly_graph 3')
sr = global$qpgraph_scorerange
scoretext = paste0('score: ', ifelse(is.null(sr), round(global$score, 2),
paste0(round(sr['mid'], 2), '\n[', round(sr['lo'], 0), ' - ', round(sr['hi'], 0), ']')),
ifelse(is.na(global$worst_residual), '', paste0('\nWR: ', round(global$worst_residual, 2))),
'\nadmix: ', nadmix,
ifelse(global$tempgraph, '\n\nPREVIEW!', ''))
gg = eg %>% mutate(rownum = 1:n()) %>%
ggplot(aes(x=x, xend=xend, y=y, yend=yend, from=from, to=to)) +
geom_segment(aes_string(linetype = 'type', col = 'as.factor(y)', text = segtext),
arrow=arrow(type = 'closed', angle = 10, length=unit(0.15, 'inches'))) +
geom_point(data = allnodes, aes(x, y, text = from), col = 'black', alpha = 0) +
annotate('text', x = min(nodes$x), y = max(nodes$yend), hjust = 0, vjust = 1,
label = scoretext) +
# annotate('text', x = min(nodes$x), y = max(nodes$yend), hjust = 0,
# label = scoretext) +
theme(panel.background = element_blank(),
axis.line = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank(),
legend.position = 'none') +
xlab('') + ylab('') +
scale_linetype_manual(values = c(admix=3, normal=1)) +
ggtitle('') +
scale_x_continuous(expand = c(0.1, 0.1))
if('highlight_unidentifiable' %in% input$plotopt) {
unid = unidentifiable_edges(global$graph)
unid2 = eg %>% right_join(unid %>% select(-type), by = c('from', 'to'))
gg = gg + geom_segment(aes_string(linetype = 'type'), col = 'red', size = 1, data = unid2)
}
gg = gg + geom_text(aes(x = (x+xend)/2, y = (y+yend)/2, label = label, text = paste(from, to, sep = ' -> ')),
size = textsize) +
geom_text(data = nodes, aes_string(label = 'name', col = 'as.factor(yend)',
from = NA, text = 'text'), size = textsize)
print('plotly_graph 4')
plt = plotly::ggplotly(gg, source = 'src', tooltip=c('text'))
print('plotly_graph 5')
plt
})
output$graphplot <- renderPlotly({
print('graphplot')
suppressWarnings({withProgress(message = 'Processing...', {plt = plotly_graph()})})
#event_register(plt, 'plotly_hover')
print('graphplot2')
plt
})
get_f2 = reactive({
print('get f2')
f2blocks = if(input$fst == 'fst') get_fstblocks() else get_f2blocks()
#req(f2blocks, p1, p2, p3, p4)
f2(f2blocks)
})
get_f3 = reactive({
print('get f3')
f2blocks = get_f2blocks()
f3(f2blocks)
})
get_f4 = reactive({
print('get f4')
input$f4_run
f2blocks = get_f2blocks()
p1 = input$f4pops1
p2 = input$f4pops2
p3 = input$f4pops3
p4 = input$f4pops4
req(f2blocks, p1, p2, p3, p4)
f4(f2blocks, p1, p2, p3, p4)
})
get_qpadm = reactive({
print('get qpadm')
input$qpadm_run
f2blocks = get_f2blocks()
p1 = input$qpadmpops1
p2 = input$qpadmpops2
p3 = input$qpadmpops3
req(f2blocks, p1, p2, p3)
#if(global$show_qpadm_rotate) {
if(input$qpadm_rotate) {
p1 = if(is.null(global$qpadmpops_left_remain)) p1 else global$qpadmpops_left_remain
p2 = union(p2, global$qpadmpops_right_add)
}
qpadm(f2blocks, p1, p2, p3)
})
plotly_f2 = reactive({
print('get f2')
f2blocks = if(input$fst == 'fst') get_fstblocks() else get_f2blocks()
pdat = f2blocks %>%
apply(1:2, mean, na.rm = T) %>%
as_tibble(rownames = 'pop1') %>%
pivot_longer(-pop1, 'pop2', values_to = 'est')
gg = pdat %>%
ggplot(aes(pop1, pop2, text = paste(pop1, pop2, round(est, 3), sep='\n'))) +
geom_tile(aes(fill = est)) +
xlab('') +
ylab('') +
theme(axis.text.x = element_text(angle = 90, hjust = 1), panel.background = element_blank()) +
scale_fill_distiller(palette = 'Purples')
plt = plotly::ggplotly(gg, source = 'src_f2', tooltip=c('text'))
plt
})
plotly_f2_indiv = reactive({
f2dat = get_f2dat_indiv()
poplist = global$poplist
ind = enframe(poplist, 'pop', 'ind') %>% unnest(ind)
dat = f2dat %>%
transmute(ind1 = pop1, ind2 = pop2, f2corr = f2, f2uncorr) %>%
left_join(ind %>% transmute(ind1 = ind, pop1 = pop), by = 'ind1') %>%
left_join(ind %>% transmute(ind2 = ind, pop2 = pop), by = 'ind2') %>%
arrange(pop1, pop2) %>%
mutate(i1 = as.numeric(factor(ind1, levels = unique(ind1[order(pop1)]))),
i2 = as.numeric(factor(ind2, levels = unique(ind2[order(pop2)]))),
f2 = ifelse(i1 > i2, f2corr, f2uncorr))
d2 = enframe(sapply(poplist, length)) %>%
arrange(name) %>%
mutate(end = cumsum(value), start = lag(end, default = 0), mn = start + value/2 + 0.5)
exp = 0
gg = dat %>%
ggplot(aes(i1, i2)) +
geom_tile(aes(fill = f2, text = paste(ind1, ind2, round(f2, 3), sep = '\n'))) +
geom_vline(aes(xintercept = end+0.5), data = d2 %>% head(-1) %>% rename(ind2 = name)) +
geom_hline(aes(yintercept = end+0.5), data = d2 %>% head(-1) %>% rename(ind2 = name)) +
scale_x_continuous(labels = d2$name, breaks = d2$mn, expand = c(exp, exp)) +
scale_y_continuous(labels = d2$name, breaks = d2$mn, expand = c(exp, exp)) +
xlab('') + ylab('') +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
panel.background = element_blank()) +
scale_fill_distiller(palette = 'Purples')
plt = plotly::ggplotly(gg, source = 'src_f2', tooltip = 'text')
plt
})
plotly_f4 = reactive({
print('get f4')
input$f4_run
#isolate({
f2blocks = get_f2blocks()
p1 = input$f4pops1[1]
p2 = input$f4pops2[1:2]
# f4pops2 needs at least two populations. how to check for this? otherwise, "Error: Only strings can be converted to symbols"
f4dat = get_f4() %>% filter(pop1 == p1, pop2 %in% p2)
req(f4dat)
#})
f4dat %<>%
mutate(pop2 = recode(pop2, !!sym(p2[1]):='x', !!sym(p2[2]):='y')) %>%
pivot_wider(names_from = pop2, values_from = c(est, se, z, p)) %>%
mutate(mul = input$f4_se)
xmin = min(f4dat$est_x - f4dat$se_x)
xmax = max(f4dat$est_x + f4dat$se_x)
ymin = min(f4dat$est_y - f4dat$se_y)
ymax = max(f4dat$est_y + f4dat$se_y)
gg = f4dat %>%
ggplot(aes(est_x, est_y, col = pop3, text = paste0('f4(', p1, ', X; ', pop3, ', ', pop4, ')'))) +
geom_smooth(method = 'lm', se = FALSE, formula=y~x-1) +
geom_hline(yintercept = 0) +
geom_vline(xintercept = 0) +
geom_point() +
theme(panel.background = element_blank()) +
xlab(paste0('f4(', p1, ', ', p2[1], '; Pop 3, Pop 4)' )) +
ylab(paste0('f4(', p1, ', ', p2[2], '; Pop 3, Pop 4)' )) +
scale_x_continuous(limits = c(xmin, xmax)) +
scale_y_continuous(limits = c(ymin, ymax))
if(input$f4_se > 0) {
gg = gg +
geom_errorbar(aes(ymin = est_y - se_y*mul, ymax = est_y + se_y*mul), width = 0) +
geom_errorbarh(aes(xmin = est_x - se_x*mul, xmax = est_x + se_x*mul), height = 0)
}
plt = plotly::ggplotly(gg, source = 'src_f4', tooltip=c('text'))
})
plotly_f42 = reactive({
print('get f4')
input$f4_run
#isolate({
f2blocks = get_f2blocks()
f4dat = get_f4() %>% filter(pop1 == input$f4pops1[1], pop3 == input$f4pops3[1], pop4 == input$f4pops4[1])
req(f4dat)
#})
# f4dat %<>%
# mutate(pop2 = recode(pop2, !!sym(p2[1]):='x', !!sym(p2[2]):='y')) %>%
# pivot_wider(names_from = pop2, values_from = c(est, se, z, p))
# xmin = min(f4dat$est_x - f4dat$se_x)
# xmax = max(f4dat$est_x + f4dat$se_x)
# ymin = min(f4dat$est_y - f4dat$se_y)
# ymax = max(f4dat$est_y + f4dat$se_y)
gg = f4dat %>%
ggplot(aes(pop2, est, col = pop2)) +
geom_hline(yintercept = 0) +
geom_point() +
theme(panel.background = element_blank()) +
xlab(paste0('' )) +
ylab(paste0('f4(', input$f4pops1[1], ', X; ', input$f4pops3[1], ', ', input$f4pops4[1],')' ))
if(input$f4_se > 0) {
mul = input$f4_se
gg = gg + geom_errorbar(aes(ymin = est - se*mul, ymax = est + se*mul), width = 0)
}
plt = plotly::ggplotly(gg, source = 'src_f42')
})
output$f2heatmap = renderPlotly({
print('f2plot')
plt = plotly_f2()
plt
})
output$f2heatmap_indiv = renderPlotly({
print('f2plot indiv')
plt = plotly_f2_indiv()
plt
})
output$f4plot = renderPlotly({
print('f4plot')
plt = plotly_f4()
plt
})
output$f4plot2 = renderPlotly({
print('f4plot')
plt = plotly_f42()
plt
})
output$qpadmcomparison = renderPlot({
print('qpadmcomparison')
withProgress(message = 'running qpadm in R...', {
qpadm_R = get_qpadm()
})
original = global$qpadmout
gg = plot_comparison(original, qpadm_R)
#plotly::ggplotly(gg, tooltip = 'label')
gg
})
output$graphcomparison = renderPlotly({
print('graphcomparison')
req(global$qpgraphout)
withProgress(message = 'running qpgraph in R...', {
qpgraph_R = get_fit()
})
#qpgraph_R$f3 = qpgraph_R$f3 %>% select(pop2, pop3, fit, est) # temporary; line should be removed
gg = plot_comparison(global$qpgraphout, qpgraph_R)
plotly::ggplotly(gg, tooltip = 'label')
})
output$downloadgraph = downloadHandler('graph.tsv', function(file) {
edges = get_fit()$edges
if(input$downloadgraph_format == 'admixtools') edges %>% admixtools:::fit_to_qpgraph_format() %>% write(file)
else if(input$downloadgraph_format == 'edgelist') edges %>% write_tsv(file)
else if(input$downloadgraph_format == 'dot') edges %>% admixtools:::write_dot(file)
else stop()
})
ex = c('Buttons', 'Scroller')
dtfun = function(...) DT::renderDataTable(..., extensions = ex, options = do, selection = 'single')
output$f2 = dtfun({format_table(get_fit()$f2)})
output$f3 = dtfun({format_table(get_fit()$f3)})
output$f4 = dtfun({if(is.null(get_fit()$f4)) as.data.frame("Activate 'return_fstats' under 'Fit'!") %>% slice(-1) else format_table(get_fit()$f4)})
output$opt = dtfun({format_table(get_fit()$opt)})
output$history = dtfun({format_table(get_history() %>% select(-graph, -edges))})
output$minus1 = dtfun({format_table(get_minus1())})
output$plus1 = dtfun(format_table(get_plus1()))
output$flipadmix = dtfun({format_table(get_flipadmix())})
output$addleaf = dtfun({format_table(get_addleaf())})
output$initresample = dtfun({format_table(get_initresample())})
output$indresample = dtfun({format_table(get_indresample())})
output$snpresample = dtfun({format_table(get_snpresample())})
output$f2stats = dtfun({format_table(get_f2())})
output$f3stats = dtfun({format_table(get_f3())})
output$f4stats = dtfun({format_table(get_f4())})
output$qpadm_weights = dtfun({format_table(get_qpadm()$weights)})
output$qpadm_f4 = dtfun({format_table(get_qpadm()$f4)})
output$qpadm_rankdrop = dtfun({format_table(get_qpadm()$rankdrop)})
output$qpadm_popdrop = dtfun({format_table(get_qpadm()$popdrop)})
output$qpadm_rot = dtfun({format_table(get_qpadm_rotate())})
onBookmark(function(state) {
names(global) %>% map(~{state$values[[.]] = global[[.]]})
})
onRestore(function(state) {
names(state$values) %>% map(~{global[[.]] = state$values[[.]]})
})
#setBookmarkExclude('sidebarItemExpanded')
#setBookmarkExclude('navbar')
observeEvent(input$tooltips, {
print('refresh tooltips')
if(input$tooltips %% 2 == 1) {
updateActionButton(session, 'tooltips', icon = icon('exclamation'))
imap(tt, ~removeTooltip(session, .y))
imap(tt, ~addTooltip(session, .y, .x))
} else {
updateActionButton(session, 'tooltips', icon = icon('question'))
imap(tt, ~removeTooltip(session, .y))
}
})
}
shinyApp(ui = ui, server = server, enableBookmarking = 'server')
#runApp(list(ui = ui, server = server), display.mode = 'showcase')
uqrmaie1/admixtools documentation built on March 20, 2024, 8:24 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
library(admixtools)
library(igraph)
library(htmlwidgets)
library(shiny)
library(shinythemes)
library(plotly)
library(tidyverse)
library(magrittr)
library(RColorBrewer)
library(shinyjs)
library(shinyBS)
library(shinyFiles)
library(shinyalert)
library(shinyWidgets)
library(shinydashboard)
#options(shiny.sanitize.errors = FALSE)
edgemul = 1000
dtstyle = 'font-size:70%'
pl = 15
do = list(pageLength = pl, lengthMenu = list(c(pl, -1), c(pl, 'All')),
buttons = c('copy', 'csv', 'excel'), scrollX = TRUE, dom = 'Blfrtip')
actionButton = function(...) actionBttn(..., style = 'unite', size = 'sm')
#cols = gg_color_hue(5, 97, 15)
cols = c("#FFF1F0", "#F7F8E4", "#E2FDF2", "#E8F9FF", "#FFF1FF")
mutfuns = c(`SPR leaves`='spr_leaves', `SPR all`='spr_all', `Swap leaves`='swap_leaves', `Move admixture edge`='move_admixedge_once', `Flip admix`='flipadmix_random', `Change root`='place_root_random', `Mutate n`='mutate_n')
box = function (..., title = NULL, footer = NULL, status = NULL, solidHeader = FALSE,
background = NULL, width = 6, height = NULL, collapsible = FALSE,
collapsed = FALSE)
{
boxClass <- "box"
if (solidHeader || !is.null(background)) {
boxClass <- paste(boxClass, "box-solid")
}
if (!is.null(status)) {
validateStatus(status)
boxClass <- paste0(boxClass, " box-", status)
}
if (collapsible && collapsed) {
boxClass <- paste(boxClass, "collapsed-box")
}
# if (!is.null(background)) {
# validateColor(background)
# boxClass <- paste0(boxClass, " bg-", background)
# }
style <- NULL
if (!is.null(height)) {
style <- paste0("background:", background, "; height: ", validateCssUnit(height))
}
titleTag <- NULL
if (!is.null(title)) {
titleTag <- h3(class = "box-title", title)
}
collapseTag <- NULL
if (collapsible) {
buttonStatus <- status %OR% "default"
collapseIcon <- if (collapsed)
"plus"
else "minus"
collapseTag <- div(class = "box-tools pull-right", tags$button(class = paste0("btn btn-box-tool"),
`data-widget` = "collapse", shiny::icon(collapseIcon)))
}
headerTag <- NULL
if (!is.null(titleTag) || !is.null(collapseTag)) {
headerTag <- div(class = "box-header", titleTag, collapseTag)
}
div(class = if (!is.null(width))
paste0("col-sm-", width), div(class = boxClass, style = if (!is.null(style))
style, headerTag, div(class = "box-body", style = paste('background:', background), ...), if (!is.null(footer))
div(class = "box-footer", footer)))
}
tt = c('data' = 'Load data here',
'dir' = 'Directory where precursors for f-statistics will be stored or have been stored.',
'popfilediv' = 'Select a PLINK or EIGENSTRAT sample metadata textfile.',
'graphfilediv' = 'A graph in ADMIXTOOLS format, or as a two column edge list.',
'genofile1' = 'Select a Packedancestrymap .geno file to extract data from. This step should only be necessary once. After that, the extracted data can be used for any further analyses.',
'adjpopdiv' = 'Assign individuals to populations. Population names can be changed, but should match population names in the graph file, if one is provided. Data will be read from disk and f2-statistics will be computed when switching to another sidebar item.',
'bookmarkdiv' = 'Click this to save the current session. This will write files to disk and generate a URL which can be used to restore the state of the progam.',
'qpgraph_fit' = 'Click this to see qpgraph fit statistics',
'qpgraph_similar_minus1' = 'Fit all graphs that result from removing one admixture edge from the current graph',
'qpgraph_similar_minusplus' = 'Fit all graphs that result from removing one admixture edge from the current graph and adding back one admixture edge',
'qpgraph_similar_plus1' = 'Fit all graphs that result from adding one admixture edge to the current graph',
'qpgraph_similar_decomposed' = 'Fit all trees that result from removing one edge from each admixture node',
'qpgraph_similar_treeneighbors' = 'Fit all trees that result from removing one edge from each admixture node and then reattaching each edge to each other edge',
'qpgraph_similar_flipadmix' = 'Fit all graphs that result from flipping the direction of each admixture edge (unless this would introduce a cycle)',
'qpgraph_update' = 'Set the current graph to the displayed graph',
'qpgraph_revert' = 'Display the original graph',
'addleaf' = 'Choose a population which should be added to the current graph',
'qpgraph_add' = 'Fit all graphs that result from adding this population to each edge',
'optim_run' = 'Run optimization. The current graph is modified in a number of different ways, and each modified graph is evaluated. Can be run iteratively across several generations.',
'optim_ngen' = 'Number of generations. In each generation the best fitting graphs are selected for the next generation.',
'optim_sec' = 'Minimum number of seconds to run',
'mutfuns' = 'SPR leaves: Subree prune and regraft leaf edges\n
SPR all: Subree prune and regraft any edges\\n
Swap leaves: Swap two random leaves\\\n
Move admixture edge: Shift one admixture edge\\\\n
Flip admix: Flip the direction of one admixture edge\\\\\n
Change root: Place the root of the tree on a random edge\\\\\\n
Mutate n: Combine two or more of the other mutation functions',
'econstraints_update_div' = 'Set limits to minimum and maximum edge lengths',
'aconstraints_update_div' = 'Set limits to admixture weights',
'qpgraphbin' = 'Choose the qpGraph executable on your computer',
'qpadmbin' = 'Choose the qpAdm executable on your computer',
'qpgraph_genofile' = 'Choose the .geno or .bed file for qpGraph to use',
'qpadm_genofile' = 'Choose the .geno or .bed file for qpAdm to use',
'usef3precomp' = 'Differences between the original and R version of qpGraph can be due to the estimation of f-statistics, or due to the fitting step. Click this to use the f-statistics generated in the original qpGraph to use in qpGraph in R.',
'init_resample' = 'Re-evaluate graph with different random intializations of starting weights. This is always done by default, but only the parameters for the starting weights that minimize the score are shown. Here, the table will show all estimated graph parameters (score, edge lengths, admixture weights) for all sets of random starting parameters.',
'ind_resample' = 'Re-evaluate the graph leaving out each sample at a time. Samples are only left out if they are not the only sample in their respective population. This can take long, as data has to be read and f2 computed separately for each sample that is being left out.',
'snp_resample' = 'Re-evaluate the graph leaving out SNP blocks. With jackknife, each SNP block is left out at a time. With bootstrap, SNP blocks are resampled with replacement to match the original number of SNP blocks (usually around 700). The number of random resamplings can be specified.',
'clear_edge' = 'Clear the current selection',
'delete_edge' = 'Delete the selected edge or node. Has to be an admixture edge or leaf node.',
'add_edge' = 'Add an admixture edge by connecting two selected non-admixture edges. Fails if it introduces cycles or attempts to connect admixture edges.',
'randgraph' = 'Generate a random admixture graph. Uses the current populations, a specified number of admixture nodes, and a specified outgroup.',
'lsqmode' = 'Least squares mode. By default, the computation of the graph likelihood score uses the inverse of the f3-statistic covariance matrix. When fitting many populations, this inverse can be unstable. Least squares mode sets the offdiagonal elements of the f3-statistic covariance matrix to zero, which makes the matrix inversion stable, but can introduce bias.',
'lambdascale' = 'Scale the estimated f2-statistics by a constant factor.',
'seed' = 'Set a random seed to always pick the same initial weights and to make the results reproducible.',
'multiprocess' = 'Enable parallel evaluation of graphs in multiple R processes. Faster, but may not work on all operating systems. See "A Future for R: A Comprehensive Overview" for more details.',
'collapse_edges' = 'Collapse all nodes which are separated by less then a certain threshold of estimated drift.',
'f2corr' = 'Untick this to compute f-statistics without bias correction factor.',
'addpop' = 'Add a population',
'removepop' = 'Remove a population',
'showseledge' = 'Selected edges or nodes are shown here',
'downloadgraph' = 'Download the fit of the current graph as a .tsv file.',
'numstart' = 'How often optimization should be initialized with random weights',
'qpgraph_diag' = 'Factor to add to diagonal elements of f3-statistic covariance matrix before inversion. Increase this value if optimization doesn\'t converge',
'qpadm_fit' = 'Return to the qpadm fit statistics',
'qpadm_diag' = 'Factor to add to diagonal elements of f4-statistic covariance matrix before inversion. Increase this value if optimization doesn\'t converge',
'qpadm_constrained' = 'Click this to constrain the admixture weights to be between 0 and 1',
'data_load' = 'Only a placeholder at the moment. Data is loaded when switching from the "Data" tab to another tab',
'b' = '',
'c' = '',
'f2' = 'These are the estimated f2 statistics',
'f3' = 'These are the estimated and fitted f3 statistics',
'opt' = 'This table contains information about the admixture weight optimization step for each iteration, such as the randomly drawn initial weights, as well as the final score ("value", the quantity that is being minimized) and convergence parameters from the R "optim" function. Run "Resample / Initial weights" to see how other estimated graph parameters depend on the initial weights.')
format_table = function(dat) {
if(is.null(dat)) return(tibble())
print('ft1')
print(names(dat))
popcols = which(str_sub(names(dat), 1, 3) == 'pop')
names = unique(c(as.matrix(dat[,popcols])))
print('ft2')
nam = shortest_unique_prefixes(names, 3)
pcols = which(names(dat) %in% c('p', 'p.value'))
print('ft3')
numcols = setdiff(which(sapply(dat, class) == 'numeric'), pcols)
print('ft4')
if(length(popcols) > 0) dat %<>% mutate_at(popcols, ~recode(., !!!purrr::set_names(nam, names)))
if('graph' %in% names(dat)) dat %<>% select(-graph)
if('edges' %in% names(dat)) dat %<>% select(-edges)
print('ft5')
out = dat %>%
mutate_at(numcols, ~round(., max(2, 2 - ceiling(log10(max(abs(na.omit(.)))))))) %>%
mutate_at(pcols, format.pval) %>%
rename_at(vars(pcols), ~'p.value')
out
}
aa = '.shiny-notification { height: 100px; width: 400px; position: fixed; top: 20%; left: calc(50% - 200px);;}'
ui = function(request) {
tagList(useShinyjs(),
useShinyalert(),
chooseSliderSkin('Modern', color = 'black'),
bsTooltip('tooltips', 'Click this to turn tooltips on or off. Needs to be reactivated after switching tabs.'),
dashboardPage(
dashboardHeader(title = 'ADMIXTOOLS 2', titleWidth = 300,
tags$li(class = 'dropdown', actionLink('tooltips', label = '', icon = icon('question'))),
tags$li(class = 'dropdown', id = 'bookmarkdiv', actionLink('._bookmark_', label = '', icon = icon('save')))),
dashboardSidebar(tags$head(tags$style(HTML(aa))), width = 300,
sidebarMenu(id = 'navbar',
menuItem('Data', tabName = 'data', expandedName = 'data', id = 'datax', icon = icon('database'), startExpanded = TRUE,
menuItem('Extract settings', tabName = 'extract_settings',
checkboxInput('fix_populations', 'Fix populations', TRUE),
numericInput('max_miss', 'max missing', value = 0, step = 0.01, min = 0, max = 1),
splitLayout(
numericInput('min_maf', 'min MAF', value = 0, step = 0.01, min = 0, max = 0.5),
numericInput('max_maf', 'max MAF', value = 0.5, step = 0.01, min = 0, max = 0.5)),
radioButtons('trans_extract', 'Mutations', choices = c('both', 'only transitions', 'only transversions')),
fileInput('keepsnps', NULL, buttonLabel = 'SNP list'),
numericInput('maxmem', 'max RAM in GB', value = 15, step = 1, min = 0)
)),
menuItem('f2', tabName = 'f2', expandedName = 'f2', icon = icon('dice-two'),
menuItem('Options', tabName = 'f2_options')),
menuItem('f3', tabName = 'f3', expandedName = 'f3', icon = icon('dice-three'),
menuItem('Options', tabName = 'f3_options')),
menuItem('f4', tabName = 'f4', expandedName = 'f4', icon = icon('dice-four'),
menuItem('Options', tabName = 'f4_options')),
menuItem('qpAdm', tabName = 'qpadm', expandedName = 'qpadm', icon = icon('balance-scale'),
actionLink('qpadm_fit', 'Fit'),
menuItem('Compare', href = 'qpadm_comparison', tabName = 'qpadm_comparison',
splitLayout(
shinyFilesButton('qpadmbin', 'qpAdm bin', 'Select qpAdm bin', FALSE),
shinyFilesButton('qpadm_genofile', 'Geno file', 'Select Packedancestrymap geno file', FALSE)),
p(HTML('<center>or</center>')),
fileInput('qpadm_out', NULL, buttonLabel = 'qpAdm output'),
hr(),
actionButton('run_qpadm', 'Run')),
menuItem('Resample', tabName = 'qpadm_resample'),
menuItem('Options', tabName = 'qpadm_options',
numericInput('qpadm_diag', 'diag', value = 0.0001, step = 0.0001),
checkboxInput('qpadm_constrained', 'Constrain weights'),
checkboxInput('qpadm_results', 'Show results', value = TRUE),
checkboxInput('qpadm_rotate', 'Rotate left', value = FALSE))),
menuItem('qpGraph', tabName = 'qpgraph', expandedName = 'qpgraph', id = 'qpgraph', icon = icon('project-diagram'),
checkboxInput('evaluate_graph', 'Evaluate graph', FALSE),
menuItem('Fit',
splitLayout(
shiny::actionButton('qpgraph_update', 'Update graph'),
shiny::actionButton('qpgraph_revert', 'Revert graph')),
#actionButton('options_update', 'Update'),
splitLayout(numericInput('qpgraph_diag', 'diag', value = 0.0001, step = 0.0001),
numericInput('numstart', '# init', value = 10, min = 1)),
splitLayout(numericInput('lambdascale', 'lambdascale', value = 1, step = 0.001),
numericInput('seed', 'Random seed', value = NULL)),
splitLayout(checkboxInput('lsqmode', 'lsqmode'),
checkboxInput('return_fstats', 'return_fstats'))),
menuItem('Modify', expandedName = 'qpgraph_modify', id = 'qpgraph_modify', tabName = 'qpgraph_modify',
hr(),
p('Selected'),
div(verbatimTextOutput('showseledge', placeholder = TRUE), tags$head(tags$style('#showseledge{max-height: 300px; background: white}'))),
splitLayout(
shiny::actionButton('clear_edge', 'Clear'),
shiny::actionButton('delete_edge', 'Delete'),
shiny::actionButton('add_edge', 'Add')),
uiOutput('qpgraph_add'),
hr(),
p('Randomize graph'),
uiOutput('nadmix'),
#splitLayout(checkboxInput('fixoutgroup', 'Fix outgroup', value = TRUE),
actionButton('randgraph', 'Randomize'),
hr()),
menuItem('Optimize', tabName = 'qpgraph_optim',
actionButton('optim_run', 'Run'),
numericInput('optim_ngen', '# generations', 10, 1, 1000),
numericInput('optim_sec', 'min seconds', 10, 1, 1000),
checkboxInput('optim_initrand', 'Start randomly'),
checkboxGroupInput('mutfuns', 'Graph modifications', choices = mutfuns, selected = mutfuns)),
menuItem('Constrain', tabName = 'qpgraph_constraints',
menuItem('Drift', tabName = 'qpgraph_econstraints',
uiOutput('econstraints')),
menuItem('Admixture', tabName = 'qpgraph_aconstraints',
uiOutput('aconstraints'))),
menuItem('Compare', tabName = 'qpgraph_comparison',
splitLayout(
shinyFilesButton('qpgraphbin', 'qpGraph bin', 'Select qpGraph bin', FALSE),
shinyFilesButton('qpgraph_genofile', 'Geno file', 'Select Packedancestrymap geno file', FALSE)),
checkboxInput('usef3precomp', 'Use qpGraph fstats in R'),
p(HTML('<center>or</center>')),
fileInput('qpgraph_out', NULL, buttonLabel = 'qpGraph output'),
hr(),
actionButton('run_qpgraph', 'Run')),
menuItem('Resample', tabName = 'qpgraph_resample',
menuItem('Initial weights', tabName = 'qpgraph_resample_init',
sliderInput('numstart2', '# init', 1, 100, 10),
actionButton('init_resample', 'Run')),
menuItem('Individuals', tabName = 'qpgraph_resample_individuals',
actionButton('ind_resample', 'Run')),
menuItem('SNPs', tabName = 'qpgraph_resample_snps',
checkboxInput('boot', 'Bootstrap'),
conditionalPanel('input.boot',
sliderInput('bootnum', '# boot', 0, 1000, 100)),
actionButton('snp_resample', 'Run'))),
menuItem('Load', fileInput('graphfile', NULL, placeholder = '', buttonLabel = 'Graph file')),
menuItem('Save', tabName = 'qpgraph_download', expandedName = 'qpgraph_download',
div(radioButtons('downloadgraph_format', 'Graph format',
choices = c(`ADMIXTOOLS`='admixtools', `Edge list`='edgelist', `DOT`='dot'),
selected = 'admixtools'),
downloadButton('downloadgraph', 'Save graph'))),
menuItem('Plot options',
checkboxGroupInput('plotopt', '',
choices = c(`Highlight unidentifiable` = 'highlight_unidentifiable',
`Reorder edges`='reorder_edges',
`Shift edges down`='shift_down',
`Collapse edges`='collapse_edges'),
selected = c('shift_down')),
sliderInput('collapse_threshold', 'Log10 collapse threshold', -6, 2, -3, step = 0.1))),
menuItem('Options', tabName = 'options', expandedName = 'options', id = 'options', icon = icon('cogs'),
checkboxInput('multiprocess', 'multiprocess', value = FALSE),
checkboxInput('f2corr', 'f2 bias correction', value = TRUE),
checkboxInput('afprod', 'allele frequency products', value = FALSE)))),
dashboardBody(uiOutput('dashboardbody'), tags$head(tags$style(HTML('.skin-black .sidebar .shiny-download-link a { color: #444; } .content {background-color: white } .main-header .logo {font-family: "GillSans", "Skia", "Avenir-Medium", "Herculanum", "Hiragino Maru Gothic Pro", "Maru Gothic Pro", "Hiragino", "Comic Sans MS"; font-weight: normal;};')))),
skin = 'black'
))
}
server = function(input, output, session) {
showLog()
dp = '~/mnt/o2/projects/admixprograms/indpairs_v42.1/'
dp = '~/Downloads/testinds/'
dp = '/Users/robert/Downloads/countdata_eas5/'
dp = '~/Downloads/countdata_iosifgraph/'
dp = ''
#dp = '~/Documents/countdata_afafam/'
#dp = '~/Documents/countdata_pavel2/'
global = reactiveValues(countdir = dp,
allinds = list.dirs(paste0(dp, '/pairs'),F,F),
poplist = list())
global$graph = NULL
home = normalizePath('~')
volumes = getVolumes()
#volumes = expr(c(home = normalizePath('~'), getVolumes()() %>% purrr::set_names()))
shinyDirChoose(input, 'dir', roots = volumes)
shinyFileChoose(input, 'genofile1', roots = volumes)
shinyFileChoose(input, 'qpgraph_genofile', roots = volumes)
shinyFileChoose(input, 'qpadm_genofile', roots = volumes)
shinyFileChoose(input, 'qpgraphbin', roots = volumes)
shinyFileChoose(input, 'qpadmbin', roots = volumes)
#global$show_qpadm_rotate = FALSE
#global$show_qpadm_results = TRUE
shinyjs::disable('qpgraph_update')
shinyjs::disable('qpgraph_revert')
output$dirout = renderText({global$countdir})
output$dirout = renderText({str_replace(global$countdir, '/Users/robert', '')})
output$show_popadjust = renderText({as.numeric(show_popadjust())})
output$show_extract = renderText({as.numeric(show_extract())})
output$show_indselect = renderText({as.numeric(show_indselect())})
output$genofile1out = renderText({parseFilePaths(eval(volumes), input$genofile1)$datapath %>% str_replace('\\.geno$', '') %>% str_replace('/Users/robert', '')})
observeEvent(input$dir, {
#if (!"path" %in% names(dir())) return()
#volumes = eval(volumes)
print('volumes')
print(volumes)
print('input$dir')
print(c(input$dir, class(input$dir)))
print('xxx')
print(parseDirPath(volumes, input$dir))
print(str(parseDirPath(volumes, input$dir)))
print(as.character(parseDirPath(volumes, input$dir)))
global$countdir = parseDirPath(volumes, input$dir)
global$iscountdata = 'indivs' %in% list.dirs(global$countdir,F,F)
global$isf2data = !'indivs' %in% list.dirs(global$countdir,F,F) &
'block_lengths_f2.rds' %in% list.files(global$countdir)
print('global$iscountdata')
print(global$iscountdata)
if(global$iscountdata) {
global$allinds = list.dirs(paste0(global$countdir, '/pairs'),F,F)
} else if(global$isf2data) {
global$allinds = list.dirs(global$countdir,F,F)
global$poplist = global$allinds %>% purrr::set_names(global$allinds)
print('global$allinds')
print(global$allinds)
} else {
}
})
observeEvent(input$textdirinput, {
nam = ifelse(!is.null(input$textdirinput) && dir.exists(input$textdirinput), 'Load', 'Create')
updateActionButton(session, "textdir", label = nam)
})
observeEvent(input$textdir, {
justcreated = FALSE
if(!dir.exists(input$textdirinput)) {
#shinyalert('Could not find directory!', "Please write the path of a directory into the box. This should be an empty directory if you haven't extracted f2-statistics/counts yet, or otherwise the directory with extracted data.")
dir.create(input$textdirinput)
justcreated = TRUE
#return()
}
global$countdir = normalizePath(input$textdirinput)
global$iscountdata = 'indivs' %in% list.dirs(global$countdir,F,F)
global$isf2data = !'indivs' %in% list.dirs(global$countdir,F,F) &
'block_lengths_f2.rds' %in% list.files(global$countdir)
print('global$iscountdata')
print(global$iscountdata)
if(global$iscountdata) {
global$allinds = list.dirs(paste0(global$countdir, '/pairs'),F,F)
} else if(global$isf2data) {
global$allinds = list.dirs(global$countdir,F,F)
global$poplist = global$allinds %>% purrr::set_names(global$allinds)
print('global$allinds')
print(global$allinds)
} else {
if(!justcreated) shinyalert("Directory does not seem to have f2 data!")
}
})
observeEvent(input$navbar, {
print('navbar detected')
print(input$navbar)
})
tolisten = reactive({
list(input$sidebarItemExpanded, input$qpadm_rotate, input$qpadm_results)
#!is.null(input$sidebarItemExpanded) && (is.null(global$exp) || input$sidebarItemExpanded != global$exp)
})
get_bod = eventReactive(tolisten(), {
print('expanded!')
exp = input$sidebarItemExpanded
print(exp)
#if(is.null(exp) || !is.null(global$exp) && exp == global$exp) return(global$bod)
if(is.null(exp)) exp = global$exp
global$exp = exp
if(exp == 'options') {
} else if(exp == 'data') {
print('navbar == data!')
global$databod = get_loaddata()
global$bod = global$databod
} else {
nam = map(paste0('pop', seq_len(length(global$poplist))), ~input[[.]])
if(is.null(nam[[1]])) {
shinyalert('Global poplist is NULL or has not been processed correctly!', global$poplist)
return()
}
global$poplist = map(nam, ~input[[.]]) %>% purrr::set_names(nam)
print('switch away from data')
if(length(unlist(global$poplist)) == 0) {
shinyalert('No samples selected!', '')
return()
}
if(exp == 'f2') {
print('switch to f2')
global$bod = fluidRow(
column(1, div(radioGroupButtons('fst', choices = c('f2', 'fst')))),
column(11, tabsetPanel(tabPanel('Pop pairs', plotlyOutput(outputId = 'f2heatmap', height = '800', width='auto')),
tabPanel('Ind pairs', plotlyOutput(outputId = 'f2heatmap_indiv', height = '800', width='auto')),
tabPanel('Pop table', dto('f2stats')))))
} else if(exp == 'f3') {
print('switch to f3')
global$bod = fluidRow(
column(12, tabsetPanel(tabPanel('Pop table', dto('f3stats')))))
} else if(exp == 'f4') {
print('switch to f4')
choices = names(global$poplist)
global$f4_left = div(map(1:4, ~selectizeInput(paste0('f4pops', .),
paste0('Population ', .),
choices = choices,
multiple = TRUE,
selected = choices[seq(.*(.-1)/2+1, .*(.-1)/2+.)])),
checkboxInput('f4_allowdup', 'Allow duplicates'),
sliderInput('f4_se', 'Standard errors', min = 0, max = 5, value = 0, step = 0.01))
#checkboxInput('f4_se', 'Show standard errors'))
popsinuse = unique(c(input$f4pops1, input$f4pops2, input$f4pops3, input$f4pops4))
# update
map(1:4, ~{
fp = paste0('f4pops', .)
sel = input[[fp]]
updateSelectizeInput(session, fp, selected = sel,
choices = union(sel, setdiff(names(global$poplist), popsinuse)))})
global$bod = fluidRow(
column(3, global$f4_left),
column(9, tabsetPanel(tabPanel('Plot', plotlyOutput(outputId = 'f4plot', height = '800', width='auto')),
tabPanel('Plot2', plotlyOutput(outputId = 'f4plot2', height = '800', width='auto')),
tabPanel('Table', dto('f4stats')))))
} else if(exp == 'qpadm') {
print('switch to qpAdm')
choices = names(global$poplist)
if(is.null(input$qpadmpops1) || is.null(input$qpadmpops2) || is.null(input$qpadmpops3)) {
global$qpadmpops1 = choices[2:3]
global$qpadmpops2 = choices[4:6]
global$qpadmpops3 = choices[1]
} else {
global$qpadmpops1 = input$qpadmpops1
global$qpadmpops2 = input$qpadmpops2
global$qpadmpops3 = input$qpadmpops3
}
#global$qpadm_leftpanel = qpadm_leftpanel()
global$bod = fluidRow(
column(3, qpadm_leftpanel()),
column(9, qpadm_rightpanel()))
} else if(exp == 'qpgraph') {
print('switch to qpGraph')
if(is.null(global$graph) || !all(get_leafnames(global$graph) %in% names(global$poplist))) {
numadm = 2
if(!is.null(global$graph)) {
shinyalert('Generating random graph because populations don\'t match',
setdiff(get_leafnames(global$graph), names(global$poplist)))
numadm = numadmix(isolate(global$graph))
}
op = if(is.null(input$outpop) || input$outpop == '< undefined >') NULL else input$outpop
global$graph = random_admixturegraph(na.omit(names(global$poplist)), numadm, TRUE, outpop = op)
}
#global$qpg_right = qpg_right_fit()
#global$qpg_right_show = 'fit'
#global$qpg_right = qpg_right()
global$bod = fluidRow(
column(8, plotlyOutput(outputId = 'graphplot', height = '800', width='auto')),
column(4, tabsetPanel(tabPanel('Fit', tabsetPanel(tabPanel('f2', dto('f2')),
tabPanel('f3', dto('f3')),
tabPanel('f4', dto('f4')),
tabPanel('opt', dto('opt')))),
tabPanel('Similar', tabsetPanel(tabPanel('-1', dto('minus1')),
tabPanel('+1', dto('plus1')),
tabPanel('flip', dto('flipadmix')))),
tabPanel('History', dto('history')),
tabPanel('Compare', plotlyOutput('graphcomparison')),
tabPanel('Resample', tabsetPanel(tabPanel('Starting values', dto('initresample')),
tabPanel('Individuals', dto('indresample')),
tabPanel('SNPs', dto('snpresample')))),
id = 'qpgraph_tabset')))
} else {
shinyalert('Error', 'implement tab')
}
}
global$bod
})
observeEvent(global$graph, {
print('observe global$graph')
updateSelectInput(session, 'addleaf',
choices = setdiff(names(global$poplist), get_leafnames(global$graph)))
})
get_graph = reactive({
print('get graph')
global$graph
})
# observeEvent(global$qpg_right, {
# print('global$qpg_right change detected!')
# global$qpgraphbod = fluidRow(
# column(8, plotlyOutput(outputId = 'graphplot', height = '800', width='auto')),
# column(4, global$qpg_right))
# })
# observeEvent(global$databod, {global$bod = global$databod})
# observeEvent(global$qpadmbod, {global$bod = global$qpadmbod})
observeEvent(global$qpgraphbod, {print('global$qpgraphbod change detected!'); global$bod = global$qpgraphbod})
observeEvent(global$bod, {print('global$bod change detected!'); b = global$bod; global$bod = b})
observeEvent(input$navbar, {
print('navbar selected')
print(input$navbar)
})
observeEvent(input$addpop, {
pops = names(global$poplist)
i = 0
repeat({i=i+1; newpop = paste0('pop', i); if(!newpop %in% pops) break})
global$poplist[newpop] = NA
})
observeEvent(input$removepop, {
global$poplist = head(global$poplist, -1)
})
observeEvent(input$popfile, {
inddat = read_table2(input$popfile$datapath, col_names = FALSE, col_types = cols())
if(ncol(inddat) == 3) colnames(inddat) = c('ind', 'sex', 'pop')
else if(ncol(inddat) == 6) colnames(inddat) = c('pop', 'ind', 'p1', 'p2', 'sex', 'pheno')
else shinyalert('Error', paste0('sample file has ',ncol(inddat),' columns!'))
global$poplist = inddat %>%
select(ind, pop) %$% split(ind, factor(pop, levels = unique(pop)))
if(length(global$allinds) == 0) { # f2data dir
global$allinds = unlist(global$poplist)
}
print('observe')
imap(names(global$poplist), ~{
bttn = paste0('group', .y)
print(bttn)
observeEvent(input[[bttn]], {
isolate({
if(input[[bttn]] == 0) return()
poplist = global$poplist
#poplistbak = global$poplistbak
countdir = global$countdir
inds = poplist[[.x]]
#if(is.null(global$poplistbak[[.x]])) {
isgrouped = length(inds) == 1 && admixtools:::is_group(countdir, inds)
if(!isgrouped) {
if(.x %in% global$allinds) {
shinyalert('Error', 'Group name cannot be identical to an existing sample or group name!')
return()
}
print('add group')
withProgress(message = paste0('grouping ', length(inds), ' samples...'), {
group_samples(global$countdir, inds, .x, overwrite = TRUE)
})
global$allinds = union(global$allinds, .x)
#global$poplistbak[[.x]] = inds
global$poplist[[.x]] = .x
updateSelectizeInput(session, .x, selected = .x, choices = global$allinds)
#updateButton(session, bttn, 'ungroup')
} else {
print('delete group')
withProgress(message = paste0('ungrouping ', length(inds), ' samples...'), {
global$poplist[[.x]] = delete_groups(global$countdir, .x)
})
#global$poplist[[.x]] = global$poplistbak[[.x]]
#global$poplistbak[[.x]] = NULL
global$allinds = setdiff(global$allinds, .x)
updateSelectizeInput(session, .x, selected = global$poplist[[.x]], choices = global$allinds)
#updateButton(session, bttn, 'group')
}
})
})
})
})
observeEvent(input$graphfile, {
gf = input$graphfile$datapath
print('read graphfile')
print(gf)
colnums = read_lines(gf) %>% str_split('\\s+') %>% map_dbl(length)
if(length(unique(colnums)) > 1) {
print('admixtools format')
global$graph = parse_qpgraph_graphfile(gf) %>% graph_from_edgelist()
} else {
print('two column format')
global$graph = read_table2(gf, col_types = cols()) %>% select(1:2) %>%
as.matrix %>% graph_from_edgelist()
}
print('read graphfile done')
})
observeEvent(input$qpgraph_out, {
global$qpgraphout = parse_qpgraph_output(input$qpgraph_out$datapath)
#global$qpg_right = plotlyOutput('graphcomparison')
})
observeEvent(input$randgraph, {
print('randg')
op = if(is.null(input$outpop) || input$outpop == '< undefined >') NULL else input$outpop
global$graph = random_admixturegraph(get_leafnames(global$graph), input$nadmix, TRUE, outpop = op)
})
observeEvent(input$qpgraph_update, {
print('qpgraph_update')
global$graph = global$selgraph
fit = get_fit()
global$edges = fit$edges
global$score = fit$score
global$worst_residual = fit$worst_residual
add_graph_to_history(global)
global$tempgraph = FALSE
shinyjs::disable('qpgraph_update')
shinyjs::disable('qpgraph_revert')
})
observeEvent(input$qpgraph_revert, {
print('qpgraph_revert')
global$selgraph = global$graph
fit = get_fit()
global$edges = fit$edges
global$score = fit$score
global$worst_residual = fit$worst_residual
add_graph_to_history(global)
global$tempgraph = FALSE
shinyjs::disable('qpgraph_update')
shinyjs::disable('qpgraph_revert')
print(global$score)
})
observeEvent(input$navbar, {
print('navbar:')
print(input$navbar)
})
observeEvent(input$optim_run, {
numgen = input$optim_ngen
stop_sec = input$optim_sec
diag = input$qpgraph_diag
lambdascale = as.numeric(input$lambdascale)
f2blocks = get_f2blocks()
pops = dimnames(f2blocks)[[1]]
f2blocks = f2blocks[pops, pops,]
outpop = get_outpop(global$graph)
nadmix = numadmix(global$graph)
g = NULL
if(!input$optim_initrand) g = global$graph
print('running opt...')
print(pops)
print(dim(f2blocks))
print(paste(numgen, nadmix))
mutfuns = sapply(input$mutfuns, get)
withProgress(message = paste('Running find_graphs...\nCheck R console to see progress'), {
# opt_results = find_graphs(f2blocks, outpop = outpop, numrep = numrep,
# numgraphs = numgraphs, numgen = numgen, numsel = numsel,
# numadmix = nadmix, initgraph = g, mutfuns = mutfuns, diag = diag)
opt_results = find_graphs(f2blocks, outpop = outpop, stop_gen = numgen, stop_sec = stop_sec,
numadmix = nadmix, initgraph = g, mutfuns = mutfuns, diag = diag)
})
print(opt_results)
opt_results %<>% filter(!is.na(score))
print('opt done')
if(nrow(opt_results) == 0) {
shinyalert('No valid graphs in output!', 'Try increasing the "diag" parameter.')
} else {
winner = opt_results %>% slice_min(score) %$% graph[[1]]
oldgraph = global$graph
oldscore = global$score
oldedges = global$edges
oldwr = global$worst_residual
withProgress(message = 'Fitting best graph...', {
global$graph = winner
fit = get_fit()
})
print('fit')
print(str(fit))
if(fit$score < oldscore) {
shinyalert('Better graph found!', paste('Old score: ', round(oldscore, 2),
'\nNew score: ', round(fit$score, 2)))
} else {
global$graph = oldgraph
global$score = oldscore
global$edges = oldedges
global$worst_residual = oldwr
add_graph_to_history(global)
fit = get_fit()
shinyalert('No better graph found!')
}
print('graph')
print(global$graph)
}
global$seed = NULL
})
qpgraphfun = reactive({
req(input$evaluate_graph)
print('qpgraphfun')
seed = isolate(as.numeric(input$seed))
if(is.na(seed)) seed = global$seed
if(input$usef3precomp && !is.null(global$precomp)) f3precomp = global$precomp
else f3precomp = NULL
#isolate({
numstart = as.numeric(input$numstart)
lambdascale = as.numeric(input$lambdascale)
if(is.null(lambdascale) || length(lambdascale) == 0) lambdascale = 1
diag = input$qpgraph_diag
lsqmode = input$lsqmode
return_fstats = input$return_fstats
print(paste('qpgraphfun:', numstart, seed, diag, lsqmode, lambdascale))
function(x, y, ...) {
args = list(...)
if(!'numstart' %in% names(args)) args[['numstart']] = numstart
args = c(list(x, y), args, diag = diag, lambdascale = lambdascale, lsqmode = lsqmode,
seed = seed, f3precomp = list(f3precomp), return_fstats = return_fstats)
do.call(quietly(qpgraph), args)$result
}
#})
})
get_fit = reactive({
print('get_fit')
if(!input$evaluate_graph) return(list(edges=global$edges, score=global$score, worst_residual=global$worst_residual))
f2blocks = get_f2blocks()
g = global$graph
input$qpgraph_revert
input$qpgraph_update
input$aconstraints_update
input$econstraints_update
global$qpgraph_ranges = NULL
global$qpgraph_scorerange = NULL
req(g, f2blocks)
replace_null = function(x, y) ifelse(is.null(x), y, x)
edges = g %>% igraph::as_edgelist() %>% as_tibble(.name_repair = ~c('from', 'to')) %>%
group_by(to) %>% mutate(type = ifelse(n() > 1, 'admix', 'normal')) %>% ungroup %>%
mutate(eid = paste(from, to, sep = ' -> '), eid2 = str_remove_all(eid, '\\.|>| '))
isolate({
if(isTRUE(global$useconstraints)) {
edges %<>% mutate(lower = map_dbl(eid2, ~replace_null(input[[.]][1], 0)),
upper = map_dbl(eid2, ~replace_null(input[[.]][2], 1e9))) %>%
mutate(upper = ifelse(type == 'normal', upper/edgemul, upper),
lower = ifelse(type == 'normal', lower/edgemul, lower)) %>% select(-type, -eid2)
#print(edges)
}
})
withProgress(message = 'Fitting graph...', {
fit = qpgraphfun()(f2blocks, edges)
})
global$edges = fit$edges
global$score = fit$score
global$worst_residual = fit$worst_residual
add_graph_to_history(global)
fit
})
get_history = reactive({
print('get_history')
isolate({
hist = global$history %>% purrr::transpose() %>% as_tibble %>%
transmute(graph, edges,
nadmix = as.numeric(nadmix),
score = as.numeric(score),
worst_residual = map_dbl(worst_residual, ~ifelse(is.null(.), NA, .))) %>%
mutate(n = 1:n(), .before = 1) %>% arrange(-n)
})
hist %>% mutate(x = list(get_fit())) %>% select(-x)
})
get_pdat = reactive({
print('get_pdat')
if(is.null(global$edges)) return()
print('get_pdat2')
init_plotdata()
})
get_f2blocks = reactive({
print('get f2')
poplist = global$poplist
req(poplist, global$countdir)
if(length(poplist) == 0) return()
pops = rep(names(poplist), sapply(poplist, length))
inds = NULL
if(global$iscountdata) inds = unlist(poplist)
withProgress(message = 'Reading data and computing f2...', {
f2blocks = f2_from_precomp(global$countdir, inds = inds, pops = pops,
apply_corr = input$f2corr, afprod = input$afprod)
})
print('get f2 done')
f2blocks
})
get_fstblocks = reactive({
print('get fst')
poplist = global$poplist
req(poplist, global$countdir)
if(length(poplist) == 0) return()
pops = rep(names(poplist), sapply(poplist, length))
inds = NULL
if(global$iscountdata) inds = unlist(poplist)
withProgress(message = 'Reading data and computing fst...', {
f2blocks = f2_from_precomp(global$countdir, inds = inds, pops = pops,
apply_corr = input$f2corr, afprod = input$afprod, fst = TRUE)
})
print('get fst done')
f2blocks
})
get_f2dat_indiv = reactive({
print('get f2 indiv')
poplist = global$poplist
req(poplist, global$countdir)
if(length(poplist) == 0) return()
withProgress(message = 'Reading data and computing f2...', {
f2dat = f2_from_precomp(global$countdir, inds = unlist(poplist), pops = unlist(poplist),
return_array = FALSE, afprod = input$afprod)
})
print('get f2 indiv done')
f2dat
})
init_plotdata = reactive({
print('init')
input$qpgraph_update
input$qpgraph_revert
input$evaluate_graph
graph = global$edges
if('igraph' %in% class(graph)) {
edges = graph %>% as_edgelist %>% as_tibble(.name_repair = ~c('from', 'to'))
} else {
edges = graph %>% as_tibble(.name_repair = make.unique)
if(nrow(edges) == 0) browser()
graph %<>% select(1:2) %>% as.matrix %>% graph_from_edgelist()
}
names(edges)[1:2] = c('from', 'to')
if('collapse_edges' %in% input$plotopt) {
withProgress(message = 'collapsing edges...', {
edges %<>% admixtools:::collapse_edges(10^input$collapse_threshold)
})
}
admixnodes = unique(edges[[2]][edges[[2]] %in% names(which(table(edges[[2]]) > 1))])
graph = edges %>% select(1:2) %>% as.matrix %>% graph_from_edgelist()
pos = data.frame(names(V(graph)), igraph::layout_as_tree(graph), stringsAsFactors = F) %>%
set_colnames(c('node', 'x', 'y'))
eg = edges %>% left_join(pos, by=c('from'='node')) %>%
left_join(pos %>% transmute(to=node, xend=x, yend=y), by='to') %>%
mutate(type = ifelse(to %in% admixnodes, 'admix', 'normal'))
if('reorder_edges' %in% input$plotopt) {
withProgress(message = 'reordering edges...', {
eg = admixtools:::fix_layout(eg %>% rename(name = from), graph) %>% rename(from = name)
})
}
if('shift_down' %in% input$plotopt) {
withProgress(message = 'reordering edges...', {
eg = admixtools:::fix_shiftdown(eg %>% rename(name = from), graph) %>% rename(from = name)
})
}
# if('simplify_graph' %in% input$plotopt) {
# withProgress(message = 'simplifying...', {
# isolate({
# global$graph = admixtools:::edges_to_igraph(eg) %>% admixtools:::simplify_graph()
# global$edges = edges = global$graph %>%
# igraph::as_edgelist() %>% as_tibble(.name_repair = ~c('from', 'to'))
# })
# eg = edges %>%
# left_join(pos, by=c('from'='node')) %>%
# left_join(pos %>% transmute(to=node, xend=x, yend=y), by='to') %>%
# mutate(type = ifelse(to %in% admixnodes, 'admix', 'normal'))
# })
# }
if('weight' %in% names(edges)) {
fa = function(x) paste0(round(x*100), '%')
fe = function(x) round(x*edgemul)
edges$label = ifelse(edges$type == 'admix', fa(edges$weight), fe(edges$weight))
if(!is.null(global$qpgraph_ranges)) {
edges %<>% left_join(global$qpgraph_ranges, by = c('from', 'to')) %>%
mutate(label = ifelse(type == 'admix',
paste(fa(mid), paste0('[', fa(lo), '-', fa(hi), ']'), sep = '\n'),
paste(fe(mid), paste0('[', fe(lo), '-', fe(hi), ']'), sep = '\n')))
}
}
if(!'label' %in% names(edges)) edges %<>% mutate(label='')
e2 = edges %>% transmute(from, to, label) %>% left_join(count(., to), by = c('from'='to'))
eg %<>% left_join(e2, by=c('from', 'to')) %>% mutate(indegree = replace_na(n, 0))
nodes = eg %>% filter(to %in% get_leafnames(graph)) %>% rename(x=xend, y=yend, xend=x, yend=y) %>%
transmute(name = to, x, y, xend, yend, to=NA, rownum = 1:n())
namedList(nodes, eg)
})
show_popadjust = reactive(length(global$poplist) > 0)
show_extract = reactive(isFALSE(global$iscountdata) && isFALSE(global$isf2data))
show_indselect = reactive(show_extract() || isTRUE(global$iscountdata))
observeEvent(input$minus1_cell_clicked, {
row = input$minus1_rows_selected
req(row)
sel = get_minus1() %>% slice(row)
global$selgraph = sel$graph[[1]]
global$edges = sel$edges[[1]]
global$score = sel$score[[1]]
global$worst_residual = sel$worst_residual[[1]]
add_graph_to_history(global)
global$tempgraph = TRUE
shinyjs::enable('qpgraph_update')
shinyjs::enable('qpgraph_revert')
print('newgraph set...')
print(numadmix(global$graph))
})
observeEvent(input$plus1_cell_clicked, {
row = input$plus1_rows_selected
req(row)
sel = get_plus1() %>% slice(row)
global$selgraph = sel$graph[[1]]
global$edges = sel$edges[[1]]
global$score = sel$score[[1]]
global$tempgraph = TRUE
shinyjs::enable('qpgraph_update')
shinyjs::enable('qpgraph_revert')
})
observeEvent(input$flipadmix_cell_clicked, {
row = input$flipadmix_rows_selected
req(row)
sel = get_flipadmix() %>% slice(row)
global$selgraph = sel$graph[[1]]
global$edges = sel$edges[[1]]
global$score = sel$score[[1]]
global$worst_residual = sel$worst_residual[[1]]
add_graph_to_history(global)
global$tempgraph = TRUE
shinyjs::enable('qpgraph_update')
shinyjs::enable('qpgraph_revert')
})
observeEvent(input$history_cell_clicked, {
row = input$history_rows_selected
req(row)
sel = get_history() %>% slice(row)
global$selgraph = sel$graph[[1]]
global$edges = sel$edges[[1]]
global$score = sel$score[[1]]
global$worst_residual = sel$worst_residual[[1]]
#add_graph_to_history(global)
global$tempgraph = TRUE
shinyjs::enable('qpgraph_update')
shinyjs::enable('qpgraph_revert')
})
observeEvent(input$addleaf_cell_clicked, {
row = input$addleaf_rows_selected
req(row)
sel = get_addleaf() %>% slice(row)
global$selgraph = sel$graph[[1]]
global$edges = sel$edges[[1]]
global$score = sel$score[[1]]
global$worst_residual = sel$worst_residual[[1]]
add_graph_to_history(global)
global$tempgraph = TRUE
shinyjs::enable('qpgraph_update')
shinyjs::enable('qpgraph_revert')
})
observeEvent(input$initresample_cell_clicked, {
row = input$initresample_rows_selected
req(row)
sel = get_initresample0() %>% slice(row)
global$edges = sel$edges[[1]]
global$score = sel$score[[1]]
global$worst_residual = sel$worst_residual[[1]]
add_graph_to_history(global)
global$tempgraph = TRUE
shinyjs::enable('qpgraph_update')
shinyjs::enable('qpgraph_revert')
})
observeEvent(input$indresample_cell_clicked, {
row = input$indresample_rows_selected
req(row)
sel = get_indresample0() %>% slice(row)
global$edges = sel$edges[[1]]
global$score = sel$score[[1]]
global$worst_residual = sel$worst_residual[[1]]
add_graph_to_history(global)
global$tempgraph = TRUE
shinyjs::enable('qpgraph_update')
shinyjs::enable('qpgraph_revert')
})
observeEvent(input$snpresample_cell_clicked, {
row = input$snpresample_rows_selected
req(row)
sel = get_snpresample0() %>% slice(row)
global$edges = sel$edges[[1]]
global$score = sel$score[[1]]
global$worst_residual = sel$worst_residual[[1]]
add_graph_to_history(global)
global$tempgraph = TRUE
shinyjs::enable('qpgraph_update')
shinyjs::enable('qpgraph_revert')
})
observeEvent(input$qpadm_rot_cell_clicked, {
row = input$qpadm_rot_rows_selected
req(row)
sel = get_qpadm_rotate() %>% slice(row)
global$qpadmpops_left_remain = sel$left[[1]]
global$qpadmpops_right_add = sel$right[[1]]
})
observeEvent(input$multiprocess, {
print('multiprocess')
if(input$multiprocess) future::plan('multiprocess')
else future::plan('sequential')
})
make_reactive_graphtabs = function(graphfun, inp) {
reactive({
print('graph neighbors')
f2blocks = get_f2blocks()
print('graph neighbors2')
#g = global$graph
g = get_graph()
withProgress(message = 'Finding graphs...', {
out = g %>% graphfun
})
print(paste('found', nrow(out), 'graphs'))
withProgress(message = paste('Evaluating', nrow(out), 'graphs...'), {
qpgfun = qpgraphfun()
out %<>% mutate(res = furrr::future_map(graph, ~qpgfun(f2blocks, .))) %>%
unnest_wider(res) %>%
select(starts_with('from'), starts_with('to'), starts_with('score'), starts_with('graph'), starts_with('edges'))
})
if(nrow(out) == 0) {
shinyalert('No graphs could be evaluated!')
return()
}
out %>% arrange(score)
})
}
get_minus1 = make_reactive_graphtabs(graph_minusone, input$qpgraph_similar_minus1)
get_plus1 = make_reactive_graphtabs(graph_plusone, input$qpgraph_similar_plus1)
get_flipadmix = make_reactive_graphtabs(admixtools::graph_flipadmix, input$qpgraph_similar_flipadmix)
get_addleaf = eventReactive(input$qpgraph_add_run, {
f2blocks = get_f2blocks()
g = global$graph
withProgress(message = 'Finding graphs...', {
out = g %>% graph_addleaf(input$addleaf)
})
withProgress(message = paste('Evaluating', nrow(out), 'graphs...'), {
qpgfun = qpgraphfun()
out %<>% mutate(res = furrr::future_map(graph, ~qpgfun(f2blocks, .))) %>%
unnest_wider(res) %>%
select(starts_with('from'), starts_with('to'), starts_with('score'), starts_with('graph'), starts_with('edges'))
})
if(nrow(out) > 0) out %<>% arrange(score)
out
})
qpgraph_scorerange = function(out, low = 0, high = 1) {
out %>%
select(score) %>%
filter(!is.na(score)) %>%
arrange(score) %>%
slice(ceiling((low+1e-10)*n()), round(n()/2), ceiling(n()*high)) %>%
add_column(type = c('lo', 'mid', 'hi'), .before = 1) %>%
deframe
}
qpgraph_ranges = function(out, low = 0, high = 1) {
out %>%
select(from, to, weight) %>%
filter(!is.na(weight)) %>%
group_by(from, to) %>%
arrange(weight) %>%
slice(ceiling((low+1e-10)*n()), round(n()/2), ceiling(n()*high)) %>%
mutate(type = c('lo', 'mid', 'hi')) %>%
ungroup %>%
select(from, to, type, weight) %>%
pivot_wider(names_from = type, values_from = weight)
}
add_graph_to_history = function(global) {
print('add_graph_to_history')
isolate({
new = list(graph = global$graph,
edges = global$edges,
nadmix = numadmix(global$graph),
score = global$score,
worst_residual = if(is.null(global$worst_residual)) NA else global$worst_residual)
global$history %<>% c(list(new))
})
}
get_initresample0 = reactive({
print('get initresample')
input$init_resample
f2blocks = get_f2blocks()
g = global$graph
numstart = as.numeric(isolate(input$numstart2))
withProgress(message = paste('Evaluating graph ', numstart, ' times...'), {
qpgfun = qpgraphfun()
out = tibble(id = seq_len(numstart)) %>%
mutate(res = furrr::future_map(id, ~qpgfun(f2blocks, g, numstart = 1))) %>%
unnest_wider(res) %>%
select(id, starts_with('score'), starts_with('opt'), starts_with('graph'), starts_with('edges'))
})
})
get_initresample = reactive({
out = get_initresample0()
if(nrow(out) > 0 && 'edges' %in% names(out)) {
global$qpgraph_scorerange = qpgraph_scorerange(out, low = 0.025, high = 0.975)
out %<>% unnest(edges)
global$qpgraph_ranges = qpgraph_ranges(out, low = 0.025, high = 0.975)
out %<>% mutate(edge = paste(from, to, sep = ' -> ')) %>%
select(id, opt, score, edge, weight) %>%
pivot_wider(names_from = edge, values_from = weight) %>%
unnest(opt)
}
print('initresample done')
out
})
get_indresample0 = reactive({
print('get indresample')
input$ind_resample
poplist = global$poplist
countdir = global$countdir
inds = unlist(poplist)
numsingle = sum(sapply(poplist, length) == 1)
pops = rep(names(poplist), sapply(poplist, length))
f2blocks = get_f2blocks()
g = global$graph
qpgfun = qpgraphfun()
fun = admixtools:::make_resample_inds_fun(~qpgfun(..., y = g))
withProgress(message = paste('Evaluating graph ', length(inds)-numsingle, ' times...'), {
out = fun(countdir, inds = inds, pops = pops, numstart = 1)
})
out
})
get_indresample = reactive({
out = get_indresample0()
if(nrow(out) > 0) {
print('head(out)')
print(head(out))
print(head(out$error))
global$qpgraph_scorerange = qpgraph_scorerange(out, low = 0.025, high = 0.975)
out %<>% unnest(edges)
global$qpgraph_ranges = qpgraph_ranges(out, low = 0.025, high = 0.975)
out %<>%
mutate(edge = paste(from, to, sep = ' -> ')) %>%
select(ind, pop, edge, score, weight) %>%
pivot_wider(names_from = edge, values_from = weight)
}
print('indresample done')
out
})
get_snpresample0 = reactive({
f2blocks = get_f2blocks()
input$snp_resample
g = global$graph
bootnum = ifelse(input$boot, input$bootnum, FALSE)
withProgress(message = paste('Evaluating graph ', ifelse(bootnum, bootnum, dim(f2blocks)[[3]]),' times...'), {
out = admixtools:::make_resample_snps_fun(~qpgraphfun()(..., y = g))(f2blocks, boot = bootnum, numstart = 1)
})
})
get_snpresample = reactive({
out = get_snpresample0()
if(nrow(out) > 0 && 'edges' %in% names(out)) {
print('head(out)')
print(head(out))
global$qpgraph_scorerange = qpgraph_scorerange(out, low = 0.025, high = 0.975)
out %<>% unnest(edges)
global$qpgraph_ranges = qpgraph_ranges(out, low = 0.025, high = 0.975)
out %<>%
mutate(edge = paste(from, to, sep = ' -> ')) %>%
select(id, edge, score, weight) %>%
pivot_wider(names_from = edge, values_from = weight)
}
print('snpres done')
out
})
get_qpadm_rightpanel = reactive({
print('get qpadm_rightpanel')
global$qpadm_rightpanel
})
output$dashboardbody = renderUI({
print('dashboardbody')
get_bod()
})
output$popselectors = renderUI({
poplist = global$poplist
indnames = unique(global$allinds)
print('popselectors')
print(indnames)
sellist = imap(names(poplist), ~{
#buttlab = ifelse(is.null(global$poplistbak[[.x]]), 'group', 'ungroup')
isgrouped = length(poplist[[.x]]) == 1 && admixtools:::is_group(global$countdir, poplist[[.x]])
buttlab = ifelse(isgrouped, 'ungroup', 'group')
fillCol(tagList(textInput(paste0('pop', .y), NULL, .x, placeholder = 'pop name'),
shiny::actionButton(paste0('group', .y), buttlab),
selectizeInput(.x, NULL, indnames,
poplist[[.x]],
multiple = TRUE, width = 'auto'),
tags$style(type="text/css", paste0('.selectize-control{white-space:pre-wrap}'))), height = '1000px')})
print('popselectors2')
do.call(splitLayout, sellist)
})
output$econstraints = renderUI({
global$useconstraints = FALSE
print('econstraints')
edgecols = isolate(get_edgecols())
edges = isolate(get_fit()$edges) %>%
mutate(nam = paste(from, to, sep = ' -> '), nam2 = str_remove_all(nam, '\\.|>| ')) %>%
left_join(edgecols, by = c('nam2'='edge')) %>%
filter(type == 'edge') %>% select(nam, nam2, weight, col)
maxlen = ceiling(max(edges$weight))*edgemul
print('econstraints 2')
div(div(id = 'econstraints_update_div', actionButton('econstraints_update', 'Update')),
sliderInput('any_econstraints', 'All edges', 0, maxlen, c(0, maxlen), dragRange = FALSE),
edges %>%
pmap(~{div(tags$style(HTML(paste0('[for=',..2,']+span>.irs-bar, [for=',..2,']+span>.irs>.irs-from, [for=',..2,']+span>.irs>.irs-to {background: ', ..4,'}'))),
sliderInput(..2, ..1, min = 0, max = maxlen, value = c(0, maxlen),
dragRange = FALSE, width = '100%'))}))
})
output$aconstraints = renderUI({
global$useconstraints = FALSE
print('aconstraints')
input$add_edge
input$delete_edge
edgecols = isolate(get_edgecols())
edges = isolate(get_fit()$edges) %>%
mutate(nam = paste(from, to, sep = ' -> '), nam2 = str_remove_all(nam, '\\.|>| ')) %>%
left_join(edgecols, by = c('nam2'='edge')) %>%
filter(type == 'admix') %>% select(nam, nam2, weight, col)
div(div(id = 'aconstraints_update_div', actionButton('aconstraints_update', 'Update')),
edges %>%
pmap(~{div(tags$style(HTML(paste0('[for=',..2,']+span>.irs-bar, [for=',..2,']+span>.irs>.irs-from, [for=',..2,']+span>.irs>.irs-to {background: ', ..4,'}'))),
sliderInput(..2, ..1, min = 0, max = 1, value = c(0, 1),
dragRange = FALSE, width = '100%'))})
)
})
output$nadmix = renderUI({
nad = 3
if(!is.null(global$graph)) nad = numadmix(global$graph)
div(sliderInput('nadmix', '# admix', 0, 10, nad),
selectizeInput('outpop', 'Outgroup', c('< undefined >', names(global$poplist)), selected = names(global$poplist)[1]))
})
output$qpgraph_add = renderUI({
print('selectInput addleaf')
choices = setdiff(names(global$poplist), get_leafnames(global$graph))
selectInput('addleaf', '', choices = choices)
})
output$showseledge = renderText({
input$clear_edge
print('showseledge1')
se = global$seledge
print('se')
print(global$seledge)
print(se)
if(!is.null(se) && (is.na(se) || se == 'NA')) {
print('is na')
se = NULL
}
print('se2')
print(se)
se
})
observeEvent(input$any_econstraints, {
print('update triggered by any_econstraints')
req(global$graph)
print('any_econstraints 2')
edges = global$graph %>% igraph::as_edgelist() %>% as_tibble(.name_repair = ~c('from', 'to')) %>%
group_by(to) %>% mutate(type = ifelse(n() > 1, 'admix', 'normal')) %>% ungroup %>%
mutate(eid = paste(from, to, sep = '->') %>% str_remove_all('\\.|>| ')) %$% eid
print('any_econstraints 3')
if(is.null(isolate(input[[edges[1]]]))) return()
edges %>% walk(~{updateSliderInput(session, ., value = c(input$any_econstraints[1],
input$any_econstraints[2]))})
print('any_econstraints 4')
})
observeEvent(input$delete_edge, {
print('del 1')
eg = get_seledge()
if(is.null(eg)) return()
gnew = global$graph
leaves = get_leafnames(gnew)
print('del 2')
print(eg)
eg %<>% str_split('\n') %>% `[[`(1)
if(str_detect(eg[1], ' -> ')) {
eg %<>% str_split(' -> ')
print(eg)
for(i in 1:length(eg)) {
#browser()
gnew = delete_admix(gnew, eg[[i]][1], eg[[i]][2])
# g <<- g
# gnew <<- gnew
if(!admixtools:::is_valid(gnew)) {
shinyalert('Error', 'Could not delete edge!')
return()
}
}
} else {
for(i in 1:length(eg)) {
gnew = delete_leaf(gnew, eg[i])
if(!admixtools:::is_valid(gnew)) {
browser()
shinyalert('Error', 'Could not delete node!')
return()
}
}
}
newleaves = setdiff(get_leafnames(gnew), leaves)
if(length(newleaves) > 0) {
# newleaves <<- newleaves
# oldgraph <<- global$graph
# newgraph <<- g
shinyalert('New leaves!', newleaves)
}
global$graph = gnew
global$seledge = NULL
global$edges = gnew %>% as_edgelist %>% as_tibble(.name_repair = ~c('from', 'to'))
global$score = NA
global$worst_residual = NA
shinyjs::disable('clear_edge')
shinyjs::disable('add_edge')
shinyjs::disable('delete_edge')
print('delete_edge done')
})
observeEvent(input$clear_edge, {
print('clear')
global$seledge = NULL
shinyjs::disable('clear_edge')
shinyjs::disable('add_edge')
shinyjs::disable('delete_edge')
})
observeEvent(input$add_edge, {
print('add')
eg = get_seledge() %>% str_split('\n') %>% `[[`(1) %>% str_split(' -> ')
g = global$graph
leaves = get_leafnames(g)
if(length(eg) == 1) {
print('input$addleaf')
print(input$addleaf)
if(is.null(input$addleaf)) {
shinyalert('Error', 'Can only add edge which connects exactly two other edges!')
return()
}
from = eg[[1]][1]
to = eg[[1]][2]
gnew = admixtools:::insert_leaf(g, input$addleaf, from, to)
if(!admixtools:::is_valid(gnew)) {
shinyalert('Error', 'Could not add population!')
return()
}
global$graph = gnew
global$seledge = NULL
global$edges = gnew %>% as_edgelist %>% as_tibble(.name_repair = ~c('from', 'to'))
global$score = NA
global$worst_residual = NA
shinyjs::disable('clear_edge')
shinyjs::disable('add_edge')
shinyjs::disable('delete_edge')
} else if(length(eg) == 2) {
#from = eg[[1]][2]
#to = eg[[2]][2]
#print(paste(from, to))
alert = function(x) shinyalert('Could not insert edge!', as.character(x))
tryCatch({
#gnew = admixtools:::insert_admix_old(g, from, to, allow_below_admix = TRUE)
gnew = admixtools:::insert_admix(g, eg[[1]][1], eg[[1]][2], eg[[2]][1], eg[[2]][2])
}, warning = alert, error = alert)
if(!exists('gnew') || is.null(gnew)) return()
if(!admixtools:::is_valid(gnew)) {
shinyalert('Error', 'Could not insert edge!')
return()
}
global$graph = gnew
global$seledge = NULL
global$edges = gnew %>% as_edgelist %>% as_tibble(.name_repair = ~c('from', 'to'))
global$score = NA
global$worst_residual = NA
shinyjs::disable('clear_edge')
shinyjs::disable('add_edge')
shinyjs::disable('delete_edge')
newleaves = setdiff(get_leafnames(gnew), leaves)
if(length(newleaves) > 0) shinyalert('New leaves!', newleaves)
} else {
shinyalert('Error', 'Too many items selected!')
return()
}
})
observeEvent(input$run_qpadm, {
print('observe run_qpadm')
global$qpadmout = get_qpadmout()
global$qpadm_rightpanel = plotOutput('qpadmcomparison')
})
observeEvent(input$run_qpgraph, {
print('observe run_qpgraph')
global$qpgraphout = get_qpgraphout()
#global$qpg_right = plotlyOutput('graphcomparison')
})
map(1:4,
~observeEvent(input[[paste0('f4pops', .)]], {
if(input$f4_allowdup) return()
popsinuse = unique(c(input$f4pops1, input$f4pops2, input$f4pops3, input$f4pops4))
map(setdiff(1:4, .), ~{
fp = paste0('f4pops', .)
sel = input[[fp]]
updateSelectizeInput(session, fp, selected = sel,
choices = union(sel, setdiff(names(global$poplist), popsinuse)))})
}))
observeEvent(input$f4_allowdup, {
map(1:4, ~{
fp = paste0('f4pops', .)
sel = input[[fp]]
updateSelectizeInput(session, fp, selected = sel,
choices = names(global$poplist))})
})
map(1:3, ~observeEvent(input[[paste0('qpadmpops', .)]], {
popsinuse = unique(c(input$qpadmpops1, input$qpadmpops2, input$qpadmpops3))
map(setdiff(1:3, .), ~{
fp = paste0('qpadmpops', .)
sel = input[[fp]]
global[[fp]] = sel
choices = union(sel, setdiff(names(global$poplist), popsinuse))
updateSelectizeInput(session, fp, selected = sel, choices = choices)})
}))
observeEvent(input$qpadm_randomize, {
nam = sample(names(global$poplist))
nl = floor(length(nam)/2)
choices = list(nam[1], nam[2:nl], nam[(nl+1):length(nam)])
map(1:3, ~{
global[[paste0('qpadmpops', .)]] = choices[[.]]
updateSelectizeInput(session, paste0('qpadmpops', .),
selected = choices[[.]], choices = nam)})
})
observeEvent(input$qpadm_rotate, {
#global$show_qpadm_rotate = input$qpadm_rotate
global$qpadmpops1 = input$qpadmpops1
global$qpadmpops2 = input$qpadmpops2
global$qpadmpops3 = input$qpadmpops3
})
observeEvent(input$qpadm_results, {
global$qpadmpops1 = input$qpadmpops1
global$qpadmpops2 = input$qpadmpops2
global$qpadmpops3 = input$qpadmpops3
#global$show_qpadm_results = input$qpadm_results
})
get_qpadm_rotate = reactive({
f2blocks = get_f2blocks()
p1 = input$qpadmpops1
p2 = input$qpadmpops2
p3 = input$qpadmpops3
req(f2blocks, p1, p2, p3)
withProgress(message = paste0('Generating qpadm models...'), {
lr = admixtools:::all_lr2(p1, length(p2))
})
withProgress(message = paste0('Evaluating ', length(lr[[1]]), ' qpadm models...'), {
out = admixtools:::qpadm_eval_rotate(f2blocks, p3, lr, p2, verbose = F) %>%
select(chisq, p, feasible, left, right)
})
out
})
get_qpadmout = reactive({
print('get_qpadmout')
get_qpadm()
volumes %<>% eval
pref = parseFilePaths(volumes, input$qpadm_genofile)$datapath
pref %<>% str_replace('\\.geno$', '') %>% normalizePath(mustWork = FALSE)
bin = parseFilePaths(volumes, input$qpadmbin)$datapath %>% normalizePath(mustWork = FALSE)
#if(length(bin) == 0) bin = '~/Downloads/admixtoools_hub_reldec19/bin/qpAdm'
#if(length(pref) == 0) pref = '~/Documents/v42.1_afafam'
target = input$qpadmpops1
left = input$qpadmpops2
right = input$qpadmpops3
withProgress(message = 'running original qpAdm...', {
dir = tempdir()
out = tryCatch(qpadm_wrapper(left, right, target, bin = bin, pref = pref, outdir = dir),
error = function(e) shinyalert('Error!', as.character(e)))
})
out
})
get_qpgraphout = reactive({
print('get_qpgraphout')
get_fit()
volumes %<>% eval
pref = parseFilePaths(eval(volumes), input$qpgraph_genofile)$datapath
pref %<>% str_replace('\\.geno$', '') %>% normalizePath(mustWork = FALSE)
bin = parseFilePaths(volumes, input$qpgraphbin)$datapath %>% normalizePath(mustWork = FALSE)
print('bin')
print(bin)
# if(length(bin) == 0) bin = '~/Downloads/admixtoools_hub_reldec19/bin/qpGraph'
# if(length(pref) == 0) pref = '~/Documents/v42.1_afafam'
g = global$graph
withProgress(message = 'running original qpGraph...', {
dir = tempdir() %>% normalizePath()
print('pref')
print(pref)
print('dir')
print(dir)
out = qpgraph_wrapper(g, bin, pref, outdir = dir, zthresh = 0)
out$f2 = NULL
print('parse fstats')
fstatsfile = paste0(dir, '/fstats.out')
if(file.exists(fstatsfile)) global$precomp = admixtools:::parse_fstats(fstatsfile)
})
out
})
observeEvent(input$econstraints_update, {global$useconstraints = TRUE})
observeEvent(input$aconstraints_update, {global$useconstraints = TRUE})
dto = function(x) div(DT::DTOutput(x), style = dtstyle)
#dtof4 = reactive({ if(input$return_fstats) dto('f4') else renderText("Activate 'return_fstats' under 'Fit'!") })
observeEvent(input$qpadm_fit, {global$qpadm_rightpanel = qpadm_rightpanel_fit()})
# observeEvent(input$qpgraph_fit, {
# print('fit!')
# global$qpgraph_ranges = NULL
# global$qpgraph_scorerange = NULL
# global$qpg_right = qpg_right_fit()
# })
qpadm_leftpanel = reactive({
print('reactive qpadm_leftpanel')
all = names(global$poplist)
choices1 = setdiff(all, c(global$qpadmpops2, global$qpadmpops3))
choices2 = setdiff(all, c(global$qpadmpops1, global$qpadmpops3))
choices3 = setdiff(all, c(global$qpadmpops1, global$qpadmpops2))
div(
selectizeInput('qpadmpops1', 'Left', choices = choices1, multiple = TRUE, selected = global$qpadmpops1),
selectizeInput('qpadmpops2', 'Right', choices = choices2, multiple = TRUE, selected = global$qpadmpops2),
selectizeInput('qpadmpops3', 'Target', choices = choices3, multiple = FALSE, selected = global$qpadmpops3),
splitLayout(actionButton('qpadm_randomize', 'Randomize')))
})
qpadm_rightpanel_fit = reactive({
tabsetPanel(id = 'qpadm_tabset',
tabPanel('Weights', dto('qpadm_weights')),
tabPanel('f4', dto('qpadm_f4')),
tabPanel('Rank drop', dto('qpadm_rankdrop')),
tabPanel('Pop drop', dto('qpadm_popdrop')))
})
qpadm_rightpanel_rotate = reactive(div(dto('qpadm_rot'), style='padding-right:50px;'))
qpadm_rightpanel = reactive({
#res = global$show_qpadm_results
#rot = global$show_qpadm_rotate
res = if(is.null(input$qpadm_results)) TRUE else input$qpadm_results
rot = if(is.null(input$qpadm_rotate)) FALSE else input$qpadm_rotate
if(rot && res) splitLayout(cellWidths = c('50%', '50%'), qpadm_rightpanel_rotate(), qpadm_rightpanel_fit())
else if(rot) qpadm_rightpanel_rotate()
else if(res) qpadm_rightpanel_fit()
else div()
})
get_loaddata = reactive({
print('load_data')
input$textdirinput
bh = '140px'
div(
fluidRow(
box(width=4, height=bh, background = cols[1], h4('Select data directory'),
div(splitLayout(textInput('textdirinput', NULL, placeholder = 'Data directory'),
actionButton('textdir', 'Create')))),
conditionalPanel('output.show_extract == "1"',
fluidRow(box(width=4, height=bh, background = cols[2], h4('Extract data'),
splitLayout(
div(
textInput('textgenofile1', NULL, placeholder = 'genotype file')
),
actionButton('extract_data', 'Extract data'))))),
conditionalPanel('output.show_indselect == "1"', (
box(width=4, height=bh, background = cols[3], h4('Select .ind or .fam file'),
div(fileInput('popfile', NULL, placeholder = '', buttonLabel = 'Ind file'), id = 'popfilediv')))),
fluidRow(column(12, conditionalPanel('input.extract_data > 0', verbatimTextOutput('console')))),
#box(width=6, textOutput('console')),
div(style = 'visibility: hidden', verbatimTextOutput('show_popadjust')),
div(style = 'visibility: hidden', verbatimTextOutput('show_extract')),
div(style = 'visibility: hidden', verbatimTextOutput('show_indselect')),
conditionalPanel('output.show_popadjust == "1"', fluidRow(column(12,
box(width=12, background = cols[5],
fluidRow(column(2, h4('Adjust populations', id = 'adjpopdiv')),
column(1, div(actionButton('removepop', '–'),
actionButton('addpop', '+')))),
fluidRow(column(11, uiOutput('popselectors')))))))))
})
observeEvent(input$extract_data, {
if(is.null(input$textgenofile1) || !file.exists(input$textgenofile1)) {
shinyalert('Error!', 'Please provide path and name for a .geno or .bed file!')
return()
}
pref = input$textgenofile1 %>% str_replace('\\.geno$|\\.bed$', '')
oldnam = names(global$poplist)
nam = map(paste0('pop', seq_len(length(global$poplist))), ~input[[.]])
if(is.null(nam[[1]])) return()
global$poplist = map(oldnam, ~input[[.]]) %>% purrr::set_names(nam)
poplist = global$poplist
inds = unlist(poplist)
pops = rep(names(poplist), sapply(poplist, length))
transitions = input$trans_extract %in% c('both', 'only transitions')
transversions = input$trans_extract %in% c('both', 'only transversions')
if(input$fix_populations) extractfun = function(...) extract_f2(pops = pops, ...)
else extractfun = extract_counts
print('inds')
print(inds)
print('global$countdir')
print(global$countdir)
withCallingHandlers({
shinyjs::html('console', '')
print(list(pref, global$countdir, inds = inds,
maxmiss = input$max_miss, minmaf = input$min_maf, maxmaf = input$max_maf,
transitions = transitions, transversions = transversions,
keepsnps = input$keepsnps, maxmem = input$maxmem*1e3))
extractfun(pref, global$countdir, inds = inds,
maxmiss = input$max_miss, minmaf = input$min_maf, maxmaf = input$max_maf,
transitions = transitions, transversions = transversions,
keepsnps = input$keepsnps, maxmem = input$maxmem*1e3)
},
message = function(m) {
shinyjs::html(id = 'console', html = m$message, add = TRUE)
})
global$iscountdata = 'indivs' %in% list.dirs(global$countdir,F,F)
global$isf2data = !'indivs' %in% list.dirs(global$countdir,F,F) &
'block_lengths_f2.rds' %in% list.files(global$countdir)
if(global$iscountdata) {
global$allinds = list.dirs(paste0(global$countdir, '/pairs'),F,F)
} else if(global$isf2data) {
global$allinds = list.dirs(global$countdir,F,F)
global$poplist = global$allinds %>% purrr::set_names(global$allinds)
}
})
observeEvent(event_data('plotly_click', source = 'src'), {
print('observe seledge')
ed = event_data('plotly_click', source = 'src')
isolate({
plt = plotly_graph()
newe = ed_to_name(ed, plt)
oldedge = global$seledge
if(!is.null(oldedge) && is.na(oldedge)) {
print('seledge is na. why?')
oldedge = NULL
}
global$seledge = newe
shinyjs::enable('clear_edge')
shinyjs::enable('add_edge')
shinyjs::enable('delete_edge')
})
if(!is.null(newe) && !is.null(oldedge) && newe != oldedge) global$seledge = paste(c(oldedge, newe), collapse = '\n')
})
get_seledge = reactive({
print('get_seledge')
print(global$seledge)
global$seledge
})
ed_to_name = function(ed, plt) {
dat = plt$x$data[[ed$curveNumber + 1]]
type = ifelse(length(dat) == 14, 'text', ifelse('hovertext' %in% names(dat), 'edge', 'node'))
pt = ed$pointNumber + 1
if (type == 'edge') {
name = dat$hovertext[pt]
} else if (type == 'node') {
col = ifelse(duplicated(dat$text)[pt], 5, 3)
name = dat$text[pt]
} else {
name = dat$text[pt]
}
print('ed_to_name22')
name
}
get_edgecols = reactive({
plt = plotly_graph()
edgecols = plt$x$data %>% map(~pluck(., 'line', 'color')) %>% discard(is.null)
txt = plt$x$data[seq_len(length(edgecols))] %>% map(~pluck(., 'text'))
tibble(edge = str_remove_all(unlist(txt), '\\.|>| '),
col = rep(as.vector(edgecols), map_dbl(txt, length))) %>%
unnest(col) %>% filter(!duplicated(.), !is.na(edge))
})
plotly_graph = reactive({
print('plotly_graph')
poplist = global$poplist
global$graph
pdat = get_pdat()
req(pdat, poplist)
nadmix = numadmix(graph_from_edgelist(as.matrix(global$edges)[,1:2]))
textsize = 2.5
print('plotly_graph 2')
eg = pdat$eg
nodelabs = imap(poplist, ~paste(.y, paste(.x, collapse = '\n'), sep = '\n')) %>%
unlist %>% enframe(name = 'name', value = 'text')
#nodelabs = poplist %>% group_by(pop) %>%
# summarize(text = paste(pop[1], paste(ind, collapse='\n'), sep = '\n'))
nodes = pdat$nodes %>% left_join(nodelabs, by = c('name'))
allnodes = eg %>% transmute(from, x, y, xend=0, yend=0, to=0) %>% filter(!duplicated(.))
segtext = "ifelse(indegree == 1, to, from)"
segtext = "paste(from, to, sep = ' -> ')"
print('plotly_graph 3')
sr = global$qpgraph_scorerange
scoretext = paste0('score: ', ifelse(is.null(sr), round(global$score, 2),
paste0(round(sr['mid'], 2), '\n[', round(sr['lo'], 0), ' - ', round(sr['hi'], 0), ']')),
ifelse(is.na(global$worst_residual), '', paste0('\nWR: ', round(global$worst_residual, 2))),
'\nadmix: ', nadmix,
ifelse(global$tempgraph, '\n\nPREVIEW!', ''))
gg = eg %>% mutate(rownum = 1:n()) %>%
ggplot(aes(x=x, xend=xend, y=y, yend=yend, from=from, to=to)) +
geom_segment(aes_string(linetype = 'type', col = 'as.factor(y)', text = segtext),
arrow=arrow(type = 'closed', angle = 10, length=unit(0.15, 'inches'))) +
geom_point(data = allnodes, aes(x, y, text = from), col = 'black', alpha = 0) +
annotate('text', x = min(nodes$x), y = max(nodes$yend), hjust = 0, vjust = 1,
label = scoretext) +
# annotate('text', x = min(nodes$x), y = max(nodes$yend), hjust = 0,
# label = scoretext) +
theme(panel.background = element_blank(),
axis.line = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank(),
legend.position = 'none') +
xlab('') + ylab('') +
scale_linetype_manual(values = c(admix=3, normal=1)) +
ggtitle('') +
scale_x_continuous(expand = c(0.1, 0.1))
if('highlight_unidentifiable' %in% input$plotopt) {
unid = unidentifiable_edges(global$graph)
unid2 = eg %>% right_join(unid %>% select(-type), by = c('from', 'to'))
gg = gg + geom_segment(aes_string(linetype = 'type'), col = 'red', size = 1, data = unid2)
}
gg = gg + geom_text(aes(x = (x+xend)/2, y = (y+yend)/2, label = label, text = paste(from, to, sep = ' -> ')),
size = textsize) +
geom_text(data = nodes, aes_string(label = 'name', col = 'as.factor(yend)',
from = NA, text = 'text'), size = textsize)
print('plotly_graph 4')
plt = plotly::ggplotly(gg, source = 'src', tooltip=c('text'))
print('plotly_graph 5')
plt
})
output$graphplot <- renderPlotly({
print('graphplot')
suppressWarnings({withProgress(message = 'Processing...', {plt = plotly_graph()})})
#event_register(plt, 'plotly_hover')
print('graphplot2')
plt
})
get_f2 = reactive({
print('get f2')
f2blocks = if(input$fst == 'fst') get_fstblocks() else get_f2blocks()
#req(f2blocks, p1, p2, p3, p4)
f2(f2blocks)
})
get_f3 = reactive({
print('get f3')
f2blocks = get_f2blocks()
f3(f2blocks)
})
get_f4 = reactive({
print('get f4')
input$f4_run
f2blocks = get_f2blocks()
p1 = input$f4pops1
p2 = input$f4pops2
p3 = input$f4pops3
p4 = input$f4pops4
req(f2blocks, p1, p2, p3, p4)
f4(f2blocks, p1, p2, p3, p4)
})
get_qpadm = reactive({
print('get qpadm')
input$qpadm_run
f2blocks = get_f2blocks()
p1 = input$qpadmpops1
p2 = input$qpadmpops2
p3 = input$qpadmpops3
req(f2blocks, p1, p2, p3)
#if(global$show_qpadm_rotate) {
if(input$qpadm_rotate) {
p1 = if(is.null(global$qpadmpops_left_remain)) p1 else global$qpadmpops_left_remain
p2 = union(p2, global$qpadmpops_right_add)
}
qpadm(f2blocks, p1, p2, p3)
})
plotly_f2 = reactive({
print('get f2')
f2blocks = if(input$fst == 'fst') get_fstblocks() else get_f2blocks()
pdat = f2blocks %>%
apply(1:2, mean, na.rm = T) %>%
as_tibble(rownames = 'pop1') %>%
pivot_longer(-pop1, 'pop2', values_to = 'est')
gg = pdat %>%
ggplot(aes(pop1, pop2, text = paste(pop1, pop2, round(est, 3), sep='\n'))) +
geom_tile(aes(fill = est)) +
xlab('') +
ylab('') +
theme(axis.text.x = element_text(angle = 90, hjust = 1), panel.background = element_blank()) +
scale_fill_distiller(palette = 'Purples')
plt = plotly::ggplotly(gg, source = 'src_f2', tooltip=c('text'))
plt
})
plotly_f2_indiv = reactive({
f2dat = get_f2dat_indiv()
poplist = global$poplist
ind = enframe(poplist, 'pop', 'ind') %>% unnest(ind)
dat = f2dat %>%
transmute(ind1 = pop1, ind2 = pop2, f2corr = f2, f2uncorr) %>%
left_join(ind %>% transmute(ind1 = ind, pop1 = pop), by = 'ind1') %>%
left_join(ind %>% transmute(ind2 = ind, pop2 = pop), by = 'ind2') %>%
arrange(pop1, pop2) %>%
mutate(i1 = as.numeric(factor(ind1, levels = unique(ind1[order(pop1)]))),
i2 = as.numeric(factor(ind2, levels = unique(ind2[order(pop2)]))),
f2 = ifelse(i1 > i2, f2corr, f2uncorr))
d2 = enframe(sapply(poplist, length)) %>%
arrange(name) %>%
mutate(end = cumsum(value), start = lag(end, default = 0), mn = start + value/2 + 0.5)
exp = 0
gg = dat %>%
ggplot(aes(i1, i2)) +
geom_tile(aes(fill = f2, text = paste(ind1, ind2, round(f2, 3), sep = '\n'))) +
geom_vline(aes(xintercept = end+0.5), data = d2 %>% head(-1) %>% rename(ind2 = name)) +
geom_hline(aes(yintercept = end+0.5), data = d2 %>% head(-1) %>% rename(ind2 = name)) +
scale_x_continuous(labels = d2$name, breaks = d2$mn, expand = c(exp, exp)) +
scale_y_continuous(labels = d2$name, breaks = d2$mn, expand = c(exp, exp)) +
xlab('') + ylab('') +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
panel.background = element_blank()) +
scale_fill_distiller(palette = 'Purples')
plt = plotly::ggplotly(gg, source = 'src_f2', tooltip = 'text')
plt
})
plotly_f4 = reactive({
print('get f4')
input$f4_run
#isolate({
f2blocks = get_f2blocks()
p1 = input$f4pops1[1]
p2 = input$f4pops2[1:2]
# f4pops2 needs at least two populations. how to check for this? otherwise, "Error: Only strings can be converted to symbols"
f4dat = get_f4() %>% filter(pop1 == p1, pop2 %in% p2)
req(f4dat)
#})
f4dat %<>%
mutate(pop2 = recode(pop2, !!sym(p2[1]):='x', !!sym(p2[2]):='y')) %>%
pivot_wider(names_from = pop2, values_from = c(est, se, z, p)) %>%
mutate(mul = input$f4_se)
xmin = min(f4dat$est_x - f4dat$se_x)
xmax = max(f4dat$est_x + f4dat$se_x)
ymin = min(f4dat$est_y - f4dat$se_y)
ymax = max(f4dat$est_y + f4dat$se_y)
gg = f4dat %>%
ggplot(aes(est_x, est_y, col = pop3, text = paste0('f4(', p1, ', X; ', pop3, ', ', pop4, ')'))) +
geom_smooth(method = 'lm', se = FALSE, formula=y~x-1) +
geom_hline(yintercept = 0) +
geom_vline(xintercept = 0) +
geom_point() +
theme(panel.background = element_blank()) +
xlab(paste0('f4(', p1, ', ', p2[1], '; Pop 3, Pop 4)' )) +
ylab(paste0('f4(', p1, ', ', p2[2], '; Pop 3, Pop 4)' )) +
scale_x_continuous(limits = c(xmin, xmax)) +
scale_y_continuous(limits = c(ymin, ymax))
if(input$f4_se > 0) {
gg = gg +
geom_errorbar(aes(ymin = est_y - se_y*mul, ymax = est_y + se_y*mul), width = 0) +
geom_errorbarh(aes(xmin = est_x - se_x*mul, xmax = est_x + se_x*mul), height = 0)
}
plt = plotly::ggplotly(gg, source = 'src_f4', tooltip=c('text'))
})
plotly_f42 = reactive({
print('get f4')
input$f4_run
#isolate({
f2blocks = get_f2blocks()
f4dat = get_f4() %>% filter(pop1 == input$f4pops1[1], pop3 == input$f4pops3[1], pop4 == input$f4pops4[1])
req(f4dat)
#})
# f4dat %<>%
# mutate(pop2 = recode(pop2, !!sym(p2[1]):='x', !!sym(p2[2]):='y')) %>%
# pivot_wider(names_from = pop2, values_from = c(est, se, z, p))
# xmin = min(f4dat$est_x - f4dat$se_x)
# xmax = max(f4dat$est_x + f4dat$se_x)
# ymin = min(f4dat$est_y - f4dat$se_y)
# ymax = max(f4dat$est_y + f4dat$se_y)
gg = f4dat %>%
ggplot(aes(pop2, est, col = pop2)) +
geom_hline(yintercept = 0) +
geom_point() +
theme(panel.background = element_blank()) +
xlab(paste0('' )) +
ylab(paste0('f4(', input$f4pops1[1], ', X; ', input$f4pops3[1], ', ', input$f4pops4[1],')' ))
if(input$f4_se > 0) {
mul = input$f4_se
gg = gg + geom_errorbar(aes(ymin = est - se*mul, ymax = est + se*mul), width = 0)
}
plt = plotly::ggplotly(gg, source = 'src_f42')
})
output$f2heatmap = renderPlotly({
print('f2plot')
plt = plotly_f2()
plt
})
output$f2heatmap_indiv = renderPlotly({
print('f2plot indiv')
plt = plotly_f2_indiv()
plt
})
output$f4plot = renderPlotly({
print('f4plot')
plt = plotly_f4()
plt
})
output$f4plot2 = renderPlotly({
print('f4plot')
plt = plotly_f42()
plt
})
output$qpadmcomparison = renderPlot({
print('qpadmcomparison')
withProgress(message = 'running qpadm in R...', {
qpadm_R = get_qpadm()
})
original = global$qpadmout
gg = plot_comparison(original, qpadm_R)
#plotly::ggplotly(gg, tooltip = 'label')
gg
})
output$graphcomparison = renderPlotly({
print('graphcomparison')
req(global$qpgraphout)
withProgress(message = 'running qpgraph in R...', {
qpgraph_R = get_fit()
})
#qpgraph_R$f3 = qpgraph_R$f3 %>% select(pop2, pop3, fit, est) # temporary; line should be removed
gg = plot_comparison(global$qpgraphout, qpgraph_R)
plotly::ggplotly(gg, tooltip = 'label')
})
output$downloadgraph = downloadHandler('graph.tsv', function(file) {
edges = get_fit()$edges
if(input$downloadgraph_format == 'admixtools') edges %>% admixtools:::fit_to_qpgraph_format() %>% write(file)
else if(input$downloadgraph_format == 'edgelist') edges %>% write_tsv(file)
else if(input$downloadgraph_format == 'dot') edges %>% admixtools:::write_dot(file)
else stop()
})
ex = c('Buttons', 'Scroller')
dtfun = function(...) DT::renderDataTable(..., extensions = ex, options = do, selection = 'single')
output$f2 = dtfun({format_table(get_fit()$f2)})
output$f3 = dtfun({format_table(get_fit()$f3)})
output$f4 = dtfun({if(is.null(get_fit()$f4)) as.data.frame("Activate 'return_fstats' under 'Fit'!") %>% slice(-1) else format_table(get_fit()$f4)})
output$opt = dtfun({format_table(get_fit()$opt)})
output$history = dtfun({format_table(get_history() %>% select(-graph, -edges))})
output$minus1 = dtfun({format_table(get_minus1())})
output$plus1 = dtfun(format_table(get_plus1()))
output$flipadmix = dtfun({format_table(get_flipadmix())})
output$addleaf = dtfun({format_table(get_addleaf())})
output$initresample = dtfun({format_table(get_initresample())})
output$indresample = dtfun({format_table(get_indresample())})
output$snpresample = dtfun({format_table(get_snpresample())})
output$f2stats = dtfun({format_table(get_f2())})
output$f3stats = dtfun({format_table(get_f3())})
output$f4stats = dtfun({format_table(get_f4())})
output$qpadm_weights = dtfun({format_table(get_qpadm()$weights)})
output$qpadm_f4 = dtfun({format_table(get_qpadm()$f4)})
output$qpadm_rankdrop = dtfun({format_table(get_qpadm()$rankdrop)})
output$qpadm_popdrop = dtfun({format_table(get_qpadm()$popdrop)})
output$qpadm_rot = dtfun({format_table(get_qpadm_rotate())})
onBookmark(function(state) {
names(global) %>% map(~{state$values[[.]] = global[[.]]})
})
onRestore(function(state) {
names(state$values) %>% map(~{global[[.]] = state$values[[.]]})
})
#setBookmarkExclude('sidebarItemExpanded')
#setBookmarkExclude('navbar')
observeEvent(input$tooltips, {
print('refresh tooltips')
if(input$tooltips %% 2 == 1) {
updateActionButton(session, 'tooltips', icon = icon('exclamation'))
imap(tt, ~removeTooltip(session, .y))
imap(tt, ~addTooltip(session, .y, .x))
} else {
updateActionButton(session, 'tooltips', icon = icon('question'))
imap(tt, ~removeTooltip(session, .y))
}
})
}
shinyApp(ui = ui, server = server, enableBookmarking = 'server')
#runApp(list(ui = ui, server = server), display.mode = 'showcase')
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.