Scripts for general data analyses (Laboratory purposes) and some algorithms for learning purposes

(Principal Component Analysis python scripts)

Assume the presence of an arbitrary 2D Gaussian distribution within a two-dimensional space.

Let us sample 1,000 random points from this 2D Gaussian distribution.

Consider each point as a set of vectors in a two-dimensional space:

$$\mathbf{X} = \left(\begin{matrix} x_1&\cdots&x_n \\ y_1&\cdots&y_n \end{matrix}\right)^\mathrm{T}\in \mathbb{R}^{n\times 2}$$

The covariance matrix for $\mathbf{X}$ is:

$$\Sigma = \begin{pmatrix} V[\mathbf{X_1}]&Cov[\mathbf{X_1},\mathbf{X_2}] \\ Cov[\mathbf{X_1},\mathbf{X_2}]& V[\mathbf{X_2}] \end{pmatrix}\in \mathbb{R}^{2\times 2}$$

, where $\mathbf{X_1} = (x_1:\cdots x_n)$, $\mathbf{X_2} = (y_1:\cdots y_n)$.

Define a projection matrix for the linear transformation $\mathbb{R}^2 \to \mathbb{R}$ as:

$$\mathbf{w} = \begin{pmatrix}w_1&w_2\end{pmatrix}^\mathrm{T}$$

The variance of the projected points onto $\mathbb{R}$ is given by: $$s^2 = \mathbf{w}^\mathrm{T}\Sigma \mathbf{w}$$

We aim to maximize the variance such that the projected points correspond to the principal component.

To prevent the variance from diverging, we fix the norm of the projection matrix at 1. Therefore, we solve the following constrained maximization problem to find the projection axis:

$$arg \max (\mathbf{w}^\mathrm{T}\Sigma \mathbf{w}), |\mathbf{w}| = 1$$

To solve this maximization problem under the given constraint, we employ the method of Lagrange multipliers. This approach introduces an auxiliary function, known as the Lagrange function, to find the extremum of a function subject to constraints. The formulation of the Lagrange multipliers method applied to this problem is:

$$\cal{L}(\mathbf{w},\lambda) = \mathbf{w}^\mathrm{T}\Sigma \mathbf{w} - \lambda(\mathbf{w}^\mathrm{T}\mathbf{w}-1)$$

At the point of maximum variance: $$\frac{\partial\cal{L}}{\partial{\mathbf{w}}} = 2\Sigma\mathbf{w}-2\lambda\mathbf{w} = 0$$

Hence, we find that:

$$ \Sigma\mathbf{w}=\lambda\mathbf{w} $$

Select the eigenvector corresponding to the largest eigenvalue (where $\lambda_1 > \lambda_2$) as the direction of the principal component.

A singular value decompositon of $\mathbf{A}$ is written as:

$$\mathbf{A}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^\mathrm{T}$$

, where

$$\mathbf{A}\in\mathbb{R}^{m\times n},\mathbf{U}\in\mathbb{R}^{m\times m},\mathbf{\Sigma}\in\mathbb{R}^{m\times n},\mathbf{V}\in\mathbb{R}^{n\times n}$$

Right Singular Vectors:

$$\mathbf{A}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^\mathrm{T}$$

$$\mathbf{A}^\mathrm{T}\mathbf{A}=(\mathbf{U}\mathbf{\Sigma}\mathbf{V}^\mathrm{T})^\mathrm{T}\mathbf{U}\mathbf{\Sigma}\mathbf{V}^\mathrm{T}$$

$$\mathbf{A}^\mathrm{T}\mathbf{A}= \mathbf{V}\Sigma^\mathrm{T}\mathbf{U}^\mathrm{T}\mathbf{U}\mathbf{\Sigma}\mathbf{V}^\mathrm{T}$$

thus

$$\mathbf{A}^\mathrm{T}\mathbf{A}=\mathbf{V}\Sigma^2\mathbf{V}^\mathrm{T}$$

Right Singular Vectors:

$$\mathbf{A}\mathbf{A}^\mathrm{T}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^\mathrm{T}(\mathbf{U}\mathbf{\Sigma}\mathbf{V}^\mathrm{T})^\mathrm{T}$$

