MakeCTFGroupCsh

#!/usr/bin/env csh
#
# This script makes CTF groups from a series of CTF parameter files according to a user-defined number
# of divisions (i.e. defocus values in Angstroms). It consists of two separate parts.
#
# The first part of this script removes astigmatism by averaging over defocusU and defocusV.
# Note that CTF-amplitude envelope and background noise information will be ignored.
# This script creates a file called sorted_defocus.dat where the average defocus
# values are sorted from near-to-focus to far_from-focus.
# Visualization of this file (e.g. in gnuplot or xmgr) helps to determine good
# values for the user-defined division values for the second part.
#
# The second part makes the actual groups based on user-defined division values.
# It also creates the corresponding CTFdat file and the CTF parameter files for the groups
#
# To execute the first part, type on the command line:
#      ./make_ctf_groups.csh 1
# To execute the second partm, type:
#      ./make_ctf_groups.csh 2
#
#
# Edit the following parameters:
#
# FIRST PART
#
# Selfile with CTF parameter files for each micrograph to be taken into account
set param_selfile="good2_ctfparams.sel"
# Selfile with the corresponding selfiles of windowed particles for each micrographs
# Note that both selfiles should have the same number of entries in the same order
set sel_selfile="good2_selfiles.sel"
# Amplitude constrast (Q0) for all groups (should be negative)
set Q0=-0.1
# Voltage (kV) for all groups
set kV=200
# Sampling rate (Ang(pix) for all groups:
set angpix=4.2
# Spherical aberration for all groups:
set sa=2.26
# Visualize results with Xmgr? (1=yes, 0=no)
set have_xmgr=1
#
# SECOND PART
# User-defined divisions (defocus values in Angstroms, from large defocus to smaller ones)
set divisions=" -60000 -50000 -40000 -30000"
# Output root name
set outname="group"
# Name for CTFdat file
set ctfdat="all_images.ctfdat"

# IMPORTANT NOTE:
# To make the CTFdat file, it is assumed that foreach
# *_Periodogramavg.ctfparam parameter file a corresponding *.raw.sel selfile exists
#
########################################
# Dont change anything below this line #
########################################

set execute_part = $1

set na=`wc -l $param_selfile|awk '{print $1}'`
set nb=`wc -l $sel_selfile|awk '{print $1}'`
if ($na != $nb) then
  echo "ERROR: " $param_selfile" and "$sel_selfile" have unequal number of entries"
  exit
endif

# FIRST PART
if ($execute_part == 1) then
 echo " -> Making all CTFs non-astigmatic and perform sorting on defocus values"
 rm -f all_defocus.dat
 foreach ctf (`cat $param_selfile |awk '{print $1}'`)
   set new=`echo ${ctf}"_sym"`
   set avg = `grep defocus ${ctf} | awk -F"=" '{avg=avg+$2} END {print avg/2}'`
   echo "defocusU= "$avg > $new
   echo "defocusV= "$avg >> $new
   echo "Q0= " $Q0 >> $new
   echo "voltage= " $kV >> $new
   echo "spherical_aberration= " $sa >> $new
   echo "sampling_rate= " $angpix >> $new
   echo "K= 1.0" >> $new
   printf '%07.0f\n' $avg >> all_defocus.dat
   xmipp_ctf_profile -i ${new} |grep -v "\#" | awk '{print $2,$4}' > ${new}.profile
 end
 sort -r all_defocus.dat >sorted_defocus.dat
 if ($have_xmgr == 1) then
     echo AUTOSCALE >tt
     echo TITLE \"Sorted defocus values \" >>tt
     echo SUBTITLE \" Use this plot to determine the division points \" >>tt
     echo YAXIS LABEL \" Defocus \(in Angstroms\)\" >>tt
   xmgr sorted_defocus.dat -param tt &
 else
   echo " Now visualize the file sorted_defocus.dat to decide on your division points"
 endif
endif

# SECOND PART
if ($execute_part == 2) then
 rm -f ${outname}_*_ctf.sel ${outname}_*_ctf.param ${outname}_*_ctf.param.profile ${ctfdat}
 set pwd=`pwd`
 set max=`head -n 1 sorted_defocus.dat |awk '{print $1-1}'`
 set min=`tail -n 1 sorted_defocus.dat |awk '{print $1+1}'`
 set oldsplit=${max}
 set gr = 0
 foreach split (`echo $divisions $min`)
  @ gr++
  echo " -> Getting CTF group number " ${gr}
  foreach ctf (`cat $param_selfile |awk '{print $1}'`)
    set new=`echo ${ctf}"_sym"`
    set def=`grep defocusU ${new} |awk -F"=" '{print $NF}'`
    set in=`echo ${def} ${oldsplit} ${split}|awk '{if ($1>=$2) {if ($1<$3) print 1}}'`
    if ($in == 1) then
      echo ${new}" 1" >>${outname}_${gr}_ctf.sel
    endif
  end
  set oldsplit=${split}
 end

 # Make the CTF parameter files for the groups
 set ngr = ${gr}
 set gr = 0
 while (${gr} < ${ngr} )
   @ gr++
   echo " -> Making CTF parameter file for CTF group number " ${gr}
   set sum = 0
   set ii = 0
   foreach mem (`cat ${outname}_${gr}_ctf.sel |awk '{print $1}'`)
     set def=`grep defocusU ${mem} |awk -F"=" '{print $2}'`
     set sum = `echo "\( ${sum} \+ ${def} \)" | bc -l`
     @ ii++
   end
   set avgdef=`echo "\( ${sum} / ${ii} \)" | bc `
   set new=`echo ${outname}_${gr}.ctfparam`
   echo "defocusU= "$avgdef > $new
   echo "defocusV= "$avgdef >> $new
   echo "Q0= " $Q0 >> $new
   echo "voltage= " $kV >> $new
   echo "spherical_aberration= " $sa >> $new
   echo "sampling_rate= " $angpix >> $new
   echo "K= 1.0" >> $new
   xmipp_ctf_profile -i ${new} | grep -v "\#" | awk '{print $2,$4}' > ${new}.profile

   # Make CTFdat file
   set npart = 0
   foreach ctf (`cat group_${gr}_ctf.sel | awk '{print $1}' |sed 's|_sym||'` )
     set nn=`grep -n ${ctf} $param_selfile |awk -F":" '{print $1}'`
     set sel=`head -n ${nn} $sel_selfile |tail -n 1 |awk '{print $1}'`
     if (-e ${sel}) then
       set part = `wc -l ${sel}|awk '{print $1}'`
       @ npart = $npart + $part
     endif
     cat ${sel} |awk -v"ctf=${pwd}/${new}" '{print $1,ctf}'>> ${ctfdat}
   end

   # Output to screen
   echo " --->  $ii micrographs, $npart particles, average defocus = $avgdef"

   # visualize
   if ($have_xmgr) then
     echo AUTOSCALE >tt
     echo TITLE \"${outname}_${gr}_ctf.sel\" >>tt
     echo SUBTITLE \" $ii members, average defocus = $avgdef, number of particles = $npart  \" >>tt
     xmgr ${new}.profile `cat ${outname}_${gr}_ctf.sel| sed 's|_sym 1|_sym.profile|' ` -param tt >& .t
     rm -f .t
   else
     echo "Visualize the following files to analyse the Ctfs of this groups:"
     cat ${outname}_${gr}_ctf.sel| sed 's|_sym 1|_sym.profile|'
     echo "###"
   endif


 end
endif
echo " exiting normally..."

-- Main.SjorsScheres - 28 Nov 2007