rrshieldscutler/splinectomer

R package version of the splinectomy longitudinal statistical analysis tools

splinectomer's Issues

Input OTU table is Relative abundance or log transformed?

Hi @RRShieldsCutler,

Thanks for a great package! I am following the tutorial on my 16S data. But I was a bit confused about the structure of the OTU table. Could you please confirm if the input OTU table consists of relative abundance or log-transformed data?

Regards
Jigyasa

Different pvalue from permuspliner from multiple runs

@RRShieldsCutler

I am following the tutorial in https://rrshieldscutler.github.io/splinectomeR/yassour_antibiotics_web.html.

When using permuspliner with same data and executing the command multiple time, the pvalue differs from each run.

Below is the command and the pvalue is different from three different runs with same command and same data

> f__Bacteroidaceae.top_family_meta <- top_family_meta %>% filter(family == 'k__Bacteria.p__Bacteroidetes.c__Bacteroidia.o__Bacteroidales.f__Bacteroidaceae')
> bacteroidaceae_result <- permuspliner(data = f__Bacteroidaceae.top_family_meta, xvar = 'timepoint',
+                        yvar = 'relative_abundance', perms = 99,
+                        category = 'antibiotics_y_n', cases = 'baby_id', retain_perm = T)

Testing between N and Y for a difference in relative_abundance

Scalpel please: performing permusplinectomy with 99 permutations...

p-value = 0.31

> f__Bacteroidaceae.top_family_meta <- top_family_meta %>% filter(family == 'k__Bacteria.p__Bacteroidetes.c__Bacteroidia.o__Bacteroidales.f__Bacteroidaceae')
> bacteroidaceae_result <- permuspliner(data = f__Bacteroidaceae.top_family_meta, xvar = 'timepoint',
+                        yvar = 'relative_abundance', perms = 99,
+                        category = 'antibiotics_y_n', cases = 'baby_id', retain_perm = T)

Testing between N and Y for a difference in relative_abundance

Scalpel please: performing permusplinectomy with 99 permutations...

p-value = 0.42

> f__Bacteroidaceae.top_family_meta <- top_family_meta %>% filter(family == 'k__Bacteria.p__Bacteroidetes.c__Bacteroidia.o__Bacteroidales.f__Bacteroidaceae')
> bacteroidaceae_result <- permuspliner(data = f__Bacteroidaceae.top_family_meta, xvar = 'timepoint',
+                        yvar = 'relative_abundance', perms = 99,
+                        category = 'antibiotics_y_n', cases = 'baby_id', retain_perm = T)

Testing between N and Y for a difference in relative_abundance

Scalpel please: performing permusplinectomy with 99 permutations...

p-value = 0.28

Using more than one category

@RRShieldsCutler
Can you please let me know if it is possible to test simultaneously for more than one category across the time points ?
Example, for the following metadata from tutorial, is it possible to use both categories "birth_mode" and "feeding_in_hospital" simultaneously across all the time points.

> head(general_metadata)
  #Sample_ID birth_year    sex        birth_mode gestational_age
1    E000823       2008   Male  Vaginal delivery            39+4
2    E001958       2008 Female  Vaginal delivery            40+1
3    E003188       2008 Female  Vaginal delivery            38+6
4    E003953       2008 Female  Vaginal delivery            39+5
5    E004628       2008   Male Ceasarean section            40+3
6    E004709       2008   Male  Vaginal delivery            40+2
  months_final_exclusive_breastfeeding months_end_breastfeeding
1                                  5.0                    13.02
2                                  0.0                    12.66
3                                  5.0                    13.28
4                                  3.0                     6.98
5                                  0.5                     4.46
6                                  0.0                     7.97
  months_started_formula_feeding feeding_in_hospital
1                             NA      OnlyBreastmilk
2                             NA               mixed
3                           10.0      OnlyBreastmilk
4                            3.0      OnlyBreastmilk
5                            0.5      OnlyBreastmilk
6                            0.0               mixed

