计算化学公社

标题: 如何将残基拟合成一条直线 [打印本页]

作者
Author:
1154975925    时间: 2024-9-5 16:07
标题: 如何将残基拟合成一条直线
(, 下载次数 Times of downloads: 4) 该图是一个微管蛋白二聚体(1JFF),其中标红的是编号(222-244)的残基,我现在想要把这图上的点拟合成一条曲线,有没有什么方法能够做到。

作者
Author:
Loading0760    时间: 2024-9-5 16:35
看起来像是同源二聚体?
既然这是一个结构那么就一定存在坐标,把残基的坐标提取出来.得到一个20*3的矩阵,接下来问题就是20个点,求一条拟合直线.
最小二乘法启动,以下代码来源自ChatGPT
  1. import numpy as np

  2. def fit_line_least_squares(points):
  3.     # 计算质心
  4.     centroid = np.mean(points, axis=0)
  5.     # 去中心化
  6.     points_centered = points - centroid
  7.     # 计算协方差矩阵
  8.     cov_matrix = np.dot(points_centered.T, points_centered)
  9.     # 对协方差矩阵进行特征值分解
  10.     _, _, vh = np.linalg.svd(cov_matrix)
  11.     # 取出最大特征值对应的特征向量作为方向向量
  12.     direction = vh[0]
  13.     # 返回质心和方向向量
  14.     return centroid, direction

  15. centroid, direction = fit_line_least_squares(points)
复制代码


你的标题有[GROMACS] 字段,如果是想分析轨迹的话,可以用MDAnalysis库
以下是读取动力学轨迹并提取CYS的代码,编号(222-244)的残基部分同理

  1. import MDAnalysis as mda
  2. u = mda.Universe("npt.gro", "md_0_2_fit_50ns.xtc")
  3. cys_selection = u.select_atoms("resname CYS")
  4. cys_selection.positions #提取坐标
复制代码

作者
Author:
student0618    时间: 2024-9-5 17:20
目的是作图初步比较一下,还是定量分析?

作图的话PyMOL settings -> cartoon -> cylindrical helices 可以把helices显示为圆柱。

分析的话如楼上所说。
作者
Author:
sobereva    时间: 2024-9-5 21:29
可以对螺旋区域计算最长主轴对应的矢量。不需要做拟合

北京科音分子动力学与GROMACS培训班(http://www.keinsci.com/workshop/KGMX_content.html)里一个VMD tcl脚本中的一段可以参考

(, 下载次数 Times of downloads: 5)

(, 下载次数 Times of downloads: 3)
作者
Author:
1154975925    时间: 2024-9-23 17:11
sobereva 发表于 2024-9-5 21:29
可以对螺旋区域计算最长主轴对应的矢量。不需要做拟合

北京科音分子动力学与GROMACS培训班(http://www. ...

那我是不是把这个top alpha改成我图中的残基就可以了,希望sob老师可以解答一下
作者
Author:
sobereva    时间: 2024-9-24 01:44
1154975925 发表于 2024-9-23 17:11
那我是不是把这个top alpha改成我图中的残基就可以了,希望sob老师可以解答一下

top alpha是指整个体系里所有alpha碳。如果你要对指定区段的残基用alpha碳来计算,加上resid选择语句
作者
Author:
1154975925    时间: 2024-10-8 21:47
sobereva 发表于 2024-9-5 21:29
可以对螺旋区域计算最长主轴对应的矢量。不需要做拟合

北京科音分子动力学与GROMACS培训班(http://www. ...

sob老师想请教一下,我用这个代码跑两次模拟,一次加电场,一次不加,能否用这两次模拟输出结果的差值来表现整个蛋白质在电场作用下扭转的角度,这个代码是针对螺旋的,不是很明白选中整个蛋白质的条件下输出结果的意义。
作者
Author:
sobereva    时间: 2024-10-9 04:26
1154975925 发表于 2024-10-8 21:47
sob老师想请教一下,我用这个代码跑两次模拟,一次加电场,一次不加,能否用这两次模拟输出结果的差值来 ...

假定蛋白质是刚性的,主轴矢量随着蛋白质旋转而相应旋转,那就可以
作者
Author:
1154975925    时间: 2024-10-9 08:59
sobereva 发表于 2024-10-9 04:26
假定蛋白质是刚性的,主轴矢量随着蛋白质旋转而相应旋转,那就可以

感谢
作者
Author:
1154975925    时间: 2024-10-11 10:11
sobereva 发表于 2024-10-9 04:26
假定蛋白质是刚性的,主轴矢量随着蛋白质旋转而相应旋转,那就可以

(, 下载次数 Times of downloads: 2) sob老师我这里蛋白质的最长主轴似乎是y轴,如果想要实现我说的功能应该是提取蛋白质y轴的转动惯量Iyy [x2 y2 z2],然后再提取其中的y分量,相应的代码应由set zcomp [lindex [lindex [lindex [measure inertia $sel] 1] 2] 2] 改为set ycomp [lindex [lindex [lindex [measure inertia $sel] 1] 1] 1],不知道是否正确,希望sob老师可以回复一下,谢谢。

作者
Author:
sobereva    时间: 2024-10-11 19:04
1154975925 发表于 2024-10-11 10:11
sob老师我这里蛋白质的最长主轴似乎是y轴,如果想要实现我说的功能应该是提取蛋白质y轴的转动惯量Iyy [x2 ...

用VMD的绘制cylinder的命令把矢量在图上画出来一看便知对不对
作者
Author:
1154975925    时间: 2024-10-17 22:55
sobereva 发表于 2024-10-11 19:04
用VMD的绘制cylinder的命令把矢量在图上画出来一看便知对不对

sob老师好,我输出了我所选择螺旋的[measure inertia $sel],{40.410728454589844 96.76710510253906 69.3785171508789} {{0.9233367443084717 0.26798388361930847 0.27501633763313293} {-0.17684435844421387 -0.33896195888519287 0.9240296483039856} {-0.34084513783454895 0.9018256068229675 0.2655846178531647}},并根据这组数据在图中显示了质心以及最后最后一行即最长惯性主轴方向向量,为了方便查看我将向量×了1000,代码如下draw cylinder {0 0 0} {-340.84513783454895 901.8256068229675 265.5846178531647} radius 2 filled yes,最后在VMD中的效果如图,可是这个最长惯性主轴似乎不能代表这个螺旋上点拟合的直线,是不是这个最长惯性主轴只能求螺旋角度差而不能求螺旋主轴的矢量,希望SOB老师能解答一下,谢谢
(, 下载次数 Times of downloads: 2)

作者
Author:
sobereva    时间: 2024-10-20 10:10
1154975925 发表于 2024-10-17 22:55
sob老师好,我输出了我所选择螺旋的[measure inertia $sel],{40.410728454589844 96.76710510253906 69. ...

一张图看不充分
从当前图上来看蓝色的柱已经合理展现出了螺旋区域的朝向




欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3