You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
137 lines
6.7 KiB
137 lines
6.7 KiB
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)
|
|
}
|