|
|
- library(dplyr)
- library(dunn.test)
- library(tibble)
- library(reshape2)
- library(ggplot2)
-
- corr_multi<-function(table, genes1, genes2, font_size=5, ruta="./", exportar=F)
- {
- cor.list<-list()
- cont<-1
- for (gene1 in genes1){
- for (gene2 in genes2){
- print(gene1)
- print(gene2)
- cor_temp<-cor.test(table[,gene1], table[,gene2])
- if (!is.na(cor_temp$p.value)){
- pval_plot<-paste("cor =",format(cor_temp$estimate, digits=2),if(cor_temp$p.value < 2.2e-16){" , p < 2.2e-16"}else{paste(", p=",format(cor_temp$p.value,digits=3))})
- }
- cor.list[[cont]]<-ggplotGrob(ggplot(data=table, aes(x=table[,gene1], y=table[,gene2]))+
- geom_point()+
- geom_smooth(method="lm")+
- geom_text(aes(x=min(table[,gene1], na.rm=T), y=(max(table[,gene2],na.rm = T)+max(table[,gene2],na.rm = T)*0.2), label=pval_plot), hjust="inward", family="serif", size=font_size)+
- labs(x=gene1, y=gene2)+
- theme_bw())
- if (exportar == T){
- png(paste0(ruta,gene1,"vs",gene2,".png"), res=900, width=3000, height=3000)
- plot(cor.list[[cont]])
- dev.off()
- }
- cont<-cont+1
- }
- }
- return(cor.list)
- }
-
- multi_stats<-function(table, value.var, x, group, stat.test, adjust="default", paired=F){
- ## Requires dplyr and tibble packets
- defaults=c("dunn"="none", "ttest"="holm", "wilcox"="holm")
- funs<-c("ttest"="pairwise.t.test", "wilcox"="pairwise.wilcox.test")
- if (adjust == "default"){adjust=defaults[stat.test]}
- stat.def<-as.data.frame(matrix(nrow=0, ncol=5))
- colnames(stat.def)<-c(x, "group1", "group2", "p.adj", "p.signif")
- for (point in unique(table[,x])){
- condition<-all(table %>% filter(table[,x] == point) %>% pull(value.var) == 0) == F
- len_group<-length(unique(table %>% filter(table[,x] == point) %>% pull(group)))
- if (condition == T & !is.na(condition) & len_group > 1){
- if(stat.test == "dunn"){
- test<-dunn.test(table %>% filter(table[,x] == point) %>% pull(value.var), table %>% filter(table[,x] == point) %>% pull(group), method=adjust)
- comp<-strsplit(test$comparisons, " - ")
- stat.temp<-data.frame(matrix(unlist(comp), nrow=length(comp), byrow=T), "p.adj"=test$P.adjusted, stringsAsFactors = F, check.names = F)
- colnames(stat.temp)[1:2]<-c("group1", "group2")
- }else if (stat.test %in% names(funs)){
- test<-get(funs[stat.test])(table %>% filter(table[,x] == point) %>% pull(value.var), table %>% filter(table[,x] == point) %>% pull(group), method=adjust, paired=paired)$p.value
- stat.temp<-melt(test)
- colnames(stat.temp)<-c("group1", "group2","p.adj")
- }
-
- stat.temp["p.signif"]<-case_when(
- stat.temp$p.adj >= 0.05 ~ "ns",
- stat.temp$p.adj < 0.0001 ~ "****",
- stat.temp$p.adj < 0.001 ~ "***",
- stat.temp$p.adj < 0.01 ~ "**",
- stat.temp$p.adj < 0.05 ~ "*"
- )
- stat.temp<-stat.temp %>% add_column(x=point, .before=T)
- colnames(stat.temp)[1]<-x
- stat.def<-rbind(stat.def, stat.temp)
- }
- }
- stat.def["Method"]<-stat.test
- return(stat.def)
- }
-
- generate_labstats<-function(table_stat, table, value.var, x, group, y="max", bracket.offset=0.05, bracket.length=0.02){
- table[,group]<-as.factor(table[,group])
- table[,x]<-as.factor(table[,x])
- se<-function(x, na.rm=F) sd(x, na.rm = na.rm)/sqrt(length(x))
- if (y == "max"){
- formula<-as.formula(paste0(colnames(table_stat)[1], "~."))
- agg<-dcast(table, formula, value.var = value.var, fun.aggregate = max, na.rm=T)
- }else if (y == "mean"){
- formula<-as.formula(paste0(colnames(table_stat)[1], "~", group))
- agg<-dcast(table, formula, value.var = value.var, fun.aggregate = mean, na.rm=T)
- agg<- data.frame(x=agg[,1], "."=apply(agg[,2:ncol(agg)], 1, max, na.rm=T))
- colnames(agg)[1]<-x
- }else if (y == "mean+sd"){
- formula<-as.formula(paste0(colnames(table_stat)[1], "~", group))
- agg<- dcast(table, formula, value.var = value.var, fun.aggregate = function(x) mean(x,na.rm=T)+sd(x,na.rm=T))
- agg<- data.frame(x=agg[,1], "."=apply(agg[,2:ncol(agg)], 1, max, na.rm=T))
- colnames(agg)[1]<-x
- }else if (y == "mean+se"){
- formula<-as.formula(paste0(colnames(table_stat)[1], "~", group))
- agg<- dcast(table, formula, value.var = value.var, fun.aggregate = function(x) mean(x,na.rm=T)+se(x,na.rm=T))
- agg<- data.frame(agg[,1], "."=apply(agg[,2:ncol(agg)], 1, max, na.rm=T))
- colnames(agg)[1]<-x
- }
-
- t<-data.frame("y1"=merge(table_stat, agg ,sort=F)[,"."]+diff(range(table[value.var], na.rm = T))*bracket.offset,
- "y2"=merge(table_stat, agg ,sort=F)[,"."]+diff(range(table[value.var], na.rm = T))*bracket.offset,
- "x1"= match(table_stat[,x], unique(table[,x]))+
- 0.75*((match(table_stat$group1, levels(table[,group]))-0.5)/length(levels(table[,group]))-0.5),
- "x2"= match(table_stat[,x], unique(table[,x]))+
- 0.75*((match(table_stat$group2, levels(table[,group]))-0.5)/length(levels(table[,group]))-0.5)
- )
- for (dia in unique(table_stat[,1])){
- t[table_stat[,x] == dia,"y1"]<-seq(t[table_stat[,x] == dia,"y1"][1],
- t[table_stat[,x] == dia,"y1"][1]+diff(range(table[,value.var], na.rm = T))*0.05*(nrow(table_stat[table_stat[,x] == dia,])-1),
- by=diff(range(table[,value.var], na.rm = T))*0.05)
- t[table_stat[,x] == dia,"y2"]<-t[table_stat[,x] == dia,"y1"]
- }
-
- t_def<-as.data.frame(matrix(ncol=4, nrow=0))
- for (row in 1:nrow(t)){
- t_def<-rbind(t_def, t[row,],
- c(t[row,"y1"]-diff(range(table[,value.var], na.rm = T))*bracket.length, t[row,"y1"], t[row,"x1"], t[row,"x1"]),
- c(t[row,"y1"]-diff(range(table[,value.var], na.rm = T))*bracket.length, t[row,"y1"], t[row,"x2"], t[row,"x2"]))
- }
- t_lab<-data.frame("x"=t$x1+(t$x2-t$x1)/2, "y"=t$y1+diff(range(table[,value.var], na.rm = T))*0.005, check.names = F)
- return(list("label"=t_lab, "brackets"=t_def))
- }
-
- secfile<-function(file){
- ext<-strsplit(file, ".", fixed = T)[[1]]
- ext<-ext[length(ext)]
- num<-1
- while(file.exists(file) == T){
- if(num == 1){
- file_tmp<-strsplit(file, ".", fixed=T)[[1]]
- file<-paste0(paste(file_tmp[1:(length(file_tmp)-1)], collapse = "."),"_",num,".",file_tmp[length(file_tmp)])
- }else{
- file_tmp<-paste(strsplit(file, "_", fixed=T)[[1]][-length(strsplit(file, "_", fixed=T)[[1]])], collapse = "_")
- file<-paste0(file_tmp, "_", num,".",ext)
- }
- num<-num+1
- }
- return(file)
- }
|
