This notebook describes the haplotype analysis of rice loci identified by means of GWAS by Martina Huber, Prof. Ronald Pierink’s group.

library(tidyr)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library("ggplot2")
library("doBy")
## 
## Attaching package: 'doBy'
## The following object is masked from 'package:dplyr':
## 
##     order_by
library("reshape2")
## 
## Attaching package: 'reshape2'
## The following object is masked from 'package:tidyr':
## 
##     smiths
library(ggbeeswarm)
library(ggpubr)
## Warning: package 'ggpubr' was built under R version 4.0.2
library(cowplot)
## 
## ********************************************************
## Note: As of version 1.0.0, cowplot does not change the
##   default ggplot2 theme anymore. To recover the previous
##   behavior, execute:
##   theme_set(theme_cowplot())
## ********************************************************
## 
## Attaching package: 'cowplot'
## The following object is masked from 'package:ggpubr':
## 
##     get_legend
library(reshape)
## 
## Attaching package: 'reshape'
## The following object is masked from 'package:cowplot':
## 
##     stamp
## The following objects are masked from 'package:reshape2':
## 
##     colsplit, melt, recast
## The following object is masked from 'package:dplyr':
## 
##     rename
## The following objects are masked from 'package:tidyr':
## 
##     expand, smiths
library(multcompView)
## Warning: package 'multcompView' was built under R version 4.0.2

upload new data

decode <- read.csv("new_haplotypes/accession index_HYDRA-ID_for_Magda.csv")
head(decode)
decode <- decode[,c(2:7)]
decode$VCF_better <- gsub("@", ".", decode$VCF_code)

pheno <- read.csv("new_haplotypes/14traits_SI_norm.csv")

# make sure that the column names are the same
colnames(pheno)
##  [1] "index.pots"        "HDRA.ID"           "NSFTV.ID"         
##  [4] "IRGC.ID"           "Variety"           "Subpop"           
##  [7] "Shoot.area..cm2."  "Hull.area..cm2."   "Solidity"         
## [10] "Perimeter..cm."    "Nr.leaves"         "Nr.tillers"       
## [13] "Plant.height..cm." "Culm.height..cm."  "Leaf.length..cm." 
## [16] "DW..g."            "Leaf.angle..º."    "Tiller.angle..º." 
## [19] "Erectness..º."     "Droopiness..º."    "SUM_rank_traits"  
## [22] "SUM_norm_traits"   "norm_SI_rank"      "SI_rank"
colnames(decode)[1] <- "HDRA.ID"


pheno2 <- merge(decode, pheno, id=c("HDRA.ID", "NSFTV.ID", "IRGC.ID"))
head(pheno2)
pheno_id <- pheno2$VCF_better
unique(pheno_id)
##   [1] "NSFTV81.006dfe9b.0"  "NSFTV39.014ab0f1.0"  "NSFTV132.02cc7c6d.0"
##   [4] "NSFTV286.044639de.0" "NSFTV643.0491646d.0" "NSFTV32.068b860d.0" 
##   [7] "NSFTV358.0764178c.0" "NSFTV275.07dac217.0" "NSFTV269.092f15e7.0"
##  [10] "NSFTV256.0be20a08.0" "NSFTV6.0d125c0e.0"   "NSFTV324.0defa551.0"
##  [13] "NSFTV287.0e787735.0" "NSFTV126.0f6a67da.0" "NSFTV110.11bf5114.0"
##  [16] "NSFTV355.122a975b.0" "NSFTV150.128dd425.0" "NSFTV72.12d82364.0" 
##  [19] "NSFTV311.137628a5.0" "NSFTV235.13ae4f27.0" "NSFTV216.149fb038.0"
##  [22] "NSFTV364.14e43b02.0" "NSFTV323.1510e4aa.0" "NSFTV277.1516d75f.0"
##  [25] "NSFTV177.15e6c437.0" "NSFTV105.16463092.0" "NSFTV278.17627a92.0"
##  [28] "NSFTV145.17c4070a.0" "NSFTV391.180a155f.0" "NSFTV253.1911c363.0"
##  [31] "NSFTV163.192e24ab.0" "NSFTV228.195567bf.0" "NSFTV27.1956dd3f.0" 
##  [34] "NSFTV219.199f4455.0" "NSFTV138.1a946dc6.0" "NSFTV133.1a95985b.0"
##  [37] "NSFTV65.1ab838a2.0"  "NSFTV25.1ce7093b.0"  "NSFTV335.1e26852c.0"
##  [40] "NSFTV328.1e4f3933.0" "NSFTV396.1eb5d579.0" "NSFTV100.1f10be3d.0"
##  [43] "NSFTV321.1f6f13c3.0" "NSFTV394.1f856ac1.0" "NSFTV94.20a4c97d.0" 
##  [46] "NSFTV49.21d3f1b3.0"  "NSFTV333.22e60af9.0" "NSFTV380.238f25f1.0"
##  [49] "NSFTV36.2593e5a7.0"  "NSFTV201.263f2baf.0" "NSFTV377.2779bba9.0"
##  [52] "NSFTV121.280279b3.0" "NSFTV241.2a7724c1.0" "NSFTV9.2a9d8d47.0"  
##  [55] "NSFTV162.2b28e441.0" "NSFTV33.2b4e06fe.0"  "NSFTV306.2c583eed.0"
##  [58] "NSFTV157.2db62b89.0" "NSFTV10.2e1c9c87.0"  "NSFTV29.30646147.0" 
##  [61] "NSFTV8.30c3073f.0"   "NSFTV221.311aaf30.0" "NSFTV144.3269952e.0"
##  [64] "NSFTV85.328163e1.0"  "NSFTV26.32a15c4d.0"  "NSFTV96.32a6808e.0" 
##  [67] "NSFTV234.32ca1759.0" "NSFTV205.348853a5.0" "NSFTV214.34df4fca.0"
##  [70] "NSFTV206.37249a74.0" "NSFTV67.39dd7feb.0"  "NSFTV119.3aa51818.0"
##  [73] "NSFTV385.3b25c24f.0" "NSFTV195.3c10a4b8.0" "NSFTV329.3e58f34c.0"
##  [76] "NSFTV143.3ea144c5.0" "NSFTV318.3eb8fab5.0" "NSFTV255.3fb62843.0"
##  [79] "NSFTV153.40f343f4.0" "NSFTV651.432387a8.0" "NSFTV17.446f6c62.0" 
##  [82] "NSFTV243.44aeff79.0" "NSFTV186.4559ff3f.0" "NSFTV397.45d3c920.0"
##  [85] "NSFTV77.45ec7861.0"  "NSFTV178.480f7505.0" "NSFTV147.48b4f132.0"
##  [88] "NSFTV639.48d3bc66.0" "NSFTV273.493f2b7e.0" "NSFTV231.49649a00.0"
##  [91] "NSFTV340.4b0a7350.0" "NSFTV345.4e35b58a.0" "NSFTV386.4f4f777a.0"
##  [94] "NSFTV161.4f643bf5.0" "NSFTV239.523d5ec4.0" "NSFTV22.52c6e2ba.0" 
##  [97] "NSFTV112.531e23fa.0" "NSFTV226.536afc14.0" "NSFTV209.5384be83.0"
## [100] "NSFTV265.544f48c0.0" "NSFTV211.5497b233.0" "NSFTV179.5505a767.0"
## [103] "NSFTV5.5533f406.0"   "NSFTV74.55d3afae.0"  "NSFTV193.5632be21.0"
## [106] "NSFTV252.56f5103e.0" "NSFTV51.57d7feea.0"  "NSFTV280.58d5bdc3.0"
## [109] "NSFTV308.593def74.0" "NSFTV282.5aa7968d.0" "NSFTV117.5b144b9c.0"
## [112] "NSFTV390.5c592759.0" "NSFTV92.5df8f871.0"  "NSFTV3.5ef1be74.0"  
## [115] "NSFTV244.5efe8622.0" "NSFTV258.60c362f1.0" "NSFTV187.60e142c3.0"
## [118] "NSFTV166.61152244.0" "NSFTV53.61b7bf53.0"  "NSFTV71.61cc3858.0" 
## [121] "NSFTV54.06334433.0"  "NSFTV125.63f298ba.0" "NSFTV279.642238d9.0"
## [124] "NSFTV87.64fa0112.0"  "NSFTV158.6516a500.0" "NSFTV56.652530bb.0" 
## [127] "NSFTV13.660f0236.0"  "NSFTV189.661eaeaa.0" "NSFTV78.66d97672.0" 
## [130] "NSFTV270.6815ec6d.0" "NSFTV46.68c2ecf8.0"  "NSFTV37.69527f35.0" 
## [133] "NSFTV98.6ab77e3e.0"  "NSFTV20.6af7c9fc.0"  "NSFTV283.6b37cb3d.0"
## [136] "NSFTV116.6cedf6aa.0" "NSFTV102.6e26f4cc.0" "NSFTV360.6f068769.0"
## [139] "NSFTV240.6f263bd1.0" "NSFTV184.7119ebee.0" "NSFTV389.7126c359.0"
## [142] "NSFTV123.714ac141.0" "NSFTV118.71bd9426.0" "NSFTV57.71d5fcf4.0" 
## [145] "NSFTV164.72497775.0" "NSFTV322.72cef953.0" "NSFTV366.73b20824.0"
## [148] "NSFTV303.73c4fdfb.0" "NSFTV90.743f41c1.0"  "NSFTV315.74c9fbbc.0"
## [151] "NSFTV66.7691a376.0"  "NSFTV128.76a1efc9.0" "NSFTV40.76f56677.0" 
## [154] "NSFTV304.773e969e.0" "NSFTV23.7741d7c1.0"  "NSFTV204.778cea6e.0"
## [157] "NSFTV113.7a723d9e.0" "NSFTV107.7b7a0d82.0" "NSFTV84.7bee5b9f.0" 
## [160] "NSFTV222.7cfe4a0c.0" "NSFTV44.7d4c6c4e.0"  "NSFTV274.7dba928e.0"
## [163] "NSFTV332.7e956f22.0" "NSFTV99.7eb7c6a8.0"  "NSFTV60.7f7fb805.0" 
## [166] "NSFTV289.7fc911f7.0" "NSFTV142.806c51cc.0" "NSFTV346.80fc89ae.0"
## [169] "NSFTV198.818143d8.0" "NSFTV4.81d03b86.0"   "NSFTV305.8288d278.0"
## [172] "NSFTV371.84ad8457.0" "NSFTV347.853e318c.0" "NSFTV140.85551f9c.0"
## [175] "NSFTV281.85959164.0" "NSFTV218.85da3a70.0" "NSFTV104.8629f76c.0"
## [178] "NSFTV137.8653bbdb.0" "NSFTV301.86556239.0" "NSFTV313.872ec0b4.0"
## [181] "NSFTV249.87621330.0" "NSFTV641.88631879.0" "NSFTV250.88e1f162.0"
## [184] "NSFTV89.897c5eef.0"  "NSFTV106.8a06320f.0" "NSFTV245.8a41a10a.0"
## [187] "NSFTV350.8adbe877.0" "NSFTV103.8c76404c.0" "NSFTV181.8cd35626.0"
## [190] "NSFTV176.8cfa8fb2.0" "NSFTV45.8d6aded4.0"  "NSFTV120.8e6220e5.0"
## [193] "NSFTV167.8e941960.0" "NSFTV129.8fafd383.0" "NSFTV155.9127f236.0"
## [196] "NSFTV217.915c0033.0" "NSFTV156.93387d96.0" "NSFTV199.958ed737.0"
## [199] "NSFTV384.9764b3c8.0" "NSFTV191.97ce6839.0" "NSFTV73.97dcf87d.0" 
## [202] "NSFTV334.09869491.0" "NSFTV642.986a3c37.0" "NSFTV172.9ac72c9e.0"
## [205] "NSFTV165.9b2223bb.0" "NSFTV248.9b37187d.0" "NSFTV285.9baa653c.0"
## [208] "NSFTV58.9c099b78.0"  "NSFTV93.9f1f4614.0"  "NSFTV101.9f782e9e.0"
## [211] "NSFTV341.9f786a86.0" "NSFTV387.a08839ce.0" "NSFTV344.a0f35768.0"
## [214] "NSFTV225.a31b4a06.0" "NSFTV284.a32e0586.0" "NSFTV348.a446becd.0"
## [217] "NSFTV349.a5a6e2dd.0" "NSFTV59.a5ec1b10.0"  "NSFTV130.a796716d.0"
## [220] "NSFTV227.a7983b03.0" "NSFTV337.a7bec464.0" "NSFTV141.a8319fc6.0"
## [223] "NSFTV342.ac7a352c.0" "NSFTV267.aceb3352.0" "NSFTV638.ad1a9c05.0"
## [226] "NSFTV326.aeed4684.0" "NSFTV353.af77442b.0" "NSFTV75.b0816082.0" 
## [229] "NSFTV242.b21f7528.0" "NSFTV213.b2320902.0" "NSFTV171.b238d197.0"
## [232] "NSFTV291.b29d3bf4.0" "NSFTV339.b3a301ea.0" "NSFTV268.b49ad0db.0"
## [235] "NSFTV343.b663b3f8.0" "NSFTV271.b67fd5f6.0" "NSFTV122.b6dc1bcc.0"
## [238] "NSFTV316.b828b757.0" "NSFTV317.b956dfb5.0" "NSFTV288.bc6b0af5.0"
## [241] "NSFTV652.bc9c6f80.0" "NSFTV139.bd0d322b.0" "NSFTV370.bd7aaa87.0"
## [244] "NSFTV83.becfe767.0"  "NSFTV188.bf18cbae.0" "NSFTV190.bf7afed8.0"
## [247] "NSFTV260.bfa16661.0" "NSFTV183.c08a5ea2.0" "NSFTV79.c0936cf1.0" 
## [250] "NSFTV307.c0d6eca3.0" "NSFTV69.c14c4e03.0"  "NSFTV373.c20dfc59.0"
## [253] "NSFTV392.c4a397e5.0" "NSFTV232.c528bce0.0" "NSFTV35.c59435a8.0" 
## [256] "NSFTV220.c5bf98cc.0" "NSFTV319.c7aeb39f.0" "NSFTV263.c7d2513e.0"
## [259] "NSFTV359.c86f37d0.0" "NSFTV200.c8a73fe7.0" "NSFTV185.c97be7ae.0"
## [262] "NSFTV76.cb5e38e6.0"  "NSFTV257.cb9c18df.0" "NSFTV18.cbd6b346.0" 
## [265] "NSFTV357.cc0ef8cf.0" "NSFTV202.ccdc6d39.0" "NSFTV261.ceb798a3.0"
## [268] "NSFTV320.cf44b628.0" "NSFTV196.cff871e8.0" "NSFTV325.d0392fe8.0"
## [271] "NSFTV131.d09d62e7.0" "NSFTV19.d22c8e33.0"  "NSFTV31.d2da857d.0" 
## [274] "NSFTV331.d438bfed.0" "NSFTV207.d628ec3a.0" "NSFTV290.d864377c.0"
## [277] "NSFTV368.db737b9b.0" "NSFTV237.dc2cf39b.0" "NSFTV369.dd2bfbfb.0"
## [280] "NSFTV30.dd6e755e.0"  "NSFTV376.de696d95.0" "NSFTV247.de7ac7bf.0"
## [283] "NSFTV151.dee26171.0" "NSFTV41.dfb13733.0"  "NSFTV154.dfdfb828.0"
## [286] "NSFTV24.e074229f.0"  "NSFTV246.e0815776.0" "NSFTV363.e32c9d62.0"
## [289] "NSFTV108.e3b049a9.0" "NSFTV152.e4533c8b.0" "NSFTV43.e5951598.0" 
## [292] "NSFTV367.e59cbfbe.0" "NSFTV208.e876c9ec.0" "NSFTV372.e8f708a5.0"
## [295] "NSFTV314.ea28671d.0" "NSFTV114.eac19fb8.0" "NSFTV236.eb17cda6.0"
## [298] "NSFTV636.ebb47c09.0" "NSFTV378.eca85a73.0" "NSFTV148.ed1beb9e.0"
## [301] "NSFTV55.ef8fd965.0"  "NSFTV70.efbb93de.0"  "NSFTV16.f0328a18.0" 
## [304] "NSFTV356.f22d2dc6.0" "NSFTV365.f2e723fc.0" "NSFTV224.f2e9bed5.0"
## [307] "NSFTV264.f2f12364.0" "NSFTV262.f30e146b.0" "NSFTV635.f3d9c979.0"
## [310] "NSFTV203.f46f03bf.0" "NSFTV644.f66e275e.0" "NSFTV223.f73fa9e9.0"
## [313] "NSFTV266.f7da2506.0" "NSFTV180.f948d7b5.0" "NSFTV336.fa180a1e.0"
## [316] "NSFTV21.fa4c4111.0"  "NSFTV327.fb2b3d82.0" "NSFTV310.fbae20bd.0"
## [319] "NSFTV276.fc26ce23.0" "NSFTV192.fdd98970.0" "NSFTV259.fe70ec33.0"
## [322] "NSFTV338.fed0b9a3.0" "NSFTV215.fed1dbae.0" "NSFTV395.ffde60f9.0"
## [325] "#N/A"

