Using Random Forest to Identity Host or Microbial Biomarkers of Non-optimal Vaginal Microbiomes
Jenna Wilson
2026-01-15
Introduction/README
Extensive literature suggests a Lactobacillus dominant cervicovaginal microbiome plays a protective role in sexual and reproductive health [1]. However, there are significant limitations to the tests currently in clinical use for identifying non-optimal microbiome states.
The aim of this project is to identify protein biomarkers of non-optimal vaginal microbiomes using a metaproteomics dataset from pregnant Australian women.
In this analysis, an “optimal” microbiome is considered to have more than 50% relative abundance of one or a combination of Lactobacillus species. “Non-optimal” includes all others, mostly Gardnerella, Bifidobacterium and diverse microbiome states.
Getting started
There are two input files required. A metaproteomics dataset that has been partially pre-processed (see below) and a file containing the samples and their case/control status for the variable of interest (in this case, an optimal (0) or non-optimal (1) vaginal microbiome).
Packages required: library(dplyr) library(caret) library(randomForest) library(kableExtra) library(knitr) library(ggplot2) (optional, for visualisation)
Computational requirements: This workflow was run on an Intel(R) Core(TM) Ultra 7 165U (2.10 GHz) laptop with 16GB RAM, 64-bit operating system, x64-based processor. These code chunks should take less than a minute to run.
Methodology and data
The metaproteomics dataset used here was generated by the Monash Proteomics and Metabolomics Platform in Clayton. This is data dependent acquisition using a TMT multiplex system (31 plexes/batches).
Data pre-processing (in brief):
Poor quality outliers were removed (using a PCA plot)
Low confidence and low abundance proteins filtered (proteins identified by only one unique peptide were excluded, and a criterion of minimum intensity of 100 across 1% of samples was applied)
Contaminants were removed (bovine, pig, mouse proteins introuced during mass spectrometry sample preparation)
These data will be uploaded to the PRIDE database https://www.ebi.ac.uk/pride/ but until then, will be made available by email.
Part 1: Normalisation, Data splitting and Imputation
First, load and inspect datasets
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(caret)## Loading required package: ggplot2## Loading required package: latticelibrary(knitr)
library(kableExtra)##
## Attaching package: 'kableExtra'## The following object is masked from 'package:dplyr':
##
## group_rowsset.seed(12345) #for reproducibility
#1.Load dataset and additional sheet that contains case/control labels
proteomics_data <- read.delim("proteomic_data_filtered_refs_included2_no_contam.csv",
header = TRUE,
sep = ",",
stringsAsFactors = FALSE)
labels <- read.csv("Lacto_Nonlacto_labelling.csv", stringsAsFactors = FALSE) #columns: SampleID, Class
labels$SampleID <- as.character(labels$SampleID)
kable(labels)| SampleID | Class |
|---|---|
| X11250 | 0 |
| X11352 | 0 |
| X11140 | 0 |
| X11273 | 0 |
| X11269 | 0 |
| X11335 | 0 |
| X11173 | 0 |
| X10993 | 0 |
| X11099 | 0 |
| X11107 | 0 |
| X11136 | 0 |
| X11089 | 0 |
| X21144 | 0 |
| X11186 | 0 |
| X11067 | 0 |
| X21096 | 0 |
| X11000 | 0 |
| X11122 | 0 |
| X11007 | 0 |
| X11791 | 0 |
| X11620 | 0 |
| X11102 | 0 |
| X21386 | 0 |
| X21205 | 0 |
| X20831 | 0 |
| X21442 | 0 |
| X11035 | 0 |
| X21116 | 0 |
| X11135 | 0 |
| X11113 | 0 |
| X21357 | 0 |
| X11055 | 0 |
| X21330 | 0 |
| X11808 | 0 |
| X11041 | 0 |
| X11353 | 0 |
| X11016 | 0 |
| X21069 | 0 |
| X20884 | 0 |
| X21198 | 0 |
| X11005 | 0 |
| X21191 | 0 |
| X10942 | 0 |
| X10975 | 0 |
| X11126 | 0 |
| X11180 | 0 |
| X11251 | 0 |
| X11158 | 0 |
| X11339 | 0 |
| X11019 | 0 |
| X11294 | 0 |
| X20931 | 0 |
| X21094 | 0 |
| X11257 | 0 |
| X11013 | 0 |
| X21445 | 0 |
| X11004 | 0 |
| X11103 | 0 |
| X11422 | 0 |
| X21435 | 0 |
| X11021 | 0 |
| X11523 | 0 |
| X11264 | 0 |
| X20837 | 0 |
| X11012 | 0 |
| X21365 | 0 |
| X11159 | 0 |
| X11146 | 0 |
| X11108 | 0 |
| X11077 | 0 |
| X11052 | 0 |
| X11028 | 0 |
| X11142 | 0 |
| X11026 | 0 |
| X10964 | 0 |
| X11425 | 0 |
| X11162 | 0 |
| X11117 | 0 |
| X11211 | 0 |
| X20964 | 0 |
| X21334 | 0 |
| X11254 | 0 |
| X21252 | 0 |
| X11178 | 1 |
| X11078 | 1 |
| X11458 | 1 |
| X11109 | 1 |
| X11023 | 1 |
| X20898 | 1 |
| X11412 | 1 |
| X11011 | 1 |
| X10946 | 1 |
| X11316 | 1 |
| X10954 | 1 |
| X11087 | 1 |
| X21162 | 1 |
| X11528 | 1 |
| X11161 | 1 |
| X11459 | 1 |
| X11643 | 1 |
| X11048 | 1 |
| X11084 | 1 |
| X11128 | 1 |
| X11134 | 1 |
| X11152 | 1 |
| X11104 | 1 |
| X21335 | 1 |
| X11133 | 1 |
| X21027 | 1 |
| X21347 | 1 |
| X20930 | 1 |
| X21495 | 1 |
| X21308 | 1 |
| X10945 | 1 |
# 2. Split metadata and quantitative data
metadata_cols <- 33
meta <- proteomics_data[, 1:metadata_cols] #this contains protein info (e.g. Accession, used later)
quant <- proteomics_data[, (metadata_cols + 1):ncol(proteomics_data)] #sample columns
head(meta, 3)## Checked Protein.FDR.Confidence..Combined Master Accession
## 1 FALSE High Master Protein D4MJJ1
## 2 FALSE High Master Protein P15636
## 3 FALSE High Master Protein A0A2I1K755
## Species.Names
## 1 [Eubacterium] siraeum V10Sc8a OX=717961
## 2 Achromobacter lyticus
## 3 Aerococcus christensenii OX=87541
## Description
## 1 Phosphoglycerate kinase OS=[Eubacterium] siraeum V10Sc8a OX=717961 GN=pgk PE=3 SV=1
## 2 SWISS-PROT:P15636 Protease I precursor Lysyl endopeptidase Achromobacter lyticus.
## 3 Peptidase OS=Aerococcus christensenii OX=87541 GN=CYJ27_05215 PE=4 SV=1
## Exp..q.value..Combined Sum.PEP.Score Coverage.... X..Peptides X..PSMs
## 1 0 5.058 7 2 92
## 2 0 60.381 20 11 7347
## 3 0 42.901 8 13 107
## X..Unique.Peptides X..AAs MW..kDa. calc..pI
## 1 2 406 43.3 5.50
## 2 11 653 68.1 7.24
## 3 11 2350 257.2 6.14
## Score.MS.Amanda.3.0..MS.Amanda.3.0 Score.Sequest.HT..Sequest.HT
## 1 8715.02 66.00
## 2 484165.67 12226.99
## 3 8429.03 159.60
## X..Peptides..by.Search.Engine...MS.Amanda.3.0
## 1 2
## 2 9
## 3 13
## X..Peptides..by.Search.Engine...Sequest.HT Contaminant Marked.as
## 1 2 FALSE Bacterial
## 2 11 TRUE
## 3 13 FALSE Bacterial
## Biological.Process Cellular.Component
## 1 other metabolic processes cytosol
## 2 protein metabolism non-structural extracellular
## 3 non-structural extracellular;other membranes
## Molecular.Function Pfam.IDs
## 1 kinase activity;other molecular function
## 2 other molecular function
## 3 Pf00746, Pf04650, Pf07554
## Entrez.Gene.ID Ensembl.Gene.ID Gene.ID Gene.Symbol
## 1 D4MJJ1; esr:ES1_08450 pgk
## 2
## 3 A0A2I1K755
## Reactome.Pathways WikiPathways X..Protein.Pathway.Groups X..Razor.Peptides
## 1 0 0
## 2 0 0
## 3 0 2head(quant, 3)## REF1 X21488 X11142 X21042 X20855 X11787 X11772 X21005 X11269 X20968 X21369
## 1 NA NA NA NA NA NA NA NA NA NA NA
## 2 1265.1 1649.2 1596.9 1012.3 974.0 1548.9 1339.7 965.1 935.9 1266.3 873.4
## 3 104.8 18.6 5.3 7.1 9.3 8.1 14.8 142.9 18.7 50.7 7.6
## X21171 X21094 X11261 X20958 X21542 X21328 X21506 REF2 X21059 X11603 X11839
## 1 NA NA NA NA NA NA NA NA NA NA NA
## 2 1152.1 1291.3 1427.8 1091.8 1211.8 1152.5 1313.3 1958.3 1545.1 2526.7 916.7
## 3 6.8 6.8 374.9 513.2 7.1 0.1 4154.3 NA NA NA NA
## X11590 X21040 X11778 X11655 X21044 X21382 X11319 X21159 X11128 X11161 X11007
## 1 NA NA NA NA NA NA NA NA NA NA NA
## 2 2000 2445.7 2271 1740.7 2059.1 1590.1 1654.6 2537.2 1172.1 1833.9 1499.1
## 3 NA NA NA NA NA NA NA NA NA NA NA
## X21335 X11667 X21266 REF3 X21347 X11087 X21064 X11113 X20842 X11794 X11555
## 1 NA NA NA NA NA NA NA NA NA NA NA
## 2 1780.4 2357.5 1989.7 2282.5 2110.0 2422 2255.3 2500.4 2152.1 2146.3 1859.3
## 3 NA NA NA 92.0 270.8 NA 207.4 2.6 4.2 3.7 13.7
## X20949 X20962 X11443 X20903 X20931 X11133 X11217 X11844 X21154 REF4 X11565
## 1 NA NA NA NA NA NA NA NA NA 183.8 114.2
## 2 2179.9 1939.1 1605.6 1834.6 1979.1 2394.3 2273.5 2367.9 1788.4 2467.3 1852.8
## 3 NA 306.6 43.6 2.3 3.2 28.7 5.5 2.9 3.4 42.4 63.3
## X21131 X11005 X11815 X11513 X11882 X20959 X21658 X21329 X21334 X11539 X21563
## 1 77.0 92.6 109.6 100.1 64.0 164.0 170.8 278.8 115.4 134.2 62.0
## 2 1539.1 1368.9 2550.0 1728.0 1997.3 2571.4 2131.2 2602.5 2687.4 2682.5 1530.5
## 3 1.7 1.9 4.9 9.5 5.2 7.9 17.5 16.0 173.6 75.6 280.1
## X21140 X11412 X11398 X11325 REF5 X11427 X21319 X20873 X21116 X11117 X21148
## 1 97.3 127.6 109.8 65.9 165.0 83.5 66.1 26.9 31.2 121.0 138.1
## 2 1943.3 2179.0 1677.8 2377.2 2059.2 1741.6 1684.3 1280.5 2153.6 2286.0 2170.3
## 3 13.1 753.6 19.4 123.4 35.1 NA 202.4 3.0 2.7 5.7 5.0
## X20877 X10964 X20880 X11377 X20884 X11136 X20905 X11409 X20820 X11587 REF6
## 1 62.0 121.7 55.1 94.9 47.2 54.2 48.4 54.7 103.4 64.3 554.2
## 2 1696.1 2461.1 1723.9 2014.7 2267.0 2525.1 2322.2 2051.1 2004.2 1981.2 1821.5
## 3 150.4 36.4 3.9 2.6 3.9 7.7 2.3 NA 53.5 125.2 NA
## X20982 X11850 X21453 X11273 X11449 X21331 X21184 X11310 X11281 X21095 X11173
## 1 962.9 436.6 374.0 448.6 366.9 684.9 398.8 569.3 619.1 845.0 320.8
## 2 1816.8 1933.5 1172.6 1765.4 2252.4 2162.3 1291.4 1482.5 1531.8 1580.6 1187.3
## 3 NA NA NA NA NA NA NA NA NA NA NA
## X11523 X11595 X11204 X11563 REF7 X11353 X11103 X21339 X10993 X11654 X21205
## 1 766.2 427.0 479.8 620.1 27.5 24.7 49.1 32.9 29.1 39.2 34.3
## 2 1511.0 2400.1 1462.1 1924.8 3001.1 2278.2 2275.9 1610.4 3431.6 2690.9 2375.3
## 3 NA NA NA NA 27.2 3.6 2.4 4.8 2.7 147.5 NA
## X20821 X11320 X21091 X11211 X21530 X21263 X11857 X11571 X20852 X10945 X11786
## 1 33.9 14.3 38.5 38.5 31.9 35.7 40.1 37.5 41.8 33.4 38.3
## 2 3095.6 1632.1 1887.9 2123.8 2632.5 1434.9 2560.6 1719.8 1956.2 2770.7 2550.6
## 3 3.6 123.4 6.2 27.6 22.9 2.9 7.2 NA 6.9 256.2 6.8
## REF8 X21120 X21164 X11484 X11474 X11126 X21191 X21429 X11669 X11899 X11104
## 1 67.3 82.8 72.1 71.9 77.6 73.0 51.4 10.0 33.0 79.6 47.8
## 2 2338.6 940.1 1032.1 1594.1 1113.0 1041.3 1097.2 505.1 639.2 1103.4 1266.9
## 3 36.3 1.9 1.4 NA NA NA 50.0 NA NA NA 454.2
## X11102 X11608 X20811 X11645 X11254 X11914 X11109 REF9 X11336 X21199 X21384
## 1 74.1 14.2 53.7 72.1 63.3 80.1 56.8 148.3 47.1 355.2 26.2
## 2 1080.0 2171.1 990.8 1180.3 1264.2 1339.3 1390.2 3463.1 582.2 2388.7 2452.4
## 3 3.0 4.6 1.9 2.3 1.9 6.2 3.1 199.4 1.6 NA 423.3
## X11396 X11771 X21089 X11695 X11518 X21069 X21492 X11471 X21099 X11534 X11041
## 1 90.2 161.5 117.8 NA 7.0 6.8 3.5 291.9 8.8 157.6 58.2
## 2 1452.8 1468.2 1216.6 1205.9 1557.0 2152.1 1476.2 1526.5 1760.2 1169.0 1430.3
## 3 14.3 0.7 1719.3 3.8 722.1 NA 677.1 10.4 348.1 11.7 5.9
## X11341 X11605 X21241 REF10 X21113 X11904 X20879 X21363 X21495 X21383 X21037
## 1 2.1 3.5 95.0 NA NA NA NA NA NA NA NA
## 2 1292.3 1447.7 1511.6 6175.8 1669.2 1719.5 2279.7 2014.3 2119.6 2350.9 2727.3
## 3 13.1 9.1 406.6 NA NA NA NA NA NA NA NA
## X11424 X21365 X11642 X11183 X11485 X11201 X21442 X21409 X11689 X21541 REF11
## 1 NA NA NA NA NA NA NA NA NA NA NA
## 2 1632.8 1668.9 1791.3 2262.4 1920.5 1332.6 1808.2 991.3 3372.9 2508 3208.4
## 3 NA NA NA NA NA NA NA NA NA NA NA
## X11052 X11108 X20835 X11825 X11077 X11408 X20887 X11762 X11327 X21261 X21502
## 1 NA NA NA NA NA NA NA NA NA NA NA
## 2 2038.3 3365.5 1217 1503 1985 1083 1524.6 2195 2212.2 2417 2333
## 3 NA NA NA NA NA NA NA NA NA NA NA
## X20930 X21445 X20923 REF12 X21278 X11388 X21004 X11159 X20799 X11788 X21322
## 1 NA NA NA 206.4 402.7 308.1 213.5 294.7 279.1 271.0 302.8
## 2 2350 2667.5 1989.1 3939.7 1356.8 1507.7 1379.1 2038.8 1099.6 1214.7 1955.0
## 3 NA NA NA NA NA NA NA NA NA NA NA
## X21311 X21616 X11078 X21219 X21375 X21428 X11528 X21221 X21510 REF13 X21028
## 1 202.2 182.9 203.0 365.1 183.0 314.1 274.4 326.6 297.1 NA NA
## 2 946.2 2401.8 1700.5 1674.1 720.9 1552.4 1317.5 1420.7 1708.1 3665.9 2096.6
## 3 NA NA NA NA NA NA NA NA NA 39.3 6.6
## X11816 X21068 X21274 X20809 X21074 X11264 X20822 X11777 X21301 X21622 X21388
## 1 NA NA NA NA NA NA NA NA NA NA NA
## 2 2200.2 3465.0 4643.6 3440.0 3273.1 2454.7 3484.1 2389.4 2410.1 2923.8 2213.7
## 3 3.0 3.3 7.6 16.8 6.4 3.8 7.8 3.8 1217.5 NA 2.6
## X11152 X20917 X20773 REF14 X11676 X21632 X21252 X11251 X21636 X21386 X20789
## 1 NA NA NA NA NA NA NA NA NA NA NA
## 2 3170.3 2473.2 546.4 4898.5 3906.4 2944.6 3883.9 4082.2 3398.1 4668.2 3202.3
## 3 4.7 2.1 NA 36.2 2.4 NA 26.7 NA NA NA NA
## X11860 X21256 X21254 X21121 X11425 X21093 X11359 X11916 X21656 REF15 X11757
## 1 NA NA NA NA NA NA NA NA NA NA NA
## 2 3054.4 3400.8 3401.2 2594.5 4107.8 4821.2 5490.7 4032.7 3818.3 2557.3 1220.1
## 3 58.8 6.5 86.6 NA 135.7 NA NA 153.5 NA 29.5 NA
## X21555 X21192 X11811 X11369 X21641 X11649 X11753 X11242 X20860 X21652 X21056
## 1 NA NA NA NA NA NA NA NA NA NA NA
## 2 1272.1 1903.3 1575.4 2420.7 1370.8 2432.2 1192.4 323.1 2021.7 1427.1 1381.3
## 3 NA NA 0.2 65.8 311.2 NA NA NA 32.8 NA 1.1
## X11566 X21336 X21049 X11422 X11352 REF16 X21248 X20782 X11019 X20831 X11548
## 1 NA NA NA NA NA NA NA NA NA NA NA
## 2 1445 2415 1312.7 1858.4 1703.2 2254.7 1220.4 1702.7 1364.5 2531 1271.4
## 3 NA NA 3.3 3.3 NA NA NA NA NA NA NA
## X11259 X20893 X11540 X11462 X11521 X21619 X20834 X21295 REF17 X11569 X20971
## 1 NA NA NA NA NA NA NA NA NA NA NA
## 2 935.4 2611.3 2105.1 1540.6 915.5 1762.1 1346.6 1264.6 2361.4 1183.1 1396.1
## 3 NA NA NA NA NA NA NA NA NA NA NA
## X21102 X11016 X11683 X11804 X10975 X21521 X11004 X11013 X11758 X11489 X11800
## 1 NA NA NA NA NA NA NA NA NA NA NA
## 2 1267.2 1327.2 1226.1 1516.2 1644.9 1649 1015.3 1548.5 2608.6 1914.4 2560.6
## 3 NA NA NA NA NA NA NA NA NA NA NA
## X21438 X11770 X11802 X21290 REF18 X21357 X11339 X21523 X21187 X21520 X20801
## 1 NA NA NA NA NA NA NA NA NA NA NA
## 2 1192.4 1939.8 1500.1 1874.8 1797.2 1031.3 1153.7 950.4 1073.7 1425.2 1198.8
## 3 NA NA NA NA NA NA NA NA NA NA NA
## X11692 X21104 X21242 X11048 X11456 X11681 X21265 X11636 X11140 X20980 X11107
## 1 NA NA NA NA NA NA NA NA NA NA NA
## 2 1139.8 1522.3 1085.6 678.8 1478.2 1250.2 784.3 1055.7 1414.1 855.5 1230.5
## 3 NA NA NA NA NA NA NA NA NA NA NA
## REF19 X11666 X11755 X21654 X11393 X11547 X21213 X11023 X11421 X21420 X21330
## 1 168.0 239.5 21.1 229.5 281.9 15.0 690.3 13.6 79.7 18.8 25.3
## 2 2746.2 1313.6 1273.6 2076.3 2312.2 2132.1 2547.9 1577.1 1347.7 1411.5 1855.6
## 3 NA NA NA NA NA NA NA NA NA NA NA
## X11620 X11012 X11915 X21260 X11791 X11618 X20851 REF20 X21166 X21144 X21186
## 1 55.4 413.8 102.4 19.4 16.3 26.5 13.0 NA NA NA NA
## 2 1537.9 1589.4 1220.9 1149.6 1395.8 2779.3 2930.3 5410.2 1476 2411.3 2327
## 3 NA NA NA NA NA NA NA NA NA NA NA
## X11519 X11122 X21518 X11774 X21012 X21436 X11546 X11884 X11905 X21308 X10946
## 1 NA NA NA NA NA NA NA NA NA NA NA
## 2 2918.2 2662.2 1583.4 2955 2636.8 1872.7 1268.2 1547.9 3226.3 1860.3 1979.8
## 3 NA NA NA NA NA NA NA NA NA NA NA
## X11055 X11247 X21130 REF21 X11646 X11165 X11360 X21362 X20781 X21207 X11680
## 1 NA NA NA 372.4 558.4 84.9 101.6 50.2 521.9 130.4 139.1
## 2 2347.9 2700.5 2690 3293.1 1779.7 1914.3 1694.1 1755.6 2279.3 2061.6 1995.5
## 3 NA NA NA NA NA NA NA NA NA NA NA
## X21378 X20920 X21136 X21650 X11533 X20803 REF22 X11647 X20859 X21628 X11820
## 1 485.7 330.5 467.7 61.8 156.3 369.9 35.4 NA NA 18.9 138.5
## 2 2175.4 1774.2 1845.7 2483.0 2236.8 2707.0 3215.5 1996.8 1509.7 2119.6 2280.7
## 3 NA NA NA NA NA NA NA NA NA NA NA
## X11898 X11186 X21326 X11067 X20837 X11021 X20898 X11364 X21317 X21337 X11448
## 1 NA NA 34.9 0.2 126.0 69.7 NA NA 1.8 29.8 127.8
## 2 2485.1 2844.2 2271.1 2762.5 2149.6 2107.1 1745.6 2138.8 2114.4 2480.8 2790.8
## 3 NA NA NA NA NA NA NA NA NA NA NA
## X21212 REF23 X21401 X11089 X20840 X21096 X11579 X20948 X11625 X20963 X11265
## 1 NA NA NA NA NA NA NA NA NA NA NA
## 2 1797.4 3417.8 2683.5 2273.5 3513.5 3112.3 3214.0 2205.1 2230.5 2234 2588.9
## 3 NA 76.7 2.8 2.9 6.9 3.1 186.6 1.1 995.3 NA 0.7
## X21559 X11809 X11257 X11146 X10954 X21053 X21162 X11394 REF24 X11922 X20753
## 1 NA NA NA NA NA NA NA NA NA NA NA
## 2 2861.0 2863.0 2996.9 2410.2 2740.3 2023.8 2868.7 2281.1 2196.1 1069.9 1568.0
## 3 4.9 31.3 NA 2.2 5.6 1.4 NA 361.2 33.2 106.7 0.4
## X11028 X21565 X21027 X11643 X11084 X20984 X11624 X21309 X11011 X11250 X11531
## 1 NA NA NA NA NA NA NA NA NA NA NA
## 2 1503.2 1574.7 1918.5 1810.2 1830.7 1661.6 1394.1 356.0 1602.1 1490.3 1524.8
## 3 NA NA 385.1 56.5 NA 117.6 153.1 4.8 NA NA NA
## X11311 X11316 X21253 X11418 REF25 X21236 X21522 X21092 X11688 X21129 X21410
## 1 NA NA NA NA NA NA NA NA NA NA NA
## 2 1436.6 1526.1 1742.7 1740.6 2252.8 1035.2 1895.7 1437.4 1293.7 1303.0 1895.2
## 3 NA 201.6 2.1 5.3 66.3 3.3 48.9 9.4 6.0 12.7 16.0
## X11659 X11335 X21663 X21087 X11178 X11162 X11545 X21440 X11510 X21421 X11134
## 1 NA NA NA NA NA NA NA NA NA NA NA
## 2 1674.1 1635.5 1275.8 2227.9 2277.1 1766.6 2276.8 1811.7 1557.1 1804.0 1368.6
## 3 17.8 10.8 4.6 5.0 4.6 16.6 12.7 1015.6 1.4 701.2 7.7
## REF26 X11284 X11234 X21413 X11158 X11458 X11469 X21446 X21653 X20987 X21353
## 1 80.0 132.4 NA 0.6 190.6 4.5 150.3 180.5 292.2 132.0 NA
## 2 2140.5 2045.4 1712.2 2102.6 1929.7 2132.1 1434.8 1698.4 1371.4 1338.6 1241
## 3 NA NA NA NA NA NA NA NA NA NA NA
## X11632 REF27 X11796 X11345 X20800 X11400 X11551 X11568 X21380 X11135 X11785
## 1 185.1 NA NA NA NA NA NA NA NA NA NA
## 2 1535.0 2275.1 2743.6 1206.9 1913.4 1316.7 1827.7 1778.8 1650.3 1571.6 1837.2
## 3 NA 78.4 8.9 18.2 28.2 18.4 18.4 21.8 22.3 22.4 33.7
## X21084 X10942 X11928 X11459 X11808 X11658 X11000 REF28 X21100 X21527 X20862
## 1 NA NA NA NA NA NA NA 25.3 15.5 20.8 15.8
## 2 1406.4 2166.9 1226.9 1304.9 1972.6 1655.0 1941.0 2474.3 1155.4 1350.9 1109.9
## 3 16.6 17.9 12.1 1654.0 41.2 3.1 32.4 NA NA NA NA
## X21511 X11892 X21372 X21408 X21169 X11212 X11822 X11180 X11630 X21664 X21416
## 1 17.3 13.9 31.8 19.6 29.8 12.5 28.0 26.2 27.4 15.6 15.2
## 2 1399.8 1132.8 1533.4 1661.9 1907.8 1330.6 1404.6 1189.7 1287.6 1162.6 993.0
## 3 NA NA NA NA NA NA NA NA NA NA NA
## X20952 X21423 X11821 REF29 X11827 X11099 X20906 X21498 X21043 X11035 X11648
## 1 14.2 29.6 13.2 43.6 12.6 47 14.2 63.2 48.2 18.6 4.8
## 2 1246.4 1286.6 751.7 2950.5 1605.4 2741 2498.7 1865.8 1994.3 2296.5 898.1
## 3 NA NA NA NA NA NA NA NA NA NA NA
## X21052 X11026 X11607 X20828 X11426 X20889 X21633 X21535 X11163 REF30 X11560
## 1 22.5 48.0 22.8 18.5 18.3 28.9 52.1 31.8 22.3 18.5 5.7
## 2 1826.7 1937.9 1992.6 1709.0 1389.7 2551.8 2015.7 1723.2 3207.1 2956.8 2309.8
## 3 NA NA NA NA NA NA NA NA NA 27.1 28.2
## X11833 X21390 X11910 X21047 X11662 X21454 X21367 X11526 X21371 X11830 X11328
## 1 2.9 1.5 19.3 2.9 39.8 4.7 15.5 3.2 6.6 5.2 7.6
## 2 1495.3 1996.5 3121.5 3293.5 2033.7 1904.7 1786.6 1962.0 1910.4 2238.2 2171.6
## 3 1.1 29.1 381.4 22.9 NA 3.1 7.1 4.7 3.2 NA 3.8
## X21621 X20964 X11434 X11294 X21198 REF31 X21175 X21307 X20978 X11464 X20747
## 1 6.1 15.1 3.3 40.7 3.1 45.6 70.8 58.4 74.9 NA NA
## 2 2272.6 2376.8 1900.9 2504.3 2294.7 4977.1 3302.9 2090.4 2927.4 3142.1 2318.8
## 3 12.6 5.7 28.3 11.9 7.0 NA NA NA NA NA NA
## X11760 X21177 X11368
## 1 NA 22.6 2.6
## 2 3287.1 3161.8 4003.2
## 3 NA NA NANext, perform IRS normalisation using the reference channels, and prepare data for imputation. The data is split into train and test datasets prior to imputation, to avoid data leakage/bias.
# 3. Identify Global REF columns (used for normalisation)
ref_cols <- grep("^REF", colnames(quant))
# 4. Compute global reference per protein (geometric mean of REF channels)
global_ref <- apply(quant[, ref_cols], 1, function(x){
x <- x[is.finite(x) & x > 0] #ignore 0/NA/Inf
if(length(x) == 0) return(NA_real_)
exp(mean(log(x)))
})
# 5. Apply IRS normalization plex-by-plex
quant_IRS <- quant
for(i in seq_along(ref_cols)){
ref_col <- ref_cols[i]
if(i < length(ref_cols)){
plex_cols <- ref_col:(ref_cols[i+1]-1)
} else {
plex_cols <- ref_col:ncol(quant)
}
ref_vals <- quant[, ref_col]
scale_factor <- global_ref / ref_vals
scale_factor[!is.finite(scale_factor)] <- NA
#Apply row-wise scaling
quant_IRS[, plex_cols] <- sweep(
quant[, plex_cols],
1,
scale_factor,
`*`
)
}
# 6. Convert to numeric and log2 transform (necessary before imputation)
quant_IRS <- as.data.frame(lapply(quant_IRS, function(x) as.numeric(as.character(x))))
quant_IRS_log <- log2(quant_IRS)
quant_IRS_log[!is.finite(as.matrix(quant_IRS_log))] <- NA
# 7. Remove REF channels (not needed downstream)
quant_clean <- quant_IRS_log[, -ref_cols, drop=FALSE]
# 8. Match labels for the remaining samples
sample_names <- colnames(quant_clean)
sample_classes <- labels$Class[match(sample_names, labels$SampleID)]
# Remove unlabeled samples (e.g. REF channels and those without matching 16S data)
valid_idx <- which(!is.na(sample_classes))
sample_names_valid <- sample_names[valid_idx]
sample_classes_valid <- sample_classes[valid_idx]
# 9. Stratified train/test split (70% train, 30% test)
set.seed(12345) #for reproducibility
train_idx <- createDataPartition(
y = sample_classes_valid,
p = 0.7,
list = FALSE
)
train_samples <- sample_names_valid[train_idx]
test_samples <- sample_names_valid[-train_idx]
train_classes <- sample_classes_valid[train_idx]
test_classes <- sample_classes_valid[-train_idx]
# Create train/test matrices
train_matrix <- quant_clean[, train_samples, drop = FALSE]
test_matrix <- quant_clean[, test_samples, drop = FALSE]Perseus-style imputation was chosen as it is a left shifted method that does not impute values from the upper distribution, which is appropriate for differential expression and random forest downstream analyses.
# 10. "Perseus-style" imputation
perseus_impute <- function(x, width=0.3, downshift=1.8){
observed <- x[!is.na(x)]
if(length(observed) < 2) return(x)
mean_obs <- mean(observed)
sd_obs <- sd(observed)
if(is.na(sd_obs) || sd_obs==0) sd_obs <- 1e-6
mu <- mean_obs - downshift*sd_obs
sigma <- sd_obs*width
x[is.na(x)] <- rnorm(sum(is.na(x)), mean=mu, sd=sigma)
return(x)
}
set.seed(12345)
imputed_train <- apply(train_matrix, 2, perseus_impute)
imputed_test <- apply(test_matrix, 2, perseus_impute)
#stop to check how many values were missing (NAs), and how many get imputed. We want these to match.
sum(is.na(train_matrix))## [1] 136160sum(is.na(test_matrix))## [1] 57583sum(is.na(train_matrix) & !is.na(imputed_train))## [1] 136160sum(is.na(test_matrix) & !is.na(imputed_test))## [1] 57583# 11. Transpose so rows = samples, columns = proteins
train_final <- as.data.frame(t(imputed_train))
test_final <- as.data.frame(t(imputed_test))
# 12. Add protein/sample labels
protein_labels <- meta$Accession
# Column names = protein labels
colnames(train_final) <- protein_labels
colnames(test_final) <- protein_labels
# Add Class and SampleID columns
train_final$Class <- train_classes
train_final$SampleID <- train_samples
test_final$Class <- test_classes
test_final$SampleID <- test_samples
# 13. Optional: save processed datasets to CSV and inspect output
write.csv(train_final, "train_data_imputed_with_labels_lacto2.csv", row.names = FALSE)
write.csv(test_final, "test_data_imputed_with_labels_lacto2.csv", row.names = FALSE)Part 2: Random Forest for predicting optimal (0) and non-optimal (1) vaginal microbiomes
We want to: A) measure the sensitivity and specificity, B) Idenitfy the top 20 most important features (potential biomarkers)
# 1. Load packages
library(randomForest)## randomForest 4.7-1.2## Type rfNews() to see new features/changes/bug fixes.##
## Attaching package: 'randomForest'## The following object is masked from 'package:ggplot2':
##
## margin## The following object is masked from 'package:dplyr':
##
## combinelibrary(dplyr)
library(ggplot2)
# 2. Prepare data for RF
train_final$Class <- as.factor(train_final$Class)
test_final$Class <- as.factor(test_final$Class)
feature_cols <- setdiff(colnames(train_final), c("SampleID", "Class"))
x_train <- train_final[, feature_cols]
y_train <- train_final$Class
x_test <- test_final[, feature_cols]
y_test <- test_final$Class
# 3. Check class distribution of train dataset
table(y_train)## y_train
## 0 1
## 55 25# 4. Train Random Forest with class balancing (two options were suggested, here I have gone with method 1)
set.seed(12345)
rf_model <- randomForest(
x = x_train,
y = y_train,
ntree = 500,
importance = TRUE,
# Method 1: class weights
classwt = c("0" = 1, "1" = 5), # give minority class higher weight
# Method 2: balanced sampling
#sampsize = c("0" = 25, "1" = 25) # each tree sees equal numbers of each class
)
# 5. Evaluate model
pred_test <- predict(rf_model, x_test)
conf_matrix <- table(Predicted = pred_test, Actual = y_test)
print(conf_matrix)## Actual
## Predicted 0 1
## 0 27 3
## 1 0 3accuracy <- mean(pred_test == y_test)
cat("Test set accuracy:", accuracy, "\n")## Test set accuracy: 0.9090909#Calculate sensitivity / specificity
sensitivity <- conf_matrix["1","1"] / sum(conf_matrix[,"1"])
specificity <- conf_matrix["0","0"] / sum(conf_matrix[,"0"])
cat("Sensitivity (cases):", sensitivity, "\n")## Sensitivity (cases): 0.5cat("Specificity (controls):", specificity, "\n")## Specificity (controls): 1# 6. Feature importance
importance_df <- data.frame(
Protein = rownames(rf_model$importance),
MeanDecreaseGini = rf_model$importance[, "MeanDecreaseGini"]
) %>% arrange(desc(MeanDecreaseGini))
#Top 20 proteins (potential biomarkers)
top20 <- head(importance_df, 20)
print(top20)## Protein MeanDecreaseGini
## A0A3E1IPY6 A0A3E1IPY6 0.3053048
## A0A9X7FES1 A0A9X7FES1 0.2426053
## A0A9X7FEU3 A0A9X7FEU3 0.2416014
## C8PAY4 C8PAY4 0.2282421
## A0A9X7FEX1 A0A9X7FEX1 0.2190282
## Q5TDH0 Q5TDH0 0.2179694
## C8PB23 C8PB23 0.2073075
## K1P574 K1P574 0.1943948
## A0AAQ2QIE5 A0AAQ2QIE5 0.1909355
## A0A9X7FEZ8 A0A9X7FEZ8 0.1900298
## C8PB65 C8PB65 0.1897678
## C8PAL6 C8PAL6 0.1846343
## C8PB25 C8PB25 0.1794131
## A0A133NWM8 A0A133NWM8 0.1723751
## C8PCZ3 C8PCZ3 0.1661992
## A0A805Z0S7 A0A805Z0S7 0.1607945
## C8PCK4 C8PCK4 0.1578538
## A0A9X7I9K3 A0A9X7I9K3 0.1573021
## A0A9X7FER2 A0A9X7FER2 0.1549394
## K1MQ15 K1MQ15 0.1378049#Plot top 20 proteins using ggplot2
ggplot(top20, aes(x = reorder(Protein, MeanDecreaseGini), y = MeanDecreaseGini)) +
geom_bar(stat="identity", fill="skyblue") +
coord_flip() +
theme_minimal() +
labs(
title = "Top 20 protein biomarkers (Balanced RF)",
y = "Mean Decrease in Gini",
x = "Protein Accession"
)
Discussion
The random forest model generated here demonstrates good overall accuracy, although the specificity is limited. This may be due to the diverse nature of the outcome group. The inclusion of inflammation or clinical information may improve the model accuracy and clinical relevance of these findings. Interestingly, a recent study by Challa et al. [2] used machine learning algorithms and multi-omic data to investigate novel diagnostic biomarkers of BV, and found metabolomic models performed better than microbial ones.
Reflection
While troubleshooting the script, the most common cause of error codes seemed to be data in the wrong format (numeric, list, character etc.). While I feel a lot more comfortable working with basic R code, including changing colours etc. I find the random forest script difficult to interpret and I’m not always sure if the AI generated code is performing what I expect it to. With more practice I hope to be less reliant on AI for the R code.
I have really enjoyed the workshop series. Workshop sessions 2 (Data and code management), 4 (Mastering Rstudio with Rmarkdown reports and version control) and 10 (Machine learning in R with caret) were especially helpful for this assignment, along with the R markdown cheat sheet https://github.com/adam-p/markdown-here/wiki/Markdown-Cheatsheet. Thank you!
Impact on professional development
This workshop series has been an enormous help for my PhD at the Burnet, which focuses on metaproteomics data analysis of vaginal samples from pregnant participants. Most of my analysis is being done in R studio, so these sessions were especially relevant and the use of R Markdown has made my work more organised and readable to my group members. Additionally, I find the figures made in R look more professional than the ones I would usually make in Excel. Many of the skills and principles taught in the workshop series (data management, reproducibility) are also transferable and employable in other fields, not just -omics related work.
Bibliography
[1]Plummer EL, Vodstrcil LA, Bradshaw CS. Unravelling the vaginal microbiome, impact on health and disease. Curr Opin Obstet Gynecol. 2024 Oct 1;36(5):338-344. doi: 10.1097/GCO.0000000000000976.
[2]Challa A, Maras JS, Nagpal S, Tripathi G, Taneja B, Kachhawa G, et al. Multi-omics analysis identifies potential microbial and metabolite diagnostic biomarkers of bacterial vaginosis. J Eur Acad Dermatol Venereol. 2024; 38: 1152–1165. https://doi.org/10.1111/jdv.19805
sessionInfo()## R version 4.5.2 (2025-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.3 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_AU.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_AU.UTF-8 LC_COLLATE=en_AU.UTF-8
## [5] LC_MONETARY=en_AU.UTF-8 LC_MESSAGES=en_AU.UTF-8
## [7] LC_PAPER=en_AU.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_AU.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Australia/Melbourne
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] randomForest_4.7-1.2 kableExtra_1.4.0 knitr_1.51
## [4] caret_7.0-1 lattice_0.22-7 ggplot2_4.0.1
## [7] dplyr_1.1.4
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.6 xfun_0.56 bslib_0.9.0
## [4] recipes_1.3.1 vctrs_0.7.0 tools_4.5.2
## [7] generics_0.1.4 stats4_4.5.2 parallel_4.5.2
## [10] tibble_3.3.1 pkgconfig_2.0.3 ModelMetrics_1.2.2.2
## [13] Matrix_1.7-4 data.table_1.18.0 RColorBrewer_1.1-3
## [16] S7_0.2.1 lifecycle_1.0.5 compiler_4.5.2
## [19] farver_2.1.2 stringr_1.6.0 textshaping_1.0.4
## [22] codetools_0.2-20 htmltools_0.5.9 class_7.3-23
## [25] sass_0.4.10 yaml_2.3.12 prodlim_2025.04.28
## [28] pillar_1.11.1 jquerylib_0.1.4 MASS_7.3-65
## [31] cachem_1.1.0 gower_1.0.2 iterators_1.0.14
## [34] rpart_4.1.24 foreach_1.5.2 nlme_3.1-168
## [37] parallelly_1.46.1 lava_1.8.2 tidyselect_1.2.1
## [40] digest_0.6.39 stringi_1.8.7 future_1.69.0
## [43] reshape2_1.4.5 purrr_1.2.1 listenv_0.10.0
## [46] labeling_0.4.3 splines_4.5.2 fastmap_1.2.0
## [49] grid_4.5.2 cli_3.6.5 magrittr_2.0.4
## [52] dichromat_2.0-0.1 survival_3.8-6 future.apply_1.20.1
## [55] withr_3.0.2 scales_1.4.0 lubridate_1.9.4
## [58] timechange_0.3.0 rmarkdown_2.30 globals_0.18.0
## [61] otel_0.2.0 nnet_7.3-20 timeDate_4051.111
## [64] evaluate_1.0.5 hardhat_1.4.2 viridisLite_0.4.2
## [67] rlang_1.1.7 Rcpp_1.1.1 glue_1.8.0
## [70] xml2_1.5.2 pROC_1.19.0.1 ipred_0.9-15
## [73] svglite_2.2.2 rstudioapi_0.18.0 jsonlite_2.0.0
## [76] R6_2.6.1 plyr_1.8.9 systemfonts_1.3.1