技术领域
[0001] 本
发明涉及
图像处理领域，尤其涉及一种冷冻电镜三维密度图分辨率检测方法。
背景技术
[0002] 作为一项新兴技术，冷冻电镜(Cryo-EM)技术近年来迅速发展并取得了卓越的成就，越来越多大分子通
过冷冻电镜技术重构得到亚
纳米级三维结构。结构
生物领域传统的三维重构技术，如x-射线三维重构技术和
核磁共振三维重构技术(NMR)，各自有其固有
缺陷，如x-射线技术只能对晶体进行研究，不能对非晶体进行三维重构；核磁共振技术不能对大分子量(超过100KDa)的分子进行三维重构。相比之下，冷冻电镜技术具有显著的优势：(1)对研究样本没有结晶要求，对晶体和非晶体都可进行三维重构；(2)对研究样本没有分子量限制，可对大分子量的生物分子进行三维重构；(3)对样品进行快速冷冻，使样品中的
水处于
玻璃态，生物分子非常接近生理状态；(4)样本中的生物分子处于生理过程的不同阶段，可研究生理过程中生物分子的动态演化过程等。冷冻电镜技术显著的优越性使其具有广阔的发展前景。
[0003] 利用冷冻电镜技术进行三维重构的主要方法是单颗粒重构(SPR)。单颗粒重构的流程如下：(1)冷冻样品制备；(2)电镜成像采集；(3)电镜图像处理，包括：CTF校正，颗粒挑选等；(4)三维重构，包括对颗粒进行分类、平均，构建初始模型，
迭代优化最终模型等。单颗粒重构得到生物分子的三维密度图后，对重构模型的
质量进行评估，即分辨率检测，具有重要的意义：有助于比较重构
算法的优劣，促进不同技术的融合交流，促进冷冻电镜技术的发展，促进结构生物学的发展等。
[0004] 目前常用的分辨率检测方法有傅里叶壳相关(Fourier shell correlation,FSC) 方法和谱
信噪比(Spectral signal-to-noise ratio,SSNR)方法。通过FSC方法检测分辨率需要将颗粒数据集随机分成两个半数据集，对这两个半数据集进行三维重构，然后计算两个半数据集重构结果在傅里叶空间的相关性，得到FSC 曲线，通过
阈值截断，得到最终的全局分辨率。使用FSC方法进行分辨率检测存在以下缺点：(1)该方法检测的实际上是重构算法的可重复性，与传统的分辨率检测意义不同，阈值的选取没有直观的含义，目前依然存在较大争议； (2)该方法易受噪声影响，由于高频部分的噪声多，该方法在计算时，往往会将噪声的相关性作为颗粒相关性来计算，出现过拟合现象；(3)该方法需要将颗粒数据集分成两个半数据集并分别进行三维重构，半数据集的重构质量低于全部数据集的重构质量，据此计算出的分辨率与全部数据集的重构结果分辨率不符，且计算过程冗余。
[0005] SSNR方法计算三维密度图在频域的
信号与噪声的功率比，然后通过阈值截断，计算全局分辨率。与FSC方法相比，计算分辨率时使用的阈值有明确的含义，如常取1作为阈值，此时信号功率等于噪声功率。现有的计算SSNR 的方法是将三维密度图进行再投影，投影空间的图像被认为是去噪后的图像，用原始图像减去投影图像得到噪声图像，然后使用由纯噪声三维密度图得到的噪声衰减因子，对投影图像的信号功率和噪声功率进行修正，最后由修正后的投影图像信号功率比修正后的噪声功率得到三维密度图的谱信噪比。使用现有方法计算SSNR是在二维空间考虑三维空间的信噪比，不具有直观性，计算结果依赖于再投影过程，且计算过程中需要原始数据集。
[0006] 冷冻电镜技术正处于快速发展阶段，且具有广阔的应用前景，而现有的冷冻电镜三维密度图分辨率检测方法都存在较大缺陷，因此需要提出一种新的分辨率检测方法。
发明内容
[0007] 针对上述
现有技术中的不足，本发明提供一种冷冻电镜三维密度图分辨率检测方法，在三维空间中直接计算三维密度图的谱信噪比，不需要再投影过程，具有简单、直观、操作便捷的优点。
[0008] 为了实现上述目的，本发明提供一种冷冻电镜三维密度图分辨率检测方法，包括步骤：
[0009] S1：输入一冷冻电镜三维密度图和所述冷冻电镜三维密度图的一纯噪声三维密度图；
[0010] S2：根据所述冷冻电镜三维密度图生成一掩膜集，所述掩膜集包括n个大小不同的掩膜，n为大于零的整数；
[0011] S3：利用各所述掩膜分别对所述纯噪声三维密度图进行分割分别获得所述纯噪声三维密度图的与各所述掩膜对应的一第一膜内数据和一第一膜外数据；
[0012] S4：根据所述第一膜内数据和所述第一膜外数据计算各所述掩膜对应的一第一膜内噪声谱功率和一第一膜外噪声谱功率；
[0013] S5：根据各组所述第一膜内噪声谱功率和所述第一膜外噪声谱功率计算获得一线性参数集λp；
[0014]
[0015] 其中，0≤p≤n，且p为整数， 为第p掩膜所对应的第一膜内噪声谱功率，为第p掩膜所对应的第一膜外噪声谱功率；λp为线性拟合参数；
[0016] S6：利用各所述掩膜分别对所述冷冻电镜三维密度图进行分割分别获得所述冷冻电镜三维密度图的与各所述掩膜对应的一第二膜内数据和一第二膜外数据；
[0017] S7：根据所述第二膜内数据、所述第二膜外数据和所述线性参数集计算获得各掩膜所对应的第二膜内数据谱信噪比曲线图；
[0018] S8：根据所述第二膜内数据谱信噪比曲线图获得一分辨率集，并根据所述分辨率集绘制一分辨率集三维曲面图；
[0019] S9：根据所述分辨率集三维曲面图计算获得所述冷冻电镜三维密度图的一全局分辨率值。
[0020] 优选地，所述掩膜采用二值掩膜，所述二值掩膜包括一掩膜区域，所述掩膜区域内的值为1，所述掩膜区域外的值为0。
[0021] 优选地，在所述S2步骤中，通过在一透射
电子显微镜图像处理
软件中设置并调整一三维密度图的密度等值面阈值参数和一向外扩展的高斯层数参数，处理所述冷冻电镜三维密度图，获得所述掩膜集。
[0022] 优选地，所述S3步骤进一步包括步骤：
[0023] 将所述纯噪声三维密度图在三维空间每个坐标上的值与当前所述掩膜相应坐标上的值相乘，获得当前掩膜所对应的所述第一膜内数据；
[0024] 将所述纯噪声三维密度图减去所述第一膜内数据，获得当前掩膜所对应的所述第一膜外数据；
[0025] 重复步骤，直至获得所有所述掩膜所对应的所述第一膜内数据和所述第一膜外数据。
[0026] 优选地，所述S4步骤进一步包括步骤：
[0027] 将当前所述掩膜所对应的所述第一膜内数据通过傅里叶变换转换到一第一
频率空间；将当前所述掩膜所对应的所述第一膜外数据通过傅里叶变换转换到一第二频率空间；
[0028] 分别在所述第一频率空间和所述第二频率空间内建立一球壳模型shell(s)， s为所述球壳模型的半径；
[0029] 根据一公式(2)计算当前所述掩膜所对应的所述第一膜内噪声谱功率和所述第一膜外噪声谱功率：
[0030]
[0031] 其中，RPS(s)为谱功率；当RPS(s)为所述第一膜内噪声谱功率时，Ks为所述第一频率空间内所述球壳模型上的点的坐标，所述Ns为所述第一频率空间内所述球壳模型上的点的数量，M(Ks)为所述第一膜内数据的傅里叶变换；
[0032] 当RPS(s)为所述第一膜外噪声谱功率时，Ks为所述第二频率空间内所述球壳模型上的点的坐标，所述Ns为所述第二频率空间内所述球壳模型上的点的数量，M(Ks)为所述第一膜外数据的傅里叶变换。
[0033] 优选地，所述S7步骤进一步包括步骤：
[0034] S71：将当前所述掩膜所对应的所述第二膜内数据通过傅里叶变换转换到一第三频率空间；将当前所述掩膜所对应的所述第二膜外数据通过傅里叶变换转换到一第四频率空间；
[0035] S72：分别在所述第三频率空间和所述第四频率空间内建立一球壳模型 shell(s)，s为所述球壳模型的半径；
[0036] S73：根据所述公式(2)计算当前所述掩膜所对应的一第二膜内信号和噪声加和的谱功率和一第二膜外噪声谱功率：
[0037] 当RPS(s)为所述第二膜内信号和噪声加和的谱功率时，Ks为所述第三频率空间内所述球壳模型上的点的坐标，所述Ns为所述第三频率空间内所述球壳模型上的点的数量，M(Ks)为所述第二膜内数据的傅里叶变换；
[0038] 当RPS(s)为所述第二膜外噪声谱功率时，Ks为所述第四频率空间内所述球壳模型上的点的坐标，所述Ns为所述第四频率空间内所述球壳模型上的点的数量，M(Ks)为所述第二膜外数据的傅里叶变换；
[0039] S74：利用一公式(3)、所述线性参数集λp和所述第二膜外噪声谱功率估算所述冷冻电镜三维密度图的一第二膜内噪声谱功率：
[0040] RPSnp(s)＝λp·RPSop(s) (3)；
[0041] 其中，RPSnp(s)为第p掩膜所对应的所述第二膜内噪声谱功率；RPSop(s) 为第p掩膜所对应的所述第二膜外噪声谱功率；
[0042] S75：将所述第二膜内信号和噪声加和的谱功率减去所述第二膜内噪声谱功率计算获得所述冷冻电镜三维密度图的一第二膜内信号谱功率RPSsp(s)：
[0043] RPSsp(s)＝RPSip(s)-RPSnp(s) (4)；
[0044] 其中，RPSip(s)为第p掩膜所对应的所述第二膜内信号和噪声加和的谱功率；RPSnp(s)为第p掩膜所对应的所述第二膜内噪声谱功率；
[0045] S76：根据所述第二膜内信号谱功率和所述第二膜内噪声谱功率计算一第二膜内数据谱信噪比，绘制当前掩膜所对应的所述第二膜内数据谱信噪比曲线图；
[0046] S77：重复步骤S71～S76直至绘制获得所有掩膜所对应的所述第二膜内数据谱信噪比曲线图。
[0047] 优选地，所述第二膜内数据谱信噪比满足公式(5)：
[0048]
[0049] 其中，SSNRp(s)为第p掩膜所对应的所述第二膜内数据谱信噪比。
[0050] 优选地，所述绘制当前掩膜所对应的所述第二膜内数据谱信噪比曲线图步骤中：将所述第二膜内数据谱信噪比取对数，获得所述第二膜内数据谱信噪比曲线图的纵坐标，将所述第二膜内数据谱信噪比所对应的所述球壳模型的半径作为所述第二膜内数据谱信噪比曲线图的横坐标，形成所述第二膜内数据谱信噪比曲线图。
[0051] 优选地，所述S8步骤进一步包括步骤：
[0052] 设定一阈值；
[0053] 在每一所述第二膜内数据谱信噪比曲线图内截取每一所述第二膜内数据谱信噪比曲线上纵坐标与所述阈值数值相等的至少一待选坐标点，并选取所述待选坐标点中横坐标最小的一所述待选坐标点作为一目标坐标点；
[0054] 计算所述目标坐标点横坐标的倒数，获得各所述掩膜所对应的所述第二膜内数据的分辨率，各所述掩膜所对应的所述第二膜内数据的分辨率形成所述分辨率集；
[0055] 将各所述掩膜所对应的所述分辨率作为z轴，将各所述掩膜所对应所述三维密度图的密度等值面阈值参数和所述向外扩展的高斯层数参数分别作为 x轴和y轴，形成所述分辨率集三维曲面图，所述分辨率集三维曲面图包括一分界线。
[0056] 优选地，所述S9步骤进一步包括步骤：
[0057] 计算所述分界线上各点的分辨率数值的均值；
[0058] 寻找所述分界线上各点中分辨率数值与所述均值最为接近的一点作为判别点；
[0059] 将所述判别点的分辨率作为所述全局分辨率值。
[0060] 本发明由于采用了以上技术方案，使其具有以下有益效果：
[0061] 本发明的冷冻电镜三维密度图分辨率检测方法，不需要用两个半数据集重构结果的相似性来估计分辨率，计算分辨率时不需要使用原始颗粒数据集，不需要进行投影过程，可直接在三维空间中计算三维密度图的谱信噪比，具有简单、直观、操作便捷的优点。
具体实施方式
[0062] 下面给出本发明的较佳
实施例，并予以详细描述，使能更好地理解本发明的功能、特点。
[0063] 本发明实施例的一种冷冻电镜三维密度图分辨率检测方法，包括步骤：
[0064] S1：输入一冷冻电镜三维密度图和冷冻电镜三维密度图的一纯噪声三维密度图。
[0065] 例如：可通过调用EMDataBank
数据库中ID code是EMD-8119的数据包。数据包包含一真实图像和两个只使用了一半数据重构得到的半数据集重构结果，其中冷冻电镜三维密度图可采用该数据包中预存的真实图像，可通过将两个半数据重构得到的半数据集重构结果相减获得纯噪声三维密度图；
[0066] S2：根据冷冻电镜三维密度图生成一掩膜集，掩膜集包括n个大小不同的掩膜，n为大于零的整数。
[0067] 其中，掩膜采用二值掩膜，二值掩膜包括一掩膜区域，掩膜区域内的值为 1，掩膜区域外的值为0。
[0068] 在S2步骤中，通过在一透射电子显微镜图像处理软件EMAN2中设置并调整一三维密度图的密度等值面阈值参数(threshold)和一向外扩展的高斯层数参数(nshells)，处理冷冻电镜三维密度图，获得掩膜集。
[0069] 当三维密度图的密度等值面阈值参数的数值越大或者向外扩展的高斯层数参数的数值越小时，生成的掩膜尺寸越小。通过使用chimera软件观察冷冻电镜三维密度图可以得到生成一个合理掩膜时三维密度图的密度等值面阈值参数的取值上下限，当颗粒图像被噪声淹没到几乎看不见时对应的值为三维密度图的密度等值面阈值参数的下限，当几乎完全看不到噪声时对应的值为三维密度图的密度等值面阈值参数的上限，步长一般选择使得上下限之间有十个取值的大小。向外扩展的高斯层数参数一般取值为2～11，步长为1。
[0070] S3：利用各掩膜分别对纯噪声三维密度图进行分割分别获得纯噪声三维密度图的与各掩膜对应的一第一膜内数据和一第一膜外数据。
[0071] 其中，S3步骤进一步包括步骤：
[0072] 将纯噪声三维密度图在三维空间每个坐标上的值与当前掩膜相应坐标上的值相乘，获得当前掩膜所对应的第一膜内数据；
[0073] 将纯噪声三维密度图减去第一膜内数据，获得当前掩膜所对应的第一膜外数据；
[0074] 重复步骤，直至获得所有掩膜所对应的第一膜内数据和第一膜外数据。
[0075] S4：根据第一膜内数据和第一膜外数据计算各掩膜对应的一第一膜内噪声谱功率和一第一膜外噪声谱功率。
[0076] 其中，S4步骤进一步包括步骤：
[0077] 将当前掩膜所对应的第一膜内数据通过傅里叶变换转换到一第一频率空间；将当前掩膜所对应的第一膜外数据通过傅里叶变换转换到一第二频率空间；
[0078] 分别在第一频率空间和第二频率空间内建立一球壳模型shell(s)，s为球壳模型的半径；
[0079] 根据一公式(2)计算当前掩膜所对应的第一膜内噪声谱功率和第一膜外噪声谱功率：
[0080]
[0081] 其中，RPS(s)为谱功率；当RPS(s)为第一膜内噪声谱功率时，Ks为第一频率空间内球壳模型上的点的坐标，Ns为第一频率空间内球壳模型上的点的数量，M(Ks)为第一膜内数据的傅里叶变换；
[0082] 当RPS(s)为第一膜外噪声谱功率时，Ks为第二频率空间内球壳模型上的点的坐标，Ns为第二频率空间内球壳模型上的点的数量，M(Ks)为第一膜外数据的傅里叶变换。
[0083] S5：根据各组第一膜内噪声谱功率和第一膜外噪声谱功率计算获得一线性参数集λp；
[0084]
[0085] 其中，0≤p≤n，且p为整数， 为第p掩膜所对应的第一膜内噪声谱功率，为第p掩膜所对应的第一膜外噪声谱功率；λp为线性拟合参数。
[0086] S6：利用各掩膜分别对冷冻电镜三维密度图进行分割分别获得冷冻电镜三维密度图的与各掩膜对应的一第二膜内数据和一第二膜外数据；
[0087] S7：根据第二膜内数据、第二膜外数据和线性参数集计算获得各掩膜所对应的第二膜内数据谱信噪比曲线图。
[0088] 其中，S7步骤进一步包括步骤：
[0089] S71：将当前掩膜所对应的第二膜内数据通过傅里叶变换转换到一第三频率空间；将当前掩膜所对应的第二膜外数据通过傅里叶变换转换到一第四频率空间；
[0090] S72：分别在第三频率空间和第四频率空间内建立一球壳模型shell(s)，s 为球壳模型的半径；
[0091] S73：根据公式(2)计算当前掩膜所对应的第二膜内信号和噪声加和的谱功率和第二膜外噪声谱功率：
[0092] 当RPS(s)为第二膜内信号和噪声加和的谱功率时，Ks为第三频率空间内球壳模型上的点的坐标，Ns为第三频率空间内球壳模型上的点的数量，M(Ks) 为第二膜内数据的傅里叶变换；
[0093] 当RPS(s)为第二膜外噪声谱功率时，Ks为第四频率空间内球壳模型上的点的坐标，Ns为第四频率空间内球壳模型上的点的数量，M(Ks)为第二膜外数据的傅里叶变换；
[0094] S74：利用一公式(3)、线性参数集λp和第二膜外噪声谱功率估算冷冻电镜三维密度图的一第二膜内噪声谱功率：
[0095] RPSnp(s)＝λp·RPSop(s) (3)；
[0096] 其中，RPSnp(s)为第p掩膜所对应的第二膜内噪声谱功率；RPSop(s)为第 p掩膜所对应的第二膜外噪声谱功率；
[0097] S75：将第二膜内信号和噪声加和的谱功率减去第二膜内噪声谱功率计算获得冷冻电镜三维密度图的一第二膜内信号谱功率RPSsp(s)：
[0098] RPSsp(s)＝RPSip(s)-RPSnp(s) (4)；
[0099] 其中，RPSip(s)为第p掩膜所对应的第二膜内信号和噪声加和的谱功率； RPSnp(s)为第p掩膜所对应的第二膜内噪声谱功率；
[0100] S76：根据第二膜内信号谱功率和第二膜内噪声谱功率计算所述第二膜内数据谱信噪比，绘制当前掩膜所对应的一第二膜内数据谱信噪比曲线图；
[0101] 具体地，将第二膜内数据谱信噪比取对数log10SSNR，获得第二膜内数据谱信噪比曲线图的纵坐标，将第二膜内数据谱信噪比所对应的球壳模型的半径作为第二膜内数据谱信噪比曲线图的横坐标，形成第二膜内数据谱信噪比曲线图；
[0102] S77：重复步骤S71～S76直至绘制获得所有掩膜所对应的第二膜内数据谱信噪比曲线图。
[0103] 第二膜内数据谱信噪比满足公式(5)：
[0104]
[0105] 其中，SSNRp(s)为第p掩膜所对应的第二膜内数据谱信噪比。
[0106] S8：根据第二膜内数据谱信噪比曲线图获得一分辨率集，并根据分辨率集绘制一分辨率集三维曲面图。
[0107] 其中，S8步骤进一步包括步骤：
[0108] 设定一阈值，本实施例中，设定阈值log100.334≈-0.48；
[0109] 在每一第二膜内数据谱信噪比曲线图内截取每一第二膜内数据谱信噪比曲线上纵坐标与阈值数值相等的至少一待选坐标点，并选取待选坐标点中横坐标最小的一待选坐标点作为一目标坐标点；
[0110] 计算目标坐标点横坐标的倒数，获得各掩膜所对应的第二膜内数据的分辨率，各掩膜所对应的第二膜内数据的分辨率形成分辨率集；
[0111] 将各掩膜所对应的分辨率作为z轴，将各掩膜所对应三维密度图的密度等值面阈值参数和向外扩展的高斯层数参数分别作为x轴和y轴，形成分辨率集三维曲面图，分辨率集三维曲面图包括一分界线。
[0112] S9：根据分辨率集三维曲面图计算获得冷冻电镜三维密度图的一全局分辨率值。
[0113] 优选地，S9步骤进一步包括步骤：
[0114] 计算分界线上各点的分辨率数值的均值；
[0115] 寻找分界线上各点中分辨率数值与均值最为接近的一点作为判别点；
[0116] 将判别点的分辨率作为全局分辨率值。
[0117] 以上结合实施例对本发明进行了详细说明，本领域中普通技术人员可根据上述说明对本发明做出种种变化例。因而，实施例中的某些细节不应构成对本发明的限定，本发明将以所附
权利要求书界定的范围作为本发明的保护范围。
