#Treatment Effect Analysis #Two tables were created for each gene one containing the log-transformed DCt value for the control fish and the immune stimulus fish, family, population, tank, weight, and experimental chip; and one containing the same parameters but for the control fish and the handling stress fish. #The analysis below was repeated for each gene per treatment and the p-values were collected, entered into a new file, and an FDR correction was completed for each treatment (all immune simulus reponses and all stress treatment responses). #The pvalues for all genes within a treatment were combined for the FDR correction. >library(lmer) >library(lmerTest) >CAL_Immune_Treatment <- read.table(file.choose().header=T,sep="\t") >CAL_Immune_Treatment.anova<- lmer(DCt.Log ~ Population*Treatment + (1|Family) + (1|Tank) + (1|Weight) + (1|Chip), REML = T, data = CAL_Immune_Treatment) >step(CAL_Immune_Treatment.anova) >CAL_Stress_Treatment <- read.table(file.choose().header=T,sep="\t") >CAL_Stress_Treatment.anova<- lmer(DCt.Log ~ Population*Treatment + (1|Family) + (1|Tank) + (1|Weight) + (1|Chip), REML = T, data = CAL_Stress_Treatment) >step(CAL_Immune_Treatment.anova) #The resulting p-values were recorded and used in the FDR correction below. One table was made for each treatment, containing all pvalues from the Treatment Effect Analysis. #FDR Correction #The tables were first ordered according to p-value from smallest to largest, then the Benjamini-Hochberg correction was applied and added to a new "FDR" column in the data table. >Immune_Treatment_Pvalues <- read.table(file.choose().header=T,sep="\t") >Stress_Treatment_Pvalues <- read.table(file.choose().header=T,sep="\t") >Immune_Treatment_Pvalues <- Immune_Treatment_Pvalues[order(Immune_Treatment_Pvalues$Pvalue),] >Immune_Treatment_Pvalues$FDR <- p.adjust(Immune_Treatment_Pvalues$Pvalue, method="BH") >Stress_Treatment_Pvalues <- Stress_Treatment_Pvalues[order(Stress_Treatment_Pvalues$Pvalue),] >Stress_Treatment_Pvalues$FDR <- p.adjust(Stress_Treatment_Pvalues$Pvalue, method="BH") #Effect Size >cohen.d(CAL_Immune_Treatment$DCt.Log, CAL_Immune_Treatment$Treatment, pooled=TRUE, paired=FALSE, na.rm=TRUE, hedges.correction=FALSE, conf.level=0.95, noncentral=FALSE) >cohen.d(CAL_Stress_Treatment$DCt.Log, CAL_Stress_Treatment$Treatment, pooled=TRUE, paired=FALSE, na.rm=TRUE, hedges.correction=FALSE, conf.level=0.95, noncentral=FALSE) #Population and Sire Effect Test #New tables were created using only one gene per treatment and the analysis was repeated for every gene at-rest (control group) and every significant treatment response gene from the Treatment Effect Analysis. #The code for hsp70 is used below as it was the only gene to exhibit a significant treatment effect to both the immune stimulus and the handling stress. >hsp70_Control <- read.table(file.choose(),header=T,sep="\t") >hsp70_Immune <- read.table(file.choose(),header=T,sep="\t") >hsp70_Stress <- read.table(file.choose(),header=T,sep="\t") >hsp70_Control.anova <- lmer(DCt.Log ~ Population + (1|Family) + (1|Tank) + (1|Weight) + (1|Chip), REML = T, data = hsp70_Control) >step(hsp70_Control.anova) >hsp70_Immune.anova <- lmer(DCt.Log ~ Population + (1|Family) + (1|Tank) + (1|Weight) + (1|Chip), REML = T, data = hsp70_Immune) >step(hsp70_Immune.anova) >hsp70_Stress.anova <- lmer(DCt.Log ~ Population + (1|Family) + (1|Tank) + (1|Weight) + (1|Chip), REML = T, data = hsp70_Stress) >step(hsp70_Stress.anova) #Tukey test for N-KEF population test (the only gene to demontrate a significant population effect). >library(multcomp) >NKEF_Control_glht <- glht(NKEF_Control.anova, linfct=mcp(Population="Tukey") >confint(NKEF_Control_glht) >cld(NKEF_Control_glht) #Population Effect per gene group (i.e. immune, growth, metabolic, and stress) #All genes for each gene group were combined by treatment and tested for a significant population effect. #The below anlaysis was repeated for the growth genes, metabolic genes, and stress genes. >Immune_Genes_Control <- read.table(file.choose(),header=T,sep="\t") >Immune_Genes_Immune_Treatment <- read.table(file.choose(),header=T,sep="\t") >Immune_Genes_Stress_Treatment <- read.table(file.choose(),header=T,sep="\t") >Immune_Genes_Control.anova <- lmet(DCt.Log ~ Population (1|Family) + (1|Tank) + (1|Weight) + (1|Target.Name), REML = T, data = Immune_Genes_Control) >step(Immune_Genes_Control.anova) >Immune_Genes_Immune_Treatment.anova <- lmet(DDCt.Log ~ Population (1|Family) + (1|Tank) + (1|Weight) + (1|Target.Name), REML = T, data = Immune_Genes_Immune_Treatment) >step(Immune_Genes_Immune_Treatment.anova) >Immune_Genes_Stress_Treatment.anova <- lmet(DDCt.Log ~ Population (1|Family) + (1|Tank) + (1|Weight) + (1|Target.Name), REML = T, data = Immune_Genes_Stress_Treatment) >step(Immune_Genes_Stress_Treatment.anova) #Pincipal Component Analysis #This anlaysis was performed for each gene group at rest and in response to the immune stimulus. The stress genes were the only gene group to exhibit a significant stress challenge effect, therefore a PCA was also done for stress genes in response to a handling stress. #Tables were created with only the log-transformed DCt or log-transformed DDCt for each gene group for each treament, with no other factors. >library(missMDA) >Stress_Genes_Control_PCA <- read.table(file.choose(),header=T,sep="\t") >Stress_Genes_ImmTrtt_PCA <- read.table(file.choose(),header=T,sep="\t") >Stress_Genes_StrTrt_PCA <- read.table(file.choose(),header=T,sep="\t") >estim_ncpPCA(Stress_Genes_Control_PCA) >estim_ncpPCA(Stress_Genes_ImmTrt_PCA) >estim_ncpPCA(Stress_Genes_StrTrt_PCA) #In the below commands, the ncp was equal to the number generated by estim_nspPCA, therefore this varied for every PCA. >Stress_Genes_Control_PCA.impute <- imputePCA(data, ncp=(results_from_estim_ncpPCA), scale = FALSE) >Stress_Genes_ImmTrt_PCA.impute <- imputePCA(data, ncp=(results_from_estim_ncpPCA), scale = FALSE) >Stress_Genes_StrTRt_PCA.impute <- imputePCA(data, ncp=(results_from_estim_ncpPCA), scale = FALSE) >Stress_Genes_Control_PCA.impute.pca <- PCA(Stress_Genes_Control_PCA.impute$completeObs) >Stress_Genes_ImmTrt_PCA.impute.pca <- PCA(Stress_Genes_ImmTrt_PCA.impute$completeObs) >Stress_Genes_StrTrt_PCA.impute.pca <- PCA(Stress_Genes_StrTrt_PCA.impute$completeObs) #Get the eigenvalues, variance, and variables. Dim1 and Dim2 were used to create Scatterplot graphs to represent the data for each PCA. >summary(Stress_Genes_Control_PCA.impute.pca) >summary(Stress_Genes_ImmTrt_PCA.impute.pca) >summary(Stress_Genes_StrTrt_PCA.impute.pca) #Retrieve standard error of each population. >Stress_Genes_Control_PCA.impute.pca$call$exart.type >Stress_Genes_ImmTrt_PCA.impute.pca$call$exart.type >Stress_Genes_StrTrt_PCA.impute.pca$call$exart.type