How to print predefined sequence in python
<p>I am trying to print predefined sequence from pdb input file in python but I am not getting expected result. I am new in python, and I have also import directory but its not working. not showing anything (unable to find error). Its just running without any output. </p>

<pre><code>import os

os.chdir('C:\Users\Vishnu\Desktop\Test_folder\Input')


for path, dirs, pdbfile in os.walk('/C:\Users\Vishnu\Desktop\Test_folder\Input'):
for line in pdbfile:
    if line[:6] != "HETATM":
        continue
    chainID = line[21:22]
    atomID = line[13:16].strip()
    if chainID not in ('A', 'B', 'C', 'D', 'E', 'F', 'G', 'H'):
        continue
    if atomID not in ('C4B', 'O4B', 'C1B', 'C2B', 'C3B'):
        continue
    with open('C:\Users\Vishnu\Desktop\Test_folder\Input', 'r') as fh:
        new = [line.rstrip() for line in fh]
    with open('C:\Users\Vishnu\Desktop\Test_folder\Output', 'w') as fh:
        [fh.write('%s\n' % line) for line in new]
        fh.write((line.rstrip()))
</code></pre>

<p>Expected output:</p>

<pre><code>HETATM 3788  C4B NAI A 302      52.695  15.486   8.535  1.00 57.28           C  
HETATM 3789  O4B NAI A 302      52.258  14.631   7.456  1.00 56.26           O  
HETATM 3794  C1B NAI A 302      53.348  13.816   7.022  1.00 53.44           C 
HETATM 3792  C2B NAI A 302      54.537  14.748   7.190  1.00 50.93           C  

HETATM 3789  O4B NAI A 302      52.258  14.631   7.456  1.00 56.26           O  
HETATM 3794  C1B NAI A 302      53.348  13.816   7.022  1.00 53.44           C 
HETATM 3792  C2B NAI A 302      54.537  14.748   7.190  1.00 50.93           C 
HETATM 3790  C3B NAI A 302      54.225  15.525   8.465  1.00 52.99           C  

HETATM 3794  C1B NAI A 302      53.348  13.816   7.022  1.00 53.44           C 
HETATM 3792  C2B NAI A 302      54.537  14.748   7.190  1.00 50.93           C 
HETATM 3790  C3B NAI A 302      54.225  15.525   8.465  1.00 52.99           C  
HETATM 3788  C4B NAI A 302      52.695  15.486   8.535  1.00 57.28           C  
HETATM 3792  C2B NAI A 302      54.537  14.748   7.190  1.00 50.93           C 
HETATM 3790  C3B NAI A 302      54.225  15.525   8.465  1.00 52.99           C  
HETATM 3788  C4B NAI A 302      52.695  15.486   8.535  1.00 57.28           C  
HETATM 3789  O4B NAI A 302      52.258  14.631   7.456  1.00 56.26           O  

HETATM 3790  C3B NAI A 302      54.225  15.525   8.465  1.00 52.99           C  
HETATM 3788  C4B NAI A 302      52.695  15.486   8.535  1.00 57.28           C  
HETATM 3789  O4B NAI A 302      52.258  14.631   7.456  1.00 56.26            O  
HETATM 3794  C1B NAI A 302      53.348  13.816   7.022  1.00 53.44           C 
</code></pre>

<p>same format for B chain also.</p>

<p>How to print predefined sequence? line [21:22] is there chain ID, chain ID may be A to H. How to define A to H chain ID? </p>

<p>I am unable to print in sequence, can any one let me know how to print predefined sequence in python?</p>

<p><strong>After Answer:</strong></p>

<p>I have updated above code with below code: </p>

<pre><code>n = 4
for chain, atoms in d.items():
    for atom, line in atoms.items():
        for i in range(len(atom)-n+1):
            for j in range(n):
                print d[chain][atomIDs[i+j]]
            print
</code></pre>

<p>I want to extend two more paragraph but not getting expected output</p>

<p>Here are my comments all combined into an answer:</p>

<pre><code>with open('1AHI.pdb') as pdbfile:
    for line in pdbfile:
        if line[:6] != "HETATM":
            continue
        chainID = line[21:22]
        atomID = line[13:16].strip()
        if chainID not in ('A', 'B'):
            continue
        if atomID not in ('C4B', 'O4B', 'C1B', 'C2B', 'C3B'):
            continue
        ## Either:
        print(line, end='')
        ## Or:
        print(line.rstrip(), end='\n')
        ## Or if Python2.x:
        print line.rstrip()
</code></pre>

<p>My first lines of code was written more than 10 years ago parsing PDB files. Don't despair. You have a long and beautiful journey ahead of you.</p>

<p>P.S. I think mmCIF is to prefer over PDB these days... Make sure you read the specifications for both file formats.</p>

<hr>

<p>I have updated the answer, but please be aware this site is for solving specific problems and not for other people to do the work for you. It is generally looked down on.</p>

<pre><code>d = {}
chainIDs = ('A', 'B',)
atomIDs = ('C4B', 'O4B', 'C1B', 'C2B', 'C3B', 'C4B')
with open('1AHI.pdb') as pdbfile:
    for line in map(str.rstrip, pdbfile):
        if line[:6] != "HETATM":
            continue
        chainID = line[21:22]
        atomID = line[13:16].strip()
        if chainID not in chainIDs:
            continue
        if atomID not in atomIDs:
            continue
        try:
            d[chainID][atomID] = line
        except KeyError:
            d[chainID] = {atomID: line}

n = 4
for chainID in chainIDs:
    for i in range(len(atomIDs)-n+1):
        for j in range(n):
            print d[chainID][atomIDs[i+j]]
        print
</code></pre>

python python-2.7 python-3.x bioinformatics biopython