Mercurial > repos > dereeper > pangenome_explorer
comparison COG/bac-genomics-scripts/prot_finder/prot_finder_pipe.sh @ 3:e42d30da7a74 draft
Uploaded
author | dereeper |
---|---|
date | Thu, 30 May 2024 11:52:25 +0000 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
2:97e4e3e818b6 | 3:e42d30da7a74 |
---|---|
1 #!/bin/bash | |
2 set -e | |
3 | |
4 ############# | |
5 # Functions # | |
6 ############# | |
7 | |
8 usage () { | |
9 cat 1>&2 << EOF # ${0##*/} parameter expansion substitution with variable '0' to get shell script filename without path | |
10 Usage: ${0##*/} [OPTION] -q query.faa -f (embl|gbk) > blast_hits.tsv | |
11 or: ${0##*/} [OPTION] -q query.faa -s subject.faa -d result_dir \\ | |
12 > result_dir/blast_hits.tsv | |
13 | |
14 Bash wrapper script to run a pipeline consisting of optional | |
15 'cds_extractor.pl' (with its options '-p -f'), BLASTP, 'prot_finder.pl', | |
16 and optional Clustal Omega. 'cds_extractor.pl' (only for shell script | |
17 option '-f') and 'prot_finder.pl' either have to be installed in the | |
18 global PATH or present in the current working directory. BLASTP is run | |
19 with disabled query filtering, locally optimal Smith-Waterman alignments, | |
20 and increasing the number of database sequences to show alignments | |
21 to 500 for BioPerl parsing (legacy: '-F F -s T -b 500', plus: '-seg | |
22 no -use_sw_tback -num_alignments 500'). | |
23 | |
24 The script ends with the STDERR message 'Pipeline finished!', if this | |
25 is not the case have a look at the log files in the result directory | |
26 for errors. | |
27 | |
28 Mandatory options: | |
29 -q <str> Path to query protein multi-FASTA file (*.faa) | |
30 with unique FASTA IDs | |
31 -f <str> File extension for files in the current working | |
32 directory to use for 'cds_extractor.pl' (e.g. | |
33 'embl' or 'gbk'); excludes shell script option '-s' | |
34 or | |
35 -s <str> Path to subject protein multi-FASTA file (*.faa) | |
36 already created with 'cds_extractor.pl' (and its | |
37 options '-p -f'), will not run 'cds_extractor.pl'; | |
38 excludes shell script option '-f' | |
39 | |
40 Optional options: | |
41 -h Print usage | |
42 -d <str> Path to result folder [default = results_i#_cq#] | |
43 -p (legacy|plus) BLASTP suite to use [default = plus] | |
44 -e <real> E-value for BLASTP [default = 1e-10] | |
45 -t <int> Number of threads to be used for BLASTP and | |
46 Clustal Omega [default = all processors on | |
47 system] | |
48 -i <int> Query identity cutoff for significant hits | |
49 [default = 70] | |
50 -c <int> Query coverage cutoff [default = 70] | |
51 -k <int> Subject coverage cutoff [default = 0] | |
52 -b Give only best hit (highest identity) for each | |
53 subject sequence | |
54 -a Multiple alignment of each multi-FASTA result | |
55 file with Clustal Omega | |
56 -o <str> Path to executable Clustal Omega binary if not | |
57 in global PATH; requires shell script option '-a' | |
58 -m Clean up all non-essential files | |
59 | |
60 Author: Andreas Leimbach <aleimba[at]gmx[dot]de> | |
61 EOF | |
62 } | |
63 | |
64 | |
65 ### Check external dependencies | |
66 check_commands () { | |
67 which "$1" > /dev/null || err "Required executable '$1' not found in global PATH, please install.$2" | |
68 } | |
69 | |
70 ### Check cutoff options input | |
71 check_cutoff_options () { | |
72 local message="Option '-$2' requires an integer number >= 0 or <= 100 as value, not '$1'!" | |
73 [[ $1 =~ ^[0-9]+$ ]] || err "$message" | |
74 (( $1 <= 100 )) || err "$message" # arithmetic expression (can only handle integer math, not float) | |
75 } | |
76 | |
77 | |
78 ### Error messages | |
79 err () { | |
80 echo -e "\n### Fatal error: $*" 1>&2 | |
81 exit 1 | |
82 } | |
83 | |
84 | |
85 ### Run status of script to STDERR instead of STDOUT | |
86 msg () { | |
87 echo -e "# $*" 1>&2 | |
88 } | |
89 | |
90 | |
91 ######## | |
92 # MAIN # | |
93 ######## | |
94 | |
95 shopt -s extglob # enable extended globs for bash | |
96 | |
97 Cmdline="$*" | |
98 | |
99 ### Getopts | |
100 Blastp_Suite="plus" | |
101 Evalue="1e-10" | |
102 Threads="$(nproc --all)" # get max number of processors on system | |
103 Ident_Cut=70 | |
104 Cov_Query_Cutoff=70 | |
105 Cov_Subject_Cutoff=0 | |
106 | |
107 while getopts ':q:f:s:d:p:e:t:i:c:k:bao:mh' opt; do # beginning ':' indicates silent mode, trailing ':' after each option requires value | |
108 case $opt in | |
109 q) Query_File=$OPTARG | |
110 [[ -r $Query_File ]] || err "Cannot read query file '$Query_File'!" | |
111 ;; | |
112 f) Subject_Ext=$OPTARG | |
113 [[ -n "$(find . -maxdepth 1 -name "*.${Subject_Ext}" -print -quit)" ]] || err "No files with the option '-f' specified file extension '$Subject_Ext' found in the current working directory!" | |
114 ;; | |
115 s) Subject_File=$OPTARG | |
116 [[ -r $Subject_File ]] || err "Cannot read subject file '$Subject_File'!" | |
117 ;; | |
118 d) Result_Dir=$OPTARG;; # checked below | |
119 p) Blastp_Suite=$OPTARG | |
120 [[ $Blastp_Suite = @(plus|legacy) ]] || err "Option '-p' only allows 'plus' for BLASTP+ or 'legacy' for legacy BLASTP as value, not '$Blastp_Suite'!" # extended glob (regex more expensive) | |
121 ;; | |
122 e) Evalue=$OPTARG | |
123 [[ $Evalue =~ ^([0-9][0-9]*|[0-9]+e-[0-9]+)$ ]] || err "Option '-e' requires a real number (either integer or scientific exponential notation) as value, not '$Evalue'!" | |
124 ;; | |
125 t) Threads=$OPTARG | |
126 [[ $Threads =~ ^[1-9][0-9]*$ ]] || err "Option '-t' requires an integer > 0 as value, not '$Threads'!" | |
127 ;; | |
128 i) Ident_Cut=$OPTARG | |
129 check_cutoff_options "$Ident_Cut" "i" | |
130 ;; | |
131 c) Cov_Query_Cutoff=$OPTARG | |
132 check_cutoff_options "$Cov_Query_Cutoff" "c" | |
133 ;; | |
134 k) Cov_Subject_Cutoff=$OPTARG | |
135 check_cutoff_options "$Cov_Subject_Cutoff" "k" | |
136 ;; | |
137 b) Opt_Best_Hit=1;; | |
138 a) Opt_Align=1;; | |
139 o) Clustal_Path=$OPTARG | |
140 [[ -x $Clustal_Path ]] || err "Option '-o' requires the path to an executable Clustal Omega binary as value, not '$Clustal_Path'!" | |
141 ;; | |
142 m) Opt_Clean_Up=1;; | |
143 h) usage; exit;; # usage function, exit code zero | |
144 \?) err "Invalid option '-$OPTARG'. See usage with '-h'!";; | |
145 :) err "Option '-$OPTARG' requires a value. See usage with '-h'!";; | |
146 esac | |
147 done | |
148 | |
149 | |
150 ### Check options and enforce mandatory options | |
151 [[ $Query_File && ($Subject_Ext || $Subject_File) ]] || err "Mandatory options '-q' and '-f' or '-s' are missing!" | |
152 | |
153 [[ $Subject_Ext && $Subject_File ]] && err "Options '-f' and '-s' given which exclude themselves. Choose either '-f' OR '-s'!" | |
154 | |
155 (( Threads <= $(nproc) )) || err "Number of threads for option '-t', '$Threads', exceeds the maximum $(nproc) processors on the system!" | |
156 | |
157 [[ ! $Opt_Align && $Clustal_Path ]] && Opt_Align=1 && msg "Option '-o' requires option '-a', forcing option '-a'!" | |
158 | |
159 | |
160 ### Check external dependencies | |
161 echo 1>&2 # newline | |
162 msg "Checking pipeline dependencies" | |
163 [[ $Opt_Align && ! $Clustal_Path ]] && check_commands "clustalo" " Or use option '-o' to give the path to the binary!" | |
164 | |
165 for exe in cds_extractor.pl formatdb blastall makeblastdb blastp prot_finder.pl; do | |
166 [[ $Subject_File && $exe == cds_extractor.pl ]] && continue | |
167 [[ $Blastp_Suite == legacy && $exe = @(makeblastdb|blastp) ]] && continue # extended glob | |
168 [[ $Blastp_Suite == plus && $exe = @(formatdb|blastall) ]] && continue | |
169 if [[ $exe = *.pl ]]; then # glob | |
170 if [[ -r "./$exe" ]]; then # present in current wd | |
171 [[ $exe =~ ^cds ]] && Cds_Extractor_Cmd="perl cds_extractor.pl" | |
172 [[ $exe =~ ^prot ]] && Prot_Finder_Cmd="perl prot_finder.pl" | |
173 continue | |
174 else | |
175 [[ $exe =~ ^cds ]] && Cds_Extractor_Cmd="cds_extractor.pl" | |
176 [[ $exe =~ ^prot ]] && Prot_Finder_Cmd="prot_finder.pl" | |
177 check_commands "$exe" " Or copy the Perl script in the current working directory." | |
178 fi | |
179 continue | |
180 fi | |
181 check_commands "$exe" | |
182 done | |
183 | |
184 msg "Script call command: ${0##*/} $Cmdline" | |
185 | |
186 | |
187 ### Create result folder | |
188 if [[ ! $Result_Dir ]]; then # can't give default before 'getopts' in case cutoffs are set by the user | |
189 Result_Dir="results_i${Ident_Cut}_cq${Cov_Query_Cutoff}" | |
190 else | |
191 Result_Dir="${Result_Dir%/}" # parameter expansion substitution to get rid of a potential '/' at the end of Result_Dir path | |
192 fi | |
193 | |
194 if [[ -d $Result_Dir ]]; then # make possible to redirect STDOUT output into result_dir (corresponding to option '-f' in 'protein_finder.pl' script) | |
195 skip=0 | |
196 for file in "$Result_Dir"/*; do | |
197 if [[ -s $file || $skip -eq 1 ]]; then # die if a file with size > 0 or more than one file already in result_dir | |
198 err "Result directory '$Result_Dir' already exists! You can use option '-d' to set a different result directory name." | |
199 fi | |
200 skip=1 | |
201 done | |
202 else | |
203 mkdir -pv "$Result_Dir" 1>&2 | |
204 fi | |
205 | |
206 | |
207 ### Run cds_extractor.pl | |
208 if [[ $Subject_Ext ]]; then | |
209 msg "Running cds_extractor.pl on all '*.$Subject_Ext' files in the current working directory" | |
210 for file in *."$Subject_Ext"; do | |
211 file_no_ext="${file%.${Subject_Ext}}.faa" # parameter expansion substitution to get rid of file extension and replace with new one (*.faa are the output files from cds_extractor) | |
212 File_Names+=("$file_no_ext") # append to array | |
213 eval "$Cds_Extractor_Cmd -i $file -p -f &>> $Result_Dir/cds_extractor.log" # '&>' instead of '/dev/null' for error catching | |
214 done | |
215 Subject_File="$Result_Dir/prot_finder.faa" # for creating BLASTP db below | |
216 cat "${File_Names[@]}" > "$Subject_File" # concatenate files stored in the array, "${array[@]}" expands to list of array elements (words) | |
217 fi | |
218 | |
219 | |
220 ### Run BLASTP | |
221 msg "Running BLASTP '$Blastp_Suite' with subject '$Subject_File', query '$Query_File', evalue '$Evalue', and $Threads threads" | |
222 Blast_Report="$Result_Dir/prot_finder.blastp" | |
223 if [[ $Blastp_Suite == legacy ]]; then | |
224 formatdb -p T -i "$Subject_File" -n prot_finder_db | |
225 blastall -p blastp -d prot_finder_db -i "$Query_File" -o "$Blast_Report" -e "$Evalue" -F F -s T -b 500 -a "$Threads" | |
226 elif [[ $Blastp_Suite == plus ]]; then | |
227 makeblastdb -in "$Subject_File" -input_type fasta -dbtype prot -out prot_finder_db &> "$Result_Dir/makeblastdb.log" # '&>' instead of '/dev/null' for error catching | |
228 blastp -db prot_finder_db -query "$Query_File" -out "$Blast_Report" -evalue "$Evalue" -seg no -use_sw_tback -num_alignments 500 -num_threads "$Threads" | |
229 fi | |
230 | |
231 | |
232 ### Run prot_finder.pl | |
233 msg "Running prot_finder.pl with identity cutoff '$Ident_Cut', query coverage cutoff '$Cov_Query_Cutoff', and subject coverage cutoff '$Cov_Subject_Cutoff'" | |
234 Cmd="$Prot_Finder_Cmd -d $Result_Dir -f -q $Query_File -s $Subject_File -r $Blast_Report -i $Ident_Cut -cov_q $Cov_Query_Cutoff -cov_s $Cov_Subject_Cutoff" | |
235 [[ $Opt_Best_Hit ]] && Cmd="$Cmd -b" # append to command | |
236 [[ $Opt_Align ]] && Cmd="$Cmd -a -t $Threads" | |
237 [[ $Clustal_Path ]] && Cmd="$Cmd -p $Clustal_Path" | |
238 eval "$Cmd" 2> "$Result_Dir/prot_finder.log" # '2>' instead of '/dev/null' for error catching | |
239 | |
240 msg "All result files stored in directory '$Result_Dir'" | |
241 | |
242 | |
243 ### Clean up non-essential files | |
244 if [[ $Opt_Clean_Up ]]; then | |
245 msg "Removing non-essential output files, option '-m'" | |
246 for file in "${File_Names[@]}"; do # remove output files from cds_extractor | |
247 rm -v "$file" 1>&2 | |
248 done | |
249 [[ $Subject_Ext ]] && rm -v "$Subject_File" 1>&2 # 'cat' from cds_extractor | |
250 if [[ $Blastp_Suite == legacy ]]; then | |
251 rm -v formatdb.log 1>&2 | |
252 [[ -r error.log ]] && rm -v error.log 1>&2 # no idea where this guy is coming from or what is its trigger | |
253 fi | |
254 rm -v prot_finder_db.p* "$Blast_Report" "$Result_Dir"/*.log "${Subject_File}.idx" 1>&2 | |
255 fi | |
256 | |
257 msg "Pipeline finished!" |