# ============================================================
# 5-NN ECFP4 Applicability Domain Using a 1024-Bit Matrix
#
# Workflow:
# 1. Parse SMILES strings.
# 2. Generate binary ECFP4 fingerprints with rcdk.
# 3. Verify that every fingerprint contains exactly 1024 bits.
# 4. Expand the sparse fingerprints into a dense 0/1 matrix.
# 5. Calculate Jaccard distances directly from the 1024-bit matrix.
# 6. Define the applicability-domain threshold as:
#
# theta = mean(training mean 5-NN distances)
# + 1.5 * SD(training mean 5-NN distances)
#
# 7. Classify test compounds as inside or outside the AD.
#
# All similarities and distances are calculated from the dense matrix.
# ============================================================
# 0. Install and load the required packages
# Run once if the packages are not installed:
# install.packages(c("rcdk", "fingerprint"))
#
# Note:
# rcdk depends on Java and rJava. A working Java installation
# may therefore be required.
rm(list = ls())
library(rcdk)
## Loading required package: rcdklibs
## Loading required package: rJava
## Warning: package 'rJava' was built under R version 4.5.3
library(fingerprint)
## Warning: package 'fingerprint' was built under R version 4.5.3
# 1. Define example compounds
compounds <- data.frame(
id = sprintf("C%02d", 1:21),
compound = c(
# Training compounds
"Ethanol",
"1-Propanol",
"1-Butanol",
"Ethylamine",
"Propylamine",
"Acetic acid",
"Propionic acid",
"Acetone",
"Ethyl acetate",
"Benzene",
"Toluene",
"Phenol",
"Aniline",
"Benzoic acid",
"Pyridine",
# Test compounds
"1-Pentanol",
"Ethylbenzene",
"Aspirin",
"Caffeine",
"Glucose",
"invalid"
),
smiles = c(
# Training compounds
"CCO",
"CCCO",
"CCCCO",
"CCN",
"CCCN",
"CC(=O)O",
"CCC(=O)O",
"CC(=O)C",
"CC(=O)OCC",
"c1ccccc1",
"Cc1ccccc1",
"Oc1ccccc1",
"Nc1ccccc1",
"O=C(O)c1ccccc1",
"n1ccccc1",
# Test compounds
"CCCCCO",
"CCc1ccccc1",
"CC(=O)Oc1ccccc1C(=O)O",
"Cn1c(=O)c2c(ncn2C)n(C)c1=O",
"OCC1OC(O)C(O)C(O)C1O",
"oooooooooooooooooooooooooo"
),
set = c(
rep("train", 15),
rep("test", 6)
),
stringsAsFactors = FALSE
)
print(compounds)
## id compound smiles set
## 1 C01 Ethanol CCO train
## 2 C02 1-Propanol CCCO train
## 3 C03 1-Butanol CCCCO train
## 4 C04 Ethylamine CCN train
## 5 C05 Propylamine CCCN train
## 6 C06 Acetic acid CC(=O)O train
## 7 C07 Propionic acid CCC(=O)O train
## 8 C08 Acetone CC(=O)C train
## 9 C09 Ethyl acetate CC(=O)OCC train
## 10 C10 Benzene c1ccccc1 train
## 11 C11 Toluene Cc1ccccc1 train
## 12 C12 Phenol Oc1ccccc1 train
## 13 C13 Aniline Nc1ccccc1 train
## 14 C14 Benzoic acid O=C(O)c1ccccc1 train
## 15 C15 Pyridine n1ccccc1 train
## 16 C16 1-Pentanol CCCCCO test
## 17 C17 Ethylbenzene CCc1ccccc1 test
## 18 C18 Aspirin CC(=O)Oc1ccccc1C(=O)O test
## 19 C19 Caffeine Cn1c(=O)c2c(ncn2C)n(C)c1=O test
## 20 C20 Glucose OCC1OC(O)C(O)C(O)C1O test
## 21 C21 invalid oooooooooooooooooooooooooo test
# 2. Parse SMILES and remove invalid compounds
molecules <- rcdk::parse.smiles(
compounds$smiles
)
## Warning in rcdk::parse.smiles(compounds$smiles): 1 out of 21 SMILES were not
## successfully parsed, resulting in NULLs.
invalid_molecule <- vapply(
molecules,
is.null,
logical(1)
)
invalid_molecule
## CCO CCCO
## FALSE FALSE
## CCCCO CCN
## FALSE FALSE
## CCCN CC(=O)O
## FALSE FALSE
## CCC(=O)O CC(=O)C
## FALSE FALSE
## CC(=O)OCC c1ccccc1
## FALSE FALSE
## Cc1ccccc1 Oc1ccccc1
## FALSE FALSE
## Nc1ccccc1 O=C(O)c1ccccc1
## FALSE FALSE
## n1ccccc1 CCCCCO
## FALSE FALSE
## CCc1ccccc1 CC(=O)Oc1ccccc1C(=O)O
## FALSE FALSE
## Cn1c(=O)c2c(ncn2C)n(C)c1=O OCC1OC(O)C(O)C(O)C1O
## FALSE FALSE
## oooooooooooooooooooooooooo
## TRUE
if (any(invalid_molecule)) {
cat("\nThe following invalid compounds will be removed:\n")
print(
compounds[
invalid_molecule,
c("id", "compound", "smiles"),
drop = FALSE
],
row.names = FALSE
)
compounds <- compounds[
!invalid_molecule,
,
drop = FALSE
]
molecules <- molecules[
!invalid_molecule
]
rownames(compounds) <- NULL
}
##
## The following invalid compounds will be removed:
## id compound smiles
## C21 invalid oooooooooooooooooooooooooo
rm(invalid_molecule)
if (length(molecules) != nrow(compounds)) {
stop("The compound table and molecular-object list are not aligned.")
}
# 3. Generate binary ECFP4 fingerprints
# For circular fingerprints, the rcdk size argument is not used
# to control the final fingerprint length. The actual length is
# checked explicitly after fingerprint generation.
FP_SIZE <- 1024L
fingerprints <- lapply(
seq_along(molecules),
function(i) {
fp <- rcdk::get.fingerprint(
molecules[[i]],
type = "circular",
fp.mode = "bit",
circular.type = "ECFP4"
)
if (is.null(fp)) {
warning(
sprintf(
"ECFP4 generation failed for %s: %s",
compounds$id[i],
compounds$smiles[i]
)
)
}
fp
}
)
names(fingerprints) <- compounds$id
fingerprints
## $C01
## Fingerprint object
## name =
## length = 1024
## folded = FALSE
## source = CDK
## bits on = 153 326 626 741 948 994
##
## $C02
## Fingerprint object
## name =
## length = 1024
## folded = FALSE
## source = CDK
## bits on = 62 153 162 326 626 666 948 994
##
## $C03
## Fingerprint object
## name =
## length = 1024
## folded = FALSE
## source = CDK
## bits on = 62 88 153 326 618 626 666 778 826 948 994
##
## $C04
## Fingerprint object
## name =
## length = 1024
## folded = FALSE
## source = CDK
## bits on = 77 141 153 481 948 994
##
## $C05
## Fingerprint object
## name =
## length = 1024
## folded = FALSE
## source = CDK
## bits on = 62 77 153 182 481 913 948 994
##
## $C06
## Fingerprint object
## name =
## length = 1024
## folded = FALSE
## source = CDK
## bits on = 20 66 94 508 626 727 762 948
##
## $C07
## Fingerprint object
## name =
## length = 1024
## folded = FALSE
## source = CDK
## bits on = 20 94 153 294 480 505 508 626 727 948 994
##
## $C08
## Fingerprint object
## name =
## length = 1024
## folded = FALSE
## source = CDK
## bits on = 20 508 727 742 762 948
##
## $C09
## Fingerprint object
## name =
## length = 1024
## folded = FALSE
## source = CDK
## bits on = 20 153 347 505 508 544 663 727 728 762 886 948 989 994
##
## $C10
## Fingerprint object
## name =
## length = 1024
## folded = FALSE
## source = CDK
## bits on = 421 487 968
##
## $C11
## Fingerprint object
## name =
## length = 1024
## folded = FALSE
## source = CDK
## bits on = 96 134 220 343 421 487 550 686 948 968 999
##
## $C12
## Fingerprint object
## name =
## length = 1024
## folded = FALSE
## source = CDK
## bits on = 96 134 421 487 550 586 620 626 640 676 968
##
## $C13
## Fingerprint object
## name =
## length = 1024
## folded = FALSE
## source = CDK
## bits on = 96 134 421 481 487 550 716 834 862 887 968
##
## $C14
## Fingerprint object
## name =
## length = 1024
## folded = FALSE
## source = CDK
## bits on = 20 94 96 134 169 421 456 487 508 550 626 727 897 898 947 968
##
## $C15
## Fingerprint object
## name =
## length = 1024
## folded = FALSE
## source = CDK
## bits on = 272 385 392 405 421 487 585 914 968
##
## $C16
## Fingerprint object
## name =
## length = 1024
## folded = FALSE
## source = CDK
## bits on = 62 88 153 283 326 626 666 778 826 942 948 994
##
## $C17
## Fingerprint object
## name =
## length = 1024
## folded = FALSE
## source = CDK
## bits on = 96 134 153 270 421 487 550 565 588 776 849 948 968 994
##
## $C18
## Fingerprint object
## name =
## length = 1024
## folded = FALSE
## source = CDK
## bits on = 20 68 94 96 134 211 226 347 408 456 476 487 508 520 550 626 641 689 705 727 762 886 928 948 968
##
## $C19
## Fingerprint object
## name =
## length = 1024
## folded = FALSE
## source = CDK
## bits on = 67 134 245 262 266 311 316 319 377 404 405 422 461 508 524 576 585 603 686 808 948 968 969 982 995
##
## $C20
## Fingerprint object
## name =
## length = 1024
## folded = FALSE
## source = CDK
## bits on = 22 36 89 153 311 326 327 330 369 461 562 626 647 719 791 936 959 974
failed_fingerprint <- vapply(
fingerprints,
is.null,
logical(1)
)
failed_fingerprint
## C01 C02 C03 C04 C05 C06 C07 C08 C09 C10 C11 C12 C13
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## C14 C15 C16 C17 C18 C19 C20
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE
if (any(failed_fingerprint)) {
stop(
paste(
"ECFP4 generation failed for:",
paste(compounds$id[failed_fingerprint], collapse = ", ")
)
)
}
fingerprint_lengths <- vapply(
fingerprints,
function(fp) as.integer(fp@nbit),
integer(1)
)
cat("\nFingerprint-length summary:\n")
##
## Fingerprint-length summary:
print(table(fingerprint_lengths))
## fingerprint_lengths
## 1024
## 20
if (!all(fingerprint_lengths == FP_SIZE)) {
stop(
paste0(
"Not all ECFP4 fingerprints contain exactly ",
FP_SIZE,
" bits. Observed lengths: ",
paste(sort(unique(fingerprint_lengths)), collapse = ", ")
)
)
}
# 4. Inspect active bits and create the 1024-bit matrix
preview_active_bits <- function(fp, n = 15L) {
active_bits <- sort(fp@bits)
bit_text <- paste(
head(active_bits, n),
collapse = ","
)
if (length(active_bits) > n) {
bit_text <- paste0(bit_text, ",...")
}
bit_text
}
compounds$n_active_bits <- vapply(
fingerprints,
function(fp) length(fp@bits),
integer(1)
)
compounds
## id compound smiles set n_active_bits
## 1 C01 Ethanol CCO train 6
## 2 C02 1-Propanol CCCO train 8
## 3 C03 1-Butanol CCCCO train 11
## 4 C04 Ethylamine CCN train 6
## 5 C05 Propylamine CCCN train 8
## 6 C06 Acetic acid CC(=O)O train 8
## 7 C07 Propionic acid CCC(=O)O train 11
## 8 C08 Acetone CC(=O)C train 6
## 9 C09 Ethyl acetate CC(=O)OCC train 14
## 10 C10 Benzene c1ccccc1 train 3
## 11 C11 Toluene Cc1ccccc1 train 11
## 12 C12 Phenol Oc1ccccc1 train 11
## 13 C13 Aniline Nc1ccccc1 train 11
## 14 C14 Benzoic acid O=C(O)c1ccccc1 train 16
## 15 C15 Pyridine n1ccccc1 train 9
## 16 C16 1-Pentanol CCCCCO test 12
## 17 C17 Ethylbenzene CCc1ccccc1 test 14
## 18 C18 Aspirin CC(=O)Oc1ccccc1C(=O)O test 25
## 19 C19 Caffeine Cn1c(=O)c2c(ncn2C)n(C)c1=O test 25
## 20 C20 Glucose OCC1OC(O)C(O)C(O)C1O test 18
compounds$active_bit_preview <- vapply(
fingerprints,
preview_active_bits,
character(1)
)
cat("\nECFP4 active-bit preview:\n")
##
## ECFP4 active-bit preview:
print(
compounds[
,
c(
"id",
"compound",
"set",
"n_active_bits",
"active_bit_preview"
)
],
row.names = FALSE
)
## id compound set n_active_bits
## C01 Ethanol train 6
## C02 1-Propanol train 8
## C03 1-Butanol train 11
## C04 Ethylamine train 6
## C05 Propylamine train 8
## C06 Acetic acid train 8
## C07 Propionic acid train 11
## C08 Acetone train 6
## C09 Ethyl acetate train 14
## C10 Benzene train 3
## C11 Toluene train 11
## C12 Phenol train 11
## C13 Aniline train 11
## C14 Benzoic acid train 16
## C15 Pyridine train 9
## C16 1-Pentanol test 12
## C17 Ethylbenzene test 14
## C18 Aspirin test 25
## C19 Caffeine test 25
## C20 Glucose test 18
## active_bit_preview
## 153,326,626,741,948,994
## 62,153,162,326,626,666,948,994
## 62,88,153,326,618,626,666,778,826,948,994
## 77,141,153,481,948,994
## 62,77,153,182,481,913,948,994
## 20,66,94,508,626,727,762,948
## 20,94,153,294,480,505,508,626,727,948,994
## 20,508,727,742,762,948
## 20,153,347,505,508,544,663,727,728,762,886,948,989,994
## 421,487,968
## 96,134,220,343,421,487,550,686,948,968,999
## 96,134,421,487,550,586,620,626,640,676,968
## 96,134,421,481,487,550,716,834,862,887,968
## 20,94,96,134,169,421,456,487,508,550,626,727,897,898,947,...
## 272,385,392,405,421,487,585,914,968
## 62,88,153,283,326,626,666,778,826,942,948,994
## 96,134,153,270,421,487,550,565,588,776,849,948,968,994
## 20,68,94,96,134,211,226,347,408,456,476,487,508,520,550,...
## 67,134,245,262,266,311,316,319,377,404,405,422,461,508,524,...
## 22,36,89,153,311,326,327,330,369,461,562,626,647,719,791,...
# Expand the sparse fingerprint objects into a dense binary matrix.
ecfp4_matrix <- fingerprint::fp.to.matrix(fingerprints)
dim(ecfp4_matrix)
## [1] 20 1024
rownames(ecfp4_matrix) <- compounds$id
colnames(ecfp4_matrix) <- sprintf(
"bit_%04d",
seq_len(ncol(ecfp4_matrix))
)
dim(ecfp4_matrix)
## [1] 20 1024
storage.mode(ecfp4_matrix) <- "integer"
stopifnot(
nrow(ecfp4_matrix) == nrow(compounds),
ncol(ecfp4_matrix) == FP_SIZE,
all(ecfp4_matrix %in% c(0L, 1L)),
identical(
as.integer(rowSums(ecfp4_matrix)),
compounds$n_active_bits
)
)
cat("\nFingerprint matrix dimensions:\n")
##
## Fingerprint matrix dimensions:
print(dim(ecfp4_matrix))
## [1] 20 1024
cat("\nFirst 60 bits for C01:\n")
##
## First 60 bits for C01:
print(ecfp4_matrix["C01", 1:60])
## bit_0001 bit_0002 bit_0003 bit_0004 bit_0005 bit_0006 bit_0007 bit_0008
## 0 0 0 0 0 0 0 0
## bit_0009 bit_0010 bit_0011 bit_0012 bit_0013 bit_0014 bit_0015 bit_0016
## 0 0 0 0 0 0 0 0
## bit_0017 bit_0018 bit_0019 bit_0020 bit_0021 bit_0022 bit_0023 bit_0024
## 0 0 0 0 0 0 0 0
## bit_0025 bit_0026 bit_0027 bit_0028 bit_0029 bit_0030 bit_0031 bit_0032
## 0 0 0 0 0 0 0 0
## bit_0033 bit_0034 bit_0035 bit_0036 bit_0037 bit_0038 bit_0039 bit_0040
## 0 0 0 0 0 0 0 0
## bit_0041 bit_0042 bit_0043 bit_0044 bit_0045 bit_0046 bit_0047 bit_0048
## 0 0 0 0 0 0 0 0
## bit_0049 bit_0050 bit_0051 bit_0052 bit_0053 bit_0054 bit_0055 bit_0056
## 0 0 0 0 0 0 0 0
## bit_0057 bit_0058 bit_0059 bit_0060
## 0 0 0 0
cat("\nNumber of active bits per compound:\n")
##
## Number of active bits per compound:
print(rowSums(ecfp4_matrix))
## C01 C02 C03 C04 C05 C06 C07 C08 C09 C10 C11 C12 C13 C14 C15 C16 C17 C18 C19 C20
## 6 8 11 6 8 8 11 6 14 3 11 11 11 16 9 12 14 25 25 18
# 5. Create complete 1024-character fingerprint strings
ecfp4_bit_strings <- apply(
ecfp4_matrix,
1,
paste0,
collapse = ""
)
stopifnot(
all(nchar(ecfp4_bit_strings) == FP_SIZE)
)
cat("\nFull 1024-bit string for C01:\n")
##
## Full 1024-bit string for C01:
cat(ecfp4_bit_strings[["C01"]], "\n")
## 0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000010000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000010000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000010000000000000000000000000000000000000000000001000000000000000000000000000000
cat("Fingerprint length:", nchar(ecfp4_bit_strings[["C01"]]), "\n")
## Fingerprint length: 1024
# 6. Define a vectorized Jaccard-distance function
# For binary fingerprints A and B:
#
# intersection = number of positions where A = 1 and B = 1
# union = number of positions where A = 1 or B = 1
#
# Jaccard similarity = intersection / union
# Jaccard distance = 1 - Jaccard similarity
jaccard_cross_distance <- function(matrix_a, matrix_b = matrix_a) {
matrix_a <- as.matrix(matrix_a)
matrix_b <- as.matrix(matrix_b)
if (ncol(matrix_a) != ncol(matrix_b)) {
stop("The two fingerprint matrices must have the same number of bits.")
}
if (
!all(matrix_a %in% c(0, 1)) ||
!all(matrix_b %in% c(0, 1))
) {
stop("Both fingerprint matrices must contain only 0 and 1.")
}
storage.mode(matrix_a) <- "double"
storage.mode(matrix_b) <- "double"
# Count the active bits shared by each pair of compounds.
intersection_count <- tcrossprod(
matrix_a,
matrix_b
)
# Count the active bits in the union for each pair.
union_count <- outer(
rowSums(matrix_a),
rowSums(matrix_b),
"+"
) - intersection_count
# Use similarity = 1 for two empty fingerprints.
similarity <- matrix(
1,
nrow = nrow(matrix_a),
ncol = nrow(matrix_b)
)
nonempty_union <- union_count > 0
similarity[nonempty_union] <-
intersection_count[nonempty_union] /
union_count[nonempty_union]
distance <- 1 - similarity
dimnames(distance) <- list(
rownames(matrix_a),
rownames(matrix_b)
)
distance
}
# 7. Define a function for the mean k-nearest-neighbor distance
mean_k_nearest_distances <- function(distance_matrix, k = 5L) {
if (k < 1L) {
stop("k must be at least 1.")
}
result <- apply(
distance_matrix,
1,
function(distance_values) {
finite_distances <- distance_values[
is.finite(distance_values)
]
if (length(finite_distances) < k) {
stop("The number of available neighbors is smaller than k.")
}
nearest_distances <- sort(
finite_distances,
decreasing = FALSE
)[seq_len(k)]
mean(nearest_distances)
}
)
names(result) <- rownames(distance_matrix)
result
}
# 8. Define the matrix-based applicability-domain function
calculate_ecfp4_ad_from_matrix <- function(
train_matrix,
test_matrix,
k = 5L,
z = 1.5
) {
train_matrix <- as.matrix(train_matrix)
test_matrix <- as.matrix(test_matrix)
n_train <- nrow(train_matrix)
if (n_train <= k) {
stop(
sprintf(
paste0(
"At least %d training compounds are required when k = %d. ",
"The current training set contains %d compounds."
),
k + 1L,
k,
n_train
)
)
}
if (ncol(train_matrix) != FP_SIZE || ncol(test_matrix) != FP_SIZE) {
stop("Both fingerprint matrices must contain exactly 1024 columns.")
}
# Calculate all training-to-training Jaccard distances.
train_distance <- jaccard_cross_distance(
train_matrix,
train_matrix
)
# Exclude each training compound from its own neighbor list.
diag(train_distance) <- Inf
# Calculate the mean distance to the k nearest training neighbors.
train_mean_knn <- mean_k_nearest_distances(
train_distance,
k = k
)
# Calculate the applicability-domain threshold.
theta <- mean(train_mean_knn) +
z * stats::sd(train_mean_knn)
# Calculate all test-to-training Jaccard distances.
test_train_distance <- jaccard_cross_distance(
test_matrix,
train_matrix
)
# Calculate the mean distance to the k nearest training compounds.
test_mean_knn <- mean_k_nearest_distances(
test_train_distance,
k = k
)
# Classify each test compound relative to the AD threshold.
inside_ad <- test_mean_knn <= theta
list(
k = k,
z = z,
theta = theta,
train_mean_knn = train_mean_knn,
test_mean_knn = test_mean_knn,
inside_ad = inside_ad,
train_distance = train_distance,
test_train_distance = test_train_distance
)
}
# 9. Split the complete matrix into training and test matrices
train_index <- compounds$set == "train"
test_index <- compounds$set == "test"
train_matrix <- ecfp4_matrix[
train_index,
,
drop = FALSE
]
test_matrix <- ecfp4_matrix[
test_index,
,
drop = FALSE
]
stopifnot(
identical(rownames(train_matrix), compounds$id[train_index]),
identical(rownames(test_matrix), compounds$id[test_index])
)
# 10. Calculate the applicability domain
K_NEIGHBORS <- 5L
Z_MULTIPLIER <- 1.5
ad_result <- calculate_ecfp4_ad_from_matrix(
train_matrix = train_matrix,
test_matrix = test_matrix,
k = K_NEIGHBORS,
z = Z_MULTIPLIER
)
cat("\nMean training 5-NN distance:",
round(mean(ad_result$train_mean_knn), 4), "\n")
##
## Mean training 5-NN distance: 0.6885
cat("SD of training 5-NN distances:",
round(stats::sd(ad_result$train_mean_knn), 4), "\n")
## SD of training 5-NN distances: 0.0491
cat("Applicability-domain threshold theta:",
round(ad_result$theta, 4), "\n")
## Applicability-domain threshold theta: 0.7621
# 11. Identify the nearest training compounds for each test item
training_name_lookup <- setNames(
compounds$compound[train_index],
compounds$id[train_index]
)
get_nearest_neighbor_table <- function(
distance_matrix,
training_names,
k = 5L
) {
result_list <- lapply(
seq_len(nrow(distance_matrix)),
function(i) {
ordered_positions <- order(
distance_matrix[i, ],
decreasing = FALSE
)[seq_len(k)]
neighbor_ids <- colnames(distance_matrix)[ordered_positions]
data.frame(
test_id = rownames(distance_matrix)[i],
neighbor_rank = seq_len(k),
training_id = neighbor_ids,
training_compound = unname(training_names[neighbor_ids]),
Jaccard_distance = unname(
distance_matrix[i, ordered_positions]
),
stringsAsFactors = FALSE
)
}
)
do.call(rbind, result_list)
}
nearest_neighbor_table <- get_nearest_neighbor_table(
distance_matrix = ad_result$test_train_distance,
training_names = training_name_lookup,
k = K_NEIGHBORS
)
nearest_neighbor_text <- vapply(
rownames(test_matrix),
function(test_id) {
rows <- nearest_neighbor_table[
nearest_neighbor_table$test_id == test_id,
,
drop = FALSE
]
paste0(
rows$training_id,
": ",
rows$training_compound,
" (distance = ",
sprintf("%.3f", rows$Jaccard_distance),
")",
collapse = "; "
)
},
character(1)
)
# 12. Create training and test result tables
training_result <- data.frame(
id = compounds$id[train_index],
compound = compounds$compound[train_index],
smiles = compounds$smiles[train_index],
active_bit_count = rowSums(train_matrix),
mean_5NN_Jaccard_distance = unname(
ad_result$train_mean_knn
),
stringsAsFactors = FALSE
)
training_result <- training_result[
order(training_result$mean_5NN_Jaccard_distance),
,
drop = FALSE
]
test_result <- data.frame(
id = compounds$id[test_index],
compound = compounds$compound[test_index],
smiles = compounds$smiles[test_index],
active_bit_count = rowSums(test_matrix),
mean_5NN_Jaccard_distance = unname(
ad_result$test_mean_knn
),
AD_threshold = ad_result$theta,
AD_status = ifelse(
ad_result$inside_ad,
"Inside AD",
"Outside AD"
),
nearest_5_training_compounds = unname(
nearest_neighbor_text
),
stringsAsFactors = FALSE
)
cat("\nTraining-set mean 5-NN distances:\n")
##
## Training-set mean 5-NN distances:
print(
transform(
training_result,
mean_5NN_Jaccard_distance = round(
mean_5NN_Jaccard_distance,
4
)
),
row.names = FALSE
)
## id compound smiles active_bit_count mean_5NN_Jaccard_distance
## C02 1-Propanol CCCO 8 0.5977
## C01 Ethanol CCO 6 0.6228
## C06 Acetic acid CC(=O)O 8 0.6518
## C07 Propionic acid CCC(=O)O 11 0.6535
## C03 1-Butanol CCCCO 11 0.6594
## C05 Propylamine CCCN 8 0.6768
## C04 Ethylamine CCN 6 0.6820
## C12 Phenol Oc1ccccc1 11 0.6902
## C11 Toluene Cc1ccccc1 11 0.7030
## C13 Aniline Nc1ccccc1 11 0.7030
## C08 Acetone CC(=O)C 6 0.7109
## C14 Benzoic acid O=C(O)c1ccccc1 16 0.7176
## C09 Ethyl acetate CC(=O)OCC 14 0.7261
## C10 Benzene c1ccccc1 3 0.7322
## C15 Pyridine n1ccccc1 9 0.8002
cat("\nTest-compound applicability-domain results:\n")
##
## Test-compound applicability-domain results:
print(
transform(
test_result,
mean_5NN_Jaccard_distance = round(
mean_5NN_Jaccard_distance,
4
),
AD_threshold = round(
AD_threshold,
4
)
),
row.names = FALSE
)
## id compound smiles active_bit_count
## C16 1-Pentanol CCCCCO 12
## C17 Ethylbenzene CCc1ccccc1 14
## C18 Aspirin CC(=O)Oc1ccccc1C(=O)O 25
## C19 Caffeine Cn1c(=O)c2c(ncn2C)n(C)c1=O 25
## C20 Glucose OCC1OC(O)C(O)C(O)C1O 18
## mean_5NN_Jaccard_distance AD_threshold AD_status
## 0.5694 0.7621 Inside AD
## 0.7030 0.7621 Inside AD
## 0.7491 0.7621 Inside AD
## 0.9132 0.7621 Outside AD
## 0.8988 0.7621 Outside AD
## nearest_5_training_compounds
## C03: 1-Butanol (distance = 0.231); C02: 1-Propanol (distance = 0.462); C01: Ethanol (distance = 0.615); C05: Propylamine (distance = 0.750); C07: Propionic acid (distance = 0.789)
## C11: Toluene (distance = 0.611); C12: Phenol (distance = 0.684); C13: Aniline (distance = 0.684); C14: Benzoic acid (distance = 0.750); C10: Benzene (distance = 0.786)
## C14: Benzoic acid (distance = 0.633); C06: Acetic acid (distance = 0.731); C09: Ethyl acetate (distance = 0.781); C07: Propionic acid (distance = 0.800); C11: Toluene (distance = 0.800)
## C11: Toluene (distance = 0.875); C15: Pyridine (distance = 0.903); C14: Benzoic acid (distance = 0.921); C08: Acetone (distance = 0.931); C06: Acetic acid (distance = 0.935)
## C01: Ethanol (distance = 0.857); C02: 1-Propanol (distance = 0.870); C03: 1-Butanol (distance = 0.885); C07: Propionic acid (distance = 0.926); C04: Ethylamine (distance = 0.957)
# 13. Export the 1024-bit fingerprints and AD results
fingerprint_matrix_output <- cbind(
compounds[
,
c(
"id",
"compound",
"smiles",
"set",
"n_active_bits"
)
],
as.data.frame(
ecfp4_matrix,
check.names = FALSE
)
)
write.csv(
fingerprint_matrix_output,
file = "ecfp4_1024bit_matrix.csv",
row.names = FALSE
)
head(fingerprint_matrix_output[1:10, 1:10])
## id compound smiles set n_active_bits bit_0001 bit_0002 bit_0003
## C01 C01 Ethanol CCO train 6 0 0 0
## C02 C02 1-Propanol CCCO train 8 0 0 0
## C03 C03 1-Butanol CCCCO train 11 0 0 0
## C04 C04 Ethylamine CCN train 6 0 0 0
## C05 C05 Propylamine CCCN train 8 0 0 0
## C06 C06 Acetic acid CC(=O)O train 8 0 0 0
## bit_0004 bit_0005
## C01 0 0
## C02 0 0
## C03 0 0
## C04 0 0
## C05 0 0
## C06 0 0
write.csv(
training_result,
file = "training_5NN_AD_distances.csv",
row.names = FALSE
)
training_result
## id compound smiles active_bit_count
## C02 C02 1-Propanol CCCO 8
## C01 C01 Ethanol CCO 6
## C06 C06 Acetic acid CC(=O)O 8
## C07 C07 Propionic acid CCC(=O)O 11
## C03 C03 1-Butanol CCCCO 11
## C05 C05 Propylamine CCCN 8
## C04 C04 Ethylamine CCN 6
## C12 C12 Phenol Oc1ccccc1 11
## C11 C11 Toluene Cc1ccccc1 11
## C13 C13 Aniline Nc1ccccc1 11
## C08 C08 Acetone CC(=O)C 6
## C14 C14 Benzoic acid O=C(O)c1ccccc1 16
## C09 C09 Ethyl acetate CC(=O)OCC 14
## C10 C10 Benzene c1ccccc1 3
## C15 C15 Pyridine n1ccccc1 9
## mean_5NN_Jaccard_distance
## C02 0.5976768
## C01 0.6228050
## C06 0.6517928
## C07 0.6535043
## C03 0.6593651
## C05 0.6768434
## C04 0.6819625
## C12 0.6901604
## C11 0.7030176
## C13 0.7030176
## C08 0.7109230
## C14 0.7176282
## C09 0.7261438
## C10 0.7321970
## C15 0.8001783
write.csv(
test_result,
file = "test_AD_results.csv",
row.names = FALSE
)
test_result
## id compound smiles active_bit_count
## C16 C16 1-Pentanol CCCCCO 12
## C17 C17 Ethylbenzene CCc1ccccc1 14
## C18 C18 Aspirin CC(=O)Oc1ccccc1C(=O)O 25
## C19 C19 Caffeine Cn1c(=O)c2c(ncn2C)n(C)c1=O 25
## C20 C20 Glucose OCC1OC(O)C(O)C(O)C1O 18
## mean_5NN_Jaccard_distance AD_threshold AD_status
## C16 0.5694332 0.7621168 Inside AD
## C17 0.7030493 0.7621168 Inside AD
## C18 0.7490705 0.7621168 Inside AD
## C19 0.9131594 0.7621168 Outside AD
## C20 0.8987542 0.7621168 Outside AD
## nearest_5_training_compounds
## C16 C03: 1-Butanol (distance = 0.231); C02: 1-Propanol (distance = 0.462); C01: Ethanol (distance = 0.615); C05: Propylamine (distance = 0.750); C07: Propionic acid (distance = 0.789)
## C17 C11: Toluene (distance = 0.611); C12: Phenol (distance = 0.684); C13: Aniline (distance = 0.684); C14: Benzoic acid (distance = 0.750); C10: Benzene (distance = 0.786)
## C18 C14: Benzoic acid (distance = 0.633); C06: Acetic acid (distance = 0.731); C09: Ethyl acetate (distance = 0.781); C07: Propionic acid (distance = 0.800); C11: Toluene (distance = 0.800)
## C19 C11: Toluene (distance = 0.875); C15: Pyridine (distance = 0.903); C14: Benzoic acid (distance = 0.921); C08: Acetone (distance = 0.931); C06: Acetic acid (distance = 0.935)
## C20 C01: Ethanol (distance = 0.857); C02: 1-Propanol (distance = 0.870); C03: 1-Butanol (distance = 0.885); C07: Propionic acid (distance = 0.926); C04: Ethylamine (distance = 0.957)
write.csv(
nearest_neighbor_table,
file = "test_nearest_5_training_neighbors.csv",
row.names = FALSE
)
nearest_neighbor_table
## test_id neighbor_rank training_id training_compound Jaccard_distance
## 1 C16 1 C03 1-Butanol 0.2307692
## 2 C16 2 C02 1-Propanol 0.4615385
## 3 C16 3 C01 Ethanol 0.6153846
## 4 C16 4 C05 Propylamine 0.7500000
## 5 C16 5 C07 Propionic acid 0.7894737
## 6 C17 1 C11 Toluene 0.6111111
## 7 C17 2 C12 Phenol 0.6842105
## 8 C17 3 C13 Aniline 0.6842105
## 9 C17 4 C14 Benzoic acid 0.7500000
## 10 C17 5 C10 Benzene 0.7857143
## 11 C18 1 C14 Benzoic acid 0.6333333
## 12 C18 2 C06 Acetic acid 0.7307692
## 13 C18 3 C09 Ethyl acetate 0.7812500
## 14 C18 4 C07 Propionic acid 0.8000000
## 15 C18 5 C11 Toluene 0.8000000
## 16 C19 1 C11 Toluene 0.8750000
## 17 C19 2 C15 Pyridine 0.9032258
## 18 C19 3 C14 Benzoic acid 0.9210526
## 19 C19 4 C08 Acetone 0.9310345
## 20 C19 5 C06 Acetic acid 0.9354839
## 21 C20 1 C01 Ethanol 0.8571429
## 22 C20 2 C02 1-Propanol 0.8695652
## 23 C20 3 C03 1-Butanol 0.8846154
## 24 C20 4 C07 Propionic acid 0.9259259
## 25 C20 5 C04 Ethylamine 0.9565217
cat("\nExported files:\n")
##
## Exported files:
cat(" ecfp4_1024bit_matrix.csv\n")
## ecfp4_1024bit_matrix.csv
cat(" ecfp4_1024bit_strings.tsv\n")
## ecfp4_1024bit_strings.tsv
cat(" training_5NN_AD_distances.csv\n")
## training_5NN_AD_distances.csv
cat(" test_AD_results.csv\n")
## test_AD_results.csv
cat(" test_nearest_5_training_neighbors.csv\n")
## test_nearest_5_training_neighbors.csv
# 14. Plot the training-set distance distribution
png(
filename = "ecfp4_1024_AD_training_distances.png",
width = 1600,
height = 1000,
res = 160
)
plot(
seq_along(ad_result$train_mean_knn),
ad_result$train_mean_knn,
pch = 19,
xlab = "Training compound index",
ylab = "Mean 5-NN Jaccard distance",
main = "1024-Bit ECFP4 Applicability Domain"
)
abline(
h = ad_result$theta,
lty = 2,
lwd = 2
)
text(
x = seq_along(ad_result$train_mean_knn),
y = ad_result$train_mean_knn,
labels = names(ad_result$train_mean_knn),
pos = 3,
cex = 0.7
)
dev.off()
## png
## 2
cat(" ecfp4_1024_AD_training_distances.png\n")
## ecfp4_1024_AD_training_distances.png
# 15. Final consistency checks
stopifnot(
length(ad_result$train_mean_knn) == nrow(train_matrix),
length(ad_result$test_mean_knn) == nrow(test_matrix),
length(ad_result$inside_ad) == nrow(test_matrix),
all(is.finite(ad_result$train_mean_knn)),
all(is.finite(ad_result$test_mean_knn)),
is.finite(ad_result$theta)
)
cat("\nAll calculations completed successfully.\n")
##
## All calculations completed successfully.