pcoa questions

Commands in mothur
#1

:slight_smile: Hi everyone，
I have some questions about “pcoa” command.

  1. The output ofthis command are two files ending in *pcoa, and *pcoa.loadings,but I can’t find the file ending in *pcoa. However, there is a file ending in *axes.
  2. Most PCoA are based on unweighted unifrac or weighted unifrac. Is this command also based on unifrac? If not, how can I do PCoA based on unifrac?
  3. Who can give me a R script to plot PCoA results
    Thanks for your help! :slight_smile:
1 Like
#2

  1. The axes file contains the plotting coordinates for graphing the results so this is what you need. Maybe the *.pcoa is a typo in the logfile or something.

  2. PCoA is just run on any distance matrix. It’s up to you what you want to run it on, you just need a phylip-formatted matrix to proceed. As you note, you can get these from the unifrac.unweighted/unifrac.weighted commands, but also from the dist.shared command (Bray-Curtis, Yue-Clayton theta, Jaccard etc).

  3. If you want anything fancy you’ll need to have a play around a bit, but for a very basic plot:

pcoa_data = read.table(file="your.pcoa.file.axes", header=T, sep='\t')
plot(pcoa_data$Axis2~pcoa_data$Axis1,xlab='Axis 1',ylab='Axis 2',main='Blah blah PCoA...')

If you have a design file that matches the groups processed in your distance matrix you can use the following to get group-based colouring or symbols:

design_file = read.table(file="your.design.file", 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 ]] 
}

Then you could plot either by colours:

plot(pcoa_data$Axis2~pcoa_data$Axis1,xlab='Axis 1',ylab='Axis 2',main='Blah blah PCoA...',pch=21,col='black',bg=col_vector)
legend(X,Y,uniq_groups,avail_cols)

Or symbols:

plot(pcoa_data$Axis2~pcoa_data$Axis1,xlab='Axis 1',ylab='Axis 2',main='Blah blah PCoA...',pch=symbol_vector)

You’ll just need to manually work out the X/Y values for placing your legend (hopefully there will be a bit of an empty spot on your graph where you can stick it).

#3

Hi dwaite,
Thank you for the quick response. Your reply give me invaluable help! :slight_smile: