Analyze chair/boat conformation from MD simulation trajectory. Part2 Scripts.

Step 1. Atom C1 to C7 are the atom of the ring.

for j in {1..100};do
rm -rf res${j}
mkdir res${j}
cd res${j}
cat > rmsd1.in << EOF
parm ../../*prmtop
trajin ../../*.dcd
EOF

for i in {1..200};do
cat >> rmsd1.in << EOF
reference ../../*.dcd ${i} [${i}]
symmrmsd :${j}@C1,C2,C3,C4,C5,C6,C7 ref [${i}] out res${j}_${i}.dat
EOF
 done

cd ..
done

 for ((j=1;j<=100;j++));do
  cd res${j}
  cpptraj -i rmsd1.in > rmsd1.log
  cd ..
 done

Step 2. Grep is much faster than compare one by one.

unset a b
for j in {1..100};do
cd res${j}
rm -f rmsd_comparison.dat
b[0]=0
for i in {1..200};do
        m=`echo ${b[*]} | grep -w $i | wc -l`
        if [ $m == 1 ];then continue; fi
        a=(`grep -E '0\.0|0\.1' res${j}_${i}.dat | awk '{print $1}'`) #rmsd<0.2
        echo $i,${a[*]} >> rmsd_comparison.dat
        b+=(${a[*]})
done
unset a b
cd ..
done

Step 3. This is a little complicated. You need to modify the PRMTOP file. And make sure every snapshot (PDB format) has he same atom order as the modified PRMTOP requires. The process convert PDB to DCD is not necessary.

mkdir total_rmsd
rm -f total_rmsd/frame_list.txt
#get standard frame of each residues
for i in {1..100};do
 cd res${i}
 a=(`cat rmsd_comparison.dat | cut -f1 -d"," | xargs echo`)
 m=`echo ${a[*]} | tr ' ' ','`
cpptraj ../../*prmtop << EOF
trajin ../../*.dcd
strip !(:${i}@C1,C2,C3,C4,C5,C6,C7)
trajout frame.pdb pdb append onlyframes $m
EOF
 cd ..
  for j in ${a[*]};do
   echo res${i}_frame${j} >> total_rmsd/frame_list.txt
  done
done

#combine pdb of ring, arrange order of atoms
for i in {1..100};do
 cd res${i}
 ((m=`grep -w C7 frame.pdb | wc -l`))
 for ((j=1;j<=$m;j++));do
  for k in C7 C6 C4 C2 C5 C3 C1;do
   grep -w ${k} frame.pdb | sed -n ${j}p >> frame2.pdb
  done
  echo ENDMDL >> frame2.pdb
 done
 cd ..
done

cat > total.pdb << EOF
CRYST1  50.000  50.000  50.000  90.00  90.00  90.00               1
EOF
for i in {1..100};do
 cat res${i}/frame2.pdb >> total.pdb
done

mv total.pdb total_rmsd/

#generate prmtop file for the ring
cd res100
cpptraj ../../*prmtop << EOF
parmstrip !:100@C1,C2,C3,C4,C5,C6,C7
parmwrite out ring.prmtop
EOF
mv ring.prmtop ../total_rmsd/
cd ..

#convert pdb to dcd for the further calculation
cd total_rmsd
cpptraj ring.prmtop << EOF
trajin total.pdb
center :1 mass
trajout total_rmsd.dcd dcd
EOF

#calculate RMSD
cat > rmsd.in << EOF
parm ring.prmtop
trajin total_rmsd.dcd
EOF
((m=`cat frame_list.txt | wc -l`))
for ((i=1;i<=$m;i++));do
cat >> rmsd.in << EOF
reference total_rmsd.dcd $i [$i]
symmrmsd :1@C1,C2,C3,C4,C5,C6,C7 ref [${i}] out rmsd_${i}.dat
EOF
done
cpptraj -i rmsd.in

Step 4. rmsd_comparison.dat and rmsd_comparison2.dat contains information between reference and related frame for each residues.

cd total_rmsd
((line=`cat frame_list.txt | wc -l`))
nl frame_list.txt > nl_frame_list.txt
unset a b
rm -f rmsd_comparison.dat rmsd_comparison2.dat
b[0]=0
for ((i=1;i<=$line;i++));do
        m=`echo ${b[*]} | grep -w $i | wc -l`
        if [ $m == 1 ];then continue; fi
        a=(`grep -E '0\.0|0\.1' rmsd_${i}.dat | awk '{print $1}'`) #rmsd<0.2
        echo $i,${a[*]} >> rmsd_comparison.dat
        echo $i,${a[*]} >> rmsd_comparison2.dat
        for j in ${a[*]};do
                grep -w ${j} nl_frame_list.txt >> rmsd_comparison.dat
        done
        b+=(${a[*]})
done
unset a b

Step 5. Extract final ref to view onlyframe.pdb in VMD. 1,2,5,9 should be the frame of ref in your systems.

parm ring.prmtop
trajin total_rmsd.dcd
trajout onlyframe.pdb pdb onlyframes 1,2,5,9 

Step 6. Analysis. Assuming 1,2 are boat conformation. 5 is chair conformation, 9 is twist boat/chair conformation. The following result will print the conformation of each residue in 50 frames.

mkdir stat_each_resid
#count total number of boat, chair and twisted-chair in each 5ns (50frames)
for k in 1 51 151;do
((l=$k+49))
for i in {1..100};do
((B=0))
((C=0))
((TC=0))
        cd res${i}; #echo res_${i}
        for ((j=$k;j<=$l;j++));do
                #echo $j
                m=`grep -w ${j} rmsd_comparison.dat | cut -f1 -d"," | head -1`
                n=`grep -w res${i}_frame${m} ../total_rmsd/nl_frame_list.txt | awk '{print $1}'`
                ((a=`grep -w $n ../total_rmsd/rmsd_comparison2.dat | cut -f1 -d"," | tail -1`))
                if [ $a == 1 ] || [ $a == 2 ] ;then ((B++));fi
                if [ $a == 5 ] ;then ((C++));fi
                if [ $a == 9 ] ;then ((TC++));fi
        done
        cd ..
        echo $k to $l >> stat_each_resid/res${i}_stat.dat
        echo B $B >> stat_each_resid/res${i}_stat.dat
        echo C $C >> stat_each_resid/res${i}_stat.dat
       echo TC $TC >> stat_each_resid/res${i}_stat.dat
done

done
Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s