Metagenomeseq package; time series analysis
1
0
Entering edit mode
pvp2101 • 0
@pvp2101-7504
Last seen 9.1 years ago
United States

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

 

metagenomeseq software error • 1.1k views
ADD COMMENT
1
Entering edit mode
@joseph-nathaniel-paulson-6442
Last seen 7.1 years ago
United States

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 COMMENT

Login before adding your answer.

Traffic: 732 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6