Click here for other works of the author on RPubs

Load package

library(knitr)

1. Calculate the copepod density for each species for each cruise-station

Import data

species <- read.table("copepod_composition.txt",header = T) #import copepod data
dens <- as.vector(read.table("cop_density.txt",header = T)[[1]]) # import density data

Data manipulation

#convert species frequency into percentage
species = species / 100

#calculate copepod density for each species for each cruise station
species.density = t(apply(species, 1, function(x) x * dens))

kable(t(head(species.density)), digits = 3, col.names = 1:6)
1 2 3 4 5 6
p1 0.000 0.000 0.000 0.000 0.000 0.000
p3 0.000 0.000 0.000 0.000 0.000 0.000
p4 0.000 0.000 0.000 0.000 0.000 0.000
p6 0.000 0.000 0.000 0.000 0.000 0.000
p13 0.000 0.000 0.000 0.000 0.000 0.000
p16 0.000 0.000 4.671 0.000 0.000 0.000
p19 0.000 0.000 27.121 0.000 0.000 0.000
p21 0.000 0.000 0.000 0.000 0.000 0.000
p23 0.000 0.000 33.910 0.000 0.000 0.000
p25 0.000 0.000 15.552 0.000 0.000 0.000
s18 0.000 0.000 0.000 0.000 119.902 0.000
s19 0.000 0.000 0.000 0.000 29.640 0.000
s20 0.000 0.000 0.000 0.000 16.286 0.000
s22 0.000 0.000 0.000 0.000 195.277 0.000
s23 0.000 0.000 0.000 0.000 417.623 0.000
s25 18.996 0.000 0.000 0.000 94.347 0.000
s27 0.000 0.000 31.935 0.000 0.000 15.968
s29 0.000 0.000 14.469 0.000 0.000 72.827
sA 0.000 0.000 48.593 4.129 4.129 0.000
sB 0.000 0.000 39.083 0.000 0.000 0.000
sC 0.000 0.000 47.988 0.000 0.000 0.000
sD 0.000 0.000 8.395 0.000 0.000 0.000
sE 0.000 0.000 160.366 0.000 17.574 0.000
sF 0.000 0.000 24.546 0.000 121.145 0.000
sG 0.000 0.000 17.780 0.000 0.000 0.000
w22 0.000 0.000 0.000 0.000 69.543 0.000
w23 0.000 0.000 0.000 0.000 135.497 0.000
w25 0.000 0.000 0.000 0.000 2.267 0.000
w27 0.000 0.000 0.000 0.000 1.773 0.000
w29 0.000 0.175 0.000 0.000 2.621 0.000
wA 0.000 0.000 0.000 0.000 0.051 0.000
wB 0.000 0.000 1.103 0.000 0.139 0.000
wC 0.000 0.000 0.291 0.000 0.000 0.000
wD 0.000 0.000 0.434 0.000 0.289 0.000

2. For each cruise-station, calculate the species richness (number of species) and Shannon diversity index

#calculate number of species in each cruise station
species_number = apply(species, 2, function(x) length(x[x > 0]))
kable(species_number, col.names = "Number of species")
Number of species
p1 6
p3 12
p4 8
p6 9
p13 31
p16 29
p19 43
p21 7
p23 46
p25 39
s18 39
s19 41
s20 32
s22 25
s23 32
s25 41
s27 40
s29 23
sA 47
sB 49
sC 44
sD 46
sE 44
sF 38
sG 25
w22 18
w23 16
w25 24
w27 33
w29 27
wA 44
wB 54
wC 64
wD 48 
#calculate Shannon diversity index for each cruise station
Shannon_index = apply(species, 2, function(x) - sum(x * log(x), na.rm = T))
kable(Shannon_index, col.names = "Shannon diversity index", digits = 2)
Shannon diversity index
p1 1.08
p3 1.26
p4 1.05
p6 1.11
p13 2.14
p16 1.41
p19 2.58
p21 1.57
p23 2.99
p25 2.48
s18 2.84
s19 2.98
s20 2.57
s22 2.57
s23 2.88
s25 3.00
s27 2.80
s29 2.12
sA 3.11
sB 2.98
sC 2.82
sD 2.94
sE 3.02
sF 2.89
sG 1.69
w22 1.98
w23 1.62
w25 1.84
w27 2.58
w29 2.36
wA 2.61
wB 3.00
wC 3.21
wD 3.01

3. Find dominant species (species >=2% of total composition in any cruise-station) and calculate the average density for the spring, summer, and winter cruise for each dominant species.

Find dominant species in spring

#find dominant species in each station during spring
p.dom = apply(species[, 1:10] >= 0.02, 2, function(x) which(x == TRUE))
p.dom.species = sort(unique(unlist(p.dom)))

