parallel_extract_core() behaviour

[Gian77]

This is what I got with a reduced size test_phyloseq.rda dataset. It does not sound right to me, the elbow method should be more stringent as described by Shade and Stopnisek. However, I am not sure if it will be for any dataset, please see below @B-Kristy @B-Kristy

I reduce the samples for speed...

test_small_phyloseq <- 
    test_phyloseq %>%
    subset_samples(site %in% "ESC") %>% 
    prune_samples(sample_sums(.) >= 1, .)

rarefied_data <- 
    parallel_rarefy(physeq = test_small_phyloseq, 
                    depth_level = 500, 
                    num_iter = 9, 
                    threads = 8, 
                    set_seed = 100)

physeq_samll_rare <- 
    do_phyloseq(physeq = test_small_phyloseq, 
                otu_rare = rarefied_data )

> core_ps_rare_inc<- parallel_extract_core(
+     physeq_samll_rare,
+     Var = "crop",
+     method = "increase",
+     increase_value = 2,
+     ncores = 10, 
+     set_seed = 321
+ ) 
Seed used: 321
✔ Input phyloseq object is valid!
ℹ Dataset is rarefied at a depth of: 500
ℹ Group: not used.
ℹ Core calculated across: crop
ℹ Calculating BC dissimilarity based on the 1st ranked OTU
✔ BC dissimilarity based on the 1st ranked OTU complete
ℹ Calculating BC dissimilarity based on ranked OTUs, starting at 2025-09-10 11:54:10.495006
ℹ Running in serial
✔ BC ranks done!                                 
ℹ Performing method 'increase'
✔ % increase method identified 4 core OTUs
✔ Analysis complete!

Got 4 core OTUs with increase method.

> core_ps_rare_elb <- parallel_extract_core(
+     physeq_samll_rare,
+     Var = "SampleType",
+     method = "elbow",
+     ncores = 10,
+     set_seed = 567
+ ) 
Seed used: 567
✔ Input phyloseq object is valid!
ℹ Dataset is rarefied at a depth of: 500
ℹ Group: not used.
ℹ Core calculated across: SampleType
ℹ Calculating BC dissimilarity based on the 1st ranked OTU
✔ BC dissimilarity based on the 1st ranked OTU complete
ℹ Calculating BC dissimilarity based on ranked OTUs, starting at 2025-09-10 11:55:11.251187
ℹ Running in serial
✔ BC ranks done!                                 
ℹ Performing method 'elbow'
✔ Elbow method identified 112 core OTUs
✔ Analysis complete!

Got 112 OTUs with elbow mehod.

[jibarozzo]

@Gian77 What version or parallel_extract_core() are you using? What branch. I tried to test in your workflow_improve branch and the version is different.

[jibarozzo]

parallel_extract_core<- function(
physeq,
Var,
method,
increase_value = NULL,
Group = NULL,
Level = NULL,
trimOTUs = TRUE,
.parallel = FALSE,
ncores = get_available_cores()
)

[jibarozzo]

@Gian77 is referring to 898a809

[jibarozzo]

@Gian77 The first thin I see is that you are using Var = "SampleType" in the elbow method vs. "crop" in increase. Apples to oranges. Not that it changes much.

core_ps_rare_elb <- parallel_extract_core(
+     physeq_samll_rare,
+     Var = "crop",
+     method = "elbow",
+     .parallel = TRUE,
+     ncores = 4,
+     set_seed = 567
+ ) 
Seed used: 567
✔ Input phyloseq object is valid!
ℹ Dataset is rarefied at a depth of: 500
ℹ Group: not used.
ℹ Core calculated across: crop
ℹ Calculating BC dissimilarity based on the 1st ranked OTU
✔ BC dissimilarity based on the 1st ranked OTU complete
ℹ Calculating BC dissimilarity based on ranked OTUs, starting at 2025-09-10 11:25:02.842172
ℹ Local environment. Using 7 core(s)/worker(s) via parallelly::availableCores().
Keeping 1 core/worker available for system.
✔ BC ranks done!
ℹ Performing method 'elbow'
✔ Elbow method identified 112 core OTUs
✔ Analysis complete!

