Day 2

Dimensionality Reduction�Coevolution and Contact Prediction

What is your favorite movie or tv show?

Conservation Review

Conservation can be displayed as a “weblogo”

MSA

Weblogo

What are we missing in our conservation/weblogo analysis?

1. Coevolution
2. Phylogeny

​

Why might “Phylogeny” be important?

1. Different groups of sequences might have different evolutionary constraints

Putting it all together

Subgroups of sequences might have different conservation patterns!

Phylogeny: Distance method Review

Distance matrix

Data matrix

Phylogenetics

Features

Samples

UPGMA/NJ

Alpha

Beta

Gamma

Delta

Epsilon

Why might we want to do this?

https://pubmed.ncbi.nlm.nih.gov/10381871/

Species tree

Gene tree

image adopted from evolution-textbook.org

gene family α

gene family β

Distance matrix

Alpha

Beta

Gamma

Delta

Epsilon

Alpha

Beta

Gamma

Delta

Epsilon

Distance matrix

Dimensionality reduction

6 Features

Samples

2 Features

X

Y

Plot

Multidimensional scaling

6 Features

2 Features

X

Y

Plot

Samples

MDS (Multidimensional scaling)�Find a low-dimensional embedding that minimizes the difference in distances.

A

dm(A)

B

dm(B)

Minimize the difference

(dm(A)-dm(B))²

LET’S TRY!

Dimensionality reduction

n (samples)

d (features)

k

n (samples)

A

B

Linear transformation

A

B

W

Linear transformation

A

B

A’

W

W^T

W

W^T

What do those “lines” mean?

W

W^T

What does this “line” mean?

A * M = B

Simple math!

2 * 3 = 6

What should be?

2

8

1

?

3

(2 * 3) + (1 * ?) = 8

?

2

8

1

2

3

(2 * 3) + (1 * 2) = 8

multiple lines and nodes

2

8

(2 * 3) + (1 * 2) = 8

(2 * 4) + (1 * 1) = 9

1

3

4

9

2

1

W

W^T

Minimize difference

(find W, so that A ≈ A’)

PCA (Principal component analysis)�W = top-k eigenvectors of the covariance matrix of A�

Transform (rotate)

A

B

W

Transform (rotate) and plot the data!

A

B

W

Y

X

MDS - multidimensional scaling

A

dm(A)

B

dm(B)

A

B

A

PCA - principal component analysis

PC2

PC1

MDS2

MDS1

PCA “rotates” the data to find axes that maximize the variance.

MDS finds a lower-dimensional manifold that preserves the distances.

LET'S TRY

https://www.pnas.org/doi/10.1073/pnas.1711913115

https://onlinelibrary.wiley.com/doi/10.1002/pro.4422

https://github.com/gelnesr/SVD-of-MSAs

Part 2

What can we learn from MSA(s)?

Review - Conservation

A

B

C

D

E

0

1

2

3

4

5

A

Sequence

Multiple Sequence Alignment

Conservation

BLAST

Review - Phylogeny

A

B

C

D

E

0

1

2

3

4

5

Phylogeny

0

4

3

2

2

4

0

3

6

2

3

3

0

3

5

2

6

3

0

4

2

2

5

4

0

Distance matrix

NJ/UPGMA

A

Sequence

Multiple Sequence Alignment

BLAST

Review - Dimensionality reduction

A

B

C

D

E

0

1

2

3

4

5

Phylogeny

0

4

3

2

2

4

0

3

6

2

3

3

0

3

5

2

6

3

0

4

2

2

5

4

0

Distance matrix

Dimensionality Reduction

NJ/UPGMA

MDS

PCA

A

Sequence

Multiple Sequence Alignment

BLAST

Today - Coevolution!

A

B

C

D

E

0

1

2

3

4

5

A

Sequence

Multiple Sequence Alignment

BLAST

Which positions appear to covary?

A

B

C

D

E

0

1

2

3

4

5

Which positions appear to covary?

A

B

C

D

E

0

2

1

3

4

5

Which positions appear to covary?

