IHC-rxfp2-stats
J Aldersey
27/01/2022
Code/meanings
2 = background
3 = rxfp2+stain
5 = neg nuclei
count = number of individual areas
area = total number of pixels
IMPORT AND ARRANGE DATA
r532hb55 <- read.csv("C://IHC/rxfp2/r532hb55_Regions_rearranged.csv")
#add meta data and reorder
r532hb55$sample <- "532"
r532hb55$tissue <- "HB"
r532hb55$section <- "55"
r532hb55$genotype <- "pp"
r532hb55$ab<- "rxfp2"
r532hb55 <- r532hb55[, c(7, 8, 9, 10, 11, 2, 1, 3, 5, 4, 6)]
r532hb73 <- read.csv("C://IHC/rxfp2/r532hb73_Regions_rearranged.csv")
#add meta data and reorder
r532hb73$sample <- "532"
r532hb73$tissue <- "HB"
r532hb73$section <- "73"
r532hb73$genotype <- "pp"
r532hb73$ab<- "rxfp2"
r532hb73 <- r532hb73[, c(7, 8, 9, 10, 11, 2, 1, 3, 5, 4, 6)]
r546hb53 <- read.csv("C://IHC/rxfp2/r546hb53_Regions_rearranged.csv")
#add meta data and reorder
r546hb53$sample <- "546"
r546hb53$tissue <- "HB"
r546hb53$section <- "53"
r546hb53$genotype <- "pp"
r546hb53$ab<- "rxfp2"
r546hb53 <- r546hb53[, c(7, 8, 9, 10, 11, 2, 1, 3, 5, 4, 6)]
r546hb57 <- read.csv("C://IHC/rxfp2/r546hb57_Regions_rearranged.csv")
#add meta data and reorder
r546hb57$sample <- "546"
r546hb57$tissue <- "HB"
r546hb57$section <- "57"
r546hb57$genotype <- "pp"
r546hb57$ab<- "rxfp2"
r546hb57 <- r546hb57[, c(7, 8, 9, 10, 11, 2, 1, 3, 5, 4, 6)]
r546hb73 <- read.csv("C://IHC/rxfp2/r546hb73_Regions_rearranged.csv")
#add meta data and reorder
r546hb73$sample <- "546"
r546hb73$tissue <- "HB"
r546hb73$section <- "73"
r546hb73$genotype <- "pp"
r546hb73$ab<- "rxfp2"
r546hb73 <- r546hb73[, c(7, 8, 9, 10, 11, 2, 1, 3, 5, 4, 6)]
r581hb35 <- read.csv("C://IHC/rxfp2/r581hb35_Regions_rearranged.csv")
#add meta data and reorder
r581hb35$sample <- "581"
r581hb35$tissue <- "HB"
r581hb35$section <- "35"
r581hb35$genotype <- "pp"
r581hb35$ab<- "rxfp2"
r581hb35 <- r581hb35[, c(7, 8, 9, 10, 11, 2, 1, 3, 5, 4, 6)]
r581hb49 <- read.csv("C://IHC/rxfp2/r581hb49_Regions_rearranged.csv")
#add meta data and reorder
r581hb49$sample <- "581"
r581hb49$tissue <- "HB"
r581hb49$section <- "49"
r581hb49$genotype <- "pp"
r581hb49$ab<- "rxfp2"
r581hb49 <- r581hb49[, c(7, 8, 9, 10, 11, 2, 1, 3, 5, 4, 6)]
r581hb63 <- read.csv("C://IHC/rxfp2/r581hb63_Regions_rearranged.csv")
#add meta data and reorder
r581hb63$sample <- "581"
r581hb63$tissue <- "HB"
r581hb63$section <- "63"
r581hb63$genotype <- "pp"
r581hb63$ab<- "rxfp2"
r581hb63 <- r581hb63[, c(7, 8, 9, 10, 11, 2, 1, 3, 5, 4, 6)]
r618hb43 <- read.csv("C://IHC/rxfp2/r618hb43_Regions_rearranged.csv")
#add meta data and reorder
r618hb43$sample <- "618"
r618hb43$tissue <- "HB"
r618hb43$section <- "43"
r618hb43$genotype <- "pp"
r618hb43$ab<- "rxfp2"
r618hb43 <- r618hb43[, c(7, 8, 9, 10, 11, 2, 1, 3, 5, 4, 6)]
r618hb63 <- read.csv("C://IHC/rxfp2/r618hb63_Regions_rearranged.csv")
#add meta data and reorder
r618hb63$sample <- "618"
r618hb63$tissue <- "HB"
r618hb63$section <- "63"
r618hb63$genotype <- "pp"
r618hb63$ab<- "rxfp2"
r618hb63 <- r618hb63[, c(7, 8, 9, 10, 11, 2, 1, 3, 5, 4, 6)]
r668hb17 <- read.csv("C://IHC/rxfp2/r668hb17_Regions_rearranged.csv")
#add meta data and reorder
r668hb17$sample <- "668"
r668hb17$tissue <- "HB"
r668hb17$section <- "17"
r668hb17$genotype <- "pp"
r668hb17$ab<- "rxfp2"
r668hb17 <- r668hb17[, c(7, 8, 9, 10, 11, 2, 1, 3, 5, 4, 6)]
r668hb51 <- read.csv("C://IHC/rxfp2/r668hb51_Regions_rearranged.csv")
#add meta data and reorder
r668hb51$sample <- "668"
r668hb51$tissue <- "HB"
r668hb51$section <- "51"
r668hb51$genotype <- "pp"
r668hb51$ab<- "rxfp2"
r668hb51 <- r668hb51[, c(7, 8, 9, 10, 11, 2, 1, 3, 5, 4, 6)]
r668hb61 <- read.csv("C://IHC/rxfp2/r668hb61_Regions_rearranged.csv")
#add meta data and reorder
r668hb61$sample <- "668"
r668hb61$tissue <- "HB"
r668hb61$section <- "61"
r668hb61$genotype <- "pp"
r668hb61$ab<- "rxfp2"
r668hb61 <- r668hb61[, c(7, 8, 9, 10, 11, 2, 1, 3, 5, 4, 6)]
r667hb17 <- read.csv("C://IHC/rxfp2/r667hb17_Regions_rearranged.csv")
#add meta data and reorder
r667hb17$sample <- "677"
r667hb17$tissue <- "HB"
r667hb17$section <- "17"
r667hb17$genotype <- "PP"
r667hb17$ab<- "rxfp2"
r667hb17 <- r667hb17[, c(7, 8, 9, 10, 11, 2, 1, 3, 5, 4, 6)]
r667hb27 <- read.csv("C://IHC/rxfp2/r667hb27_Regions_rearranged.csv")
#add meta data and reorder
r667hb27$sample <- "677"
r667hb27$tissue <- "HB"
r667hb27$section <- "27"
r667hb27$genotype <- "PP"
r667hb27$ab<- "rxfp2"
r667hb27 <- r667hb27[, c(7, 8, 9, 10, 11, 2, 1, 3, 5, 4, 6)]
r667hb49 <- read.csv("C://IHC/rxfp2/r667hb49_Regions_rearranged.csv")
#add meta data and reorder
r667hb49$sample <- "677"
r667hb49$tissue <- "HB"
r667hb49$section <- "49"
r667hb49$genotype <- "PP"
r667hb49$ab<- "rxfp2"
r667hb49 <- r667hb49[, c(7, 8, 9, 10, 11, 2, 1, 3, 5, 4, 6)]
r694fs23 <- read.csv("C://IHC/rxfp2/r694fs23_Regions_rearranged.csv")
#add meta data and reorder
r694fs23$sample <- "694"
r694fs23$tissue <- "FS"
r694fs23$section <- "23"
r694fs23$genotype <- "PP"
r694fs23$ab<- "rxfp2"
r694fs23 <- r694fs23[, c(7, 8, 9, 10, 11, 2, 1, 3, 5, 4, 6)]
r709hb39 <- read.csv("C://IHC/rxfp2/r709hb39_Regions_rearranged.csv")
#add meta data and reorder
r709hb39$sample <- "709"
r709hb39$tissue <- "HB"
r709hb39$section <- "39"
r709hb39$genotype <- "PP"
r709hb39$ab<- "rxfp2"
r709hb39 <- r709hb39[, c(7, 8, 9, 10, 11, 2, 1, 3, 5, 4, 6)]
#replace with 709hb7 and 709hb27 as they were randomly selected over 709hb49
r709hb27 <- read.csv("C://IHC/rxfp2/r709hb27_Regions_rearranged.csv")
#add meta data and reorder
r709hb27$sample <- "709"
r709hb27$tissue <- "HB"
r709hb27$section <- "27"
r709hb27$genotype <- "PP"
r709hb27$ab<- "rxfp2"
r709hb27 <- r709hb27[, c(7, 8, 9, 10, 11, 2, 1, 3, 5, 4, 6)]
r709hb7 <- read.csv("C://IHC/rxfp2/r709hb7_Regions_rearranged.csv")
#add meta data and reorder
r709hb7$sample <- "709"
r709hb7$tissue <- "HB"
r709hb7$section <- "7"
r709hb7$genotype <- "PP"
r709hb7$ab<- "rxfp2"
r709hb7 <- r709hb7[, c(7, 8, 9, 10, 11, 2, 1, 3, 5, 4, 6)]
r736hb13 <- read.csv("C://IHC/rxfp2/r736hb13_Regions_rearranged.csv")
#add meta data and reorder
r736hb13$sample <- "736"
r736hb13$tissue <- "HB"
r736hb13$section <- "13"
r736hb13$genotype <- "pp"
r736hb13$ab<- "rxfp2"
r736hb13 <- r736hb13[, c(7, 8, 9, 10, 11, 2, 1, 3, 5, 4, 6)]
r736hb33 <- read.csv("C://IHC/rxfp2/r736hb33_Regions_rearranged.csv")
#add meta data and reorder
r736hb33$sample <- "736"
r736hb33$tissue <- "HB"
r736hb33$section <- "33"
r736hb33$genotype <- "pp"
r736hb33$ab<- "rxfp2"
r736hb33 <- r736hb33[, c(7, 8, 9, 10, 11, 2, 1, 3, 5, 4, 6)]
r736hb63 <- read.csv("C://IHC/rxfp2/r736hb63_Regions_rearranged.csv")
#add meta data and reorder
r736hb63$sample <- "736"
r736hb63$tissue <- "HB"
r736hb63$section <- "63"
r736hb63$genotype <- "pp"
r736hb63$ab<- "rxfp2"
r736hb63 <- r736hb63[, c(7, 8, 9, 10, 11, 2, 1, 3, 5, 4, 6)]
raw_data <- rbind(r532hb55, r532hb73, r546hb53, r546hb57, r546hb73, r581hb35,
r581hb49, r581hb63, r618hb43, r618hb63, r668hb17, r668hb51,
r668hb61, r736hb13, r736hb33, r736hb63, r694fs23, r709hb39,
r709hb7, r709hb27, r667hb17, r667hb27, r667hb49) Checking for Normality of raw data
Conclusion: Not all data is normal
##
## Shapiro-Wilk normality test
##
## data: rxfp2_raw_data$count_2
## W = 0.88001, p-value = 0.01001
##
## Shapiro-Wilk normality test
##
## data: rxfp2_raw_data$count_3
## W = 0.77424, p-value = 0.0001525
##
## Shapiro-Wilk normality test
##
## data: rxfp2_raw_data$count_5
## W = 0.93953, p-value = 0.1756
##
## Shapiro-Wilk normality test
##
## data: rxfp2_raw_data$area_2
## W = 0.94807, p-value = 0.2667
##
## Shapiro-Wilk normality test
##
## data: rxfp2_raw_data$area_3
## W = 0.97321, p-value = 0.7651
##
## Shapiro-Wilk normality test
##
## data: rxfp2_raw_data$area_5
## W = 0.8343, p-value = 0.001432
####averaging####
1 = rxfp2+stain, 2 = background, 5 = neg nuclei
#bind data from the same sample
r532HB <-rbind(r532hb55,
r532hb73)
#calculate the average and standard deviation of each measure
mean_532HB <- data.frame (sample = "532HB",
genotype = "pp",
ab = "rxfp2",
count_3 = mean(r532HB$count_3),
count_2 = mean(r532HB$count_2),
count_5 = mean(r532HB$count_5),
area_3 = mean(r532HB$area_3),
area_2 = mean(r532HB$area_2),
area_5 = mean(r532HB$area_5),
sd_count_3 = sd(r532HB$count_3),
sd_count_2 = sd(r532HB$count_2),
sd_count_5 = sd(r532HB$count_5),
sd_area_3 = sd(r532HB$area_3),
sd_area_2 = sd(r532HB$area_2),
sd_area_5 = sd(r532HB$area_5))
r546HB <-rbind(r546hb53,
r546hb57,
r546hb73)
mean_546HB <- data.frame (sample = "546HB",
genotype = "pp",
ab = "rxfp2",
count_3 = mean(r546HB$count_3),
count_2 = mean(r546HB$count_2),
count_5 = mean(r546HB$count_5),
area_3 = mean(r546HB$area_3),
area_2 = mean(r546HB$area_2),
area_5 = mean(r546HB$area_5),
sd_count_3 = sd(r546HB$count_3),
sd_count_2 = sd(r546HB$count_2),
sd_count_5 = sd(r546HB$count_5),
sd_area_3 = sd(r546HB$area_3),
sd_area_2 = sd(r546HB$area_2),
sd_area_5 = sd(r546HB$area_5))
r581HB <-rbind(r581hb35,
r581hb49,
r581hb63)
mean_581HB <- data.frame (sample = "581HB",
genotype = "pp",
ab = "rxfp2",
count_3 = mean(r581HB$count_3),
count_2 = mean(r581HB$count_2),
count_5 = mean(r581HB$count_5),
area_3 = mean(r581HB$area_3),
area_2 = mean(r581HB$area_2),
area_5 = mean(r581HB$area_5),
sd_count_3 = sd(r581HB$count_3),
sd_count_2 = sd(r581HB$count_2),
sd_count_5 = sd(r581HB$count_5),
sd_area_3 = sd(r581HB$area_3),
sd_area_2 = sd(r581HB$area_2),
sd_area_5 = sd(r581HB$area_5))
r618HB <-rbind(r618hb43, r618hb63, r668hb17)
mean_618HB <- data.frame (sample = "618HB",
genotype = "pp",
ab = "rxfp2",
count_3 = mean(r618HB$count_3),
count_2 = mean(r618HB$count_2),
count_5 = mean(r618HB$count_5),
area_3 = mean(r618HB$area_3),
area_2 = mean(r618HB$area_2),
area_5 = mean(r618HB$area_5),
sd_count_3 = sd(r618HB$count_3),
sd_count_2 = sd(r618HB$count_2),
sd_count_5 = sd(r618HB$count_5),
sd_area_3 = sd(r618HB$area_3),
sd_area_2 = sd(r618HB$area_2),
sd_area_5 = sd(r618HB$area_5))
r668HB <-rbind(r668hb17, r668hb51,
r668hb61)
mean_668HB <- data.frame (sample = "668HB",
genotype = "pp",
ab = "rxfp2",
count_3 = mean(r668HB$count_3),
count_2 = mean(r668HB$count_2),
count_5 = mean(r668HB$count_5),
area_3 = mean(r668HB$area_3),
area_2 = mean(r668HB$area_2),
area_5 = mean(r668HB$area_5),
sd_count_3 = sd(r668HB$count_3),
sd_count_2 = sd(r668HB$count_2),
sd_count_5 = sd(r668HB$count_5),
sd_area_3 = sd(r668HB$area_3),
sd_area_2 = sd(r668HB$area_2),
sd_area_5 = sd(r668HB$area_5))
r736HB <-rbind(r736hb13, r736hb33, r736hb63)
mean_736HB <- data.frame (sample = "736HB",
genotype = "pp",
ab = "rxfp2",
count_3 = mean(r736HB$count_3),
count_2 = mean(r736HB$count_2),
count_5 = mean(r736HB$count_5),
area_3 = mean(r736HB$area_3),
area_2 = mean(r736HB$area_2),
area_5 = mean(r736HB$area_5),
sd_count_3 = sd(r736HB$count_3),
sd_count_2 = sd(r736HB$count_2),
sd_count_5 = sd(r736HB$count_5),
sd_area_3 = sd(r736HB$area_3),
sd_area_2 = sd(r736HB$area_2),
sd_area_5 = sd(r736HB$area_5))
#only one sample for 694fs
mean_694FS <- data.frame(sample = "694FS",
genotype = "PP",
ab = "rxfp2",
count_3 = r694fs23$count_3,
count_2 = r694fs23$count_2,
count_5 = r694fs23$count_5,
area_3 = r694fs23$area_3,
area_2 = r694fs23$area_2,
area_5 = r694fs23$area_5,
sd_count_3 = 0,
sd_count_2 = 0,
sd_count_5 = 0,
sd_area_3 = 0,
sd_area_2 = 0,
sd_area_5 = 0)
r667HB <-rbind(r667hb17, r667hb27, r667hb49)
mean_667HB <- data.frame (sample = "667HB",
genotype = "PP",
ab = "rxfp2",
count_3 = mean(r667HB$count_3),
count_2 = mean(r667HB$count_2),
count_5 = mean(r667HB$count_5),
area_3 = mean(r667HB$area_3),
area_2 = mean(r667HB$area_2),
area_5 = mean(r667HB$area_5),
sd_count_3 = sd(r667HB$count_3),
sd_count_2 = sd(r667HB$count_2),
sd_count_5 = sd(r667HB$count_5),
sd_area_3 = sd(r667HB$area_3),
sd_area_2 = sd(r667HB$area_2),
sd_area_5 = sd(r667HB$area_5))
r709HB <-rbind(r709hb39,
r709hb7, r709hb27)
mean_709HB <- data.frame (sample = "709HB",
genotype = "PP",
ab = "rxfp2",
count_3 = mean(r709HB$count_3),
count_2 = mean(r709HB$count_2),
count_5 = mean(r709HB$count_5),
area_3 = mean(r709HB$area_3),
area_2 = mean(r709HB$area_2),
area_5 = mean(r709HB$area_5),
sd_count_3 = sd(r709HB$count_3),
sd_count_2 = sd(r709HB$count_2),
sd_count_5 = sd(r709HB$count_5),
sd_area_3 = sd(r709HB$area_3),
sd_area_2 = sd(r709HB$area_2),
sd_area_5 = sd(r709HB$area_5))
rxfp2_meanforsample <-rbind(mean_532HB, mean_546HB, mean_581HB, mean_736HB, mean_668HB, mean_667HB, mean_694FS, mean_709HB)Checking for normality of data averaged per sample.
Conclusion: Data is normal as p-value is >0.05
##
## Shapiro-Wilk normality test
##
## data: rxfp2_meanforsample$count_2
## W = 0.93574, p-value = 0.5697
##
## Shapiro-Wilk normality test
##
## data: rxfp2_meanforsample$count_3
## W = 0.67885, p-value = 0.001296
##
## Shapiro-Wilk normality test
##
## data: rxfp2_meanforsample$count_5
## W = 0.95635, p-value = 0.7747
##
## Shapiro-Wilk normality test
##
## data: rxfp2_meanforsample$area_2
## W = 0.98961, p-value = 0.9945
##
## Shapiro-Wilk normality test
##
## data: rxfp2_meanforsample$area_3
## W = 0.86883, p-value = 0.1468
##
## Shapiro-Wilk normality test
##
## data: rxfp2_meanforsample$area_5
## W = 0.97314, p-value = 0.9214Calculating proportions
ggboxplot(rxfp2_meanforsample, x = "genotype", y = "TA", fill = "genotype") +
theme_classic()
ggboxplot(rxfp2_meanforsample, x = "genotype", y = "pTA_2", fill = "genotype") +
theme_classic() +
labs(y = "Background/total area (proportion)")
rxfp2_meanforsample$pTA_3_percent <- rxfp2_meanforsample$pTA_3*100
ggboxplot(rxfp2_meanforsample, x = "genotype", y = "pTA_3_percent", fill = "genotype") +
theme_classic(base_size = 18) +
labs(y = "% area RXFP2 positive nerves")
ggboxplot(rxfp2_meanforsample, x = "genotype", y = "pTA_5", fill = "genotype") +
theme_classic() +
labs(y = "unstained nuclei area/total area (proportion)")
rxfp2hsum <- subset(rxfp2_meanforsample, genotype=="horned")
rxfp2psum <- subset(rxfp2_meanforsample, genotype=="polled")
write.csv(rxfp2hsum, "rxfp2hsum.csv", row.names = FALSE)
write.csv(rxfp2psum, "rxfp2psum.csv", row.names = FALSE)Fisher’s Exact Test
Testing to see if there is a significant difference in positive areas (nerves) between horned and polled tissues
horn <-rbind(mean_532HB, mean_546HB, mean_581HB, mean_668HB, mean_736HB)
polled <- rbind(mean_694FS, mean_709HB, mean_667HB)horn$A2 <- horn$area_2 + horn$area_5 #calculate negative pixels
horn_means <- as.data.frame(round(colMeans(horn[, c(7, 16)]))) # find mean
horn_means## round(colMeans(horn[, c(7, 16)]))
## area_3 93136
## A2 8656556polled$A2 <- polled$area_2 + polled$area_5
polled_means <- as.data.frame(round(colMeans(polled[, c(7, 16)])))
polled_means## round(colMeans(polled[, c(7, 16)]))
## area_3 39236
## A2 4552424area_3_mean <-mean(polled$area_3)#create new data frame with following format for each variable
# | + | - | Total
# pp | | |
# PP | | |
#Total| | |
dat <- data.frame(
"Positive" = c(horn_means[1,], polled_means[1,]),
"Negative" = c(horn_means[2,], polled_means[2,]),
row.names = c("Horned", "Polled"),
stringsAsFactors = FALSE)
dat## Positive Negative
## Horned 93136 8656556
## Polled 39236 4552424mosaicplot(dat,
main = "Mosaic plot",
color = TRUE)
test <- fisher.test(dat)
test##
## Fisher's Exact Test for Count Data
##
## data: dat
## p-value < 2.2e-16
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 1.233573 1.263244
## sample estimates:
## odds ratio
## 1.248286test$p.value## [1] 6.578718e-303write.csv(dat, "C:/IHC/rxfp2_nerve.csv")Based on the Fisher Exact Test, the number of positive areas (pixels) in horned samples is significantly different to polled samples.