R语言--果蝇攀爬--ImageJ测量数据处理

发表于 分类于 阅读次数： 本文字数： 757 阅读时长 ≈ 3 分钟

目的

处理ImageJ测量的果蝇攀爬数据

实例

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
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
library(cluster) 
# 设定工作目录
setwd("E:/Desktop/")

# 导入数据
X20180609_WF_IF <- read.csv("E:/Desktop/20180609-WF-IF.csv")
X20180609_CM_AM <- read.csv("E:/Desktop/20180609-CM-AM.csv")
X20180609_CF_AF <- read.csv("E:/Desktop/20180609-CF-AF.csv")

clusGap(X20180609_WF_IF$X[which(X20180609_WF_IF$Slice==1)],FUN=kmeans,nstart=25,K.max=10,B=500)
## 自定义函数，获取每管果蝇的攀爬指数----
ppindex <- function(X){
g1 <- c()
g2 <- c()
g3 <- c()
g4 <- c()
g5 <- c()
g6 <- c()
for (i in 1:max(X$Slice)) {
  slice_x <- X$X[which(X$Slice==i)]
  slice_y <- X$Y[which(X$Slice==i)]
  slice_max <- (max(slice_x)+100)/6
  par(mfrow = c(1,6))
  # 管1
  slice_x <- slice_x[which(slice_x < slice_max)]
  slice_y <- slice_y[which(slice_x < slice_max)]
  plot(slice_x,slice_y)
  slice_y_fq0 <- length(slice_y[which(slice_y < 200)])
  slice_y_fq1 <- length(slice_y[which(slice_y > 200&slice_y < 400)])
  slice_y_fq2 <- length(slice_y[which(slice_y > 400)])-2
  index_1 <- (slice_y_fq0*0+slice_y_fq1*1+slice_y_fq2*2)/(slice_y_fq0+slice_y_fq1+slice_y_fq2)
  g1 <- c(g1,index_1)
  # 管2
  slice_x2 <- slice_x[which(slice_max < slice_x&slice_x < slice_max*2)]
  slice_y2 <- slice_y[which(slice_max < slice_x&slice_x < slice_max*2)]
  plot(slice_x2,slice_y2)
  slice_y2_fq0 <- length(slice_y2[which(slice_y2 < 200)])
  slice_y2_fq1 <- length(slice_y2[which(slice_y2 > 200&slice_y2 < 400)])
  slice_y2_fq2 <- length(slice_y2[which(slice_y2 > 400)])-2
  index_2 <- (slice_y2_fq0*0+slice_y2_fq1*1+slice_y2_fq2*2)/(slice_y2_fq0+slice_y2_fq1+slice_y2_fq2)
  g2 <- c(g2,index_2)
  # 管3
  slice_x3 <- slice_x[which(slice_max*2 < slice_x&slice_x < slice_max*3)]
  slice_y3 <- slice_y[which(slice_max*2 < slice_x&slice_x < slice_max*3)]
  plot(slice_x3,slice_y3)
  slice_y3_fq0 <- length(slice_y3[which(slice_y3 < 200)])
  slice_y3_fq1 <- length(slice_y3[which(slice_y3 > 200&slice_y3 < 400)])
  slice_y3_fq2 <- length(slice_y3[which(slice_y3 > 400)])-2
  index_3 <- (slice_y3_fq0*0+slice_y3_fq1*1+slice_y3_fq2*2)/(slice_y3_fq0+slice_y3_fq1+slice_y3_fq2)
  g3 <- c(g3,index_3)
  # 管4
  slice_x4 <- slice_x[which(slice_max*3 < slice_x&slice_x < slice_max*4)]
  slice_y4 <- slice_y[which(slice_max*3 < slice_x&slice_x < slice_max*4)]
  plot(slice_x4,slice_y4)
  slice_y4_fq0 <- length(slice_y4[which(slice_y4 < 200)])
  slice_y4_fq1 <- length(slice_y4[which(slice_y4 > 200&slice_y4 < 400)])
  slice_y4_fq2 <- length(slice_y4[which(slice_y4 > 400)])-2
  index_4 <- (slice_y4_fq0*0+slice_y4_fq1*1+slice_y4_fq2*2)/(slice_y4_fq0+slice_y4_fq1+slice_y4_fq2)
  g4 <- c(g4,index_4)
  # 管5
  slice_x5 <- slice_x[which(slice_max*4 < slice_x&slice_x < slice_max*5)]
  slice_y5 <- slice_y[which(slice_max*4 < slice_x&slice_x < slice_max*5)]
  plot(slice_x5,slice_y5)
  slice_y5_fq0 <- length(slice_y5[which(slice_y5 < 200)])
  slice_y5_fq1 <- length(slice_y5[which(slice_y5 > 200&slice_y5 < 400)])
  slice_y5_fq2 <- length(slice_y5[which(slice_y5 > 400)])-2
  index_5 <- (slice_y5_fq0*0+slice_y5_fq1*1+slice_y5_fq2*2)/(slice_y5_fq0+slice_y5_fq1+slice_y5_fq2)
  g5 <- c(g5,index_5)
  # 管6
  slice_x6 <- slice_x[which(slice_max*5 < slice_x&slice_x < slice_max*6)]
  slice_y6 <- slice_y[which(slice_max*5 < slice_x&slice_x < slice_max*6)]
  plot(slice_x6,slice_y6)
  slice_y6_fq0 <- length(slice_y6[which(slice_y6 < 200)])
  slice_y6_fq1 <- length(slice_y6[which(slice_y6 > 200&slice_y6 < 400)])
  slice_y6_fq2 <- length(slice_y6[which(slice_y6 > 400)])-2
  index_6 <- (slice_y6_fq0*0+slice_y6_fq1*1+slice_y6_fq2*2)/(slice_y6_fq0+slice_y6_fq1+slice_y6_fq2)
  g6 <- c(g6,index_6)
  }
  flyindex <- data.frame(g1,g2,g3,g4,g5,g6)
  return(flyindex)
}