A

B

C

D

E

0

2

3

4

1

5

48

P(■■) = ?

0

1

49

Position-specific scoring matrix (PSSM)

■�■

0

1

■�■

Position-specific scoring matrix (PSSM)

0

1

Markov Random Fields

for Sequences

Alan S. Lapedes

Position-specific scoring matrix (PSSM)

■�■

P(■■)=exp(V₁[■]+V₂[■]+W₁₂[■][■])/Z = 0.36

​

Z = exp(V₁[■]+V₂[■]+W₁₂[■][■])

+exp(V₁[■]+V₂[■]+W₁₂[■][■])� +exp(V₁[■]+V₂[■]+W₁₂[■][■])

+exp(V₁[■]+V₂[■]+W₁₂[■][■])

Markov Random Fields

Cons:

• Likelihood calculations VERY slow
• 20^Length combinations for proteins
• For 100 aa proteins thats: 1.3x10¹³⁰
• For scale:�There are 1.0x10⁸² Atoms in the UNIVERSE!

+ reg(w,v)

pll - Besag et al 1975, 1977

GREMLIN - Balakrishnan et al. 2009�GREMLIN - Balakrishnan et al. 2011�GREMLIN_v2 - Kamisetty et al. 2013�plmDCA - Ekeberg et al. 2013

GREMLIN (Generative REgularized ModeLs of proteINs)

Let’s pretend there are:

• 3 positions
• 3 characters:

0

1

2

conservation

coevolution

0

1

2

x = one_hot(MSA)

x

0

1

2

x' = x@w

w

x

0

1

2

[w]eights = coevolution

x'

x' = x@w + b

b

x

0

1

2

[w]eights = coevolution

[b]ias = conservation

w

x' = softmax(x@w + b)

softmax(x) = exp(x)/sum(exp(x))

�To make sure probabilities at each position sum to 1.0

b

x

x'

0

1

2

w

Fully connected layer (Dense layer)�P(x₁,x₂,x₃) ≈ P(x₁|x₁,x₂,x₃) * P(x₂|x₁,x₂,x₃) * P(x₃|x₁,x₂,x₃)

x' = softmax(x@w + b)

