Protein Aptamer Deconvolution Simulation
Alex Abbas
2024-05-21
Introduction
This report is an exercise in simulating a toy version of a protein assay inspired by the aptamer-based Nautilus device. The proteome used is the collection of human blood-secreted proteins from the Protein Atlas. All aptamers are tripeptides, and 10% of possible tripeptides here actually have an aptamer to be used for detection.
Data
Proteins
Human proteins from the Protein Atlas known to be secreted to blood:
hpa = read.csv("data/proteinatlas.csv")
hpa = hpa %>% dplyr::select(Gene, Ensembl, Gene.description, Uniprot, Secretome.location)
hpa = hpa %>% mutate(SecretedBlood = Secretome.location == "Secreted to blood")
secretome = hpa %>% dplyr::filter(SecretedBlood) %>% dplyr::select(-Secretome.location,-SecretedBlood)
# build secretome library and aptamer library
edb <- EnsDb.Hsapiens.v86
prts1 <- proteins(edb, filter = GeneNameFilter(secretome$Gene), columns = c("gene_id","protein_sequence")) %>% data.frame
prts = prts1 %>% group_by(gene_name) %>% dplyr::slice(1)
meanConc = read.csv("data/meanConc.txt",sep="\t",header=F) # from
colnames(meanConc) = c("gene_name","description","X40.g.L.1","x")
meanConc = meanConc %>% dplyr::select(-x)
meanConc <- meanConc %>% separate(X40.g.L.1, c('MeanConc', 'MeanUnit'), sep = " ") %>% mutate(MeanConc = as.numeric(MeanConc))
multHash = c(
"g/L" = 1e9,
"mg/L" = 1e6,
"µg/L" = 1e3,
"ng/L" = 1
)
meanConc = meanConc %>% mutate(MeanVal = MeanConc * multHash[MeanUnit])
prts = inner_join(prts,meanConc)
datatable(prts %>% dplyr::select(-protein_sequence, -MeanVal))Aptamers
Here we create the aptamers. 1% of possible aptamers are available for use. Binding affinities of aptamers are assumed to be experimentally measured with a standard deviation of log10 0.6. A differential affinity of specific to non-specific aptamer binding is 1 * 104.
# set parameters
onTargetAffinity = 1e-5
offTargetAffinity = 1e-1
conc = 1e-3
pctAptamerCoverage = 0.1
affinityStdDev = 0.6
get_kmers <- function(string, k=3) {
if (k > nchar(string)) { stop("k cannot be greater than the length of the string") } # Check if k is greater than the length of the string
kmers <- character(0) # Initialize an empty vector to store the kmers
for (i in 1:(nchar(string) - k + 1)) { # Loop through the string to extract kmers
kmer <- substr(string, i, i + k - 1)
kmers <- c(kmers, kmer)
}
kmers
}
trimers = sapply(prts$protein_sequence, get_kmers) %>% unlist %>% unique
aptamers = sample(trimers, size = floor(length(trimers)/2))
onTargetAffinities = rnorm(length(aptamers), mean = log10(onTargetAffinity), sd = affinityStdDev)
names(onTargetAffinities) = aptamers
offTargetAffinities = rnorm(length(aptamers), mean = log10(offTargetAffinity), sd = affinityStdDev)
names(offTargetAffinities) = aptamers3949 aptamers are made.
Plot the distribution of observed aptamer affinities for all tripeptides. The large peak is low affinity non-specific binding, and the small peak around -6 is the high affinity on-target binding.
affinityMatrix = do.call(cbind,lapply(prts$protein_sequence, function(aseq) {
aseq_trimer_table = table(get_kmers(aseq))
aseq_trimer_apt = aseq_trimer_table[names(aseq_trimer_table) %in% aptamers]
aseq_trimer_apt[setdiff(aptamers,names(aseq_trimer_apt))] = 0
all_aptamers_counts = aseq_trimer_apt[aptamers]
affinities = offTargetAffinities + onTargetAffinities * all_aptamers_counts
affinities
}))
colnames(affinityMatrix) = prts$gene_name
probMatrix = 1 / (1 + exp(1) ^ (affinityMatrix - log10(conc)))
plot(density(affinityMatrix),xlim=c(-9,1))
Probability fit
Let’s fit a two-group gaussian mixed model to the binding probabilities. Fit on a 10000-aptamer random subset for efficiency.
dataToFit = sample(unlist(probMatrix),1e4)
fit = Mclust(dataToFit, G=2, model="V")
summary(fit)----------------------------------------------------
Gaussian finite mixture model fitted by EM algorithm
----------------------------------------------------
Mclust V (univariate, unequal variance) model with 2 components:
log-likelihood n df BIC ICL
10605 10000 5 21163.96 20774.05
Clustering table:
1 2
9369 631 plot(fit, what="density", xlab = "Probability")
Analysis
Make test sample
Now let’s create one test sample of 100 proteins on 100 tiles (100% occupancy!), with proteins selected at random according to typical normal healthy blood serum protein relative abundances. Mostly albumin!
numOfSpots = 1e2
testProteins = sample(prts$gene_name, size = numOfSpots, prob = prts$MeanVal / max(prts$MeanVal), replace = T)
aptamerHits = do.call(cbind, lapply(testProteins, function(testProtein) {
sapply(aptamers, function(aptamer) {
p = probMatrix[aptamer,testProtein]
sample(c(1,0), 1, prob = c(p,1-p))
})
}))
table(testProteins)testProteins
ALB APOA1 APOA2 APOB APOC2 APOH C4BPB C9
54 3 1 1 1 1 1 1
CFH FGG HP ORM1 SERPINA1 SERPINA3 TF TTR
1 2 2 1 1 1 25 2
VTN
2 Assess test sample
Finally, let’s run our virtual assay and empirically identify the proteins by the discrete values of aptamers that bound or did not bind to each protein, by computing a joint probability of each observed aptamer binding event for each possible protein identity. The raw score is just the pearson correlation of the observed aptamer binding states and the aptamer binding probabilities of the top match. The marginal p-value is the log10 difference between the probability of the top hit and the probability of the second place hit (so, the higher the better).
cors = cor(aptamerHits, probMatrix)
testResults = bind_rows(lapply(c(1:100),function(test_i){
aCors = cors[test_i,]
maxCor = max(aCors)
secondCor = sort(aCors,decreasing = T)[2]
fit = fitdistr(aCors, "normal")
para = fit$estimate
p1 = pnorm(maxCor, mean = para["mean"], sd = para["sd"], lower.tail=F)
p2 = pnorm(secondCor, mean = para["mean"], sd = para["sd"], lower.tail=F)
names(p2) = NULL
assignment = prts$gene_name[aCors==maxCor]
c("True_protein" = testProteins[test_i], "Inferred_protein" = assignment, "Score" = maxCor, "Marginal_Pval" = -log10(p1) + log10(p2))
}))
datatable(testResults)And a confusion table of the results.
confusion = confusionMatrix(data = factor(testResults$Inferred_protein), reference = factor(testResults$True_protein))
confusionConfusion Matrix and Statistics
Reference
Prediction ALB APOA1 APOA2 APOB APOC2 APOH C4BPB C9 CFH FGG HP ORM1 SERPINA1
ALB 54 0 0 0 0 0 0 0 0 0 0 0 0
APOA1 0 3 0 0 0 0 0 0 0 0 0 0 0
APOA2 0 0 1 0 0 0 0 0 0 0 0 0 0
APOB 0 0 0 1 0 0 0 0 0 0 0 0 0
APOC2 0 0 0 0 1 0 0 0 0 0 0 0 0
APOH 0 0 0 0 0 1 0 0 0 0 0 0 0
C4BPB 0 0 0 0 0 0 1 0 0 0 0 0 0
C9 0 0 0 0 0 0 0 1 0 0 0 0 0
CFH 0 0 0 0 0 0 0 0 1 0 0 0 0
FGG 0 0 0 0 0 0 0 0 0 2 0 0 0
HP 0 0 0 0 0 0 0 0 0 0 2 0 0
ORM1 0 0 0 0 0 0 0 0 0 0 0 1 0
SERPINA1 0 0 0 0 0 0 0 0 0 0 0 0 1
SERPINA3 0 0 0 0 0 0 0 0 0 0 0 0 0
TF 0 0 0 0 0 0 0 0 0 0 0 0 0
TTR 0 0 0 0 0 0 0 0 0 0 0 0 0
VTN 0 0 0 0 0 0 0 0 0 0 0 0 0
Reference
Prediction SERPINA3 TF TTR VTN
ALB 0 0 0 0
APOA1 0 0 0 0
APOA2 0 0 0 0
APOB 0 0 0 0
APOC2 0 0 0 0
APOH 0 0 0 0
C4BPB 0 0 0 0
C9 0 0 0 0
CFH 0 0 0 0
FGG 0 0 0 0
HP 0 0 0 0
ORM1 0 0 0 0
SERPINA1 0 0 0 0
SERPINA3 1 0 0 0
TF 0 25 0 0
TTR 0 0 2 0
VTN 0 0 0 2
Overall Statistics
Accuracy : 1
95% CI : (0.9638, 1)
No Information Rate : 0.54
P-Value [Acc > NIR] : < 2.2e-16
Kappa : 1
Mcnemar's Test P-Value : NA
Statistics by Class:
Class: ALB Class: APOA1 Class: APOA2 Class: APOB
Sensitivity 1.00 1.00 1.00 1.00
Specificity 1.00 1.00 1.00 1.00
Pos Pred Value 1.00 1.00 1.00 1.00
Neg Pred Value 1.00 1.00 1.00 1.00
Prevalence 0.54 0.03 0.01 0.01
Detection Rate 0.54 0.03 0.01 0.01
Detection Prevalence 0.54 0.03 0.01 0.01
Balanced Accuracy 1.00 1.00 1.00 1.00
Class: APOC2 Class: APOH Class: C4BPB Class: C9 Class: CFH
Sensitivity 1.00 1.00 1.00 1.00 1.00
Specificity 1.00 1.00 1.00 1.00 1.00
Pos Pred Value 1.00 1.00 1.00 1.00 1.00
Neg Pred Value 1.00 1.00 1.00 1.00 1.00
Prevalence 0.01 0.01 0.01 0.01 0.01
Detection Rate 0.01 0.01 0.01 0.01 0.01
Detection Prevalence 0.01 0.01 0.01 0.01 0.01
Balanced Accuracy 1.00 1.00 1.00 1.00 1.00
Class: FGG Class: HP Class: ORM1 Class: SERPINA1
Sensitivity 1.00 1.00 1.00 1.00
Specificity 1.00 1.00 1.00 1.00
Pos Pred Value 1.00 1.00 1.00 1.00
Neg Pred Value 1.00 1.00 1.00 1.00
Prevalence 0.02 0.02 0.01 0.01
Detection Rate 0.02 0.02 0.01 0.01
Detection Prevalence 0.02 0.02 0.01 0.01
Balanced Accuracy 1.00 1.00 1.00 1.00
Class: SERPINA3 Class: TF Class: TTR Class: VTN
Sensitivity 1.00 1.00 1.00 1.00
Specificity 1.00 1.00 1.00 1.00
Pos Pred Value 1.00 1.00 1.00 1.00
Neg Pred Value 1.00 1.00 1.00 1.00
Prevalence 0.01 0.25 0.02 0.02
Detection Rate 0.01 0.25 0.02 0.02
Detection Prevalence 0.01 0.25 0.02 0.02
Balanced Accuracy 1.00 1.00 1.00 1.00