Locus 1 Os03g0843700 RAP2

list.files(pattern=".csv")
##  [1] "L1_Haplotype_Os03g0841800_RAP2_gene_and_pheno.csv"
##  [2] "L1_Haplotype_Os03g0843700_RAP2_gene.csv"          
##  [3] "L1_Haplotype_Os03g0845000_RAP2_gene.csv"          
##  [4] "L1_Haplotype_Os03g0845700_RAP2_gene.csv"          
##  [5] "L1_Haplotype_Os03g0845800_RAP2_gene.csv"          
##  [6] "L1_Haplotype_Os03g0848700_RAP2_gene.csv"          
##  [7] "L1_Haplotype_Os03g62400_gene_and_pheno.csv"       
##  [8] "L2_Haplotype_Os05g0420500_RAP2_gene.csv"          
##  [9] "L2_Haplotype_Os05g0420600_RAP2_gene.csv"          
## [10] "L2_Haplotype_Os05g0420900_RAP2_gene.csv"          
## [11] "L3_Haplotype_Os06g0346300_RAP2_gene.csv"          
## [12] "L4_Haplotype_Os07g0623200_RAP2_gene.csv"          
## [13] "L4_Haplotype_Os07g0623501_RAP2_gene.csv"          
## [14] "L4_Haplotype_Os07g0623600_RAP2_gene.csv"          
## [15] "L4_Haplotype_Os07g43040_gene.csv"                 
## [16] "L6_Haplotype_Os11g0216000_RAB2_gene.csv"          
## [17] "Locus_1_final_file.csv"                           
## [18] "Locus_2_final_file.csv"                           
## [19] "Locus_3_final_file.csv"                           
## [20] "Locus_4_final_file.csv"                           
## [21] "Locus_6_final_file.csv"
temp <- read.csv("L1_Haplotype_Os03g0843700_RAP2_gene.csv")
head(temp)
haplotype2 <- temp[,c(1, 14)]

ultimate2 <- merge(pheno2, haplotype2, by=c("VCF_better"))
head(ultimate2)
unique(ultimate2$haplotype)
## [1] "AAAAAAA" "AAAAATT" "TAAAATA" "AAAAATA" "ATATTTA"
# let's see which haplotype is most popular among accessions:

ultimate2 %>% group_by(haplotype) %>%
  tally()

In this case - we have 6 original haplotypes - but only 1 of them is represented by 3 or more accessions

ultimate2$haplotype <- gsub("AAAAATT", "hap01", ultimate2$haplotype)
ultimate2$haplotype <- gsub("AAAAAAA", "hap02", ultimate2$haplotype)
ultimate2$haplotype <- gsub("AAAAATA", "hap03", ultimate2$haplotype)
ultimate2$haplotype <- gsub("TAAAATA", "hap04", ultimate2$haplotype)
ultimate2$haplotype <- gsub("ATATTTA", "hap05", ultimate2$haplotype)

ultimate2$haplotype <- factor(ultimate2$haplotype, levels = c("hap01", "hap02", "hap03", "hap04", "hap05"))
colnames(ultimate2)
##  [1] "VCF_better"        "HDRA.ID"           "NSFTV.ID"         
##  [4] "IRGC.ID"           "Subpop"            "Accession.Name"   
##  [7] "VCF_code"          "index.pots"        "Variety"          
## [10] "Shoot.area..cm2."  "Hull.area..cm2."   "Solidity"         
## [13] "Perimeter..cm."    "Nr.leaves"         "Nr.tillers"       
## [16] "Plant.height..cm." "Culm.height..cm."  "Leaf.length..cm." 
## [19] "DW..g."            "Leaf.angle..º."    "Tiller.angle..º." 
## [22] "Erectness..º."     "Droopiness..º."    "SUM_rank_traits"  
## [25] "SUM_norm_traits"   "norm_SI_rank"      "SI_rank"          
## [28] "haplotype"

OK - now let’s graph them:

Output <- TukeyHSD(aov(aov(Shoot.area..cm2. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Shoot.area..cm2., colour = haplotype)) 
my_box_plot <- my_box_plot + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot <- my_box_plot + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot <- my_box_plot + theme(legend.position = "none")
my_box_plot <- my_box_plot + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3", "firebrick3", "firebrick3"))
my_box_plot <- my_box_plot + xlab("") + ylab(expression(paste("Shoot area (", cm^2, ")", sep = "")))
my_box_plot <- my_box_plot + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot <- my_box_plot + stat_pvalue_manual(test, label = "Tukey", y.position = 950)
my_box_plot

Output <- TukeyHSD(aov(aov(Hull.area..cm2. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot2 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Hull.area..cm2., colour = haplotype)) 
my_box_plot2 <- my_box_plot2 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot2 <- my_box_plot2 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot2 <- my_box_plot2 + theme(legend.position = "none")
my_box_plot2 <- my_box_plot2 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3", "firebrick3", "firebrick3"))
my_box_plot2 <- my_box_plot2 + xlab("") + ylab(expression(paste("Hull area (", cm^2, ")", sep = "")))
my_box_plot2 <- my_box_plot2 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot2 <- my_box_plot2 + stat_pvalue_manual(test, label = "Tukey", y.position = 6000)
my_box_plot2

Output <- TukeyHSD(aov(aov(Solidity ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot3 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Solidity, colour = haplotype)) 
my_box_plot3 <- my_box_plot3 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot3 <- my_box_plot3 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot3 <- my_box_plot3 + theme(legend.position = "none")
my_box_plot3 <- my_box_plot3 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3", "firebrick3", "firebrick3"))
my_box_plot3 <- my_box_plot3 + xlab("") + ylab("Solidity (a.u.)")
my_box_plot3 <- my_box_plot3 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot3 <- my_box_plot3 + stat_pvalue_manual(test, label = "Tukey", y.position = 0.21)
my_box_plot3

Output <- TukeyHSD(aov(aov(Perimeter..cm. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot4 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Perimeter..cm., colour = haplotype)) 
my_box_plot4 <- my_box_plot4 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot4 <- my_box_plot4 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot4 <- my_box_plot4 + theme(legend.position = "none")
my_box_plot4 <- my_box_plot4 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3", "firebrick3", "firebrick3"))
my_box_plot4 <- my_box_plot4 + xlab("") + ylab("Perimeter (cm)")
my_box_plot4 <- my_box_plot4 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot4 <- my_box_plot4 + stat_pvalue_manual(test, label = "Tukey", y.position = 2500)
my_box_plot4

Output <- TukeyHSD(aov(aov(Nr.leaves ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot5 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Nr.leaves, colour = haplotype)) 
my_box_plot5 <- my_box_plot5 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot5 <- my_box_plot5 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot5 <- my_box_plot5 + theme(legend.position = "none")
my_box_plot5 <- my_box_plot5 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3", "firebrick3", "firebrick3"))
my_box_plot5 <- my_box_plot5 + xlab("") + ylab("Leaf number")
my_box_plot5 <- my_box_plot5 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot5 <- my_box_plot5 + stat_pvalue_manual(test, label = "Tukey", y.position = 55)
my_box_plot5

Output <- TukeyHSD(aov(aov(Nr.tillers ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot6 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Nr.tillers, colour = haplotype)) 
my_box_plot6 <- my_box_plot6 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot6 <- my_box_plot6 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot6 <- my_box_plot6 + theme(legend.position = "none")
my_box_plot6 <- my_box_plot6 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3", "firebrick3", "firebrick3"))
my_box_plot6 <- my_box_plot6 + xlab("") + ylab("Tiller number")
my_box_plot6 <- my_box_plot6 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot6 <- my_box_plot6 + stat_pvalue_manual(test, label = "Tukey", y.position = 17)
my_box_plot6

Output <- TukeyHSD(aov(aov(Plant.height..cm. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot7 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Plant.height..cm., colour = haplotype)) 
my_box_plot7 <- my_box_plot7 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot7 <- my_box_plot7 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot7 <- my_box_plot7 + theme(legend.position = "none")
my_box_plot7 <- my_box_plot7 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3", "firebrick3", "firebrick3"))
my_box_plot7 <- my_box_plot7 + xlab("") + ylab("Plant height (cm)")
my_box_plot7 <- my_box_plot7 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot7 <- my_box_plot7 + stat_pvalue_manual(test, label = "Tukey", y.position = 95)
my_box_plot7

Output <- TukeyHSD(aov(aov(Culm.height..cm. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot8 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Culm.height..cm., colour = haplotype)) 
my_box_plot8 <- my_box_plot8 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot8 <- my_box_plot8 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot8 <- my_box_plot8 + theme(legend.position = "none")
my_box_plot8 <- my_box_plot8 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3", "firebrick3", "firebrick3"))
my_box_plot8 <- my_box_plot8 + xlab("") + ylab("Culm height (cm)")
my_box_plot8 <- my_box_plot8 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot8 <- my_box_plot8 + stat_pvalue_manual(test, label = "Tukey", y.position = 35)
my_box_plot8

Output <- TukeyHSD(aov(aov(Leaf.length..cm. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot9 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Leaf.length..cm., colour = haplotype)) 
my_box_plot9 <- my_box_plot9 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot9 <- my_box_plot9 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot9 <- my_box_plot9 + theme(legend.position = "none")
my_box_plot9 <- my_box_plot9 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3", "firebrick3", "firebrick3"))
my_box_plot9 <- my_box_plot9 + xlab("") + ylab("Leaf length (cm)")
my_box_plot9 <- my_box_plot9 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot9 <- my_box_plot9 + stat_pvalue_manual(test, label = "Tukey", y.position = 65)
my_box_plot9

Output <- TukeyHSD(aov(aov(DW..g. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot10 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = DW..g., colour = haplotype)) 
my_box_plot10 <- my_box_plot10 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot10 <- my_box_plot10 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot10 <- my_box_plot10 + theme(legend.position = "none")
my_box_plot10 <- my_box_plot10 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3", "firebrick3", "firebrick3"))
my_box_plot10 <- my_box_plot10 + xlab("") + ylab("Dry weight (g)")
my_box_plot10 <- my_box_plot10 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot10 <- my_box_plot10 + stat_pvalue_manual(test, label = "Tukey", y.position = 4)
my_box_plot10

Output <- TukeyHSD(aov(aov(Leaf.angle..º. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot11 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Leaf.angle..º., colour = haplotype)) 
my_box_plot11 <- my_box_plot11 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot11 <- my_box_plot11 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot11 <- my_box_plot11 + theme(legend.position = "none")
my_box_plot11 <- my_box_plot11 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3", "firebrick3", "firebrick3"))
my_box_plot11 <- my_box_plot11 + xlab("") + ylab("Leaf angle (º)")
my_box_plot11 <- my_box_plot11 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot11 <- my_box_plot11 + stat_pvalue_manual(test, label = "Tukey", y.position = 65)
my_box_plot11

Output <- TukeyHSD(aov(aov(Tiller.angle..º. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot12 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Tiller.angle..º., colour = haplotype)) 
my_box_plot12 <- my_box_plot12 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot12 <- my_box_plot12 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot12 <- my_box_plot12 + theme(legend.position = "none")
my_box_plot12 <- my_box_plot12 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3", "firebrick3", "firebrick3"))
my_box_plot12 <- my_box_plot12 + xlab("") + ylab("Tiller angle (º)")
my_box_plot12 <- my_box_plot12 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot12 <- my_box_plot12 + stat_pvalue_manual(test, label = "Tukey", y.position = 37)
my_box_plot12

Output <- TukeyHSD(aov(aov(Erectness..º. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot13 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Erectness..º., colour = haplotype)) 
my_box_plot13 <- my_box_plot13 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot13 <- my_box_plot13 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot13 <- my_box_plot13 + theme(legend.position = "none")
my_box_plot13 <- my_box_plot13 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3", "firebrick3", "firebrick3"))
my_box_plot13 <- my_box_plot13 + xlab("") + ylab("Erectness (º)")
my_box_plot13 <- my_box_plot13 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot13 <- my_box_plot13 + stat_pvalue_manual(test, label = "Tukey", y.position = 165)
my_box_plot13

Output <- TukeyHSD(aov(aov(Droopiness..º. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot14 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Droopiness..º., colour = haplotype)) 
my_box_plot14 <- my_box_plot14 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot14 <- my_box_plot14 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot14 <- my_box_plot14 + theme(legend.position = "none")
my_box_plot14 <- my_box_plot14 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3", "firebrick3", "firebrick3"))
my_box_plot14 <- my_box_plot14 + xlab("") + ylab("Droopiness (º)")
my_box_plot14 <- my_box_plot14 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot14 <- my_box_plot14 + stat_pvalue_manual(test, label = "Tukey", y.position = 170)
my_box_plot14

pdf("Locus_01_Haplotypes_Os03g0843700_RAP2_NEW.pdf", height = 30, width = 12)
plot_grid(my_box_plot,  my_box_plot2, my_box_plot3, my_box_plot4, my_box_plot5, my_box_plot6, my_box_plot7, my_box_plot8, my_box_plot9, my_box_plot10, my_box_plot11, my_box_plot12, my_box_plot13, my_box_plot14,
          ncol=2,
          align = "hv", labels=c("AUTO"), 
          label_size = 24)
dev.off()
## quartz_off_screen 
##                 2

Locus 1 Os03g0845000 RAP2

list.files(pattern=".csv")
##  [1] "L1_Haplotype_Os03g0841800_RAP2_gene_and_pheno.csv"
##  [2] "L1_Haplotype_Os03g0843700_RAP2_gene.csv"          
##  [3] "L1_Haplotype_Os03g0845000_RAP2_gene.csv"          
##  [4] "L1_Haplotype_Os03g0845700_RAP2_gene.csv"          
##  [5] "L1_Haplotype_Os03g0845800_RAP2_gene.csv"          
##  [6] "L1_Haplotype_Os03g0848700_RAP2_gene.csv"          
##  [7] "L1_Haplotype_Os03g62400_gene_and_pheno.csv"       
##  [8] "L2_Haplotype_Os05g0420500_RAP2_gene.csv"          
##  [9] "L2_Haplotype_Os05g0420600_RAP2_gene.csv"          
## [10] "L2_Haplotype_Os05g0420900_RAP2_gene.csv"          
## [11] "L3_Haplotype_Os06g0346300_RAP2_gene.csv"          
## [12] "L4_Haplotype_Os07g0623200_RAP2_gene.csv"          
## [13] "L4_Haplotype_Os07g0623501_RAP2_gene.csv"          
## [14] "L4_Haplotype_Os07g0623600_RAP2_gene.csv"          
## [15] "L4_Haplotype_Os07g43040_gene.csv"                 
## [16] "L6_Haplotype_Os11g0216000_RAB2_gene.csv"          
## [17] "Locus_1_final_file.csv"                           
## [18] "Locus_2_final_file.csv"                           
## [19] "Locus_3_final_file.csv"                           
## [20] "Locus_4_final_file.csv"                           
## [21] "Locus_6_final_file.csv"
temp <- read.csv("L1_Haplotype_Os03g0845000_RAP2_gene.csv")
head(temp)
haplotype2 <- temp[,c(1, 14)]

ultimate2 <- merge(pheno2, haplotype2, by=c("VCF_better"))
head(ultimate2)
unique(ultimate2$haplotype)
## [1] "AA" "TA" "TT"
# let's see which haplotype is most popular among accessions:

ultimate2 %>% group_by(haplotype) %>%
  tally()
ultimate2$haplotype <- gsub("TA", "hap01", ultimate2$haplotype)
ultimate2$haplotype <- gsub("AA", "hap02", ultimate2$haplotype)
ultimate2$haplotype <- gsub("TT", "hap03", ultimate2$haplotype)

ultimate2$haplotype <- factor(ultimate2$haplotype, levels = c("hap01", "hap02", "hap03"))
colnames(ultimate2)
##  [1] "VCF_better"        "HDRA.ID"           "NSFTV.ID"         
##  [4] "IRGC.ID"           "Subpop"            "Accession.Name"   
##  [7] "VCF_code"          "index.pots"        "Variety"          
## [10] "Shoot.area..cm2."  "Hull.area..cm2."   "Solidity"         
## [13] "Perimeter..cm."    "Nr.leaves"         "Nr.tillers"       
## [16] "Plant.height..cm." "Culm.height..cm."  "Leaf.length..cm." 
## [19] "DW..g."            "Leaf.angle..º."    "Tiller.angle..º." 
## [22] "Erectness..º."     "Droopiness..º."    "SUM_rank_traits"  
## [25] "SUM_norm_traits"   "norm_SI_rank"      "SI_rank"          
## [28] "haplotype"

OK - now let’s graph them:

Output <- TukeyHSD(aov(aov(Shoot.area..cm2. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Shoot.area..cm2., colour = haplotype)) 
my_box_plot <- my_box_plot + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot <- my_box_plot + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot <- my_box_plot + theme(legend.position = "none")
my_box_plot <- my_box_plot + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot <- my_box_plot + xlab("") + ylab(expression(paste("Shoot area (", cm^2, ")", sep = "")))
my_box_plot <- my_box_plot + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot <- my_box_plot + stat_pvalue_manual(test, label = "Tukey", y.position = 950)
my_box_plot

Output <- TukeyHSD(aov(aov(Hull.area..cm2. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot2 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Hull.area..cm2., colour = haplotype)) 
my_box_plot2 <- my_box_plot2 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot2 <- my_box_plot2 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot2 <- my_box_plot2 + theme(legend.position = "none")
my_box_plot2 <- my_box_plot2 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot2 <- my_box_plot2 + xlab("") + ylab(expression(paste("Hull area (", cm^2, ")", sep = "")))
my_box_plot2 <- my_box_plot2 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot2 <- my_box_plot2 + stat_pvalue_manual(test, label = "Tukey", y.position = 6000)
my_box_plot2

Output <- TukeyHSD(aov(aov(Solidity ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot3 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Solidity, colour = haplotype)) 
my_box_plot3 <- my_box_plot3 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot3 <- my_box_plot3 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot3 <- my_box_plot3 + theme(legend.position = "none")
my_box_plot3 <- my_box_plot3 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot3 <- my_box_plot3 + xlab("") + ylab("Solidity (a.u.)")
my_box_plot3 <- my_box_plot3 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot3 <- my_box_plot3 + stat_pvalue_manual(test, label = "Tukey", y.position = 0.21)
my_box_plot3

Output <- TukeyHSD(aov(aov(Perimeter..cm. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot4 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Perimeter..cm., colour = haplotype)) 
my_box_plot4 <- my_box_plot4 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot4 <- my_box_plot4 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot4 <- my_box_plot4 + theme(legend.position = "none")
my_box_plot4 <- my_box_plot4 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot4 <- my_box_plot4 + xlab("") + ylab("Perimeter (cm)")
my_box_plot4 <- my_box_plot4 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot4 <- my_box_plot4 + stat_pvalue_manual(test, label = "Tukey", y.position = 2500)
my_box_plot4

Output <- TukeyHSD(aov(aov(Nr.leaves ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot5 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Nr.leaves, colour = haplotype)) 
my_box_plot5 <- my_box_plot5 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot5 <- my_box_plot5 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot5 <- my_box_plot5 + theme(legend.position = "none")
my_box_plot5 <- my_box_plot5 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot5 <- my_box_plot5 + xlab("") + ylab("Leaf number")
my_box_plot5 <- my_box_plot5 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot5 <- my_box_plot5 + stat_pvalue_manual(test, label = "Tukey", y.position = 55)
my_box_plot5

Output <- TukeyHSD(aov(aov(Nr.tillers ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot6 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Nr.tillers, colour = haplotype)) 
my_box_plot6 <- my_box_plot6 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot6 <- my_box_plot6 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot6 <- my_box_plot6 + theme(legend.position = "none")
my_box_plot6 <- my_box_plot6 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot6 <- my_box_plot6 + xlab("") + ylab("Tiller number")
my_box_plot6 <- my_box_plot6 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot6 <- my_box_plot6 + stat_pvalue_manual(test, label = "Tukey", y.position = 17)
my_box_plot6

Output <- TukeyHSD(aov(aov(Plant.height..cm. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot7 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Plant.height..cm., colour = haplotype)) 
my_box_plot7 <- my_box_plot7 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot7 <- my_box_plot7 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot7 <- my_box_plot7 + theme(legend.position = "none")
my_box_plot7 <- my_box_plot7 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot7 <- my_box_plot7 + xlab("") + ylab("Plant height (cm)")
my_box_plot7 <- my_box_plot7 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot7 <- my_box_plot7 + stat_pvalue_manual(test, label = "Tukey", y.position = 95)
my_box_plot7

Output <- TukeyHSD(aov(aov(Culm.height..cm. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot8 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Culm.height..cm., colour = haplotype)) 
my_box_plot8 <- my_box_plot8 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot8 <- my_box_plot8 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot8 <- my_box_plot8 + theme(legend.position = "none")
my_box_plot8 <- my_box_plot8 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot8 <- my_box_plot8 + xlab("") + ylab("Culm height (cm)")
my_box_plot8 <- my_box_plot8 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot8 <- my_box_plot8 + stat_pvalue_manual(test, label = "Tukey", y.position = 35)
my_box_plot8

Output <- TukeyHSD(aov(aov(Leaf.length..cm. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot9 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Leaf.length..cm., colour = haplotype)) 
my_box_plot9 <- my_box_plot9 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot9 <- my_box_plot9 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot9 <- my_box_plot9 + theme(legend.position = "none")
my_box_plot9 <- my_box_plot9 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot9 <- my_box_plot9 + xlab("") + ylab("Leaf length (cm)")
my_box_plot9 <- my_box_plot9 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot9 <- my_box_plot9 + stat_pvalue_manual(test, label = "Tukey", y.position = 65)
my_box_plot9

Output <- TukeyHSD(aov(aov(DW..g. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot10 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = DW..g., colour = haplotype)) 
my_box_plot10 <- my_box_plot10 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot10 <- my_box_plot10 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot10 <- my_box_plot10 + theme(legend.position = "none")
my_box_plot10 <- my_box_plot10 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot10 <- my_box_plot10 + xlab("") + ylab("Dry weight (g)")
my_box_plot10 <- my_box_plot10 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot10 <- my_box_plot10 + stat_pvalue_manual(test, label = "Tukey", y.position = 4) 
my_box_plot10

Output <- TukeyHSD(aov(aov(Leaf.angle..º. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot11 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Leaf.angle..º., colour = haplotype)) 
my_box_plot11 <- my_box_plot11 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot11 <- my_box_plot11 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot11 <- my_box_plot11 + theme(legend.position = "none")
my_box_plot11 <- my_box_plot11 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot11 <- my_box_plot11 + xlab("") + ylab("Leaf angle (º)")
my_box_plot11 <- my_box_plot11 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot11 <- my_box_plot11 + stat_pvalue_manual(test, label = "Tukey", y.position = 65)
my_box_plot11

Output <- TukeyHSD(aov(aov(Tiller.angle..º. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot12 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Tiller.angle..º., colour = haplotype)) 
my_box_plot12 <- my_box_plot12 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot12 <- my_box_plot12 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot12 <- my_box_plot12 + theme(legend.position = "none")
my_box_plot12 <- my_box_plot12 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot12 <- my_box_plot12 + xlab("") + ylab("Tiller angle (º)")
my_box_plot12 <- my_box_plot12 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot12 <- my_box_plot12 + stat_pvalue_manual(test, label = "Tukey", y.position = 37)
my_box_plot12

Output <- TukeyHSD(aov(aov(Erectness..º. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot13 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Erectness..º., colour = haplotype)) 
my_box_plot13 <- my_box_plot13 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot13 <- my_box_plot13 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot13 <- my_box_plot13 + theme(legend.position = "none")
my_box_plot13 <- my_box_plot13 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot13 <- my_box_plot13 + xlab("") + ylab("Erectness (º)")
my_box_plot13 <- my_box_plot13 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot13 <- my_box_plot13 + stat_pvalue_manual(test, label = "Tukey", y.position = 165)
my_box_plot13

Output <- TukeyHSD(aov(aov(Droopiness..º. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot14 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Droopiness..º., colour = haplotype)) 
my_box_plot14 <- my_box_plot14 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot14 <- my_box_plot14 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot14 <- my_box_plot14 + theme(legend.position = "none")
my_box_plot14 <- my_box_plot14 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot14 <- my_box_plot14 + xlab("") + ylab("Droopiness (º)")
my_box_plot14 <- my_box_plot14 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot14 <- my_box_plot14 + stat_pvalue_manual(test, label = "Tukey", y.position = 170)
my_box_plot14

pdf("Locus_01_Haplotypes_Os03g0845000_RAP2_NEW.pdf", height = 30, width = 12)
plot_grid(my_box_plot,  my_box_plot2, my_box_plot3, my_box_plot4, my_box_plot5, my_box_plot6, my_box_plot7, my_box_plot8, my_box_plot9, my_box_plot10, my_box_plot11, my_box_plot12, my_box_plot13, my_box_plot14,
          ncol=2,
          align = "hv", labels=c("AUTO"), 
          label_size = 24)
dev.off()
## quartz_off_screen 
##                 2

Locus 1 Os03g0845700 RAP2

list.files(pattern=".csv")
##  [1] "L1_Haplotype_Os03g0841800_RAP2_gene_and_pheno.csv"
##  [2] "L1_Haplotype_Os03g0843700_RAP2_gene.csv"          
##  [3] "L1_Haplotype_Os03g0845000_RAP2_gene.csv"          
##  [4] "L1_Haplotype_Os03g0845700_RAP2_gene.csv"          
##  [5] "L1_Haplotype_Os03g0845800_RAP2_gene.csv"          
##  [6] "L1_Haplotype_Os03g0848700_RAP2_gene.csv"          
##  [7] "L1_Haplotype_Os03g62400_gene_and_pheno.csv"       
##  [8] "L2_Haplotype_Os05g0420500_RAP2_gene.csv"          
##  [9] "L2_Haplotype_Os05g0420600_RAP2_gene.csv"          
## [10] "L2_Haplotype_Os05g0420900_RAP2_gene.csv"          
## [11] "L3_Haplotype_Os06g0346300_RAP2_gene.csv"          
## [12] "L4_Haplotype_Os07g0623200_RAP2_gene.csv"          
## [13] "L4_Haplotype_Os07g0623501_RAP2_gene.csv"          
## [14] "L4_Haplotype_Os07g0623600_RAP2_gene.csv"          
## [15] "L4_Haplotype_Os07g43040_gene.csv"                 
## [16] "L6_Haplotype_Os11g0216000_RAB2_gene.csv"          
## [17] "Locus_1_final_file.csv"                           
## [18] "Locus_2_final_file.csv"                           
## [19] "Locus_3_final_file.csv"                           
## [20] "Locus_4_final_file.csv"                           
## [21] "Locus_6_final_file.csv"
temp <- read.csv("L1_Haplotype_Os03g0845700_RAP2_gene.csv")
head(temp)
haplotype2 <- temp[,c(1, 14)]

ultimate2 <- merge(pheno2, haplotype2, by=c("VCF_better"))
head(ultimate2)
unique(ultimate2$haplotype)
## [1] "A" "T"
# let's see which haplotype is most popular among accessions:

ultimate2 %>% group_by(haplotype) %>%
  tally()
ultimate2$haplotype <- gsub("T", "hap01", ultimate2$haplotype)
ultimate2$haplotype <- gsub("A", "hap02", ultimate2$haplotype)

ultimate2$haplotype <- factor(ultimate2$haplotype, levels = c("hap01", "hap02"))
colnames(ultimate2)
##  [1] "VCF_better"        "HDRA.ID"           "NSFTV.ID"         
##  [4] "IRGC.ID"           "Subpop"            "Accession.Name"   
##  [7] "VCF_code"          "index.pots"        "Variety"          
## [10] "Shoot.area..cm2."  "Hull.area..cm2."   "Solidity"         
## [13] "Perimeter..cm."    "Nr.leaves"         "Nr.tillers"       
## [16] "Plant.height..cm." "Culm.height..cm."  "Leaf.length..cm." 
## [19] "DW..g."            "Leaf.angle..º."    "Tiller.angle..º." 
## [22] "Erectness..º."     "Droopiness..º."    "SUM_rank_traits"  
## [25] "SUM_norm_traits"   "norm_SI_rank"      "SI_rank"          
## [28] "haplotype"

OK - now let’s graph them:

Output <- TukeyHSD(aov(aov(Shoot.area..cm2. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
P4
## [1] 7.428502e-13
test2 <- test[c(1,3),]

my_box_plot <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Shoot.area..cm2., colour = haplotype)) 
my_box_plot <- my_box_plot + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot <- my_box_plot + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot <- my_box_plot + theme(legend.position = "none")
my_box_plot <- my_box_plot + scale_colour_manual(values = c("steelblue", "firebrick3"))
my_box_plot <- my_box_plot + xlab("") + ylab(expression(paste("Shoot area (", cm^2, ")", sep = "")))
my_box_plot <- my_box_plot + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot <- my_box_plot + stat_pvalue_manual(test2, label = "Tukey", y.position = 950)
my_box_plot

Output <- TukeyHSD(aov(aov(Hull.area..cm2. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
P4
## [1] 6.97653e-12
test2 <- test[c(1,3),]

my_box_plot2 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Hull.area..cm2., colour = haplotype)) 
my_box_plot2 <- my_box_plot2 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot2 <- my_box_plot2 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot2 <- my_box_plot2 + theme(legend.position = "none")
my_box_plot2 <- my_box_plot2 + scale_colour_manual(values = c("steelblue", "firebrick3"))
my_box_plot2 <- my_box_plot2 + xlab("") + ylab(expression(paste("Hull area (", cm^2, ")", sep = "")))
my_box_plot2 <- my_box_plot2 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot2 <- my_box_plot2 + stat_pvalue_manual(test2, label = "Tukey", y.position = 6000)
my_box_plot2

Output <- TukeyHSD(aov(aov(Solidity ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
P4
## [1] 8.019962e-11
test2 <- test[c(1,3),]

my_box_plot3 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Solidity, colour = haplotype)) 
my_box_plot3 <- my_box_plot3 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot3 <- my_box_plot3 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot3 <- my_box_plot3 + theme(legend.position = "none")
my_box_plot3 <- my_box_plot3 + scale_colour_manual(values = c("steelblue", "firebrick3"))
my_box_plot3 <- my_box_plot3 + xlab("") + ylab("Solidity (a.u.)")
my_box_plot3 <- my_box_plot3 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot3 <- my_box_plot3 + stat_pvalue_manual(test2, label = "Tukey", y.position = 0.21)
my_box_plot3

Output <- TukeyHSD(aov(aov(Perimeter..cm. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
P4
## [1] 7.610579e-13
my_box_plot4 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Perimeter..cm., colour = haplotype)) 
my_box_plot4 <- my_box_plot4 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot4 <- my_box_plot4 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot4 <- my_box_plot4 + theme(legend.position = "none")
my_box_plot4 <- my_box_plot4 + scale_colour_manual(values = c("steelblue", "firebrick3"))
my_box_plot4 <- my_box_plot4 + xlab("") + ylab("Perimeter (cm)")
my_box_plot4 <- my_box_plot4 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot4 <- my_box_plot4 + stat_pvalue_manual(test2, label = "Tukey", y.position = 2500)
my_box_plot4

Output <- TukeyHSD(aov(aov(Nr.leaves ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
P4
## [1] 7.456258e-13
my_box_plot5 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Nr.leaves, colour = haplotype)) 
my_box_plot5 <- my_box_plot5 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot5 <- my_box_plot5 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot5 <- my_box_plot5 + theme(legend.position = "none")
my_box_plot5 <- my_box_plot5 + scale_colour_manual(values = c("steelblue", "firebrick3"))
my_box_plot5 <- my_box_plot5 + xlab("") + ylab("Leaf number")
my_box_plot5 <- my_box_plot5 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot5 <- my_box_plot5 + stat_pvalue_manual(test2, label = "Tukey", y.position = 55)
my_box_plot5

Output <- TukeyHSD(aov(aov(Nr.tillers ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
P4
## [1] 7.444045e-13
my_box_plot6 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Nr.tillers, colour = haplotype)) 
my_box_plot6 <- my_box_plot6 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot6 <- my_box_plot6 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot6 <- my_box_plot6 + theme(legend.position = "none")
my_box_plot6 <- my_box_plot6 + scale_colour_manual(values = c("steelblue", "firebrick3"))
my_box_plot6 <- my_box_plot6 + xlab("") + ylab("Tiller number")
my_box_plot6 <- my_box_plot6 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot6 <- my_box_plot6 + stat_pvalue_manual(test2, label = "Tukey", y.position = 17)
my_box_plot6

Output <- TukeyHSD(aov(aov(Plant.height..cm. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
P4
## [1] 0.1882969
test3 <- test2
test3$Tukey <- gsub("b", "a", test3$Tukey)

my_box_plot7 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Plant.height..cm., colour = haplotype)) 
my_box_plot7 <- my_box_plot7 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot7 <- my_box_plot7 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot7 <- my_box_plot7 + theme(legend.position = "none")
my_box_plot7 <- my_box_plot7 + scale_colour_manual(values = c("steelblue", "firebrick3"))
my_box_plot7 <- my_box_plot7 + xlab("") + ylab("Plant height (cm)")
my_box_plot7 <- my_box_plot7 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot7 <- my_box_plot7 + stat_pvalue_manual(test3, label = "Tukey", y.position = 95)
my_box_plot7

Output <- TukeyHSD(aov(aov(Culm.height..cm. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
P4
## [1] 0.01431405
my_box_plot8 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Culm.height..cm., colour = haplotype)) 
my_box_plot8 <- my_box_plot8 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot8 <- my_box_plot8 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot8 <- my_box_plot8 + theme(legend.position = "none")
my_box_plot8 <- my_box_plot8 + scale_colour_manual(values = c("steelblue", "firebrick3"))
my_box_plot8 <- my_box_plot8 + xlab("") + ylab("Culm height (cm)")
my_box_plot8 <- my_box_plot8 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot8 <- my_box_plot8 + stat_pvalue_manual(test2, label = "Tukey", y.position = 35)
my_box_plot8

Output <- TukeyHSD(aov(aov(Leaf.length..cm. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
P4
## [1] 0.009384355
my_box_plot9 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Leaf.length..cm., colour = haplotype)) 
my_box_plot9 <- my_box_plot9 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot9 <- my_box_plot9 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot9 <- my_box_plot9 + theme(legend.position = "none")
my_box_plot9 <- my_box_plot9 + scale_colour_manual(values = c("steelblue", "firebrick3"))
my_box_plot9 <- my_box_plot9 + xlab("") + ylab("Leaf length (cm)")
my_box_plot9 <- my_box_plot9 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot9 <- my_box_plot9 + stat_pvalue_manual(test2, label = "Tukey", y.position = 65)
my_box_plot9

Output <- TukeyHSD(aov(aov(DW..g. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
P4
## [1] 7.456258e-13
my_box_plot10 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = DW..g., colour = haplotype)) 
my_box_plot10 <- my_box_plot10 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot10 <- my_box_plot10 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot10 <- my_box_plot10 + theme(legend.position = "none")
my_box_plot10 <- my_box_plot10 + scale_colour_manual(values = c("steelblue", "firebrick3"))
my_box_plot10 <- my_box_plot10 + xlab("") + ylab("Dry weight (g)")
my_box_plot10 <- my_box_plot10 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot10 <- my_box_plot10 + stat_pvalue_manual(test2, label = "Tukey", y.position = 4) 
my_box_plot10

Output <- TukeyHSD(aov(aov(Leaf.angle..º. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
P4
## [1] 0.02075252
my_box_plot11 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Leaf.angle..º., colour = haplotype)) 
my_box_plot11 <- my_box_plot11 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot11 <- my_box_plot11 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot11 <- my_box_plot11 + theme(legend.position = "none")
my_box_plot11 <- my_box_plot11 + scale_colour_manual(values = c("steelblue", "firebrick3"))
my_box_plot11 <- my_box_plot11 + xlab("") + ylab("Leaf angle (º)")
my_box_plot11 <- my_box_plot11 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot11 <- my_box_plot11 + stat_pvalue_manual(test2, label = "Tukey", y.position = 65)
my_box_plot11

Output <- TukeyHSD(aov(aov(Tiller.angle..º. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
P4
## [1] 1.198577e-06
my_box_plot12 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Tiller.angle..º., colour = haplotype)) 
my_box_plot12 <- my_box_plot12 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot12 <- my_box_plot12 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot12 <- my_box_plot12 + theme(legend.position = "none")
my_box_plot12 <- my_box_plot12 + scale_colour_manual(values = c("steelblue", "firebrick3"))
my_box_plot12 <- my_box_plot12 + xlab("") + ylab("Tiller angle (º)")
my_box_plot12 <- my_box_plot12 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot12 <- my_box_plot12 + stat_pvalue_manual(test2, label = "Tukey", y.position = 37)
my_box_plot12

Output <- TukeyHSD(aov(aov(Erectness..º. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
P4
## [1] 4.694578e-12
my_box_plot13 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Erectness..º., colour = haplotype)) 
my_box_plot13 <- my_box_plot13 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot13 <- my_box_plot13 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot13 <- my_box_plot13 + theme(legend.position = "none")
my_box_plot13 <- my_box_plot13 + scale_colour_manual(values = c("steelblue", "firebrick3"))
my_box_plot13 <- my_box_plot13 + xlab("") + ylab("Erectness (º)")
my_box_plot13 <- my_box_plot13 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot13 <- my_box_plot13 + stat_pvalue_manual(test2, label = "Tukey", y.position = 165)
my_box_plot13

Output <- TukeyHSD(aov(aov(Droopiness..º. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
P4
## [1] 4.694578e-12
my_box_plot14 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Droopiness..º., colour = haplotype)) 
my_box_plot14 <- my_box_plot14 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot14 <- my_box_plot14 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot14 <- my_box_plot14 + theme(legend.position = "none")
my_box_plot14 <- my_box_plot14 + scale_colour_manual(values = c("steelblue", "firebrick3"))
my_box_plot14 <- my_box_plot14 + xlab("") + ylab("Droopiness (º)")
my_box_plot14 <- my_box_plot14 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot14 <- my_box_plot14 + stat_pvalue_manual(test2, label = "Tukey", y.position = 170)
my_box_plot14

pdf("Locus_01_Haplotypes_Os03g0845700_RAP2_NEW.pdf", height = 30, width = 12)
plot_grid(my_box_plot,  my_box_plot2, my_box_plot3, my_box_plot4, my_box_plot5, my_box_plot6, 
          my_box_plot7, my_box_plot8, my_box_plot9, my_box_plot10, my_box_plot11, my_box_plot12, 
          my_box_plot13, my_box_plot14,
          ncol=2,
          align = "hv", labels=c("AUTO"), 
          label_size = 24)
dev.off()
## quartz_off_screen 
##                 2

Locus 2 Os5g0420600

list.files(pattern=".csv")
##  [1] "L1_Haplotype_Os03g0841800_RAP2_gene_and_pheno.csv"
##  [2] "L1_Haplotype_Os03g0843700_RAP2_gene.csv"          
##  [3] "L1_Haplotype_Os03g0845000_RAP2_gene.csv"          
##  [4] "L1_Haplotype_Os03g0845700_RAP2_gene.csv"          
##  [5] "L1_Haplotype_Os03g0845800_RAP2_gene.csv"          
##  [6] "L1_Haplotype_Os03g0848700_RAP2_gene.csv"          
##  [7] "L1_Haplotype_Os03g62400_gene_and_pheno.csv"       
##  [8] "L2_Haplotype_Os05g0420500_RAP2_gene.csv"          
##  [9] "L2_Haplotype_Os05g0420600_RAP2_gene.csv"          
## [10] "L2_Haplotype_Os05g0420900_RAP2_gene.csv"          
## [11] "L3_Haplotype_Os06g0346300_RAP2_gene.csv"          
## [12] "L4_Haplotype_Os07g0623200_RAP2_gene.csv"          
## [13] "L4_Haplotype_Os07g0623501_RAP2_gene.csv"          
## [14] "L4_Haplotype_Os07g0623600_RAP2_gene.csv"          
## [15] "L4_Haplotype_Os07g43040_gene.csv"                 
## [16] "L6_Haplotype_Os11g0216000_RAB2_gene.csv"          
## [17] "Locus_1_final_file.csv"                           
## [18] "Locus_2_final_file.csv"                           
## [19] "Locus_3_final_file.csv"                           
## [20] "Locus_4_final_file.csv"                           
## [21] "Locus_6_final_file.csv"
temp <- read.csv("L2_Haplotype_Os05g0420600_RAP2_gene.csv")
head(temp)
haplotype2 <- temp[,c(1, 14)]

ultimate2 <- merge(pheno2, haplotype2, by=c("VCF_better"))
head(ultimate2)
unique(ultimate2$haplotype)
## [1] "AAAA" "ATAA" "TTAA"
# let's see which haplotype is most popular among accessions:

ultimate2 %>% group_by(haplotype) %>%
  tally()

In this case - we have 6 original haplotypes - but only 1 of them is represented by 3 or more accessions

ultimate2$haplotype <- gsub("AAAA", "hap01", ultimate2$haplotype)
ultimate2$haplotype <- gsub("ATAA", "hap02", ultimate2$haplotype)
ultimate2$haplotype <- gsub("TTAA", "hap03", ultimate2$haplotype)

ultimate2$haplotype <- factor(ultimate2$haplotype, levels = c("hap01", "hap02", "hap03"))
colnames(ultimate2)
##  [1] "VCF_better"        "HDRA.ID"           "NSFTV.ID"         
##  [4] "IRGC.ID"           "Subpop"            "Accession.Name"   
##  [7] "VCF_code"          "index.pots"        "Variety"          
## [10] "Shoot.area..cm2."  "Hull.area..cm2."   "Solidity"         
## [13] "Perimeter..cm."    "Nr.leaves"         "Nr.tillers"       
## [16] "Plant.height..cm." "Culm.height..cm."  "Leaf.length..cm." 
## [19] "DW..g."            "Leaf.angle..º."    "Tiller.angle..º." 
## [22] "Erectness..º."     "Droopiness..º."    "SUM_rank_traits"  
## [25] "SUM_norm_traits"   "norm_SI_rank"      "SI_rank"          
## [28] "haplotype"
Output <- TukeyHSD(aov(aov(Shoot.area..cm2. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Shoot.area..cm2., colour = haplotype)) 
my_box_plot <- my_box_plot + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot <- my_box_plot + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot <- my_box_plot + theme(legend.position = "none")
my_box_plot <- my_box_plot + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot <- my_box_plot + xlab("") + ylab(expression(paste("Shoot area (", cm^2, ")", sep = "")))
my_box_plot <- my_box_plot + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot <- my_box_plot + stat_pvalue_manual(test, label = "Tukey", y.position = 950)
my_box_plot

Output <- TukeyHSD(aov(aov(Hull.area..cm2. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot2 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Hull.area..cm2., colour = haplotype)) 
my_box_plot2 <- my_box_plot2 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot2 <- my_box_plot2 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot2 <- my_box_plot2 + theme(legend.position = "none")
my_box_plot2 <- my_box_plot2 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot2 <- my_box_plot2 + xlab("") + ylab(expression(paste("Hull area (", cm^2, ")", sep = "")))
my_box_plot2 <- my_box_plot2 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot2 <- my_box_plot2 + stat_pvalue_manual(test, label = "Tukey", y.position = 6000)
my_box_plot2

Output <- TukeyHSD(aov(aov(Solidity ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot3 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Solidity, colour = haplotype)) 
my_box_plot3 <- my_box_plot3 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot3 <- my_box_plot3 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot3 <- my_box_plot3 + theme(legend.position = "none")
my_box_plot3 <- my_box_plot3 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot3 <- my_box_plot3 + xlab("") + ylab("Solidity (a.u.)")
my_box_plot3 <- my_box_plot3 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot3 <- my_box_plot3 + stat_pvalue_manual(test, label = "Tukey", y.position = 0.21)
my_box_plot3

Output <- TukeyHSD(aov(aov(Perimeter..cm. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot4 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Perimeter..cm., colour = haplotype)) 
my_box_plot4 <- my_box_plot4 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot4 <- my_box_plot4 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot4 <- my_box_plot4 + theme(legend.position = "none")
my_box_plot4 <- my_box_plot4 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot4 <- my_box_plot4 + xlab("") + ylab("Perimeter (cm)")
my_box_plot4 <- my_box_plot4 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot4 <- my_box_plot4 + stat_pvalue_manual(test, label = "Tukey", y.position = 2500)
my_box_plot4

Output <- TukeyHSD(aov(aov(Nr.leaves ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot5 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Nr.leaves, colour = haplotype)) 
my_box_plot5 <- my_box_plot5 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot5 <- my_box_plot5 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot5 <- my_box_plot5 + theme(legend.position = "none")
my_box_plot5 <- my_box_plot5 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot5 <- my_box_plot5 + xlab("") + ylab("Leaf number")
my_box_plot5 <- my_box_plot5 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot5 <- my_box_plot5 + stat_pvalue_manual(test, label = "Tukey", y.position = 55)
my_box_plot5

Output <- TukeyHSD(aov(aov(Nr.tillers ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot6 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Nr.tillers, colour = haplotype)) 
my_box_plot6 <- my_box_plot6 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot6 <- my_box_plot6 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot6 <- my_box_plot6 + theme(legend.position = "none")
my_box_plot6 <- my_box_plot6 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot6 <- my_box_plot6 + xlab("") + ylab("Tiller number")
my_box_plot6 <- my_box_plot6 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot6 <- my_box_plot6 + stat_pvalue_manual(test, label = "Tukey", y.position = 17)
my_box_plot6

Output <- TukeyHSD(aov(aov(Plant.height..cm. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot7 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Plant.height..cm., colour = haplotype)) 
my_box_plot7 <- my_box_plot7 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot7 <- my_box_plot7 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot7 <- my_box_plot7 + theme(legend.position = "none")
my_box_plot7 <- my_box_plot7 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot7 <- my_box_plot7 + xlab("") + ylab("Plant height (cm)")
my_box_plot7 <- my_box_plot7 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot7 <- my_box_plot7 + stat_pvalue_manual(test, label = "Tukey", y.position = 95)
my_box_plot7

Output <- TukeyHSD(aov(aov(Culm.height..cm. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot8 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Culm.height..cm., colour = haplotype)) 
my_box_plot8 <- my_box_plot8 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot8 <- my_box_plot8 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot8 <- my_box_plot8 + theme(legend.position = "none")
my_box_plot8 <- my_box_plot8 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot8 <- my_box_plot8 + xlab("") + ylab("Culm height (cm)")
my_box_plot8 <- my_box_plot8 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot8 <- my_box_plot8 + stat_pvalue_manual(test, label = "Tukey", y.position = 35)
my_box_plot8

Output <- TukeyHSD(aov(aov(Leaf.length..cm. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot9 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Leaf.length..cm., colour = haplotype)) 
my_box_plot9 <- my_box_plot9 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot9 <- my_box_plot9 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot9 <- my_box_plot9 + theme(legend.position = "none")
my_box_plot9 <- my_box_plot9 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot9 <- my_box_plot9 + xlab("") + ylab("Leaf length (cm)")
my_box_plot9 <- my_box_plot9 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot9 <- my_box_plot9 + stat_pvalue_manual(test, label = "Tukey", y.position = 65)
my_box_plot9

Output <- TukeyHSD(aov(aov(DW..g. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot10 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = DW..g., colour = haplotype)) 
my_box_plot10 <- my_box_plot10 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot10 <- my_box_plot10 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot10 <- my_box_plot10 + theme(legend.position = "none")
my_box_plot10 <- my_box_plot10 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot10 <- my_box_plot10 + xlab("") + ylab("Dry weight (g)")
my_box_plot10 <- my_box_plot10 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot10 <- my_box_plot10 + stat_pvalue_manual(test, label = "Tukey", y.position = 4) 
my_box_plot10

Output <- TukeyHSD(aov(aov(Leaf.angle..º. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot11 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Leaf.angle..º., colour = haplotype)) 
my_box_plot11 <- my_box_plot11 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot11 <- my_box_plot11 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot11 <- my_box_plot11 + theme(legend.position = "none")
my_box_plot11 <- my_box_plot11 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot11 <- my_box_plot11 + xlab("") + ylab("Leaf angle (º)")
my_box_plot11 <- my_box_plot11 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot11 <- my_box_plot11 + stat_pvalue_manual(test, label = "Tukey", y.position = 65)
my_box_plot11

Output <- TukeyHSD(aov(aov(Tiller.angle..º. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot12 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Tiller.angle..º., colour = haplotype)) 
my_box_plot12 <- my_box_plot12 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot12 <- my_box_plot12 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot12 <- my_box_plot12 + theme(legend.position = "none")
my_box_plot12 <- my_box_plot12 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot12 <- my_box_plot12 + xlab("") + ylab("Tiller angle (º)")
my_box_plot12 <- my_box_plot12 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot12 <- my_box_plot12 + stat_pvalue_manual(test, label = "Tukey", y.position = 37)
my_box_plot12

Output <- TukeyHSD(aov(aov(Erectness..º. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot13 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Erectness..º., colour = haplotype)) 
my_box_plot13 <- my_box_plot13 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot13 <- my_box_plot13 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot13 <- my_box_plot13 + theme(legend.position = "none")
my_box_plot13 <- my_box_plot13 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot13 <- my_box_plot13 + xlab("") + ylab("Erectness (º)")
my_box_plot13 <- my_box_plot13 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot13 <- my_box_plot13 + stat_pvalue_manual(test, label = "Tukey", y.position = 165)
my_box_plot13

Output <- TukeyHSD(aov(aov(Droopiness..º. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot14 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Droopiness..º., colour = haplotype)) 
my_box_plot14 <- my_box_plot14 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot14 <- my_box_plot14 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot14 <- my_box_plot14 + theme(legend.position = "none")
my_box_plot14 <- my_box_plot14 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot14 <- my_box_plot14 + xlab("") + ylab("Droopiness (º)")
my_box_plot14 <- my_box_plot14 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot14 <- my_box_plot14 + stat_pvalue_manual(test, label = "Tukey", y.position = 170)
my_box_plot14

pdf("Locus_02_Haplotypes_Os05g0420600_RAP2_NEW.pdf", height = 30, width = 12)
plot_grid(my_box_plot,  my_box_plot2, my_box_plot3, my_box_plot4, my_box_plot5, my_box_plot6, my_box_plot7, my_box_plot8, my_box_plot9, my_box_plot10, my_box_plot11, my_box_plot12, my_box_plot13, my_box_plot14,
          ncol=2,
          align = "hv", labels=c("AUTO"), 
          label_size = 24)
dev.off()
## quartz_off_screen 
##                 2

Locus 2 Os05g0420900 RAP2

list.files(pattern=".csv")
##  [1] "L1_Haplotype_Os03g0841800_RAP2_gene_and_pheno.csv"
##  [2] "L1_Haplotype_Os03g0843700_RAP2_gene.csv"          
##  [3] "L1_Haplotype_Os03g0845000_RAP2_gene.csv"          
##  [4] "L1_Haplotype_Os03g0845700_RAP2_gene.csv"          
##  [5] "L1_Haplotype_Os03g0845800_RAP2_gene.csv"          
##  [6] "L1_Haplotype_Os03g0848700_RAP2_gene.csv"          
##  [7] "L1_Haplotype_Os03g62400_gene_and_pheno.csv"       
##  [8] "L2_Haplotype_Os05g0420500_RAP2_gene.csv"          
##  [9] "L2_Haplotype_Os05g0420600_RAP2_gene.csv"          
## [10] "L2_Haplotype_Os05g0420900_RAP2_gene.csv"          
## [11] "L3_Haplotype_Os06g0346300_RAP2_gene.csv"          
## [12] "L4_Haplotype_Os07g0623200_RAP2_gene.csv"          
## [13] "L4_Haplotype_Os07g0623501_RAP2_gene.csv"          
## [14] "L4_Haplotype_Os07g0623600_RAP2_gene.csv"          
## [15] "L4_Haplotype_Os07g43040_gene.csv"                 
## [16] "L6_Haplotype_Os11g0216000_RAB2_gene.csv"          
## [17] "Locus_1_final_file.csv"                           
## [18] "Locus_2_final_file.csv"                           
## [19] "Locus_3_final_file.csv"                           
## [20] "Locus_4_final_file.csv"                           
## [21] "Locus_6_final_file.csv"
temp <- read.csv("L2_Haplotype_Os05g0420900_RAP2_gene.csv")
head(temp)
haplotype2 <- temp[,c(1, 14)]

ultimate2 <- merge(pheno2, haplotype2, by=c("VCF_better"))
head(ultimate2)
unique(ultimate2$haplotype)
## [1] "AAAAAA" "AATAAA" "ATTAAA"
# let's see which haplotype is most popular among accessions:

ultimate2 %>% group_by(haplotype) %>%
  tally()

In this case - we have 6 original haplotypes - but only 1 of them is represented by 3 or more accessions

ultimate2$haplotype <- gsub("AAAAAA", "hap01", ultimate2$haplotype)
ultimate2$haplotype <- gsub("AATAAA", "hap02", ultimate2$haplotype)
ultimate2$haplotype <- gsub("ATTAAA", "hap03", ultimate2$haplotype)

ultimate2$haplotype <- factor(ultimate2$haplotype, levels = c("hap01", "hap02", "hap03"))
colnames(ultimate2)
##  [1] "VCF_better"        "HDRA.ID"           "NSFTV.ID"         
##  [4] "IRGC.ID"           "Subpop"            "Accession.Name"   
##  [7] "VCF_code"          "index.pots"        "Variety"          
## [10] "Shoot.area..cm2."  "Hull.area..cm2."   "Solidity"         
## [13] "Perimeter..cm."    "Nr.leaves"         "Nr.tillers"       
## [16] "Plant.height..cm." "Culm.height..cm."  "Leaf.length..cm." 
## [19] "DW..g."            "Leaf.angle..º."    "Tiller.angle..º." 
## [22] "Erectness..º."     "Droopiness..º."    "SUM_rank_traits"  
## [25] "SUM_norm_traits"   "norm_SI_rank"      "SI_rank"          
## [28] "haplotype"
Output <- TukeyHSD(aov(aov(Shoot.area..cm2. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Shoot.area..cm2., colour = haplotype)) 
my_box_plot <- my_box_plot + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot <- my_box_plot + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot <- my_box_plot + theme(legend.position = "none")
my_box_plot <- my_box_plot + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot <- my_box_plot + xlab("") + ylab(expression(paste("Shoot area (", cm^2, ")", sep = "")))
my_box_plot <- my_box_plot + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot <- my_box_plot + stat_pvalue_manual(test, label = "Tukey", y.position = 950)
my_box_plot

Output <- TukeyHSD(aov(aov(Hull.area..cm2. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot2 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Hull.area..cm2., colour = haplotype)) 
my_box_plot2 <- my_box_plot2 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot2 <- my_box_plot2 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot2 <- my_box_plot2 + theme(legend.position = "none")
my_box_plot2 <- my_box_plot2 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot2 <- my_box_plot2 + xlab("") + ylab(expression(paste("Hull area (", cm^2, ")", sep = "")))
my_box_plot2 <- my_box_plot2 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot2 <- my_box_plot2 + stat_pvalue_manual(test, label = "Tukey", y.position = 6000)
my_box_plot2

Output <- TukeyHSD(aov(aov(Solidity ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot3 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Solidity, colour = haplotype)) 
my_box_plot3 <- my_box_plot3 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot3 <- my_box_plot3 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot3 <- my_box_plot3 + theme(legend.position = "none")
my_box_plot3 <- my_box_plot3 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot3 <- my_box_plot3 + xlab("") + ylab("Solidity (a.u.)")
my_box_plot3 <- my_box_plot3 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot3 <- my_box_plot3 + stat_pvalue_manual(test, label = "Tukey", y.position = 0.21)
my_box_plot3

Output <- TukeyHSD(aov(aov(Perimeter..cm. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot4 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Perimeter..cm., colour = haplotype)) 
my_box_plot4 <- my_box_plot4 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot4 <- my_box_plot4 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot4 <- my_box_plot4 + theme(legend.position = "none")
my_box_plot4 <- my_box_plot4 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot4 <- my_box_plot4 + xlab("") + ylab("Perimeter (cm)")
my_box_plot4 <- my_box_plot4 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot4 <- my_box_plot4 + stat_pvalue_manual(test, label = "Tukey", y.position = 2500)
my_box_plot4

Output <- TukeyHSD(aov(aov(Nr.leaves ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot5 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Nr.leaves, colour = haplotype)) 
my_box_plot5 <- my_box_plot5 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot5 <- my_box_plot5 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot5 <- my_box_plot5 + theme(legend.position = "none")
my_box_plot5 <- my_box_plot5 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot5 <- my_box_plot5 + xlab("") + ylab("Leaf number")
my_box_plot5 <- my_box_plot5 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot5 <- my_box_plot5 + stat_pvalue_manual(test, label = "Tukey", y.position = 55)
my_box_plot5

Output <- TukeyHSD(aov(aov(Nr.tillers ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot6 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Nr.tillers, colour = haplotype)) 
my_box_plot6 <- my_box_plot6 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot6 <- my_box_plot6 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot6 <- my_box_plot6 + theme(legend.position = "none")
my_box_plot6 <- my_box_plot6 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot6 <- my_box_plot6 + xlab("") + ylab("Tiller number")
my_box_plot6 <- my_box_plot6 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot6 <- my_box_plot6 + stat_pvalue_manual(test, label = "Tukey", y.position = 17)
my_box_plot6

Output <- TukeyHSD(aov(aov(Plant.height..cm. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot7 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Plant.height..cm., colour = haplotype)) 
my_box_plot7 <- my_box_plot7 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot7 <- my_box_plot7 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot7 <- my_box_plot7 + theme(legend.position = "none")
my_box_plot7 <- my_box_plot7 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot7 <- my_box_plot7 + xlab("") + ylab("Plant height (cm)")
my_box_plot7 <- my_box_plot7 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot7 <- my_box_plot7 + stat_pvalue_manual(test, label = "Tukey", y.position = 95)
my_box_plot7

Output <- TukeyHSD(aov(aov(Culm.height..cm. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot8 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Culm.height..cm., colour = haplotype)) 
my_box_plot8 <- my_box_plot8 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot8 <- my_box_plot8 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot8 <- my_box_plot8 + theme(legend.position = "none")
my_box_plot8 <- my_box_plot8 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot8 <- my_box_plot8 + xlab("") + ylab("Culm height (cm)")
my_box_plot8 <- my_box_plot8 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot8 <- my_box_plot8 + stat_pvalue_manual(test, label = "Tukey", y.position = 35)
my_box_plot8

Output <- TukeyHSD(aov(aov(Leaf.length..cm. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot9 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Leaf.length..cm., colour = haplotype)) 
my_box_plot9 <- my_box_plot9 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot9 <- my_box_plot9 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot9 <- my_box_plot9 + theme(legend.position = "none")
my_box_plot9 <- my_box_plot9 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot9 <- my_box_plot9 + xlab("") + ylab("Leaf length (cm)")
my_box_plot9 <- my_box_plot9 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot9 <- my_box_plot9 + stat_pvalue_manual(test, label = "Tukey", y.position = 65)
my_box_plot9

Output <- TukeyHSD(aov(aov(DW..g. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot10 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = DW..g., colour = haplotype)) 
my_box_plot10 <- my_box_plot10 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot10 <- my_box_plot10 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot10 <- my_box_plot10 + theme(legend.position = "none")
my_box_plot10 <- my_box_plot10 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot10 <- my_box_plot10 + xlab("") + ylab("Dry weight (g)")
my_box_plot10 <- my_box_plot10 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot10 <- my_box_plot10 + stat_pvalue_manual(test, label = "Tukey", y.position = 4) 
my_box_plot10

Output <- TukeyHSD(aov(aov(Leaf.angle..º. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot11 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Leaf.angle..º., colour = haplotype)) 
my_box_plot11 <- my_box_plot11 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot11 <- my_box_plot11 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot11 <- my_box_plot11 + theme(legend.position = "none")
my_box_plot11 <- my_box_plot11 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot11 <- my_box_plot11 + xlab("") + ylab("Leaf angle (º)")
my_box_plot11 <- my_box_plot11 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot11 <- my_box_plot11 + stat_pvalue_manual(test, label = "Tukey", y.position = 65)
my_box_plot11

Output <- TukeyHSD(aov(aov(Tiller.angle..º. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot12 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Tiller.angle..º., colour = haplotype)) 
my_box_plot12 <- my_box_plot12 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot12 <- my_box_plot12 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot12 <- my_box_plot12 + theme(legend.position = "none")
my_box_plot12 <- my_box_plot12 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot12 <- my_box_plot12 + xlab("") + ylab("Tiller angle (º)")
my_box_plot12 <- my_box_plot12 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot12 <- my_box_plot12 + stat_pvalue_manual(test, label = "Tukey", y.position = 37)
my_box_plot12

Output <- TukeyHSD(aov(aov(Erectness..º. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot13 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Erectness..º., colour = haplotype)) 
my_box_plot13 <- my_box_plot13 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot13 <- my_box_plot13 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot13 <- my_box_plot13 + theme(legend.position = "none")
my_box_plot13 <- my_box_plot13 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot13 <- my_box_plot13 + xlab("") + ylab("Erectness (º)")
my_box_plot13 <- my_box_plot13 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot13 <- my_box_plot13 + stat_pvalue_manual(test, label = "Tukey", y.position = 165)
my_box_plot13

Output <- TukeyHSD(aov(aov(Droopiness..º. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot14 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Droopiness..º., colour = haplotype)) 
my_box_plot14 <- my_box_plot14 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot14 <- my_box_plot14 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot14 <- my_box_plot14 + theme(legend.position = "none")
my_box_plot14 <- my_box_plot14 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3"))
my_box_plot14 <- my_box_plot14 + xlab("") + ylab("Droopiness (º)")
my_box_plot14 <- my_box_plot14 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot14 <- my_box_plot14 + stat_pvalue_manual(test, label = "Tukey", y.position = 170)
my_box_plot14

pdf("Locus_02_Haplotypes_Os05g0420900_RAP2_NEW.pdf", height = 30, width = 12)
plot_grid(my_box_plot,  my_box_plot2, my_box_plot3, my_box_plot4, my_box_plot5, my_box_plot6, my_box_plot7, my_box_plot8, my_box_plot9, my_box_plot10, my_box_plot11, my_box_plot12, my_box_plot13, my_box_plot14,
          ncol=2,
          align = "hv", labels=c("AUTO"), 
          label_size = 24)
dev.off()
## quartz_off_screen 
##                 2

Locus 4 Os07g0623200

list.files(pattern=".csv")
##  [1] "L1_Haplotype_Os03g0841800_RAP2_gene_and_pheno.csv"
##  [2] "L1_Haplotype_Os03g0843700_RAP2_gene.csv"          
##  [3] "L1_Haplotype_Os03g0845000_RAP2_gene.csv"          
##  [4] "L1_Haplotype_Os03g0845700_RAP2_gene.csv"          
##  [5] "L1_Haplotype_Os03g0845800_RAP2_gene.csv"          
##  [6] "L1_Haplotype_Os03g0848700_RAP2_gene.csv"          
##  [7] "L1_Haplotype_Os03g62400_gene_and_pheno.csv"       
##  [8] "L2_Haplotype_Os05g0420500_RAP2_gene.csv"          
##  [9] "L2_Haplotype_Os05g0420600_RAP2_gene.csv"          
## [10] "L2_Haplotype_Os05g0420900_RAP2_gene.csv"          
## [11] "L3_Haplotype_Os06g0346300_RAP2_gene.csv"          
## [12] "L4_Haplotype_Os07g0623200_RAP2_gene.csv"          
## [13] "L4_Haplotype_Os07g0623501_RAP2_gene.csv"          
## [14] "L4_Haplotype_Os07g0623600_RAP2_gene.csv"          
## [15] "L4_Haplotype_Os07g43040_gene.csv"                 
## [16] "L6_Haplotype_Os11g0216000_RAB2_gene.csv"          
## [17] "Locus_1_final_file.csv"                           
## [18] "Locus_2_final_file.csv"                           
## [19] "Locus_3_final_file.csv"                           
## [20] "Locus_4_final_file.csv"                           
## [21] "Locus_6_final_file.csv"
temp <- read.csv("L4_Haplotype_Os07g0623200_RAP2_gene.csv")
head(temp)
haplotype2 <- temp[,c(1, 14)]

ultimate2 <- merge(pheno2, haplotype2, by=c("VCF_better"))
head(ultimate2)
unique(ultimate2$haplotype)
## [1] "AAAA" "ATAA" "ATAT" "AATA"
# let's see which haplotype is most popular among accessions:

ultimate2 %>% group_by(haplotype) %>%
  tally()

In this case - we have 6 original haplotypes - but only 1 of them is represented by 3 or more accessions

ultimate2$haplotype <- gsub("AAAA", "hap01", ultimate2$haplotype)
ultimate2$haplotype <- gsub("ATAA", "hap02", ultimate2$haplotype)
ultimate2$haplotype <- gsub("ATAT", "hap03", ultimate2$haplotype)
ultimate2$haplotype <- gsub("AATA", "hap04", ultimate2$haplotype)

ultimate2$haplotype <- factor(ultimate2$haplotype, levels = c("hap01", "hap02", "hap03", "hap04"))
colnames(ultimate2)
##  [1] "VCF_better"        "HDRA.ID"           "NSFTV.ID"         
##  [4] "IRGC.ID"           "Subpop"            "Accession.Name"   
##  [7] "VCF_code"          "index.pots"        "Variety"          
## [10] "Shoot.area..cm2."  "Hull.area..cm2."   "Solidity"         
## [13] "Perimeter..cm."    "Nr.leaves"         "Nr.tillers"       
## [16] "Plant.height..cm." "Culm.height..cm."  "Leaf.length..cm." 
## [19] "DW..g."            "Leaf.angle..º."    "Tiller.angle..º." 
## [22] "Erectness..º."     "Droopiness..º."    "SUM_rank_traits"  
## [25] "SUM_norm_traits"   "norm_SI_rank"      "SI_rank"          
## [28] "haplotype"
Output <- TukeyHSD(aov(aov(Shoot.area..cm2. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Shoot.area..cm2., colour = haplotype)) 
my_box_plot <- my_box_plot + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot <- my_box_plot + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot <- my_box_plot + theme(legend.position = "none")
my_box_plot <- my_box_plot + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3", "firebrick3"))
my_box_plot <- my_box_plot + xlab("") + ylab(expression(paste("Shoot area (", cm^2, ")", sep = "")))
my_box_plot <- my_box_plot + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot <- my_box_plot + stat_pvalue_manual(test, label = "Tukey", y.position = 950)
my_box_plot

Output <- TukeyHSD(aov(aov(Hull.area..cm2. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot2 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Hull.area..cm2., colour = haplotype)) 
my_box_plot2 <- my_box_plot2 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot2 <- my_box_plot2 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot2 <- my_box_plot2 + theme(legend.position = "none")
my_box_plot2 <- my_box_plot2 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3", "firebrick3"))
my_box_plot2 <- my_box_plot2 + xlab("") + ylab(expression(paste("Hull area (", cm^2, ")", sep = "")))
my_box_plot2 <- my_box_plot2 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot2 <- my_box_plot2 + stat_pvalue_manual(test, label = "Tukey", y.position = 6000)
my_box_plot2

Output <- TukeyHSD(aov(aov(Solidity ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot3 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Solidity, colour = haplotype)) 
my_box_plot3 <- my_box_plot3 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot3 <- my_box_plot3 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot3 <- my_box_plot3 + theme(legend.position = "none")
my_box_plot3 <- my_box_plot3 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3", "firebrick3"))
my_box_plot3 <- my_box_plot3 + xlab("") + ylab("Solidity (a.u.)")
my_box_plot3 <- my_box_plot3 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot3 <- my_box_plot3 + stat_pvalue_manual(test, label = "Tukey", y.position = 0.21)
my_box_plot3

Output <- TukeyHSD(aov(aov(Perimeter..cm. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot4 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Perimeter..cm., colour = haplotype)) 
my_box_plot4 <- my_box_plot4 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot4 <- my_box_plot4 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot4 <- my_box_plot4 + theme(legend.position = "none")
my_box_plot4 <- my_box_plot4 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3", "firebrick3"))
my_box_plot4 <- my_box_plot4 + xlab("") + ylab("Perimeter (cm)")
my_box_plot4 <- my_box_plot4 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot4 <- my_box_plot4 + stat_pvalue_manual(test, label = "Tukey", y.position = 2500)
my_box_plot4

Output <- TukeyHSD(aov(aov(Nr.leaves ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot5 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Nr.leaves, colour = haplotype)) 
my_box_plot5 <- my_box_plot5 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot5 <- my_box_plot5 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot5 <- my_box_plot5 + theme(legend.position = "none")
my_box_plot5 <- my_box_plot5 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3", "firebrick3"))
my_box_plot5 <- my_box_plot5 + xlab("") + ylab("Leaf number")
my_box_plot5 <- my_box_plot5 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot5 <- my_box_plot5 + stat_pvalue_manual(test, label = "Tukey", y.position = 55)
my_box_plot5

Output <- TukeyHSD(aov(aov(Nr.tillers ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot6 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Nr.tillers, colour = haplotype)) 
my_box_plot6 <- my_box_plot6 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot6 <- my_box_plot6 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot6 <- my_box_plot6 + theme(legend.position = "none")
my_box_plot6 <- my_box_plot6 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3", "firebrick3"))
my_box_plot6 <- my_box_plot6 + xlab("") + ylab("Tiller number")
my_box_plot6 <- my_box_plot6 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot6 <- my_box_plot6 + stat_pvalue_manual(test, label = "Tukey", y.position = 17)
my_box_plot6

Output <- TukeyHSD(aov(aov(Plant.height..cm. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot7 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Plant.height..cm., colour = haplotype)) 
my_box_plot7 <- my_box_plot7 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot7 <- my_box_plot7 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot7 <- my_box_plot7 + theme(legend.position = "none")
my_box_plot7 <- my_box_plot7 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3", "firebrick3"))
my_box_plot7 <- my_box_plot7 + xlab("") + ylab("Plant height (cm)")
my_box_plot7 <- my_box_plot7 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot7 <- my_box_plot7 + stat_pvalue_manual(test, label = "Tukey", y.position = 95)
my_box_plot7

Output <- TukeyHSD(aov(aov(Culm.height..cm. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot8 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Culm.height..cm., colour = haplotype)) 
my_box_plot8 <- my_box_plot8 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot8 <- my_box_plot8 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot8 <- my_box_plot8 + theme(legend.position = "none")
my_box_plot8 <- my_box_plot8 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3", "firebrick3"))
my_box_plot8 <- my_box_plot8 + xlab("") + ylab("Culm height (cm)")
my_box_plot8 <- my_box_plot8 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot8 <- my_box_plot8 + stat_pvalue_manual(test, label = "Tukey", y.position = 35)
my_box_plot8

Output <- TukeyHSD(aov(aov(Leaf.length..cm. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot9 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Leaf.length..cm., colour = haplotype)) 
my_box_plot9 <- my_box_plot9 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot9 <- my_box_plot9 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot9 <- my_box_plot9 + theme(legend.position = "none")
my_box_plot9 <- my_box_plot9 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3", "firebrick3"))
my_box_plot9 <- my_box_plot9 + xlab("") + ylab("Leaf length (cm)")
my_box_plot9 <- my_box_plot9 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot9 <- my_box_plot9 + stat_pvalue_manual(test, label = "Tukey", y.position = 65)
my_box_plot9

Output <- TukeyHSD(aov(aov(DW..g. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot10 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = DW..g., colour = haplotype)) 
my_box_plot10 <- my_box_plot10 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot10 <- my_box_plot10 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot10 <- my_box_plot10 + theme(legend.position = "none")
my_box_plot10 <- my_box_plot10 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3", "firebrick3"))
my_box_plot10 <- my_box_plot10 + xlab("") + ylab("Dry weight (g)")
my_box_plot10 <- my_box_plot10 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot10 <- my_box_plot10 + stat_pvalue_manual(test, label = "Tukey", y.position = 4) 
my_box_plot10

Output <- TukeyHSD(aov(aov(Leaf.angle..º. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot11 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Leaf.angle..º., colour = haplotype)) 
my_box_plot11 <- my_box_plot11 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot11 <- my_box_plot11 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot11 <- my_box_plot11 + theme(legend.position = "none")
my_box_plot11 <- my_box_plot11 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3", "firebrick3"))
my_box_plot11 <- my_box_plot11 + xlab("") + ylab("Leaf angle (º)")
my_box_plot11 <- my_box_plot11 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot11 <- my_box_plot11 + stat_pvalue_manual(test, label = "Tukey", y.position = 65)
my_box_plot11

Output <- TukeyHSD(aov(aov(Tiller.angle..º. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot12 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Tiller.angle..º., colour = haplotype)) 
my_box_plot12 <- my_box_plot12 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot12 <- my_box_plot12 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot12 <- my_box_plot12 + theme(legend.position = "none")
my_box_plot12 <- my_box_plot12 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3", "firebrick3"))
my_box_plot12 <- my_box_plot12 + xlab("") + ylab("Tiller angle (º)")
my_box_plot12 <- my_box_plot12 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot12 <- my_box_plot12 + stat_pvalue_manual(test, label = "Tukey", y.position = 37)
my_box_plot12

Output <- TukeyHSD(aov(aov(Erectness..º. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot13 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Erectness..º., colour = haplotype)) 
my_box_plot13 <- my_box_plot13 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot13 <- my_box_plot13 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot13 <- my_box_plot13 + theme(legend.position = "none")
my_box_plot13 <- my_box_plot13 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3", "firebrick3"))
my_box_plot13 <- my_box_plot13 + xlab("") + ylab("Erectness (º)")
my_box_plot13 <- my_box_plot13 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot13 <- my_box_plot13 + stat_pvalue_manual(test, label = "Tukey", y.position = 165)
my_box_plot13

Output <- TukeyHSD(aov(aov(Droopiness..º. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

my_box_plot14 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Droopiness..º., colour = haplotype)) 
my_box_plot14 <- my_box_plot14 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot14 <- my_box_plot14 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot14 <- my_box_plot14 + theme(legend.position = "none")
my_box_plot14 <- my_box_plot14 + scale_colour_manual(values = c("steelblue", "firebrick3", "firebrick3", "firebrick3"))
my_box_plot14 <- my_box_plot14 + xlab("") + ylab("Droopiness (º)")
my_box_plot14 <- my_box_plot14 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot14 <- my_box_plot14 + stat_pvalue_manual(test, label = "Tukey", y.position = 170)
my_box_plot14

pdf("Locus_04_Haplotypes_Os07g0623200_RAP2_NEW.pdf", height = 30, width = 12)
plot_grid(my_box_plot,  my_box_plot2, my_box_plot3, my_box_plot4, my_box_plot5, my_box_plot6, my_box_plot7, my_box_plot8, my_box_plot9, my_box_plot10, my_box_plot11, my_box_plot12, my_box_plot13, my_box_plot14,
          ncol=2,
          align = "hv", labels=c("AUTO"), 
          label_size = 24)
dev.off()
## quartz_off_screen 
##                 2

Locus 6 Os11g0216000

list.files(pattern=".csv")
##  [1] "L1_Haplotype_Os03g0841800_RAP2_gene_and_pheno.csv"
##  [2] "L1_Haplotype_Os03g0843700_RAP2_gene.csv"          
##  [3] "L1_Haplotype_Os03g0845000_RAP2_gene.csv"          
##  [4] "L1_Haplotype_Os03g0845700_RAP2_gene.csv"          
##  [5] "L1_Haplotype_Os03g0845800_RAP2_gene.csv"          
##  [6] "L1_Haplotype_Os03g0848700_RAP2_gene.csv"          
##  [7] "L1_Haplotype_Os03g62400_gene_and_pheno.csv"       
##  [8] "L2_Haplotype_Os05g0420500_RAP2_gene.csv"          
##  [9] "L2_Haplotype_Os05g0420600_RAP2_gene.csv"          
## [10] "L2_Haplotype_Os05g0420900_RAP2_gene.csv"          
## [11] "L3_Haplotype_Os06g0346300_RAP2_gene.csv"          
## [12] "L4_Haplotype_Os07g0623200_RAP2_gene.csv"          
## [13] "L4_Haplotype_Os07g0623501_RAP2_gene.csv"          
## [14] "L4_Haplotype_Os07g0623600_RAP2_gene.csv"          
## [15] "L4_Haplotype_Os07g43040_gene.csv"                 
## [16] "L6_Haplotype_Os11g0216000_RAB2_gene.csv"          
## [17] "Locus_1_final_file.csv"                           
## [18] "Locus_2_final_file.csv"                           
## [19] "Locus_3_final_file.csv"                           
## [20] "Locus_4_final_file.csv"                           
## [21] "Locus_6_final_file.csv"
temp <- read.csv("L6_Haplotype_Os11g0216000_RAB2_gene.csv")
head(temp)
haplotype2 <- temp[,c(1, 14)]

ultimate2 <- merge(pheno2, haplotype2, by=c("VCF_better"))
head(ultimate2)
unique(ultimate2$haplotype)
## [1] "AAAAAAAAA" "AAAATAAAA"
# let's see which haplotype is most popular among accessions:

ultimate2 %>% group_by(haplotype) %>%
  tally()
ultimate2$haplotype <- gsub("AAAAAAAAA", "hap01", ultimate2$haplotype)
ultimate2$haplotype <- gsub("AAAATAAAA", "hap02", ultimate2$haplotype)

ultimate2$haplotype <- factor(ultimate2$haplotype, levels = c("hap01", "hap02"))
colnames(ultimate2)
##  [1] "VCF_better"        "HDRA.ID"           "NSFTV.ID"         
##  [4] "IRGC.ID"           "Subpop"            "Accession.Name"   
##  [7] "VCF_code"          "index.pots"        "Variety"          
## [10] "Shoot.area..cm2."  "Hull.area..cm2."   "Solidity"         
## [13] "Perimeter..cm."    "Nr.leaves"         "Nr.tillers"       
## [16] "Plant.height..cm." "Culm.height..cm."  "Leaf.length..cm." 
## [19] "DW..g."            "Leaf.angle..º."    "Tiller.angle..º." 
## [22] "Erectness..º."     "Droopiness..º."    "SUM_rank_traits"  
## [25] "SUM_norm_traits"   "norm_SI_rank"      "SI_rank"          
## [28] "haplotype"

OK - now let’s graph them:

Output <- TukeyHSD(aov(aov(Shoot.area..cm2. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
P4
## [1] 3.952704e-09
my_box_plot <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Shoot.area..cm2., colour = haplotype)) 
my_box_plot <- my_box_plot + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot <- my_box_plot + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot <- my_box_plot + theme(legend.position = "none")
my_box_plot <- my_box_plot + scale_colour_manual(values = c("steelblue", "firebrick3"))
my_box_plot <- my_box_plot + xlab("") + ylab(expression(paste("Shoot area (", cm^2, ")", sep = "")))
my_box_plot <- my_box_plot + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot <- my_box_plot + stat_pvalue_manual(test2, label = "Tukey", y.position = 950)
my_box_plot

Output <- TukeyHSD(aov(aov(Hull.area..cm2. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
P4
## [1] 7.073674e-08
my_box_plot2 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Hull.area..cm2., colour = haplotype)) 
my_box_plot2 <- my_box_plot2 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot2 <- my_box_plot2 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot2 <- my_box_plot2 + theme(legend.position = "none")
my_box_plot2 <- my_box_plot2 + scale_colour_manual(values = c("steelblue", "firebrick3"))
my_box_plot2 <- my_box_plot2 + xlab("") + ylab(expression(paste("Hull area (", cm^2, ")", sep = "")))
my_box_plot2 <- my_box_plot2 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot2 <- my_box_plot2 + stat_pvalue_manual(test2, label = "Tukey", y.position = 6000)
my_box_plot2

Output <- TukeyHSD(aov(aov(Solidity ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
P4
## [1] 0.0452482
my_box_plot3 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Solidity, colour = haplotype)) 
my_box_plot3 <- my_box_plot3 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot3 <- my_box_plot3 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot3 <- my_box_plot3 + theme(legend.position = "none")
my_box_plot3 <- my_box_plot3 + scale_colour_manual(values = c("steelblue", "firebrick3"))
my_box_plot3 <- my_box_plot3 + xlab("") + ylab("Solidity (a.u.)")
my_box_plot3 <- my_box_plot3 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot3 <- my_box_plot3 + stat_pvalue_manual(test2, label = "Tukey", y.position = 0.21)
my_box_plot3

Output <- TukeyHSD(aov(aov(Perimeter..cm. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
P4
## [1] 1.671474e-11
my_box_plot4 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Perimeter..cm., colour = haplotype)) 
my_box_plot4 <- my_box_plot4 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot4 <- my_box_plot4 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot4 <- my_box_plot4 + theme(legend.position = "none")
my_box_plot4 <- my_box_plot4 + scale_colour_manual(values = c("steelblue", "firebrick3"))
my_box_plot4 <- my_box_plot4 + xlab("") + ylab("Perimeter (cm)")
my_box_plot4 <- my_box_plot4 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot4 <- my_box_plot4 + stat_pvalue_manual(test2, label = "Tukey", y.position = 2500)
my_box_plot4

Output <- TukeyHSD(aov(aov(Nr.leaves ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
P4
## [1] 0
my_box_plot5 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Nr.leaves, colour = haplotype)) 
my_box_plot5 <- my_box_plot5 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot5 <- my_box_plot5 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot5 <- my_box_plot5 + theme(legend.position = "none")
my_box_plot5 <- my_box_plot5 + scale_colour_manual(values = c("steelblue", "firebrick3"))
my_box_plot5 <- my_box_plot5 + xlab("") + ylab("Leaf number")
my_box_plot5 <- my_box_plot5 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot5 <- my_box_plot5 + stat_pvalue_manual(test2, label = "Tukey", y.position = 55)
my_box_plot5

Output <- TukeyHSD(aov(aov(Nr.tillers ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
P4
## [1] 0
my_box_plot6 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Nr.tillers, colour = haplotype)) 
my_box_plot6 <- my_box_plot6 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot6 <- my_box_plot6 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot6 <- my_box_plot6 + theme(legend.position = "none")
my_box_plot6 <- my_box_plot6 + scale_colour_manual(values = c("steelblue", "firebrick3"))
my_box_plot6 <- my_box_plot6 + xlab("") + ylab("Tiller number")
my_box_plot6 <- my_box_plot6 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot6 <- my_box_plot6 + stat_pvalue_manual(test2, label = "Tukey", y.position = 17)
my_box_plot6

Output <- TukeyHSD(aov(aov(Plant.height..cm. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
P4
## [1] 0.5464665
my_box_plot7 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Plant.height..cm., colour = haplotype)) 
my_box_plot7 <- my_box_plot7 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot7 <- my_box_plot7 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot7 <- my_box_plot7 + theme(legend.position = "none")
my_box_plot7 <- my_box_plot7 + scale_colour_manual(values = c("steelblue", "firebrick3"))
my_box_plot7 <- my_box_plot7 + xlab("") + ylab("Plant height (cm)")
my_box_plot7 <- my_box_plot7 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot7 <- my_box_plot7 + stat_pvalue_manual(test3, label = "Tukey", y.position = 95)
my_box_plot7

Output <- TukeyHSD(aov(aov(Culm.height..cm. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
P4
## [1] 0.01020817
my_box_plot8 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Culm.height..cm., colour = haplotype)) 
my_box_plot8 <- my_box_plot8 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot8 <- my_box_plot8 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot8 <- my_box_plot8 + theme(legend.position = "none")
my_box_plot8 <- my_box_plot8 + scale_colour_manual(values = c("steelblue", "firebrick3"))
my_box_plot8 <- my_box_plot8 + xlab("") + ylab("Culm height (cm)")
my_box_plot8 <- my_box_plot8 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot8 <- my_box_plot8 + stat_pvalue_manual(test2, label = "Tukey", y.position = 35)
my_box_plot8

Output <- TukeyHSD(aov(aov(Leaf.length..cm. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
P4
## [1] 0.9102221
my_box_plot9 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Leaf.length..cm., colour = haplotype)) 
my_box_plot9 <- my_box_plot9 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot9 <- my_box_plot9 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot9 <- my_box_plot9 + theme(legend.position = "none")
my_box_plot9 <- my_box_plot9 + scale_colour_manual(values = c("steelblue", "firebrick3"))
my_box_plot9 <- my_box_plot9 + xlab("") + ylab("Leaf length (cm)")
my_box_plot9 <- my_box_plot9 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot9 <- my_box_plot9 + stat_pvalue_manual(test3, label = "Tukey", y.position = 65)
my_box_plot9

Output <- TukeyHSD(aov(aov(DW..g. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
P4
## [1] 1.336687e-08
my_box_plot10 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = DW..g., colour = haplotype)) 
my_box_plot10 <- my_box_plot10 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot10 <- my_box_plot10 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot10 <- my_box_plot10 + theme(legend.position = "none")
my_box_plot10 <- my_box_plot10 + scale_colour_manual(values = c("steelblue", "firebrick3"))
my_box_plot10 <- my_box_plot10 + xlab("") + ylab("Dry weight (g)")
my_box_plot10 <- my_box_plot10 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot10 <- my_box_plot10 + stat_pvalue_manual(test2, label = "Tukey", y.position = 4) 
my_box_plot10

Output <- TukeyHSD(aov(aov(Leaf.angle..º. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
P4
## [1] 0.6813148
my_box_plot11 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Leaf.angle..º., colour = haplotype)) 
my_box_plot11 <- my_box_plot11 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot11 <- my_box_plot11 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot11 <- my_box_plot11 + theme(legend.position = "none")
my_box_plot11 <- my_box_plot11 + scale_colour_manual(values = c("steelblue", "firebrick3"))
my_box_plot11 <- my_box_plot11 + xlab("") + ylab("Leaf angle (º)")
my_box_plot11 <- my_box_plot11 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot11 <- my_box_plot11 + stat_pvalue_manual(test3, label = "Tukey", y.position = 65)
my_box_plot11

Output <- TukeyHSD(aov(aov(Tiller.angle..º. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
P4
## [1] 2.522885e-05
my_box_plot12 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Tiller.angle..º., colour = haplotype)) 
my_box_plot12 <- my_box_plot12 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot12 <- my_box_plot12 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot12 <- my_box_plot12 + theme(legend.position = "none")
my_box_plot12 <- my_box_plot12 + scale_colour_manual(values = c("steelblue", "firebrick3"))
my_box_plot12 <- my_box_plot12 + xlab("") + ylab("Tiller angle (º)")
my_box_plot12 <- my_box_plot12 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot12 <- my_box_plot12 + stat_pvalue_manual(test2, label = "Tukey", y.position = 37)
my_box_plot12

Output <- TukeyHSD(aov(aov(Erectness..º. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
P4
## [1] 0.2157439
my_box_plot13 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Erectness..º., colour = haplotype)) 
my_box_plot13 <- my_box_plot13 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot13 <- my_box_plot13 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot13 <- my_box_plot13 + theme(legend.position = "none")
my_box_plot13 <- my_box_plot13 + scale_colour_manual(values = c("steelblue", "firebrick3"))
my_box_plot13 <- my_box_plot13 + xlab("") + ylab("Erectness (º)")
my_box_plot13 <- my_box_plot13 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot13 <- my_box_plot13 + stat_pvalue_manual(test3, label = "Tukey", y.position = 165)
my_box_plot13

Output <- TukeyHSD(aov(aov(Droopiness..º. ~ haplotype, data = ultimate2)))
P4 = Output$haplotype[,'p adj']
P4
## [1] 0.2157439
my_box_plot14 <- ggplot(data = ultimate2, mapping = aes(x = haplotype, y = Droopiness..º., colour = haplotype)) 
my_box_plot14 <- my_box_plot14 + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot14 <- my_box_plot14 + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: `fun.y` is deprecated. Use `fun` instead.
my_box_plot14 <- my_box_plot14 + theme(legend.position = "none")
my_box_plot14 <- my_box_plot14 + scale_colour_manual(values = c("steelblue", "firebrick3"))
my_box_plot14 <- my_box_plot14 + xlab("") + ylab("Droopiness (º)")
my_box_plot14 <- my_box_plot14 + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
my_box_plot14 <- my_box_plot14 + stat_pvalue_manual(test3, label = "Tukey", y.position = 170)
my_box_plot14

pdf("Locus_06_Haplotypes_Os11g216000_RAP2_NEW.pdf", height = 30, width = 12)
plot_grid(my_box_plot,  my_box_plot2, my_box_plot3, my_box_plot4, my_box_plot5, my_box_plot6, 
          my_box_plot7, my_box_plot8, my_box_plot9, my_box_plot10, my_box_plot11, my_box_plot12, 
          my_box_plot13, my_box_plot14,
          ncol=2,
          align = "hv", labels=c("AUTO"), 
          label_size = 24)
dev.off()
## quartz_off_screen 
##                 2