[Gian77]

Good catch @jibarozzo . I change from GlobalPatterns to physeq_samll_rare, forgot to change the Var. I rerun with crop on both, this below what I got, and seems more reasonable. I will rerun with the full dataset later.

> core_ps_rare_inc<- parallel_extract_core(
+     physeq_samll_rare,
+     Var = "crop",
+     method = "increase",
+     increase_value = 2,
+     ncores = 10, 
+     set_seed = 888
+ ) 
Seed used: 888
✔ Input phyloseq object is valid!
ℹ Dataset is rarefied at a depth of: 500
ℹ Group: not used.
ℹ Core calculated across: crop
ℹ Calculating BC dissimilarity based on the 1st ranked OTU
✔ BC dissimilarity based on the 1st ranked OTU complete
ℹ Calculating BC dissimilarity based on ranked OTUs, starting at 2025-09-10 12:30:05.58101
ℹ Running in serial
✔ BC ranks done!                                 
ℹ Performing method 'increase'
! No OTUs meet the specified increase threshold. Using top OTU only.
✔ % increase method identified 1 core OTUs
✔ Analysis complete!

> core_ps_rare_elb <- parallel_extract_core(
+     physeq_samll_rare,
+     Var = "crop",
+     method = "elbow",
+     ncores = 10,
+     set_seed = 345
+ ) 
Seed used: 345
✔ Input phyloseq object is valid!
ℹ Dataset is rarefied at a depth of: 500
ℹ Group: not used.
ℹ Core calculated across: crop
ℹ Calculating BC dissimilarity based on the 1st ranked OTU
✔ BC dissimilarity based on the 1st ranked OTU complete
ℹ Calculating BC dissimilarity based on ranked OTUs, starting at 2025-09-10 12:30:22.652519
ℹ Running in serial
✔ BC ranks done!                                 
ℹ Performing method 'elbow'
✔ Elbow method identified 1 core OTUs
✔ Analysis complete!

[jibarozzo]

This issue has also highlighted how the parallel helper function get_available_cores() is overriding the cores specidified by the user. I will tackle this when I implement changes mentioned in #46. I specified .parallel = TRUE, ncores = 4 and received this message: Local environment. Using 7 core(s)/worker(s) via parallelly::availableCores(). Keeping 1 core/worker available for system.

[jibarozzo]

@Gian77

> core_ps_rare_elb <- parallel_extract_core(
+     physeq_samll_rare,
+     Var = "crop",
+     method = "elbow",
+     .parallel = TRUE,
+     ncores = 4,
+     set_seed = 567
+ ) 
Seed used: 567
✔ Input phyloseq object is valid!
ℹ Dataset is rarefied at a depth of: 500
ℹ Group: not used.
ℹ Core calculated across: crop
ℹ Calculating BC dissimilarity based on the 1st ranked OTU
✔ BC dissimilarity based on the 1st ranked OTU complete
ℹ Calculating BC dissimilarity based on ranked OTUs, starting at 2025-09-10 11:38:32.123688
ℹ Local environment. Using 7 core(s)/worker(s) via parallelly::availableCores().
Keeping 1 core/worker available for system.
✔ BC ranks done!
ℹ Performing method 'elbow'
✔ Elbow method identified 1 core OTUs
✔ Analysis complete!

[Gian77]

With the updated version (with chatGTP help) to remove tibble indexing facilitating macos to work, I still have this weird results that elbow generates a larger core compared to increase and I do not think is realistic.

> core_ps_rare_inc<- parallel_extract_core(
+     physeq_samll_rare,
+     Var = "crop",
+     method = "increase",
+     increase_value = 2,
+     ncores = 10, 
+     set_seed = 888
+ ) 
Seed used: 888
✔ Input phyloseq object is valid!
ℹ Dataset is rarefied at a depth of: 500
ℹ Group: not used.
ℹ Core calculated across: crop
ℹ Calculating BC dissimilarity based on the 1st ranked OTU
✔ BC dissimilarity based on the 1st ranked OTU complete
ℹ Calculating BC dissimilarity based on ranked OTUs, starting at 2025-09-11 14:54:26.040625
ℹ Running in serial
✔ BC ranks done!                                 
ℹ Performing method 'increase'
! No OTUs meet the specified increase threshold. Using top OTU only.
✔ % increase method identified 4 core OTUs
✔ Analysis complete!


