forked from theaidenlab/AGWG-merge
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathrun-asm-pipeline-post-review.sh
225 lines (185 loc) · 9.16 KB
/
run-asm-pipeline-post-review.sh
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
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
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
#!/bin/bash
##########
#The MIT License (MIT)
#
# Copyright (c) 2018 Aiden Lab
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
##########
#
# 3D DNA de novo genome assembly pipeline: 180114 version.
#
echo `readlink -f $0`" "$*
USAGE_short="
*****************************************************
3D de novo assembly: version 180114
USAGE: ./run-asm-pipeline-post-review.sh [options] -r <review.assembly> <path_to_input_fasta> <path_to_input_mnd>
DESCRIPTION:
This is a script to finalize assemblies (represented in input by draft fasta and deduplicated list of alignments of Hi-C reads to this fasta as produced by the Juicer pipeline) into chromosome-length fasta sequences, after review in Juicebox Assembly Tools Module (represented by review.assembly file). The script will produce an output fasta file, a Hi-C map of the final assembly, and a few supplementary annotation files to help review the result in Juicebox.
ARGUMENTS:
path_to_input_fasta Specify file path to draft assembly fasta file.
path_to_input_mnd Specify path to deduplicated list of alignments of Hi-C reads to the draft assembly fasta as produced by the Juicer pipeline: the merged_nodups file (mnd).
OPTIONS:
-r|--review Path to review \".assembly\" file.
--sort-output Sorts output scaffolds by size (prior to adding gaps), default is false.
-i|--input input_size Specifies threshold input contig/scaffold size (default is 15000). Contigs/scaffolds smaller than input_size are going to be ignored. Only matters if running including seal. Should be the same as used for running the original script.
-s|--stage stage Assembly steps to run on top of the reviewed assembly, can be seal and finalize. Default is finalize.
-q|--mapq mapq Mapq threshold for final map visualization, default is 1.
-h|--help Shows this help. Type --help for a full set of options.
*****************************************************
"
pipeline=`cd "$( dirname $0)" && pwd`
## default parameter setup
input_size=15000 # contigs/scaffolds smaller than input_size are ignored
mapq=1 # read mapping quality threshold for Hi-C scaffolder
stage="finalize" # by default run only final pipeline
sort_output=false
############### HANDLE OPTIONS ###############
while :; do
case $1 in
-h|--help)
echo "$USAGE_short" >&1
exit 0
;;
-r|--review) OPTARG=$2
if [[ -f $OPTARG ]]; then
echo " -r|--review flag was triggered, treating file $OPTARG as a JB4A review file for draft fasta in arguments." >&1
review_assembly=$OPTARG
else
echo ":( File not found in the suggested review assembly file path. Exiting!" >&2
exit 1
fi
shift
;;
-i|--input) OPTARG=$2
re='^[0-9]+$'
if [[ $OPTARG =~ $re ]]; then
echo " -i|--input flag was triggered, filtering draft contigs/scaffolds smaller than $OPTARG." >&1
input_size=$OPTARG
else
echo ":( Wrong syntax for input size threshold. Using the default value ${input_size}." >&2
fi
shift
;;
-s|--stage) OPTARG=$2
if [ "$OPTARG" == "seal" ] || [ "$OPTARG" == "finalize" ]; then
echo " -s|--stage flag was triggered, fast-forwarding to \"$OPTARG\" pipeline section." >&1
stage=$OPTARG
else
echo ":( Wrong syntax for pipeline stage. Exiting!" >&2
fi
shift
;;
-q|--mapq) OPTARG=$2 ##TODO: check that propagates consistently, not tested sufficiently
re='^[0-9]+$'
if [[ $OPTARG =~ $re ]]; then
echo " -q|--mapq flag was triggered, scaffolding using reads with at least $OPTARG mapping quality." >&1
mapq=$OPTARG
else
echo ":( Wrong syntax for mapping quality. Using the default value ${mapq}." >&2
fi
shift
;;
## long menu options
--sort-output)
echo " --sort-output was triggered, will sort output scaffolds by size." >&1
sort_output=true
;;
--) # End of all options
shift
break
;;
-?*)
echo ":| WARNING: Unknown option. Ignoring: ${1}" >&2
;;
*) # Default case: If no more options then break out of the loop.
break
esac
shift
done
############### HANDLE EXTERNAL DEPENDENCIES ###############
## GNU Parallel Dependency
parallel="false"
if hash parallel 2>/dev/null; then
ver=`parallel --version | awk 'NR==1{print \$3}'`
[ $ver -ge 20150322 ] && parallel="true"
fi
[ $parallel == "false" ] && echo ":| WARNING: GNU Parallel version 20150322 or later not installed. We highly recommend to install it to increase performance. Starting pipeline without parallelization!" >&2
############### HANDLE ARGUMENTS ###############
[ -z ${review_assembly} ] || [ -z $1 ] || [ -z $2 ] && echo >&2 "Not sure how to parse your input: files not listed or not found at expected locations. Exiting!" && echo >&2 "$USAGE_short" && exit 1
[ ! -s ${review_assembly} ] || [ ! -s $1 ] || [ ! -s $2 ] && echo >&2 "Not sure how to parse your input: files not listed or not found at expected locations. Exiting!" && echo >&2 "$USAGE_short" && exit 1
## TODO: check file format
if [ "$#" -ne 2 ]; then
echo >&2 "Illegal number of arguments. Please double check your input. Exiting!" && echo >&2 "$USAGE_short" && exit 1
fi
orig_fasta=$1
orig_mnd=$2
genomeid=$(basename "$orig_fasta" .fasta)
genomeid=$(basename "$genomeid" .fna)
genomeid=$(basename "$genomeid" .fa)
if [ ${orig_mnd} != ${genomeid}.mnd.txt ]; then
if [ -f ${genomeid}.mnd.txt ]; then
cmp --silent ${orig_mnd} ${genomeid}.mnd.txt || { echo >&2 ":( Please remove or rename file ${genomeid}.mnd.txt. Exiting!" && exit 1; }
fi
ln -sf ${orig_mnd} ${genomeid}.mnd.txt
fi
orig_mnd=${genomeid}.mnd.txt
## TODO: add some checks to make more user-proof against swapping
[ -f ${genomeid}.cprops ] || awk -f ${pipeline}/utils/generate-sorted-cprops-file.awk ${orig_fasta} > ${genomeid}.cprops
orig_cprops=${genomeid}.cprops
[ ! -f ${orig_cprops} ] && echo >&2 ":( No cprops file found. Please rerun the pipeline from scratch. Exiting!" && exit 1
## calculate zoom
# TODO: move this to mismatch detector, pass only scale
totlength=`awk '{total+=$3}END{print total}' ${orig_cprops}`
scale=$(( 1 + $totlength / 2100000000 ))
if [ $scale -ne 1 ]; then
editor_coarse_resolution=$((editor_coarse_resolution/scale))
editor_coarse_region=$((editor_coarse_region/scale))
editor_fine_resolution=$((editor_fine_resolution/scale))
polisher_coarse_resolution=$((polisher_coarse_resolution/scale))
polisher_coarse_region=$((polisher_coarse_region/scale))
polisher_fine_resolution=$((polisher_fine_resolution/scale))
splitter_coarse_resolution=$((splitter_coarse_resolution/scale))
splitter_coarse_region=$((splitter_coarse_region/scale))
splitter_fine_resolution=$((splitter_fine_resolution/scale))
fi
if [ "$stage" != "finalize" ]; then
############### SEALING ###############
awk -v cprops=${genomeid}.split.cprops -v asm=${genomeid}.split.asm '$1~/^>/{$1=substr($1,2); print > cprops;next}{print > asm}' ${review_assembly}
if [ "$sort_output" == "true" ]; then
bash ${pipeline}/run-asm-pipeline.sh -s seal -i ${input_size} --sort-output ${orig_fasta} ${orig_mnd}
else
bash ${pipeline}/run-asm-pipeline.sh -s seal -i ${input_size} ${orig_fasta} ${orig_mnd}
fi
else
############### FINALIZING ###############
echo "###############" >&1
echo "Finilizing output:" >&1
awk -v cprops=${genomeid}.final.cprops -v asm=${genomeid}.final.asm '$1~/^>/{$1=substr($1,2); print > cprops;next}{print > asm}' ${review_assembly}
if [ "$sort_output" == "true" ]; then
awk -v input_size=${input_size} -f ${pipeline}/utils/sort-asm-by-size.awk ${genomeid}.final.cprops ${genomeid}.final.asm > ${genomeid}.final.asm.tmp && mv ${genomeid}.final.asm.tmp ${genomeid}.final.asm
fi
# build final map
bash ${pipeline}/edit/edit-mnd-according-to-new-cprops.sh ${genomeid}.final.cprops ${orig_mnd} > ${genomeid}.final.mnd.txt
bash ${pipeline}/visualize/run-asm-visualizer.sh -p ${parallel} -q ${mapq} -i ${genomeid}.final.cprops ${genomeid}.final.asm ${genomeid}.final.mnd.txt
rm ${genomeid}.final.mnd.txt temp.${genomeid}.final.asm_mnd.txt
# build final fasta
awk -f ${pipeline}/edit/edit-fasta-according-to-new-cprops.awk ${genomeid}.final.cprops ${orig_fasta} > ${genomeid}.final.fasta
bash ${pipeline}/finalize/finalize-output.sh -s ${input_size} -l ${genomeid} ${genomeid}.final.cprops ${genomeid}.final.asm ${genomeid}.final.fasta final
fi