Question: Metagenomeseq package; time series analysis
0
gravatar for pvp2101
4.3 years ago by
pvp21010
United States
pvp21010 wrote:

I was going to implement code to look across all genus to see what genus of bacteria differed between time periods.  We harvested oysters over the period of 4 months from 2 sites in WA. 

 

The code I modified for this analysis is here:

 

timeSeriesFits = lapply(genusi,function(i){

  fitTimeSeries(obj=obj2,

                feature=i,

                class= ??,

                id="Location",

                time="Collection_Date_cat",

                lvl='genus',

                B=10) # This should be set to a higher value

})

 

In this case, I converted the dates to numbers into the time variable.  The ID I understood as should be corresponding to the two sites over time, so I put the ID as the location.  However, I’m not sure what to put down for class, since it is a bit different from the mousedata example, where you used status which seems to represent nutritional status of the mouse.  I put the ID variable as the individual sample names and the class I changed as location, but its spitting out

Error in env$fun(x, wk, env$env) : gss error in rk: unknown factor levels

 

Thank you

 

ADD COMMENTlink modified 4.3 years ago by Joseph Nathaniel Paulson270 • written 4.3 years ago by pvp21010
Answer: Metagenomeseq package; time series analysis
1
gravatar for Joseph Nathaniel Paulson
4.3 years ago by
United States
Joseph Nathaniel Paulson270 wrote:

Hi pvp2101,

Thanks for the interest. Since it sounds like you're interested in comparing locations across time your class should be the "Location" variable. To ensure that the id's are unique (sample ids) I created a new variable called 'id' given the column names. Make sure to upgrade to version 1.9.30 as we've made a number of under the hood tweaks to improve stability.

genusi = fData(obj2)$genus
pData(obj2)$id = colnames(obj2)

timeSeriesFits = lapply(genusi,function(i){
  fitTimeSeries(obj=obj2,
                feature=i,
                class= "Location",
                id="id",
                time="Collection_Date_cat",
                lvl='genus',
                B=10) # This should be set to a higher value
})
res = sapply(timeSeriesFits,function(i)i[1]);
names(res) = genusi
res = res[-grep("No statistically",res)]

Best wishes,

ADD COMMENTlink modified 4.3 years ago • written 4.3 years ago by Joseph Nathaniel Paulson270
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 310 users visited in the last hour