Tcl统计配体附近5埃内的氨基酸，在200帧内出现的次数

lxhpdx

1. set sel [atomselect top "same residue as {protein within 5 of resname MOL}"]  
   
2. for {set i 0} {$i<200} {incr i} {  
   
3.     $sel frame $i
   
4.     $sel update  
   
5.     set residue_ids [$sel get resid]   
   
6.     set unique_residue_ids [lsort -unique $residue_ids]   
   
7.     if {[llength $unique_residue_ids] > 0} {  
   
8.         puts "Frame: $i   unique residue IDs: $unique_residue_ids"  
   
9.     } else {  
   
10.         puts "Frame: $i   No residues selected"  
    
11.     }   
    
12. }

复制代码

这个脚本可以正确统计每帧出现的残基ID如下图

，于是乎我想统计一下每个残基在200帧中出现了几次，每帧只计数一次。

1. array set residue_counts {}   
   
2. set sel [atomselect top "same residue as {protein within 5 of resname MOL}"]  
   
3. for {set i 0} {$i<200} {incr i} {  
   
4.     $sel frame $i
   
5.     $sel update   
   
6.     set residue_ids [$sel get resid]  
   
7.     set unique_residue_ids [lsort -unique $residue_ids]  
   
8.       
   
9.     foreach resid $unique_residue_ids {  
   
10.         if {[info exists residue_counts($resid)]} {  
    
11.          
    
12.             incr residue_counts($resid)  
    
13.         } else {  
    
14.             
    
15.             set residue_counts($resid) 1  
    
16.         }  
    
17.     }  
    
18. }  
    
19. foreach resid [array names residue_counts] {  
    
20.     puts "Residue ID: $resid, Total Count: $residue_counts($resid)"  
    
21. }  
    
22.   

复制代码

然后统计出来的残基次数就出现了大于200的情况。

不知道是代码哪里写错了，请各位老师指点。

lxhpdx

1. # 选择与名为"MOL"的残基距离在5埃之内的蛋白质原子
   
2. set sel [atomselect top "same residue as {protein within 5 of resname MOL}"]
   
3. 
   
4. # 循环处理200帧的数据
   
5. for {set i 0} {$i < 200} {incr i} {
   
6.     # 设置当前帧
   
7.     $sel frame $i
   
8.     # 更新选择
   
9.     $sel update
   
10.    
    
11.     # 获取当前帧中选中的残基ID
    
12.     set residue_ids [$sel get resid]
    
13.    
    
14.     # 对残基ID进行去重，得到当前帧中的唯一残基ID
    
15.     set unique_residue_ids [lsort -unique $residue_ids]
    
16.    
    
17.     # 遍历唯一的残基ID
    
18.     foreach resid $unique_residue_ids {
    
19.         # 如果该残基ID已经存在于residue_frames数组中，则添加当前帧号
    
20.         if {[info exists residue_frames($resid)]} {
    
21.             lappend residue_frames($resid) $i
    
22.         } else {
    
23.             # 否则，初始化该残基ID的帧集合并添加当前帧号
    
24.             set residue_frames($resid) [list $i]
    
25.         }
    
26.     }
    
27. }
    
28. 
    
29. # 输出每个残基ID以及其出现的帧数
    
30. foreach resid [array names residue_frames] {
    
31.     set frame_count [llength [lsort -unique $residue_frames($resid)]]
    
32.     puts "Residue ID: $resid, Total Frames: $frame_count"
    
33. }

复制代码

lxhpdx

lxhpdx 发表于 2024-6-12 18:53

问题已解决