Dear all,
I'm trying to add P.adj values from dunnTest to the ggplot. I tried stat_pvalue_manual() but it doesn't seem to work.
please see example below.
I appriciate your help,
BW,
Eteri
setwd("/Users/Etuna/Dropbox/Etuna PHD/Etuna/17q21 paper")
library(tidyverse)
#> Warning: package 'tidyverse' was built under R version 3.5.2
#> Warning: package 'ggplot2' was built under R version 3.5.2
#> Warning: package 'tibble' was built under R version 3.5.2
#> Warning: package 'tidyr' was built under R version 3.5.2
#> Warning: package 'purrr' was built under R version 3.5.2
#> Warning: package 'dplyr' was built under R version 3.5.2
#> Warning: package 'stringr' was built under R version 3.5.2
#> Warning: package 'forcats' was built under R version 3.5.2
library(psych)
#>
#> Attaching package: 'psych'
#> The following objects are masked from 'package:ggplot2':
#>
#> %+%, alpha
library(readxl)
#> Warning: package 'readxl' was built under R version 3.5.2
library(FSA)
#> ## FSA v0.8.22. See citation('FSA') if used in publication.
#> ## Run fishR() for related website and fishR('IFAR') for related book.
#>
#> Attaching package: 'FSA'
#> The following object is masked from 'package:psych':
#>
#> headtail
library(reprex)
#> Warning: package 'reprex' was built under R version 3.5.2
age16_rv_combined<-read.csv(file="age16_RV_Rawdata.csv", header = TRUE, sep=",")
A17q21SNPs_selected<-read.csv(file="17q21SNPs_selected.csv", header = TRUE, sep=",")
age16_RV_SNP_Rawdata<-merge.data.frame(age16_rv_combined,A17q21SNPs_selected,by="number")
Res<-matrix(0,6,5)
for (j in 8:12){
for(i in 2:7){ ##change index according to the variables (cytokines) in the dataframe of reference
x<-age16_RV_SNP_Rawdata[,i]
y<-age16_RV_SNP_Rawdata[,j]
Res[i-1,j-7]<-kruskal.test(x ~ y, data = age16_RV_SNP_Rawdata)$p.value ##test on SNPs: change SNP
}
}
adj_pv<-p.adjust(Res, method = "BH") #multiple testing adjustment
as.data.frame(cbind(colnames(age16_RV_SNP_Rawdata[,2:7]), as.numeric(Res),adj_pv)) #print the results
#> V1 V2 adj_pv
#> 1 IFN_Alpha_RV16 0.206968022022279 0.344946703370466
#> 2 IFN_beta_RV16 0.066593189156671 0.285399382100019
#> 3 CXCL10IP10_RV16 0.536832106426285 0.619421661261098
#> 4 IFN_Alpha_RV1B 0.150520597767185 0.318732203929838
#> 5 IFN_beta_RV1B 0.00561800792956939 0.0421350594717704
#> 6 CXCL10IP10_RV1B 0.740672521398965 0.740672521398965
#> 7 IFN_Alpha_RV16 0.236249864163047 0.373026101310074
#> 8 IFN_beta_RV16 0.082378788461456 0.287711461646697
#> 9 CXCL10IP10_RV16 0.417456427305632 0.521820534132041
#> 10 IFN_Alpha_RV1B 0.143279199996726 0.318732203929838
#> 11 IFN_beta_RV1B 0.00762956551620507 0.0457773930972304
#> 12 CXCL10IP10_RV1B 0.602449506445157 0.645481614048383
#> 13 IFN_Alpha_RV16 0.180442233246166 0.318732203929838
#> 14 IFN_beta_RV16 0.0486374093867026 0.243187046933513
#> 15 CXCL10IP10_RV16 0.444764114055887 0.533716936867065
#> 16 IFN_Alpha_RV1B 0.119701680822495 0.299254202056236
#> 17 IFN_beta_RV1B 0.00296086471821577 0.0421350594717704
#> 18 CXCL10IP10_RV1B 0.592849337496398 0.645481614048383
#> 19 IFN_Alpha_RV16 0.301290137273629 0.43041448181947
#> 20 IFN_beta_RV16 0.0959038205488989 0.287711461646697
#> 21 CXCL10IP10_RV16 0.394245728398496 0.514233558780647
#> 22 IFN_Alpha_RV1B 0.180614915560242 0.318732203929838
#> 23 IFN_beta_RV1B 0.00546512278214013 0.0421350594717704
#> 24 CXCL10IP10_RV1B 0.659235659146462 0.681967923254961
#> 25 IFN_Alpha_RV16 0.177109477698654 0.318732203929838
#> 26 IFN_beta_RV16 0.0890734929594249 0.287711461646697
#> 27 CXCL10IP10_RV16 0.363721374266969 0.49598369218223
#> 28 IFN_Alpha_RV1B 0.112715672629601 0.299254202056236
#> 29 IFN_beta_RV1B 0.00032481756114297 0.00974452683428912
#> 30 CXCL10IP10_RV1B 0.257614331287415 0.386421496931123
fsPT = dunnTest(IFN_beta_RV1B ~as.factor(rs7216389_T),data=age16_RV_SNP_Rawdata,method="bh") ##post-hoc test on SNPs: change SNP
#> Warning: Some rows deleted from 'x' and 'g' because missing data.
fsPT
#> Dunn (1964) Kruskal-Wallis multiple comparison
#> p-values adjusted with the Benjamini-Hochberg method.
#> Comparison Z P.unadj P.adj
#> 1 0 - 1 1.945059 0.051767836 0.077651754
#> 2 0 - 2 3.217925 0.001291214 0.003873642
#> 3 1 - 2 1.702576 0.088647417 0.088647417
snp_rs7216389<-as.factor(age16_RV_SNP_Rawdata$rs7216389_T)
levels(snp_rs7216389)<-c("CC","CT","TT")
p2<-ggplot(aes(x=as.factor(rs7216389_T),y=IFN_beta_RV1B),data = age16_RV_SNP_Rawdata)+theme(axis.text.y=element_text(size = 12),axis.text.x = element_text(size = 12))+geom_boxplot(outlier.size = -1)+xlab("rs7216389 genotype in" ~italic(GSDMB))+ylab("RV-A1 induced IFN-β\n(pg/mL)")+
theme(axis.title.y = element_text(size = rel(1.4), angle = 90))+ theme(axis.title.x = element_text(size = rel(1.4), angle = 00))+
scale_x_discrete(labels=c("0" = "CC", "1" = "CT", "2" = "TT"))
A2<-p2+geom_jitter(position = position_jitter(0.15),aes(color=snp_rs7216389),show.legend=FALSE)
A2
#> Warning: Removed 35 rows containing non-finite values (stat_boxplot).
#> Warning: Removed 35 rows containing missing values (geom_point).
Created on 2020-04-03 by the reprex package (v0.3.0)