forked from atarashansky/SAMap
-
Notifications
You must be signed in to change notification settings - Fork 0
/
map_genes.sh
executable file
·145 lines (130 loc) · 4.63 KB
/
map_genes.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
#!/bin/bash
#parse arguments
die() {
printf '%s\n' "$1" >&2
exit 1
}
_help() {
printf "Usage: map_genes.sh [--tr1] [--t1] [--n1] [--tr2] [--t2] [--n2]\n\t[--tr1]: path to transcriptome/proteome 1\n\t[--t1]: is 1 a transcriptome [nucl] or proteome [prot]\n\t[--n1]: two character identifier of 1\n\t[--tr2]: path to transcriptome/proteome 2\n\t[--t2]: is 2 a transcriptome [nucl] or proteome [prot]\n\t[--n2]: two character identifier of 2"
}
if [ $# -eq 0 ]; then
_help
exit 1
fi
while :; do
case $1 in
-h|-\?|--help)
_help
exit
;;
--tr1) # Takes an option argument; ensure it has been specified.
if [ "$2" ]; then
tr1=$2
shift
else
die 'ERROR: "--tr1" requires a non-empty option argument.'
fi
;;
--tr2) # Takes an option argument; ensure it has been specified.
if [ "$2" ]; then
tr2=$2
shift
else
die 'ERROR: "--tr2" requires a non-empty option argument.'
fi
;;
--n1) # Takes an option argument; ensure it has been specified.
if [ "$2" ]; then
n1=$2
shift
else
die 'ERROR: "--n1" requires a non-empty option argument.'
fi
;;
--n2) # Takes an option argument; ensure it has been specified.
if [ "$2" ]; then
n2=$2
shift
else
die 'ERROR: "--n2" requires a non-empty option argument.'
fi
;;
--t1) # Takes an option argument; ensure it has been specified.
if [ "$2" ]; then
t1=$2
shift
else
die 'ERROR: "--t1" requires a non-empty option argument.'
fi
;;
--t2) # Takes an option argument; ensure it has been specified.
if [ "$2" ]; then
t2=$2
shift
else
die 'ERROR: "--t2" requires a non-empty option argument.'
fi
;;
--threads) # Takes an option argument; ensure it has been specified.
if [ "$2" ]; then
n_threads=$2
shift
fi
;;
--) # End of all options.
shift
break
;;
-?*)
printf 'WARN: Unknown option (ignored): %s\n' "$1" >&2
;;
*) # Default case: No more options, so break out of the loop.
break
esac
shift
done
if [ -z ${n_threads+x} ]; then
n_threads=8
fi
#read -e -p 'Transcriptome 1: ' tr1
#read -e -p 'Transcriptome type: ' t1
#read -e -p 'Transcriptome ID: ' n1
#read -e -p 'Transcriptome 2: ' tr2
#read -e -p 'Transcriptome type: ' t2
#read -e -p 'Transcriptome ID: ' n2
n1="${n1:0:2}"
n2="${n2:0:2}"
mkdir -p "maps"
mkdir -p "maps/${n1}${n2}"
if [ ! -f "${tr1}.nhr" ] && [ ! -f "${tr1}.phr" ]
then
makeblastdb -in $tr1 -dbtype $t1
fi
if [ ! -f "${tr2}.nhr" ] && [ ! -f "${tr2}.phr" ]
then
makeblastdb -in $tr2 -dbtype $t2
fi
if [[ "$t1" == "nucl" && "$t2" == "nucl" ]]
then
echo "Running tblastx in both directions"
tblastx -query $tr1 -db $tr2 -outfmt 6 -out "maps/${n1}${n2}/${n1}_to_${n2}.txt" -num_threads ${n_threads} -max_hsps 1 -evalue 1e-6
tblastx -query $tr2 -db $tr1 -outfmt 6 -out "maps/${n1}${n2}/${n2}_to_${n1}.txt" -num_threads ${n_threads} -max_hsps 1 -evalue 1e-6
fi
if [[ "$t1" == "nucl" && "$t2" == "prot" ]]
then
echo "Running blastx from 1 to 2 and tblastn from 2 to 1"
blastx -query $tr1 -db $tr2 -outfmt 6 -out "maps/${n1}${n2}/${n1}_to_${n2}.txt" -num_threads ${n_threads} -max_hsps 1 -evalue 1e-6
tblastn -query $tr2 -db $tr1 -outfmt 6 -out "maps/${n1}${n2}/${n2}_to_${n1}.txt" -num_threads ${n_threads} -max_hsps 1 -evalue 1e-6
fi
if [[ "$t1" == "prot" && "$t2" == "nucl" ]]
then
echo "Running tblastn from 1 to 2 and blastx from 2 to 1"
tblastn -query $tr1 -db $tr2 -outfmt 6 -out "maps/${n1}${n2}/${n1}_to_${n2}.txt" -num_threads ${n_threads} -max_hsps 1 -evalue 1e-6
blastx -query $tr2 -db $tr1 -outfmt 6 -out "maps/${n1}${n2}/${n2}_to_${n1}.txt" -num_threads ${n_threads} -max_hsps 1 -evalue 1e-6
fi
if [[ "$t1" == "prot" && "$t2" == "prot" ]]
then
echo "Running blastp in both directions"
blastp -query $tr1 -db $tr2 -outfmt 6 -out "maps/${n1}${n2}/${n1}_to_${n2}.txt" -num_threads ${n_threads} -max_hsps 1 -evalue 1e-6
blastp -query $tr2 -db $tr1 -outfmt 6 -out "maps/${n1}${n2}/${n2}_to_${n1}.txt" -num_threads ${n_threads} -max_hsps 1 -evalue 1e-6
fi