> core_ps_rare_elb <- parallel_extract_core(
+     physeq_samll_rare,
+     Var = "crop",
+     method = "elbow",
+     ncores = 10,
+     set_seed = 345
+ ) 
Seed used: 345
✔ Input phyloseq object is valid!
ℹ Dataset is rarefied at a depth of: 500
ℹ Group: not used.
ℹ Core calculated across: crop
ℹ Calculating BC dissimilarity based on the 1st ranked OTU
✔ BC dissimilarity based on the 1st ranked OTU complete
ℹ Calculating BC dissimilarity based on ranked OTUs, starting at 2025-09-11 14:55:03.913443
ℹ Running in serial
✔ BC ranks done!                                 
ℹ Performing method 'elbow'
✔ Elbow method identified 245 core OTUs
✔ Analysis complete!

[Gian77]

ok, why this is happening? @jibarozzo @jibarozzo I downloaded my original extract_core_function from here and compared with the parallel_extract_core and run both. Need to go step by step and see what's going on... it runs also run so much slower than the original.

new code finds 1 core OTU

> rarefied_physeq
phyloseq-class experiment-level object
otu_table()   OTU Table:          [ 11975 taxa and 26 samples ]:
sample_data() Sample Data:        [ 26 samples by 7 sample variables ]:
tax_table()   Taxonomy Table:     [ 11975 taxa by 7 taxonomic ranks ]:
taxa are rows

>  sample_sums(rarefied_physeq)
  AQC1cm   AQC4cm   AQC7cm      CC1      CL3    Even1    Even2    Even3  F21Plmr LMEpi24M  M11Fcsw  M11Plmr  M11Tong  M31Fcsw  M31Plmr 
    5000     5000     5000     5000     5000     5000     5000     5000     5000     5000     5000     5000     5000     5000     5000 
 M31Tong      NP2      NP3      NP5 SLEpi20M      SV1  TRRsed1  TRRsed2  TRRsed3     TS28     TS29 
    5000     5000     5000     5000     5000     5000     5000     5000     5000     5000     5000 

>  # extract core
>  system.time(
+  GP_core_rare_inc<- parallel_extract_core(
+      rarefied_physeq,
+      Var = "SampleType",
+      method = "increase",
+      increase_value = 2,
+      .parallel = FALSE,
+      ncores = 10, 
+      set_seed=154
+  ) 
+  )
Seed used: 154
✔ Input phyloseq object is valid!
ℹ Dataset is rarefied at a depth of: 5000
ℹ Group: not used.
ℹ Core calculated across: SampleType
ℹ Calculating BC dissimilarity based on the 1st ranked OTU
✔ BC dissimilarity based on the 1st ranked OTU complete
ℹ Calculating BC dissimilarity based on ranked OTUs, starting at 2025-09-11 16:47:28.453326
ℹ Running in serial
✔ BC ranks done!                                 
ℹ Performing method 'increase'
! No OTUs meet the specified increase threshold. Using top OTU only.
✔ % increase method identified 1 core OTUs
✔ Analysis complete!
   user  system elapsed 
239.959   0.547 240.230 

>  system.time(
+  GP_core_rare_elb <- parallel_extract_core(
+      rarefied_physeq,
+      Var = "SampleType",
+      method = "elbow",
+      ncores = 10,
+      set_seed=312
+  )
+  )
Seed used: 312
✔ Input phyloseq object is valid!
ℹ Dataset is rarefied at a depth of: 5000
ℹ Group: not used.
ℹ Core calculated across: SampleType
ℹ Calculating BC dissimilarity based on the 1st ranked OTU
✔ BC dissimilarity based on the 1st ranked OTU complete
ℹ Calculating BC dissimilarity based on ranked OTUs, starting at 2025-09-11 16:56:18.06286
ℹ Running in serial
✔ BC ranks done!                                 
ℹ Performing method 'elbow'
[1] "549656"
✔ Elbow method identified 1 core OTUs
✔ Analysis complete!
   user  system elapsed 
