The envfit function needs two pieces of data, your initial data matrix and a table of metadata that groups the individual samples according to some environmental groupings. Usually in R the hardest part of an analysis is getting your data into R, once you have the data the actual analysis is very simple.

Since you've done your CCA already, you should have a matrix of your data, probably an OTU table. This should be in the format of samples as rows, and OTU abundances as columns. You then need a table that lists your environmental data (pH, temperature, salinity or whatever), where the rows of the table correspond to the rows of your OTU table, and each column reports a different environmental factor.

As a completely made-up example, this code will make a quick OTU table called 'myTable'. In this table the first 10 rows contain values between 1 and 10, and rows 11 - 20 contain between 15 and 25. This is meaningless data, but gives a strong difference between these row groupings.

myTable = matrix(nrow=20,ncol=10)

for(i in 1:10) {myTable[i,] = sample(1:10,10)}

for(i in 11:20) {myTable[i,] = sample(15:25,10)}

Then this code will make a data frame that contains some metadata about the OTU table. In this, the first column of my metadata table groups the rows into two groups, 1-10 and 11-20. The second column groups all the odd-numbered rows together and all the even-numbered rows together. Again, this is meaningless but gives us one metadata grouping that should provide a very strong signal and one with a very weak, insignificant signal.

myMetadata = data.frame(FirstColumn = character(), SecondColumn = character(), stringsAsFactors=FALSE)

myMetadata[1:10,1] = 'G1'

myMetadata[11:20,1] = 'G2'

myMetadata[seq(1,20,2),2] = 'G3'

myMetadata[seq(2,20,2),2] = 'G4'

To run the envfit, all you need to do is call the command(s):

library(vegan)

envfit(myTable ~ myMetadata$FirstColumn, data = myMetadata, perm=1000)

This will perform the envfit function against your OTU table, grouping the individual samples (rows) according to how they are grouped in the first column of myMetadata. The perm parameter specifies how many random permutations to test for assigning significance to the vectors. You should see a high R2 value with strong statistical significance (p-value). If you test the data using the other column:

envfit(myTable ~ myMetadata$SecondColumn, data = myMetadata, perm=1000)

You should see an extremely weak (~0) R2 with no significance.

With your own data, you'll just need to step through each metadata column in order to record the fit/significance of each parameter. I hope that helps. If you need anything else just let me know.

The envfit function needs two pieces of data, your initial data matrix and a table of metadata that groups the individual samples according to some environmental groupings. Usually in R the hardest part of an analysis is getting your data into R, once you have the data the actual analysis is very simple.

Since you've done your CCA already, you should have a matrix of your data, probably an OTU table. This should be in the format of samples as rows, and OTU abundances as columns. You then need a table that lists your environmental data (pH, temperature, salinity or whatever), where the rows of the table correspond to the rows of your OTU table, and each column reports a different environmental factor.

As a completely made-up example, this code will make a quick OTU table called 'myTable'. In this table the first 10 rows contain values between 1 and 10, and rows 11 - 20 contain between 15 and 25. This is meaningless data, but gives a strong difference between these row groupings.

myTable = matrix(nrow=20,ncol=10)

for(i in 1:10) {myTable[i,] = sample(1:10,10)}

for(i in 11:20) {myTable[i,] = sample(15:25,10)}

Then this code will make a data frame that contains some metadata about the OTU table. In this, the first column of my metadata table groups the rows into two groups, 1-10 and 11-20. The second column groups all the odd-numbered rows together and all the even-numbered rows together. Again, this is meaningless but gives us one metadata grouping that should provide a very strong signal and one with a very weak, insignificant signal.

myMetadata = data.frame(FirstColumn = character(), SecondColumn = character(), stringsAsFactors=FALSE)

myMetadata[1:10,1] = 'G1'

myMetadata[11:20,1] = 'G2'

myMetadata[seq(1,20,2),2] = 'G3'

myMetadata[seq(2,20,2),2] = 'G4'

To run the envfit, all you need to do is call the command(s):

library(vegan)

envfit(myTable ~ myMetadata$FirstColumn, data = myMetadata, perm=1000)

This will perform the envfit function against your OTU table, grouping the individual samples (rows) according to how they are grouped in the first column of myMetadata. The perm parameter specifies how many random permutations to test for assigning significance to the vectors. You should see a high R2 value with strong statistical significance (p-value). If you test the data using the other column:

envfit(myTable ~ myMetadata$SecondColumn, data = myMetadata, perm=1000)

You should see an extremely weak (~0) R2 with no significance.

With your own data, you'll just need to step through each metadata column in order to record the fit/significance of each parameter. I hope that helps. If you need anything else just let me know.

Hi David

Thank you for taking time to help! Despite you very detailed answer, I'm not sure I understood everything.

Yes I have 2 tables, one with the individual and environmental factors (10 columns), end one with the OTUs (63 columns). Each table refers to the same rows = individuals. What I don't understand is how you create your metadata table. In your example, you choose arbitrary to group 1 to 10 and 11 to 20. But in my case, how am I suppose to choose? based on CCA analysis? Here are my results:

> CCorA(Factors, Order)