Thanks

How does splinectomeR handle missing values for a subject?

I am working with a dataset with 4 timepoints. Not every subject has a sample for each timepoint. How does splinectomeR handle missing values for a subject?

x-axis ticks in trendyspliner

Hi @RRShieldsCutler
Thank you for this very useful package.

I am wondering if there is a reason for the limits of the x-axis ticks in trendyspliner.

I ran trendyspliner with the ChickWeight example data like this:

data("ChickWeight")
set.seed(100)
result <- trendyspliner(data = ChickWeight, x = 'Time', y = 'weight', cases = 'Chick', category = 'Diet', group = '1', perms = 999)

and then made a plot like this:

result_plot <- trendyspliner.plot.perms(result)
result_plot

Then I looked for the x-axis limits of the plot:

layer_scales(result_plot)$x$range$range

This yields a min of 1.0500 and a max of 20.0025

However, range(ChickWeight$Time) is from 0 to 21.

Is there a reason that this does not match?

I had the same issue with my own data, and when I replace the default x-axis ticks, the first and last label are not covered by the plot.

Thanks in advance for your help.

trendyspliner returns pvalues either 0.001 or 0.002

@RRShieldsCutler
I am using trendyspliner to test a response in one group (here, feeding_in_hospital == "OnlyFormula") over time for each of the bacterium family. When i run through the method using the commands as shown below, i get pvalue either 0.001 or 0.002 only. When i change the feeding_in_hospital == "mixed", i do get pvalues either 0.001 or 0.002. I am little confused on how to interpret the results. Not sure if the input or parameters are correct. Can you please help to understand the output.
Thanks

top20_family_meta <- top20_family_meta[which(top20_family_meta$feeding_in_hospital == "OnlyFormula"),]


f_top20_abx_pvals <- list()
for (f in unique(top20_family_meta$family)) {
  f__df <- top20_family_meta %>% filter(family == f)
  
  set.seed(123456)
  f.result <- trendyspliner(data = f__df,
                               xvar = 'timepoint', quiet = T, cases = 'baby_id',
                               mean_center = T, yvar = 'relative_abundance', perms = 999)
  message(f, " : ",f.result$pval)
  f_top20_abx_pvals <- append(f_top20_abx_pvals, f.result$pval)
}

Output

k__Archaea.p__Euryarchaeota.c__Methanobacteria.o__Methanobacteriales.f__Methanobacteriaceae : 0.002
k__Bacteria.p__Actinobacteria.c__Actinobacteria.o__Actinomycetales.f__Actinomycetaceae : 0.002
k__Bacteria.p__Actinobacteria.c__Actinobacteria.o__Actinomycetales.f__Corynebacteriaceae : 0.001
k__Bacteria.p__Actinobacteria.c__Actinobacteria.o__Actinomycetales.f__Microbacteriaceae : 0.002
k__Bacteria.p__Actinobacteria.c__Actinobacteria.o__Actinomycetales.f__Micrococcaceae : 0.002
k__Bacteria.p__Actinobacteria.c__Actinobacteria.o__Bifidobacteriales.f__Bifidobacteriaceae : 0.001
k__Bacteria.p__Actinobacteria.c__Coriobacteriia.o__Coriobacteriales.f__Coriobacteriaceae : 0.001
k__Bacteria.p__Bacteroidetes.c__Bacteroidia.o__Bacteroidales.f__Bacteroidaceae : 0.001
k__Bacteria.p__Bacteroidetes.c__Bacteroidia.o__Bacteroidales.f__Porphyromonadaceae : 0.001
k__Bacteria.p__Bacteroidetes.c__Bacteroidia.o__Bacteroidales.f__Prevotellaceae : 0.002

sliding spliner not significant any time point

On doing permuspliner test between groups, I got a p value of 0.034,

