552 lines
22 KiB
R
Executable file
552 lines
22 KiB
R
Executable file
#!/usr/bin/env Rscript
|
|
#########################################################
|
|
# TASK: KS test for PS/DUET lineage distributions
|
|
#=======================================================================
|
|
#!/usr/bin/env Rscript
|
|
#source("~/git/LSHTM_analysis/config/alr.R")
|
|
source("~/git/LSHTM_analysis/config/embb.R")
|
|
#source("~/git/LSHTM_analysis/config/katg.R")
|
|
#source("~/git/LSHTM_analysis/config/gid.R")
|
|
#source("~/git/LSHTM_analysis/config/pnca.R")
|
|
#source("~/git/LSHTM_analysis/config/rpob.R")
|
|
|
|
# get plottting dfs
|
|
source("~/git/LSHTM_analysis/scripts/plotting/get_plotting_dfs.R")
|
|
|
|
source("~/git/LSHTM_analysis/scripts/plotting/plotting_colnames.R")
|
|
#=============
|
|
# Output
|
|
#=============
|
|
outdir_images = paste0("~/git/Writing/thesis/images/results/", tolower(gene), "/")
|
|
outdir_stats = paste0(outdir_images,"stats/")
|
|
|
|
# ks test by lineage
|
|
#ks_lineage = paste0(outdir, "/KS_lineage_all_muts.csv")
|
|
|
|
###########################
|
|
# Data for stats
|
|
# you need df2 or df2_comp
|
|
# since this is one-many relationship
|
|
###########################
|
|
|
|
# REASSIGNMENT
|
|
df2 = merged_df2[, colnames(merged_df2)%in%plotting_cols]
|
|
|
|
# quick checks
|
|
colnames(df2)
|
|
str(df2)
|
|
|
|
########################################################################
|
|
table(df2$lineage); str(df2$lineage)
|
|
|
|
# subset only lineages1-4
|
|
sel_lineages = c("L1"
|
|
, "L2"
|
|
, "L3"
|
|
, "L4")
|
|
|
|
# subset selected lineages
|
|
df_lin = subset(df2, subset = lineage %in% sel_lineages)
|
|
|
|
table(df_lin$lineage)
|
|
table(df_lin$sensitivity)
|
|
|
|
table(df_lin$lineage, df_lin$sensitivity)
|
|
|
|
#ensure lineage is a factor
|
|
#str(df_lin$lineage); str(df_lin$lineage_labels)
|
|
#df_lin$lineage = as.factor(df_lin$lineage)
|
|
#df_lin$lineage_labels = as.factor(df_lin$lineage)
|
|
table(df_lin$lineage); table(df_lin$lineage_labels)]
|
|
|
|
#==============================
|
|
# Stats for average stability
|
|
#=============================
|
|
# individual: CHECKS
|
|
lin1 = df_lin[df_lin$lineage == "L1",]$avg_stability_scaled
|
|
lin2 = df_lin[df_lin$lineage == "L2",]$avg_stability_scaled
|
|
lin3 = df_lin[df_lin$lineage == "L3",]$avg_stability_scaled
|
|
lin4 = df_lin[df_lin$lineage == "L4",]$avg_stability_scaled
|
|
|
|
ks.test(lin1, lin4)
|
|
|
|
ks.test(df_lin$avg_stability_scaled[df_lin$lineage == "L1"]
|
|
, df_lin$avg_stability_scaled[df_lin$lineage == "L4"])
|
|
|
|
#=======================================================================
|
|
my_lineages = levels(factor(df_lin$lineage)); my_lineages
|
|
#=======================================================================
|
|
# Loop
|
|
#0 : < 2.2e-16
|
|
#=====================
|
|
# Lineage 1 comparisons
|
|
#=====================
|
|
my_lin1 = "L1"
|
|
#my_lineages_comp_l1 = c("L2", "L3", "L4")
|
|
my_lineages_comp_l1 = my_lineages[-match(my_lin1, my_lineages)]
|
|
|
|
ks_df_l1 = data.frame()
|
|
|
|
for (i in my_lineages_comp_l1){
|
|
cat(i)
|
|
l1_df = data.frame(comparison = NA, method = NA, ks_statistic = NA, ks_pvalue = NA
|
|
, n_samples = NA
|
|
, n_samples_total = NA)
|
|
|
|
lineage_comp = paste0(my_lin1, " vs ", i)
|
|
n_samples_lin = nrow(df_lin[df_lin$lineage == my_lin1,])
|
|
n_samples_i = nrow(df_lin[df_lin$lineage == i,])
|
|
n_samples_all = paste0(my_lin1, "(", n_samples_lin, ")"
|
|
, ", "
|
|
, i, "(", n_samples_i, ")")
|
|
|
|
n_samples_total = n_samples_lin + n_samples_i
|
|
|
|
ks_method = ks.test(df_lin$avg_stability_scaled[df_lin$lineage == my_lin1]
|
|
, df_lin$avg_stability_scaled[df_lin$lineage == i])$method
|
|
|
|
ks_statistic = ks.test(df_lin$avg_stability_scaled[df_lin$lineage == my_lin1]
|
|
, df_lin$avg_stability_scaled[df_lin$lineage == i])$statistic
|
|
|
|
ks_pvalue = ks.test(df_lin$avg_stability_scaled[df_lin$lineage == my_lin1]
|
|
, df_lin$avg_stability_scaled[df_lin$lineage == i])$p.value
|
|
|
|
# print(c(lineage_comp, ks_method, ks_statistic[[1]], ks_pval))
|
|
l1_df$comparison = lineage_comp
|
|
l1_df$method = ks_method
|
|
l1_df$ks_statistic = ks_statistic[[1]]
|
|
l1_df$ks_pvalue = ks_pvalue
|
|
l1_df$n_samples = n_samples_all
|
|
l1_df$n_samples_total= n_samples_total
|
|
|
|
print(l1_df)
|
|
ks_df_l1 = rbind(ks_df_l1,l1_df)
|
|
|
|
}
|
|
|
|
ks_df_l1
|
|
|
|
# adjusted p-value: bonferroni
|
|
ks_df_l1$p_adj_bonferroni = p.adjust(ks_df_l1$ks_pvalue, method = "bonferroni")
|
|
ks_df_l1$signif_bon = ks_df_l1$p_adj_bonferroni
|
|
ks_df_l1 = dplyr::mutate(ks_df_l1
|
|
, signif_bon = case_when(signif_bon == 0.05 ~ "."
|
|
, signif_bon <=0.0001 ~ '****'
|
|
, signif_bon <=0.001 ~ '***'
|
|
, signif_bon <=0.01 ~ '**'
|
|
, signif_bon <0.05 ~ '*'
|
|
, TRUE ~ 'ns'))
|
|
|
|
# adjusted p-value:fdr
|
|
ks_df_l1$p_adj_fdr = p.adjust(ks_df_l1$ks_pvalue, method = "fdr")
|
|
ks_df_l1$signif_fdr = ks_df_l1$p_adj_fdr
|
|
ks_df_l1 = dplyr::mutate(ks_df_l1
|
|
, signif_fdr = case_when(signif_fdr == 0.05 ~ "."
|
|
, signif_fdr <=0.0001 ~ '****'
|
|
, signif_fdr <=0.001 ~ '***'
|
|
, signif_fdr <=0.01 ~ '**'
|
|
, signif_fdr <0.05 ~ '*'
|
|
, TRUE ~ 'ns'))
|
|
|
|
ks_df_l1
|
|
|
|
#####################################################################
|
|
|
|
#=====================
|
|
# Lineage 2 comparisons
|
|
#=====================
|
|
my_lin2 = "L2"
|
|
#my_lineages_comp_l2 = c("L1", lineage3", "L4")
|
|
my_lineages_comp_l2 = my_lineages[-match(my_lin2, my_lineages)]
|
|
|
|
ks_df_l2 = data.frame()
|
|
|
|
for (i in my_lineages_comp_l2){
|
|
cat(i)
|
|
l2_df = data.frame(comparison = NA, method = NA, ks_statistic = NA, ks_pvalue = NA
|
|
, n_samples = NA
|
|
, n_samples_total = NA)
|
|
|
|
lineage_comp = paste0(my_lin2, " vs ", i)
|
|
n_samples_lin = nrow(df_lin[df_lin$lineage == my_lin2,])
|
|
n_samples_i = nrow(df_lin[df_lin$lineage == i,])
|
|
n_samples_all = paste0(my_lin2, "(", n_samples_lin, ")"
|
|
, ", "
|
|
, i, "(", n_samples_i, ")")
|
|
|
|
n_samples_total = n_samples_lin + n_samples_i
|
|
|
|
ks_method = ks.test(df_lin$avg_stability_scaled[df_lin$lineage == my_lin2]
|
|
, df_lin$avg_stability_scaled[df_lin$lineage == i])$method
|
|
|
|
ks_statistic = ks.test(df_lin$avg_stability_scaled[df_lin$lineage == my_lin2]
|
|
, df_lin$avg_stability_scaled[df_lin$lineage == i])$statistic
|
|
|
|
ks_pvalue = ks.test(df_lin$avg_stability_scaled[df_lin$lineage == my_lin2]
|
|
, df_lin$avg_stability_scaled[df_lin$lineage == i])$p.value
|
|
|
|
# print(c(lineage_comp, ks_method, ks_statistic[[1]], ks_pval))
|
|
l2_df$comparison = lineage_comp
|
|
l2_df$method = ks_method
|
|
l2_df$ks_statistic = ks_statistic[[1]]
|
|
l2_df$ks_pvalue = ks_pvalue
|
|
l2_df$n_samples = n_samples_all
|
|
l2_df$n_samples_total = n_samples_total
|
|
|
|
print(l2_df)
|
|
ks_df_l2 = rbind(ks_df_l2, l2_df)
|
|
|
|
}
|
|
|
|
|
|
# adjusted p-value: bonferroni
|
|
ks_df_l2$p_adj_bonferroni = p.adjust(ks_df_l2$ks_pvalue, method = "bonferroni")
|
|
ks_df_l2$signif_bon = ks_df_l2$p_adj_bonferroni
|
|
ks_df_l2 = dplyr::mutate(ks_df_l2
|
|
, signif_bon = case_when(signif_bon == 0.05 ~ "."
|
|
, signif_bon <=0.0001 ~ '****'
|
|
, signif_bon <=0.001 ~ '***'
|
|
, signif_bon <=0.01 ~ '**'
|
|
, signif_bon <0.05 ~ '*'
|
|
, TRUE ~ 'ns'))
|
|
|
|
# adjusted p-value:fdr
|
|
ks_df_l2$p_adj_fdr = p.adjust(ks_df_l2$ks_pvalue, method = "fdr")
|
|
ks_df_l2$signif_fdr = ks_df_l2$p_adj_fdr
|
|
ks_df_l2 = dplyr::mutate(ks_df_l2
|
|
, signif_fdr = case_when(signif_fdr == 0.05 ~ "."
|
|
, signif_fdr <=0.0001 ~ '****'
|
|
, signif_fdr <=0.001 ~ '***'
|
|
, signif_fdr <=0.01 ~ '**'
|
|
, signif_fdr <0.05 ~ '*'
|
|
, TRUE ~ 'ns'))
|
|
|
|
ks_df_l2
|
|
|
|
#=====================
|
|
# Lineage 3 comparisons
|
|
#=====================
|
|
my_lin3 = "L3"
|
|
#my_lineages_comp_l3 = c("L1", lineage2", "L4")
|
|
my_lineages_comp_l3 = my_lineages[-match(my_lin3, my_lineages)]
|
|
|
|
ks_df_l3 = data.frame()
|
|
|
|
for (i in my_lineages_comp_l3){
|
|
cat(i)
|
|
l3_df = data.frame(comparison = NA, method = NA, ks_statistic = NA, ks_pvalue = NA
|
|
, n_samples = NA
|
|
, n_samples_total = NA)
|
|
|
|
lineage_comp = paste0(my_lin3, " vs ", i)
|
|
n_samples_lin = nrow(df_lin[df_lin$lineage == my_lin3,])
|
|
n_samples_i = nrow(df_lin[df_lin$lineage == i,])
|
|
n_samples_all = paste0(my_lin3, "(", n_samples_lin, ")"
|
|
, ", "
|
|
, i, "(", n_samples_i, ")")
|
|
n_samples_total = n_samples_lin + n_samples_i
|
|
|
|
ks_method = ks.test(df_lin$avg_stability_scaled[df_lin$lineage == my_lin3]
|
|
, df_lin$avg_stability_scaled[df_lin$lineage == i])$method
|
|
|
|
ks_statistic = ks.test(df_lin$avg_stability_scaled[df_lin$lineage == my_lin3]
|
|
, df_lin$avg_stability_scaled[df_lin$lineage == i])$statistic
|
|
|
|
ks_pvalue = ks.test(df_lin$avg_stability_scaled[df_lin$lineage == my_lin3]
|
|
, df_lin$avg_stability_scaled[df_lin$lineage == i])$p.value
|
|
|
|
# print(c(lineage_comp, ks_method, ks_statistic[[1]], ks_pval))
|
|
l3_df$comparison = lineage_comp
|
|
l3_df$method = ks_method
|
|
l3_df$ks_statistic = ks_statistic[[1]]
|
|
l3_df$ks_pvalue = ks_pvalue
|
|
l3_df$n_samples = n_samples_all
|
|
l3_df$n_samples_total = n_samples_total
|
|
|
|
print(l3_df)
|
|
ks_df_l3 = rbind(ks_df_l3, l3_df)
|
|
|
|
}
|
|
|
|
# adjusted p-value: bonferroni
|
|
ks_df_l3$p_adj_bonferroni = p.adjust(ks_df_l3$ks_pvalue, method = "bonferroni")
|
|
ks_df_l3$signif_bon = ks_df_l3$p_adj_bonferroni
|
|
ks_df_l3 = dplyr::mutate(ks_df_l3
|
|
, signif_bon = case_when(signif_bon == 0.05 ~ "."
|
|
, signif_bon <=0.0001 ~ '****'
|
|
, signif_bon <=0.001 ~ '***'
|
|
, signif_bon <=0.01 ~ '**'
|
|
, signif_bon <0.05 ~ '*'
|
|
, TRUE ~ 'ns'))
|
|
|
|
# adjusted p-value:fdr
|
|
ks_df_l3$p_adj_fdr = p.adjust(ks_df_l3$ks_pvalue, method = "fdr")
|
|
ks_df_l3$signif_fdr = ks_df_l3$p_adj_fdr
|
|
ks_df_l3 = dplyr::mutate(ks_df_l3
|
|
, signif_fdr = case_when(signif_fdr == 0.05 ~ "."
|
|
, signif_fdr <=0.0001 ~ '****'
|
|
, signif_fdr <=0.001 ~ '***'
|
|
, signif_fdr <=0.01 ~ '**'
|
|
, signif_fdr <0.05 ~ '*'
|
|
, TRUE ~ 'ns'))
|
|
|
|
ks_df_l3
|
|
|
|
####################################################################
|
|
# combine all 3 ks_dfs
|
|
n_dfs = 3
|
|
if ( all.equal(nrow(ks_df_l1), nrow(ks_df_l2), nrow(ks_df_l3)) &&
|
|
all.equal(ncol(ks_df_l1), ncol(ks_df_l2), ncol(ks_df_l3)) ){
|
|
cat("\nPASS: Calculating expected rows and cols for sanity checks on combined_dfs")
|
|
expected_rows = nrow(ks_df_l1) * n_dfs
|
|
expected_cols = ncol(ks_df_l1)
|
|
ks_df_combined = rbind(ks_df_l1, ks_df_l2, ks_df_l3)
|
|
if ( nrow(ks_df_combined) == expected_rows && ncol(ks_df_combined) == expected_cols ){
|
|
cat("\nPASS: combined df successfully created"
|
|
, "\nnrow combined_df:", nrow(ks_df_combined)
|
|
, "\nncol combined_df:", ncol(ks_df_combined))
|
|
}
|
|
else{
|
|
cat("\nFAIL: Dim mismatch"
|
|
, "\nExpected rows:", expected_rows
|
|
, "\nGot:", nrow(ks_df_combined)
|
|
, "\nExpected cols:", expected_cols
|
|
, "\nGot:", ncol(ks_df_combined))
|
|
}
|
|
}else{
|
|
cat("\nFAIL: Could not generate expected_rows and/or expected_cols"
|
|
, "\nCheck hardcoded value of n_dfs")
|
|
}
|
|
|
|
#----------------------------
|
|
# ADD extra cols: formatting
|
|
#----------------------------
|
|
# adding pvalue_signif
|
|
ks_df_combined$pvalue_signif = ks_df_combined$ks_pvalue
|
|
str(ks_df_combined$pvalue_signif)
|
|
|
|
ks_df_combined = dplyr::mutate(ks_df_combined
|
|
, pvalue_signif = case_when(pvalue_signif == 0.05 ~ "."
|
|
, pvalue_signif <=0.0001 ~ '****'
|
|
, pvalue_signif <=0.001 ~ '***'
|
|
, pvalue_signif <=0.01 ~ '**'
|
|
, pvalue_signif <0.05 ~ '*'
|
|
, TRUE ~ 'ns'))
|
|
|
|
# Remove duplicates
|
|
rows_to_remove = c("L2 vs L1", "L3 vs L1", "L3 vs L2")
|
|
|
|
ks_df_combined_f = ks_df_combined[-match(rows_to_remove, ks_df_combined$comparison),]
|
|
#################################################################################
|
|
#================================
|
|
# R vs S distribution: Overall
|
|
#=================================
|
|
df_lin_R = df_lin[df_lin$sensitivity == "R",]
|
|
df_lin_S = df_lin[df_lin$sensitivity == "S",]
|
|
|
|
overall_RS_df = data.frame(comparison = NA, method = NA, ks_statistic = NA, ks_pvalue = NA
|
|
, n_samples = NA
|
|
, n_samples_total = NA)
|
|
|
|
comp_type = "overall R vs S distributions"
|
|
|
|
n_samples_R = nrow(df_lin_R)
|
|
n_samples_S = nrow(df_lin_S)
|
|
n_samples_all = paste0("R(n=", n_samples_R, ")" , " ", "S(n=", n_samples_S, ")")
|
|
n_samples_total = n_samples_R+n_samples_S
|
|
|
|
ks.test(df_lin_R$avg_stability_scaled
|
|
, df_lin_S$avg_stability_scaled)
|
|
|
|
ks_method = ks.test(df_lin_R$avg_stability_scaled
|
|
, df_lin_S$avg_stability_scaled)$method
|
|
|
|
ks_statistic = ks.test(df_lin_R$avg_stability_scaled
|
|
, df_lin_S$avg_stability_scaled)$statistic
|
|
|
|
ks_pvalue = ks.test(df_lin_R$avg_stability_scaled
|
|
, df_lin_S$avg_stability_scaled)$p.value
|
|
|
|
overall_RS_df$comparison = comp_type
|
|
overall_RS_df$method = ks_method
|
|
overall_RS_df$ks_statistic = ks_statistic
|
|
overall_RS_df$ks_pvalue = ks_pvalue
|
|
overall_RS_df$n_samples = n_samples_all
|
|
overall_RS_df$n_samples_total= n_samples_total
|
|
|
|
#----------------------------
|
|
# ADD extra cols
|
|
#----------------------------
|
|
# adjusted p-value: bonferroni
|
|
overall_RS_df$p_adj_bonferroni = p.adjust(overall_RS_df$ks_pvalue, method = "bonferroni")
|
|
overall_RS_df$signif_bon = overall_RS_df$p_adj_bonferroni
|
|
overall_RS_df = dplyr::mutate(overall_RS_df
|
|
, signif_bon = case_when(signif_bon == 0.05 ~ "."
|
|
, signif_bon <=0.0001 ~ '****'
|
|
, signif_bon <=0.001 ~ '***'
|
|
, signif_bon <=0.01 ~ '**'
|
|
, signif_bon <0.05 ~ '*'
|
|
, TRUE ~ 'ns'))
|
|
|
|
# adjusted p-value:fdr
|
|
overall_RS_df$p_adj_fdr = p.adjust(overall_RS_df$ks_pvalue, method = "fdr")
|
|
overall_RS_df$signif_fdr = overall_RS_df$p_adj_fdr
|
|
overall_RS_df = dplyr::mutate(overall_RS_df
|
|
, signif_fdr = case_when(signif_fdr == 0.05 ~ "."
|
|
, signif_fdr <=0.0001 ~ '****'
|
|
, signif_fdr <=0.001 ~ '***'
|
|
, signif_fdr <=0.01 ~ '**'
|
|
, signif_fdr <0.05 ~ '*'
|
|
, TRUE ~ 'ns'))
|
|
# unadjusted p-values
|
|
overall_RS_df$pvalue_signif = overall_RS_df$ks_pvalue
|
|
overall_RS_df = dplyr::mutate(overall_RS_df
|
|
, pvalue_signif = case_when(pvalue_signif == 0.05 ~ "."
|
|
, pvalue_signif <=0.0001 ~ '****'
|
|
, pvalue_signif <=0.001 ~ '***'
|
|
, pvalue_signif <=0.01 ~ '**'
|
|
, pvalue_signif <0.05 ~ '*'
|
|
, TRUE ~ 'ns'))
|
|
#####################################################################
|
|
if (all(colnames(ks_df_combined_f) == colnames(overall_RS_df))){
|
|
|
|
cat("\nPASS:combining KS test results")
|
|
}else{
|
|
stop("\nAbort: Cannot combine results for KS test")
|
|
}
|
|
|
|
ks_df_combined_f2 = rbind(ks_df_combined_f, overall_RS_df)
|
|
|
|
# ADD extra cols
|
|
ks_df_combined_f2$ks_comp_type = "between_lineages"
|
|
ks_df_combined_f2$gene_name = tolower(gene)
|
|
|
|
ks_df_combined_f2
|
|
|
|
###########################################################################
|
|
#=================================
|
|
# Within lineage R vs S: MANUAL
|
|
#=================================
|
|
lin1 = df_lin[df_lin$lineage == my_lin1,]
|
|
ks.test(lin1$avg_stability_scaled[lin1$sensitivity == "R"]
|
|
, lin1$avg_stability_scaled[lin1$sensitivity == "S"])
|
|
|
|
lin2 = df_lin[df_lin$lineage == my_lin2,]
|
|
ks.test(lin2$avg_stability_scaled[lin2$sensitivity == "R"]
|
|
, lin2$avg_stability_scaled[lin2$sensitivity == "S"])
|
|
|
|
|
|
lin3 = df_lin[df_lin$lineage == "L3",]
|
|
ks.test(lin3$avg_stability_scaled[lin3$sensitivity == "R"]
|
|
, lin3$avg_stability_scaled[lin3$sensitivity == "S"])
|
|
|
|
lin4 = df_lin[df_lin$lineage == "L4",]
|
|
ks.test(lin4$avg_stability_scaled[lin4$sensitivity == "R"]
|
|
, lin4$avg_stability_scaled[lin4$sensitivity == "S"])
|
|
###################################################################
|
|
#=======================
|
|
# Within lineage R vs S: LOOP
|
|
#=======================
|
|
all_within_lin_df = data.frame()
|
|
|
|
for (i in c(my_lin1
|
|
, my_lin2
|
|
, my_lin3
|
|
)){
|
|
#cat(i)
|
|
|
|
within_lin_df = data.frame(comparison = NA, method = NA, ks_statistic = NA, ks_pvalue = NA
|
|
, n_samples = NA
|
|
, n_samples_total = NA)
|
|
|
|
within_lineage_comp = paste0(i, ": R vs S")
|
|
my_lin_df = df_lin[df_lin$lineage == i,]
|
|
|
|
lin_R_n = nrow(my_lin_df[my_lin_df$sensitivity == "R",])
|
|
lin_S_n = nrow(my_lin_df[my_lin_df$sensitivity == "S",])
|
|
lin_total_n = lin_R_n+lin_S_n
|
|
n_samples_all = paste0("R(n=", lin_R_n, ")" , ", ", "S(n=", lin_S_n, ")")
|
|
|
|
ks_method = ks.test(my_lin_df$avg_stability_scaled[my_lin_df$sensitivity == "R"]
|
|
, my_lin_df$avg_stability_scaled[my_lin_df$sensitivity == "S"])$method
|
|
|
|
ks_statistic = ks.test(my_lin_df$avg_stability_scaled[my_lin_df$sensitivity == "R"]
|
|
, my_lin_df$avg_stability_scaled[my_lin_df$sensitivity == "S"])$statistic
|
|
|
|
ks_pvalue =ks.test(my_lin_df$avg_stability_scaled[my_lin_df$sensitivity == "R"]
|
|
, my_lin_df$avg_stability_scaled[my_lin_df$sensitivity == "S"])$p.value
|
|
|
|
# print(c(lineage_comp, ks_method, ks_statistic[[1]], ks_pval))
|
|
within_lin_df$comparison = within_lineage_comp
|
|
within_lin_df$method = ks_method
|
|
within_lin_df$ks_statistic = ks_statistic[[1]]
|
|
within_lin_df$ks_pvalue = ks_pvalue
|
|
within_lin_df$n_samples = n_samples_all
|
|
within_lin_df$n_samples_total=lin_total_n
|
|
|
|
print(my_lin_df)
|
|
all_within_lin_df = rbind(all_within_lin_df, within_lin_df)
|
|
|
|
}
|
|
all_within_lin_df
|
|
|
|
#----------------------------
|
|
# ADD extra cols: formatting
|
|
#----------------------------
|
|
# adjusted p-value: bonferroni
|
|
all_within_lin_df$p_adj_bonferroni = p.adjust(all_within_lin_df$ks_pvalue, method = "bonferroni")
|
|
all_within_lin_df$signif_bon = all_within_lin_df$p_adj_bonferroni
|
|
all_within_lin_df = dplyr::mutate(all_within_lin_df
|
|
, signif_bon = case_when(signif_bon == 0.05 ~ "."
|
|
, signif_bon <=0.0001 ~ '****'
|
|
, signif_bon <=0.001 ~ '***'
|
|
, signif_bon <=0.01 ~ '**'
|
|
, signif_bon <0.05 ~ '*'
|
|
, TRUE ~ 'ns'))
|
|
|
|
# adjusted p-value:fdr
|
|
all_within_lin_df$p_adj_fdr = p.adjust(all_within_lin_df$ks_pvalue, method = "fdr")
|
|
all_within_lin_df$signif_fdr = all_within_lin_df$p_adj_fdr
|
|
all_within_lin_df = dplyr::mutate(all_within_lin_df
|
|
, signif_fdr = case_when(signif_fdr == 0.05 ~ "."
|
|
, signif_fdr <=0.0001 ~ '****'
|
|
, signif_fdr <=0.001 ~ '***'
|
|
, signif_fdr <=0.01 ~ '**'
|
|
, signif_fdr <0.05 ~ '*'
|
|
, TRUE ~ 'ns'))
|
|
|
|
# unadjusted p-value
|
|
all_within_lin_df$pvalue_signif = all_within_lin_df$ks_pvalue
|
|
str(all_within_lin_df$pvalue_signif)
|
|
all_within_lin_df = dplyr::mutate(all_within_lin_df
|
|
, pvalue_signif = case_when(pvalue_signif == 0.05 ~ "."
|
|
, pvalue_signif <=0.0001 ~ '****'
|
|
, pvalue_signif <=0.001 ~ '***'
|
|
, pvalue_signif <=0.01 ~ '**'
|
|
, pvalue_signif <0.05 ~ '*'
|
|
, TRUE ~ 'ns'))
|
|
|
|
# ADD info cols
|
|
all_within_lin_df$ks_comp_type = "within_lineages"
|
|
all_within_lin_df$gene_name = tolower(gene)
|
|
|
|
##################################################################
|
|
|
|
if (all(colnames(ks_df_combined_f2) == colnames(all_within_lin_df))){
|
|
|
|
cat("\nPASS:combining KS test results")
|
|
}else{
|
|
stop("\nAbort: Cannot combine results for KS test")
|
|
}
|
|
|
|
ks_df_combined_all = rbind(ks_df_combined_f2, all_within_lin_df)
|
|
|
|
#--------------------
|
|
# write output file: KS test within grpup
|
|
#----------------------
|
|
Out_ks_all = paste0(outdir_stats
|
|
, tolower(gene)
|
|
, "_ks_lineage_all_comp.csv")
|
|
|
|
cat("Output of KS test all comparisons:", Out_ks_all )
|
|
write.csv(ks_df_combined_all, Out_ks_all, row.names = FALSE)
|