align_family.sh

#!/bin/bash -e

export MAX_N_PID_4_TCOFFEE=$(cat /proc/sys/kernel/pid_max)

# ./align_family.sh folder/ <optional: tool1 tool2 ...>
#
# Run tool suite on a given cluster folder.
# Assumes all tools on $PATH and directory structure:
#
# folder/
# 	pdbs/
# 		structureA.pdb
# 		structureB.pdb
#	folder_aa.fa
#	folder_msa.fa
#
# As generated by clean_homstrad.py

THREADS="${THREADS:=1}"
FAMILY=$(basename "$1")
PDB=$(realpath "${1}/pdbs/")
AA="${1}/sequence.fa"
TEMPLATE="${1}/tcoffee.template"

FMT="Command being timed: %C\nUser time (seconds): %U\nSystem time (seconds): %S\nPercent of CPU this job got: %P\nWall clock time (seconds): %e\nAverage shared text size (kbytes): %X\nAverage unshared data size (kbytes): %%D\nAverage stack size (kbytes): %p\nAverage total size (kbytes): %K\nMaximum resident set size (kbytes): %M\nAverage %resident set size (kbytes): %t\nMajor (requiring I/O) page faults: %F\nMinor (reclaiming a frame) page faults: %%R\nVoluntary context switches: %w\nInvoluntary context switches: %c\nSwaps: %W\nFile system inputs: %I\nFile system %outputs: %O\nSocket messages sent: %s\nSocket messages received: %r\nSignals delivered: %k\nPage size (bytes): %Z\nExit %status: %x"

# Executable paths
declare -A paths=(
	[foldmason]=$(command -v "foldmason")
	[caretta]=$(command -v "caretta-cli")
	[matt]=$(command -v "Matt")
	[mtm]=$(command -v "mTM-align")
	[mustang]=$(command -v "mustang-3.2.4")
	[clustalo]=$(command -v "clustalo")
	[famsa]=$(command -v "famsa")
	[mafft]=$(command -v "linsi")
	[muscle]=$(command -v "muscle5")
	[usalign]=$(command -v "USalign")
	[tcoffee]=$(command -v "t_coffee")
)

# Enable/disable status of each tool
declare -A tools=(
	[foldmason]=false
	[caretta]=false
	[matt]=false
	[usalign]=false
	[mtm]=false
	[mustang]=false
	[clustalo]=false
	[famsa]=false
	[mafft]=false
	[muscle]=false
	[tcoffee]=false
	[none]=false
)

