library(gplots)
library(RColorBrewer)
library(ggplot2)
library(pheatmap)
library(dplyr)
library(ape)
library(cluster)
#library(reshape)
#library(ggtree)
library(data.table)

Attaching package: ‘gplots’


The following object is masked from ‘package:stats’:

    lowess



Attaching package: ‘dplyr’


The following objects are masked from ‘package:stats’:

    filter, lag


The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union



Attaching package: ‘data.table’


The following objects are masked from ‘package:dplyr’:

    between, first, last

allele.data <- read.csv("Diversity_of_CRISPR_loci_in_Escherichia.tsv", sep = '\t')#, fileEncoding="UTF-16LE")

rnames <- allele.data[,1]                            # assign labels in column 1 to "rnames"
mat_allele.data <- data.matrix(allele.data[,2:ncol(allele.data)])  # transform column 2-5 into a matrix
rownames(mat_allele.data) <- rnames                  # assign row names

#rnames <- baps.data[,1]                            # assign labels in column 1 to "rnames"
#mat_baps.data <- data.matrix(baps.data[,2:ncol(baps.data)])  # transform column 2-5 into a matrix
#rownames(mat_baps.data) <- rnames                  # assign row names

head(mat_allele.data)
dim(mat_allele.data)

A matrix: 6 × 7 of type int

	X2.1	X2.2	X2.3	X2.2.3	X4.1	X4.2	X4.1.2
SMS-3-5	23	0	0	0	0	0	1
101.1	7	3	16	0	0	0	2
H10407	4	3	8	0	0	0	2
K12	14	3	7	0	0	0	2
ATCC 8739	22	3	29	0	0	0	4
53638	10	3	5	0	0	0	4

1. 100
2. 7

#allelic profiles heatmap
#ap.hmap <- pheatmap(mat_allele.data, clustering_method = "complete", annotation_row = as.data.frame(mat_mic.data), cluster_rows=1, cluster_cols = 0, show_colnames=0, show_rownames=0)
ap.hmap <- pheatmap(mat_allele.data, clustering_method = "complete", cluster_rows=1, cluster_cols = 0, show_colnames=0, show_rownames=0)

cgMLST_dissim_mat <- as.matrix(daisy(mat_allele.data, metric='gower'))
head(cgMLST_dissim_mat)

Warning message in daisy(mat_allele.data, metric = "gower"):
“binary variable(s) 2 treated as interval scaled”

A matrix: 6 × 100 of type dbl

	SMS-3-5	101.1	H10407	K12	ATCC 8739	53638	IAI1	HS	55989	SE11	⋯	EC63 (B2)	EC64 (B2)	EC65 (B2)	EC66 (B2)	EC67 (B1)	EC68 (B1)	EC69 (B1)	EC70 (B1)	EC71 (B1)	EC72 (B1)
SMS-3-5	0.0000000	0.32643678	0.30131363	0.24876847	0.3761905	0.31510673	0.3274220	0.2904762	0.34548440	0.2796388	⋯	0.2055556	0.1809524	0.3277778	0.1523810	0.30164204	0.32216749	0.22939245	0.27766831	0.306239737	0.27290640
101.1	0.3264368	0.00000000	0.05369458	0.07766831	0.1926108	0.12561576	0.1152709	0.4359606	0.01904762	0.1055829	⋯	0.3796114	0.2978654	0.5018336	0.2692939	0.03431856	0.03382594	0.09704433	0.04876847	0.048768473	0.05353038
H10407	0.3013136	0.05369458	0.00000000	0.05254516	0.2463054	0.10049261	0.1689655	0.4108374	0.04417077	0.1592775	⋯	0.3259168	0.2441708	0.4481390	0.2155993	0.01937603	0.08752053	0.07192118	0.03349754	0.004926108	0.03825944
K12	0.2487685	0.07766831	0.05254516	0.00000000	0.2036125	0.08604269	0.1262726	0.3963875	0.09671593	0.1165846	⋯	0.3686097	0.2868637	0.4908320	0.2582923	0.05287356	0.07339901	0.01937603	0.02889984	0.057471264	0.02413793
ATCC 8739	0.3761905	0.19261084	0.24630542	0.20361248	0.0000000	0.17536946	0.0773399	0.4857143	0.21165846	0.1060755	⋯	0.5722222	0.4333333	0.6944444	0.4619048	0.22692939	0.15878489	0.20394089	0.21280788	0.241379310	0.20804598
53638	0.3151067	0.12561576	0.10049261	0.08604269	0.1753695	0.00000000	0.1551724	0.3103448	0.14466338	0.2026273	⋯	0.3968528	0.2579639	0.5190750	0.2865353	0.10082102	0.13087028	0.08571429	0.07684729	0.105418719	0.08160920

gd.hmap <- pheatmap(cgMLST_dissim_mat, cluster_rows=1, cluster_cols = 1, show_rownames=0, show_colnames=0)

#d <- dist(mat_allele.data, method = "manhattan") #gave 6 clusters
#h_clust <- hclust(d, method = "ward.D") #ward clustering gave 6 clusters
#h_clust <- hclust(d, method = "complete") #complete clustering gave too many clusters
cgMLST_dissim_mat_for_clust <- daisy(mat_allele.data, metric='gower')
h_clust <- hclust(cgMLST_dissim_mat_for_clust, method = "ward.D") #ward clustering gave 6 clusters
cgMLST_tree.newClustering <- as.phylo(hclust(cgMLST_dissim_mat_for_clust, method = "ward.D"))
plot(h_clust)
rect.hclust(h_clust,k=12)

Warning message in daisy(mat_allele.data, metric = "gower"):
“binary variable(s) 2 treated as interval scaled”

#extract clusters
clusters <- cutree(h_clust,k=12)
#save to file
write.csv(file = "Diversity_of_CRISPR_loci_in_Escherichia.gower.ward.D.results.12groups.csv",clusters)