#汇总----
WI_F_Index <- ppindex(X20180609_WF_IF) 
names(WI_F_Index) <- c("WF1","WF2","WF3","IF1","IF2","IF3")
WI_F_Index <- ppindex(X20180609_CF_AF)
names(CA_F_Index) <- c("CF1","CF2","CF3","AF1","AF2","AF3")
CA_M_Index <- ppindex(X20180609_CM_AM)
names(CA_M_Index) <- c("CM1","CM2","CM3","AM1","AM2","AM3")

# 保存
library(xlsx)
write.xlsx2(WI_F_Index, sheetName = "WI_F_Index", file='flyIndex.xlsx',append = TRUE)
write.xlsx2(CA_F_Index, sheetName = "CA_F_Index", file='flyIndex.xlsx',append = TRUE)
write.xlsx2(CA_M_Index, sheetName = "CA_M_Index", file='flyIndex.xlsx',append = TRUE)
dev.off()

###########################################################################################
#单管----
# 导入数据
X20180609_IM3 <- read.csv("E:/Desktop/20180609-IM3.csv")
X20180609_IM2 <- read.csv("E:/Desktop/20180609-IM2.csv")
X20180609_IM1 <- read.csv("E:/Desktop/20180609-IM1.csv")
X20180609_WM3 <- read.csv("E:/Desktop/20180609-WM3.csv")
X20180609_WM2 <- read.csv("E:/Desktop/20180609-WM2.csv")
X20180609_WM1 <- read.csv("E:/Desktop/20180609-WM1.csv")

pp2index <- function(x) {
  index <- c()
  for (i in 1:max(x$Slice)) {
    slice_x <- x$X[which(x$Slice == i)]
    slice_y <- x$Y[which(x$Slice == i)]
    slice_y_fq0 <- length(slice_y[which(slice_y < 340)])
    slice_y_fq1 <- length(slice_y[which(slice_y > 340 &
                                          slice_y < 600)])
    if ((length(slice_y[which(slice_y > 600 &
                              slice_y < 900)]) - 1) < 0) {
      slice_y_fq2 <- 0
    } else {
      slice_y_fq2 <- length(slice_y[which(slice_y > 600 &
                                            slice_y < 900)]) - 1
    }
    if ((length(slice_y[which(slice_y > 900 &
                              slice_y < 1200)]) - 1) < 0) {
      slice_y_fq3 <- 0
    } else {
      slice_y_fq3 <- length(slice_y[which(slice_y > 900 &
                                            slice_y < 1200)]) - 1
    }
    if ((length(slice_y[which(slice_y > 1200)]) - 1) < 0) {
      slice_y_fq4 <- 0
    } else {
      slice_y_fq4 <- length(slice_y[which(slice_y > 1200)]) - 1
    }
    
    index_1 <- (slice_y_fq0 * 0 + slice_y_fq1 * 1 + slice_y_fq2 * 2 + slice_y_fq3 *
                  3 + slice_y_fq4 * 4) / (slice_y_fq0 + slice_y_fq1 + slice_y_fq2 + slice_y_fq3 +
                                            slice_y_fq4)
    index <- c(index,index_1)
  }
  return(index)
}
IM3 <- pp2index(X20180609_IM3)
IM2 <- pp2index(X20180609_IM2)
IM1 <- pp2index(X20180609_IM1)
WM1 <- pp2index(X20180609_WM1)
WM2 <- pp2index(X20180609_WM3)
WM3 <- pp2index(X20180609_WM2)
WI_M_Index <- cbind(WM1=WM1,WM2=WM2,WM3=WM3,IM1=IM1,IM2=IM2,IM3=IM3)
write.xlsx2(WI_M_Index, sheetName = "WI_M_Index", file='flyIndex.xlsx',append = TRUE)