thus

$$\mathbf{A}\mathbf{A}^\mathrm{T}=\mathbf{U}\Sigma^2\mathbf{U}^\mathrm{T}$$

Hence,

$$\mathbf{A}\mathbf{v}=\sigma \mathbf{u}$$ $$\mathbf{A}^\mathrm{T}\mathbf{u}=\sigma\mathbf{v}$$

Then, we can derive the following:

$$\mathbf{A}^\mathrm{T}\mathbf{A}\mathbf{v}=\sigma^2\mathbf{I}\mathbf{v}$$

We can solve this equation as follows:

$$|\mathbf{A}^\mathrm{T}\mathbf{A}\mathbf{v} - \sigma^2\mathbf{I}\mathbf{v}| = 0$$

Therefore, the singular value is described as:

$$\sqrt{\lambda}$$

where

$$\det (\mathbf{B}-\lambda\mathbf{I})=0$$

Assume that we have two series of data that describe expression patterns of stress response genes.

data1:

data2:

Using the raw data, it is impossible to visualize as the data dimension is higher than 3. Thus, try to reduce the dimensionality.

As the first step, let us visualize the singular values for each dimension.

Since the cumulative contribution rate exceeds 90% up to three dimensions, we will project the data into a three-dimensional space for analysis.

Upon projecting the data into a three-dimensional space, it becomes evident that the gene expression patterns of data1 and data2 are distinctively separated.

The specific growth rate is a coefficient that represents how much one cell divides per unit of time.

Under constant conditions, with the cell mass denoted as X (g), the specific growth rate as μ ($\text{time}^{-1}$), and time as t($h$), the following equation holds:

$$\frac{dX}{dt}=\mu X$$

thus

$$\ln \frac{X}{X_0}=\mu \Delta t$$

where

$X_0$ (g) is the initial cell mass.

Here, we calculate the specific growth rate ($\text{time}^{-1}$) of cells using the cell concentration (OD600) at each time point during the exponential growth phase.

Here is an example of sequential OD600 data of an Escherichia coli strain for 20 hours.

The exponential phase in this case is where OD600 is in between 0.01 and 0.10, therefore we only focus on the area.

Assume that the cell math icreases exponentially at the phase, we set a fitting model given by:

$$X = e^{\mu t + \theta}$$

is also written as:

$$\ln X = \mu t + \theta$$

thus we use a normal equation:

$$\begin{pmatrix}\mu&\theta \end{pmatrix} = (\mathbf{W}^\mathrm{T}\mathbf{W})^{-1} \mathbf{W}^\mathrm{T}\mathbf{f}\in \mathbb{R}^2$$

where

$$\mathbf{f} = (\ln X_1\cdots \ln X_n)^\mathrm{T}\in \mathbb{R}^n$$

$$\mathbf{W}=\begin{pmatrix} t_1&1&1\\ \vdots &\vdots&\vdots\\ t_n&1&1 \end{pmatrix}$$

With the model, we obtained the result where the specific growth rate is:

$$1.53 (\text{h}^{-1})$$

and the fitted growth curve is as shown below.

In this section, we calculate Oligonucleotide Melting Temperature, denoted as Tm, with Nearest Neighbors method.

With Nearest Neighbors method, Tm is calculated as:

$$T_{m} = \frac{\Delta H}{A + \Delta S + R\ln \frac{C}{4}} - 273.15+16.6\log[Na^{+}]$$

, where

Tm = Melting temperature in ℃

$\Delta H$ = enthalpy change in $\frac{kcal}{mol}$

A = constant of -0.0108 $\frac{kcal}{K・mol}$

$\Delta S$ = entropy change in $\frac{kcal}{K・mol}$

R = gas constant of 0.00199 $\frac{kcal}{K・mol}$

C = oligonucleotide concentration in $M$ (as $mol/L$)

$[Na^{+}]$ = Sodium ion conc. in $M$ (as $mol/L$)

and the parameters are already emprically determined as shown:

1. Create a virtual environment for python3

python3 -m venv venv

2. Activate the venv

source venv/bin/activate

3. Leave the environment

deactivate

pip freeze > requirements.txt

pip install -r requirements.txt