assembly.sh

#!/bin/bash 

usage()
{
	echo "Usage: $0 [-b 1|2][-D nr|u90][-d Fasta_Dir]|[-f Fasta_file][-o Output_dir][-I Work_ID][-m model][-N Decoy_num][-n Output_num]"
	echo "	-b	rounds of blast when using Psipred - 1 or 2 (default: 1)"
	echo "	-D	the sequence database - nr or u90 (default: u90)"
	echo "	-d	the directory as input"
	echo "	-f	a fasta file as input"
	echo "	-I	work ID (when using -f option)"
	echo "	-m	the model type used for model evaluation - G(General) or D(Dwan) (default: G)"
	echo "	-N	number of PDB files generated (threading method) (default: 50)"
	echo "	-n	number of PDB files for output (default: 5)"
	echo "	-o	the directory for output (default: $PWD)"
	exit 2
}

if [ "$#" -eq "0" ]; then
	usage
fi
#put 'path' file and this script in the same directory
pathfile=${0%/*}/softwarePath.txt
if [ ! -r $pathfile ]; then
	echo "Cannot read file: $pathfile"
	exit 2
fi
IFS_OLD=$IFS
IFS=$'\n'
for line in `cat $pathfile`
do
	case "$line" in
	"Modeller")
		count=1;;
	"nr database")
		count=2;;
	"u90 database")
		count=3;;
	*)
		line=`echo $line | sed 's/^[[:space:]]*//'`
		case "$count" in
		"1")	## Modeller
			Modeller=$line/bin;;
		"2")	## nr database
			dbname1=$line/nr;;
		"3")	## u90 database
			dbname2=$line/u90;;
		esac;;
	esac
done
unset line
IFS=$IFS_OLD

