0
$\begingroup$

I have some data from mass spectrometry of some plant samples (11 morphospecies and 4 treatments for each). Each species/treatment yields a graphic like the one below.

Graphic from mass spectrometry

Each column in the graphic indicates how much of that chemical substance is present in the sample. All the data was assembled in a single table, as shown below.

sp  technique   mz  abundance
sp1 ESIneg  118.89  3.01
sp1 ESIneg  172.72  3.20
sp1 ESIneg  202.94  3.80
sp1 ESIpos  118.30  2.59
sp1 ESIpos  170.68  3.13
sp1 ESIpos  257.97  3.28
sp2 ESIneg  132.33  22.22
sp2 ESIneg  211.84  3.87

The table has 1587 rows and 4 columns. Now, I would like to obtain a dendrogram like the following:

enter image description here

which was obtained from this simple example in R:

d <- dist(as.matrix(mtcars))   # find distance matrix
hc <- hclust(d)                # apply hirarchical clustering
plot(hc)                       # plot the dendrogram

The problem is: mtcars in this example is a table with no repeating rows. While in my case, there each species/treatment spans many rows (one for each column in the first graphic), I'm in a complete lost about how to proceed. In a related question, I've been told to use dummy variables through model.matrix, but this still don't give me the sort of data I need.

How can i get a dendrogram showing the relatedness of all morphospecies, from my initial data?

$\endgroup$
  • $\begingroup$ What are objects to cluster? ESIneg is one, ESIpos is another, etc.? If yes how many such do you have? Different rows for ESIneg - how are they different, can you just avarage them somehow? How do you see the similarity between ESIneg and ESIpos? $\endgroup$ – ttnphns Jan 11 '17 at 17:35
  • $\begingroup$ @ttnphns The objects may be either the species or species+treatment (species+technique). I intend to do both analysis (only species and species together with treatment). Not sure how I'll "average" all treatments yet. I have 11 species and 4 such treatments (APCI + and -, ESI + and -). $\endgroup$ – Rodrigo Jan 12 '17 at 23:42
  • $\begingroup$ You will need to compute a similarity measure based on which the dendrogramm can be computed. However you do it, you will need a pairwaise distance matrix to do hierarchical clustering (--> dendrogramm). Which similarity measure/distance metric would be fit is a question that can only be answered by domain experts. $\endgroup$ – Nikolas Rieble Jan 18 '17 at 9:04
1
+50
$\begingroup$

Whereas the example you provide uses a distance matrix between all observations dist(as.matrix(mtcars)), you currently have a matrix of raw data, which is not formatted properly.

What you need to do in the first step is: Format your table so that

  • for each object to be clustered (morphospecies?), all measurements are in a single row
  • each type of measurement is in a single column.

In the next step, you will want to compute the similarity (or as you say: relatedness) between each pair of objects. This is the distance matrix d. Therefore you will have to choose a distance metric / measure of similarity. The choice depends on which characteristics you want to capture. You could start with Euclidean Distance and later use more sophisticated measures.

If you have this matrix, then you can simply use the commands in R that you presented:

hc <- hclust(d)                # apply hirarchical clustering
plot(hc)                       # plot the dendrogram

and the dendrogramm will be yours. Note here: I do not have domain knowledge, nor did you explain the abbreviations. Therefore this example might be wrong, yet it could help to illustrate. I will try to explain this using the table you provided:

Assuming sp means species and sp1, sp2, ... are the objects to be clustered.

sp  technique   mz  abundance
sp1 ESIneg  118.89  3.01
sp1 ESIneg  172.72  3.20
sp1 ESIneg  202.94  3.80
sp1 ESIpos  118.30  2.59
sp1 ESIpos  170.68  3.13
sp1 ESIpos  257.97  3.28
sp2 ESIneg  132.33  22.22
sp2 ESIneg  211.84  3.87

Now you want to have a single row for each sp such as:

The feature n_mz1 here means technique: ESIneg, measure: mz, measurement: 1

sp  n_mz1  n_mz2  n_mz3  p_mz1  p_mz2  p_mz3  n_ab_1 n_ab_2  ...
sp1 118.89 172.72 202.94 118.30 170.68 257.97 3.01   3.20    ...
sp2 ....

Or you could instead compute features based on the all mz and ab measurements for a single species such as average mz and average ab:

sp  avg_mz avg_ab 
sp1 118.89 172.72
sp2 ....

Based on such a feature matrix, you can compute the distance exactly as

d <- dist(as.matrix(mtcars))
$\endgroup$
  • $\begingroup$ Thank you, it worked perfectly! I'm gonna post an answer just to show exactly how I achieved your suggestion. $\endgroup$ – Rodrigo Jan 18 '17 at 15:11
  • $\begingroup$ That is a very good idea. $\endgroup$ – Nikolas Rieble Jan 18 '17 at 15:12
0
$\begingroup$

Following Nikolas answer, I did the following.

The mz values (mass/electric charge) range from 102 up to 999 (atomic units). So I have ~900 values for each treatment * 4 treatments = 3600 columns.

db <- read.csv('my_data.csv',sep='\t')
db$mzR <- round(db$mz) # rounded the values of mz to better group them in columns
db$an2 <- NA #simplifying 
db$an2[which(db$analysis == 'ESI negativo')] <- 'ESIn'
db$an2[which(db$analysis == 'ESI positivo')] <- 'ESIp'
db$an2[which(db$analysis == 'APCI negativo')] <- 'APCIn'
db$an2[which(db$analysis == 'APCI positivo')] <- 'APCIp'
m <- matrix(NA,length(unique(db$sp)),3600)
rownames(m) <- unique(db$sp)
colnames(m) <- c(paste('ESIn',101:1000,sep=''),
                 paste('ESIp',101:1000,sep=''),
                 paste('APCIn',101:1000,sep=''),
                 paste('APCIp',101:1000,sep=''))
# Example for the first row: m['sp1','ESIn119'] <- 3.01
for (i in 1:nrow(db)) {
  m[db$sp[i],paste(db$an2[i],db$mzR[i],sep='')] <- db$abundance[i]
}
d <- dist(m)
hc <- hclust(d) # apply hirarchical clustering
plot(hc) # dendro1

This produced a matrix m with lots of NAs. So I used the Bray-Curtis distance, which is more appropriate for tables with lots of empty spaces (first I turned all NAs to zeros).

install.packages('vegan')
library(vegan)
m2 <- m
m2[is.na(m2)] <- 0
d2 <- vegdist(m2,method='bray')
hc2 <- hclust(d2)
plot(hc2) # dendro2

Another approach was to group columns in groups of 10. So all mz from 110.0 up to 119.99 became 11.

m3 <- matrix(0,length(unique(db$sp)),364)
rownames(m3) <- unique(db$sp)
colnames(m3) <- c(paste('ESIn',10:100,sep=''),
                 paste('ESIp',10:100,sep=''),
                 paste('APCIn',10:100,sep=''),
                 paste('APCIp',10:100,sep=''))
for (i in 1:nrow(db)) {
  m3[db$sp[i],paste(db$an2[i],trunc(db$mz[i]/10),sep='')] <- 
    m3[db$sp[i],paste(db$an2[i],trunc(db$mz[i]/10),sep='')] + db$abundance[i]
}
d3 <- vegdist(m3,method='bray')
hc3 <- hclust(d3)
plot(hc3) # dendro3
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.