Study note on differential privacy

2026-03-06T00:00:00+00:00

Linking attacks in genomic studies</a></h2>
Several studies have shown that releasing gene expression data can leak private information about individuals (Harmanci & Gerstein 2016</a>; Walker et al. 2024</a>).
This is how an adversary would do it:

A study releases anonymized expression data with sensitive disease labels (e.g., HIV status) </li>

Gain access to eQTL summary statistics by examining how genotypes correlate with expression </li>

Predict genotypes from expression vectors by quantizing eQTL-predictable values to match possible 0 (homozygous reference), 1 (heterozygous), 2 (homozygous alternative) alleles. </li>

Match the predicted genotypes (bins) to a reference genotype database </li>

A successful match over many SNPs will probabilistically identify a targeted individual </li> </ol>
Walker et al. showed 93%+ linking accuracy on single-cell data, even with all the inherent interference.
Can we release a bit noisy yet safer eQTL summary statistics?</a></h2>
One way to mitigate this is to release differentially private eQTL summary statistics instead of exact ones. Can we add enough noise to prevent linking while keeping the statistics useful?
Let's build a toy example:

$n = 100$ individuals</li>
One SNP with minor allele frequency 0.3</li>
Gene expression as a function of genotype: $E = \beta \cdot G + \epsilon$ (an eQTL)</li> </ul>
`library(ggplot2) set.seed(42) </code></pre>`
n <- 100 maf <- 0.3 ## genotype (private) and expression (released in the study) G <- rbinom(n, size = 2, prob = maf) E <- 0.5 G + rnorm(n) dat <- data.frame(id = 1:n, genotype = G, expression = round(E, 3)) dat </code></pre>
## id genotype expression ## 1 1 2 1.322 ## 2 2 2 0.216 ## 3 3 0 1.576 ## 4 4 1 1.143 ## 5 5 1 0.590 ## 6 6 1 0.777 ## 7 7 1 1.179 ## 8 8 0 0.090 ## 9 9 1 -2.493 ## 10 10 1 0.785 ## 11 11 0 -0.367 ## 12 12 1 0.685 ## 13 13 2 1.582 ## 14 14 0 1.400 ## 15 15 0 -0.727 ## 16 16 2 2.303 ## 17 17 2 1.336 ## 18 18 0 1.039 ## 19 19 0 0.921 ## 20 20 1 1.221 ## 21 21 1 -0.543 ## 22 22 0 -0.090 ## 23 23 2 1.624 ## 24 24 2 0.046 ## 25 25 0 -0.543 ## 26 26 1 1.081 ## 27 27 0 0.768 ## 28 28 1 0.964 ## 29 29 0 -0.886 ## 30 30 1 -0.600 ## 31 31 1 2.013 ## 32 32 1 0.758 ## 33 33 0 0.088 ## 34 34 1 0.379 ## 35 35 0 -1.194 ## 36 36 1 1.112 ## 37 37 0 -0.217 ## 38 38 0 -0.183 ## 39 39 1 1.433 ## 40 40 1 1.322 ## 41 41 0 1.392 ## 42 42 0 -0.476 ## 43 43 0 0.650 ## 44 44 2 2.391 ## 45 45 0 -1.111 ## 46 46 2 0.139 ## 47 47 1 -0.632 ## 48 48 1 -0.959 ## 49 49 2 1.080 ## 50 50 1 1.153 ## 51 51 0 1.201 ## 52 52 0 1.045 ## 53 53 0 -1.003 ## 54 54 1 2.348 ## 55 55 0 -0.667 ## 56 56 1 0.606 ## 57 57 1 0.078 ## 58 58 0 -0.122 ## 59 59 0 0.188 ## 60 60 1 0.619 ## 61 61 1 0.475 ## 62 62 2 1.108 ## 63 63 1 0.015 ## 64 64 1 -0.004 ## 65 65 1 -1.161 ## 66 66 0 -0.382 ## 67 67 0 -0.513 ## 68 68 1 3.202 ## 69 69 1 -0.862 ## 70 70 0 0.137 ## 71 71 0 -1.494 ## 72 72 0 -1.470 ## 73 73 0 0.125 ## 74 74 0 -0.997 ## 75 75 0 -0.002 ## 76 76 1 0.072 ## 77 77 0 -0.614 ## 78 78 0 -2.025 ## 79 79 1 -0.725 ## 80 80 0 0.180 ## 81 81 1 1.068 ## 82 82 0 -0.493 ## 83 83 0 0.000 ## 84 84 1 1.623 ## 85 85 1 1.940 ## 86 86 1 -0.597 ## 87 87 0 -0.117 ## 88 88 0 1.201 ## 89 89 0 -0.470 ## 90 90 0 -0.052 ## 91 91 1 0.414 ## 92 92 0 -0.888 ## 93 93 0 -0.445 ## 94 94 2 0.971 ## 95 95 2 0.586 ## 96 96 1 1.613 ## 97 97 0 -0.481 ## 98 98 1 0.067 ## 99 99 1 1.197 ## 100 100 1 -0.556 </code></pre> The study releases the expression data (anonymized) along with disease labels. The genotypes are private. But the eQTL summary statistic, i.e., the regression coefficient $\hat{\beta}$ from $E \sim G$, captures the genotype-expression relationship:
`beta_hat <- coef(lm(E ~ G)) cat("eQTL intercept:", beta_hat[1], "\n") </code></pre>`
`## eQTL intercept: -0.1238356 </code></pre>`
cat("eQTL effect: ", beta_hat[2], "\n") </code></pre> ## eQTL effect: 0.6449423 </code></pre> An attacker with access to these exact eQTL statistics can predict genotypes from expression. Each eQTL is a separate marginal regression, one SNP and one gene: $$ E_k = \beta_k G_k + \epsilon_k, \quad k = 1, \ldots, p $$ With many such pairs, the attacker builds a genotype profile for each individual: for each gene $k$, predict the most likely genotype $G_k$ given the observed expression $E_k$ and the released $\hat{\beta}_k$. More eQTLs = more bits of genotype information = more accurate linking. ## p independent eQTLs, each a marginal regression p <- 20 G_multi <- matrix(rbinom(n p, size = 2, prob = maf), nrow = n, ncol = p) E_multi <- matrix(nrow = n, ncol = p) beta_true_vec <- rep(0.5, p) for (k in 1:p) { E_multi[, k] <- beta_true_vec[k] G_multi[, k] + rnorm(n) } ## attacker runs marginal regression per eQTL betahats <- sapply(1:p, function(k) coef(lm(E_multi[, k] ~ G_multi[, k]))[2]) ## predict genotype per eQTL: round(E / beta) clipped to {0, 1, 2} G_predicted <- sapply(1:p, function(k) { pmin(pmax(round(E_multi[, k] / betahats[k]), 0), 2) }) ## how many genotypes does the attacker get right? accuracy <- mean(apply(G_predicted == G_multi, 2, mean)) cat("Genotype prediction accuracy:", round(accuracy, 3), "\n") </code></pre> ## Genotype prediction accuracy: 0.425 </code></pre> ## linking: match each individual's predicted genotype profile to the reference ## using Hamming distance link_correct <- 0 for (i in 1:n) { dists <- apply(G_multi, 1, function(g_ref) sum(g_ref != G_predicted[i, ])) link_correct <- link_correct + (which.min(dists) == i) } cat("Linking accuracy:", link_correct, "/", n, "\n") </code></pre> ## Linking accuracy: 3 / 100 </code></pre> Each eQTL gives a noisy estimate of one SNP's genotype. But combined across many eQTLs, the predicted genotype profile becomes unique enough to link individuals to a reference database. This exactly how we do identity attack. Protecting eQTL statistics with differential privacy</a></h2> Before releasing $\hat{\beta}$, we can add random noise while calibrating the vulnerability of the data release. For a simple linear regression, changing one person's record shifts $\hat{\beta}$ by at most some $\Delta$. A common option is to add random jitter sampled from a Laplace distribution, as used for the sparse coding, L1-penalty, or LASSO. Namely, we achieve $\epsilon$-DP by adding noise sampled from $\mathsf{Laplace}(0,\Delta/\epsilon)$ set.seed(42) ## sensitivity: how much can one person change a marginal beta? ## for bounded expression in [-B, B]: B <- 3 ## Laplace noise generator rlaplace <- function(n, scale) { u <- runif(n, -0.5, 0.5) -scale sign(u) * log(1 - 2 * abs(u)) } ## conservative sensitivity for a marginal regression coefficient sensitivity <- 2 B / n ## compare linking accuracy under different epsilon values epsilons <- c(0.1, 0.5, 1.0, 10) par(mfrow = c(2, 2)) for (eps in epsilons) { ## add Laplace noise to each marginal beta beta_noisy <- betahats + rlaplace(p, scale = sensitivity / eps) ## attacker predicts genotypes using noisy betas G_pred_noisy <- sapply(1:p, function(k) { pmin(pmax(round(E_multi[, k] / beta_noisy[k]), 0), 2) }) ## linking via Hamming distance linked <- 0 for (i in 1:n) { dists <- apply(G_multi, 1, function(g) sum(g != G_pred_noisy[i, ])) linked <- linked + (which.min(dists) == i) } acc <- mean(G_pred_noisy == G_multi) plot(1:p, betahats, pch = 19, col = "black", ylim = range(beta_noisy, betahats), xlab = "eQTL index", ylab = "Effect size", main = bquote(epsilon == .(eps) ~ " (linked " ~ .(linked) ~ "/" ~ .(n) ~ ")")) points(1:p, beta_noisy, pch = 4, col = "red") legend("topright", legend = c("true", "noisy"), pch = c(19, 4), col = c("black", "red"), cex = 0.8) } </code></pre> At $\epsilon = 0.1$ (strong privacy), the attacker's predictions are essentially random. But the coefficients deviate too much from the true one.</li> At $\epsilon = 10$ (weak privacy), the predictions are nearly as good as the exact case. But the data security breach may occur.</li> The researcher must choose $\epsilon$ based on the acceptable trade-off between utility and privacy</li> </ul> The cost: can we still detect true eQTLs?</a></h2> Privacy is only useful if the released statistics are still scientifically meaningful. Let's simulate a more realistic scenario: many genes, some with real eQTLs and some without, and ask whether we can still distinguish them after adding DP noise. set.seed(123) n_genes <- 200 n_causal <- 20 # true eQTLs ## one SNP, same genotype vector for simplicity G_util <- rbinom(n, size = 2, prob = maf) ## simulate expression for each gene beta_sim <- rep(0, n_genes) beta_sim[1:n_causal] <- rnorm(n_causal, mean = 0.5, sd = 0.2) E_mat <- sapply(beta_sim, function(b) b G_util + rnorm(n)) ## compute true eQTL betas beta_est <- apply(E_mat, 2, function(e) coef(lm(e ~ G_util))[2]) ## add DP noise at different epsilon levels epsilons <- c(0.1, 0.5, 1.0, 10) results <- do.call(rbind, lapply(epsilons, function(eps) { if (is.finite(eps)) { noisy <- beta_est + rlaplace(n_genes, scale = sensitivity / eps) } else { noisy <- beta_est } data.frame( gene = 1:n_genes, beta_noisy = noisy, causal = ifelse(1:n_genes <= n_causal, "true eQTL", "null"), epsilon = if (is.finite(eps)) paste0("epsilon == ", eps) else "epsilon == Inf" ) })) ggplot(results, aes(x = beta_noisy, fill = causal)) + geom_histogram(bins = 40, alpha = 0.7, position = "identity") + facet_wrap(~epsilon, scales = "free_y", labeller = label_parsed) + geom_vline(xintercept = 0, linetype = "dashed", color = "grey40") + labs(x = "Released eQTL effect size", y = "Count", title = "True vs. null eQTLs under differential privacy", fill = "") + theme_minimal() </code></pre> What $\epsilon$ actually promises: It's all about membership attack/protection</a></h2> So what's the theory behind it? We don't want to reveal whether an individual's data point is included by adding some noise with an $\epsilon$-DP guarantee. $$\Pr{\text{stat with Laplace noise} \mid \text{you're in DB}} \le e^{\epsilon} \times \Pr{\text{stat with Laplace noise} \mid \text{you're out of DB}}$$ The intuition is that an attacker who sees the output can't tell whether your data point is in or out since the output would have looked essentially the same regardless. Small $\epsilon$ (say 0.1): the two versions are nearly indistinguishable, so your privacy is strong, but we had to add a lot of noise, so the statistics are less accurate.</li> Large $\epsilon$ (say 10): little noise, accurate statistics, but now your presence visibly shifts the output, so privacy is weak.</li> </ul> We can make this concrete. Suppose before seeing the data the attacker gives some odds that you're in the study. After seeing the released statistics, they update those odds. The $\epsilon$-DP guarantee says this update can only be so large: the odds can grow (or shrink) by at most a factor of $e^{\epsilon}$: $$\frac{\text{odds you're in, after seeing the output}}{\text{odds you're in, before}} \le e^{\epsilon}$$ So $\epsilon$ has a direct reading as a worst-case multiplier on the attacker's betting odds: $\epsilon$</th> $e^{\epsilon}$</th> plain reading</th></tr></thead> 0.1</td> ≈ 1.1</td> odds of guessing membership rise by at most ~10%</td></tr> 0.5</td> ≈ 1.6</td> at most ~1.6×, namely a 50/50 attacker reaches at most ~62%</td></tr> 1.0</td> ≈ 2.7</td> at most ~2.7×, namely a 50/50 attacker reaches at most ~73%</td></tr> </tbody></table> Here is a concrete example. Say the attacker starts at even odds, a coin flip, 50/50 on whether you're in the study. With the exact statistic released, suppose the evidence is strong enough to push their odds up by 100 fold, to 100-to-1. That's a jump from 50% to about 99% confidence that you're in, essentially identified. Now suppose we release the statistic under $\epsilon$-DP instead. The odds can move by at most $e^{\epsilon}$, so: At $\epsilon = 1$: the 100 fold increase of the security breach is capped at $e^1 \approx 2.7\times$, i.e. odds of at most 2.7-to-1, or about 73% confidence, still uncertain.</li> </ul> u- At $\epsilon = 0.5$: capped at $e^{0.5} \approx 1.6\times$, odds of at most 1.6-to-1, or about 62% confidence, only a nudge above the coin flip. At $\epsilon = 0.1$: capped at $e^{0.1} \approx 1.1\times$, odds of at most 1.1-to-1, or about 52%, barely better than the coin flip they started with.</li> </ul> The protocol isn't really making the attack impossible. Instead, we make its guesswork much less likely. American Society of Human Genetics meeting 2025 2025-10-16T00:00:00+00:00 TL;DR: Biased memos from ASHG 2025 talks/posters. AI/ML and new omics technologies always create new interests in the genetics community, but certain aspects of core genetic questions aren't easily changed. Large-scale biobank GWAS data became the norm; then what's next? Can we think of anything beyond polygenic risk predictions? Contextual information matters, but how so? Oct 15</a></h2> Non-linear GWAS methods</a></h3> Refine additive association models with dominant and recessive models</li> Concerns over Type I error calibration</li> How is it different from model checking?</li> </ul> Context-specific polygenic score models</a></h3> Fonseca, Andrew Dahl, U of Chicago</li> Improve prediction performance Can contexts be anything?</li> </ul> </li> Straightforward idea: calculate interaction GWAS</li> Build PGS models base off interaction coefficients</li> Can we come up with summary-based extensions?</li> </ul> How can include gene-gene interactions in MR?</a></h3> No satisfactory answers</li> </ul> Oct 16</a></h2> Squeezing AI models into genetic</li> </ul> Distinguishing causal variants underlying natural selection</a></h3> Arya Rao, Sabeti Lab (Harvard), Steven Riley (Yale), DeepSweep </li> Selective sweep: SLC45A2, causal adaptive variants </li> Distinguish causal and neural-hitchhiker variants </li> Simulated data (coalescence model). Summary statistics. </li> Transformer-based prediction (do we need label data?) </li> GWAS and eQTL enrichment </li> But where is causal guarantee? </li> Question from audience: Have you compared with colocalized or finemapping results?</li> Do ages of variants matter?</li> </ul> </li> Takeaway: Probably useful... run simulation and train black box models. </li> </ul> PubMind</code> language model trained on PUBMED</a></h3> Train LLM on millions of abstracts</li> Existing resources: ClinVar, HGMD, LitVar2, ...</li> inference module: SNV, SV, CNV, Gene fusion</li> Focus on genetic evidence of pathogenecity</li> What are the queries?</li> LLaMa3.3, DeepSeek</li> Standardize disease name</li> https://pubmind.wglab.org</code></a></li> This is still a database... but the benefit of LLM is interactive API</li> We could find more, but how should we use the method?</li> </ul> Spatial 3D reconstruction of H&E data</a></h3> Hu, Yu, Emory</li> Use long range dependency between snapshots</li> Multi-resolution alignment</li> Metamorphism (gradient flow)</li> Cool data: CyCIF + H&E on colorectal cancer</li> Cool 3D image... What did we learn? Can we not find the same cell type identity mapping?</li> </ul> </li> A question from moderator: Any downstream analysis?</li> </ul> Xenotransplantation</a></h3> Eloi Schmauch, Keating, NYU</li> Pig kidney transplantation: pig vs. human cells?</li> Too much technology details? Clustering?</li> Visualization by Celldega</code></a></li> Cell type colocalization over the course of transplantation</li> Questions: How good is the pig genome?</li> Suggestions: Try out different technologies for adjacent slides?</li> </ul> MtDNA-based linage tracing in spatial data</a></h3> Rong, Zhang, UPenn</li> mtDNA far more abundant in cells</li> Can linage information inferred from mtDNA implicate spatial changes?</li> Send cDNA to (1) spatial slide-seq (2) PCR amplification for genotyping</li> Input: MT Loci with cell type deconvolution (spots)</li> Question: Which variant is significantly over-representing cell types?</li> Assign MT variants to cell types</li> Barrett's Esophagus</li> Bracht and Rong et al. (2025)</li> Spot-level localization (testing the purity of nearest neighbors)</li> </ul> Can we simply use MT expression to construct cell lineage? Why not using expression-based CNV method? Spatial neighbour analysis</a></h3> Jun Inamo, U of Colorado</li> Juvenile Idiopathic Arthritis synovium</li> 10x Xenium, n=9</li> co-varying neighborhood analysis (Reshef YA, Nat Biotech 2022</a>)</li> Can we relate this with locally-linear embedding?</li> Show dynamics of many fibroblast marker genes as a function of distance from endothelial cells</li> spatialCoocur</code> (by permutation tests)</li> Compute colocalization scores between two cell types (e.g., T and myeloid, CXCL9-CXCR3)</li> Tertiary Lymphoid Structure (TLS): T-cell and B-cell together</li> </ul> Clonal hematopoiesis, non-coding putative driver mutations</a></h3> Josh Weinstock, Emory</li> How does DNA changes over time?</li> Clonal hematopoiesis of indeterminate potential</li> Clonal extensions lack a mutation in canonical genes.</li> 490k UKBB whole genome</li> frequency of alternate allele as a function of age at blood draw</li> clonal hematopoiesis genes recapitulated!</li> AML risk was also captured in phewas on the variants</li> CH point mutations... how much variance was explained?</li> somatic.emory.edu</a></li> preprint</a></li> Questions from audience: How do we know this is somatic or rare germline?</li> Most of them are germline point mutations.</li> Do cell types matter?</li> Copy number variants?</li> </ul> </li> We may be able to ascertain "age" as a factor for somatic vs. germline... Does it matter?</li> </ul> After the ASHG meeting</a></h2> Sishir Subedi 2025-09-30T00:00:00+00:00 A NEW WEBSITE</a> I am a Ph.D. student in the Department of Bioinformatics. I am interested in developing computational methods to understand the dynamics of cancer genomes. Mainly, I am fascinated by the interplay between the genomic architecture and the adaptive mechanisms of tumor cells in their heterogeneous microenvironment. Previously, I worked as a bioinformatician where I was involved in developing bioinformatics platforms for oncology diagnostics assays. During this time, pondering at clinical NGS data, I appreciated the importance of untangling DNA sequences and how it translates to helping a hopeful cancer patient. Also, I realized that each patient is unique, so is each cell, and I am excited to learn more about them- one cell at a time. I couldn't agree more with what Aristotle famously wrote, "The more you know, the more you realize you don't know." While offline, I am up for anything outdoor. I like visiting new places, hiking hills, and kayaking. STAT 548 PhD Qualifying Papers 2025-08-28T00:00:00+00:00 Introduction</a></h2> I am interested in almost all problems in computational biology and genomics. I expect a student to propose novel statistical approaches that can address challenges in data analysis and modelling of high-dimensional, large-volume biological problems. Feel free to contact me (ypp@stat.ubc.ca</code>). Format</a></h2> You may organize your report including the following sections. Problem definition (1 page): Extract mathematical/statistical problems from the paper and organize them. What are the input data? What is the expected output? </li> Significance (1 paragraph): Why is this an interesting problem? What can be learned by studying this problem? Why is it exciting for you? Author contribution: How did the author(s) find the solution? What was a novel contribution beyond traditional approaches? </li> Limitations/challenges (1 paragraph): What are the assumptions? Are they realistic? What are the technical limitations that the authors acknowledge or not? </li> Novel idea/methods (1-2 pages): Propose your idea and statistical methods. You could interpret the underlying problem in a different formulation. What are related problems/frameworks, but not adopted by the authors? </li> Results (1-2 pages): Include one figure that sketches your approaches. Show tables and figures that clearly demonstrate your methods. </li> Discussion (1 page): Briefly discuss what you have learned and what you would achieve if you were to develop this to a full paper. How would you validate your findings in independent studies, including wet-lab experiments? </li> </ul> Available Papers</a></h2> Smith, L. H., & VanderWeele, T. J. (2019). Bounding bias due to selection. Epidemiology (Cambridge, Mass.), 30(4), 509–516.</a> </li> Brown, B. C., Tokolyi, A., Morris, J. A., Lappalainen, T., & Knowles, D. A. (2025). Large-scale causal discovery using interventional data sheds light on gene network structure in k562 cells. Nature Communications, 16(1), 9628. </a> </li> Deshpande, S. K. (2025). FlexBART: Flexible Bayesian regression trees with categorical predictors. Journal of Computational and Graphical Statistics: A Joint Publication of American Statistical Association, Institute of Mathematical Statistics, Interface Foundation of North America, 34(3), 1117–1126. </li> Marton, S., Lüdtke, S., Bartelt, C., & Stuckenschmidt, H. (2023, October 13). GRANDE: Gradient-Based Decision Tree Ensembles for Tabular Data. The Twelfth International Conference on Learning Representations.</a> </li> Bakhtiari, M., Bonn, S., Theis, F., Zolotareva, O., & Baumbach, J. (2025). FedscGen: privacy-preserving federated batch effect correction of single-cell RNA sequencing data. Genome Biology, 26(1), 216. </li> TAKEN Rautenstrauch, P., & Ohler, U. (2025). Shortcomings of silhouette in single-cell integration benchmarking. Nature Biotechnology, 1–5. </li> TAKEN Zhang, L., Liu, L., Ji, J., Yan, R., Guo, P., Gong, W., Xue, F., Zhou, X., & Yuan, Z. (2025). Efficient Mendelian randomization analysis with self-adaptive determination of sample structure and multiple pleiotropic effects. The American Journal of Human Genetics, 0(0). https://doi.org/10.1016/j.ajhg.2025.06.002 </li> </ol> Excerpts from Ideaflow on creativity - part 3 2025-07-25T00:00:00+00:00 These are quotes and memos from Ideaflow written by Jeremy Utley and Perry Klebahn. It is the final post on Ideaflow. How can we widen the channel of Ideaflow? Flows from many sources will guarantee consistency and much-needed volume. Interdisciplinary collaboration:</a></h3> The beautify of collaboration is that we each have different blind spots. </blockquote> Collaboration doesn't have to be limited to someone with whom I can have a chat in person. It can happen via email and idea exchange in threads. As Jeremy Utley pointed out, it can even happen with AI agents. Since our goal is to generate the widest possible array of perspectives, it's crucial to actively inhibit this instinct, [which is peer pressure preventing them from being opinionated.] </blockquote> Curiosity drives innovation:</a></h3> A good question is specific. </blockquote> So, we may draw our attention to a problem, not a presenter. Even if you start with a provocative question, there’s only so much you can do by looking at a problem from one angle. When ideaflow ebbs, move your frame. </blockquote> The whole point of having a meeting is not to bolster what I believe is true. On the contrary, we need to shake things up. Nearly every good frame begins the same way: How might we...? A good HMW question allows for plenty of exploration while leaving enough structure to keep the discussion focused. </blockquote> Even as we approach the end of a project, we should still raise a question, "how might we...?" Feed the flow:</a></h3> We connect what we have, bringing together two or more elements in a new way. Abundant ideaflow requires enormous amounts of raw material to make more of these unexpected combinations. </blockquote> You don’t innovate in the pickle business by eating cucumbers all day. The more 'distant' the origins of your inputs, the more valuable and interesting the resulting combinations will be. </blockquote> "Wonder Wander":</a></h3> Walking alone is a powerful aid to creativity. </blockquote> “Walking outside produced the most novel and highest quality analogies,” the researchers concluded, but added that moving the body at all “opens up the free flow of ideas.” </blockquote> On a Wonder Wander, volume is key. Write down questions and connections that come to mind and keep walking. </blockquote> The idea is to provoke insights in an uninhibited way. </blockquote> Create space for a difficult problem:</a></h3> Consciously surrendering the difficult problem invited Einstein’s subconscious mind to step in. The next day, Einstein visited Besso again: “Thank you,” he told his friend. “I’ve completely solved the problem.” </blockquote> Daniel Kahneman and Amos Tversky coinvented behavioral economics in large part by taking long walks and joking around with each other. </blockquote> “The secret to doing good research,” Tversky once said, “is always to be a little underemployed. You waste years by not being able to waste hours.” </blockquote> Rather than beat himself [Poincare] up, he spent a few days by the shore. On a morning walk along the bluffs, a solution came to him with “brevity, suddenness and immediate certainty.” Recognizing a roadblock and making the effort to let go of the problem at hand is a counterintuitive but essential creative skill. </blockquote> If we don’t simultaneously carve away less important uses of our time to create space for reflection and contemplation—distance from the problem at hand—we only undermine the effort to boost ideaflow. </blockquote> How to tactically withdraw to create necessary space: Water (e.g., taking a shower or bath) </li> Task switching </li> Turn to a hobby </li> Nap </li> Talk </li> Move </li> </ul> When you withdraw, you are investing in the highest value you bring to your work. If you give yourself permission to step away when you’re stuck, your patience will be rewarded. </blockquote> Creativity is not an option. It's a must in my field. The more I contemplate the sentences that stood out to me, the more obvious they become. The gist of research in my field is far from what leading groups in my own field could have done. I picked up this book with some sense of urgency. If I keep on doing what other people might expect me to do, would I be able to keep my research program? Creativity is a discipline. Let's walk the walk. Excerpts from Ideaflow on creativity - part 2 2025-07-21T00:00:00+00:00 These are quotes and memos from Ideaflow written by Jeremy Utley and Perry Klebahn. TL;DR: Test to choose the right idea out of many. "Even if you're an acknowledged expert in the field, you simply aren't qualified to decide which ideas to pursue in the absence of real-world data. Nobody is!" </li> "Go through stages: test, analyze results, refine, test again." Refining steps make a solid product. My very first implementation of a hidden Markov Model worked just fine. On the benchmark data, the performance was slightly off, compared to the one TA in Machine Learning class showcased, it was just off by a very minute level. The points were taken just because of slight mistakes. If the algorithm were deployed to a real-world example and repeatedly used, the errors could have snowballed into a big mistake. If the test hints at something, it's important to refine and test again. Most people hate this. </li> "Scientists don't win Nobel Prizes for publishing a lack of positive results." A great discovery could be hidden until a sufficient heap of negative results was piled up. </li> "A good test rules out lots of options that won't work while homing in on the ones that might, vastly reducing the risk of failure." </li> "Sir James Dyson tested at least one variation a day for four years, fastidiously documenting the impact of each incremental change to the design." A great product is the result of assiduous fine-tuning. With the help of AI tools, I think such a process can be a lot less painful. </li> "Getting an unbiased outside perspective is too important to skip." This is so difficult not to. </li> "Winnow by excitement" </li> "If you're not curious to find out whether an idea will work, don't try to find out!" </li> </ul> I've struggled so hard to reclaim my right to have "little joy of making one discovery." My reading is definitely biased toward achieving that goal. Just as Jeremy Utley mentioned, AI tools are never an obstacle to constant idea flows. On the contrary, it is a great tool for generating ideas and conducting thought experiments to figure out the right problem, which will keep me awake at night with excitement. Excerpts from Ideaflow on creativity - part 1 2025-05-10T00:00:00+00:00 These are quotes and memos from Ideaflow written by Jeremy Utley and Perry Klebahn. TL;DR: Quantity drives quality. Creativity is learned. "Creativity is doing more than the first thing that comes to your mind." - a seventh-grader in Ohio </li> Creativity isn't a gift reserved for a precious few. It's learned. </li> In other words, it's the capacity to keep generating ideas after the first one that's "good enough." </li> Your brain never makes anything truly from scratch. </li> "Inspiration is a discipline." - Lacrae </li> Breakthrough results don't happen by picking the "right" idea from a list but rather by winnowing ideas down through real-world experimentation and iteration until a clear winner emerges. </li> In most cases, the flow of ideas gets disrupted not by the individuals on a team but by the well-meaning, misguided leader who is so fixated on feasibility and relevance that every divergence from the status quo gets immediately quashed. </li> To be disruptive, the new way of doing something must upend the old one. Designing the transistor, which ultimately makes vacuum tubes obsolete, is a disruptive innovation. To be clear, this is scary--if you happen to be a vacuum-tube manufacturer. </li> Fear is what actually stands in the way of improvement. </li> The riskiest move is no move at all. Conditions are always changing. </li> </ul> Embedding, attention, and transformers 2025-04-10T00:00:00+00:00 We hear lots of excitement about a transformer network. However, I did not find an easy explanation (at least to me) on why and how it works. 1. How should we "embed" (locate) genes (words)? And why do we need embedding?</a></h2> Embedding is like mapping a gene expression value in some latent space well spread between the values. In other words, from a point in the embedding space, we can recover each gene's name and its expression value (or ranking). It's important to understand what the embedding of a gene onto $d$-dimensional space does. Since my mindset has been stuck in a matrix factorization-like model for many years, I find this embedding step quite counter-intuitive. My other complaint is that "embedding" is not a good terminology at all. I think the whole step is more closely related to "stochastic" registration, dictionary, etc (just my complaint). A document (or cell) as a bag of words</a></h3> While not doing explicit embedding, in a typical fully connected neural network layer, we mark (or implicitly embed) a gene to a unique index (a variable or a visible unit). However, gene expression counts are not uniquely handled. In other words, the same gene can be "replicated" as if we keep on adding its contribution to the next hidden layer. The underlying concept seems very similar to Replicated Softmax: an Undirected Topic Model</a>. If we have a softmax model that maps for each word $d=1,\ldots,D$, we would have the energy function (to be minimized): $$E = -\sum_{d=1}^{D} \sum_{g} \sum_{k} W_{gk}^{(d)} v_{g}^{(d)} h_{k}$$ where we denote $g$ for a gene $g$ and $k$ for a factor/topic $k$; $v_{g}^{(d)} = 1$ if and only if a gene $g$ is present in a document $d$; each word $d$'s model have weights $W_{gk}^{(d)}$ between a topic $k$ and gene $g$. Assuming all the models share the same weights, namely $W_{gk}^{(d)}=W_{gk}, \forall d \in [D]$, and $x_{g}$ counts the frequency of a gene $g$ across all the words in this document, we can simplify: $$E = -\sum_{g} \sum_{k} W_{gk} x_{g} h_{k}$$ This makes sense if we take each document as "a bag of words," where each gene is independently drawn within a document, conditioning on the topic proportion and latent topic membership of the word. What have we missed with the bag-of-words assumption?</a></h4> Since bag-of-words models have been so successful for many years and seem good enough for single-cell modelling, it is hard to justify why we would need a more sophisticated language model. So, this is not an exhaustive survey on what the bag-of-words missed: Monotone contribution: We cannot map gene (word) counts in a non-linear fashion, meaning that a higher count of a gene $g$ will always have more weight on the same gene $g$ toward the next layer, i.e., $W_{gk}$ will increase for a particular $h_{k}$. In NLP, double or triple negation expressions are difficult to understand with a bag-of-words model. If a certain housekeeping gene (ribosomal or mitochondrial) is just highly expressed, its influence will dominate and mask out subtle yet significant context-dependent changes. </li> Lack of dependency within a layer: We need many layers to represent word-to-word non-linear relationships. A single layer of the softmax model cannot represent ordering or any sort of positional information across words. Using a softmax function or similar can have some effect of "lateral inhibition"-like mechanisms so that a sparse set of keywords can be selected for each topic. However, it can also make known or unknown stopwords appear across multiple topics without apparent meaning. There is no direct influence across topics that having some high-frequency word in one topic prevents the other topic from picking up the same words, except for some statistical pressure. </li> </ol> Gene embedding is like word embedding in NLP</a></h3> An embedding can bring more flexibility on interpreting the number of the same word (gene) occurred within a document (cell). First of all, what is "word embedding?" In PyTorch</code>'s documentation for Embedding</code></a> class, it says: Input: (), IntTensor or LongTensor of arbitrary shape containing the indices to extract</code> and</li> Output: (,H), where * is the input shape and H = embedding_dim</code></li> </ul> For scGPT</code></a>, this embedding class is directly applied to integer-valued expression tensor. In a nutshell, each cell is like a fixed length (number of words/tokens) document (or sentence) with tokens (genes). Discretized/binned gene expression levels on each gene (word) uniquely mapped onto some $d$-dimensional latent space. For instance, $x_{g}=1$ and $x_{g}=10$ are independently located apart from each other in some embedding space. Upon training, we may be able to tell whether $x_{g}=1$ and $x_{g}=10$ mean semantically similar but very different from $x_{g}=5$. A gene embedding step clearly brings non-monotone expressivity in a language model. Anecdotally, "not" and "not, not, not" can be located in nearby locations apart from "not, not" if that's easier to picture. We do not specify their positions but learn semantic similarity from examples. Okay, what about gene/word identity? Many single-cell foundation models include a positional embedding. If only the highly variable/expressed genes were considered, and subsequent gene-by-gene attention modules could have been mapped out in GPU memory, we wouldn't need to consider this. As long as we can feed genes in the same order, we are uniquely mapping gene symbols onto identity space. A token of the g-th position is always the same gene. Positional embedding would have been redundant. However, if the length of accessible genes from one cell to another greatly varies, and padding unexpressed genes/words with zeros is highly inefficient, we need to think of a cell as a series of tokens, just like a sentence of words dynamically rearranged. Here, a token $\neq$ a word. 2. What attention mechanisms can do with gene embedding vectors?</a></h2> Now, let's discuss the elephant in the room: What in the world is attention? Why is attention all you need?</a> First, we need to know what else has been considered so essential in previous work besides attention mechanisms. In NLP, a sentence of words has been modelled as a sequence of hidden states, such as recurrent neural network (RNN) models. One of the most recent RNN architectures is Long Short-Term Memory (LSTM), where the layer keeps the memory of previous words' hidden states to figure out the state of a current word. Since "memory" is not a part of the model parameters, it is somehow wasted, computationally cumbersome and not ideal for dealing with long sentences. How can we represent a sentence of words in a meaningful way?</a></h3> Okay. We need to model word-to-word dependency structures without inferring what each word's state could be. Why did we want to assign some latent states to words in a sentence? From the perspective of a hidden Markov model (a grandfather of a language model), knowing previous states, we can identify blocks of consecutive words, and then we can aggregate information within each block without too much loss of information. Keeping memory in the model is useful to parse out the meaning of a sentence. It might be an effective way to explain what the meaning of an overall sentence would be for a non-native English reader/writer, segmenting a sentence into multiple blocks (clauses), and so on. At the same time, such a pedagogy could become an obstacle when students want to generate a new sentence by themselves. An attention mechanism simplifies this sentence-generating process.</a></h3> It is like the way a baby learns a new language to me. A stream of words comes as a sentence. The first word will be understood together with other words that arrive later in the same sentence. It might be the imminent one or much later. As a neural model accumulates examples of sentences, associating the first word with others, the second with others, and so on, we can follow the probability of such word-to-word associations. Provided that the order of genes within a cell sentence/document is fixed, we can keep track of how frequently genes are jointly associated (or disassociated). For each gene, we can fit a softmax function (probability) to learn about frequent "neighbours" in the embedded space. Well, very roughly. For this type of process, there is no memory of word-specific states. In other words, if we consider a single cell as a scene of a movie or play, we don't need to know whether a gene plays an antagonist or protagonist in the scene. We just need to know who else appears here and how they interact with one another. Hashing, dictionary, and (query, key, and value)</a></h3> Interactions between genes matter, and keeping track of interactions is our only interest. For simplicity, although we don't need to do so, let's assume that gene embedding coordinates are fixed. We need to bring two genes closer if they talk to each other more frequently than others. Borrowing the concept of locality sensitive hashing</a>, we use a query of one gene to locate neighbouring genes according to their key in the embedding space. For closer gene pairs key and query should be similar, while distant pairs have largely independent key and query patterns. We could revise the embedding coordinates directly based on semantic similarity scores. However, if we want to stack up this gene-to-gene matching process within many different contexts, defining additional transformation functions for key and query vectors can render a more general modelling capability. More specifically, for the embedding vector of a gene $g_{1}$, namely $x_{g_{1}}$, and other $g_{2}$, namely $x_{g_{2}}$, we transform them into $[Q_{g_{1}1}, \ldots, Q_{g_{1}d}]$ and $[K_{g_{2}1},\ldots,K_{g_{2}d}]$, $d$-vectors. Then, we test how well the query and key vectors agree with each other by taking the dot product; this proximity information can be scaled and converted to a probably by the softmax transformation. $$\textsf{Attention}(Q,K,V) = \underbrace{ \textsf{softmax}\left( \frac{QK^{\top}}{\sqrt{d}} \right) }_{\color{blue}\textsf{proximity between genes}} V.$$ The value vectors, $V$, map how each gene can work. Roughly, we may think of a traditional fully connected neural network layer as a diagonal attention matrix and only the V matrix. There, we encode gene expression patterns into some hidden codes. Putting them all together, in terms of the scene analogy, genes act on some value $V$ here, while interacting with their partners and influencing one another. In the next scene, we see the consolidated net effects of the entire scene. Sure, we may rule out some genes/actors in the current scene as they fit in the story. Moreover, instead of considering all pairwise interactions, we may restrict/focus on feasible interactions, potentially set by physical constraints, such as protein-protein interactions and cellular locations. 3. Can attention-based models be foundational for single-cell RNA-seq modelling?</a></h2> I think gene embedding is a good idea. Well, I also like the way we represent gene-gene interactions as a stack of attention layers. But could we derive a foundational model? For a model to be foundational, well, what features, characteristics, or properties should a model possess and demonstrate? Very roughly... If we roughly consider each transformer layer modelled by attention networks as a scene or a set of rules for generating a scene, we may be able to run a realistic show with fifty of them stacked together, so is true with diffusion-based models. Well, is it foundational? We may need to define what the "foundation" model exactly means, though. I will leave all these sorts of questions for other posts. Transparent parameterization for a partially identified model 2025-04-07T00:00:00+00:00 TL;DR: Consider transparent parameterization for a partially identified model if a strong dependency with prior distributions was found in Bayesian causal inference. Transparent Parametrizations of Models for Potential Outcomes</a></li> </ul> Set-up: $Z \in {0,1}$: a randomized assignment or "instrument"; $Z$ is a treatment, presumed to be randomized, e.g., the assigned treatment.</li> $X \in {0,1}$: an exposure measure; $X$ is an exposure subsequent to treatment assignment.</li> $Y \in {0,1}$: a final response; $Y$ is the response.</li> There might be confounding on $X$ and $Y$</li> </ul> Goal: What is the effect of $X$ on $Y$? Identifiability matters, so we should make which parts are identifiable and which parts not.</a></h2> It is often argued that identifiability is of secondary importance in a Bayesian analysis provided that the prior and likelihood lead to a proper joint posterior for all the parameters in the model. [...] we argue that partially identified models should be re-parameterized so that the complete parameter vector may be divided into point-identified and entirely non-identified subvectors. Such an approach facilitates "transparency", allowing a reader to see clearly which parts of the analysis have been informed by the data </blockquote> A potential outcome problem</a></h3> We further define: $X^{(z=0)}$ or $X^{(z=1)}$ for the treatment a patient would receive if assigned to $Z=z$</li> $Y^{(x,z)}$ for the outcome for a given patient if they were to be assigned to $Z = z$, and then were exposed to $X = x$.</li> </ul> Why is it a potential outcome problem? A person may or may not take a pill (exposure) in relations with the assignment $Z$: compliance type</th> $X^{(z=0)}$</th> $X^{(z=1)}$</th></tr></thead> Never Taker (NT)</td> 0</td> 0</td></tr> Complier (CO)</td> 0</td> 1</td></tr> Defier (DE)</td> 1</td> 0</td></tr> Always Taker (AT)</td> 1</td> 1</td></tr> </tbody></table> What are the potential responses regarding the exposure and response variables? response type</th> $Y^{(x=0)}$</th> $Y^{(x=1)}$</th></tr></thead> Never Recover (NR)</td> 0</td> 0</td></tr> Helped (HE)</td> 0</td> 1</td></tr> Hurt (HU)</td> 1</td> 0</td></tr> Always Recover (AR)</td> 1</td> 1</td></tr> </tbody></table> Pearl (2000) considering a subset of data</a></h4> If we only consider "possible" (or partially factual) cases, $Z=0$ implies $X=0$, which leaves us only NT and CO. We can assume $Z \to X \to Y$ but no direct path $Z \to Y$ that skips over $X$. $$\textsf{ACE}(X \to Y) \equiv \mathbb{E}\left[Y^{(x=1,z)} - Y^{(x=0,z)}\right]$$ Pearl proposes analyzing the model by placing a prior distribution over $p(tX , tY)$ and then using Gibbs sampling to sample from the resulting posterior distribution for ACE($X \to Y$). </blockquote> Sensitivity to the choice of prior distribution over compliance and response types: If the model were identified, we would expect such a change in the prior to have little effect (the smallest observed count is 12). However, as the plot shows, this perturbation makes a considerable difference to the posterior. </blockquote> $$\textsf{ACE}(X \to Y) = \pi_{\textsf{CO}}\left(\gamma_{\textsf{CO}}^{(1,\cdot)} - \gamma_{\textsf{CO}}^{(0,\cdot)}\right) + \pi_{\textsf{NT}}\left(\gamma_{\textsf{NT}}^{(1,\cdot)} - \gamma_{\textsf{NT}}^{(0,\cdot)}\right)$$ Why is it not fully identified? The following factor graph (or PGM-like) well summarizes. We need to calculate two types of differences: $\gamma_{\textsf{CO}}^{(1,\cdot)}$ vs. $\gamma_{\textsf{CO}}^{(1,\cdot)}$</li> $\gamma_{\textsf{NT}}^{(1,\cdot)}$ vs. $\gamma_{\textsf{NT}}^{(1,\cdot)}$</li> </ul> Here, $\gamma_{\textsf{NT}}^{(1,\cdot)}$ has to rely on Bayesian sampling over some prior distribution. Okay, we could average out unidentified random variables, but if they were not attached to any part of the data, we would resort to the shape of the prior distribution. From a saturated model to the model considered by Hirano et al. (2000)</a></h4> We need to consider more data to be able to identify all the parameters: Or at least, we should make each parameter attached to measurable sufficient statistics with a bit more assumptions: Monotonicity of $X$ of $Z$: There are no defiers, meaning $X^{(0)} \le X^{(1)}$. </li> Stochastic exclusion of NT: $\gamma_{NT}^{(0,1)} = \gamma_{NT}^{(0,0)}$ </li> Stochastic exclusion of AT: $\gamma_{AT}^{(1,1)} = \gamma_{AT}^{(1,0)}$ </li> </ul> Intent-to-treat effect (Hirano et al. 2000): $$\textsf{ITT}(\textsf{compliance type}) = \mathbb{E}\left[ Y^{(X(1),1)} - Y^{(X(0),0)} | \textsf{compliance type}\right].$$ We can parameterize this by: $\textsf{ITT}(\textsf{CO}) = \gamma_{CO}^{(x=1,z=1)} - \gamma_{CO}^{(x=0,z=0)}$ </li> $\textsf{ITT}(\textsf{NT}) = \gamma_{NT}^{(x=0,z=1)} - \gamma_{NT}^{(x=0,z=0)} = 0$ (by the above exclusion assumption 2) </li> $\textsf{ITT}(\textsf{AT}) = \gamma_{AT}^{(x=1,z=1)} - \gamma_{AT}^{(x=1,z=0)} = 0$ (by the above exclusion assumption 3) </li> </ul> My takeaway</a></h2> I didn't go over all the inequalities laid out by the authors, so there is more clarity to be found. </li> Working out some sort of factor graph representation will be useful to see which parameters are identifiable. </li> I wonder if we can design an algorithm that detects partially or fully identified model parameters. </li> A Bayesian inference engine geared toward "integrate out" will not always yield identifiable results. </li> </ul> A stochastic block model and decades later 2025-03-30T00:00:00+00:00 Generative models</a></h2> Discrete (non-mixed) stochastic block models</a></h3> One of the most well-known examples in early social network study is this: An Information Flow Model for Conflict and Fission in Small Groups</a> Zachary (1977). At the beginning of the study there was an incipient conflict between the club president, John A., and Mr. Hi over the price of karate lessons. Mr. Hi, who wished to raise prices, claimed the authority to set his own lesson fees, since he was the instructor. John A., who wished to stabilize prices, claimed the </blockquote> Letting information flow from the president (John A, or vertex 1</code>) to a new faction leader (Mr. Hi, or vertex 34</code>), the min-cut, max-flow algorithm results in two groups in the social network data: In classical stochastic block modelling, we assign a vertex (node) to a certain group or block (cluster). Each vertex may belong to only one of the blocks. An algorithm for clustering relational data with applications to social network analysis and comparison with multidimensional scaling</a> (1975) </li> Learning systems of concepts with an infinite relational model</a> (2006) or use this link</a> </li> </ul> We can define a generative model for a binary relationship between vertices $i$ and $j$: $$Y_{ij} | B \cdot \sim \prod_{k,l} B_{kl}^{z_{i}^{(k)} z_{j}^{(l)}}$$ where $z_{i}^{(k)} = 1$ if and only if a vertex $i$ belongs to a cluster $k$; otherwise, $z_{i}^{(k)} = 0$. Here, $z_{i}^{k}$ can be modelled by Chinese Restaurant Process</a> a priori to accommodate an arbitrary number of groups/blocks/clusters. Mixed-membership stochastic block models</a></h3> Modelling all the group structures by a block matrix $B$ (Bernoulli or Poisson rate) is somewhat restrictive and would require unnecessarily many groups to model detailed structures. Instead of giving vertices discrete membership assignments, we can relax that their membership can be partial or more continuous in some latent space. Latent Space Approaches to Social Network Analysis</a> (2002)</li> </ul> $$p(Y) = \sigma\left( \beta_{0} + \sum_{a} \beta_{a} x_{ija} + \left. \sum_{k} z_{ik} z_{jk} \middle/ \sqrt{ \sum_{k} z_{jk}^{2} } \right. \right)$$ where $x_{ija}$ denotes an attribute variable for a pair $(i,j)$ and $\sigma(x)=1/(1+\exp(-x))$ for brevity and $z_{ik} \in \mathbb{R}$. Similarly, we can give a partial membership of a vertex by assigning discrete membership of each vertex's incoming and outgoing interactions (edges/pairs adjacent to a vertex). Mixed Membership Stochastic Blockmodels</a> (2008)</li> </ul> $$p(Y_{ij} =1 |\cdot) = \sum_{k,l} z_{i\to{j}}^{(k)} B_{kl} z_{j\to{i}}^{(l)}$$ where we assign a discrete membership on the (directed) edges (all pairs). Here, each edge can belong to only one of the cluster, namely, $z_{i \to {j}}^{(k)} \in (0,1)$ and $\sum_{k} z_{i \to{j}}^{(k)} = 1$, but a vertex can be adjacent to multiple edge clusters without constraints. Nonparametric Latent Feature Models for Link Prediction</a> (2009)</li> </ul> We can relax the latent membership matrix $z_{ik}$ to allow mixed membership on nodes more explicitly: $$p(Y_{ij} = 1 | \cdot) = \sigma\left( \sum_{k,l} z_{ik}^{\top} W_{kl} z_{jl} \right)$$ where $z_{ik} \in {0,1}$ without any row and column constraints. Efficient and principled method for detecting communities in networks</a> (2011)</li> </ul> Another similar formulation can be found here with Poisson: $$p(Y_{ij}| \cdot) = \frac{\left(\sum_{k} \theta_{ik} \theta_{jk}\right)^{Y_{ij}}}{Y_{ij}!} \exp\left( \sum_{k} \theta_{ik} \theta_{jk} \right),$$ where $\theta > 0$. Here, we can partially introduce latent variables for edges $$\log p(Y_{ij}| \cdot) \ge \sum_{k} \left[ Y_{ij} z_{ij}^{(k)} \log\frac{\theta_{ik}\theta_{jk}}{z_{ij}^{(k)}} - \theta_{ik} \theta_{jk} \right] - \log Y_{ij}$$ where the bound was derived by Jansen's inequality. Let's put a stochastic block model onto a VAE framework</a></h2> Stochastic Blockmodels meet Graph Neural Networks</a> (2019) It is therefore desirable to have a model that retains the basic spirit to OSBM/LFRM</a> (e.g., easy interpretability and strong link prediction performance), but with greater expressiveness, and a simpler and scalable inference procedure. </blockquote> The decoder model (the left) is almost the same as the latent feature relational model: $$p(A_{ij}| \cdot) = \sigma(f_{i}(z_{i})^{\top}f_{j}(z_{j}))$$ where $f_{i}$ is a neural network-like model for a vertex $i$ that maps a latent state $z_{i}$ to latent feature vectors. I think the encoder part is more interesting: $q(v_{ik}) = \textsf{Beta}(c_{ik}, d_{ik})$ </li> $q(b_{ik}) = \textsf{Bernoulli}(c_{ik}, d_{ik})$ </li> $q(r_{i}) = \mathcal{N}(\mu_{i}, \textsf{diag}(\sigma_{i}^{2}))$ </li> </ul> All these "variational" parameters are modelled by a neural network model mapping a full adjacency matrix $A$ and a feature matrix $X$. We can make sampling steps differentiable by approximate reparameterization tricks. ML/AI cheat sheet: Encoder 2025-03-27T00:00:00+00:00 Single-cell (pseudotime) trajectory analysis 2025-03-24T00:00:00+00:00 TL;DR: The goal is to give a "pseudotime" for each cell. Use minimum spanning tree (MST) to estimate "independent" lineages. Sort the cells within each lineage to assign "time" information. Pseudotime trajectory inference</a></h3> monocle</code> The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells</a></a></h4> monocle</code> deals with pairs of cells directly:</li> </ul> We developed Monocle to informatically order the cells by their progress through differentiation rather than by the time they were collected, maximizing the transcriptional similarity between successive pairs of cells. </blockquote> Independent Component Analysis to reduce dimensions:</li> </ul> It reduces the dimensionality of this space using independent component analysis17. Dimensionality reduction transforms the cell data from a high-dimensional space into a low-dimensional one that preserves essential relationships between cell populations but is much easier to visualize and interpret. </blockquote> Fit a minimum spanning tree on the cells and identify the longest path:</li> </ul> Monocle constructs a minimum spanning tree (MST) on the cells [...] The algorithm finds the longest path through the MST, corresponding to the longest sequence of transcriptionally similar cells. </blockquote> Assign "pseudo time" to each cell along the path:</li> </ul> Finally, Monocle uses this sequence to produce a 'trajectory' of an individual cell's progress through differentiation. </blockquote> TSCAN</code> TSCAN: Pseudo-time reconstruction and evaluation in single-cell RNA-seq analysis</a></a></h4> Gene clustering to reduce "drop out" genes:</li> </ul> Before pseudo-time reconstruction, [...] in order to alleviate the effect of drop-out events on the subsequent analyses, genes with similar expression patterns are grouped into clusters by hierarchical clustering (using Euclidean distance and complete linkage). [...] For each cluster and each cell, the expression measurements of all genes in the cluster are averaged to produce a cluster-level expression which will be used for subsequent MST construction. </blockquote> Then, they performed PCA on the clustered matrix.</li> </ul> After gene clustering, single-cell transcriptome for cell $i$ becomes a $H$ dimensional vector $E_{i}$ Here, $H$ is the number of gene clusters. $E_{i}$ still has high dimension, and many components in this vector are still correlated. The dimensionality makes visualization and statistical modeling difficult. For this reason, TSCAN further reduces the dimension of using principal component analysis (PCA). [...] After PCA, the H dimensional vector is mapped to a lower dimensional space and becomes a $K$ dimensional vector $\tilde{E}_{i}$. Here, $K$ is much smaller than $H$. </blockquote> Cell clustering by fitting a mixture of multivariate Gaussian</li> </ul> The clustering is performed using the mclust (22) package in R which fits a mixture of multivariate normal distributions to the data $\tilde{E}_{i}$. </blockquote> Minimum spanning tree</li> </ul> Next, TSCAN constructs a minimum spanning tree to connect all cluster centers. In a connected and undirected graph, a spanning tree is a subgraph that is a tree and connects all the vertices (or 'nodes'). [...] Unlike the MST approach used by Monocle where the tree is constructed to connect individual cells, the MST in TSCAN is constructed to connect clusters of cells. </blockquote> Identify "main" branch and pick the root randomly</li> </ul> A tree may have multiple branches. By default, we define the main path of the tree (solid lines in Figure 1D) as the path with the largest number of clusters. If more than one path has the same largest number of clusters, the path with the largest number of cells becomes the main path. The main path has two ends. Without other information, one end will be randomly picked up as the origin of the path. </blockquote> Cell ordering and pseudo-time calculation; we can learn an average direction between clusters $i$ and $j$ and apply each cell in the cluster $i$ (the preceding one) with respect to this vectorial information:</li> </ul> All cells in these clusters will be ordered along the path as follows. Let $C_{i}$ ($i$ = 1, 2, $\ldots$, $M$) indicate the ordered clusters, where $M$ is the number of clusters on the ordered path. </blockquote> Sort all the cells by these vectorial information:</li> </ul> Cell orderings are determined in three steps. First, for cells which are in the same cluster and are projected onto the same edge, their order is determined by the projected values on the edge. Second, within each cluster, the order of cells projected onto different edges is determined by the order of edges, which is given by the cluster-level ordering. Third, the order of cells in different clusters is determined by the order of clusters. In this way, all cells can be placed in order. </blockquote> slingshot</code> Slingshot: cell lineage and pseudotime inference for single-cell transcriptomics</a></a></h4> Estimate cluster-level MST</li> </ul> Slingshot identifies lineages by treating clusters of cells as nodes in a graph and drawing a minimum spanning tree (MST) between the nodes, [...] We have found that a Mahalanobis-like distance, i.e., a covariance-scaled Euclidean distance, that accounts for cluster shape, works well in practice, but users have the option of specifying any type of distance measure (e.g., Euclidean, Manhattan). </blockquote> Estimate pseudotime by principal curves</li> </ul> The second stage of Slingshot is concerned with assigning pseudotimes to individual cells. For this purpose, we make use of principal curves to draw a path through the gene expression space of each lineage </blockquote> Project all data points onto the curve and calculate the arc length from the beginning of the curve to each point’s projection. Setting the lowest value to zero, this produces pseudotimes. </li> For each dimension $j$, $j \in {1,…,J'}$, use the cells' pseudotimes to predict their coordinates, typically with a smoothing spline. This produces a set of J' functions which collectively map pseudotime values, thereby defining a smooth curve in J' dimensions. </li> Repeat this process until convergence. We use the sum of squared distances between cells’ actual coordinates and their projections on the curves to determine convergence. </li> </ol> monocle3</code> The single-cell transcriptional landscape of mammalian organogenesis</a></a></h4> I think monocle3</code> made a wrong decision. Okay, the original monocle</code> might have resulted in a noisy MST backbone, or the method is infeasible for a million cells. However, there is no rationale behind "projecting data onto UMAP and performing clustering," and the trajectories monocle3</code> report report never looks like a trajectory (because UMAP internally optimizes for cluster-like patterns).</li> </ul> Monocle3</code> first projects cells onto a low-dimensional space encoding transcriptional state using UMAP. It then groups mutually similar cells using the Louvain community detection algorithm, and merges adjacent groups into ‘supergroups’. Finally, it resolves the paths or trajectories that individual cells can take during development, identifying the locations of branches and convergences within each supergroup. Finally, it resolves the paths or trajectories that individual cells can take during development, identifying the locations of branches and convergences within each supergroup. </blockquote> Downstream analysis to statistically test "dynamic" gene programs</a></h3> Trajectory inference across multiple conditions with condiments</a> </li> Enabling single-cell trajectory network enrichment</a> </li> A statistical framework for differential pseudotime analysis with multiple single-cell RNA-seq samples</a> </li> </ul> Brain Age Gap GWAS 2025-03-22T00:00:00+00:00 Genetically supported targets and drug repurposing for brain aging: A systematic study in the UK Biobank</a> Brain age gap (BAG), the difference between chronological and biological/estimated ages, is an interesting observation to think about. $\Delta$ age sounds more natural to me. A colleague and collaborator in my previous lab also attempted to estimate biological aging using molecular profiles of postmortem data: Rate of brain aging and APOE ε4 are synergistic risk factors for Alzheimer’s disease</a></li> </ul> As noted in the paper, several related works have been done to measure the difference between biological and chronological ages in multiple organs. Heterogeneous aging across multiple organ systems and prediction of chronic disease and mortality</a></li> </ul> Although existing studies emphasize the polygenic architecture of brain aging, genetically supported drug target studies of brain aging, as a tool to anticipate the effect of drug action on brain aging, are incredibly lacking [regarding the above]. </blockquote> Brain imaging data confer an invaluable window to look into how old (biologically) person's cognitive ability is. Train $f$: image $\to$ age (and $\Delta$ age)</li> Some genetic variant $\to$ $\Delta$ age</li> The same genetic variant $\to$ gene expressions (GTEx v8)</li> Sort out druggable targets overlapping TWAS genes (eQTL and GWAS co-localized) with drug-gene interaction networks</li> Mendelian Randomization to build up a causal model: genetic variant $\to$ druggable protein/gene $\to$ $\Delta$ brain age.</li> </ol> It's great to know that 3d-vision transformer works great: A promising finding in our study was the clear advantage of using a state-of-the-art deep learning model, namely, 3D-ViT, for brain age estimation. [...] We performed early stopping to prevent potential overfitting problems if the validation error did not improve within 10 epochs. </blockquote> However, it also raises my question: What if an optimal predictor can perfectly estimate a person's age based on image or molecular profiles? Could we ever see $\Delta$ age? How do we know that enough is enough? Finding genetic underpinnings of the druggable disease mechanisms is probably one of the most important contributions of this work (as we can get around potential confounding issues). But can we really? To accurately identify druggable targets for brain aging, it is crucial to mitigate the confounding impact of brain disorders on brain aging. </blockquote> Potentially, I think ageing processes are naturally heteroscedastic, so measuring $\Delta$, not taking into account over/under-dispersion, seems like a serious issue. Finding a cancer-related gene can be both interesting and worrying. In addition to confirming the associations between well-known genes like RUNX2, CRHR1, and INPP5A with BAG, our study has also identified genes, including TP53, that are associated with BAG in the context of brain aging. </blockquote> What if the whole analysis is just confounded by a person's age? Even so, the study has some value since we characterized what would have happened during ageing and potential drug targets and repurposing, etc. Notes on 'How altered cancer metabolism influences cancer progression' 2025-03-17T00:00:00+00:00 Matthew Vander Heiden</a> gave this talk at the BC Cancer Seminar. This is just my note jotted down during the seminar. TL;DR Potentially interesting previous work is already out. Cancer tissue of origin constrains the growth and metabolism of metastases</a> </li> A CRISPRi/a screening platform to study cellular nutrient transport in diverse microenvironments</a> </li> </ul> How do cancer cells adapt to make metabolic alterations? What factors affect how cancer cell metabolism changes? (1) lineage, (2) genetic Mutation, and (3) tissue environment. The cells definitely need to express many enzymes to overcome thermodynamic barriers. Well, nutrients found in normal tissues can tell us about those found in tumours. It makes sense to think that cancer needs to adapt to the metabolic environment (shaped by resident tissues/cells) rather than changing it. A "seed and soil" model</a></h4> Cancer tissue of origin constrains the growth and metabolism of metastases</a> </li> PDAC tumour cells prefer to grow in environments similar to the tissue of origin. </li> How much effect can be attributed to the soil? How much variance can be explained by metabolite changes? </li> DepMap</code> project also profiles cancer metabolism MetMap</code> A metastasis map of human cancer cell lines</a> </li> Preliminary work reveal the complexity of underlying regulatory mechanisms. Specific nutrient deficiencies not necessarily predict metastatic tumour growth. They tested their hypothesis by measuring growth of auxotrophs. </li> Is this an opportunity to study how cells and genes respond to differential metabolic environments? </li> </ul> Energy consumption vs. budget vs. cellular growth</a></h4> Nutrients are both the source of ATP and budget (protein, nucleic acids, lipids, etc.) and eventually NAD+ regeneration. Mitosis aligned with oxygen consumption, NAD+ regeneration, and less protein synthesis, where cell growth was reduced. Respiration-deficient cells grow slowly, why? Cells that underwent unsuccessful mitosis die more frequently. </li> Aurora kinase can shed lights into mechanisms as cellular phenotype are similar. Aneuploidy confirmed by single-cell CNV profiles. </li> Basically Cancer is a "reduced" state. These results emphasize that somatic mutation is necessary but not sufficient to drive cancer. </li> Nucleotide imbalance leads to stress and lineage-dependent changes, and eventually, gene expression program changes (epigenetic states and enzymatic activities to reach a steady state?). </li> </ul> Immune and metabolism</a></h4> Adenosine nucleotide metabolic pathways can influence immune response (suppression). </li> Interesting paper: The adenosine pathway in immuno-oncology</a> </li> Adenosine uptake blockage rescued CD8+ T-cells proliferation. Blocking A2A receptor might have a therapeutic value. </li> Synergy with immune checkpont blockade therapy. </li> </ul> `.R/Makevars` on silicon Mac 2025-01-22T00:00:00+00:00 TL;DR: include OBJCXX=clang++ -std=gnu++17</code> in .R/Makevars</code> to help install.pacakges</code> to help Object-C code compilation. Error messages</a></h3> I wanted to install tidyverse</code> on a new Mac. Almost everything worked out of the box, except for one thing: R</code>'s systemfonts</code> triggered an hour of pain. clang++</code> stopped parsing when it encountered mac/FontManagerMac.mm</code> while complaining that it was unable to understand c++11</code> features, e.g., constexpr</code>. clang++ -I"/opt/homebrew/Cellar/r/4.4.2_2/lib/R/include" -DNDEBUG -I/opt/homebrew/opt/freetype/include/freetype2 -I/opt/homebrew/opt/libpng/include/libpng16 -I'/opt/homebrew/lib/R/4.4/site-library/cpp11/include' -I/opt/homebrew/opt/gettext/include -I/opt/homebrew/opt/readline/include -I/opt/homebrew/opt/xz/include -I/opt/homebrew/include -fobjc-arc -fPIC -c mac/FontManagerMac.mm -o mac/FontManagerMac.o In file included from mac/FontManagerMac.mm:5: In file included from mac/../FontDescriptor.h:11: In file included from mac/../utils.h:9: In file included from /opt/homebrew/lib/R/4.4/site-library/cpp11/include/cpp11/protect.hpp:10: /opt/homebrew/lib/R/4.4/site-library/cpp11/include/cpp11/R.hpp:52:1: error: unknown type name 'constexpr' 52 | constexpr R_xlen_t operator"" _xl(unsigned long long int value) { return value; } | ^ /opt/homebrew/lib/R/4.4/site-library/cpp11/include/cpp11/R.hpp:52:19: error: expected ';' after top level declarator 52 | constexpr R_xlen_t operator"" _xl(unsigned long long int value) { return value; } | ^ | ; </code></pre> Solution: make your Clang++</code> understand new features</a></h3> I believe CRAN packages are generally well prepared to handle C/C++ compiling, except for Object-C! We need to include the following line to .R/Makevars</code> for a silicon Mac: OBJCXX=clang++ -std=gnu++17 </code></pre> We could additionally tweak so that compilation processes can take place on many files with an additional flag: MAKEFLAGS = -j10 </code></pre> Here, we run 10 jobs simultaneously. Another issue on openmp</code></a></h3> In order to compile Rcpp</code> codes with <omp.h></code>, we also need libomp</code> for multi-threading. We can add more lines to .R/Makevars</code> and here is an example: LDFLAGS += -L/opt/homebrew/opt/llvm/lib LDFLAGS += -L/opt/homebrew/opt/libomp/lib -lomp CXXFLAGS += -I/opt/homebrew/opt/llvm/include CXXFLAGS += -I/opt/homebrew/opt/libomp/include -Xclang -fopenmp CPPFLAGS += -I/opt/homebrew/opt/llvm/include CPPFLAGS += -I/opt/homebrew/opt/libomp/include -Xclang -fopenmp PKG_LDFLAGS += $(LDFLAGS) PKG_CXXFLAGS += $(CXXFLAGS) PKG_CPPFLAGS += $(CPPFLAGS) OBJCXX = clang++ -std=gnu++17 MAKEFLAGS = -j10 </code></pre> Alcohol can cause cancers? 2025-01-03T00:00:00+00:00 Perhaps based on this prospective study</a> and many other articles, the US Surgeon General identified heavy alcohol consumption as a potential cause of many cancers</a>. There are several biological mechanisms. See Four ways alcohol can cause cancer</a>. All of them remind me of cancer biology courses I took a decade ago. I agree alcohol can cause cancer in some ways, but proving the causality and measuring causal effect sizes will probably be another research question. Showing tobacco causes cancer took nearly a half-century from the earliest wakeup calls until the US Surgeon General finally acknowledged the fact, and perhaps longer until the tobacco companies took meaningful action. I am wondering if there are any (genetic) instrumental variables to establish some sense of causality. Well... In fact, there are some: Genetic architecture of alcohol consumption identified by a genotype-stratified GWAS and impact on esophageal cancer risk in Japanese people</a> No wonder several alcohol intake GWAS variants are located in biochemically-relevant genes, including ALDH1B and ALDH2; they are also "colocalized" in interaction analysis with cancer GWAS studies.</li> </ul> </li> Beyond GWAS of Colorectal Cancer: Evidence of Interaction with Alcohol Consumption and Putative Causal Variant for the 10q24.2 Region</a> "Another possibility is that light-to-moderate drinking has a protective effect on risk of colorectal cancer, even though heavier consumption is detrimental." This clearly sets their questions and why they designed the study to show in which direction effects could go. </li> Also, G x E interaction analysis with stratification--non-drinking, light-drinking, and heavy-drinking. We see some significant hits (why?), and they end up with somewhat bizarre conclusion: "Our study suggests that the association with colorectal cancer in 10q24.2/COX15 observed in genome-wide association study is strongest in nondrinkers." What about the other variants? </li> What if drinking behaviour is a mediator variable (genetic variant to drinking to cancer)? The stratification step could create artificial associations between genetic variants and cancer phenotypes. </li> </ul> </li> The causal effects of genetically predicted alcohol consumption on endometrial cancer risk from a Mendelian randomization study</a> Summary-based MR analysis. </li> Also, a bit weird conclusion: "This MR study suggests that genetically predicted alcohol consumption is a protective factor for EC, particularly for EEC, and this protective effect may be mediated through the reduction of HCG and IGF1." </li> What if people abstained from drinking because of health concerns, such as a family history of endometrial cancer? IV assumptions would be easily invalid. EC GWAS is also sampled from the same UK Biobank cohorts, which could put the study design in a precarious position. </li> </ul> </li> </ul> Okay, I may want to believe that alcohol is a cause of cancer. But I don't think we have a definitive answer to the underlying causal questions. Studies often disagree with each other. It's like opening up another decade of debate among researchers, policymakers, and the general public. However, I don't think we are missing a chance to potential protective effects of alcohol by not cutting down on drinking. A more constructive way forward will be looping in molecular mechanisms in causal inference exercises. We identified the mechanisms in laboratory settings. How does it work at a population level? Are there alcohol-interaction molecular QTL studies? Shouldn't we consider gene-gene interactions and compensatory mechanisms? Lots of "could've done" items. Large Language Model vs. Human Experts 2024-11-28T00:00:00+00:00 I've randomly stumbled upon this article, Large language models surpass human experts in predicting neuroscience results</a>. Should we worry about this? Tests were drawn from BrainBench, testing/checking neuroscience knowledge.</li> </ol> Co-authors (Supplementary Table 5) and GPT-4 (Azure OpenAI API; version 2023-05-15) created test cases that formed BrainBench. All test cases were sourced from Journal of Neuroscience abstracts published in 2023 under the Creative Commons Attribution 4.0 International License (CC-BY). The abstracts are organized into five sections, namely, behavioural/cognitive, systems/circuits, neurobiology of disease, development/plasticity/repair and cellular/molecular. </blockquote> LLMs outperformed human experts. LLM ($>$ 0.8 accuracy) vs. predoctoral student ($\approx$ 0.65), doctoral student ($\approx$ 0.6), postdoctoral researcher ($\approx$ 0.65), faculty ($\approx$ 0.65). Shouldn't we (humans) be more critical in paper reading? </li> Should we update neuroscience textbooks? </li> Not so significant. But what happened to the doctoral students? </li> </ul> </li> Human expertise not necessarily aligns with LLM expertise. </li> </ol> Perplexity measures how surprising a text passage is to an LLM. Using these measures (Supplementary Fig. 6), the mean Spearman correlation between an LLM and human experts was 0.15 ($\pm$0.03), whereas the mean Spearman correlation between LLMs was 0.75 ($\pm$0.08). </blockquote> No indication of brute-force memorization...</li> </ol> We found no indication that BrainBench was memorized by LLMs (Supplementary Fig. 7). ... As a final check (Methods and Supplementary Fig. 8), we confirmed that LLMs do not perform better on items published earlier in 2023 (for example, January 2023 versus October 2023), which addresses the concern that early items are more likely to have a preprint or other precursor appear in the training set that affected BrainBench performance. </blockquote> Github repo: https://huggingface.co/BrainGPT</code></li> </ol> Why we need to consider learning a new language 2024-11-25T00:00:00+00:00 TL;DR: Consider learning a new programming language. Some key concepts are only available in certain languages. When I was a graduate student, Python</code> was not as popular as they are now. Matlab was the ultimate sweet spot for much-needed speed, ease of programming, powerful visualization, and simulation. It was complete liberation. I no longer needed to deal with layers and layers of inheritance of Java classes. Even silly features like that I can type "why" on Matlab's terminal to get a nonsense answer were charm. However, as soon as I embarked on projects involving graph-based traversal and classical graph algorithms, Matlab was no longer an ideal tool for me. Everything underneath quick and easy API needs to be implemented with C++</code> codes within a notorious MEX</code>, trying hard to package with well-designed template libraries. I first turned to Python 2.7</code> and survived many years with it, but I was just not so happy with its indentation rule; often, my codes were barely running (at least, I felt like that). Later, I explored O'CaML</code> to bring back the ultimate joy of functional programming. Back then, I had many versions of the same MCMC algorithm written in multiple languages. Arguably, being bilingual or polylingual is beneficial in many ways. Also, arguably, learning a new programming language is like an abridged version of learning a new natural language. Learning a new language is perhaps far from mapping one to another language. If I had to think of everything in Korean and translate it into English sentence by sentence and word by word, well, then I would not have taken full advantage of every aspect of learning a new language. Likewise, if I had to code up a prototype in R</code> and convert it to C++</code>, then I should seriously consider picking up C++</code> instead of attempting quick fixes by asking any chatbots. Key concepts like memory management are simply hidden underneath the rug, so-called convenience. If any chatbots do a wonderful job, it would brutally deprive us of any possible true learning opportunities. There is no ultimate programming language, but many modern languages were designed to address diverse computing problems. For an applied scientist, having many things in their tool belts can make a big difference. It could shorten tasks from days to hours or even minutes. Learning awk</code> will be an eye-opening experience if one needs to deal with tens of thousands of lines of compressed text files. However, I would not want to use awk</code> (although I tried once) to transpose row-wise tabular data into column-wise ones. Although I heavily depend on well-designed and fast routines of data.table</code> in R</code>, I would never attempt to implement recursion in R</code> and often get uncomfortable with messy and unsafe parallelization. I think my brain has well past the peak age of agility and flexibility. It takes a substantially longer time for me to catch up with a new one. Nevertheless, I still feel some urgency. Language shapes the way I think. In many ways, a programming language is a language of science. It can teach us how to abstract away the underlying computing problems of real-world scientific problems. Some people might do computational thinking with pseudocode, but I do with some language I am currently into. Since different languages have been developed to serve different types of computing problems, it is definitely worth diversifying our thoughts. STAT 548 PhD Qualifying Papers 2024-09-03T00:00:00+00:00 Introduction</a></h2> I am interested in almost all problems in computational biology and genomics. I expect a student to propose novel statistical approaches that can address challenges in data analysis and modelling of high-dimensional, large-volume biological problems. Feel free to contact me (ypp@stat.ubc.ca</code>). Format</a></h2> You may organize your report including the following sections. Problem definition (1 page): Extract mathematical/statistical problems from the paper and organize them. What are the input data? What is the expected output? </li> Significance (1 paragraph): Why is this an interesting problem? What can be learned by studying this problem? Why is it exciting for you? Author contribution: How did the author(s) find the solution? What was a novel contribution beyond traditional approaches? </li> Limitations/challenges (1 paragraph): What are the assumptions? Are they realistic? What are the technical limitations that the authors acknowledge or not? </li> Novel idea/methods (1-2 pages): Propose your idea and statistical methods. You could interpret the underlying problem in a different formulation. What are related problems/frameworks, but not adopted by the authors? </li> Results (1-2 pages): Include one figure that sketches your approaches. Show tables and figures that clearly demonstrate your methods. </li> Discussion (1 page): Briefly discuss what you have learned and what you would achieve if you were to develop this to a full paper. How would you validate your findings in independent studies, including wet-lab experiments? </li> </ul> Available Papers</a></h2> Madrigal, A., Lu, T., Soto, L. M., & Najafabadi, H. S. (2024). A unified model for interpretable latent embedding of multi-sample, multi-condition single-cell data</a>. Nature Communications, 15(1), 6573. </li> Qiu, Y., Sun, J., & Zhou, X.-H. (2023). Unveiling the unobservable: Causal inference on multiple derived outcomes</a>. Journal of the American Statistical Association, 1–12. </li> </ol> Papers claimed by other students</a></h2> Bridgeford, E. W., Powell, M., Kiar, G., Noble, S., Chung, J., Panda, S., Lawrence, R., Xu, T., Milham, M., Caffo, B., & Vogelstein, J. T. (2024). When no answer is better than a wrong answer: a causal perspective on batch effects</a>. In bioRxiv (p. 2021.09.03.458920). </li> Demirel, I., Alaa, A., Philippakis, A., & Sontag, D. (2024). Prediction-powered Generalization of Causal Inferences</a>. International Conference on Machine Learning, 10385–10408. </li> </ol> Some configuration updates 2024-08-30T00:00:00+00:00 Emacs for Rust with rustic</code> and lsp</code></a></h3> Learning Rust is my personal goal for this year. It is very strict in type checking and memory ownership. So, I needed some help until I could fully grasp its key concepts and syntax. I use Emacs for almost everything, including coding. After an hour of searching, I think I found a pretty good setup (a part of .emacs.d/init.el</code>): ;; Ensure that rust-mode is installed (use-package rust-mode :ensure t :config ;; Set the key binding for rust-format-buffer (define-key rust-mode-map (kbd "C-c C-f") #'rust-format-buffer)) (use-package lsp-mode :init (add-hook 'rust-mode-hook #'lsp) :custom (lsp-rust-server 'rust-analyzer)) (use-package rustic :ensure :bind (:map rustic-mode-map ("M-j" . lsp-ui-imenu) ("M-?" . lsp-find-references) ("C-c C-c a" . lsp-execute-code-action) ("C-c C-c r" . lsp-rename) ("C-c C-c q" . lsp-workspace-restart) ("C-c C-c Q" . lsp-workspace-shutdown) ("C-c C-c s" . lsp-rust-analyzer-status)) :config (add-hook 'rustic-mode-hook (lambda () (flyspell-mode -1))) ) ;; Stop reformatting every time on save (rust-disable-format-on-save) </code></pre> We can peek into documentation and fn</code> and struct</code> definitions. (setq lsp-inlay-hint-enable t) (setq lsp-ui-peek-enable t) (setq lsp-format-on-save nil) (setq rustic-format-on-save nil) (setq rust-format-on-save nil) (setq special-display-buffer-names '("rustic-compilation")) (setq special-display-function (lambda (buffer &optional args) (split-window) (switch-to-buffer buffer) (get-buffer-window buffer 0))) (use-package company :ensure :custom (company-idle-delay 0.25) ;; how long to wait until popup ;; (company-begin-commands nil) ;; uncomment to disable popup :bind (:map company-active-map ("C-n". company-select-next) ("C-p". company-select-previous) ("M-<". company-select-first) ("M->". company-select-last))) (use-package yasnippet :ensure :config (yas-reload-all) (add-hook 'prog-mode-hook 'yas-minor-mode) (add-hook 'text-mode-hook 'yas-minor-mode)) </code></pre> When we deal with multiple projects, unused and leftover workspace caused performance drawback. So, the following wrapper is pretty handy. (defun cleanup-lsp () "Remove all the workspace folders from LSP" (interactive) (let ((folders (lsp-session-folders (lsp-session)))) (while folders (lsp-workspace-folders-remove (car folders)) (setq folders (cdr folders))))) </code></pre> Overall, coding Rust is pretty fun, and working with lsp</code>, company</code>, yasnippet</code> can certainly teach me my mistakes interactively. They are quite resource-heavy, but I will need "a tutor" for a while. Fix unstable Internet</a></h3> Working from home for a couple of weeks, I learned that my network is very unstable, especially with the ethernet connection! The problem was my internet service provider is not top-notch, so the download and upload speed shouldn't exceed 100mb/s, but my ethernet device kept on trying to find better speed but couldn't do it. I found a nice solution that the network speed can be simply capped by disabling autoneg</code>: sudo lshw -C network </code></pre> sudo ethtool --change ${eth} speed 100 autoneg off </code></pre> Fix display issues</a></h3> Amazon Prime Day got me a new monitor, which is spacious and bulky enough to force me to fold up my laptop screen on my small desk. But I found Ubuntu keeps on activating my laptop screen. I can exclusively switch on my HDMI screen with this command: xrandr --output eDP-1 --off --output HDMI-1 --primary --auto </code></pre> A new config for blogging 2024-08-01T00:00:00+00:00 I've kept on telling my people that researchers should have a way to constantly write out what they think and learn, even though things may turn out to be wrong at some point. I posted this</a> (from Nathan Pyle's comic, "Strange Planet") on my door to remind me of the importance of conveying ideas with evidence. I also believe that paper writing starts at the beginning of any project with a small step at a time. Idea conception, writing a paper, and doing research</a>. Although I keep my markdown notes and scanned copies of scribbles, I have been searching for a quick and easy way to keep my ideas. First of all, I was fed up enough with R</code>'s bookdown</code> and Hugo</code> setup. It was more complex than I would ever need. I have some experience with old web programming, but Hugo</code> was too hard for me. After searching and searching, I found Zola</code> was simple enough, command line friendly, and editor friendly. It also took me some time to figure out, though. First, install zola</a> locally. More details on installation can be found: https://www.getzola.org/documentation/getting-started/installation/</code> and some web search. Then, we can set up a new zola</code> directory: zola init ${BLOG} </code></pre> Here, ${BLOG}</code> can be any working directory. I was too lazy to tweak HTML pages and styles, so I decide to use Apollo theme</a>. git submodule add https://github.com/not-matthias/apollo themes/apollo </code></pre> or simply clone: git clone https://github.com/not-matthias/apollo themes/apollo </code></pre> Out of the box, things work fine, but it's worth customizing config.toml</code>: theme = "apollo" taxonomies = [{ name = "tags" }] # The URL the site will be built for base_url = "CHANGE ME" # Whether to automatically compile all Sass files in the sass directory compile_sass = true # Whether to build a search index to be used later on by a JavaScript library build_search_index = true output_dir = "CHANGE ME" generate_feeds = true [markdown] # Whether to do syntax highlighting # Theme can be customised by setting the `highlight_theme` variable to a theme supported by Zola highlight_code = true highlight_theme = "ayu-light" [extra] # See this for more options: https://github.com/not-matthias/apollo/blob/main/config.toml#L14 toc = true use_cdn = false theme = "toggle" socials = [ ] menu = [ ] mathjax = true mathjax_dollar_inline_enable = true fancy_code = true insert_anchor_links = "heading" favicon = "/icon/favicon.ico" </code></pre> Finally, we can create _index.md</code> as a front page: emacs content/_index.md </code></pre> _index.md</code> could be something like this: +++ title= "my blog" template = "homepage.html" +++ Welcome to my blog </code></pre> We can ask zola</code> to generate HTML and js</code> pages: zola build </code></pre> A new site will be located under where the output_dir</code> is pointing to. We can simply push the entire directory to GitHub</code> or a similar kind. Every time we need to update, we can create a new 'md' file under 'content' where different. To have a handle to a new set of posts, we need to add a new menu item in config.toml</code>: menu = [ { name = "/posts", url = "/posts", weight = 1 }, ] </code></pre> In content/posts/</code>, we can add _index.md</code> to help navigate posts: +++ paginate_by = 7 title = "Posts" sort_by = "date" insert_anchor_links = "heading" +++ </code></pre> Happy blogging! Pattie Ye 2024-05-31T00:00:00+00:00 I was a Master’s student in Bioinformatics starting in Fall 2021 under the co-supervision of Dr. Yongjin Park and Dr. Ramon Klein Geltink</a>. I completed my BSc from the University of British Columbia in 2021 with a combined major in Computer Science and Statistics. My main research interests pertain to studying the underlying biological mechanisms of autoimmune and inflammatory disease using large scale data analysis. Outside of academia, I am a member of the Varsity Men’s Rowing team at UBC as a coxswain. I served as rowing team captain for two years and have won two national championships. I also enjoy cooking, reading, and taking care of houseplants. Research interests: Computational Biology</li> Statistical Genomics</li> Complex Disease</li> Autoimmune and Inflammatory Disease</li> Immunometabolism</li> </ul> Yichen Zhang 2024-05-31T00:00:00+00:00 I was a Ph.D. student in the Department of Statistics. I employ techniques from deep learning and bayesian statistics to shed light on the mechanism of complex human diseases. I am generally interested in developing interpretable and scalable algorithms to uncover the underlying structure of biological data. Before that, I obtained my Master’s Degree in Biostatistics from the University of Washington, working on matrix completion and cell type deconvolution problems with Noah Simon. Besides research, I am a food enthusiast. I love cooking, exploring local restaurants, watching food & cooking documentaries and brewing coffee and tea. I do weight training from time to time so that I can enjoy more food. Research interests: Computational Biology</li> Bayesian Statistics</li> Deep Learning</li> Graphical Model</li> </ul> Yaoyue (Yolanda) Feng 2024-05-01T00:00:00+00:00 I hope to use the combination of bioinformatics and traditional biological experiments to explore the disease progression and drug resistance in Acute Myeloid Leukemia, in order to help more clinical patients. I obtained my Master of Research Degree in Cancer Biology from Imperial College London, working on Hepatocellular Carcinoma under hypoxia condition. And I studied cisplatin resistance in Ovarian Cancer at Sichuan University during my undergraduate stage. In my spare time, I prefer to stay at home watching movies or cooking instead of outdoor activities. I only go out with my friends when the weather is really nice. I also like basketball, but I watch it more often than I play it myself. Maxwell Douglas 2024-01-15T00:00:00+00:00 Hi! I'm a computational biologist passionate about the intersection of software, analytics, and biologic/medical data. I love talking shop about the how, when, and why of deriving insights from data. My research background is heavy in the various 'omics' technologies with a particular focus on genomics. interests: Computational Biology</li> Computational Oncology</li> Open-source software</li> Causal Inference</li> Data-Science</li> Translational Science</li> Public Health</li> </ul> Hobbies: Trail Running</li> Soccer</li> Skiing</li> Climbing</li> Backpacking/hiking</li> DND</li> Volunteering & community-building/outreach</li> Struggling to keep house-plants alive</li> </ul> Adjusting batch effect in pseudobulk scRNA 2023-06-19T00:00:00+00:00 Inspiration</a></h3> Our previous work demonstrates that fiding causal difference at a pseudobulk level is easier than dealing with noisy cell-level data: Counterfactual inference for single-cell gene expression analysis</a> </li> Interestingly, other people considered identifying batch effects should be treated as a causal effect estimation problem: Batch Effects are Causal Effects: Applications in Human Connectomics</a> </li> </ul> A generative scheme for a single-cell count matrix with multiplicative batch effects</a></h3> We encounter single-cell expression data consisting of multiple batches. One of the primary goals is to identify cell types (clusters/factors) and cell-type-specific gene expression patterns. However, distinguishing batch-specific and cell-type-specific genes only by a factorization method is challenging and often not identifiable from data alone. For each gene $g$ and cell $j$, the gene expression $Y_{gj}$ were sampled from Poisson distribution with the rate parameter: $$\lambda_{gj} = \lambda_{gj}^{\textsf{unbiased}} \times \prod_{k} \delta_{gk}^{X_{kj}},$$ affected by the batch effects $\delta_{gk}$. More formally, letting $X_{kj}$ be a batch membership matrix, assigning a cell $j$ to a batch $k$ if and only if $X_{kj}=1$, we assume the average gene expression rates are linearly affected by in the log-transformed space: $$\mathbb{E}\left[\ln Y_{gj}\right] = \ln \left( \sum_{t} \beta_{gt} \theta_{jt} \right) + \sum_{k} \ln\delta_{gk} X_{kj}.$$ set.seed(1331) m <- 500 # genes n <- 1000 # cells nb <- 2 # batches ## 1. batch membership X <- matrix(0, n, nb) batch <- sample(nb, n, replace = TRUE) for(b in 1:nb){ X[batch == b, b] <- 1 } ## 2. batch effects W.true <- matrix(rnorm(mnb), m, nb) ln.delta <- apply(W.true %% t(X), 2, scale) ## 3. true effects K <- 5 .beta <- matrix(rgamma(m K, 1), m, K) .theta <- matrix(rgamma(n K, 1), n, K) lambda.true <- .beta %% t(.theta) lambda <- lambda.true exp(ln.delta) yy <- apply(lambda, 2, function(l) sapply(l, rpois, n=1)) oo <- order(apply(t(.theta), 2, which.max)) </code></pre> If we can accurately estimate a true batch effect matrix, say $\delta_{gk}$, it is straightforward to adjust the difference between batches. How can we identify the true batch effect $\delta_{gk}$ for all the genes $g$ specifically expressed in the batch $k$? If we match cells $i$ and $j$ sampled from the batches $a$ and $b$, respectively, we expect the batch-specific difference $\delta_{ga} \neq \delta_{gk}$ will persist and even amplify, but the difference originated from cell types will vanish. This problem is equivalent to estimating the potential outcome of gene expressions in each batch $k$, $\mathbb{E}\left[Y_{gj}^{(k)}\right]$. A causal inference approach to identify batch effects</a></h3> To dissect batch-specific effect in a causal inference (potential outcome) framework, we assume our confounding variables $Q$ are well-distributed across different batches: Overlap: $0 < p(X_{kj}=1|Q) < 1$ for all $k$.</li> </ul> Moreover, we assume these covariates are sufficient enough to induce conditional dependence between potential (imputed) gene expression and batch assignment mechanisms: Strong ignorability: $(Y(k), Y(k')) \perp\perp X | Q$ for all $k,k'$ pairs.</li> </ul> Estimation of the batch effects by matching</a></h3> Suppose we can counterfactually estimate gene expressions of a certain cell $j$ if the cell was measured in different batches other than the observed batch $k$. $$Z_{gj} = \frac{ \sum_{i} (1 - X_{ik}) w_{ji} Y_{gi} }{ \sum_{i} (1 - X_{ik}) w_{ji} }$$ Like many other batch correction methods invented for single-cell RNA-seq analysis, we will assume $Z_{gj}$ reliably contain biologically-relevant cell state information while excluding the batch-specific effects to which the cell $j$ belong. Observed log-likelihood: $$\prod_{j} p(Y_{gj}|\mu_{gs},\delta_{gk},X_{jk}) =\prod_{j} \operatorname{Poisson}(Y_{gj}|\mu_{gs} \sum_{k} \delta_{gk} X_{jk})$$ Counterfactual log-likelihood: $$\prod_{j} p(Z_{gj}|\mu_{gs}, \gamma_{gs}) = \prod_{j} \operatorname{Poisson}(Z_{gj}|\mu_{gs} \gamma_{gs})$$ Local update: Maximize batch $s$-specific parameters</a></h4> Let's update $\mu_{gs}$ for a gene $g$ in a sample $s$: $$\mathbb{E}\left[\mu_{gs}\right] \approx \frac{ \sum_{j \in C_{s}} Y_{gj} + \sum_{j \in C_{s}} Z_{gj} }{\sum_{k} \delta_{gk} n_{sk} + n_{s} \gamma_{gs}}$$ Letting $p_{sk} = n_{sk} / n_{s}$, $$\mu_{gs} \gets \frac{ Y_{gs} + Z_{gs}}{\sum_{k} \delta_{gk} p_{sk} + \gamma_{gs}}$$ If $\delta_{gk} \to 0$ and $p_{sk}=1$, meaning that this sample $s$ is just sampled from the batch $k$ only, $\mu_{gs} \to Y_{gs} + Z_{gs}$ and $Y_{gs} \to Y_{gsk} = 0$. Therefore, $\mu_{gs} \to Z_{gs}$. Global update</a></h4> $$\mathbb{E}\left[\delta_{gk}\right] \approx \frac{\sum_{s} \sum_{j \in C_{s}} X_{kj} Y_{gj}}{\sum_{s} \mu_{gs} \sum_{j \in C_{s}} X_{kj}}$$ $$\delta_{gk} \gets \frac{\sum_{s} Y_{gsk} n_{sk}}{\sum_{s} \mu_{gs} n_{sk}}$$ If $Y_{gsk} \to \mu_{gs}$ for all $s$, $\delta_{gk} \to 1$. If $Y_{gsk} < \mu_{gs}$ in all $s$, $\delta_{gk} < 1$. If $Y_{gsk} \to 0$ for all $s$, $\delta_{gk} \to 0$. Algorithm</a></h4> Initialize batch effect $\delta_{gk} \gets 1$ for each gene $g$ and batch $k$ </li> Initialize $\gamma_{gs} \gets 1$ for each sample $s$ </li> Static global stat: $S_{gk} \gets 0$ </li> For each pseudo-bulk sample $s$ with cells $C_{s}$, $n_{sk} \gets \sum_{j \in C_{s}} X_{kj}$, $n_{s} \gets \sum_{k} n_{sk}$, $p_{sk} \gets n_{sk}/n_{s}$ </li> $Y_{gs} \gets \sum_{j \in C_{s}} Y_{gj} / n_{s}$ </li> $Y_{gsk} \gets \sum_{j \in C_{s}} Y_{gj} X_{kj} / n_{s}$ </li> $Z_{gs} \gets \sum_{j \in C_{s}} Z_{gj} / n_{s}$ after matching and imputation </li> $S_{gk} \gets S_{gk} + Y_{gsk} n_{sk}$ </li> </ul> </li> Iterative-updated global stat: $T_{gk} \gets 0$ </li> (Local step) For each PB sample $s$: $\delta_{gs} \gets \sum_{k} \delta_{gk} p_{sk}$ </li> $\mu_{gs} \gets (Y_{gs} + Z_{gs}) / (\gamma_{gs} + \delta_{gs})$ </li> $\gamma_{gs} \gets (Y_{gs})/(\mu_{gs})$ </li> For each $k$: $T_{gk} \gets T_{gk} + \mu_{gs} n_{sk}$ </li> </ul> </li> (Global step) For each batch $k$: $\delta_{gk} \gets S_{gk} / T_{gk}$</li> </ul> </li> Repeat the previous three steps (5-7) until convergence </li> </ol> A toy example</a></h3> ## 1. project K <- 5 R <- matrix(rnorm(m K), K, m) Q.raw <- R %% yy # K x n </code></pre> Before we adjust batch membership in the random projection matrix: cor(t(Q.raw), X) </code></pre> ## [,1] [,2] ## [1,] 0.7617260 -0.7617260 ## [2,] 0.8283630 -0.8283630 ## [3,] 0.8099248 -0.8099248 ## [4,] -0.7250199 0.7250199 ## [5,] 0.6651915 -0.6651915 </code></pre> ## 2. regress out ## ## X theta = X inv(X'X) X' Y ## = U D V' V inv(D^2) V' (U D V')' Y ## = U inv(D) V' V D U' Y ## = U U' Y x.svd <- svd(X) U <- x.svd$u U.t <- t(x.svd$u) Q.t <- t(Q.raw) Q.t <- Q.t - U %% U.t %% Q.t Q <- t(Q.t) </code></pre> After we adjust the batch effects: cor(Q.t, X) </code></pre> ## [,1] [,2] ## [1,] -3.700647e-16 3.700647e-16 ## [2,] -6.145885e-16 6.145885e-16 ## [3,] 1.390839e-16 -1.390839e-16 ## [4,] 3.551708e-16 -3.551708e-16 ## [5,] -7.981593e-17 7.981593e-17 </code></pre> q.svd <- svd(Q) ## 3. sorting B <- (sign(q.svd$v) + 1)/2 ss <- apply(sweep(B, 2, 2^(seq(0,K-1)), ``), 1, sum) + 1 </code></pre> feat.dn <- apply(Q, 2, function(x) x / sqrt(sum(x^2))) knn <- 3 d <- nrow(feat.dn) library(RcppAnnoy) ## a. construct dictionary for each batch dict.list <- lapply(sort(unique(batch)), function(b) { new(AnnoyAngular, d) }) for(j in 1:length(batch)){ b <- batch[j] dict.list[[b]]$addItem(j, feat.dn[,j]) } for(dd in dict.list){ dd$build(50) } ## b. a simplified routine to retrieve and estimate counterfactual y .counterfactual <- function(j){ v <- feat.dn[,j] nn <- c() dd <- c() for(k in 1:nb){ if(k == batch[j]) next .nn <- dict.list[[k]]$getNNsByVector(v, knn) .dd <- apply(feat.dn[, .nn], 2, function(u) sum((u - v)^2)) nn <- c(nn, .nn) dd <- c(dd, .dd) } w <- exp(-(dd - max(dd))) w <- w/sum(w) yy[, nn, drop = FALSE] %% matrix(w, ncol=1) } </code></pre> ngene <- nrow(yy) nbatch <- ncol(X) nsample <- max(ss) .delta.db <- matrix(1, ngene, nbatch) # gene x batch effects .delta.num.db <- matrix(0, ngene, nbatch) # gene x batch numerators .delta.denom.db <- matrix(0, ngene, nbatch) # gene x batch denominators .prob.bs <- matrix(0, nbatch, nsample) # batch x sample probabilities .size.bs <- matrix(0, nbatch, nsample) # batch x sample freq .ybar.ds <- matrix(0, ngene, nsample) # gene x sample observed average .zbar.ds <- matrix(0, ngene, nsample) # gene x sample imputed average .mu.ds <- matrix(1, ngene, nsample) # gene x sample adjusted average ## Precalculate some statistics for(s in 1:nsample){ if(sum(ss == s) < 1) next .yy <- yy[, ss == s, drop = FALSE] .zz <- do.call(cbind, lapply(which(ss == s), .counterfactual)) .ybar.ds[,s] <- apply(.yy, 1, mean) .zbar.ds[,s] <- apply(.zz, 1, mean) .prob.bs[,s] <- colMeans(X[ss == s, ]) .size.bs[,s] <- colSums(X[ss == s, ]) .y.dsb <- yy[, ss == s, drop = FALSE] %% X[ss == s, , drop = FALSE] .delta.num.db <- .delta.num.db + .y.dsb } .gamma.ds <- matrix(1, ngene, nsample) for(iter in 1:100){ .mu.ds <- (.ybar.ds + .zbar.ds) / (.delta.db %% .prob.bs + .gamma.ds + 1e-8) .gamma.ds <- .zbar.ds / (.mu.ds + 1e-8) .delta.db <- .delta.num.db / (.mu.ds %% t(.size.bs) + 1e-8) } </code></pre> Can we recover the original batch effects? par(mfrow=c(1,2)) plot(.delta.db[,1], W.true[,1], xlab="estimated delta (batch 1)", ylab="true delta effect") plot(.delta.db[,2], W.true[,2], xlab="estimated delta (batch 2)", ylab="true delta effect") </code></pre> Are they independent of the cell type effects? par(mfrow=c(1,2)) y.true <- sweep(lambda.true %% X, 2, colSums(X), `/`) plot(.delta.db[,1], y.true[,1], xlab="estimated delta 1", ylab="true y mean 1", main = "batch1") plot(.delta.db[,2], y.true[,2], xlab="estimated delta 2", ylab="true y mean 2", main = "batch2") </code></pre> While adjusting the estimated batch effects, can we recover the unbiased cell type effects? The following is before adjustment: ybar <- sweep(yy %% X, 2, colSums(X), `/`) par(mfrow=c(1,2)) plot(log1p(ybar[,1]), y.true[,1], xlab="log1p sample mean", ylab="true y mean #1", main = "batch1") plot(log1p(ybar[,2]), y.true[,2], xlab="log1p sample mean", ylab="true y mean #2", main = "batch2") </code></pre> Here, we adjusted the batch effects: ybar.adj <- sweep((yy / .delta.db[, batch]) %% X, 2, colSums(X), `/`) par(mfrow=c(1,2)) plot(log1p(ybar.adj[,1]), y.true[,1], xlab="adjusted sample mean", ylab="true y mean", main = "batch1") plot(log1p(ybar.adj[,2]), y.true[,2], xlab="adjusted sample mean", ylab="true y mean", main = "batch2") </code></pre> par(mfrow=c(1,2)) .tsne <- Rtsne::Rtsne(log(1 + t(yy)), num_threads=4)$Y plot(.tsne[,1], .tsne[,2], xlab = "tsne1", ylab = "tsne2", col = batch, main = "before") .tsne <- Rtsne::Rtsne(log(1 + t(yy/.delta.db[,batch])), num_threads=4)$Y plot(.tsne[,1], .tsne[,2], xlab = "tsne1 adj", ylab = "tsne2 adj", col = batch, main = "after") </code></pre> Ming Yuan 2022-09-01T00:00:00+00:00 I am a PhD student in Biostatistics at UBC. Before this ongoing MSc/PhD experience, I completed my BSc degree at Sun Yat-sen University, China, majoring in Statistics in 2020. My main research interests are causal inference and its application to genomic data. interests: Variable selection</li> Causal inference</li> </ul> Sam Khalilitousi 2022-01-31T00:00:00+00:00 I have always been intrigued by the ever-growing biological insights that are made possible with computers. This led me to pursue Biomedical Engineering (specializing in Cellular Bioengineering) at UBC, where I have been gaining experience in machine learning and bioinformatics analysis. With a passion centred in machine learning and computational biology, I want to use creativity to help make novel discoveries. Looking forward to connecting and collaborating with the team at the Park Lab! https://samktt.github.io/sam/</a> Kevin Lam 2022-01-15T00:00:00+00:00 I am a Statistics PhD student, supervised by Dr. Yongjin Park and Dr. Benjamin Bloem-Reddy</a>. I am mainly interested in applying and developing methods that help characterize underlying structures in biological data sets, for example using causal inference. Before the start of the PhD, I spent a year as a research assistant at the MRC Centre for Global Infectious Disease Analysis at Imperial College London where I collaborated with the Malawi Ministry of Health. Our team provided HIV modelling tools to support stakeholders with their policy decisions such as selecting new locations for health facilities for people living with HIV. During my MSc studies, I also worked on identifying the origins of drug resistant malaria in South East Asia using community detection in network analysis. Outside of my studies, I enjoy watching movies or discovering new restaurants. I also play the piano from time to time when the opportunity arises. Research interests: Causal Inference</li> Computational Biology</li> Machine Learning</li> Public Health</li> </ul> Yongjin P Park 2020-09-20T00:00:00+00:00 I have been an assistant professor in Pathology and Laboratory Medicine, as well as Statistics at UBC</a>. I hold Canada Research Chair Tier-2</a>. My research primarily focuses on developing scalable probabilistic inference methods to elucidate hidden causal mechanisms of human diseases, such as cancer and other common/complex disorders. Note: I started adding "P" as a middle initial since 2020 to avoid confusion with many other researchers with the same first and last name. "P" may be anything among "Probably," "Peter," "Projecting," etc. Besides research, I love cooking while trying out different types of cuisines and recipes. I also like to watch history documentaries and films and read books on history. When I was young, I used to play Starcraft and practice hip-hop dance, but I gave them up once I realized that the joy of coding and doing science exceeded and was more rewarding to me. Lightweight tools for integrative single-cell data analysis 2020-09-06T00:00:00+00:00 There are already many tools.</a></h3> Ever since single-cell technology was introduced to the world, every year, every month, we see more tools. See scRNA-tools</a>, for example. Single-cell data analysis often involves many steps of data processing, quantification, and statistical inference. The technology is imperfect, and each step deserves careful consideration of different types of approaches. We, bioinformaticians, test a different combination of small tools and tailor a pipeline that works best for each project. Perhaps, there is no master algorithm, and the competitions among computational biologists are indeed fierce. The history of genomics repeats itself. For a large part of current single-cell data analysis, we probably see what had happened in RNA-seq, even gene expression microarray analysis. Once we knew how to normalize scRNA-seq data, we would like to know how to test differentially-expressed genes. Then, we would want to build more audacious models, such as gene regulatory networks. Are there any points that we may contribute from a unique perspective? Can existing tools address all the challenges of single-cell data analysis, e.g., "eleven grand challenges in single-cell data science"</a>? I believe not. We know there has always been enough room for innovation and creativity. Okay, that sounds meek and quite passive. Why do we dare to jump into the battlefield? We are not a software-engineering lab, but we can become impatient enough to feel the need for a better tool. Moreover, we want to understand why and how a certain method works, and when that method doesn't work as expected. The underlying data matrix is sparse, really sparse, but most methods make it dense.</a></h3> Yes, the data matrix is sparse, but it will easily blow up memory if you fill in all the zeros. Our first impression was that memory overhead was never an issue for most people. Lots of methods focus on replacing the zeros with some other fractional values. Not undermining the validity of imputation methods, we worried more about the scalability of the strategy. Can such a strategy work for millions of cells with billions of non-zero elements? Most methods are tested on a rather small data set, including thousands of cells. If a data matrix is sparse, why not we build a method, taking advantage of the sparsity? What if most zero values did not occur by missingness? What if droplet scRNA-seq is not zero-inflated</a>? We develop a scalable approach to index millions of cells for fast random access. We like to work on a sparse data matrix formatted as a list of triplets</a>, and compressed by gzip</code> or bgzip</code>. The rows correspond to features/genes; the columns correspond to UMI barcodes/cells. We assume the triplets (row, column, value)</code> are sorted by the columns (cells) in ascending order. To mark the line where the columns change, we modify the indexing scheme of TABIX</a> and create a separate .index</code> file to store all the random access points. We implement commonly-used matrix operations to manipulate data more swiftly. In some cases, it is worthwhile to re-implement low-level functions, since we may not need all the steps implemented in a general-purpose library. We were just too impatient to wait a long I/O time until we finish reading a full 10x matrix (.mtx.gz</code>) using either scipy.io.mmread</code> in Python, or Matrix::readMM</code> in R. As the size of the 10x matrix increase more and more, we also discovered there is a numerical limit for R's Matrix::readMM</code> (less than two billions non-zero elements) and had to face a similar issue in a low-memory machine, such as a modest 8Gb laptop. So, we did the dirty job of C++</code> implementation. Combined with the indexing scheme, now we can efficiently take a subset of a large matrix stored in a file, not worrying about memory capacity, minimizing the overhead of character parsing and I/O. We want to focus on integrative data analysis across multiple data types and modalities</a></h3> With these tools, including the ones under development, we want to do science. We want others to do science using our tools. We will continuously update a bit more complicated routines for convenience, including cell type annotation, data aggregation, cell type deconvolution, batch-balancing k-nearest neighborhood, etc. Causal Inference in Genomics 2020-09-04T00:00:00+00:00 TL;DR. We reinvent traditional bioinformatics methods to ascertain the causality of discoveries made in high-dimensional omics data. Massive data generation in genomics has transformed the methodology and practice of researches in medicine and human biology. A common practice of research is genomic and functional genomic profiling of hundreds of individuals. If desired data is not available, large consortium-level projects and collaborations can often meet the needs. In terms of both volume and dimensionality, the sizes of genomics data have exploded in the past decade, and it is not hard to predict that the speed of data accumulation will be further accelerated. Genomics has made a breakthrough in medicine and science, armed with two powerful scientific methods--hypothesis testing and machine learning (ML). Hypothesis testing with well-designed null distributions proves meaningful scientific discoveries are made on the observed data. Likewise, a traditional ML method is designed to efficiently search the most plausible models from a particular class of models, mainly focusing on minimizing generalization errors on the observed data. However, biology's ultimate question is perhaps more about unobserved principles, asking counterfactual "what if" questions, rather than simply describing phenomena. We desire to uncover causal mechanisms and laws underlying observations and ultimately provide testable cause-effect relationships. Causal inference based on a causal structure model</a></h3> Most causal inference concepts based on a structural equation model are already adopted in a carefully designed biological experiment. We can trace back to Sewall Wright's path diagram model. For instance, an experimental proof demonstrating the effect of a gene "X" to some trait "Y" would require two types of interventions that disruption of the gene X makes changes to the trait Y, and recovery of the gene can attenuate the disrupted functions. However, a common practice of high-throughput data analysis often lack such intervention steps, but hurriedly arrive at some conclusion based on observed correlation patterns. Many of these correlations may turn out non-causal but somewhat confounded by unobserved variables. Identification of causal effects by an instrumental variable (mediation)</a></h4> The key concept is an intervention. However, introducing an intervention is often impossible for technical and ethical reasons, and there are so many things to perturb to test the effects. Moreover, the effects of causal mechanisms might be too small to quantify within a reasonable statistical error reliably. How do we ascertain causal effects by intervention? Stepping back, let us think about two principal axes of the data generation process: the breadth of population-level measurements across tens of thousands of samples and the depth of regulatory genomics across multiple layers of causal mechanisms, seeking a precise answer up to a single-cell resolution. {{< figure library="true" src="project/cancer_genomics_research.png" title="Integrative cancer research builds on two principal axes of genomics data. (A) Population-level studies investigate variation across large cohorts. (B) Regulatory genomics researches delve into multiple layers of molecular networks." >}} We can recognize a causal direction from genetic information to phenotypic variation if we consider that genotypes are shaped by nature’s randomized controlled trial (RCT), or natural interventions. Given that, the goal is to identify causal mediators located in the middle of the above causal diagram. In the context of genome-wide association studies (GWAS), the mediator variables include relevant cell types and target genes derived from tissue-level or cell-type-level eQTL data. Removing unwanted variability by "control" data</a></h3> Many scientific discoveries are best represented as a set of contrastive statements, such as "A rather than B" (Peter Lipton</a>) because such a contrastive explanation clarifies our claim's scope and basis. The validity of causal argument can be empirically justified by triangulation of contrastive arguments (Lipton 1991). Therefore, an essential aspect of our research should focus on finding a suitable set of null hypotheses (non-causal facts) to be contrasted against our theory of interest to uphold. What are non-causal effects that may obfuscate the activity of scientific research? Are they generated uniformly at random? In practice, it is hard to characterize non-causal effects without considering confounding variables, and confounders create unwanted non-zero correlations between cause and effect variables. Confounder correction is a causal inference problem. So, we need to ask the question of identifiability, such as, "How do we know this variable is a confounder?" or "What is a legitimate criterion that distinguishes between causal and non-causal effects?". We could start our journey by asking: Are there any control data on which we can safely claim that there is no causal effect whatsoever? </li> Are they already available/observed? Or, should we construct/estimate them? </li> </ul> Like many research problems, control data (or variables) are hidden, yet to be estimated. In the mediation problem, for instance, we need to think about how to characterize confounding variables that relate two mediators (genes, M1 and M0) (A) and the other type of confounders between different phenotypes (B). Unless we adjust these confounding effects, mediation effects are not identifiable by statistical inference. Fortunately, we can construct control data that only include non-genetic correlation structures by exploiting human genetics's nature. Since genetic correlations are confined within each linkage disequilibrium (LD) block, we can defuse putative causal effects by projecting summary statistics onto a different, independent LD block (C). We actively construct control data to identify confounding variables. {{< figure library="true" src="project/mediation_confounder.png" title="Confounder correction in action. The figure describes essential steps of confounding factor correction in GWAS mediation analysis. (A) Type I confounding variable (S1) between mediators (M's) (B) Type II confounding variable (S2) between phenotypes (Y's). (C) Joint projection of mediator and phenotype variables makes interventions to remove mediating effect but leave confounding effects.">}} On this matter, we are broadly interested in general genomics problems, not just statistical genetics ones. Although the definition of a confounding variable is problem-specific, we seek to design probabilistic models and ML algorithms that semi-automate the overall process. The methods will indicate a set of putative confounding variables with a level of uncertainty by taking input data and prior knowledge (e.g., known causal and non-causal variables). A potential outcome framework for causally-differential expression analysis</a></h3> Recently, our interest in Donald Rubin's potential outcome framework</a> grew, and started thinking about a type of causal problems that can be better approached by Rudin's causal model (RCM). "Bayesian" causal inference for single-cell data analysis</a></h4> As a first step, we revisit commonly-practised bioinformatics problems and reformulate them in RCMs. The ultimate goal of RCM is to ask and estimate counterfactual questions. Causality is tied to an action (doing), rather than an observation (seeing), and counterfactual causal inference seeks to take an impossible action and estimate its effects. Suppose we want to test the causal statement: A gene expression of Y was changed because of a disease X, i.e., $X \to Y$ with a directed edge (going beyond $Y \sim X$). We define $X$ to take 1 or 0 for "yes" or "no" and the resulting gene expression $Y$ to take a real number. In observational studies, including most high-througput assays, we observe a pair of $X_{i}$ and $Y_{i}$ for an individual $i$; there, the observed $Y_{i}$ value stems from the observed $X_{i}$, either 1 or 0, not both. E.g., $i$</th> X</th> Y(0)</th> Y(1)</th></tr></thead> 1</td> 0</td> $Y_1$</td> ?</td></tr> 2</td> 1</td> ?</td> $Y_2$</td></tr> 3</td> 1</td> ?</td> $Y_3$</td></tr> 4</td> 0</td> $Y_4$</td> ?</td></tr> </td> </td> $\vdots$</td> </td></tr> </tbody></table> In counterfactual inference, we want to infer the other unobserved and not observable cases. For $X_{i}=1$, we want to know the potential outcome of $Y_{i}$ if we had $X_{i} = 0$, and vice versa. RCMs generally translate the underlying problem as a Bayesian imputation problem, or a matrix completion problem, which contains at least 50% of missing values. Provided an accurate Bayesian imputation method, we can quantify individual-level causal effects, $Y_{i}(1) - Y_{i}(0)$, as well as average causal effects, $\mathbb{E}\left[Y(1) - Y(0)\right]$. This is truly substantial advancement because we no longer worry about identifiability but focus on a statistical inference problem. Well, there is a catch. In traditional omics data analysis, such an imputation problem is pretty challenging unless we have a good understanding of the data-generating process, which is often a black box to us. For one individual, we only have one data point $Y_{i}$ unless we borrow information from related genes. We might have to rely on faithful pretreatment variables $W$ as in a causal path, $W \to X \to Y$, and construct propensity scores to reweight observed samples. However, single-cell data provide a unique opportunity for modelling high-dimensional omics profiles. A clear advantage over bulk data is simply that we have more data points per sample (individual). Assuming that single-cell profiling's sensitivity will become better each year, it is not too far until we can achieve reliable and accurate Bayesian imputation methods--not only the observed values but also the counterfactual values. Bayesian modelling of summary statistics data 2020-09-03T00:00:00+00:00 TL;DR. We develop Bayesian machine learning methods that excavate a latent structure or underlying joint models of summary statistics data. What are summary statistics data?</a></h3> There is no clear boundary to define that data are made available at the summary or observation level. Suppose the size of a data matrix/tensor $Z$ is smaller than the underlying sample size (number of observations). Each element in the data matrix $Z$ is generated as a result of some transformation of raw data. In that case, we often consider that the data $Z$ is observed at the summary level. Summary statistics data could be first and second-order moments (mean and variance), p-values, z-scores, effect sizes, and standard errors, etc. As the most well-known summary statistics data, we could think of the results of genome-wide association studies (GWAS). Why do we bother summary statistics?</a></h3> The foremost reason is that people feel pretty comfortable with sharing summary statistics. Unlike we have to deal with strict regulations to access individual-level data, a user (researcher) may find easily-downloadable GWAS summary statistics vectors from a public database. For the other party, those who provide data, there is almost no need to anonymize it since each data point does not correspond to an individual. For most clinical and public health studies, research groups share their final results as a giant summary statistics table. Other researchers gain essential knowledge by using them in subsequent research. Great, but it sounds like everything is done and remains nothing to be done. I just said that summary statistics data are often final to some studies. Yes, to some people, but not to those who care about a generative model involving lots of parameters. Most summary data only investigate relationships between one variable with the other variable, but not many to many. Only a set of marginal probabilities, but not the joint probability of all. Can we recover a joint, multivariate model from a set of univariate models? Taking it further, can we estimate a joint probability of multiple studies conducted independently by different groups by combining summary data collected from each of them? Summary-based regression modelling</a></h3> Zhu and Stephens (2017) PMID: 29399241</a> proposed an interesting idea, termed regression with summary statistics (RSS) model. RSS recovers the high-dimensional parameters of a multivariate (linear) regression model, highlighting a correspondence between the multivariate and univariate association statistics. We also found the RSS type of modelling can be effective enough to capture causal SNPs in high-dimensional variable selection problems. Here, we show simulation results on the reference panel genotype matrix (the 1000 genomes project</a>; on chromosome 1, European samples, n=503). We implemented a suite of related summary-based methods zQTL</a>. {{< figure library="true" src="project/summary_regression.png" title="A summary-based learning method can estimate multivariate parameters as accurate as a traditional Bayesian learning method for a sparse regression problem. AUPRC</code>: Area Under Precision Recall Curve. Column panels: simulations varying the number of causal SNPs. h2</code>: the proportion of phenotypic variance explained by causal SNPs. BSLMM</code>: Bayesian sparse linear mixed effect model; yqtl</code>: a sparse regression model with spike-slab prior, trained by stochastic variational inference; zqtl</code>: a summary-based multivariate regression modelling.">}} Summary statistics-based machine learning</a></h3> Our next challenge is to widen the scope of computational biology, genomics applications that can be benefited from summary-based learning methods. About a decade ago, the ML community also realized the value of exploiting the spectral structure of the summary data matrix/tensor. Hsu et al. (2012) Hsu+Kakade+Zhang</a> beautifully demonstrates an optimal state of a hidden Markov model can be easily recovered by the spectral decomposition of marginal probability matrices.

