Hi all,
I got the pcoa results from mothur, named “xx.pcoa.axes”. I would like to plot the results as the 3d pcoa plot. Thanks to the script provided by dwaite, I was able to plot 2d pcoa. The code as below, in case it might be useful to someone else
pcoa_data = read.table(file="10_horsepick.thetayc.0.03.lt.pcoa.axes", header=T, sep='\t')
design_file = read.table(file="10_horse.design", sep='\t',stringsAsFactors=FALSE)
uniq_groups = unique(design_file[,2])
num_uniq_groups = length(uniq_groups)
avail_cols = rainbow(num_uniq_groups)
avail_symbols = 1:num_uniq_groups
colour_map = list()
symbol_map = list()
for(i in 1:num_uniq_groups) {
current_group = uniq_groups[i]
colour_map[[ current_group ]] = avail_cols[i]
symbol_map[[ current_group ]] = avail_symbols[i]
}
col_vector = rep(NA,nrow(design_file))
symbol_vector = rep(NA,nrow(design_file))
for(i in 1:nrow(design_file)) {
group = design_file[i,2]
col_vector[i] = colour_map[[ group ]]
symbol_vector [i] = symbol_map [[ group ]]
}
plot(pcoa_data$axis2~pcoa_data$axis1,xlab='Axis 1',ylab='Axis 2',main='PCoA',col='black', pch=symbol_vector)
legend(X,Y,uniq_groups,avail_cols)
It would be nice to assign different colors to each group as well.
Thanks many!!