packagePath=${0%/*}
cd $packagePath
packagePath=$PWD
Jpredbin=$packagePath/FR-t5/bin
Jpreddata=$packagePath/FR-t5/data
Fea37Pack=$packagePath/MEFTop

if [ ! -f $Modeller/modpy.sh ]; then
	echo "Path error! - Modeller: ${Modeller%/*}">&2
	exit 2
fi
if [ ! -d $Jpredbin ]||[ ! -d $Jpreddata ]; then
	echo "Path error! - Jpreditor: ${Jpredbin%/*}">&2
	exit 2
fi
if [ ! -d $Fea37Pack ]; then
	echo "Path error! - Fea37 Package: $Fea37Pack">&2
	exit 2
fi
if [ ! -d ${dbname1%/*} ]&&[ ! -d ${dbname2%/*} ]; then
	echo -e "Path error! - Sequence database:\n\t$dbname1\n\t$dbname2">&2
	exit 2
fi

NNconExecSh=$Fea37Pack/Program/nncon/bin/predict_sssa_cm.sh
if [ ! -x $NNconExecSh ]; then
	echo "File error! - nncon: $NNconExecSh">&2
	exit 2
fi

DSSPexec=$Fea37Pack/Program/DSSP/dsspcmbi
if [ ! -x $DSSPexec ]; then
	echo "File error! - dssp: $DSSPexec">&2
	exit 2
fi

HCpath=$Fea37Pack/Program/HC
HCexec=$Fea37Pack/Program/HC/HCsta
if [ ! -x $HCexec ]; then
	echo "File error! - Hydrophobic Core: $HCexec">&2
	exit 2
fi

Fea37exec=$Fea37Pack/Program/Fea37/GenFea
if [ ! -x $Fea37exec ]; then
	echo "File error! - Fea37: $Fea37exec">&2
	exit 2
fi

svmscaleExec=$Fea37Pack/Program/libsvm-weights-3.11/svm-scale
svmpredictExec=$Fea37Pack/Program/libsvm-weights-3.11/svm-predict
if [ ! -x $svmscaleExec ]||[ ! -x $svmpredictExec ]; then
	echo -e "File error! - libsvm:\n\t$svmscaleExec\n\t$svmpredictExec">&2
	exit 2
fi

##default
[ -n "$dbname2" ]&&dbname=$dbname2||dbname=$dbname1
mainpath=$PWD
pwd=$PWD
cutoff=6
dnum=50
onum=5
fflag=0		# read one file(1);read a directory(2)
bflag=1		# rounds of BLAST when using Psipred
mflag="General"		# the model type you would use

while getopts b:D:d:f:I:m:N:n:o: OPT
do
	case "$OPT" in
	#[-b 1|2] rounds of BLAST when using Psipred
	"b")
		bflag=$OPTARG
		if [ $bflag -ne 1 ]&&[ $bflag -ne 2 ]; then
			echo "!argError--please use '-b 1' or '-b 2'">&2
			exit 2
		fi;;
	#[-D nr|u90] select sequence library (NR database OR u90 database)
	#	default: u90 database
	"D")
		case "$OPTARG" in
		"nr"|"NR")
			if [ -z "$dbname1" ]; then
				echo "!argError--can't find nr database">&2
				exit 2
			fi
			dbname=$dbname1;;
		"u90"|"U90")
			if [ -z "$dbname2" ]; then
				echo "!argError--can't find u90 database">&2
				exit 2
			fi
			dbname=$dbname2;;
		*)	
			echo "!argError--please use '-D nr' or '-D u90'">&2
			exit 2;;
		esac;;
	#[-d directory] read a directory which contains fasta files
	"d")
		if [ -d "$OPTARG" ] && [ -n "$(ls $OPTARG/*.fasta)" ]; then	
			cd $OPTARG
			in_dir=$PWD
			fflag=2
			cd $pwd
		else
			echo "!argError--invalid file path: $OPTARG">&2
			exit 2
		fi;;

	#[-f file] the name of one fasta file
	"f")
		if [ -f "$OPTARG" ]; then
			fflag=1
			in_ID=${OPTARG##*/}
			in_ID=${in_ID%.fasta}
			in_id=$OPTARG
		elif [ -f "$OPTARG.fasta" ]; then
			fflag=1
			in_ID=${OPTARG##*/}
			in_id=$OPTARG.fasta
		else
			echo "!argError--invalid work ID: $OPTARG">&2
			exit 2
		fi;;
	#[-I ID] work ID of one fasta file
	"I")
		ID_tmp=$OPTARG;;
	#[-m model] the model type used for model evaluation
	"m")
		case "$OPTARG" in
		"G"|"General")
			mflag="General";;
		"D"|"Dwan")
			mflag="Dwan";;
		*)	
			echo "!argError--please use '-m G' or '-m D'">&2
			exit 2;;
		esac;;
	#[-N dnum] number of PDB files generated by THREADING method (default:50)
	"N")
		dnum=$OPTARG;;
	#[-n onum] number of PDB files for output (default:5)
	"n")
		onum=$OPTARG;;
	#[-o Output_dir]
		# Three directories will be made
		# Input/  fasta files
		# Work/   all intermediate file
		# Output/ the final result
	"o")
		mkdir -p $OPTARG
		if [ $? -ne 0 ]; then
			echo "!argError--invalid output path: $OPTARG">&2
			exit 2
		else
			cd $OPTARG
			mainpath=$PWD
		fi
		cd $pwd;;
	"?")	
		usage;;
	esac
done

#model & par path
Model=`ls $Fea37Pack/Model/$mflag/*.model | head -n1`
Par=`ls $Fea37Pack/Model/$mflag/*.par | head -n1`
if [ ! -r $Model ]||[ ! -r $Par ]; then
	echo -e "File error! - Fea37:\n\t$Model\n\t$Par">&2
	exit 2
fi
##===model===## :
if [ "$mflag" == "General" ]; then
	k=2.15 #( Z-score <= $cutoff )
else
	k=1.60 #( Z-score <= $cutoff )
fi

mkdir -p $mainpath/Input
mkdir -p $mainpath/Work
mkdir -p $mainpath/Output

