Recipe: DistMap0

created by jeff101

Profile

• Name: DistMap0
• ID: 101659
• Shared with: Public
• Parent: None
• Children: None
• Created on: December 04, 2015 at 18:31 PM UTC
• Updated on: December 04, 2015 at 18:31 PM UTC
• Description:

  DistMap0c.txt 12/4/15 1220pm noon code

Code

-- jeff101 began DistMap0 12/3/15 354pm
-- last updated 12/4/15 1220pm noon

recipename='DistMap0'

function Main()
local tot=structure.GetCount()
local nbands=band.GetCount()
print(' ')
print(recipename..': Puzzle '..puzzle.GetName())
print('   with ID '..puzzle.GetPuzzleID()..' at '..os.date()) 
print('   has score '..trunc3(current.GetScore())..', '..nbands..' bands, and '..tot..' residues.')
print(' ')
--
-- next part calculates Ca-Ca distance map in dists and dds
--
local dist,dd,dcuts,j,j1,j2,minj1,minj2,maxj1,maxj2
local sss={}
local aas={}
local dists={}
local dds={}
local mindist= -1+0 -- will be min nonzero Ca-Ca distance found
local maxdist= -1+0 -- will be max nonzero Ca-Ca distance found
dcuts={2, 4, 6, 8, 10, 15, 20, 30, 50}
for j1=1,tot do
   dists[j1]={}
   dds[j1]={}
   sss[j1]=structure.GetSecondaryStructure(j1)
   aas[j1]=structure.GetAminoAcid(j1)
end -- for j1
for j1=1,tot do
   for j2=j1,tot do
      if j1==j2 then
         dist=0+0
         dd=getdd(dcuts,dist)
         dists[j1][j2]=dist
         dds[j1][j2]=dd
      else
         dist=getdist(j1,j2,0,0) 
         -- dist should now be distance between alpha-carbons j1 and j2
         if dist>0 then
            if maxdist== -1 or dist>maxdist then
               maxdist=dist
               maxj1=j1
               maxj2=j2
            end -- if maxdist
            if mindist== -1 or dist<mindist then
               mindist=dist
               minj1=j1
               minj2=j2
            end -- if mindist
         end -- if dist
         dd=getdd(dcuts,dist)
         dists[j1][j2]=dist
         dists[j2][j1]=dist
         dds[j1][j2]=dd
         dds[j2][j1]=dd      
      end -- if j1==j2
   end -- for j2
end -- for j1
--
-- next part prints the distance map in dds
--
print('In the alpha-carbon distance map(s) below:')
print('   0 means distance is in range (-Inf,'..dcuts[1]..').')
for j=1,#dcuts-1 do
   print('   '..j..' means distance is in range ['..dcuts[j]..','..dcuts[j+1]..').')