#calculate average density for each spring dominant species during each season
p.dom.spring = apply(species.density[p.dom.species, 1:10 ], 1, mean)
p.dom.summer = apply(species.density[p.dom.species, 11:25], 1, mean)
p.dom.winter = apply(species.density[p.dom.species, 26:34], 1, mean)
p.dom.result = cbind(p.dom.spring, p.dom.summer, p.dom.winter)
rownames(p.dom.result) = p.dom.species

#show density of spring dominant species in each season
kable(p.dom.result, col.names = c("Spring", "Summer", "Winter"), caption = "Density of spring dominant species", digits = 3)
Density of spring dominant species
Spring Summer Winter
3 8.125 26.210 0.203
14 86.466 2.500 4.603
16 4.501 1.904 0.137
35 11.927 63.691 0.672
40 11.446 23.416 0.681
54 4.873 40.485 0.229
60 28.353 41.598 23.154
72 10.711 156.407 0.996
76 5.550 3.362 0.000
84 54.093 18.917 9.400
85 626.171 414.205 18.136
88 11.223 419.718 0.038
111 9.938 4.661 0.005
123 28.781 83.105 0.991
126 115.089 38.427 0.000
135 7.362 197.810 0.000
142 161.229 6.446 0.000
158 3.378 16.906 0.011
161 21.429 0.000 0.000
169 33.968 341.827 1.254

Find dominant species in summer

#find dominant species in each station during summer
s.dom=apply(species[,11:25]>=0.02,2,function(x) which(x==TRUE))
s.dom.species=sort(unique(unlist(s.dom)))

#calculate average density for each summer dominant species during each season
s.dom.spring=apply(species.density[s.dom.species,1:10],1,mean)
s.dom.summer=apply(species.density[s.dom.species,11:25],1,mean)
s.dom.winter=apply(species.density[s.dom.species,26:34],1,mean)
s.dom.result=cbind(s.dom.spring, s.dom.summer, s.dom.winter)
rownames(s.dom.result) = s.dom.species

#show density of summer dominant species in each season
kable(s.dom.result, col.names = c("Spring", "Summer", "Winter"), caption = "Density of summer dominant species", digits = 3)
Density of summer dominant species
Spring Summer Winter
3 8.125 26.210 0.203
5 0.000 67.728 23.576
15 3.426 204.692 1.267
20 1.274 42.367 0.170
35 11.927 63.691 0.672
40 11.446 23.416 0.681
51 0.000 28.013 0.864
54 4.873 40.485 0.229
60 28.353 41.598 23.154
72 10.711 156.407 0.996
73 0.307 31.915 0.054
79 1.459 26.637 0.000
80 3.744 9.886 0.000
81 1.760 13.302 0.005
84 54.093 18.917 9.400
85 626.171 414.205 18.136
86 0.000 107.269 0.025
88 11.223 419.718 0.038
106 0.768 20.458 4.463
112 3.691 267.765 0.970
117 0.473 134.511 0.000
118 0.237 53.422 0.000
120 1.248 16.453 1.225
123 28.781 83.105 0.991
126 115.089 38.427 0.000
135 7.362 197.810 0.000
145 1.753 97.267 0.086
147 0.000 75.731 1.564
148 0.973 21.038 0.023
151 0.000 39.107 0.186
158 3.378 16.906 0.011
164 1.123 116.336 0.950
165 1.915 14.383 0.154
169 33.968 341.827 1.254

Find dominant species in Winter

#find dominant species in each station during winter
w.dom=apply(species[,26:34]>=0.02,2,function(x) which(x==TRUE))
w.dom.species=sort(unique(unlist(w.dom)))

#calculate average density for each winter dominant species during each season
w.dom.spring=apply(species.density[w.dom.species,1:10],1,mean)
w.dom.summer=apply(species.density[w.dom.species,11:25],1,mean)
w.dom.winter=apply(species.density[w.dom.species,26:34],1,mean)
w.dom.result=cbind(w.dom.spring, w.dom.summer, w.dom.winter)
rownames(w.dom.result) = w.dom.species

#show density of winter dominant species in each season
kable(w.dom.result, col.names = c("Spring", "Summer", "Winter"), caption = "Density of winter dominant species", digits = 3)
Density of winter dominant species
Spring Summer Winter
3 8.125 26.210 0.203
5 0.000 67.728 23.576
14 86.466 2.500 4.603
15 3.426 204.692 1.267
16 4.501 1.904 0.137
35 11.927 63.691 0.672
40 11.446 23.416 0.681
51 0.000 28.013 0.864
52 0.000 5.324 1.412
54 4.873 40.485 0.229
55 0.301 0.727 0.421
60 28.353 41.598 23.154
72 10.711 156.407 0.996
84 54.093 18.917 9.400
85 626.171 414.205 18.136
106 0.768 20.458 4.463
112 3.691 267.765 0.970
120 1.248 16.453 1.225
123 28.781 83.105 0.991
147 0.000 75.731 1.564
166 3.540 4.435 0.263
169 33.968 341.827 1.254