org-annotate.sh 3.27 KB
Newer Older
1 2 3 4
#!/bin/bash
#
#
#
alain viari's avatar
alain viari committed
5
#                           Annotate Organelle 
6 7 8 9 10 11
#
#========================================================================================
#
#
#========================================================================================

alain viari's avatar
alain viari committed
12
# -- CAUTION -- Works as long as the script 
13
#               is not called through a symlink
alain viari's avatar
alain viari committed
14 15
THIS_DIR="$(dirname ${BASH_SOURCE[0]})"
source "${THIS_DIR}/scripts/bash_init.sh"
16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

pushTmpDir ORG.organnot

	if [[ ! "$1" =~ ^/ ]]; then
		QUERY="${CALL_DIR}/$1"
	else
		QUERY="$1"
	fi

	RESULTS=$(basename ${QUERY/.*/})
	LOG="${CALL_DIR}/${RESULTS}.log"
	
	rm -f ${LOG}
	openLogFile ${LOG}

31 32 33 34
	loginfo "Normalizing the structure of the Chloroplast sequence..."
		loginfo "   LSC + IRB + SSC + IRA"
		${PROG_DIR}/detectors/normalize/bin/go_normalize.sh ${QUERY} > "${RESULTS}.norm.fasta"
	loginfo "Done."
35
	
36
	loginfo "Annotating the Inverted repeats and Single copies (LSC and SSC)..."
37
		${PROG_DIR}/detectors/ir/bin/go_ir.sh "${RESULTS}.norm.fasta" > "${RESULTS}.annot"		
38
	loginfo "Done."
39
	
40
	loginfo "Annotating the tRNA..."
41
		${PROG_DIR}/detectors/trna/bin/go_trna.sh "${RESULTS}.norm.fasta" >> "${RESULTS}.annot"
42
	loginfo "Done."
43
	
44
	loginfo "Annotating the rRNA genes..."
45
		${PROG_DIR}/detectors/rrna/bin/go_rrna.sh "${RESULTS}.norm.fasta" >> "${RESULTS}.annot"
46
	loginfo "Done."
alain viari's avatar
alain viari committed
47 48

	loginfo "Annotating the CDS..."
49
		${PROG_DIR}/detectors/cds/bin/go_cds.sh "${RESULTS}.norm.fasta" >> "${RESULTS}.annot"
alain viari's avatar
alain viari committed
50
	loginfo "Done."
51 52 53 54 55 56 57
	
	loginfo "Printing annotations header..."
		echo "XX"
    	echo "FH   Key             Location/Qualifiers"
	loginfo "Done."
	
	loginfo "Ordering annotations..."
58
		awk '/^.....(misc|repeat|rRNA|tRNA|gene)/ { \
59 60
		        match($3,"[0-9][0-9]*"); \
		        pos=substr($3,RSTART,RLENGTH)*1000 + 1; \
61 62 63
		        print pos,$0;    \
		        next} \
		      { pos++; \
64
		        print pos,$0}' "${RESULTS}.annot" | \
65
		sort -nk1 |\
66 67
		awk '{ \
		        match($0,"^[0-9]* ");\
68
		        line=substr($0,RLENGTH+1);\
69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98
		        print line}' 
	loginfo "Done."
	
	loginfo "Closing annotations table..."
		echo "XX"
	loginfo "Done."
	
	loginfo "Computing statistics on nucleotide usage..."
		awk '! /^>/ { \
			    seq=toupper($0); \
				gsub(" ","",seq); \
			    lseq=length(seq); \
				for (i=0; i < lseq; i++) { \
					freq[substr(seq,i,1)]++}\
					} \
			 END { \
			 	other=0; \
			 	for (i in freq) { \
			 		if (i!="A" && i!="C" && i!="G" && i!="T") {\
			 			other+=freq[i] \
			 			} \
			 		}; \
			 		print "SQ   Sequence "\
			 		      (freq["A"]+freq["C"]+freq["G"]+freq["T"]+other) \
			 		      " BP; "\
			 		      freq["A"]" A; "\
			 		      freq["C"]" C; "\
			 		      freq["G"]" G; "\
			 		      freq["T"]" T; "\
			 		      other" other;" \
99
			 }' "${RESULTS}.norm.fasta"
100 101 102
	loginfo "Done."
	
	loginfo "Reformating sequences..."
103
		lines=$(wc -l "${RESULTS}.norm.fasta" | awk '{print $1}')
104 105 106 107 108 109 110 111 112 113 114 115 116 117
		awk -v lines=$lines ' \
			! /^>/ { \
					seq=tolower($0); \
					gsub(" ","",seq); \
					printf("     ") ;\
					for (i=0; i < 6; i++) { \
						f=substr(seq,i * 10, 10); \
						pos+=length(f); \
						f = f  substr("          ",1,10-length(f)); \
						printf("%s ",f) \
					}; \
					if (NR==lines) \
					  {pos-=1}; \
					printf("   %6d\n",pos) \
118
			   }' "${RESULTS}.norm.fasta"
119 120 121 122 123 124
	loginfo "Done."
	
	loginfo "Closing sequence part..."
		echo "//"
	loginfo "Done."

125 126
popTmpDir