-
Notifications
You must be signed in to change notification settings - Fork 2
/
StrepLab-JanOw_GBS-Typer.sh
executable file
·198 lines (172 loc) · 6.79 KB
/
StrepLab-JanOw_GBS-Typer.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
#!/bin/bash -l
temp_path=$(pwd)
export PATH=$PATH:$temp_path
## -- begin embedded SGE options --
read -a PARAM <<< $(/bin/sed -n ${SGE_TASK_ID}p $1/job-control.txt)
## -- end embedded SGE options --
###Load Modules###
#. /usr/share/Modules/init/bash
module load perl/5.22.1
module load ncbi-blast+/2.2.29
module load BEDTools/2.17.0
module load freebayes/0.9.21
module load prodigal/2.60
module load cutadapt/1.8.3
module load srst2/0.1.7
###This script is called for each job in the qsub array. The purpose of this code is to read in and parse a line of the job-control.txt file
###created by 'StrepLab-JanOw_GAS-wrapr.sh' and pass that information, as arguments, to other programs responsible for various parts of strain
###characterization (MLST, emm type and antibiotic drug resistance prediction).
readPair_1=${PARAM[0]}
readPair_2=${PARAM[1]}
allDB_dir=${PARAM[2]}
batch_out=${PARAM[3]}
sampl_out=${PARAM[4]}
delete='true'
###Start Doing Stuff###
cd "$sampl_out"
batch_name=$(echo "$readPair_1" | awk -F"/" '{print $(NF-4)}')
out_name=$(echo "$readPair_1" | awk -F"/" '{print $(NF-4)"--"$(NF)}' | sed 's/_S[0-9]\+_L[0-9]\+_R[0-9]\+.*//g') ###Use This For Batches off the MiSeq###
#out_name=$(echo "$readPair_1" | awk -F"/" '{print $(NF-1)"--"$(NF)}' | sed 's/_S[0-9]\+_L[0-9]\+_R[0-9]\+.*//g') ###Otherwise Use This###
just_name=$(echo "$readPair_1" | awk -F"/" '{print $(NF)}' | sed 's/_S[0-9]\+_L[0-9]\+_R[0-9]\+.*//g')
out_nameMLST=MLST_"$just_name"
###Pre-Process Paired-end Reads###
fastq1_trimd=cutadapt_"$just_name"_S1_L001_R1_001.fastq
fastq2_trimd=cutadapt_"$just_name"_S1_L001_R2_001.fastq
cutadapt -b AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC -q 20 --minimum-length 50 --paired-output temp2.fastq -o temp1.fastq $readPair_1 $readPair_2
cutadapt -b AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT -q 20 --minimum-length 50 --paired-output $fastq1_trimd -o $fastq2_trimd temp2.fastq temp1.fastq
rm temp1.fastq
rm temp2.fastq
module load fastqc/0.11.5
mkdir "$just_name"_R1_cut
mkdir "$just_name"_R2_cut
fastqc "$fastq1_trimd" --outdir=./"$just_name"_R1_cut
fastqc "$fastq2_trimd" --outdir=./"$just_name"_R2_cut
module unload fastqc/0.11.5
###Call MLST###
srst2 --samtools_args "\-A" --mlst_delimiter '_' --input_pe "$readPair_1" "$readPair_2" --output "$out_nameMLST" --save_scores --mlst_db "$allDB_dir/Streptococcus_agalactiae.fasta" --mlst_definitions "$allDB_dir/sagalactiae.txt" --min_coverage 99.999
###Check and extract new MLST alleles###
MLST_allele_checkr.pl "$out_nameMLST"__mlst__Streptococcus_agalactiae__results.txt "$out_nameMLST"__*.Streptococcus_agalactiae.sorted.bam "$allDB_dir/Streptococcus_agalactiae.fasta"
###Call GBS Serotype###
GBS_Serotyper.pl -1 "$readPair_1" -2 "$readPair_2" -r "$allDB_dir/GBS_seroT_Gene-DB_Final.fasta" -n "$just_name"
###Call GBS bLactam Resistances###
module unload perl/5.22.1
module load perl/5.16.1-MT
PBP-Gene_Typer.pl -1 "$readPair_1" -2 "$readPair_2" -r "$allDB_dir/GBS_bLactam_Ref.fasta" -n "$just_name" -s GBS -p 1A,2B,2X
module unload perl/5.16.1-MT
module load perl/5.22.1
###Call GBS Misc. Resistances###
GBS_Res_Typer.pl -1 "$readPair_1" -2 "$readPair_2" -d "$allDB_dir" -r GBS_Res_Gene-DB_Final.fasta -n "$just_name"
GBS_Target2MIC.pl TEMP_Res_Results.txt "$just_name" TEMP_pbpID_Results.txt
###Type Surface Proteins###
GBS_Surface_Typer.pl -1 "$readPair_1" -2 "$readPair_2" -r "$allDB_dir" -p GBS_Surface_Gene-DB_Final.fasta -n "$just_name"
###Output the emm type/MLST/drug resistance data for this sample to it's results output file###
tabl_out="TABLE_Isolate_Typing_results.txt"
bin_out="BIN_Isolate_Typing_results.txt"
printf "$just_name\t" >> "$tabl_out"
printf "$just_name," >> "$bin_out"
###Serotype Output###
sero_out="NF"
while read -r line
do
if [[ -n "$line" ]]
then
justTarget=$(echo "$line" | awk -F"\t" '{print $3}')
if [[ "$sero_out" == "NF" ]]
then
sero_out="$justTarget"
else
sero_out="$sero_out;$justTarget"
fi
fi
done <<< "$(sed 1d TEMP_SeroType_Results.txt)"
printf "$sero_out\t" >> "$tabl_out"
printf "$sero_out," >> "$bin_out"
###MLST OUTPUT###
sed 1d "$out_nameMLST"__mlst__Streptococcus_agalactiae__results.txt | while read -r line
do
MLST_tabl=$(echo "$line" | cut -f2-9)
echo "MLST line: $MLST_tabl\n";
printf "$MLST_tabl\t" >> "$tabl_out"
MLST_val=$(echo "$line" | awk -F" " '{print $2}')
printf "$MLST_val," >> "$bin_out"
done
###PBP_ID Output###
justPBPs="NF"
sed 1d TEMP_pbpID_Results.txt | while read -r line
do
if [[ -n "$line" ]]
then
justPBPs=$(echo "$line" | awk -F"\t" '{print $2}' | tr ':' '\t')
justPBP_BIN=$(echo "$line" | awk -F"\t" '{print $2}' | tr ':' ',')
fi
printf "$justPBPs\t" >> "$tabl_out"
printf "$justPBP_BIN," >> "$bin_out"
done
###Resistance Targets###
#while read -r line
#do
# RES_targ=$(echo "$line" | cut -f2)
# printf "$RES_targ\t" >> "$tabl_out"
#done < TEMP_Res_Results.txt
#cat BIN_Res_Results.txt | sed 's/$/,/g' >> "$bin_out"
###Resistance Targets###
while read -r line
do
#RES_targ=$(echo "$line" | cut -f2)
#printf "$RES_targ\t" >> "$tabl_out"
printf "$line\t" | tr ',' '\t' >> "$tabl_out"
done < RES-MIC_"$just_name"
###Surface Targets###
while read -r line
do
SURF_targ=$(echo "$line" | cut -f2)
printf "$SURF_targ\t" >> "$tabl_out"
done < TEMP_Surface_Results.txt
cat BIN_Surface_Results.txt >> "$bin_out"
if [[ -e $(echo ./velvet_output/*_Logfile.txt) ]]
then
vel_metrics=$(echo ./velvet_output/*_Logfile.txt)
print "velvet metrics file: $vel_metrics\n";
velvetMetrics.pl -i "$vel_metrics";
line=$(cat velvet_qual_metrics.txt | tr ',' '\t')
printf "$line\t" >> "$tabl_out"
printf "$readPair_1\t" >> "$tabl_out";
pwd | xargs -I{} echo {}"/velvet_output/contigs.fa" >> "$tabl_out"
else
printf "NA\tNA\tNA\tNA\t$readPair_1\tNA\n" >> "$tabl_out"
fi
#printf "\n" >> "$tabl_out"
#printf "\n" >> "$tabl_out"
#printf "\n" >> "$bin_out"
###Remove Unneeded Files###
if ${delete};
then
echo "The delete flag is true. Will delete all unneeded files"
rm *.pileup
rm *.*sorted.*am
rm *.scores
rm *ARGannot_r1.*
rm *__genes__*
rm *log
rm *Results.txt
rm ARG-RESFI_fullgenes_results.txt BIN_Isolate_Typing_results.txt BLACTAM_MIC_RF_with_SIR.txt Serotype_Extraction_Sequence.fna velvet_qual_metrics.txt ./velvet_output/Sequences ./velvet_output/stats.txt ./velvet_output/Log BLACTAM_MIC_RF.txt
rm RES-MIC*
rm cutadapt_*_R*
rm FOLP_target_*
rm OUT_*
rm TEMP_*
rm velvet_output/*Graph*
rm -r PBP_to_MIC_temp/
rm -r *_R1_cut
rm -r *_R2_cut
else
echo "The delete flag is false. Will keep all unneeded files"
fi
###Unload Modules###
module unload perl/5.22.1
module unload ncbi-blast+/2.2.29
module unload BEDTools/2.17.0
module unload freebayes/0.9.21
module unload prodigal/2.60
module unload cutadapt/1.8.3
module unload srst2/0.1.7