-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathabundances.R
43 lines (33 loc) · 1.77 KB
/
abundances.R
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
#
# Supp step 2: Just look at some fun facts and figures
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# 'normalized_order_counts.csv' is generated by 'normalize.R'
normalizedData <- read.csv('data/normalized_order_counts.csv', row.names=1)
exoticData <- subset(normalizedData, Plants=="exotic")
nativeData <- subset(normalizedData, Plants=="native")
# Per-treatment total counts of each order
exoticCounts <- apply(exoticData[,c(7:dim(exoticData)[2])], MARGIN=2, FUN=sum)
nativeCounts <- apply(nativeData[,c(7:dim(nativeData)[2])], MARGIN=2, FUN=sum)
counts <- cbind(exoticCounts, nativeCounts)
colnames(counts) <- c("Exotic", "Native")
# Counts sorted in descending order (whichever treatment has a higher count is used in the sorting)
counts <- counts[order(apply(counts, MARGIN=1, FUN=max), decreasing=T),]
# Visualize the 10 highest count orders
bigCounts <- head(counts, n=10)
colrs <- rainbow(dim(bigCounts)[1])
barplot(bigCounts, col=colrs, legend=rownames(bigCounts))
# 'gtools' provides the function 'foldChanges'
library(gtools)
# Visualize the degree of difference in order counts between treatments
noZeroCounts <- counts[(apply(counts, MARGIN=1, FUN=min) > 0),]
foldChanges <- foldchange(noZeroCounts[,"Exotic"], noZeroCounts[,"Native"])
foldChanges <- foldChanges[order(abs(foldChanges))]
fcColr <- c("red", "blue")[(foldChanges > 0) + 1]
barplot(foldChanges, col=fcColr)
legend(0, 6, legend=c("More prevalent in native", "More prevalent in exotic"), col=c("red", "blue"), pch=c(19,19))
# Visualize again, with a minor adjustment so that orders with no difference appear to have no difference (originally, the difference looked substantial)
adj <- c(1, -1)[(foldChanges > 0) + 1]
adjFoldChanges <- foldChanges + adj
barplot(adjFoldChanges, col=fcColr)
q(save="no")