-
Notifications
You must be signed in to change notification settings - Fork 0
/
Heatmap+Mantel_FST_SOX_BS.R
135 lines (121 loc) · 3.83 KB
/
Heatmap+Mantel_FST_SOX_BS.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
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
# Plot heatmap of Bathymodiolus SOX symbiont FST values and test correlation with sampling year & host species 2020-01-15
# Merle Ücker
# Load libraries
library(gplots)
library(RColorBrewer)
library(vegan)
# Load data
fst <- read.table("Pairwise_mean_FST_SOX_symbionts_Broken_Spur.csv", sep = "\t", header = TRUE, check.names = FALSE)
# Set row names to bin names
rownames <- fst[,1]
row.names(fst) <- rownames
fst[,1] <- NULL
fst_matrix <- as.matrix(fst)
# Load metadata
metadata <- read.table("metadata2", sep="\t", header = TRUE)
# metadata2:
# Species Year
# Bhyb 2001
# Bhyb 2001
# Bput 2001
# Bhyb 1997
# Bput 2001
# Bhyb 1997
# Bput 1997
# Bput 2001
# Bhyb 1997
# Bhyb 1997
# Bhyb 1997
# Bput 1997
# Bhyb 1997
# Bhyb 2001
# Bput 2001
# Bput 2001
# Bput 2001
# Bhyb 2001
# Bhyb 2001
# Bput 2001
# Bhyb 2001
# Bhyb 2001
# Bput 2001
# Bput 2001
# Bhyb 2001
# Bhyb 2001
# Bhyb 2001
# Bhyb 2001
# Bput 2001
# Bput 2001
year = as.matrix(vegdist(metadata$Year, method = "euclidean"))
species <- read.table("Pairwise_distances_host_SNPs_FST.txt")
# Run Mantel test to test for correlation
fst_year = mantel(fst_matrix, year, method = "spearman", na.rm = TRUE, permutations = 9999)
fst_year
fst_species = mantel(fst_matrix, species, method = "spearman", na.rm = TRUE, permutations = 9999)
fst_species
# Plotting
my_palette <- colorRampPalette(c("white", "gray88", "#00bcd4"))(n = 30)
# Save as pdf image
pdf("Heatmap_FST_BS_nodensity.pdf",
width = 10,
height = 10,
pointsize = 12)
# Make cluster for dendogram
distance = dist(fst_matrix, method = "euclidean")
cluster = hclust(distance, method = "complete")
# Plot heatmap
heatmap.2(fst_matrix,
main = "", # heat map title
keysize=1.2,
key.title = "" ,
key.xlab = expression('F'[ST]),
key.ylab = "Density",
density.info="none", # turns off density plot inside color legend
trace="none", # turns off trace lines inside the heat map
margins =c(10,10), # widens margins around plot
col=my_palette, # use on color palette defined earlier
#breaks=col_breaks, # enable color transition at specified limits
dendrogram="both", # only draw a row dendrogram
Rowv = as.dendrogram(cluster),
Colv= rev(as.dendrogram(cluster)), # turn off column clustering
cexRow = 1.2,
cexCol = 1.2,
RowSideColors = c( # can be used to have coloured bars according to some feature, e.g. location or species, on the side
rep("#ffb219", 2),
rep("#cc0000", 2),
rep("#ffb219", 1),
rep("#cc0000", 1),
rep("#ffb219", 4),
rep("#cc0000", 1),
rep("#ffb219", 2),
rep("#cc0000", 3),
rep("#ffb219", 2),
rep("#cc0000", 1),
rep("#ffb219", 3),
rep("#cc0000", 2),
rep("#ffb219", 4),
rep("#cc0000", 2)),
ColSideColors = c(
rep("#788497", 1),
rep("#afb6c2", 1),
rep("#788497", 1),
rep("#afb6c2", 9),
rep("#788497", 18)))
# Adding a colour legend for the categories
# plot.new () # enable this to just plot the legend without the plot
par(lend = 2) # makes colour blocks in legend edgy
# Create a legend for the colours used in RowSideColors. Syntax: legend = c("species/location"), col = c("colour for this species/location")
legend("bottomleft",
legend = c("Bathymodiolus hybrid", "B. puteoserpentis"),
col = c("#ffb219","#cc0000"),
lty = 2,
lwd = 12,
title = "Symbionts of",
cex = 1)
legend("bottomright",
legend = c("1997", "2001"),
col = c("#afb6c2","#788497"),
lty = 2,
lwd = 12,
title = "Sampling year",
cex = 1)
dev.off() # close the PDF device