Information - Matrix1: rank=10 < order11

Information - Matrix2: rank=20 < order63

Canonical Correlation Analysis

Call:

CCorA(Y = Factors, X = Order)

                           Y X

Matrix Ranks 10 20

Pillai's trace: 3.996365

Significance of Pillai's trace:

from F-distribution: 0.0086079

                                         CanAxis1 CanAxis2 CanAxis3 CanAxis4 CanAxis5 CanAxis6 CanAxis7 CanAxis8 CanAxis9 CanAxis10

Canonical Correlations 0.85302 0.81898 0.78536 0.72082 0.63851 0.57999 0.48636 0.47677 0.42233 0.2744

                                    Y | X X | Y

RDA R squares 0.45066 0.1790

adj. RDA R squares 0.17599 0.0148

Usually it would be based on some pre-defined information about your samples. For example, in the paper you linked the researchers took a patient history for each individual and then created groupings based on these (antibiotic usage, vaccination, sex, infection history etc). If you give me some context about the data you're working with I could give you some specific examples, but ideas you could try include pH, temperature, salinity (if environmental samples), or diet, sex, age (if animal samples). Some factors require sample measurements, others you might be able to deduce from the literature.

If you don't have any of this information on hand, there are other techniques you can use to try and let your data cluster naturally, and then look for commonalities between the samples in the clusters.

But in your example you don't group columns, you group rows. I have about 200 hundred surveyed people and I don't see how am I supposed to group them.

The idea is to keep all the variables sex, age, infection... and to make a selection ONLY on the OTUs (to keep a dozen instead of 100). The second step would be a multivariate analysis with the selected variables, as in the paper. I'm still in the first step, i.e. preselection of relevant OTUs

Ah, sorry. I (clearly) misunderstood the question. 

The easiest way I can think of to extract the top XXX OTUs would be to calculate the total or average abundance of each OTU (your choice), then to filter the OTU table based on the top XXX. Assuming that 'Order' is your OTU table, you could get the top 12 OTUs (by summed abundance) using the following

Temp_Order = addmargins(Order, FUN = list(Total = sum), margin=1)

top_twelve_columns = order(Temp_Order [5,],decreasing=TRUE)[1:12]

New_Order = Order[,top_twelve_columns]

Where 'New_Order' contains all the rows, but only the 12 most abundant columns. You can change "sum" in the first line to mean, median etc if you want to use a different metric to sort on.

For the second part, you can then put this matrix straight into envfit with your existing Factors table. I actually didn't realise this, but if you just dump the entire Factors table into envfit it will work out the correlation for each column in the Factors table:

Env_Results = envfit(New_Order, Factors)

If you then type "Env_Results$vectors" into the console you'll get back a table with the R2 and p-value for each of your Factor columns. Since you've already done your analysis with Factors/Order the column/row relationship should be correct already, but you can always transpose them using

Order = t(Order)

If need be.

Oops, there's a typo in the second line.

top_twelve_columns = order(Temp_Order [5,],decreasing=TRUE)[1:12]

Should be

 top_twelve_columns = order(Temp_Order [nrow(Temp_Order),],decreasing=TRUE)[1:12]

yes, but I wanted to select the OTUs regarding their importance in terms of correlation, and not only in terms of relative frequency/abundance in the population. In the envfit function:

envfit(ord,env,permutations=999, strata=NULL,choices=c(1,2),display="sites")

I could not figure out what was the ord and env

but simply trying with my 2 groups of data, i.e. envfit(Factors,Order,permutations=999, strata=NULL,choices=c(1:7),display="sites") it seems that it works!!! and I can simply select the OTUs with a p-value

why is the permutations always 999?

Because the real dataset (the original one, from which the permutations are drawn) is the 1000th ;)

David W Waite thanks for your response. I know the steps to perform the analysis, but I confused the order of rows and collumns. I want to fit two treatments (warmer vs actual temp) that could affect germination speed (days) of certain groups of seedlings. Is it valid to repeat the days or the species (as the case may be) if each have germination in both treatments in each row to be able to develop a fit with envfit function?

Regards,

Criss

Dear David W Waite ,

I just ran the codes of the your 'December 1, 2014' I get the same results.

However I ran into a problem with my own data.

Although the code seems to analyse all the different OTU's/columns it is actually only analyzing the first one collumn of 'myTable'. (in my own dataset it is only analyzing the first two OTU's).

You can check this by changing column 2 to 10 in 'myTable' to values between 15 and 25

>for(i in 1:10) {myTable[i,2:10] = sample(15:25,9)}

Now OTU1/colums 1 G1 has values ranging 1 to 10

Now OTU1/colums 1 G2 has values ranging 15 to 25

All other OTUs/columns have values ranging 15 to 25 (so there is no difference between G1 and G2 in these columns)

When you now run the envfit code

>envfit(myTable ~ myMetadata$FirstColumn, data = myMetadata, perm=1000)

it will give you the exact same result, the envfit is only analysing the first collumn of myTable

(if you change all values of OTU1 to 15:25 instead of collumn 2 to 10 you won't find a significant difference.)

Do you know how to analyse all the OTU's at the same time?

Cheers, Rik