American Society of Human Genetics meeting 2025

After the ASHG meeting</a></h2>

Sishir Subedi

STAT 548 PhD Qualifying Papers

Excerpts from Ideaflow on creativity - part 3

Excerpts from Ideaflow on creativity - part 2

Excerpts from Ideaflow on creativity - part 1

Embedding, attention, and transformers

Transparent parameterization for a partially identified model

A stochastic block model and decades later

ML/AI cheat sheet: Encoder

Single-cell (pseudotime) trajectory analysis

Brain Age Gap GWAS

Notes on 'How altered cancer metabolism influences cancer progression'

`.R/Makevars` on silicon Mac

Alcohol can cause cancers?

Large Language Model vs. Human Experts

Why we need to consider learning a new language

STAT 548 PhD Qualifying Papers

Some configuration updates

A new config for blogging

Pattie Ye

Yichen Zhang

Yaoyue (Yolanda) Feng

Maxwell Douglas

Adjusting batch effect in pseudobulk scRNA

Ming Yuan

Sam Khalilitousi

Kevin Lam

Yongjin P Park

Lightweight tools for integrative single-cell data analysis

Causal Inference in Genomics

Bayesian modelling of summary statistics data

Wanxin Li

Yi-Ting(Tina) Hsu

Study note on differential privacy