�loss = -x*log(x')

loss

b

x

x'

0

1

2

w

Minimize difference using categorical-crossentropy��Can learn an identity mapping (self-connections)

Remove self-connections (Pseudo-likelihood)�P(x₁,x₂,x₃) ≈ P(x₁|x₂,x₃) * P(x₂|x₁,x₃) * P(x₃|x₁,x₂)

x' = softmax(x@w + b)

�loss = -x*log(x')

loss

b

x

x'

0

1

2

w

Minimize difference using categorical-crossentropy

Balakrishnan, S., Kamisetty, H., Carbonell, J.G. and Langmead, C.J., 2009. �Structure Learning for Generative Models of Protein Fold Families.

maximize PLL

+ reg(w,v)

pll - Besag et al 1975, 1977

GREMLIN - Balakrishnan et al. 2009�GREMLIN - Balakrishnan et al. 2011�GREMLIN_v2 - Kamisetty et al. 2013�plmDCA - Ekeberg et al. 2013

0

1

2

GREMLIN (Generative REgularized ModeLs of proteINs)

x' = softmax(x@w + b)

b

x

x'

0

1

2

w

x' = x@w

w

x

0

1

2

[w]eights = coevolution

x'

w

x

0

1

2

[w]eights = coevolution

loss

x'

x' = x@w

�loss = (x' - x)^2�w = I

analytical solution: Identity matrix

w

x

0

1

2

[w]eights = coevolution

loss

x'

x' = x@w

�loss = (x' - x)² /N - 2Tr(w)^�w = cov(x)^-1 + I

Direct Coupling Analysis

(analytical solution: inverse covariance)

Dauparas, J., Wang, H., Swartz, A., Koo, P., Nitzan, M. and Ovchinnikov, S., 2019. Unified framework for modeling multivariate distributions in biological sequences. arXiv��Morcos, F.,… Weigt, M., 2011. Direct-coupling analysis of residue coevolution captures native contacts across many protein families. PNAS

w

x

0

1

2

x'

x' = x@w@w_T

�loss = (x' - x)²�w = PCA(x).components_

Autoencoder�(analytical solution: principle component analysis)

loss

w^T

w

x

0

1

2

x'

x' = x@w@w_T

�loss = (x' - x)²�w = PCA(x).components_

Autoencoder�(analytical solution: principle component analysis)

loss

w^T

PC₁

PC₂

0

1

2

No connections (Conservation)

P(x₁,x₂,x₃) ≈ P(x₁) * P(x₂) * P(x₃)

loss

b

x

x'

0

1

2

x' = softmax(b)

�loss = -x*log(x')

x' = softmax(x@w + b)�loss = -x*log(x') + λb² + λw²

loss

GREMLIN (Generative REgularized ModeLs of proteINs)

b

x

x'

0

1

2

w

(W)eights of MRF

0

1

2

0

1

2

0

1

2

w

(W)eights of MRF can be “reduced” to a Contact map

0

1

2

0

1

2

0

1

2

0

1

2

0

1

2

norm

w

0

1

2

0

1

2

0

1

2

0

1

2

0

1

2

w

0

1

2

0

1

2

0

1

2

0

1

2

0

1

2

Low rank correction often required

Sparse matrix

low-rank matrix

Decomposition

APC = raw - AP

Dunn et al. 2008

Mutual information without the influence of phylogeny or entropy dramatically improves residue contact prediction.

Contact map

50S ribosomal protein L6

N

C

C

XRAY Contacts

How to read a contact map

50S ribosomal protein L6

N

C

C

XRAY Contacts

Overlay of predicted contacts on real contacts

XRAY Contacts

​

Predicted Contacts

LETS TRY

A

B

C

D

E

0

1

2

3

4

5

Sequences

Positions

Co-evolution

0

1

2

3

4

5

0

1

2

3

4

5

0

1

2

3

4

5

Tomorrow

MRF/Potts

(markov random fields)

MG

(multivariate gaussian)

AE

(autoencoders)

PSSM

(position-specific scoring matrix)

BERT/�Transformer

AlphaFold/�RoseTTAFold

All these models can be expressed as “Autoencoders”

PSSM

MRF

MG

AE

BERT

Bonus

Find reduced dimensionality that only preserves local distances (or neighbors)?

Why might we want to do this?

SM = Gaussian(DM)

Similarity Matrix

Distance Matrix

Similarity Matrix

A

dm

1

2

3

4

5

6

7

8

x

y

z

Samples

1

2

3

4

5

6

7

8

Samples

1

2

3

4

5

6

7

8

Samples

Features

sm

dm_ij = exp(-dm_ij² / 2σ²)

1

2

3

4

5

6

7

8

Samples

1

2

3

4

5

6

7

8

Samples

D_ij = (x_i-x_j)² + (y_i-y_j)² + (z_i-z_j)²

A

dm

B

dm

A

dm

B

dm

sm

sm

MDS

UMAP/TSNE

A

B

A

PCA

TSNE (t-distributed Stochastic Neighbor Embedding)

UMAP (Uniform Manifold Approximation and Projection)

MDS

UMAP/TSNE

Phylogeny

PCA

Encode integers by which prime-factors they share

8 million integers

78628D to 2D

https://johnhw.github.io/umap_primes/index.md.html

CHECK IT OUT

https://johnhw.github.io/umap_primes/three_umap.html

Cool interactive resources

1. PCA �https://setosa.io/ev/principal-component-analysis/
2. TSNE �https://distill.pub/2016/misread-tsne/
3. UMAP�https://pair-code.github.io/understanding-umap/�https://johnhw.github.io/umap_primes/index.md.html

​

w

x

0

1

2

x'

x' = softmax(x@w)@w_T

�loss = (x' - x)²�w = KMeans(x).means_

Categorical Autoencoder�(solution: kmeans)

loss

w^T