#code
result <- permuspliner(data = data2, xvar = 'Time_point', yvar = 'y', cases = 'Subject_ID', category = 'Residence', groups = c('0','1'))

Groups detected: 0 and 1.

Now testing between variables 0 and 1 for a difference in the response labeled y

Scalpel please: performing permusplinectomy with 999 permutations...

Area between groups successfully calculated, now spinning up permutations...
...permutations completed...

p-value = 0.033_

On doing the slidingspliner, I got the plot showing not significant at any time point. I was expecting the difference at 4 week, but the pvalue showed no significance. How do I explain that,?

#code
slide_result <- sliding_spliner(data = data2,
xvar = 'Time_point', yvar = 'y',
category = 'Residence', groups =c('0','1'), cases = 'Subject_ID', cut_low = 3,
ints = 3)

E1 = sliding_spliner.plot.pvals(data = slide_result, xvar = 'Time_point')+
scale_x_continuous(breaks = c(0.2,1,1.8), labels = c("1 week", "4 weeks", "24 weeks"))

Looking forward to your guidance

Thanks
Sanuja

Help understanding error message in permuspliner

Hi,

I'm trying to follow the protocol presented in your vignette with some of my own data and I get the following error:

permuspliner(data = sphingobium, xvar = 'Date', yvar = 'Abundance', category = 'Category', cases = 'Identifier', groups=c('SS-25_well', 'Bkg_surr'), cut_low = 4, perms = 99)

Groups detected: SS-25_well and Bkg_surr .

Now testing between variables SS-25_well and Bkg_surr for a difference in the response labeled Abundance 

Scalpel please: performing permusplinectomy with 99 permutations...

Area between groups successfully calculated, now spinning up permutations...
Error in smooth.spline(x = randy_v1[, xvar], y = randy_v1[, yvar], spar = set_spar,  : 
  need at least four unique 'x' values

I'm not sure here what the "x values" are. In my original data I believe I have more than 4 dates, which is what I'm passing in as the xvar to permuspliner

For reference, I'm attaching my data:
sphingobium.rds.gz

I'm sure it's something simple I've misunderstood.

Best,
Connor

issue with p values from trendyspliner

Hello,

When I feed trendyspliner a vector of 0's it gives me a significant p value; in fact, it gives me a significant p value in all cases I have tried. I am using the latest version installed from github. Permuspliner will give a non significant value if fed a vector of 0's (no variation between two groups).

Error in seq.default(from = x.min, to = x.max, by = ((x.max - x.min)/(ints - : 'from' must be a finite number

Trying to do a longitudinal analysis of Shannon diversity over time. But only 3 time points- week 1, week 4, week 24

slide_result <- sliding_spliner(data = div.df,
xvar = 'Time', yvar = 'diversity_shannon',
category = 'Residence', groups =c('Farm','Control'), cases = 'Subject_ID')

Error in seq.default(from = x.min, to = x.max, by = ((x.max - x.min)/(ints - :
'from' must be a finite number

Time-series data doubt in category

Hi,
Thank you very much for this wonderful package. I don't have any specific query related to something not working properly in the package.

But I have a query of a time-series data from subjects followed for 3 months at different time intervals. And as a category we have a positive and negative report for bacterial infection in which the patient can be positive for infection at some time-point while at other negative.

So there is a within category variation in and across patients. I even tried running and as expected ran with error saying ""Please ensure that the cases are uniquely labeled between the two groups"".

I would appreciate if you can provide any suggestion on how can we use this kind of data for analysis in this package.

Thanks
Shriram

Error with sliding_spliner()

When I use sliding_spliner(), I get the following error:

Warning in min(as.numeric(df[, xvar])) :
no non-missing arguments to min; returning Inf
Warning in max(as.numeric(df[, xvar])) :
no non-missing arguments to max; returning -Inf
Error in seq.default(from = x.min, to = x.max, by = ((x.max - x.min)/(ints - :
'from' must be a finite number

Thank you in advance for your help!