-
Notifications
You must be signed in to change notification settings - Fork 0
/
CWOC_chrysoscelis.Rmd
175 lines (129 loc) · 7.24 KB
/
CWOC_chrysoscelis.Rmd
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
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
---
title: "Counts Without Complexes Analysis for *Hyla chrysoscelis*"
author: "Lauren White"
date: "April 12, 2018"
output: html_document
---
## Load libraries, dataset, and edgelist
```{r}
library(igraph)
#edgelist2<-read.csv("frogedgelist.csv")
edgelist<-read.csv("frogedgelist_cwoc.csv")
#cwc<-read.csv("CountsWithComplexes.csv")
cwoc<-read.csv("CountsWithoutComplexes.csv")
```
## Translate into a weighted affiliation matrix
Using sparse matrices
```{r}
library('Matrix')
A <- spMatrix(nrow=length(unique(edgelist$Species)),
ncol=length(unique(edgelist$Location)),
i = as.numeric(factor(edgelist$Species)),
j = as.numeric(factor(edgelist$Location)),
x = rep(1, length(as.numeric(edgelist$Species))) )
row.names(A) <- levels(factor(edgelist$Species))
colnames(A) <- levels(factor(edgelist$Location))
A
#Frogs have edges if they were at same transect & stop on same day
Arow <- tcrossprod(A)
Arow
```
## Global network for species that co-call with *Hyla chrysoscelis*
A global network with all the frogs that co-call with *Hyla chrysoscelis*. Unweighted edges.
```{r}
#which row in the association matrix corresponds to Hyla chrysoscelis
species<-which(rownames(Arow)=="Hyla chrysoscelis")
cocall<-Arow[,species] #Row in matrix containing Hylea chrysoscelis
#Which species co-call with Hyla cinera (edge weight > 0)
cocall<-which(cocall>0)
#Generate sub containing only species that co-call with Hyla chrysoscelis
mat_all<-Arow[which(rownames(Arow)%in% names(cocall)), which(colnames(Arow) %in% names(cocall))]
mat_all[mat_all > 1] <- 1
#Create igraph object
g<-graph_from_adjacency_matrix(mat_all, mode = "undirected", weighted = NULL, diag = FALSE, add.colnames = NULL, add.rownames = NA)
summary(g)
V(g)$number=as.character(cwoc$Species.ID[match(V(g)$name,cwoc$Species)]) # This code says to create a vertex attribute called "Sex" by extracting the value of the column "Sex" in the attributes file when the Bird ID number matches the vertex name.
#Optional: assign the "IsHyla" attribute as the vertex color (Just comment out if you want everything the same color)
attrib_all<-data.frame(Species=names(cocall), IsHyla=names(cocall)=="Hyla chrysoscelis")
V(g)$IsHyla=as.character(attrib_all$IsHyla[match(V(g)$name, attrib_all$Species)])
V(g)$color=V(g)$IsHyla
V(g)$color=gsub(FALSE,"white",V(g)$color)
V(g)$color=gsub(TRUE,"white",V(g)$color)
#Look at degree distribution for glboal network
hist(degree(g))
library(ggplot2)
tiff("Chrysoscelis_Global_dd_cwoc.tiff", height =10 , width =12, units = "cm", compression = "lzw", res = 1200)
qplot(degree(g), geom="histogram", main = "Degree distribution: Co-calling with Hyla chrysoscelis", xlab = "degree", fill=I("blue"), col=I("black"), alpha=I(.2), xlim=c(0,60))
dev.off()
ps.options(fonts = "serif") #default font for igraph package is serif; not immediately recognized by postscript function
postscript("Chrysoscelis_Global_dd_cwoc.eps", height =10 , width =12)
qplot(degree(g), geom="histogram", main = "Degree distribution: Co-calling with Hyla chrysoscelis", xlab = "degree", fill=I("blue"), col=I("black"), alpha=I(.2), xlim=c(0,60))
dev.off()
#Reduce labeling relative to node size
V(g)$label.cex = .40
V(g)$vertex.size =6
l<-layout.fruchterman.reingold(g)
#Note: #21 corresponds to the Species.ID for *Hyla chrysoscelis*
tiff("Chrysoscelis_Global_network_cwoc_labeled.tiff", height =10 , width =10, units = "cm", compression = "lzw", res = 1200)
plot.igraph(g,vertex.label=V(g)$number, layout=l, vertex.size=12, edge.width=.2, edge.color="black")
dev.off()
postscript("Chrysoscelis_Global_network_cwoc_labeled.eps", height =10 , width =10)
plot.igraph(g,vertex.label=V(g)$number, layout=l, vertex.size=12, edge.width=.2, edge.color="black")
dev.off()
tiff("Chrysoscelis_Global_network_cwoc_unlabeled.tiff", height =10 , width =10, units = "cm", compression = "lzw", res = 1200)
plot.igraph(g,vertex.label=NA, layout=l, vertex.size=12, edge.width=.2, edge.color="black")
dev.off()
postscript("Chrysoscelis_Global_network_cwoc_unlabeled.eps")
plot.igraph(g,vertex.label=NA, layout=l, vertex.size=12, edge.width=.2, edge.color="black")
dev.off()
#Note: you can adjust node location manually using this command:
#tkplot(g)
```
## Top ten co-callers
A network with *Hyla chrysoscelis* and the 10 frogs that most commonly co-call. With weighted edges.
```{r}
species<-which(rownames(Arow)=="Hyla chrysoscelis")
cocall<-Arow[,species] #Row in matrix containing Hylea chrysoscelis
cocall<-sort(cocall, decreasing = TRUE)
topten<-cocall[1:11] #includeing Hyla chrysoscelis (adjust last number in sequence to get top x number of co-calling species)
mat10<-Arow[which(rownames(Arow)%in% names(topten)), which(colnames(Arow) %in% names(topten))] #subnetwork containing only top ten species that co-call with Hyla chrysoscelis
h<-graph_from_adjacency_matrix(mat10, mode = "undirected", weighted = TRUE, diag = FALSE, add.colnames = NULL, add.rownames = NA)
summary(h)
E(h)$weight #weighted edges this time
attrib<-data.frame(Species=names(topten), Weight=topten, IsHyla=names(topten)=="Hyla chrysoscelis") #create an attribute data frame- is the species the focal species, or not?
#assign as vertex attribute- is species focal species or not?
V(h)$IsHyla=as.character(attrib$IsHyla[match(V(h)$name,attrib$Species)])
V(h)$color=V(h)$IsHyla #assign the "IsHyla" attribute as the vertex color
V(h)$color=gsub(FALSE,"white",V(h)$color)
V(h)$color=gsub(TRUE,"white",V(h)$color)
l <- layout_with_fr(h)
l2 <- layout_as_star(h, center = V(h)[4], order = NULL)
l3 <-layout.circle(h)
#reduce vertex label size
V(h)$label.cex = .30
tiff("Chrysoscelis_Topten_star_cwoc_labeled.tiff", height =10 , width =10, units = "cm", compression = "lzw", res = 1200)
plot.igraph(h,vertex.label=V(h)$name,layout=l2,edge.width=sqrt(E(h)$weight)/15, edge.color="black")
dev.off()
ps.options(fonts = "serif") #default font for igraph package is serif; not immediately recognized by postscript function
postscript("Chrysoscelis_Topten_star_cwoc_labeled.eps")
plot.igraph(h,vertex.label=V(h)$name,layout=l2,edge.width=sqrt(E(h)$weight)/15, edge.color="black")
dev.off()
tiff("Chrysoscelis_Topten_star_cwoc_unlabeled.tiff", height =10 , width =10, units = "cm", compression = "lzw", res = 1200)
plot.igraph(h,vertex.label=NA,layout=l2,edge.width=sqrt(E(h)$weight)/15, edge.color="black")
dev.off()
postscript("Chrysoscelis_Topten_star_cwoc_unlabeled.eps", height =10 , width =10)
plot.igraph(h,vertex.label=NA,layout=l2,edge.width=sqrt(E(h)$weight)/15, edge.color="black")
dev.off()
tiff("Chrysoscelis_Topten_circle_cwoc_labeled.tiff", height =10 , width =10, units = "cm", compression = "lzw", res = 1200)
plot.igraph(h,vertex.label=V(h)$name,layout=l3,edge.width=sqrt(E(h)$weight)/15, edge.color="black")
dev.off()
postscript("Chrysoscelis_Topten_circle_cwoc_labeled.eps", height =10 , width =10)
plot.igraph(h,vertex.label=V(h)$name,layout=l3,edge.width=sqrt(E(h)$weight)/15, edge.color="black")
dev.off()
tiff("Chrysoscelis_Topten_circle_cwoc_unlabeled.tiff", height =10 , width =10, units = "cm", compression = "lzw", res = 1200)
plot.igraph(h,vertex.label=NA,layout=l3,edge.width=sqrt(E(h)$weight)/15, edge.color="black")
dev.off()
postscript("Chrysoscelis_Topten_circle_cwoc_unlabeled.eps")
plot.igraph(h,vertex.label=NA,layout=l3,edge.width=sqrt(E(h)$weight)/15, edge.color="black")
dev.off()
```