-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathR script
71 lines (58 loc) · 1.88 KB
/
R script
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
# Load required libraries
library(stats)
library(ggplot2)
library(car)
library(nlme)
# Read data
data <- read.csv("edar_xedar_expression.csv")
# 1. Data Preparation and Quality Control ----
# Calculate CV for actin
cv_actin <- function(data) {
cv <- (sd(data$actin)/mean(data$actin)) * 100
return(cv)
}
# 2. Statistical Analysis ----
# Shapiro-Wilk test for normality
shapiro.test(data$edar)
shapiro.test(data$xedar)
# Two-way ANOVA for Edar
edar_anova <- aov(edar ~ species * day, data = data)
summary(edar_anova)
# Two-way ANOVA for Xedar
xedar_anova <- aov(xedar ~ species * day, data = data)
summary(xedar_anova)
# Tukey post-hoc tests
TukeyHSD(edar_anova)
TukeyHSD(xedar_anova)
# 3. Correlation Analysis ----
# Separate by species
mouse_data <- subset(data, species == "mouse")
rat_data <- subset(data, species == "rat")
# Pearson correlation
mouse_cor <- cor.test(mouse_data$edar, mouse_data$xedar, method = "pearson")
rat_cor <- cor.test(rat_data$edar, rat_data$xedar, method = "pearson")
# 4. Polynomial Regression ----
# Edar models
mouse_edar_model <- lm(edar ~ poly(day, 2), data = mouse_data)
rat_edar_model <- lm(edar ~ poly(day, 2), data = rat_data)
# Xedar models
mouse_xedar_model <- lm(xedar ~ poly(day, 3), data = mouse_data)
rat_xedar_model <- lm(xedar ~ poly(day, 2), data = rat_data)
# Calculate R-squared values
summary(mouse_edar_model)$r.squared
summary(rat_edar_model)$r.squared
summary(mouse_xedar_model)$r.squared
summary(rat_xedar_model)$r.squared
# 5. Visualization ----
# Expression patterns plot
ggplot(data, aes(x = day, y = value, color = species)) +
geom_point() +
geom_line() +
facet_wrap(~gene, scales = "free_y") +
theme_minimal() +
labs(x = "Day", y = "Expression Level")
# Save results
save(edar_anova, xedar_anova, mouse_cor, rat_cor,
mouse_edar_model, rat_edar_model,
mouse_xedar_model, rat_xedar_model,
file = "analysis_results.RData")