tests.R
In vitkl/MItools: Protein interaction data tools: retrieving data from IntAct and most other databases using PSICQUIC service
library(roxygen2)
library(devtools)
document()
install()
library(PItools)
library(GenomicRanges)
library(R.utils)
data = fread("../viral_project/processed_data_files/viral_human_net_w_domains", sep = "\t", stringsAsFactors = F)
file_BioPlex3 = "../viral_project/processed_data_files/human_net_w_domains"
gunzip(paste0(file_BioPlex3,".gz"), remove = F)
data = fread(file_BioPlex3, sep = "\t", stringsAsFactors = F)
data = data[IDs_interactor_human_B != "UNKNOWN"]
unlink(file_BioPlex3)
res = permutationPval(interactions2permute = IDs_interactor_viral ~ IDs_interactor_human,
associations2test = IDs_interactor_viral ~ IDs_domain_human,
node_attr = list(IDs_interactor_viral ~ IDs_interactor_viral_degree,
IDs_domain_human ~ domain_count,
IDs_interactor_viral + IDs_domain_human ~ domain_frequency_per_IDs_interactor_viral),
data = data,
statistic = IDs_interactor_viral + IDs_domain_human ~ .N,
select_nodes = IDs_domain_human ~ domain_count >= 1,
N = 100, clustermq_jobs = 20, split_comp_inner_N = 2,
cores = NULL, seed = 2, clustermq = T, clustermq_mem = 2000)
plot(res)
res
library(lineprof)
library(PItools)
library(GenomicRanges)
library(R.utils)
data = fread("../viral_project/processed_data_files/viral_human_net_w_domains", sep = "\t", stringsAsFactors = F)
data_num = copy(data)
data_num[, IDs_interactor_human := as.integer(as.factor(IDs_interactor_human))]
data_num[, IDs_domain_human := as.integer(as.factor(IDs_domain_human))]
data_num[, IDs_interactor_viral := as.integer(as.factor(IDs_interactor_viral))]
l = microbenchmark::microbenchmark({
res <<- permutationPval(interactions2permute = IDs_interactor_viral ~ IDs_interactor_human,
associations2test = IDs_interactor_viral ~ IDs_domain_human,
node_attr = list(IDs_interactor_viral ~ IDs_interactor_viral_degree,
IDs_domain_human ~ domain_count,
IDs_interactor_viral + IDs_domain_human ~ domain_frequency_per_IDs_interactor_viral),
data = data,
statistic = IDs_interactor_viral + IDs_domain_human ~ .N,
select_nodes = IDs_domain_human ~ domain_count >= 1,
N = 10,
cores = 1, seed = 2, clustermq = F)
}, {
res_num <<- permutationPval(interactions2permute = IDs_interactor_viral ~ IDs_interactor_human,
associations2test = IDs_interactor_viral ~ IDs_domain_human,
node_attr = list(IDs_interactor_viral ~ IDs_interactor_viral_degree,
IDs_domain_human ~ domain_count,
IDs_interactor_viral + IDs_domain_human ~ domain_frequency_per_IDs_interactor_viral),
data = data_num,
statistic = IDs_interactor_viral + IDs_domain_human ~ .N,
select_nodes = IDs_domain_human ~ domain_count >= 1,
N = 10,
cores = 1, seed = 2, clustermq = F)
}, times = 12)
l
microbenchmark::microbenchmark({res <- permutationPval(interactions2permute = IDs_interactor_human_A ~ IDs_interactor_human_B, # first set of interacting pairs (XY) that are to be permuted
associations2test = IDs_interactor_human_A ~ IDs_domain_human_B, # set of interacting pairs to be tested (XZ), YZ interactions are assumed
node_attr = list(IDs_interactor_human_A ~ IDs_interactor_human_A_degree,
IDs_domain_human_B ~ domain_count),
data = data,
statistic = IDs_interactor_human_A + IDs_domain_human_B ~ .N,
select_nodes = IDs_domain_human_B ~ domain_count >= 1,
N = 180,
cores = NULL, seed = 2, also_permuteYZ = F,
clustermq = T, clustermq_mem = 20000,
split_comp_inner_N = 3, clustermq_jobs = 30,
clustermq_log_worker = F)}, times = 1)
plot(res, IDs_interactor_human_A + IDs_domain_human_B ~ log10(not_missing))
res$data_with_pval[!is.na(IDs_domain_human_B)]
{permutationPval(interactions2permute = IDs_interactor_viral ~ IDs_interactor_human,
associations2test = IDs_interactor_viral ~ IDs_domain_human,
node_attr = list(IDs_interactor_viral ~ IDs_interactor_viral_degree,
IDs_domain_human ~ domain_count,
IDs_interactor_viral + IDs_domain_human ~ domain_frequency_per_IDs_interactor_viral),
data = data,
statistic = IDs_interactor_viral + IDs_domain_human ~ .N,
select_nodes = IDs_domain_human ~ domain_count >= 1,
N = 1000,
cores = 3, seed = 2, clustermq = F)}
set.seed(1)
random = randomInteractome(n_prot = 200, degree_dist = NULL, taxid = "9606", database = "imex", protein_only = TRUE)
res_g = res
all.equal(res_g, res)
res2 = permutationPval(interactions2permute = IDs_interactor_viral ~ IDs_interactor_human,
associations2test = IDs_interactor_viral ~ IDs_domain_human,
node_attr = list(IDs_interactor_viral ~ IDs_interactor_viral_degree,
IDs_domain_human ~ domain_count,
IDs_interactor_viral + IDs_domain_human ~ domain_frequency_per_IDs_interactor_viral),
data = data,
statistic = IDs_interactor_viral + IDs_domain_human ~ .N / IDs_interactor_viral_degree,
select_nodes = IDs_domain_human ~ domain_count > 16,
N = 10,
cores = NULL, seed = NULL)
microbenchmark::microbenchmark({res = permutationPval(interactions2permute = IDs_interactor_viral ~ IDs_interactor_human,
associations2test = IDs_interactor_viral ~ IDs_domain_human,
node_attr = list(IDs_interactor_viral ~ IDs_interactor_viral_degree,
IDs_domain_human ~ domain_count,
IDs_interactor_viral + IDs_domain_human ~ domain_frequency_per_IDs_interactor_viral),
data = data,
statistic = IDs_interactor_viral + IDs_domain_human ~ .N / IDs_interactor_viral_degree,
select_nodes = IDs_domain_human ~ domain_count >= 1,
N = 10,
cores = NULL, seed = 1)}, times = 10)
profvis::profvis({resEnv = permutationPval(interactions2permute = IDs_interactor_viral ~ IDs_interactor_human,
associations2test = IDs_interactor_viral ~ IDs_domain_human,
node_attr = list(IDs_interactor_viral ~ IDs_interactor_viral_degree,
IDs_domain_human ~ domain_count,
IDs_interactor_viral + IDs_domain_human ~ domain_frequency_per_IDs_interactor_viral),
data = data,
statistic = IDs_interactor_viral + IDs_domain_human ~ .N / IDs_interactor_viral_degree,
select_nodes = IDs_domain_human ~ domain_count >= 1,
N = 10,
cores = NULL, seed = 1)})
# Fisher test
microbenchmark::microbenchmark({resFISHER = permutationPval(interactions2permute = IDs_interactor_viral ~ IDs_interactor_human,
associations2test = IDs_interactor_viral ~ IDs_domain_human,
node_attr = list(IDs_interactor_viral ~ IDs_interactor_viral_degree,
IDs_domain_human ~ domain_count + N_prot_w_interactors,
IDs_interactor_viral + IDs_domain_human ~ domain_count_per_IDs_interactor_viral),
data = data,
statistic = IDs_interactor_viral + IDs_domain_human ~ fisher.test(matrix(c(unique(domain_count), unique(N_prot_w_interactors) - unique(domain_count), unique(domain_count_per_IDs_interactor_viral), unique(IDs_interactor_viral_degree) - unique(domain_count_per_IDs_interactor_viral)),2,2), alternative = "greater", conf.int = F)$p.value,
select_nodes = IDs_domain_human ~ domain_count >= 1,
N = 100,
cores = NULL, seed = 1)}, times = 10)
qplot(x = resFISHER$data_with_pval[p.value < 0.5, IDs_interactor_viral_degree], y = resFISHER$data_with_pval[p.value < 0.5, domain_count], geom = "bin2d") + scale_x_log10() + scale_y_log10()
qplot(x = res$data_with_pval[p.value < 0.5, IDs_interactor_viral_degree], y = res$data_with_pval[p.value < 0.5, domain_count], geom = "bin2d") + scale_x_log10() + scale_y_log10()
qplot(x = resFISHER$data_with_pval[p.value < 0.01, IDs_interactor_viral_degree], y = resFISHER$data_with_pval[p.value < 0.01, domain_count], geom = "bin2d") + scale_x_log10() + scale_y_log10()
qplot(x = res$data_with_pval[p.value < 0.01, IDs_interactor_viral_degree], y = res$data_with_pval[p.value < 0.01, domain_count], geom = "bin2d") + scale_x_log10() + scale_y_log10()
# Unit: seconds (without inner and outer replicate)
#expr
#{ resFISHER = permutationPval(interactions2permute = IDs_interactor_viral ~ IDs_interactor_human, associations2test = IDs_interactor_viral ~ IDs_domain_human, node_attr = list(IDs_interactor_viral ~ IDs_interactor_viral_degree, IDs_domain_human ~ domain_count + N_prot_w_interactors, IDs_interactor_viral + IDs_domain_human ~ domain_count_per_IDs_interactor_viral), data = data, statistic = IDs_interactor_viral + IDs_domain_human ~ fisher.test(matrix(c(unique(domain_count), unique(N_prot_w_interactors) - unique(domain_count), unique(domain_count_per_IDs_interactor_viral), unique(IDs_interactor_viral_degree) - unique(domain_count_per_IDs_interactor_viral)), 2, 2), alternative = "greater", conf.int = F)$p.value, select_nodes = IDs_domain_human ~ domain_count >= 1, N = 100, cores = NULL, seed = 1) }
#min lq mean median uq max neval
#34.89687 35.35519 35.3463 35.39729 35.42503 35.49184 10
# Unit: seconds (with inner and outer replicate)
# expr
# { resFISHER = permutationPval(interactions2permute = IDs_interactor_viral ~ IDs_interactor_human, associations2test = IDs_interactor_viral ~ IDs_domain_human, node_attr = list(IDs_interactor_viral ~ IDs_interactor_viral_degree, IDs_domain_human ~ domain_count + N_prot_w_interactors, IDs_interactor_viral + IDs_domain_human ~ domain_count_per_IDs_interactor_viral), data = data, statistic = IDs_interactor_viral + IDs_domain_human ~ fisher.test(matrix(c(unique(domain_count), unique(N_prot_w_interactors) - unique(domain_count), unique(domain_count_per_IDs_interactor_viral), unique(IDs_interactor_viral_degree) - unique(domain_count_per_IDs_interactor_viral)), 2, 2), alternative = "greater", conf.int = F)$p.value, select_nodes = IDs_domain_human ~ domain_count >= 1, N = 100, cores = NULL, seed = 1) }
# min lq mean median uq max neval
# 38.02992 38.33805 38.7546 38.73306 38.86937 39.60011 10
library(PItools)
library(rtracklayer)
library(ggplot2)
data = fread("../viral_project/processed_data_files/viral_human_net_w_domains", sep = "\t", stringsAsFactors = F)
time = proc.time()
res = permutationPval(interactions2permute = IDs_interactor_viral ~ IDs_interactor_human,
associations2test = IDs_interactor_viral ~ IDs_domain_human,
node_attr = list(IDs_interactor_viral ~ IDs_interactor_viral_degree,
IDs_domain_human ~ domain_count,
IDs_interactor_viral + IDs_domain_human ~ domain_frequency_per_IDs_interactor_viral),
data = data,
statistic = IDs_interactor_viral + IDs_domain_human ~ .N / IDs_interactor_viral_degree,
select_nodes = IDs_domain_human ~ domain_count >= 1,
N = 10000,
cores = NULL, seed = 2)
proc.time() - time
time = proc.time()
resFISHER = permutationPval(interactions2permute = IDs_interactor_viral ~ IDs_interactor_human,
associations2test = IDs_interactor_viral ~ IDs_domain_human,
node_attr = list(IDs_interactor_viral ~ IDs_interactor_viral_degree,
IDs_domain_human ~ domain_count + N_prot_w_interactors,
IDs_interactor_viral + IDs_domain_human ~ domain_count_per_IDs_interactor_viral),
data = data,
statistic = IDs_interactor_viral + IDs_domain_human ~ fisher.test(matrix(c(unique(domain_count), unique(N_prot_w_interactors) - unique(domain_count), unique(domain_count_per_IDs_interactor_viral), unique(IDs_interactor_viral_degree) - unique(domain_count_per_IDs_interactor_viral)),2,2), alternative = "greater", conf.int = F)$p.value,
select_nodes = IDs_domain_human ~ domain_count >= 1,
N = 10000,
cores = NULL, seed = 1)
resFISHER$data_with_pval[, p.value := 1 - p.value]
proc.time() - time
resFISHER$IDs_interactor_viral_degreeVSdomain_count = qplot(x = resFISHER$data_with_pval[order(p.value)[1:250], IDs_interactor_viral_degree], y = resFISHER$data_with_pval[order(p.value)[1:250], domain_count], geom = "bin2d") + scale_x_log10() + scale_y_log10()
res$IDs_interactor_viral_degreeVSdomain_count = qplot(x = res$data_with_pval[order(p.value)[1:250], IDs_interactor_viral_degree], y = res$data_with_pval[order(p.value)[1:250], domain_count], geom = "bin2d") + scale_x_log10() + scale_y_log10()
interactiondomains = fread("http://elm.eu.org/interactiondomains.tsv")
interactiondomains[, pfam_id := `Interaction Domain Id`]
domains_known = interactiondomains[, unique(pfam_id)]
"../viral_project/processed_data_files/InterProScan_domains_nonredundant.gff3"
"../viral_project/processed_data_files/all_human_viral_protein_domains.gff3.gz"
InterProScan_domains_nonred = import(con = "../viral_project/processed_data_files/InterProScan_domains_nonredundant.gff3", format = "gff3")
domains_mapping = unique(data.table(any_id = as.character(InterProScan_domains_nonred$Name), interpro_id = as.character(InterProScan_domains_nonred$Dbxref)))
domains_known_mapped = unique(domains_mapping[any_id %in% domains_known, interpro_id])
domains_not_mapped = unique(domains_known[!domains_known %in% domains_mapping$any_id])
test_enrichment = function(N, res, domains_known_mapped){
res$data_pval = unique(res$data_with_pval[,.(IDs_interactor_viral, IDs_domain_human, p.value, domain_type, domain_count, IDs_interactor_viral_degree)])
res$data_pval[, pval_fdr := p.adjust(p.value, method = "fdr")]
hist(res$data_pval[, pval_fdr], breaks = seq(0,1,0.01))
domains_found = res$data_pval[order(p.value)[1:N], unique(IDs_domain_human)]
alldomains = res$data_pval[, unique(IDs_domain_human)]
known = factor(alldomains %in% domains_known_mapped, levels = c("TRUE", "FALSE"))
found = factor(alldomains %in% domains_found, levels = c("TRUE", "FALSE"))
table_res = table(known, found)
test = fisher.test(table(known, found), alternative = "greater", conf.int = T)
return(c(test$p.value, test$estimate, table_res["TRUE", "TRUE"]))
}
enrichment = sapply(seq(25, 500, 25), test_enrichment, res, domains_known_mapped)
colnames(enrichment) = seq(25, 500, 25)
enrichmentFISHER = sapply(seq(25, 500, 25), test_enrichment, resFISHER, domains_known_mapped)
colnames(enrichmentFISHER) = seq(25, 500, 25)
plot(colnames(enrichment), enrichment[2,], ylab = "Fisher test odds ratio", xlab = "top N viral protein - domain pairs selected", col = "red", type = "l", ylim = c(0,18))
lines(x = colnames(enrichment), y = enrichmentFISHER[2,], col = "blue", type = "l")
legend(x = 80, y = 17.5, c("statictic used in permutation test:","domain frequency among interactors of a viral protein", "Fisher test pval: domain overrepresentation over the background"), col = c("white","red", "blue"), lty = 1 ,merge = TRUE)
plot(colnames(enrichment), enrichment[3,], ylab = "known domain found", xlab = "top N viral protein - domain pairs selected", col = "red", type = "l", ylim = c(0,length(domains_known_mapped)+1))
lines(x = colnames(enrichment), y = enrichmentFISHER[3,], col = "blue", type = "l")
abline(h = length(domains_known_mapped), col = "green")
legend(x = 80, y = 50, c("statictic used in permutation test:","domain frequency among interactors of a viral protein", "Fisher test pval: domain overrepresentation over the background", "domains known to interact with linear motifs"), col = c("white","red", "blue", "green"), lty = 1 , merge = TRUE)
plot(colnames(enrichment), enrichment[1,], ylab = "Fisher test pvalue", xlab = "top N viral protein - domain pairs selected", col = "red", type = "l", ylim = c(0,0.004))
lines(x = colnames(enrichment), y = enrichmentFISHER[1,], col = "blue", type = "l")
legend(x = 80, y = 0.0041, c("statictic used in permutation test:","domain frequency among interactors of a viral protein", "Fisher test pval: domain overrepresentation over the background"), col = c("white","red", "blue"), lty = 1 ,merge = TRUE)
big_jobs = sapply(list.files(), function(file){length(readLines(file))}) == 1
error_paths = paste0("/hps/nobackup/research/petsalaki/users/vitalii/vitalii/viral_project/qslimfinder.Full_IntAct4.FALSE/log_dir/error/", gsub("\\.sh","", names(big_jobs)[big_jobs]))
sapply(error_paths, function(error_path) {
system(paste0("cat ", error_path," | grep Terminated"), intern=T)
})
log_paths = paste0("/hps/nobackup/research/petsalaki/users/vitalii/vitalii/viral_project/qslimfinder.Full_IntAct4.FALSE/log_dir/log/", gsub("\\.sh","", names(big_jobs)[big_jobs]))
job_status = sapply(log_paths, function(log_path) {
system(paste0("cat ", log_path," | grep TERM_MEMLIMIT"), intern=T)
})
# how many terminated because of memory
sum(job_status == "TERM_MEMLIMIT: job killed after reaching LSF memory usage limit.") # 9
sum(sapply(job_status, function(element) {
length(element) == 0
})) # 65
all.jobs = sapply(list.files(), function(file){length(readLines(file))})
error_paths = paste0("/hps/nobackup/research/petsalaki/users/vitalii/vitalii/viral_project/qslimfinder.Full_IntAct4.FALSE/log_dir/error/", gsub("\\.sh","", names(all.jobs)))
errors = sapply(error_paths, function(error_path) {
system(paste0("cat ", error_path," | grep Terminated"), intern=T)
})
log_paths = paste0("/hps/nobackup/research/petsalaki/users/vitalii/vitalii/viral_project/qslimfinder.Full_IntAct4.FALSE/log_dir/log/", gsub("\\.sh","", names(all.jobs)))
successes = sapply(log_paths, function(log_path) {
system(paste0("cat ", log_path," | grep \"Successfully completed\""), intern=T)
})
## TEST for benchmarking
occurence_file = "../viral_project/qslimfinder.Full_IntAct3.FALSE/result/occurence.txt"
main_file = "../viral_project/qslimfinder.Full_IntAct3.FALSE/result/main_result.txt"
domain_res_file = "../viral_project/processed_data_files/domain_res_count_20171019.RData"
motif_setup = "../viral_project/processed_data_files/QSLIMFinder_instances_h2v_qslimfinder.Full_IntAct3.FALSE_clust201802.RData"
neg_set = c("all_instances", "all_proteins", "random")[1]
domain_results_obj = "res_count"
motif_input_obj = "forSLIMFinder_Ready"
motif_setup_obj2 = NULL
occurence_filt = NULL
one_from_cloud = T
dbfile_main = "../viral_project/data_files/instances_all.gff"
dburl_main = "http://elm.eu.org/instances.gff?q=None&taxon=Homo%20sapiens&instance_logic="
dbfile_query = "../viral_project/data_files/instances_query.gff"
dburl_query = "http://elm.eu.org/instances.gff?q=all&taxon=irus&instance_logic="
query_res_query_only = T
motif_types = c("DOC", "MOD", "LIG")
all_res_excl_query = T
merge_motif_variants = F
seed = 21
N = 100
replace = T
within1sequence = T
query_predictor_col = "Sig"
all_predictor_col = "Sig"
normalise = T
minoverlap = 2
minoverlap_redundant = 5
merge_domain_data = T
merge_by_occurence_mcols = c("query", "interacts_with")
merge_by_domain_res_cols = c("IDs_interactor_viral", "IDs_interactor_human", "IDs_domain_human", "Taxid_interactor_human","Taxid_interactor_viral")
merge_by_non_query_domain_res_cols = c("IDs_interactor_human_A", "IDs_interactor_human_B", "IDs_domain_human_B", "Taxid_interactor_human_A","Taxid_interactor_human_B")
count_ranges_by = list(by = "IDs_domain_human", name = "motif_occ_per_domain",
normalise_by = "domain_count", normalised_name = "normalised_motif_occ_per_domain")
filter_by_domain_data = "p.value < 0.5"
motif_pval_cutoff = 1
select_predictor_per_range = max
non_query_domain_res_file = "../viral_project/processed_data_files/predict_domain_human_clust20180819.RData"
non_query_domain_results_obj = "res_count_all" # NULL res_count_all
non_query_domains_N = 0
non_query_set_only = F
query_domains_only = F
min_non_query_domain_support = 3
select_top_domain = F
# creating a website
# Install release version from CRAN
install.packages("pkgdown")
setwd("/Users/vk7/Desktop/ebi_projects/PItools/")
## benchmarking ways to combine p-values for domains
prior = mcol[p.value<1,mean(1-p.value)]
mcol[, combined_p.value := updatePvalue(prior, p.value), by = .(query, Pattern, interacts_with, domain)]
dom_elm = merge(mcol, elm_interactions,
by.x = c("query", "interacts_with", "domain"),
by.y = c("interactorElm", "interactorDomain", "InterProID"),
all.x = T, all.y = F)
dom_elm[!is.na(Elm), Correct_SLiM_binding := "yes"]
dom_elm[is.na(Elm), Correct_SLiM_binding := "no"]
# filter and select unique to have one row per protein_with_motif - domain pair
dom3 = unique(dom_elm[,.(query, domain, p.value = combined_p.value, Correct_SLiM_binding)])
dom3[, both_SLIM_interaction := uniqueN(Correct_SLiM_binding) >= 2,
by = .(query)]
dom3 = dom3[both_SLIM_interaction == T]
# calculate N domain protein pairs
dom3[, N_pairs := uniqueN(paste0(query, domain)), by = .(Correct_SLiM_binding)]
dom3[, network := paste0(unique(N_pairs[Correct_SLiM_binding == "yes"]),
" correct domains ",
"/ ", unique(N_pairs[Correct_SLiM_binding == "no"]),
" total")]
ggplot(dom3, aes(p.value, color = Correct_SLiM_binding, fill = Correct_SLiM_binding)) +
geom_histogram(alpha = 0.3, aes(y = ..density..), bins = 50, position = "identity") +
theme_bw() +
facet_grid(network~., scales = "free_y")+
theme(legend.position = "none",
strip.text.y = element_text(angle = 0),
strip.text = element_text(size = 14),
legend.title = element_text(size = 14),
legend.text = element_text(size = 14),
axis.title = element_text(size = 14),
axis.text = element_text(size = 14),
text = element_text(size = 14))
### benchmark end
prior = mcol[p.value<1,mean(1-p.value)]
prior
updatePvalue = function(p.value1, p.value2, N = 1, incr = 1e-10){
N = max(length(p.value1), length(p.value2))
# add/substract small value to p-value to avoid 0 likelihood
p.value2[p.value2 == 1] = 1 - incr
p.value2[p.value2 == 0] = incr
# calculate likelihood
likelihood_of_2_1 = dbeta(p.value2, p.value1 * N, (1 - p.value1) * N)
likelihood_of_2_not_1 = dbeta(p.value2, (1 - p.value1) * N, p.value1 * N)
p.value = sum(likelihood_of_2_1 * p.value1) /
(sum(likelihood_of_2_1 * p.value1) +
sum(likelihood_of_2_not_1 * (1 - p.value1)))
p.value
}
probs = c(prior, 0.999, 0, 0.979)
b = updatePvalue(probs[1], probs[2:4], N = 1)
b
up2 = updatePvalue(up1$p.value, probs[3], N = up1$N)
up3 = updatePvalue(up2$p.value, probs[4], N = up2$N)
up3
combn(c(0.999, 0, 0.979),2)
fit = beta.mle(c(0.999, 0.998, 0.979), tol = 1e-09)
hist(rbeta(10000, 0.4999 + 100, 0.4999 + 100))
mean(fit$param[1] / (fit$param[1] + fit$param[2]))
### draft end
pkgdown::build_site()
vitkl/MItools documentation built on May 29, 2019, 2:55 p.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(roxygen2)
library(devtools)
document()
install()
library(PItools)
library(GenomicRanges)
library(R.utils)
data = fread("../viral_project/processed_data_files/viral_human_net_w_domains", sep = "\t", stringsAsFactors = F)
file_BioPlex3 = "../viral_project/processed_data_files/human_net_w_domains"
gunzip(paste0(file_BioPlex3,".gz"), remove = F)
data = fread(file_BioPlex3, sep = "\t", stringsAsFactors = F)
data = data[IDs_interactor_human_B != "UNKNOWN"]
unlink(file_BioPlex3)
res = permutationPval(interactions2permute = IDs_interactor_viral ~ IDs_interactor_human,
associations2test = IDs_interactor_viral ~ IDs_domain_human,
node_attr = list(IDs_interactor_viral ~ IDs_interactor_viral_degree,
IDs_domain_human ~ domain_count,
IDs_interactor_viral + IDs_domain_human ~ domain_frequency_per_IDs_interactor_viral),
data = data,
statistic = IDs_interactor_viral + IDs_domain_human ~ .N,
select_nodes = IDs_domain_human ~ domain_count >= 1,
N = 100, clustermq_jobs = 20, split_comp_inner_N = 2,
cores = NULL, seed = 2, clustermq = T, clustermq_mem = 2000)
plot(res)
res
library(lineprof)
library(PItools)
library(GenomicRanges)
library(R.utils)
data = fread("../viral_project/processed_data_files/viral_human_net_w_domains", sep = "\t", stringsAsFactors = F)
data_num = copy(data)
data_num[, IDs_interactor_human := as.integer(as.factor(IDs_interactor_human))]
data_num[, IDs_domain_human := as.integer(as.factor(IDs_domain_human))]
data_num[, IDs_interactor_viral := as.integer(as.factor(IDs_interactor_viral))]
l = microbenchmark::microbenchmark({
res <<- permutationPval(interactions2permute = IDs_interactor_viral ~ IDs_interactor_human,
associations2test = IDs_interactor_viral ~ IDs_domain_human,
node_attr = list(IDs_interactor_viral ~ IDs_interactor_viral_degree,
IDs_domain_human ~ domain_count,
IDs_interactor_viral + IDs_domain_human ~ domain_frequency_per_IDs_interactor_viral),
data = data,
statistic = IDs_interactor_viral + IDs_domain_human ~ .N,
select_nodes = IDs_domain_human ~ domain_count >= 1,
N = 10,
cores = 1, seed = 2, clustermq = F)
}, {
res_num <<- permutationPval(interactions2permute = IDs_interactor_viral ~ IDs_interactor_human,
associations2test = IDs_interactor_viral ~ IDs_domain_human,
node_attr = list(IDs_interactor_viral ~ IDs_interactor_viral_degree,
IDs_domain_human ~ domain_count,
IDs_interactor_viral + IDs_domain_human ~ domain_frequency_per_IDs_interactor_viral),
data = data_num,
statistic = IDs_interactor_viral + IDs_domain_human ~ .N,
select_nodes = IDs_domain_human ~ domain_count >= 1,
N = 10,
cores = 1, seed = 2, clustermq = F)
}, times = 12)
l
microbenchmark::microbenchmark({res <- permutationPval(interactions2permute = IDs_interactor_human_A ~ IDs_interactor_human_B, # first set of interacting pairs (XY) that are to be permuted
associations2test = IDs_interactor_human_A ~ IDs_domain_human_B, # set of interacting pairs to be tested (XZ), YZ interactions are assumed
node_attr = list(IDs_interactor_human_A ~ IDs_interactor_human_A_degree,
IDs_domain_human_B ~ domain_count),
data = data,
statistic = IDs_interactor_human_A + IDs_domain_human_B ~ .N,
select_nodes = IDs_domain_human_B ~ domain_count >= 1,
N = 180,
cores = NULL, seed = 2, also_permuteYZ = F,
clustermq = T, clustermq_mem = 20000,
split_comp_inner_N = 3, clustermq_jobs = 30,
clustermq_log_worker = F)}, times = 1)
plot(res, IDs_interactor_human_A + IDs_domain_human_B ~ log10(not_missing))
res$data_with_pval[!is.na(IDs_domain_human_B)]
{permutationPval(interactions2permute = IDs_interactor_viral ~ IDs_interactor_human,
associations2test = IDs_interactor_viral ~ IDs_domain_human,
node_attr = list(IDs_interactor_viral ~ IDs_interactor_viral_degree,
IDs_domain_human ~ domain_count,
IDs_interactor_viral + IDs_domain_human ~ domain_frequency_per_IDs_interactor_viral),
data = data,
statistic = IDs_interactor_viral + IDs_domain_human ~ .N,
select_nodes = IDs_domain_human ~ domain_count >= 1,
N = 1000,
cores = 3, seed = 2, clustermq = F)}
set.seed(1)
random = randomInteractome(n_prot = 200, degree_dist = NULL, taxid = "9606", database = "imex", protein_only = TRUE)
res_g = res
all.equal(res_g, res)
res2 = permutationPval(interactions2permute = IDs_interactor_viral ~ IDs_interactor_human,
associations2test = IDs_interactor_viral ~ IDs_domain_human,
node_attr = list(IDs_interactor_viral ~ IDs_interactor_viral_degree,
IDs_domain_human ~ domain_count,
IDs_interactor_viral + IDs_domain_human ~ domain_frequency_per_IDs_interactor_viral),
data = data,
statistic = IDs_interactor_viral + IDs_domain_human ~ .N / IDs_interactor_viral_degree,
select_nodes = IDs_domain_human ~ domain_count > 16,
N = 10,
cores = NULL, seed = NULL)
microbenchmark::microbenchmark({res = permutationPval(interactions2permute = IDs_interactor_viral ~ IDs_interactor_human,
associations2test = IDs_interactor_viral ~ IDs_domain_human,
node_attr = list(IDs_interactor_viral ~ IDs_interactor_viral_degree,
IDs_domain_human ~ domain_count,
IDs_interactor_viral + IDs_domain_human ~ domain_frequency_per_IDs_interactor_viral),
data = data,
statistic = IDs_interactor_viral + IDs_domain_human ~ .N / IDs_interactor_viral_degree,
select_nodes = IDs_domain_human ~ domain_count >= 1,
N = 10,
cores = NULL, seed = 1)}, times = 10)
profvis::profvis({resEnv = permutationPval(interactions2permute = IDs_interactor_viral ~ IDs_interactor_human,
associations2test = IDs_interactor_viral ~ IDs_domain_human,
node_attr = list(IDs_interactor_viral ~ IDs_interactor_viral_degree,
IDs_domain_human ~ domain_count,
IDs_interactor_viral + IDs_domain_human ~ domain_frequency_per_IDs_interactor_viral),
data = data,
statistic = IDs_interactor_viral + IDs_domain_human ~ .N / IDs_interactor_viral_degree,
select_nodes = IDs_domain_human ~ domain_count >= 1,
N = 10,
cores = NULL, seed = 1)})
# Fisher test
microbenchmark::microbenchmark({resFISHER = permutationPval(interactions2permute = IDs_interactor_viral ~ IDs_interactor_human,
associations2test = IDs_interactor_viral ~ IDs_domain_human,
node_attr = list(IDs_interactor_viral ~ IDs_interactor_viral_degree,
IDs_domain_human ~ domain_count + N_prot_w_interactors,
IDs_interactor_viral + IDs_domain_human ~ domain_count_per_IDs_interactor_viral),
data = data,
statistic = IDs_interactor_viral + IDs_domain_human ~ fisher.test(matrix(c(unique(domain_count), unique(N_prot_w_interactors) - unique(domain_count), unique(domain_count_per_IDs_interactor_viral), unique(IDs_interactor_viral_degree) - unique(domain_count_per_IDs_interactor_viral)),2,2), alternative = "greater", conf.int = F)$p.value,
select_nodes = IDs_domain_human ~ domain_count >= 1,
N = 100,
cores = NULL, seed = 1)}, times = 10)
qplot(x = resFISHER$data_with_pval[p.value < 0.5, IDs_interactor_viral_degree], y = resFISHER$data_with_pval[p.value < 0.5, domain_count], geom = "bin2d") + scale_x_log10() + scale_y_log10()
qplot(x = res$data_with_pval[p.value < 0.5, IDs_interactor_viral_degree], y = res$data_with_pval[p.value < 0.5, domain_count], geom = "bin2d") + scale_x_log10() + scale_y_log10()
qplot(x = resFISHER$data_with_pval[p.value < 0.01, IDs_interactor_viral_degree], y = resFISHER$data_with_pval[p.value < 0.01, domain_count], geom = "bin2d") + scale_x_log10() + scale_y_log10()
qplot(x = res$data_with_pval[p.value < 0.01, IDs_interactor_viral_degree], y = res$data_with_pval[p.value < 0.01, domain_count], geom = "bin2d") + scale_x_log10() + scale_y_log10()
# Unit: seconds (without inner and outer replicate)
#expr
#{ resFISHER = permutationPval(interactions2permute = IDs_interactor_viral ~ IDs_interactor_human, associations2test = IDs_interactor_viral ~ IDs_domain_human, node_attr = list(IDs_interactor_viral ~ IDs_interactor_viral_degree, IDs_domain_human ~ domain_count + N_prot_w_interactors, IDs_interactor_viral + IDs_domain_human ~ domain_count_per_IDs_interactor_viral), data = data, statistic = IDs_interactor_viral + IDs_domain_human ~ fisher.test(matrix(c(unique(domain_count), unique(N_prot_w_interactors) - unique(domain_count), unique(domain_count_per_IDs_interactor_viral), unique(IDs_interactor_viral_degree) - unique(domain_count_per_IDs_interactor_viral)), 2, 2), alternative = "greater", conf.int = F)$p.value, select_nodes = IDs_domain_human ~ domain_count >= 1, N = 100, cores = NULL, seed = 1) }
#min lq mean median uq max neval
#34.89687 35.35519 35.3463 35.39729 35.42503 35.49184 10
# Unit: seconds (with inner and outer replicate)
# expr
# { resFISHER = permutationPval(interactions2permute = IDs_interactor_viral ~ IDs_interactor_human, associations2test = IDs_interactor_viral ~ IDs_domain_human, node_attr = list(IDs_interactor_viral ~ IDs_interactor_viral_degree, IDs_domain_human ~ domain_count + N_prot_w_interactors, IDs_interactor_viral + IDs_domain_human ~ domain_count_per_IDs_interactor_viral), data = data, statistic = IDs_interactor_viral + IDs_domain_human ~ fisher.test(matrix(c(unique(domain_count), unique(N_prot_w_interactors) - unique(domain_count), unique(domain_count_per_IDs_interactor_viral), unique(IDs_interactor_viral_degree) - unique(domain_count_per_IDs_interactor_viral)), 2, 2), alternative = "greater", conf.int = F)$p.value, select_nodes = IDs_domain_human ~ domain_count >= 1, N = 100, cores = NULL, seed = 1) }
# min lq mean median uq max neval
# 38.02992 38.33805 38.7546 38.73306 38.86937 39.60011 10
library(PItools)
library(rtracklayer)
library(ggplot2)
data = fread("../viral_project/processed_data_files/viral_human_net_w_domains", sep = "\t", stringsAsFactors = F)
time = proc.time()
res = permutationPval(interactions2permute = IDs_interactor_viral ~ IDs_interactor_human,
associations2test = IDs_interactor_viral ~ IDs_domain_human,
node_attr = list(IDs_interactor_viral ~ IDs_interactor_viral_degree,
IDs_domain_human ~ domain_count,
IDs_interactor_viral + IDs_domain_human ~ domain_frequency_per_IDs_interactor_viral),
data = data,
statistic = IDs_interactor_viral + IDs_domain_human ~ .N / IDs_interactor_viral_degree,
select_nodes = IDs_domain_human ~ domain_count >= 1,
N = 10000,
cores = NULL, seed = 2)
proc.time() - time
time = proc.time()
resFISHER = permutationPval(interactions2permute = IDs_interactor_viral ~ IDs_interactor_human,
associations2test = IDs_interactor_viral ~ IDs_domain_human,
node_attr = list(IDs_interactor_viral ~ IDs_interactor_viral_degree,
IDs_domain_human ~ domain_count + N_prot_w_interactors,
IDs_interactor_viral + IDs_domain_human ~ domain_count_per_IDs_interactor_viral),
data = data,
statistic = IDs_interactor_viral + IDs_domain_human ~ fisher.test(matrix(c(unique(domain_count), unique(N_prot_w_interactors) - unique(domain_count), unique(domain_count_per_IDs_interactor_viral), unique(IDs_interactor_viral_degree) - unique(domain_count_per_IDs_interactor_viral)),2,2), alternative = "greater", conf.int = F)$p.value,
select_nodes = IDs_domain_human ~ domain_count >= 1,
N = 10000,
cores = NULL, seed = 1)
resFISHER$data_with_pval[, p.value := 1 - p.value]
proc.time() - time
resFISHER$IDs_interactor_viral_degreeVSdomain_count = qplot(x = resFISHER$data_with_pval[order(p.value)[1:250], IDs_interactor_viral_degree], y = resFISHER$data_with_pval[order(p.value)[1:250], domain_count], geom = "bin2d") + scale_x_log10() + scale_y_log10()
res$IDs_interactor_viral_degreeVSdomain_count = qplot(x = res$data_with_pval[order(p.value)[1:250], IDs_interactor_viral_degree], y = res$data_with_pval[order(p.value)[1:250], domain_count], geom = "bin2d") + scale_x_log10() + scale_y_log10()
interactiondomains = fread("http://elm.eu.org/interactiondomains.tsv")
interactiondomains[, pfam_id := `Interaction Domain Id`]
domains_known = interactiondomains[, unique(pfam_id)]
"../viral_project/processed_data_files/InterProScan_domains_nonredundant.gff3"
"../viral_project/processed_data_files/all_human_viral_protein_domains.gff3.gz"
InterProScan_domains_nonred = import(con = "../viral_project/processed_data_files/InterProScan_domains_nonredundant.gff3", format = "gff3")
domains_mapping = unique(data.table(any_id = as.character(InterProScan_domains_nonred$Name), interpro_id = as.character(InterProScan_domains_nonred$Dbxref)))
domains_known_mapped = unique(domains_mapping[any_id %in% domains_known, interpro_id])
domains_not_mapped = unique(domains_known[!domains_known %in% domains_mapping$any_id])
test_enrichment = function(N, res, domains_known_mapped){
res$data_pval = unique(res$data_with_pval[,.(IDs_interactor_viral, IDs_domain_human, p.value, domain_type, domain_count, IDs_interactor_viral_degree)])
res$data_pval[, pval_fdr := p.adjust(p.value, method = "fdr")]
hist(res$data_pval[, pval_fdr], breaks = seq(0,1,0.01))
domains_found = res$data_pval[order(p.value)[1:N], unique(IDs_domain_human)]
alldomains = res$data_pval[, unique(IDs_domain_human)]
known = factor(alldomains %in% domains_known_mapped, levels = c("TRUE", "FALSE"))
found = factor(alldomains %in% domains_found, levels = c("TRUE", "FALSE"))
table_res = table(known, found)
test = fisher.test(table(known, found), alternative = "greater", conf.int = T)
return(c(test$p.value, test$estimate, table_res["TRUE", "TRUE"]))
}
enrichment = sapply(seq(25, 500, 25), test_enrichment, res, domains_known_mapped)
colnames(enrichment) = seq(25, 500, 25)
enrichmentFISHER = sapply(seq(25, 500, 25), test_enrichment, resFISHER, domains_known_mapped)
colnames(enrichmentFISHER) = seq(25, 500, 25)
plot(colnames(enrichment), enrichment[2,], ylab = "Fisher test odds ratio", xlab = "top N viral protein - domain pairs selected", col = "red", type = "l", ylim = c(0,18))
lines(x = colnames(enrichment), y = enrichmentFISHER[2,], col = "blue", type = "l")
legend(x = 80, y = 17.5, c("statictic used in permutation test:","domain frequency among interactors of a viral protein", "Fisher test pval: domain overrepresentation over the background"), col = c("white","red", "blue"), lty = 1 ,merge = TRUE)
plot(colnames(enrichment), enrichment[3,], ylab = "known domain found", xlab = "top N viral protein - domain pairs selected", col = "red", type = "l", ylim = c(0,length(domains_known_mapped)+1))
lines(x = colnames(enrichment), y = enrichmentFISHER[3,], col = "blue", type = "l")
abline(h = length(domains_known_mapped), col = "green")
legend(x = 80, y = 50, c("statictic used in permutation test:","domain frequency among interactors of a viral protein", "Fisher test pval: domain overrepresentation over the background", "domains known to interact with linear motifs"), col = c("white","red", "blue", "green"), lty = 1 , merge = TRUE)
plot(colnames(enrichment), enrichment[1,], ylab = "Fisher test pvalue", xlab = "top N viral protein - domain pairs selected", col = "red", type = "l", ylim = c(0,0.004))
lines(x = colnames(enrichment), y = enrichmentFISHER[1,], col = "blue", type = "l")
legend(x = 80, y = 0.0041, c("statictic used in permutation test:","domain frequency among interactors of a viral protein", "Fisher test pval: domain overrepresentation over the background"), col = c("white","red", "blue"), lty = 1 ,merge = TRUE)
big_jobs = sapply(list.files(), function(file){length(readLines(file))}) == 1
error_paths = paste0("/hps/nobackup/research/petsalaki/users/vitalii/vitalii/viral_project/qslimfinder.Full_IntAct4.FALSE/log_dir/error/", gsub("\\.sh","", names(big_jobs)[big_jobs]))
sapply(error_paths, function(error_path) {
system(paste0("cat ", error_path," | grep Terminated"), intern=T)
})
log_paths = paste0("/hps/nobackup/research/petsalaki/users/vitalii/vitalii/viral_project/qslimfinder.Full_IntAct4.FALSE/log_dir/log/", gsub("\\.sh","", names(big_jobs)[big_jobs]))
job_status = sapply(log_paths, function(log_path) {
system(paste0("cat ", log_path," | grep TERM_MEMLIMIT"), intern=T)
})
# how many terminated because of memory
sum(job_status == "TERM_MEMLIMIT: job killed after reaching LSF memory usage limit.") # 9
sum(sapply(job_status, function(element) {
length(element) == 0
})) # 65
all.jobs = sapply(list.files(), function(file){length(readLines(file))})
error_paths = paste0("/hps/nobackup/research/petsalaki/users/vitalii/vitalii/viral_project/qslimfinder.Full_IntAct4.FALSE/log_dir/error/", gsub("\\.sh","", names(all.jobs)))
errors = sapply(error_paths, function(error_path) {
system(paste0("cat ", error_path," | grep Terminated"), intern=T)
})
log_paths = paste0("/hps/nobackup/research/petsalaki/users/vitalii/vitalii/viral_project/qslimfinder.Full_IntAct4.FALSE/log_dir/log/", gsub("\\.sh","", names(all.jobs)))
successes = sapply(log_paths, function(log_path) {
system(paste0("cat ", log_path," | grep \"Successfully completed\""), intern=T)
})
## TEST for benchmarking
occurence_file = "../viral_project/qslimfinder.Full_IntAct3.FALSE/result/occurence.txt"
main_file = "../viral_project/qslimfinder.Full_IntAct3.FALSE/result/main_result.txt"
domain_res_file = "../viral_project/processed_data_files/domain_res_count_20171019.RData"
motif_setup = "../viral_project/processed_data_files/QSLIMFinder_instances_h2v_qslimfinder.Full_IntAct3.FALSE_clust201802.RData"
neg_set = c("all_instances", "all_proteins", "random")[1]
domain_results_obj = "res_count"
motif_input_obj = "forSLIMFinder_Ready"
motif_setup_obj2 = NULL
occurence_filt = NULL
one_from_cloud = T
dbfile_main = "../viral_project/data_files/instances_all.gff"
dburl_main = "http://elm.eu.org/instances.gff?q=None&taxon=Homo%20sapiens&instance_logic="
dbfile_query = "../viral_project/data_files/instances_query.gff"
dburl_query = "http://elm.eu.org/instances.gff?q=all&taxon=irus&instance_logic="
query_res_query_only = T
motif_types = c("DOC", "MOD", "LIG")
all_res_excl_query = T
merge_motif_variants = F
seed = 21
N = 100
replace = T
within1sequence = T
query_predictor_col = "Sig"
all_predictor_col = "Sig"
normalise = T
minoverlap = 2
minoverlap_redundant = 5
merge_domain_data = T
merge_by_occurence_mcols = c("query", "interacts_with")
merge_by_domain_res_cols = c("IDs_interactor_viral", "IDs_interactor_human", "IDs_domain_human", "Taxid_interactor_human","Taxid_interactor_viral")
merge_by_non_query_domain_res_cols = c("IDs_interactor_human_A", "IDs_interactor_human_B", "IDs_domain_human_B", "Taxid_interactor_human_A","Taxid_interactor_human_B")
count_ranges_by = list(by = "IDs_domain_human", name = "motif_occ_per_domain",
normalise_by = "domain_count", normalised_name = "normalised_motif_occ_per_domain")
filter_by_domain_data = "p.value < 0.5"
motif_pval_cutoff = 1
select_predictor_per_range = max
non_query_domain_res_file = "../viral_project/processed_data_files/predict_domain_human_clust20180819.RData"
non_query_domain_results_obj = "res_count_all" # NULL res_count_all
non_query_domains_N = 0
non_query_set_only = F
query_domains_only = F
min_non_query_domain_support = 3
select_top_domain = F
# creating a website
# Install release version from CRAN
install.packages("pkgdown")
setwd("/Users/vk7/Desktop/ebi_projects/PItools/")
## benchmarking ways to combine p-values for domains
prior = mcol[p.value<1,mean(1-p.value)]
mcol[, combined_p.value := updatePvalue(prior, p.value), by = .(query, Pattern, interacts_with, domain)]
dom_elm = merge(mcol, elm_interactions,
by.x = c("query", "interacts_with", "domain"),
by.y = c("interactorElm", "interactorDomain", "InterProID"),
all.x = T, all.y = F)
dom_elm[!is.na(Elm), Correct_SLiM_binding := "yes"]
dom_elm[is.na(Elm), Correct_SLiM_binding := "no"]
# filter and select unique to have one row per protein_with_motif - domain pair
dom3 = unique(dom_elm[,.(query, domain, p.value = combined_p.value, Correct_SLiM_binding)])
dom3[, both_SLIM_interaction := uniqueN(Correct_SLiM_binding) >= 2,
by = .(query)]
dom3 = dom3[both_SLIM_interaction == T]
# calculate N domain protein pairs
dom3[, N_pairs := uniqueN(paste0(query, domain)), by = .(Correct_SLiM_binding)]
dom3[, network := paste0(unique(N_pairs[Correct_SLiM_binding == "yes"]),
" correct domains ",
"/ ", unique(N_pairs[Correct_SLiM_binding == "no"]),
" total")]
ggplot(dom3, aes(p.value, color = Correct_SLiM_binding, fill = Correct_SLiM_binding)) +
geom_histogram(alpha = 0.3, aes(y = ..density..), bins = 50, position = "identity") +
theme_bw() +
facet_grid(network~., scales = "free_y")+
theme(legend.position = "none",
strip.text.y = element_text(angle = 0),
strip.text = element_text(size = 14),
legend.title = element_text(size = 14),
legend.text = element_text(size = 14),
axis.title = element_text(size = 14),
axis.text = element_text(size = 14),
text = element_text(size = 14))
### benchmark end
prior = mcol[p.value<1,mean(1-p.value)]
prior
updatePvalue = function(p.value1, p.value2, N = 1, incr = 1e-10){
N = max(length(p.value1), length(p.value2))
# add/substract small value to p-value to avoid 0 likelihood
p.value2[p.value2 == 1] = 1 - incr
p.value2[p.value2 == 0] = incr
# calculate likelihood
likelihood_of_2_1 = dbeta(p.value2, p.value1 * N, (1 - p.value1) * N)
likelihood_of_2_not_1 = dbeta(p.value2, (1 - p.value1) * N, p.value1 * N)
p.value = sum(likelihood_of_2_1 * p.value1) /
(sum(likelihood_of_2_1 * p.value1) +
sum(likelihood_of_2_not_1 * (1 - p.value1)))
p.value
}
probs = c(prior, 0.999, 0, 0.979)
b = updatePvalue(probs[1], probs[2:4], N = 1)
b
up2 = updatePvalue(up1$p.value, probs[3], N = up1$N)
up3 = updatePvalue(up2$p.value, probs[4], N = up2$N)
up3
combn(c(0.999, 0, 0.979),2)
fit = beta.mle(c(0.999, 0.998, 0.979), tol = 1e-09)
hist(rbeta(10000, 0.4999 + 100, 0.4999 + 100))
mean(fit$param[1] / (fit$param[1] + fit$param[2]))
### draft end
pkgdown::build_site()
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.