build-pdb-hhm-db

#!/bin/sh

set -e

function num_cpus {
  if [ -r /proc/cpuinfo ]; then
    NPROC=$(cat /proc/cpuinfo | grep '^processor' | wc -l | tr -d ' ')
  else
    NPROC=1
  fi
  echo $NPROC
}

function msg {
  echo -e $* >&2
}

function usage {
  msg "Usage: `basename $0` [flags] new-db-name hhsuite-db pdb-dir"
  msg
  msg "--help"
  msg "\tShow this help message."
  msg "--cpu num-cpus"
  msg "\tThe number of CPUs for HHsuite to use."
  exit 1
}

num_cpus=`num_cpus`

while true; do
  case "$1" in
    -cpu|--cpu)
      num_cpus=$2
      shift 2
      ;;
    -h|-help|--help)
      usage
      ;;
    -*|--*)
      msg "Invalid flag $1."
      msg
      usage
      ;;
    *)
      break
      ;;
  esac
done

if [ $# != 3 ]; then
  usage
fi

pdb_hhm_db="$1" # The new DB we're creating.
hhsuite_db="$2"  # The DB we're using to generate MSAs.
pdb_dir="$3"     # The directory containing the PDB files.
transient="$pdb_dir/transient" # A temp dir to store fasta/a3m/hhm files.

if [ -d "$pdb_hhm_db" ]; then
  msg "$pdb_hhm_db already exists. Quitting..."
  exit 1
fi
mkdir -p $pdb_hhm_db/pdb

rm -rf "$transient"
mkdir "$transient"

msg "Generating a FASTA file for each PDB file..."
for f in "$pdb_dir"/*.pdb "$pdb_dir"/*.ent "$pdb_dir"/*.ent.gz; do
  if [ ! -f "$f" ]; then
    continue
  fi

  basef=$(basename "$f")
  name=""
  if [[ "$f" = *.pdb ]]; then
    name="${basef%*.pdb}"
  elif [[ "$f" = *.ent ]]; then
    name="${basef%*.ent}"
  elif [[ "$f" = *.ent.gz ]]; then
    name="${basef%*.ent.gz}"
  else
    msg "Unknown PDB extension: $f"
    exit 1
  fi
  pdb2fasta --chain ${name:4} "$f" "$transient/$name.fasta" 2>&1
done

msg "Copying PDB files..."
for f in "$pdb_dir"/*.{pdb,ent,ent.gz}; do
  if [ ! -f "$f" ]; then
    continue
  fi
  cp -a "$f" "$pdb_hhm_db"/pdb/
done

msg "Building a multiple-sequence alignment for each sequence."
msg "This may take a while. Go grab a drink."
glob="$transient"/'*.fasta'
multithread.pl \
  "$glob" \
  'hhblits -i $file -d '"$hhsuite_db"' -n 2 -mact 0.35 -oa3m $name.a3m' \
  --cpu $num_cpus 2>&1

msg "Performing secondary structure annotation on each MSA."
glob="$transient"/'*.a3m'
multithread.pl \
  "$glob" \
  'addss.pl $file' \
  --cpu $num_cpus 2>&1

msg "Cleaning up messes made by hhsuite's 'addss.pl' script."
for f in "$transient"/*.a3m; do
  clean-a3m "$f"
done

msg "Generating an HMM for each MSA."
glob="$transient"/'*.a3m'
multithread.pl \
  "$glob" \
  'hhmake -i $file -pcm 4 -pca 2.5 -pcb 0.5 -pcc 1.0 -gapb 1.0 -gapd 0.15 -gape 1.0 -gapf 0.6 -gapg 0.6 -gapi 0.6' \
  --cpu $num_cpus 2>&1

msg "Building the HHblits database."
hhblitsdb.pl \
  -o "$pdb_hhm_db/$(basename "$pdb_hhm_db")" \
  -ia3m "$transient" \
  -ihhm "$transient" \
  -cpu $num_cpus