Early Prediction of Anti-PD-1 Therapy Response in Hepatocellular Carcinoma Using Gut Microbiota Biomarkers

Early Prediction of Anti-PD-1 Therapy Response in Hepatocellular Carcinoma Using Gut Microbiota Biomarkers | medRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-P4HH5NV'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search Early Prediction of Anti-PD-1 Therapy Response in Hepatocellular Carcinoma Using Gut Microbiota Biomarkers Phing Chian Chai , View ORCID Profile Oleg V. Moskvin , Rohan Williams , Damien Keogh doi: https://doi.org/10.1101/2025.05.15.25327728 Phing Chian Chai 1 BluMaiden Biosciences , Singapore , 139965 Find this author on Google Scholar Find this author on PubMed Search for this author on this site Oleg V. Moskvin 1 BluMaiden Biosciences , Singapore , 139965 Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Oleg V. Moskvin For correspondence: oleg{at}blumaiden.com damien{at}blumaiden.com Rohan Williams 1 BluMaiden Biosciences , Singapore , 139965 Find this author on Google Scholar Find this author on PubMed Search for this author on this site Damien Keogh 1 BluMaiden Biosciences , Singapore , 139965 Find this author on Google Scholar Find this author on PubMed Search for this author on this site For correspondence: oleg{at}blumaiden.com damien{at}blumaiden.com Abstract Full Text Info/History Metrics Data/Code Preview PDF Abstract Hepatocellular carcinoma (HCC) remains a leading cause of cancer-related mortality, and response rates to anti-PD-1 therapy are suboptimal. Previous models predicting therapy response often fail to provide reliable early predictions, limiting their clinical utility. To develop a robust predictive model based on gut microbiota features that can accurately discriminate between responders (R) and non-responders (NR) to anti-PD-1 therapy at baseline, enhancing early clinical decision-making. We reanalyzed microbiome abundance data from a cohort of HCC patients receiving anti-PD-1 therapy, originally published by Wu et al. 1 . By incorporating microbial taxa from phylum to genus levels, we addressed taxonomic assignment uncertainty. Simultaneously, we leveraged low-abundance, informative features via data transformation. Using statistical and machine learning tools, we rectified the informative features, reducing them to 8 key taxa. A neural network model was developed and performance-assessed. Our predictive model, utilizing 8 key microbial taxa, demonstrated an AUROC of 0.784 ± 0.024 on test data. Also, unlike the negative result of the original publication, our model could also successfully discriminate responders from non-responders at baseline, providing a critical advantage at the early intervention stage. This study presents a clinically valuable predictive model that enhances the ability to forecast responses to anti-PD-1 therapy in HCC patients using gut microbiota features. The model’s simplicity and high predictive accuracy offer advantages in personalized treatment planning, contributing to the improvement of patient outcomes in HCC immunotherapy. Introduction Hepatocellular carcinoma (HCC) is a major global health concern, being the sixth most common cancer and the third leading cause of cancer-related deaths worldwide. The prognosis for advanced HCC remains poor, with limited effective treatment options historically available 2 . Immune checkpoint inhibitor (ICI) therapies have emerged as a promising treatment modality for HCC. ICIs, such as those targeting programmed cell death protein 1 (PD-1) and its ligand (PD-L1), have shown significant potential in reactivating the immune system to target and destroy cancer cells 3 . Despite the promise of ICIs, a significant proportion of HCC patients do not respond to these therapies, exhibiting either primary or acquired resistance 4 . This resistance poses a major challenge, as it limits the overall efficacy of ICIs and complicates treatment planning. Identifying patients who are likely to be resistant to ICIs before starting therapy is crucial for optimizing treatment strategies and improving outcomes 5 . Currently, host-derived biomarkers such as genomic variants, expressed proteins or enrichment of specific host cell types in the tumor microenvironment have been used to predict responses to ICIs. Examples include tumor mutational burden (TMB), microsatellite instability and genetic variants in a number of specific pathways that may be associated with either therapy response or absence of it 6 , the expression level of PD-L1 itself 7 and the abundance of CD8+ T cells 8 . However, these biomarkers often have either limited predictive power or/are invasive to obtain, requiring tissue biopsies. The role of microbiota in response to anti-PD-1 therapy was noted as early as in 2018, when higher abundances of Fecalibacterium and Bacteroidales were associated with longer progression-free survival 9 . Other studies have demonstrated that the diversity and specific composition of gut microbiota can significantly impact the effectiveness of ICIs by modulating host immune responses. Notably, enrichment of particular commensals such as Faecalibacterium or Akkermansia muciniphila , as well as overall microbial diversity, has been correlated with enhanced response to PD-1 blockade in several cancers 10 11 . Currently, the gut microbiome is being considered as a potential source of biomarkers for predicting responses to immune checkpoint inhibitors 12 The biological mechanisms of microbiota-host interactions that affect the immunotherapy responses are overviewed in 13 and 14 , providing foundational knowledge for future predictive model development. Moreover, stool sample collection and analysis to assess gut microbiota offer several practical benefits. Unlike tissue biopsies required for genomic biomarker analysis, stool sample collection is noninvasive, low-cost and readily available. Thus, the inclusion of microbiome-based assays for response prediction could add significant value in this clinical setting. Research has identified microbiome-based therapy response predictors for various cancers, underscoring their potential in oncology. For instance, a network-based microbiome feature survival predictor outperformed traditional host-related factors such as PD-L1 levels, antibiotic use, and age in predicting outcomes in non-small cell lung cancer (NSCLC) 15 . Similarly, in melanoma, metagenomics studies have highlighted the role of specific bacterial taxa in influencing therapy response 16 . These findings illustrate the broader applicability of microbiome biomarkers across cancer types. The strategic vision of the Society for Immunotherapy of Cancer (SITC) expressed last year 17 highlights the importance of pre-clinical models in advancing microbiome-based predictors, emphasizing the need for continued research and innovation in this area. Importantly, the existing research already points out the causal links between gut microbiota and cancer immunotherapy, beyond statistical associations. This conclusion may be derived on the basis of either Mendelian randomization studies that, for example, demonstrate a causal link between gut microbiota composition and the risk of biliary tract cancer 18 or an improvement of therapy response after direct supplementation of live bacterial supplementation 19 . In hepatocellular carcinoma, microbiome-based predictors are still in the research stage, but early findings are promising. While studies have noted the involvement of gut microbiota in anti-PD-1 therapy response in hepatobiliary cancers, including HCC 20 , the development of robust predictive models remains ongoing. The relationship between gut microbiota and HCC prognosis is complex, necessitating further investigation to validate potential biomarkers and refine prediction models. As a result, as of now, the field of HCC prognostic markers based on the microbiome continues to stay at the research stage 21 . The main reason for the absence of a clinically-ready model is the non-performance of the reported models at baseline. In this study, we reanalyze publicly available microbiome and clinical data from hepatocellular carcinoma (HCC) patients who were undergoing anti-PD-1-based immunotherapy 1 . Our re-analysis, based on custom feature engineering and machine learning methodology, demonstrates the feasibility of a robust, rigorously derived model that can predict the response to anti-PD-1 therapy based on samples obtained before the onset of the drug intervention, a finding that was not achieved in the original analysis. In this paper, we successfully addressed the known outstanding issues of predicting the anti-PD-1 therapy response from gut microbiome data that impede clinical applications: a) the granularity of the microbial features that are sparsely distributed between the individuals, b) the large number of features in the models, c) absence of the predictive signal at the baseline. This was achieved by a) considering a set of features across taxonomy ranks, b) incorporating features with lower abundance but high prevalence via data transformation and c) designing a custom machine learning architecture. Our analysis offers a substantive advance on the findings of the source study that only detected a response-related signature for subjects after commencement of the course of the drug treatment. In contrast, our custom modeling study resulted in a qualitative shift in the clinical utility of the result by enabling the response prediction before the subjects are exposed to any dose of the drug. This sets the stage for meaningful translation of microbiome research in this area to clinical practice, contributing to the improvement of the quality of life of cancer patients and using the critically limited time to search for alternative treatments in the case of predicted non-responsiveness. Methods Dataset We used the publicly available data from Wu et al. 1 , a study of 35 patients (31 males, 4 females, mean age of 53.8 years) with hepatocellular carcinoma (HCC) who undergo anti-PD-1-based immunotherapy. The majority of the patients were in the advanced stage of HCC, classified as Barcelona Clinic Liver Cancer stage C (32/35, 91.45%). The distinction between Responders (R, 46% of the cohort) and Non-responders (NR, 54%) was performed in the original publication on the basis of a combination of complete response (CR), a partial response (PR), or stable disease (SD) labels and the follow-up time 1 . This 46%/54% distribution shows a relatively balanced dataset, allowing for robust comparative analyses. Notably, the R and NR groups were also balanced by age, gender, liver function status, and underlying etiologies, with the exception of baseline alpha-fetoprotein (AFP) levels, immature granulocyte counts, and lymphocyte counts. The raw 16S rRNA gene sequencing data representing this study in the National Genomics Data Center (NGDC) Genome Sequence Archive (GSA) database (accession: CRA005195) is accompanied by a limited metadata file with only two columns containing the patients’ data: “Public description” with R / NR status, and “Collection date”, apart from other non-informative in the given context or sequencing-specific data. Applying the knowledge of the three-week interval between fecal sample collections, we were able to match the samples to all 35 patients (Supp Data 1). The number of anti-PD-1 injections and corresponding fecal samples varied among patients, as illustrated in Figure S1 . Smaller number of samples belong to subjects dropped out of the study for various reasons. Download figure Open in new tab Figure S1: Number of Anti-PD-1 doses received by each patient. Barplot showing the number of anti-PD-1 doses received by each patient. Pink bars represent responders (R), and blue bars represent non-responder (NR). The plot highlights the variability in the number of doses administered to each patient. The identifiers of the patients on the x-axis are arbitrary and based on the order of samples given in the metadata. Microbiome data analysis Raw sequencing data (FASTQ files) representing the above dataset were downloaded from GSA and re-processed via BluMaiden standard 16S amplicon analysis workflow, using the year 2023 update of the taxonomy databases. Low-level preprocessing of the data followed the BluMaiden KEYSTONE™ pipeline; modifications relevant to the optimization of the predictive model’s performance are described in the Results. Data analysis was performed using a combination of R statistical environment 22 and Python programming language, using a combination of custom scripting with open-source packages. Calculation of community diversity measures and PERMANOVA analysis were performed using the scikit-bio package (ver. 0.6.0). Differential abundance of microbiota features relevant to the clinical phenotypes was assessed in a multivariate analysis setting using Maaslin2 23 version 1.14.1. The analysis was run with min_abundance = 0 and min_prevalence = 0.1 parameters, setting factors of interest (Responder Status and Doses Taken) as target variables. Due to the use of a custom normalization approach, the internal standardization of the algorithm was switched off by setting standardize = FALSE . Predictive modeling Custom ML model construction was performed using Python programming language, open-source packages and in-house developments of BluMaiden Biosciences. Feature importance assessment was performed using Shapley Additive exPlanations (SHAP) values. We conducted 300 bootstrap trials with optimized hyperparameters and selected the top quartile of trials where features consistently exhibited high absolute SHAP values in the validation set. The model was evaluated on the test dataset using AUROC scores to measure predictive accuracy. To visualize the prediction output on individual samples, the binary classification of the responder status was generated based on a probability score (below or above 0.5) mapped by a sigmoid function. The biomarkers identified as a result of predictive modeling are claimed in the United States (US) Patent Application No: 63/762204. Results and Discussion Microbiome compositional differences in clinical groups Conventional characterization of gut microbiome compositional differences between clinical groups involves both ecological measures (species diversity, richness and evenness) and associations of the abundance levels of the microbial features with clinical parameters of interest via analysis of variance and mixed linear models. The original paper applied those methods and confirmed their insufficiency in deriving clinically applicable biomarkers, as opposed to non-linear machine learning methods that demonstrated practical potential. Since our findings confirm this observation (and also identified flaws and unused potential of the original study, making our predictive model - see section 4.2 - the main outcome of this re-analysis, we only briefly discuss the traditional descriptive approaches to a) establish a common baseline with the published study, b) explore deeper taxonomy representation that includes multiple taxonomy ranks to balance generalizability and precision and c) leverage the updated status of taxonomy interpretation of the same sequencing data. Figure N1A shows the results of ecological parameter analysis, which included Shannon diversity, Chao1 species richness and Simpson’s evenness. Apparently, responders and non-responders exhibited similar levels of baseline microbial diversity when assessed at the genus level. While the Shannon Diversity index indicated that responders had slightly higher diversity, none of the differences in all ecological measures between responders and non-responders were statistically significant (Mann Whitney-U test - Shannon, p =0.253; Chao1, p =0.619; Simpson’s, p =0.456), consistent with the findings of the original study performed at OTU level. The non-informativeness of the overall ecological measures of the microbial communities in predicting the therapy response status was further supported by the lack of differences in Shannon diversity between the two groups after a more detailed look that included: a) consideration of 5 taxonomy ranks (from phylum to genus) and b) or separate assessment of both the entire cohort of samples and a subset of samples describing the baseline microbiome state, before receiving the first dose ( Figure NS2 ). The latter subset is the most clinically interesting, since any signal of intervention response identified there, before the actual drug was ever applied, would provide the earliest and the most valuable patient stratification insights. Download figure Open in new tab Figure S2: Shannon diversity analysis of microbial communities. Boxplots showing Shannon diversity indices for microbial communities at different taxonomic levels (phylum, class, order, family, genus) for responder (R) and non-responder (NR). The upper row represents samples taken at baseline, while the lower row includes all samples. Pink boxes represent responders, and blue boxes represent non-responders. No significant differences in Shannon diversity were observed between R and NR at any taxonomic level. Download figure Open in new tab Figure 1: Baseline microbial diversity and composition analysis. A) Alpha diversity metrics comparing responders (R) and non-responders (NR) at baseline. The boxplots display three different measures of alpha diversity: Shannon Diversity, Chao1 Species Richness, and Simpson’s Evenness; B) Principal Coordinate Analysis (PCoA) plot of microbial community composition at baseline, using Bray-Curtis dissimilarity. Each dot represents a sample, with blue dots indicating NR and pink dots indicating R; C) Stacked bar plot showing the proportions of microbial genera in each patient at baseline. The top eight genera are displayed, with all other genera aggregated into the category ‘Others’. The color-coded bars represent the relative abundance of each genus, highlighting the overall microbial composition differences between the patients. Beta diversity was next assessed to probe the differences in microbial community composition between responders and non-responders. The principal coordinate analysis (PCoA) plot based on Bray-Curtis dissimilarity ( Figure N1B ) indicated no significant differences in microbial communities between the two groups at baseline (PERMANOVA: F=0.988, p=0.391). The first principal coordinate (PCo1) explains 20.8% of the variance, while the second principal coordinate (PCo2) explains 15.3%. Other dissimilarity metrics, including Jaccard distance and weighted Unifrac, were also unhelpful in revealing the difference at baseline ( Figure NS3 ). Also, compositional representations of neither the top genera ( Figure N1C ) nor the other 4 taxonomy ranks ( Figure NS4 ) were able to reveal the microbiome compositional patterns of association with the therapy response. Interestingly, the initial study 1 reported differences prior to the third dose of anti-PD-1 injection (including baseline) using OTU-based beta diversity measures. This discrepancy may be attributed to the different taxonomic resolution used in our analysis (genus level, as opposed to OTU). The latter representation, while expected to capture finer-scale variations in both signal and noise in the microbial composition, still was not helpful to derive a robust predictive model in the original study, presumably, because of the limited generalizability / higher sparseness of the OTU features. Taken together, both alpha and beta diversity measures at baseline failed to meaningfully stratify the patients into responders and non-responders. Download figure Open in new tab Figure S3: Principal Coordinate Analysis (PCoA) of microbial community composition. Principal Coordinate Analysis (PCoA) plots depicting beta diversity measures using Bray-Curtis dissimilarity, Jaccard distance, and weighted UniFrac for responders (R) and non-responders (NR). Pink dots represent R and blue dots represent NR. The upper row shows the PCoA plots at baseline. The principal coordinates, PCo1 and PCo2 on the x- and y-axis, respectively, represent the top two axes of variance in the microbial community composition. The lower row shows the PCoA plots for all samples. Download figure Open in new tab Figure S4: Baseline microbiota composition of the study subjects Stacked bar plots showing the proportions of microbiota at four taxonomic levels - phylum, class, order, and family - before anti-PD-1 intervention in non-responders (NR) and responders (R). For all plots, each bar represents an individual patient, with colors indicating different taxa. Only the top 8 taxa by overall abundance are annotated, with aggregation of the lower-abundant taxa as “Others”. Similar to the observation in the original publication on the accumulation of the microbiome differences between R and NR groups following the course of intervention (that allowed the authors to build the R/NR classification (not prediction, de facto) model based on the data from subjects deep into the treatment), we also could detect an association between the responder status and microbiome compositions of the drug-exposed individuals. In particular, significant differences in microbial community composition between R and NR groups were detected in the entire pool of samples using either Bray-Curtis dissimilarity (PERMANOVA: F=2.793, p=0.004) or the weighted Unifrac distance (PERMANOVA: F=4.251, p=0.004). These results indicate that the microbial communities of responders and non-responders diverge more noticeably when the entire course of treatment is considered ( Figure NS3 ). Besides group-level differences, we checked if a more granular approach of associating individual microbial features with therapy response can reveal a statistically acceptable signal. For this purpose, we applied generalized linear modeling implemented in the MaAsLin2 package. This analysis revealed that none of the genera were significantly associated with either the responder status of the patients or the number of anti-PD-1 doses taken ( Figure N2 ). The absence of significant associations in this multivariate analysis (q: 0.428-0.998 (responder status), 0.456-0.962 (doses taken)) adds another level of microbiome compositional representation (individual microbial features analyzed in a conventional way) that, similar to more generalized parameters like alpha and beta diversity is insufficient to detect a signal of response to anti-PD-1 therapy. This highlights the complexity of the fecal microbiome and its interaction with clinical outcomes, suggesting that more sophisticated models may be required to uncover meaningful biomarkers for patient stratification, especially in the absence of other useful clinical metadata. We addressed the challenge of building such a non-conventional well-performing model below. Download figure Open in new tab Figure 2: Phylogenetic tree with a heatmap of multivariate analysis. The phylogenetic tree displays 100 genera, filtered to remove low-abundance taxa, with branches and terminal leaves colored according to their respective phyla. The outer circular heatmap shows q-values from the multivariate MaAsLin2 analysis, which assessed associations between “Responder” status, the number of anti-PD-1 “Doses”, and the abundance of each genus. Genera are arranged based on their phylogenetic relationships, with q-values indicating the strength of association (darker shades represent lower q-values, which are more significant). Microbiome-based predictor of treatment response As mentioned in Methods, the set of study subjects had a favorable characteristic of being balanced by the majority of clinical parameters. This balance sets the stage to discover robust clinical biomarkers of the therapy response that won’t be biased by other clinical differences. If the challenge of invisibility of the association signal (item 4.1) can be successfully resolved using less conventional approaches, the resulting markers are expected to have the potential to be unbiased and applicable for patient stratification. To reveal the clinically applicable signal, we addressed the following challenges: 1) conservative partitioning of the subjects into train and test sets while preventing any data leakage, 2) overcoming high granularity of taxonomic features when observed across subjects 3) finding the right balance between feature’s within-sample abundance and across-sample frequency of the detection and 4) optimizing the learning architecture. Partitioning Samples into Train and Test Sets For building an ML-based predictive model, partitioning the patients into suitable train and test sets was critical, given the high variance in the number of anti-PD-1 doses taken by each patient and the respective variance in the number of fecal samples per person. In the original article by Wu et al., the authors do not specify how they conducted their train-test split beyond mentioning a 70:30 train-validation split, presumably done randomly. This approach poses a risk of data leakage, as samples from the same patients at different time points could end up in both train and test sets. Since fecal microbiome profiles tend to be more invariant within subjects than between subjects (PERMANOVA, p=0.001), it is important to place all samples from a single patient in either train or test set exclusively. To adequately represent both responders and non-responders in the training and test sets while preventing data leakage, we used heuristic optimization (see Methods) to split the samples into 109 samples from 24 patients in the training set, and 46 samples from a distinct group of 11 patients in the test set. The resulting 70:30 train-test split showed the balanced distribution by both response, overall number of samples and other metrics ( Figure N3A-B ). This robust train-test split mitigates the risk of data leakage and maintains a balanced representation, ensuring the overall effectiveness and real-world applicability of the generated model. Download figure Open in new tab Figure 3: Sample and subject counts from heuristic-optimized train-test dataset. A) Pie charts displaying the distribution of samples and subject counts in the train and test sets for the total cohort, responders (R) and non-responders (NR). Blue sections represent the counts in the train set, while pink sections represent the counts in the test set; B) Box plots showing the highest number of anti-PD-1 doses recorded for R and NR in both the train (blue) and test (pink) sets. The distribution illustrates the balance achieved in the train-test split, with both sets maintaining similar dose distribution profiles for R and NR. Even though Wu et al. achieved a high AUROC performance score of 0.837 on validation samples for their microbiome-based random forest model, we believe that might have been artificially inflated due to data leakage. To establish a more realistic benchmark for model performance, we leveraged various established auto-ML libraries to serve as benchmarks. As expected, AutoGluon, TPOT, auto-sklearn, and mljar-supervised only produced the best models with AUROCs between 0.479 and 0.696 on our leak-free test set with 408 features ( Figure N4A , see next section). Improving this out-of-the-box performance required custom decisions on both the model architecture and the microbiome data representation per se . Download figure Open in new tab Figure 4: AutoML library performance with heuristic optimized train-test dataset and workflow for modelling from 16S Fastq data. A) Bar plot showing the AUROC performance of the best models built using four different AutoML libraries: auto-sklearn, AutoGluon, mljar-supervised, and TPOT. The validation and test AUROC scores are presented, illustrating the benchmark performance established by these libraries. The best models were selected based on the highest validation scores achieved; B) Conceptual overview of the BM pipelines for 16S rRNA data processing. The bioinformatics pipeline handles initial data processing steps, while the modelling pipeline encompasses data transformation, feature selection, hyperparameter tuning and model evaluation. Predictive modeling To develop robust predictive models, we carefully processed the microbiome abundance data to extract meaningful features while addressing common challenges in fecal microbiome analysis. First, we utilised microbial taxa from phylum to genus ranks as features. Although this method introduces overlapping features, it is necessary due to the inherent uncertainty in taxonomic assignment and partial mutual support of less informative (but more stable) features of higher taxonomic ranks and potentially more informative (but sparse across individuals) features of lower taxonomy ranks. Moreover, due to resolution limitations of the 16S metabarcoding data, sequencing reads often remain unclassified at lower taxonomic ranks (for example, at genus level) while being represented by aggregated counts at higher ranks. By incorporating features from multiple taxonomic levels, we captured a more complete picture of the microbial community structure with reduced sparsity issues. This approach resulted in a pool of 408 taxonomic features: 14 from phylum, 21 from class, 55 from order, 92 from family, and 226 from genus ( Figure N4B ). Because of the notorious interpersonal variability of microbiome profiles, each individual in the training set had, on average, only 25.5 ± 8.0% of the total features represented. This sparsity is a common issue in microbiome data and can lead to significant data loss if not managed appropriately. A popular practice of filtering out the low-abundance features risks discarding taxa that might still provide valuable information. To leverage this opportunity, we applied sample-wise log transformation to stabilize variance, normalize the distribution of feature counts, mitigate the influence of extreme values and ensure the retention of low-abundance but potentially informative taxa ( Figure N4B ). Despite these preprocessing steps, many microbial features may still be affected by noise from sequencing or unavoidably uncertain taxonomic assignments. To address this, we employed our proprietary suite of statistical and machine learning tools to identify and rank features based on their association with and informativeness on the predictive outcome, specifically the R and NR status in the training set. Together, these techniques reduced the number of features from 408 to 8, which constituted the core microbial taxa that could serve as biomarkers for detecting responders vs. non-responders to anti-PD-1 therapy among HCC patients. Our final model was a neural network utilizing these 8 key taxa as features. After hyperparameter tuning, the performance of this model, as illustrated in the ROC curve ( Figure N5A ), achieved an AUROC of 0.784 ± 0.024 on test data. This result indicates a good predictive capability, significantly outperforming all auto-ML benchmarks (auto-sklearn: 0.479, AutoGluon: 0.604, mljar-supervised: 0.696, TPOT: 0.689). As opposed to Wu et al.’s model, which utilized 18 biomarkers, our model’s use of 8 highly informative taxa offers several practical benefits. Additionally, our leakage-free approach ensures that the model’s performance is a true reflection of its predictive power, free from biases introduced by data leakage, which can artificially inflate performance metrics. Our simpler model with less than half of biomarkers reduces the complexity and cost of clinical testing while enhancing interpretability. Even though our test AUROC appeared to lag Wu et al.’s internal validation score of 0.837, we could achieve a validation AUROC of 0.915 ± 0.059 while closely reproducing the original paper’s analysis, employing their non-conservative choices (data not shown). Download figure Open in new tab Figure 5: Neural network model performance and responder status predictions. A) Receiver operating characteristic (ROC) curve showing the performance of the 8-taxa neural network model. The blue line represents the mean performance across 100 trials, while the cyan shaded area shows the standard deviation. The model achieved an AUROC of 0.784 ± 0.024 on the test data, significantly outperforming random chance; B) Violin plots illustrating the distribution of AUROC scores across 100 models for different sample groups: all samples, samples before intervention (baseline), samples before the third dose of anti-PD-1, and samples from the third dose onwards. The horizontal lines within the violins indicate the median performance; C) Heatmap of the predicted R/NR status based on the probability score provided by the model with a threshold of 0.5. The probabilities are the mean of 100 trials. The heatmap includes 155 samples from all 35 patients in both the train and test sets. The bottom horizontal heatmap shows the ground truth responder status based on clinical evaluation, with pink denoting responders and blue denoting non-responders. The patient identifiers are listed along the x-axis, and the number of doses is represented on the y-axis. The qualitative advantage of our model over the one from the original publication is the ability to detect the response signal before starting the drug intervention. The model by Wu et al. was not able to see this signal at all (random chance AUROC=0.5 level of prediction) despite their looser standards of data handling, especially the absence of control for data leakage. The ability to detect the therapy response signal in the microbiome data at baseline provides a tremendous advantage in clinical utility of such models, especially in combination with the parsimonious feature composition allowing for simpler and cheaper implementation of the necessary test. Insights into the informative taxa and clinical utility projection Our 8 biomarkers included 3 taxonomy features at the family level ( Acidaminococcaceae , Pasteurellaceae and Tannerellaceae ) and 5 - at the genus level ( Ligilactobacillus , Dorea , Barnesiella , Ruminococcus and Atopobium ). In contrast, Wu et al.’s model exclusively utilized biomarkers from the genus level. All three of our family-level biomarkers have corresponding genus-level biomarkers (one for each family) used in Wu et al.’s model ( Figure N6A ). These family-level markers have the potential to capture a broader range of bacterial species with similar biological properties and predictive values, overcoming the granularity of genus-level features across individuals. Notably, the most informative marker based on SHAP (Shapley Additive exPlanations) values in our study, Acidaminococcaceae , maps to two taxa at the genus level, Acidaminococcus and Phascolarctobacterium , with only the latter identified as a key microbial biomarker (though the least important) in Wu et al.’s study. This indicates that our all-taxonomic-rank approach could take advantage of features encompassing a broader microbial diversity. Download figure Open in new tab Figure 6: Biomarkers importance derived from SHAP analysis. A) Comparison of the 8 biomarkers identified in this study with the 18 biomarkers from Wu et al. Biomarkers are arranged in decreasing importance, with exact matches, cross-rank matches, and high phylogenetic closeness (with a possibility of the de facto exact match given the phylogenetic uncertainty and evolving classification) indicated. Our 8 biomarkers include Acidaminococcaceae, Ligilactobacillus, Pasteurellaceae, Dorea, Barnesiella, Tannerellaceae, Ruminococcus , and Atopobium , highlighting overlapping and distinct taxa compared to Wu et al.; B) SHAP analysis beeswarm plot illustrating the impact of the 8 biomarkers on the model’s output. Each point represents a SHAP value for a single prediction, with color indicating biomarker abundance after transformation and min-max normalization (purple for high abundance and blue for low abundance). Higher positive SHAP values indicate a stronger contribution towards predicting responder status, while higher negative SHAP values indicate a stronger contribution towards predicting non-responder status. The points are an aggregation of predictions from the validation samples, generated from the top 37 models out of 150 models, independently trained with heuristic optimized bootstrapped train dataset. For our genus-level biomarkers, three ( Dorea , Barnesiella and Ruminococcus ) overlapped with Wu et al.’s key microbial biomarkers, highlighting their significance and providing convergence in the identification of crucial biomarkers. Interestingly, our method identified Ligilactobacillus as the second most informative biomarker, whereas Wu et al.’s model identified Lactobacillus as the most important biomarker. Ligilactobacillus genus is taxonomically close to Lactobacillus and was historically reclassified from the latter. Similar to Lactobacillus , it also has a substantial track record of probiotic properties 24 . Taxonomic closeness and the respective functional overlap between the two genera suggest that both models could be capturing the same underlying biological signal. The SHAP value plot provides insights into the contribution and impact of each biomarker on our model’s output ( Figure N6B ). The top contributors to the model’s predictions, Acidaminococcaceae and Ligilactobacillus , are shown to be negatively associated with responder status. Tannerellaceae and Atopobium also show a negative association with responders while biomarkers, as evidenced by their negative SHAP values. Conversely, Pasteurellaceae, Dorea, Barnesiella and Ruminococcus positively contributed to the model’s responder status. The detection of Acidaminococcaceae as a family-level biomarker of non-response makes sense if we recall that one genus of this family - Phascolarctobacterium – was captured in the same context by Wu et al., while another genus of this family – Acidaminococcus – was reported to be associated with adverse effects of immune checkpoint inhibitor therapy, being enriched by nearly 200 times in the Adverse Events group 25 . The underlying functional nature of the similar direction of the influence of different genera of this family and the ability to generalize the biomarker to the family taxonomy cannot be easily tracked in the literature and requires a focused investigation. Despite the fact that selection of this family-level marker by our procedure finds an instant echo in the above-mentioned relevant phenotypic associations involving the underlying genera, we prefer to avoid the temptation of over-analyzing the composition of our biomarker signature on the basis of the existing literature sentiment on health- or disease-associated role of individual taxonomy features. With the nature of the neural network model that highly relies on the complex non-linear relationships between the marker features abundance and clinical phenotypes, involvement of network-level relationships, other than “healthy” versus “disease-associated” traits of individual features is reasonably expected. This makes functional interpretation of the signature taxa a longer path that requires specific focused studies. The practical advantage of the unveiled biomarker combination, despite its unclear functional background, comes from solving two outstanding problems: 1) overcoming biomarker sparseness and 2) enabling early detection of the therapy response signal, before the first doze of the drug. Indeed, the entire set of our 8 biomarkers was found in at least 12.3% of the patients, with individual marker detection showing the median prevalence of 67.8% across all patients. In comparison, the median prevalence of biomarkers from Wu et al.’s model was substantially lower, at 51%, with some of the markers (like Aggregatibacter and Parasutterella ) dropping all the way to 1.94% of the patients, which is over 6 times less than the minimal prevalence observed for our model’s markers. This sparseness of the markers of the original model is expected to decrease the consistency of predictive power across different patient cohorts. The robustness of our model was additionally confirmed by its performance across different stages of anti-PD-1 intervention: it showed a useful degree of predictive power even at early stages ( Figure N5B ), such as the baseline (AUROC: 0.681 ± 0.056) and before the third dose of anti-PD-1 (AUROC: 0.688 ± 0.058). This qualitatively differs from the original model that performed randomly with an AUROC of around 0.5 at any time before the third dose. The ability to stratify patients at these early stages, especially before the very first dose of the drug, is invaluable for clinicians to make timely and informed decisions about the course of treatment. This result suggests that our model is able to detect the response signal that is not the result of drug-induced microbiome perturbations but rather belongs to the persistent individual microbiome signature. To provide additional evidence in support of detecting individualized response signal in microbiomes by our predictor, independently of the drug-induced microbiome-restructuring effects (captured by other studies, including Wu et al.), we provide per-sample visualization of the prediction model’s calls (probabilities of therapy response), with coloring split at the probability 0.5 ( Figure N5C ), i.e. the lower values (shades of teal) signify samples suggested the absence of therapy response, while higher values (shades of pink) mark samples giving the “responder” predictive signature of various confidence. Those per-sample prediction visualizations are stacked by subject on top of the ground truth response indication to allow analysis of per-subject dynamics of the response prediction signals. First, it is easily noticeable that in the absolute majority of the cases, the sign of the response prediction (color of the stacked squares), despite changes in the confidence level, stayed the same, with one prominent exception of subject #6 of the test set. This suggest that our model captures the properties of the individual microbiomes per se, as opposed to respond to the drug-induced microbiota perturbations as seen in the original publication 1 . At the same time, our model still captures the latter as a part of the overall signal. This is visible in the evolution of the confidence of the model calls following the increasing number of the drug doses, as seen in patients 24 and 31 (train set), as well as 6 and 13 (test set). This phenomenon contributes to overall increase in the classification performance after the 3 rd doze in our case as well ( Figure N5B ), resulting in correct prediction of the responder status for the majority of the patients based on their last fecal sample (87.5% and 72.7% on train and test sets, respectively). The above-mentioned case of subject #6 shows a rare instance when the initial microbiome composition of a non-responder subject was of the “responder” type as assessed by the model and this phenomenon even exacerbated after the first dose; however, subsequent doses resulted in the evolution of subsequent treatment-induced microbiome composition to converge to the correct – responder – signature. While it is early to speculate about the possibility of changing the responder status during the drug treatment course, given the limited data and restriction of its nature to the predictive model’s calls (that have their unavoidable degrees of uncertainty), such phenomenon, suggested by a remarkable switch of the prediction direction for the subject #6, may be theoretically considered and is worth investigating by plugging in clinical and metabolomic data in an integrative microbiota-host analysis. The latter may shed light on the nature of the possible dynamic therapy response alterations in the light of, for example, three distinct known HCC metabolic subtypes 26 that may interact with gut microbiota in specific ways and result in either stable or evolving response to the drug intervention. Despite the exception above, we should underline again the impressive overall stability of the R/NR binary calls by the model for individual subjects, capturing the individualized microbiome signatures as the main part of the detected signal. To add an intuitive numerical insight to this statement, please note that the subjects who received the longest treatment followed up by microbiome samples (14 dozes in the training set’s subject #14 or 12 doses in the test set’s subject #25) had, as the absolute majority of the subjects represented by multiple samples, single-direction response prediction in all the samples that belong to a subject, despite having a potential of 2 14 =16,384 or 2 12 =4,096 combinations of the binary classes if the classification was performed randomly. Conclusion We have shown that the issue of large inter-subject variation in the gut microbiome composition that usually limits the possibilities of creating a performant drug response predictor (that should balance informativeness and prevalence of individual features) can be overcome by selecting a representative mix of features across taxonomy ranks and processing them as inputs to a custom-optimized neural network. The resulting model is advantageous compared to other models reported in the area because of a combination of 1) its capability to stratify the subjects before the drug intervention, relying largely on the innate composition of the individual microbiome rather than drug-induced microbiome perturbations and 2) doing so in a parsimonious manner, without requiring a long list of biomarkers or expensive strain-level data resolution. This combination of traits opens the doors to implementing the results of our original data processing decisions and the resulting model as a simple and affordable clinical test for subject stratification before planning the application of the ICP inhibitor treatment decision. Data Availability The data that support the findings of this study (deposited by the authors of the original study 1 ) are openly available in NCBI Sequence Read Archive, BioProject PRJCA006814. Code Availability The BluMaiden KEYSTONE™ pipeline is described in Singapore Patent Application 10202500185R. Author contributions PCC has designed and executed analysis, interpreted results and drafted the manuscript. OVM designed and supervised the study, interpreted the results and wrote the manuscript. RW advised on study design and microbiome analysis procedures and edited the manuscript. DK designed and supervised the project, interpreted the results and edited the manuscript. Competing Interests All the authors are either current or former employees or shareholders of BluMaiden Biosciences. The authors declare no additional competing interests. References ↵ Wu , H. et al. Dynamic microbiome and metabolome analyses reveal the interaction between gut microbiota and anti-PD-1 based immunotherapy in hepatocellular carcinoma . International journal of cancer. Journal international du cancer 151 , 1321 – 1334 ( 2022 ). doi: 10.1002/ijc.34118 OpenUrl CrossRef PubMed ↵ Bicer , F. , Kure , C. , Ozluk , A. A. , El-Rayes , B. F. & Akce , M . Advances in Immunotherapy for Hepatocellular Carcinoma (HCC) . Curr Oncol 30 , 9789 – 9812 ( 2023 ). doi: 10.3390/curroncol30110711 OpenUrl CrossRef PubMed ↵ Ruli , T. M. , Jr. , et al. Immune Checkpoint Inhibitors in Hepatocellular Carcinoma and Their Hepatic-Related Side Effects: A Review . Cancers 16 ( 2024 ). doi: 10.3390/cancers16112042 OpenUrl CrossRef ↵ De Lorenzo , S. , Tovoli , F. & Trevisani , F . Mechanisms of Primary and Acquired Resistance to Immune Checkpoint Inhibitors in Patients with Hepatocellular Carcinoma . Cancers 14 ( 2022 ). doi: 10.3390/cancers14194616 OpenUrl CrossRef ↵ Kou , L. et al. The progress of research on immune checkpoint inhibitor resistance and reversal strategies for hepatocellular carcinoma . Cancer Immunol Immunother 72 , 3953 – 3969 ( 2023 ). doi: 10.1007/s00262-023-03568-3 OpenUrl CrossRef PubMed ↵ Bai , R. , Lv , Z. , Xu , D. & Cui , J . Predictive biomarkers for cancer immunotherapy with immune checkpoint inhibitors . Biomark Res 8 , 34 ( 2020 ). doi: 10.1186/s40364-020-00209-0 OpenUrl CrossRef PubMed ↵ Hao , L. et al. The current status and future of PD-L1 in liver cancer . Frontiers in immunology 14 , 1323581 ( 2023 ). doi: 10.3389/fimmu.2023.1323581 OpenUrl CrossRef PubMed ↵ Tumeh , P. C. et al. PD-1 blockade induces responses by inhibiting adaptive immune resistance . Nature 515 , 568 – 571 ( 2014 ). doi: 10.1038/nature13954 OpenUrl CrossRef PubMed Web of Science ↵ Gopalakrishnan , V. et al. Gut microbiome modulates response to anti-PD-1 immunotherapy in melanoma patients . Science 359 , 97 – 103 ( 2018 ). doi: 10.1126/science.aan4236 OpenUrl Abstract / FREE Full Text ↵ Routy , B. et al. Gut microbiome influences efficacy of PD-1-based immunotherapy against epithelial tumors . Science 359 , 91 – 97 ( 2018 ). doi: 10.1126/science.aan3706 OpenUrl Abstract / FREE Full Text ↵ Matson , V. et al. The commensal microbiome is associated with anti-PD-1 efficacy in metastatic melanoma patients . Science 359 , 104 – 108 ( 2018 ). doi: 10.1126/science.aao3290 OpenUrl Abstract / FREE Full Text ↵ Zhang , M. , Liu , J. & Xia , Q . Role of gut microbiome in cancer immunotherapy: from predictive biomarker to therapeutic target . Exp Hematol Oncol 12 , 84 ( 2023 ). doi: 10.1186/s40164-023-00442-x OpenUrl CrossRef ↵ Li , Z. et al. Critical role of the gut microbiota in immune responses and cancer immunotherapy . J Hematol Oncol 17 , 33 ( 2024 ). doi: 10.1186/s13045-024-01541-w OpenUrl CrossRef PubMed ↵ Gao , X. & Jiang , J . Exploring the regulatory mechanism of intestinal flora based on PD-1 receptor/ligand targeted cancer immunotherapy . Frontiers in immunology 15 , 1359029 ( 2024 ). doi: 10.3389/fimmu.2024.1359029 OpenUrl CrossRef PubMed ↵ Derosa , L. et al. Custom scoring based on ecological topology of gut microbiota associated with cancer immunotherapy outcome . Cell 187 , 3373 – 3389 e3316 ( 2024 ). doi: 10.1016/j.cell.2024.05.029 OpenUrl CrossRef PubMed ↵ Lee , K. A. et al. Cross-cohort gut microbiome associations with immune checkpoint inhibitor response in advanced melanoma . Nature medicine 28 , 535 – 544 ( 2022 ). doi: 10.1038/s41591-022-01695-5 OpenUrl CrossRef PubMed ↵ Emens , L. A. et al. Challenges and opportunities in cancer immunotherapy: a Society for Immunotherapy of Cancer (SITC) strategic vision . J Immunother Cancer 12 ( 2024 ). doi: 10.1136/jitc-2024-009063 OpenUrl Abstract / FREE Full Text ↵ Wang , K. et al. The causal relationship between gut microbiota and biliary tract cancer: comprehensive bidirectional Mendelian randomization analysis . Frontiers in cellular and infection microbiology 14 , 1308742 ( 2024 ). doi: 10.3389/fcimb.2024.1308742 OpenUrl CrossRef ↵ Ebrahimi , H. et al. Cabozantinib and nivolumab with or without live bacterial supplementation in metastatic renal cell carcinoma: a randomized phase 1 trial . Nature medicine ( 2024 ). doi: 10.1038/s41591-024-03086-4 OpenUrl CrossRef ↵ Mao , J. et al. Gut microbiome is associated with the clinical response to anti-PD-1 based immunotherapy in hepatobiliary cancers . J Immunother Cancer 9 ( 2021 ). doi: 10.1136/jitc-2021-003334 OpenUrl Abstract / FREE Full Text ↵ Zhang , N. et al. Biomarkers and prognostic factors of PD-1/PD-L1 inhibitor-based therapy in patients with advanced hepatocellular carcinoma . Biomark Res 12 , 26 ( 2024 ). doi: 10.1186/s40364-023-00535-z OpenUrl CrossRef PubMed ↵ R Team: R, A Language and Environment for Statistical Computing . ( 2015 ). ↵ Morgan , X. C. et al. Dysfunction of the intestinal microbiome in inflammatory bowel disease and treatment . Genome biology 13 , R79 ( 2012 ). doi: 10.1186/gb-2012-13-9-r79 OpenUrl CrossRef PubMed ↵ Yang , Y. et al. Understanding Ligilactobacillus salivarius from Probiotic Properties to Omics Technology: A Review . Foods 13 ( 2024 ). doi: 10.3390/foods13060895 OpenUrl CrossRef PubMed ↵ Hamada , K. et al. Turicibacter and Acidaminococcus predict immune-related adverse events and efficacy of immune checkpoint inhibitor . Frontiers in immunology 14 , 1164724 ( 2023 ). doi: 10.3389/fimmu.2023.1164724 OpenUrl CrossRef PubMed ↵ Bidkhori , G. et al. Metabolic network-based stratification of hepatocellular carcinoma reveals three distinct tumor subtypes . Proceedings of the National Academy of Sciences of the United States of America 115 , E11874 – E11883 ( 2018 ). doi: 10.1073/pnas.1807305115 OpenUrl Abstract / FREE Full Text View the discussion thread. Back to top Previous Next Posted May 16, 2025. Download PDF Data/Code Email Thank you for your interest in spreading the word about medRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Early Prediction of Anti-PD-1 Therapy Response in Hepatocellular Carcinoma Using Gut Microbiota Biomarkers Message Subject (Your Name) has forwarded a page to you from medRxiv Message Body (Your Name) thought you would like to see this page from the medRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Early Prediction of Anti-PD-1 Therapy Response in Hepatocellular Carcinoma Using Gut Microbiota Biomarkers Phing Chian Chai , Oleg V. Moskvin , Rohan Williams , Damien Keogh medRxiv 2025.05.15.25327728; doi: https://doi.org/10.1101/2025.05.15.25327728 Share This Article: Copy Citation Tools Early Prediction of Anti-PD-1 Therapy Response in Hepatocellular Carcinoma Using Gut Microbiota Biomarkers Phing Chian Chai , Oleg V. Moskvin , Rohan Williams , Damien Keogh medRxiv 2025.05.15.25327728; doi: https://doi.org/10.1101/2025.05.15.25327728 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Health Informatics Subject Areas All Articles Addiction Medicine (568) Allergy and Immunology (863) Anesthesia (299) Cardiovascular Medicine (4425) Dentistry and Oral Medicine (443) Dermatology (382) Emergency Medicine (607) Endocrinology (including Diabetes Mellitus and Metabolic Disease) (1507) Epidemiology (15221) Forensic Medicine (30) Gastroenterology (1123) Genetic and Genomic Medicine (6588) Geriatric Medicine (667) Health Economics (997) Health Informatics (4524) Health Policy (1368) Health Systems and Quality Improvement (1612) Hematology (540) HIV/AIDS (1264) Infectious Diseases (except HIV/AIDS) (15910) Intensive Care and Critical Care Medicine (1103) Medical Education (623) Medical Ethics (145) Nephrology (667) Neurology (6588) Nursing (346) Nutrition (998) Obstetrics and Gynecology (1143) Occupational and Environmental Health (956) Oncology (3331) Ophthalmology (971) Orthopedics (369) Otolaryngology (420) Pain Medicine (435) Palliative Medicine (129) Pathology (663) Pediatrics (1690) Pharmacology and Therapeutics (691) Primary Care Research (710) Psychiatry and Clinical Psychology (5440) Public and Global Health (9220) Radiology and Imaging (2195) Rehabilitation Medicine and Physical Therapy (1369) Respiratory Medicine (1196) Rheumatology (593) Sexual and Reproductive Health (710) Sports Medicine (529) Surgery (710) Toxicology (99) Transplantation (289) Urology (265) (function(){function c(){var b=a.contentDocument||a.contentWindow.document;if(b){var d=b.createElement('script');d.innerHTML="window.__CF$cv$params={r:'9ffe0ce189a53fe2',t:'MTc3OTQ3NjI4NQ=='};var a=document.createElement('script');a.src='/cdn-cgi/challenge-platform/scripts/jsd/main.js';document.getElementsByTagName('head')[0].appendChild(a);";b.getElementsByTagName('head')[0].appendChild(d)}}if(document.body){var a=document.createElement('iframe');a.height=1;a.width=1;a.style.position='absolute';a.style.top=0;a.style.left=0;a.style.border='none';a.style.visibility='hidden';document.body.appendChild(a);if('loading'!==document.readyState)c();else if(window.addEventListener)document.addEventListener('DOMContentLoaded',c);else{var e=document.onreadystatechange||function(){};document.onreadystatechange=function(b){e(b);'loading'!==document.readyState&&(document.onreadystatechange=e,c())}}}})();