Wanxin Li

Yi-Ting(Tina) Hsu

Study note on differential privacy

American Society of Human Genetics meeting 2025

Oct 15</a></h2> Non-linear GWAS methods</a></h3> Refine additive association models with dominant and recessive models</li> Concerns over Type I error calibration</li>

Sishir Subedi

STAT 548 PhD Qualifying Papers

Excerpts from Ideaflow on creativity - part 3

Excerpts from Ideaflow on creativity - part 2

Excerpts from Ideaflow on creativity - part 1

Embedding, attention, and transformers

Transparent parameterization for a partially identified model

A stochastic block model and decades later

ML/AI cheat sheet: Encoder

Single-cell (pseudotime) trajectory analysis

Brain Age Gap GWAS

Notes on 'How altered cancer metabolism influences cancer progression'

`.R/Makevars` on silicon Mac

Alcohol can cause cancers?

Large Language Model vs. Human Experts

Why we need to consider learning a new language

STAT 548 PhD Qualifying Papers

Some configuration updates

A new config for blogging

Pattie Ye

Yichen Zhang

Yaoyue (Yolanda) Feng

Maxwell Douglas

Adjusting batch effect in pseudobulk scRNA

Ming Yuan

Sam Khalilitousi

Kevin Lam

Yongjin P Park

Lightweight tools for integrative single-cell data analysis

Causal Inference in Genomics

Bayesian modelling of summary statistics data