Predic() {
	
	ID=$1
	timelog=$mainpath/$ID.time
	exec 3>$timelog
	#Working Dir
	
	#workdir, nearly all temporary files are in this folder
	subpath=$mainpath/Work/$ID
	decoypath=$mainpath/Work/$ID/decoy
	#final path, final output files( 5 pdb model) will be copied to this folder
	finalpath=$mainpath/Output/$ID
	
	mkdir $subpath
	mkdir $decoypath
	mkdir $finalpath	
	
	#Get 1 line sequence
	cd $mainpath
	echo `date +%T`"	Start...$ID">&3
	echo -n `date +%T`"	Get 1 line sequence...">&3
	
	#Make fasta file one line
	inFile=$mainpath/Input/$ID.fasta
	if [ ! -r $infile ]; then
		echo "Cannot read fasta file: $infile">&2
		return 1
	fi
	lock=0
	IFS_OLD=$IFS
	IFS=$'\n'
	for line in `cat $inFile`
	do	
		if [ "${line:0:1}" == ">" ]; then
			echo "$line"
			lock=1
		elif [ -n "$line" ]; then
			echo -n "$line" | tr [a-z] [A-Z]
			[ "$lock" -ge "1" ] && lock=$[$lock+1]
		elif [ -z "$line" ]&&[ "$lock" -ge "2" ]; then
			break
		fi
	done > $decoypath/$ID.fasta
	IFS=$IFS_OLD

	if [ ! -f $decoypath/$ID.fasta ]; then
		echo "Cannot find $ID.fasta in working directory">&2
		return 1
	fi
	
	if [ "$lock" -lt "2" ]; then
		echo "$ID.fasta is not a regular fasta file">&2
		rm -f $decoypath/$ID.fasta
		return 1
	fi
	unset lock
	echo "done">&3
	
	echo -n `date +%T`"	Psipred: Predict the secondary structure...">&3
	cd $decoypath
	
	#Psipred
	if [ $bflag -eq 1 ]; then
		$Jpredbin/blastpgp -Q $ID.psi -b 0 -j 3 -h 0.001 -d $dbname -i $ID.fasta -C temp.chk > temp.blast
		#echo `date +%T`"	blast-round...done">&3
	else
		$Jpredbin/blastpgp -Q $ID.psi -j 3 -h 0.001 -d $dbname -i $ID.fasta -C temp1.chk > temp1.blast
		#echo `date +%T`"	blast-round1...done">&3
		$Jpredbin/blastpgp -v 5000 -b 0 -j 3 -h 0.001 -d $dbname -i $ID.fasta -C temp.chk > temp.blast
		#echo `date +%T`"	blast-round2...done">&3
	fi
	echo temp.chk > $ID.psitmp.pn
	echo $ID.fasta > $ID.psitmp.sn
	$Jpredbin/makemat -P $ID.psitmp
	cp temp.mtx $ID.psitmp.mtx
	$Jpredbin/psipred $ID.psitmp.mtx $Jpreddata/weights.dat $Jpreddata/weights.dat2 $Jpreddata/weights.dat3 > $ID.ss
	$Jpredbin/psipass2  $Jpreddata/weights_p2.dat 1 1.0 1.0 $ID.ss2 $ID.ss > $ID.horiz
	rm -f *psitmp* temp* *.ss *.ss2
	echo "done">&3
	Threading
}