end -- for j
print('   '..#dcuts..' means distance is in range ['..dcuts[#dcuts]..',Inf).')
print(string.format('All nonzero distances range from %.2f to %.2f',mindist,maxdist))
print(string.format('   for the Ca-Ca pairs %d-%d and %d-%d respectively.',minj1,minj2,maxj1,maxj2))    
print(' ')
makechart(sss,aas,dds,dcuts,mindist,maxdist,minj1,minj2,maxj1,maxj2,tot,tot)
print(' ')
if tot>74 then
   makechart(sss,aas,dds,dcuts,mindist,maxdist,minj1,minj2,maxj1,maxj2,tot,74)
   print(' ')
end
if tot>34 then
   makechart(sss,aas,dds,dcuts,mindist,maxdist,minj1,minj2,maxj1,maxj2,tot,34)
   print(' ')
end
if tot>18 then
   makechart(sss,aas,dds,dcuts,mindist,maxdist,minj1,minj2,maxj1,maxj2,tot,18)
   print(' ')
end
makechart(sss,aas,dds,dcuts,mindist,maxdist,minj1,minj2,maxj1,maxj2,tot,9)
end -- Main()

function makechart(sss,aas,dds,dcuts,mindist,maxdist,minj1,minj2,maxj1,maxj2,tot,jdim)
local tmpstr,j100s,j10s,j1s,alleq,ssstr,aastr
local js,j,j1,j2,j100,j10,jrem
js={}
j1s='      '
j10s=j1s
j100s=j1s
alleq=j1s
ssstr=j1s
aastr=j1s
for j=1,jdim do
   js[j]=1+round0((j-1)*(tot-1)/(jdim-1)) -- want j=1 to jdim to give js=1 to tot
   j100=trunc0(js[j]/100)
   jrem=js[j]-j100*100
   j10=trunc0(jrem/10)
   j1=jrem-j10*10
   j100s=(j100s..j100)
   j10s=(j10s..j10)
   j1s=(j1s..j1)
   alleq=(alleq..'|')
   ssstr=(ssstr..sss[js[j]])
   aastr=(aastr..aas[js[j]])
end -- for j
for j1=jdim,1,-1 do
   if jdim==9 then
      tmpstr=string.format('%03d-',js[j1])
   else
      tmpstr=string.format('%s%03d%s-',sss[js[j1]],js[j1],aas[js[j1]])
   end -- if jdim
   for j2=1,jdim do
      tmpstr=string.format('%s%d',tmpstr,dds[js[j1]][js[j2]])
   end -- for j2
   if jdim==9 then
      if j1==9 then
         tmpstr=string.format('%s %s: for all %d amino acids',tmpstr,recipename,tot)
      elseif j1==8 then
         tmpstr=string.format('%s    min Ca-Ca dist=%.2f for %d-%d',tmpstr,mindist,minj1,minj2)
      elseif j1==7 then
         tmpstr=string.format('%s    max Ca-Ca dist=%.2f for %d-%d',tmpstr,maxdist,maxj1,maxj2)
      elseif j1==6 then
         tmpstr=string.format('%s Ca-Ca dist ranges chart key:',tmpstr)
      elseif j1==5 then
         tmpstr=string.format('%s    0 for 0-%d, 1 for %d-%d',tmpstr,dcuts[1],dcuts[1],dcuts[2])
      elseif j1==1 then
         tmpstr=string.format('%s    8 for %d-%d, 9 for %d-Inf',tmpstr,dcuts[8],dcuts[9],dcuts[9])
      else
         j=2*(5-j1) -- j1==4 -> j=2, j1==3 -> j=4, j1==2 -> j=6 
         tmpstr=string.format('%s    %d for %d-%d, %d for %d-%d',tmpstr,j,dcuts[j],dcuts[j+1],j+1,dcuts[j+1],dcuts[j+2])
      end -- if j1
   end -- if jdim
   print(tmpstr)
end -- for j1
if jdim~=9 then
   print(alleq)
   print(aastr)
   print(j100s)
   print(j10s)
   print(j1s)
   print(ssstr)
end -- if jdim
end -- makechart()

-- below returns an integer val from 0 to #dcuts
-- representing what range the distance dist lies in
-- val should be 0 if dist<dcuts[1]             or          (-Inf,dcuts[1])
-- val should be j if dcuts[j]<=dist<dcuts[j+1] or      [dcuts[j],dcuts[j+1])
-- val should be #dcuts if dcuts[#dcuts]<=dist  or [dcuts[#dcuts],Inf)
-- ideally #dcuts is 9 or less
function getdd(dcuts,dist)
local val=0+0
local j
for j=1,#dcuts do
   if dist>=dcuts[j] then
      val=j
   end -- if dist
end -- for j
return val
end -- getdd()

-- below finds the distance between residue r1 atom a1 and residue r2 atom a2
function getdist(r1,r2,a1,a2)
local tot,nuband
local dist= -1+0 -- default value
if a1==a2 and a1==0 then
   dist=structure.GetDistance(r1,r2) 
   -- above finds distance between atom 0's, the alpha-carbons on each residue
else
   -- below finds distance between general atoms on each residue
   tot=band.GetCount()
   nuband=band.AddBetweenSegments(r1,r2,a1,a2)
   if nuband>tot then
      dist=band.GetLength(nuband)
      band.Delete(nuband)
   else
      -- 4/9/14 below happened for a1==a2 and a1==0 when r1 and r2 differed by 1
      -- 5/26/14 below happened for a1==1 and a2==4 when r1 and r2 differed by 1 or less
      print(string.format("ERROR: getdist(%d,%d,%d,%d) failed so returning -1.",r1,r2,a1,a2))
   end -- if nuband
end -- if a1==a2
return dist
end -- getdist()

-- below rounds down to 3 decimal places (nearest 0.001)
function trunc3(numi)
local numo=trunc0(numi*1000)/1000
return numo
end -- trunc3()

-- below rounds to 2 decimal places (nearest 0.01)
function round2(numi)
local numo=round0(numi*100)/100
return numo
end -- round2()

-- below rounds to 1 decimal place (nearest 0.1)
function round1(numi)
local numo=round0(numi*10)/10
return numo
end -- round1()

-- below rounds val to the nearest integer to get res
function round0(val)
local res=math.floor(val)
if val-res>=0.5 then
   res=res+1
end -- if val
return res
end -- round0()

-- below rounds val down to the nearest integer to get res
function trunc0(val)
local res=math.floor(val)
return res
end -- trunc0()

Main()

Comments

[jeff101]

DistMap0 finds the distances between every pair of alpha-carbons in the protein
and makes several maps of these distances. It gives output like below:

DistMap0: Puzzle 973: Post-CASP Tc853: Server Model Contacts with ID 998376
   at 12/03/15 20:33:23 has score 11250.722
   0 bands and 152 residues.
 
In the alpha-carbon distance map below:
   0 means distance is in range (-Inf,2).
   1 means distance is in range [2,4).
   2 means distance is in range [4,6).
   3 means distance is in range [6,8).
   4 means distance is in range [8,10).
   5 means distance is in range [10,15).
   6 means distance is in range [15,20).
   7 means distance is in range [20,30).
   8 means distance is in range [30,50).
   9 means distance is in range [50,Inf).
All nonzero distances range from 3.75 to 44.71
   for the Ca-Ca pairs 35-36 and 58-121 respectively.

L152d-6666778756777645787776554776567760
L147d-7776667567777774566677775254566506
L143k-8876556678778876654567776556665057
L138g-8887778778888887777567775564550567
L134t-7777778778888877776345665676305666
E129v-7766677767777766776555552675035655
E125v-8877778777788876777577775350564646
L120n-8887778788888887777778887305776557
E115n-8877677677778875666677775033665527
E111t-7766677667777765777566650575255654
L106d-5767778767777757887765405787567775
L102g-5545677767777666777754046787567775
H097k-7755567777777777776530456787557776
H093t-7766567778777777775403566777546677
H088e-8877677778778877775045775675535567
L083t-8876456788777887630556777677667467
L079d-8876545688667887503777787677777568
E074d-8776655477567775056777787677777657
E070d-7655657355566650577777675576677645
H065k-5556778624776505788777656787678774
E060p-5146777655664056788877677888788876
H056f-7455667566450466777877777888788877
L051n-7654555577405676667777777788788777
L047s-7655555266044675567777777787788777
E042y-6667778620676545788887777787788876
E038n-6556777602676525788777666787677765
E033t-7655545066255663467777776677777657
H028e-8776540578557787556777787788788678
L024i-7763204477556775545766777777777567
H019i-7753025577556776654655676677677567
H015k-7540336567545665666765576777677666
H010l-5304567556555455777765466787678776
E006t-5035777656664156788877577888778876
L001g-0557778766777557888877557888778876
      ||||||||||||||||||||||||||||||||||
      gtlkiietnysnfpkdddtetkgdtnnvvtgkdd
      0000000000000000000000111111111111
      0011122334455667778899001122233445
      1605948382716050493837261505948372
      LEHHHLHEEELLHEHEELLHHHLLEELEELLLLL
 
152-675785650 152 amino acids in all
133-767884605 min Ca-Ca dist=3.75 for 35-36
114-867876066 max Ca-Ca dist=44.71 for 58-121
095-757770645 Ca-Ca dist ranges chart key:
077-868807788 0 for 0-2, 1 for 2-4
058-677087887 2 for 4-6, 3 for 6-8
039-670787775 4 for 8-10, 5 for 10-15
020-707765667 6 for 15-20, 7 for 20-30
001-076687876 8 for 30-50, 9 for 50-Inf

The final map above is small enough to appear in the Recipe Output,
so it can be included in a snapshot of the protein when a puzzle ends.

Below is a sample snapshot made right after the DistMap0 run above.
You can click on the below image for a more detailed view.

[jeff101]

[jeff101]

It would be neat if we could use html code within the Recipe Output window.
Then I could use characters like █ ▓ ▒ ░ or █ ▀ ▄ ▌ ▐ ▙ ▛ ▜ ▟ ▚ ▞ ▖ ▘ ▝ ▗
from http://www.w3schools.com/charsets/ref_utf_block.asp
to make the above distance maps.

If you know how to use html code within the Recipe Output window,
please post here or at http://fold.it/portal/node/2001554 about it.

Thanks!

[jeff101]

Distance and Contact Maps show the protein in a way that does not depend on the protein's orientation in 3D space.

If two structures for the same protein look different in 3D, it could be (1) because they are identical structures viewed from different directions or (2) because they are truly different structures. Looking at their Distance or Contact Maps can tell if the two structures are actually the same or truly different.

[jeff101]

Below is a Matlab image made from the full Distance Map made by the DistMap0 run above.
This map is a 152 x 152 grid of integers from 0 to 8. Note the resemblance between this
map and the Contact Map shown above. Click on the image below for more detail.

[jeff101]

p.247 of the book "Proteins: Structures and Molecular Properties"
by Thomas E Creighton, Macmillan, 1993, ISBN 071677030X or 9780716770305
shows Contact Maps for the homologous proteins alpha-chymotrypsin and elastase
that look very similar to each other.

https://books.google.com/books?id=hu8T_kI1LrkC&pg=PA247
https://books.google.com/books?id=hu8T_kI1LrkC

[jeff101]

DistMap1 (http://fold.it/portal/recipe/101679)
is a more elaborate version of DistMap0 that
can make *.eps image files like the Matlab image above,
uses a variety of shaded character sets to show distances
in the Recipe Output window and scriptlog.default.xml file,
and can make statistical charts suitable for input into
spreadsheet programs like Excel.