Background
These R code chunks are related to the following research article which has been published recently. The article can be found via the following link: manuscript link. As the article is sufficiently detailed about the objective and results of the study, this will only serve as a tutorial to perform the analysis and draw the graphs included in the article. Authors believe in reproducible research and this is their way of contributing to the current knowledge. If this code helps you in any way we kindly ask you to cite the aforementioned article (that’s the best way of saying thanks). All regression codes will not be made available in this tutorial because the authors will give a full commented tutorial on the topic.
Prerequisites
We first load all packages and data sets that we will be needing for the analysis. The full data sets will be provided on a reasonable request. However, we will use the head() function to give you an idea about the data structure.
##load data
##======= regression =======##
regre <- read.csv("D:/Chris/ETUDES1/PhD thesis/R files/DATABASE (excel)/in ovo 5th trial/RTs model/regre.csv")
regre1 <- read.csv("D:/Chris/ETUDES1/PhD thesis/R files/DATABASE (excel)/in ovo 5th trial/RTs model/regre1.csv")
head(regre)## days rt trt
## 1 0 39.3 CON
## 2 0 40.0 CON
## 3 0 40.3 CON
## 4 0 40.0 CON
## 5 0 39.7 CON
## 6 0 39.7 CON+HShead(regre1)## days rt trt
## 1 1 41.7 CON
## 2 1 41.7 CON
## 3 1 41.3 CON
## 4 1 41.3 CON
## 5 1 40.8 CON
## 6 1 42.5 CON+HS##======= heat map =======##
rt1<- read.csv("D:/Chris/ETUDES1/PhD thesis/R files/DATABASE (excel)/in ovo 5th trial/RTs model/heat2.csv",row.names = 1)
rt2<- read.csv("D:/Chris/ETUDES1/PhD thesis/R files/DATABASE (excel)/in ovo 5th trial/RTs model/heat3.csv",row.names = 1)
head(rt1)## CT0 CT4 CT8 CT12 CT16 CT24 CT28 HS1 HS3 HS5 HS7
## CON1 39.3 39.8 41.7 41.4 41.5 41.6 41.3 41.7 41.3 41.6 41.2
## CON2 40.0 41.3 41.8 41.4 41.6 42.0 42.2 41.7 41.5 41.5 41.3
## CON3 40.3 40.7 41.5 41.3 41.2 41.3 41.3 41.3 41.5 41.5 41.3
## CON4 40.0 41.1 41.3 41.4 41.6 41.5 41.6 41.3 41.8 41.7 41.6
## CON5 39.7 41.0 41.1 41.3 41.2 41.1 41.2 40.8 41.3 40.8 41.3
## CON+HS1 39.7 41.0 41.5 41.6 41.7 41.6 41.7 42.5 42.4 42.7 42.5head(rt2)## Temperature Challenge
## CON1 CON NT
## CON2 CON NT
## CON3 CON NT
## CON4 CON NT
## CON5 CON NT
## CON+HS1 CON+HS HS##======= PCA =======##
allpca <- read.csv("D:/Chris/ETUDES1/PhD thesis/R files/DATABASE (excel)/in ovo 5th trial/RTs model/heat.csv",row.names = 1)
head(allpca)## CT0 CT4 CT8 CT12 CT16 CT24 CT28 HS1 HS3 HS5 HS7 Treatment
## CON1 39.3 39.8 41.7 41.4 41.5 41.6 41.3 41.7 41.3 41.6 41.2 CON
## CON2 40.0 41.3 41.8 41.4 41.6 42.0 42.2 41.7 41.5 41.5 41.3 CON
## CON3 40.3 40.7 41.5 41.3 41.2 41.3 41.3 41.3 41.5 41.5 41.3 CON
## CON4 40.0 41.1 41.3 41.4 41.6 41.5 41.6 41.3 41.8 41.7 41.6 CON
## CON5 39.7 41.0 41.1 41.3 41.2 41.1 41.2 40.8 41.3 40.8 41.3 CON
## CON+HS1 39.7 41.0 41.5 41.6 41.7 41.6 41.7 42.5 42.4 42.7 42.5 CON+HS
## Challenge
## CON1 NT
## CON2 NT
## CON3 NT
## CON4 NT
## CON5 NT
## CON+HS1 HS##======= survival analysis =======##
sdata <- read.csv("D:/Chris/ETUDES1/PhD thesis/R files/DATABASE (excel)/in ovo 5th trial/chronic_2.csv")
head(sdata)## Treatment id time status Temperature Humidity Bodyweight THI
## 1 CON+HS B1 358 0 22.6 65 1495 28.96
## 2 CON+HS B2 301 0 22.6 65 1140 28.96
## 3 CON+HS B3 358 0 22.6 65 1765 28.96
## 4 CON+HS B4 358 0 22.6 65 1600 28.96
## 5 CON+HS B5 284 1 28.6 39 1635 30.16
## 6 CON+HS B7 358 0 22.6 65 1375 28.96##load packages
library(tidyverse) ## for data wrangling and ggplot2 plots
library(FactoMineR) ## for performing principal component analysis
library(factoextra) ## for improving PCA plots
library(survminer) ## for performing survival analysis
library(survival) ## for performing survival analysis
library(patchwork) ## to combine plots together
library(ComplexHeatmap) ## for the heatmap
library(circlize) ## for the color palette in the heat mapFigure: regression plots
This figure is a plot describing the association between age and rectal temperature of broilers reared under both thermoneutral and heat stress conditions. we previously conducted polynomial regression analysis and detected that the quadratic trend was a better fit for the data. we will not discuss the regression analysis here since it will be a whole topic of another code.
##define a color palette for the graphs
Mycol <- c("CON" = "black", "CON+HS" = "red" , "G10+HS" = "green", "TM+HS" = "blue", "G10+TM+HS" = "coral")
Mycol2 <- c("CON+HS" = "red" , "G10+HS" = "green", "TM+HS" = "blue", "G10+TM+HS" = "coral")
##plot the thermoneutral figure
con_plot <- ggplot(regre)+ aes(x=days, y = rt)+
geom_jitter(aes(x=days, y=rt, color = trt))+
stat_smooth(aes(y = rt), method = "lm", formula = y ~ x + I(x^2), color = "black")+ #### quadratic regression
labs(color = "Treatment", x = "Rearing days", y = "Rectal temperature (°C)")+ ###legend and axis name
scale_colour_manual(values = Mycol)+ ### define custom colors
ylim(38, 43)+
theme_classic()
##plot the heat stress figure
hs_temp <- regre1 %>%
filter(trt != "CON") ## exclude CON group
hs_plot <- ggplot(hs_temp)+ aes(x=days, y = rt)+
geom_jitter(aes(x=days, y=rt, color = trt))+
stat_smooth(aes(y = rt), method = "lm", formula = y ~ x + I(x^2), color = "black")+ #### quadratic regression
labs(color = "Treatment", x = "Heat stress days", y = "Rectal temperature (°C)")+ ###legend and axis name
scale_colour_manual(values = Mycol2)+ ### define custom colors
ylim(41, 46)+
theme_classic()
##combine both plot in a single one
(con_plot/ (hs_plot + theme (legend.position = "none"))) + ## put the con_plot above the hs_plot
plot_annotation(tag_levels = 'a') + ##adding a, and b as plot indices
plot_layout(guides = "collect") ## combine legends in one
Figure: clustered heatmap
This figure is a clustered heat map showing the overall trend of the broiler’s rectal temperature during the whole trial. we will use the K-mean algorithm to specify the number of the cluster we would like to separate rectal temperatures into. we will minimize the code to its strict minimum, feel free to reach out if you have any further questions.
## transform data frame in matrix
rtonly1<-as.matrix(rt1)
## transpose matrix
rtonly1<- t(rtonly1)
##colors with circlize (values 38 to 46 are max and min value from dataset)
col_fun = colorRampPalette(c("Yellow" , "orange", "red"))(100)
##create annotation on the to pf the heat map
#### if you have more than one annotation example of challenge and treatment ###
ha2 = HeatmapAnnotation(Treatment = rt2$Temperature,
Challenge = rt2$Challenge,
col = list(Treatment =c("CON"="black","CON+HS"="red","G10+HS"="green","TM+HS"="blue", "G10+TM+HS"="orange" )
, Challenge = c("NT" = "gray", "HS" = "purple")))
##plot the heat map with custom annotations
Heatmap(rtonly1, name = "RT (°C)", #name of the values drawn
col = col_fun, #define color palette for values
rect_gp = gpar(col = "white", lwd = 2), #separate heat map tiles with a blank space
column_dend_height = unit(2, "cm"), #adjust columns height
row_dend_width = unit(3, "cm"),
row_km =3, #split the heat map in 3 clusters with the k-means algorithm
top_annotation = ha2 ) #add the annotation on the top of the heat map
Figure: principal component analysis
This figure is a PCA plot showing patterns in the broiler’s rectal temperature during the whole trial. we will use rectal temperatures as active variables and metadata as supplementary variables.
res <- PCA(allpca, quali.sup = 12:13, graph = F) ## performing the PCA (provide standard plots)
#plot individual tables (with both treatment and challenge as supplementary variables)
plot1 <- fviz_pca_ind(res, axes = c(1, 2),
label = "none", # hide individual labels
habillage = "Challenge",# color by groups
palette = c("red","blue"),
repel = TRUE, # avoid text overlap
addEllipses = TRUE, ellipse.type = "confidence",
ggtheme = theme_minimal())
plot2 <- fviz_pca_ind(res, axes = c(1, 2),
label = "none", # hide individual labels
habillage = "Treatment",# color by groups
palette = c("black","red","green","coral","blue"),
repel = TRUE,
addEllipses = TRUE, ellipse.type = "confidence",
ggtheme = theme_minimal())
#plot variable tables with both
plot3 <- fviz_pca_var(res, col.var="contrib", repel = TRUE) +
scale_color_gradient2(low="dark blue", mid="blue",
high="red", midpoint=7) +
theme_minimal()
plot4 <- fviz_pca_var(res, col.var="cos2", repel = TRUE) +
scale_color_gradient2(low="dark blue", mid="blue",
high="red", midpoint=0.5) + theme_minimal()
#combine all plots
(plot1+plot3)/(plot2+plot4) +
plot_layout(guides = 'collect') +
plot_annotation(tag_levels = 'a')
Figure: survival curves
This figure will present the Kaplan-Meir survival curves and risk table of broilers exposed to heat stress.
#we fit the Kaplan-meir survival model
fit2 <- survfit(Surv(time, status) ~ Treatment, data = sdata)
#plot the survival curve and risk table
ggsurvplot(fit2, pval = FALSE, ##remove p value
risk.table.col = "strata" , #color the groups differently (in the bottom table)
xlab="Time (min)",
break.x.by = 50,
risk.table=TRUE, #add the bottom table
ggtheme = theme_bw(),
font.subtitle=12, font.legend=12,
palette = Mycol2 , #use my custom predefined colors
risk.table.fontsize=4,
risk.table.y.text = FALSE,
legend.title = "Treatments",
legend="right",
legend.labs = c("CON+HS", "G10+HS","G10+TM+HS","TM+HS"), #reorder legend
data=sdata)