227.777   0.637 227.984 

old code find

> length(test_origin_elbow[[1]])
[1] 10
> length(test_origin_increase[[1]])
[1] 48

and the code below...

> system.time(
+     test_origin_increase <- ExtractCore(physeq = rarefied_physeq,
+                                         Var = "SampleType",
+                                         method = "increase")
+ )
[1] "Dataset is rarefied at a depth of:"
[1] 5000
phyloseq-class experiment-level object
otu_table()   OTU Table:          [ 11975 taxa and 26 samples ]:
sample_data() Sample Data:        [ 26 samples by 7 sample variables ]:
tax_table()   Taxonomy Table:     [ 11975 taxa by 7 taxonomic ranks ]:
taxa are rows
[1] 11975     7
[1] "Dataset not subsetted"
[1] "Check: dimension of datframe and metadata"
[1] 11975    26
[1] 26  8
`summarise()` has grouped output by 'otu'. You can override using the `.groups` argument.
Joining with `by = join_by(rank)`
 [1] "549656" "317182" "94166"  "331820" "189047" "329744" "360229" "550960" "98605"  "279599" "534609" "158660" "298875" "171551"
[15] "322235" "471122" "193761" "108747" "252398" "317658" "64396"  "327536" "340216" "114612" "82772"  "228764" "361496" "248140"
[29] "319505" "12812"  "471141" "561077" "235567" "557211" "141709" "213853" "200741" "484436" "469778" "259569" "306180" "513590"
[43] "359171" "518438" "357795" "470973" "114821" "326977"
   user  system elapsed 
  4.914   0.031   4.964 
Warning message:
`add_rownames()` was deprecated in dplyr 1.0.0.
ℹ Please use `tibble::rownames_to_column()` instead.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated. 

> test_origin_increase[[1]]
 [1] "549656" "317182" "94166"  "331820" "189047" "329744" "360229" "550960" "98605"  "279599" "534609" "158660" "298875" "171551"
[15] "322235" "471122" "193761" "108747" "252398" "317658" "64396"  "327536" "340216" "114612" "82772"  "228764" "361496" "248140"
[29] "319505" "12812"  "471141" "561077" "235567" "557211" "141709" "213853" "200741" "484436" "469778" "259569" "306180" "513590"
[43] "359171" "518438" "357795" "470973" "114821" "326977"
> system.time(
+     test_origin_elbow <- ExtractCore(physeq = rarefied_physeq,
+                                      Var = "SampleType",
+                                      method = "elbow")
+ )
[1] "Dataset is rarefied at a depth of:"
[1] 5000
phyloseq-class experiment-level object
otu_table()   OTU Table:          [ 11975 taxa and 26 samples ]:
sample_data() Sample Data:        [ 26 samples by 7 sample variables ]:
tax_table()   Taxonomy Table:     [ 11975 taxa by 7 taxonomic ranks ]:
taxa are rows
[1] 11975     7
[1] "Dataset not subsetted"
[1] "Check: dimension of datframe and metadata"
[1] 11975    26
[1] 26  8
`summarise()` has grouped output by 'otu'. You can override using the `.groups` argument.
Joining with `by = join_by(rank)`
 [1] "549656" "317182" "94166"  "331820" "189047" "329744" "360229" "550960" "98605"  "279599"
   user  system elapsed 
  5.176   0.009   5.213 
Warning message:
`add_rownames()` was deprecated in dplyr 1.0.0.
ℹ Please use `tibble::rownames_to_column()` instead.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated. 

> test_origin_elbow[[1]]
 [1] "549656" "317182" "94166"  "331820" "189047" "329744" "360229" "550960" "98605"  "279599"

[jibarozzo]

They way we are now using 'pivot_longer', which.max, etc etc. how we are manipulating the DF inside, I believe could be that cause of this.

[jibarozzo]

We need to take a step back, and like you said go line by line.

[jibarozzo]

Scratch ChatGPt and go oldschool haha