荧光寿命显微镜的原理及其相图的Python绘制（Phasor Plot for FLIM-FRET）

卢旺收录于类别生物信息和系列 Python Scripts

2022-10-13 2022-10-13 约 1858 字预计阅读 4 分钟

首先简单说明荧光寿命

荧光寿命是指荧光分子在被激发后，能量从激发态回落至基态所需时间的一半（$t_{1/2}$），是荧光分子的固有物理性质（在没有其他因素干扰时，通常是一个常数），是荧光分子除了在波长、偏振、强度之外的另一可以用以挖掘的角度

荧光寿命的测量需要使用荧光寿命显微镜（FLIM），目前最常用的测量方法是时间相关单光子记数法（Time-Correlated Single-Photon Counting ，TCSPC），该方法是在短时间激活荧光分子后，连续一段时间测量光强（光子个数），从而绘制出荧光强度的衰减曲线，并在此基础上拟合指数函数

激活的频率是需要适应测量荧光分子的寿命，因为需要等待一段时间进行测量，这个时间太长和太短都不好，生物荧光分子的寿命通常在12ns以下，因而多使用80MHz的频闪激光作为光源
在实际的图像采集过程中，是在每个像素格子内激活进行测量，为了得到足够的光子数进行数据拟合，通常会Frame Repetitiors 10次

对数据的拟合使用指数衰减的模型，其中$\tau$为荧光寿命，$i$是为了区分不同像素位置，$n$指代10次重复，$A$是最大强度，$Bkgr$是背景强度，$IRF$是仪器响应函数 $$ y(t) = ( \sum_{i=0}^{n} A[i] e^{-\frac{t}{\tau[i]}} + Bkgr ) ※ IRF $$ 最后，数据拟合后得到的效果如下图所示，每个像素点都有一个对应的寿命值，需要划定一个信号区域（ROI）求取平均来得到最终的寿命值，该值需要通过卡方检测

如果生物样品中含有多种荧光分子（混合着不同的荧光寿命），也可以修改模型为多个指数衰减，拟合多个寿命值

关于相图（Phasor Plot）

简单说，相图就是对荧光寿命的数据做傅里叶变换，拆分出频域空间的信号，然后将每一个点的频率信号值画在单位圆上。假设其衰减曲线的信号为指数函数，则可以推到出其在频域空间中横纵轴的计算公式为：

$$ g_{i,j}(\omega) = \frac {\int_0^\infty I_{i,j}(t) cos(\omega t) dt}{\int_0^\infty I_{i,j}(t)dt} $$

$$ s_{i,j}(\omega) = \frac {\int_0^\infty I_{i,j}(t) sin(\omega t) dt}{\int_0^\infty I_{i,j}(t)dt} $$

两者间也必然存在如下关系，说明标准的样品，其荧光信号应该落在这个半圆上

$$ s_{i,j}^2 + (g_{i,j} - \frac12)^2 = \frac14 $$

在频域空间，信号的叠加符合向量加法。因而当荧光拍摄出多个组分时，其信号可能落在单位半圆内部，此时可以将其做向量分解，找到组成其信号的可能组分。

使用Python手绘荧光寿命的Phasor Plot

因为FLIM的数据分析会占据获取有效数据的一多半时间，而显微镜配套的分析处理软件是商业化不公开的，为了节省机时，故将FLIM采集的数据导出为标准格式.ptu，自己回来处理数据；通过指数拟合获得寿命的代码网络上已经有很多了，甚至有imageJ插件直接使用，但更为合理的Phasor分析方法却没有可行的方案，看着算法并不是很困难就自己动手了

1、引包，需要引入FLIM标准格式.ptu读取模块、计算模块numpy、绘图模块matplotlib和seaborn，PTUreader是来自于Github上pull下来的库后微调得到的（原代码也可以直接使用）

1
2
3
4
5
6
7


from readPTU_FLIM import PTUreader
from matplotlib import pyplot as plt
from numpy import fft
import numpy as np
import seaborn as sns
from tqdm import tqdm
sns.set(style='white', palette='muted', color_codes=True)

2、读入.ptu格式文件，并展示图像；

1
2
3
4
5


fpath = '211015Image007.ptu'
f  = PTUreader(fpath, print_header_data = False)
flim_data_stack, intensity_image = f.get_flim_data_stack()
plt.imshow(intensity_image)
plt.colorbar()

使用flim_data_stack.shape可以看到数据的结构，本次数据为(512, 512, 2, 138)，即图像采集尺寸为512×512，包含两个通道（mNG和DAPI），以及在两次脉冲激活间记录了138个时间点信息（80MHz的激光脉冲理论上应该间隔12.5ns，在数据的头文件中f.head['MeasDesc_Resolution']会包含这一数值）

3、可以将所有像素的光子数加起来，观察这个指数衰减图像

1
2
3


plt.figure(figsize=(15,6))
x = np.sum(np.sum(flim_data_stack[:,:,0,:], axis=1), axis=0)
plt.plot(x)

4、通过快速傅里叶变换，计算绘制Phasor图

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40


K = flim_data_stack.shape[3]+1
T = f.head['MeasDesc_Resolution']
omega = 2*np.pi/(K*T)

semicircle = [[], []]
for i in np.arange(0, 100, 0.01):
    semicircle[0].append(1/(1+(i)**2))
    semicircle[1].append(i/(1+(i)**2))

result = []
for i,x in tqdm(enumerate(flim_data_stack)):
    for j,y in enumerate(x):
        d = y[0][20:]
        dd = d/sum(d)
        if sum(d)<=100:
            continue
        ft = np.fft.fft(dd)
        result.append([ft.real[1], -ft.imag[1], intensity_image[i, j]])

step = flim_data_stack.shape[3] - 20
t = np.linspace(0,step*T, step)

ruler = []
for tau in range(1,11):
    tau = tau*1e-9
    y = np.exp(t/tau)
    ynorm = y/sum(y)
    ft = np.fft.fft(ynorm)
    ruler.append([ft.real[1], ft.imag[1]])

g = [i[0] for i in result]
s = [i[1] for i in result]
gr = [i[0] for i in ruler]
sr = [i[1] for i in ruler]
plt.scatter(g, s)
plt.scatter(gr, sr, marker='o', color='r')
# plt.xlim([0,1])
# plt.ylim([0,1])
plt.gca().set_aspect(1)
plt.plot(semicircle[0], semicircle[1], color='black', linewidth=1)

使用seaborn绘制出的、好看一点的密度图，看起来和软件输出的有点像了

1
2
3
4
5
6
7
8


result = np.array(result)
plt.figure(figsize=(10,10))
plt.xlim([0,1])
plt.ylim([0,1])
plt.gca().set_aspect(1)
plt.plot(semicircle[0], semicircle[1], color='black', linewidth=1)
sns.kdeplot(result[:,0], result[:,1], cmap='Reds', shade=True)
plt.scatter(gr, sr, marker='o', color='b')