Threading(){
	
	#Threading
	echo "@@ Threading Method @@">&3
	echo -n `date +%T`"	Threading prediction...">&3
	$Jpredbin/JTHREADER_MIAOZHICHAO $ID.horiz $ID.psi $dnum $ID.aln
	echo "done">&3
	
	#Modelling
	echo -n `date +%T`"	Run modeller...">&3
	for i in `seq 1 $dnum`
	do
		Tid=`cat $ID.aln | sed -n "$(($i*3-2))p" | sed 's/^>//g'`
		Tid=${Tid%% *}
		seq1=`cat $ID.aln | sed -n "$(($i*3-1))p"`
		seq2=`cat $ID.aln | sed -n "$(($i*3))p"`
		
		Ali=$ID-$i.ali
		touch $Ali
		echo -ne ">P1;$Tid\n" > $Ali
		echo -ne "structureX:$Tid:FIRST:@:LAST:@: : : :\n" >> $Ali
		echo -ne "$seq1*\n\n" >> $Ali
		echo -ne ">P1;$ID\n" >> $Ali
		echo -ne "sequence: $ID: : : : : : : :\n" >> $Ali
		echo -ne "$seq2*\n\n" >> $Ali
		
		Py=$ID-$i.py
		touch $Py
		echo -ne "from modeller import *\nfrom modeller.automodel import *\n\n" > $Py
		echo -ne "log.verbose()\nenv = environ()\n\n" >> $Py
		echo -ne "env.io.atom_files_directory =		 ['.','$decoypath']\n" >> $Py
		echo -ne "a = automodel(env,\n" >> $Py
		echo -ne "alnfile  = '$Ali',\n" >> $Py
		echo -ne "knowns   = '$Tid',\n" >> $Py
		echo -ne "sequence = '$ID')\n" >> $Py
		echo -ne "a.starting_model= 1\na.ending_model  = 1\na.make()" >> $Py
		
		cp $Jpreddata/ent/$Tid.ent $Tid.pdb
		$Modeller/modpy.sh python $Py
		mv *.B99990001.pdb $ID-$i.pdb
		rm $Tid.pdb *ali *py *ini *rsr *sch *V9999* *D00000*
	done
	echo "done">&3
	
	ls $ID-*.pdb | sed 's/.pdb//g' > $subpath/"$ID"_t.id
	
	#Add pdb head
	echo -n `date +%T`"	Add PDB head...">&3
	for id in `cat $subpath/"$ID"_t.id`
	do
		num=${id#*-}
		model=`sed -n $(($num*3-2))'p' $decoypath/$ID.aln | sed 's/^>//g' | sed 's/ .*$//g'`
		if [ "${#model}" -eq "5" ]; then
			Parent=`echo ${model:0:4}_${model:4:1}`
		else
			Parent=$model
		fi
		sed -i '1i PFRMAT TS\nTARGET '$ID'\nAUTHOR 3986-3318-3807\nMETHOD Jiang_Server:threading\nMODEL \nPARENT '$Parent $decoypath/$id.pdb
	done
	echo "done">&3
	
	maxZ=`sed -n '1p' $decoypath/$ID.aln | sed 's/^.* //g'`
	if [ `echo "$maxZ>$cutoff"|bc` -eq 1 ];then
		echo "Max(Z-score) > 6">&3
		# Get the best five decoys:
		echo -n `date +%T`"	Select the best five models...">&3
		for id in `cat $subpath/"$ID"_t.id`
		do
			num=${id#*-}
			model=`sed -n $(($num*3-2))'p' $decoypath/$ID.aln | sed 's/^>//g' | sed 's/ .*$//g'`
			Zscore=`sed -n $(($num*3-2))'p' $decoypath/$ID.aln | sed 's/^.* //g'`
			echo -e "$id\t$model\t$Zscore" >> $subpath/."$ID"_t.score
		done
		cat $subpath/."$ID"_t.score | sort -grk 3 > $subpath/"$ID"_t.scoreZ
		rm $subpath/."$ID"_t.score 
		Output $subpath/"$ID"_t.scoreZ $decoypath $finalpath
		echo "done">&3
		#echo `date +%T`"	Finish...$ID">&3
		return 0
	else
		echo "Max(Z-score) <= 6">&3
	fi
	Evaluate
	# Get the best five decoys:
	echo -n `date +%T`"	Select the best five models...">&3
	local line=1
	for id in `cat $subpath/"$ID"_t.id`
	do
		num=${id#*-}
		model=`sed -n $(($num*3-2))'p' $decoypath/$ID.aln | sed 's/^>//g' | sed 's/ .*$//g'`
		Zscore=`sed -n $(($num*3-2))'p' $decoypath/$ID.aln | sed 's/^.* //g'`
		Pscore=`sed -n $line'p' $subpath/"$ID"_t.out`
		Mscore=`echo "scale=4;($Zscore+$k*$Pscore)/1" | bc`
		echo -e "$id\t$model\t$Zscore\t$Pscore\t$Mscore" >> $subpath/."$ID"_t.score
		line=$[$line+1]
	done
	cat $subpath/."$ID"_t.score | sort -grk 5 > $subpath/"$ID"_t.scoreM
	rm $subpath/."$ID"_t.score
	Output $subpath/"$ID"_t.scoreM $decoypath $finalpath
	echo "done">&3
	#echo `date +%T`"	Finish...$ID">&3
	return 0
}


Evaluate(){
	
	#output dir
	NNoutpath=$subpath/NNcm
	SSoutpath=$subpath/sssa
	dssppath=$subpath/dssp
	HCoutpath=$subpath/HCdecoy
	
	mkdir $NNoutpath
	mkdir $SSoutpath
	mkdir $dssppath
	mkdir $HCoutpath

	local idfile=$subpath/"$ID"_t.id
	local outfile=$subpath/"$ID"_t.out

	echo -n `date +%T`"	NNcon: calculate the ss/sa & cm...">&3
	$NNconExecSh $decoypath/$ID.fasta $NNoutpath
	echo "done">&3

	mv $NNoutpath/$ID'_fasta.sssa' $SSoutpath/

	echo -n `date +%T`"	Generate dssp files from pdb...">&3
	for id in `cat $idfile`
	do
		cd $dssppath
		$DSSPexec $decoypath/$id.pdb $dssppath/$id.dssp
	done
	echo "done">&3

	echo -n `date +%T`"	Get Hydrophobic Core information...">&3
	cd $HCpath
	for id in `cat $idfile`
	do
		mkdir $HCoutpath/$id
		$HCexec $id $decoypath/ $dssppath/ $HCoutpath/$id/
		if [ -n "`ls -A $HCoutpath/$id/`" ]; then
			ls $HCoutpath/$id/*inf | sed 's/^.*\///' | sed 's/.inf$//' >> $HCoutpath/$id/$id.hl
		fi
	done
	echo "done">&3

	echo -n `date +%T`"	Static features for model selection...">&3
	cd $Fea37Pack/Program/Fea37/
	$Fea37exec $idfile $decoypath/$ID.fasta $NNoutpath/$ID"_fasta.cm8a" $NNoutpath/$ID"_fasta.cm12a" $SSoutpath/$ID.sssa $decoypath/ $dssppath/ $subpath/HCdecoy/ $subpath/$ID.scale
	echo "done">&3

	echo -n `date +%T`"	Get the score of decoy...">&3
	$svmscaleExec -r $Par $subpath/$ID.scale > $subpath/$ID.scale1
	$svmpredictExec $subpath/$ID.scale1 $Model $outfile
	echo "done">&3
}

Output(){
	local num=`cat $1 | wc -l`
	[ $onum -gt $num ] && onum=$num
	for i in `seq 1 $onum`
	do
	id=`cat $1 | sed -n $i'p' | cut -f 1`
	cat $2/$id.pdb | awk -v line=$i 'BEGIN{OFS=" "}\
			{
			if($1~/^MODEL/)
				print $1,line;
			else if($1~/^TER/)
				print $1;
			else if($1!~/^EXPDTA/ && $1!~/^REMARK/)
				print $0;
			}' > $3/$ID-$i.pdb
	done

	cat /dev/null > $3/../$ID.txt

	for mo in `ls $3`
	do
		cat $3/$mo >> $3/../$ID.txt
	done
}

Summary(){
	if [ $2 -eq 0 ]; then
		echo `date +%T`"	Finish...$1" | tee -a $timelog
		return 0
	else
		echo `date +%T`"	ERROR...$1" | tee -a $timelog
		return 1
	fi
}
########################### Main ##############################
if [ $fflag -eq 0 ]; then
	usage
elif [ $fflag -eq 1 ]; then
	[ -n "$ID_tmp" ] && in_ID=$ID_tmp
	cp $in_id $mainpath/Input/$in_ID.fasta
	Predic $in_ID
	Summary $in_ID $?
	exit $?
else
	cp $in_dir/*.fasta $mainpath/Input/
	cd $in_dir
	file_num=`ls *.fasta | wc -l`
	error_num=0
	for in_ID in `ls *.fasta | sed 's/.fasta$//g'`
	do
		Predic $in_ID
		Summary $in_ID $?
		error_num=$[$error_num+$?]
	done
	if [ $error_num -eq 0 ]; then
		exit 0
	elif [ $error_num -eq $file_num ]; then
		exit 2
	else
		exit 1
	fi
fi