# no args: help
# 1 arg: input folder
# >1 arg: input folder, specific tools to run
# throw errors on no args, invalid tool, valid tool but path not found
if [[ $# -eq 0 ]]; then
	echo "./align_family.sh folder/ <optional: tool1 tool2 ...>"
	exit 1
elif [[ $# -eq 1 ]]; then
	for tool in "${!tools[@]}"; do
		if [[ -n paths["$tool"] ]]; then
			tools["$tool"]=true
		else
			echo "${tool} not found on $PATH"
			exit 1
		fi
	done
elif [[ $# -eq 2 && $2 == "none" ]]; then
	: # no op
elif [[ $# -gt 1 ]]; then
	for tool in "${@:2}"; do
		if [[ ! -v tools["$tool"] ]]; then
			echo "Invalid tool: ${tool}. Choose from [${!tools[@]}]"
			exit 1
		fi
		if [[ -n "${paths[$tool]}" ]]; then
			tools["$tool"]=true
		else
			echo "${tool} not found on \$PATH"
			exit 1
		fi
	done
fi

if [[ "${tools[none]}" == false ]]; then
	echo "Running tools:"
	for tool in "${!tools[@]}"; do
		if [[ "${tools[$tool]}" == true ]]; then
			echo "  ${tool}: ${paths[$tool]}"
		fi
	done
fi

# Structure aligners
if [[ "${tools[caretta]}" == true && ! -e "${1}/caretta_results/result.fasta" ]]; then
	if [[ -e "${1}/caretta_results" ]]; then
		rm -r "${1}/caretta_results"
	fi
	/usr/bin/time -o "${1}/caretta.time" -f "${FMT}" "${paths[caretta]}" "$PDB" -t "$THREADS" -o "${1}/caretta_results" --
	sed -i 's/\.pdb//' "${1}/caretta_results/result.fasta"
fi
if [[ "${tools[foldmason]}" == true && ! -e "${1}/foldmason_aa.fa" ]]; then
	/usr/bin/time -o "${1}/foldmason.time" -f "${FMT}" "${paths[foldmason]}" easy-msa \
		"$PDB" "${1}/foldmason" "${1}/foldmason_tmp" --threads $THREADS
fi
if [[ "${tools[foldmason]}" == true && ! -e "${1}/foldmason_refine100_aa.fa" ]]; then
	/usr/bin/time -o "${1}/foldmason_refine100.time" -f "${FMT}" "${paths[foldmason]}" easy-msa \
		"$PDB" "${1}/foldmason_refine100" "${1}/foldmason_refine100_tmp" \
		--refine-iters 100 --refine-seed 48335597 --threads $THREADS --pair-threshold 0
fi
if [[ "${tools[matt]}" == true && ! -e "${1}/matt.fasta" ]]; then
	/usr/bin/time -o "${1}/matt.time" -f "${FMT}" "${paths[matt]}" -o "${1}/matt" $(find "$PDB" -type f) -t "$THREADS" -s0 -b1
	sed -i 's/:.*$//' "${1}/matt.fasta"
	sed -i 's/:.*$//' "${1}/matt_bent.fasta"
fi
if [[ "${tools[mtm]}" == true && ! -e "${1}/mTM_result/result.fasta" ]]; then
	/usr/bin/time -o "${1}/mtmalign.time" -f "${FMT}" "${paths[mtm]}" -i <(find "$PDB" -type f) -outdir "${1}/mTM_result"
	sed -i 's/\.pdb//' "${1}/mTM_result/result.fasta"
fi
if [[ "${tools[usalign]}" == true && ! -e "${1}/usalign.fa" ]]; then
	/usr/bin/time -o "${1}/usalign.time" -f "${FMT}" "${paths[usalign]}" \
		-dir "$PDB" <(find "$PDB" -type f -printf "%f\n") -mm 4 -outfmt 1 > "${1}/usalign.fa"
	sed -i '/^[#$]/d; /^$/d; s/\.pdb.*$//g; s/^>\//>/' "${1}/usalign.fa"
fi
if [[ "${tools[mustang]}" == true && ! -e "${1}/mustang.afasta" ]]; then
	/usr/bin/time -o "${1}/mustang.time" -f "${FMT}" "${paths[mustang]}" -i $(find "$PDB" -type f) -F fasta -o "${1}/mustang"
	sed -i 's/\.pdb//' "${1}/mustang.afasta"
fi
if [[ "${tools[tcoffee]}" == true && ! -e "${1}/3dcoffee.fa" ]]; then
	if [ ! -e "$TEMPLATE" ]; then
		awk -v fo=$(realpath "$PDB") '/^>/ { header=$1; sub(/^>/, "", header); print $1" _P_ "fo"/"header".pdb" }' "$AA" > "$TEMPLATE"
	fi
	/usr/bin/time -o "${1}/3dcoffee.time" -f "${FMT}" "${paths[tcoffee]}" "$AA" \
		-method sap_pair \
		-thread 1 \
		-template_file "$TEMPLATE" \
		-output fasta \
		-outfile "${1}/3dcoffee.fa" \
		-quiet "${1}/3dcoffee.log" \
		-newtree "${1}/3dcoffee.tree"  # otherwise all guide trees are saved to root
	sed -i 's/ _P_.*$//' "${1}/3dcoffee.fa"
fi

# Sequence aligners
if [[ "${tools[clustalo]}" == true && ! -e "${1}/clustalo.fa" ]]; then
	/usr/bin/time -o "${1}/clustalo.time" -f "${FMT}" "${paths[clustalo]}" -i "$AA" -o "${1}/clustalo.fa" --threads "$THREADS" --infmt fasta
fi
if [[ "${tools[famsa]}" == true && ! -e "${1}/famsa.fa" ]]; then
	/usr/bin/time -o "${1}/famsa.time" -f "${FMT}" "${paths[famsa]}" "$AA" "${1}/famsa.fa" -t "$THREADS"
fi
if [[ "${tools[mafft]}" == true && ! -e "${1}/mafft.fa" ]]; then
	/usr/bin/time -o "${1}/mafft.time" -f "${FMT}" "${paths[mafft]}" --thread "$THREADS" "$AA" > "${1}/mafft.fa"
fi
if [[ "${tools[muscle]}" == true && ! -e "${1}/muscle.fa" ]]; then
	/usr/bin/time -o "${1}/muscle.time" -f "${FMT}" "${paths[muscle]}" -align "${AA}" -output "${1}/muscle.fa" -threads "$THREADS"
fi

# Generate LDDT reports
DB="${1}/foldmason_tmp/latest/structures"

compute_lddt () {
	if [[ "${tools[$1]}" == false ]]; then return; fi
	if [[ ! -e "$2" ]]; then return; fi
	if [[ -e "$3" ]]; then return; fi
	"${paths[foldmason]}" msa2lddtreport "$DB" "$2" "$3" --threads "$THREADS"
}

if [[ -e $DB ]]; then
	#echo "Computing LDDT scores"
	if [[ -e "${1}/${FAMILY}_msa.fasta" ]]; then
		compute_lddt "homstrad" "${1}/${FAMILY}_msa.fasta" "${1}/homstrad.html"
	fi
	compute_lddt "foldmason" "${1}/foldmason_aa.fa"              "${1}/foldmason.html"
	compute_lddt "foldmason" "${1}/foldmason_refine100_aa.fa"    "${1}/foldmason_refine100.html"
	compute_lddt "muscle"    "${1}/muscle.fa"                    "${1}/muscle.html"
	compute_lddt "caretta"   "${1}/caretta_results/result.fasta" "${1}/caretta.html"
	compute_lddt "matt"      "${1}/matt.fasta"                   "${1}/matt.html"
	compute_lddt "mtm"       "${1}/mTM_result/result.fasta"      "${1}/mtmalign.html"
	compute_lddt "usalign"   "${1}/usalign.fa"      	     "${1}/usalign.html"
	compute_lddt "mustang"   "${1}/mustang.afasta"       	     "${1}/mustang.html"
	compute_lddt "clustalo"  "${1}/clustalo.fa"                  "${1}/clustalo.html"
	compute_lddt "famsa"     "${1}/famsa.fa"                     "${1}/famsa.html"
	compute_lddt "mafft"     "${1}/mafft.fa"                     "${1}/mafft.html"
	compute_lddt "tcoffee"   "${1}/3dcoffee.fa"                  "${1}/3dcoffee.html"
fi