-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathmotifDynAna
executable file
·257 lines (232 loc) · 13.1 KB
/
motifDynAna
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
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
#!/bin/bash
#PBS -l nodes=1:ppn=4
GENOME="mm9"
PROCESSOR=1
NMOTIFS=50
LENGTH="8,10,12" ## recommended 7,8,9,10,11,12,13,14 (Andersson et al. 2014)
REDUCE_THRESHOLD="0.6" ## argument in compareMotifs.pl (-reduceThreshold) ## recommended 0.75 (Andersson et al. 2014)
#INFO="1.5" ## argument in compareMotifs.pl (-info)
#MIN_T=50 ## argument in compareMotifs.pl (-minT)
PVALUE="0.01" ## recommended 1e-15 (Andersson et al. 2014)
MIN_P_T=0 ## recommended 3 (Andersson et al. 2014)
B=100 ## recommended 30 (Andersson et al. 2014)
S=25 ## recommended 100 (Andersson et al. 2014)
SIZE=200
#### usage ####
usage() {
echo Program: "motifDynAna (compute dynamics in the enrichment of motifs across multiple samples)"
echo Author: BRIC, University of Copenhagen, Denmark
echo Version: 1.0
echo Contact: [email protected]
echo "Usage: motifDynAna -i <dir> -o <dir> [OPTIONS]"
echo "Options:"
echo " -i <dir> [directory containing result from prior run of motifAna script]"
echo " [if multiple, seperate them by a comma]"
echo " -o <dir> [output directory]"
echo "[OPTIONS]"
echo " -g <string> [genome for which to perform the analysis (mm9, hg19, mm10, hg38, danRer7; default: mm9)]"
echo " -p <int> [number of processors to use (default: 1)]"
echo " -c <int> [number of top motifs to analyze from each sample (default: 50)]"
echo " -t <string> [motif length (default=8,10,12)]"
echo " -u <float> [similarity threshold used to remove similar motifs (default: 0.6; **OBSOLETE**)]"
echo " -v <float> [remove motifs with information content less than # (default: not used; **OBSOLETE**)]"
echo " -w <int> [remove motifs with less than # number of target instances (default: not used; **OBSOLETE**)]"
echo " -x <float> [p-value cutoff (default: 0.01)]"
echo " -y <float> [remove motifs with target percentage less than # (default: 0)]"
echo " -z <float> [remove motifs with background percentage greater than # (default: 100)]"
echo " -n <int> [number of motifs to optimize for each motif length (default: 25)]"
echo " -b <file> [custom background file]"
echo " -d <int> [size of region (default: 200bp)]"
echo " [The size of the region used for motif finding is important. If analyzing ChIP-Seq peaks from a transcription factor,]"
echo " [Chuck would recommend 50 bp for establishing the primary motif bound by a given transcription factor and 200 bp for ]"
echo " [finding both primary and 'co-enriched' motifs for a transcription factor. When looking at histone marked regions, ]"
echo " [500-1000 bp is probably a good idea (i.e. H3K4me or H3/H4 acetylated regions). In theory, HOMER can work with very ]"
echo " [large regions (i.e. 10kb), but with the larger the regions comes more sequence and longer execution time. These ]"
echo " [regions will be based off the center of the peaks. If you prefer an offset, you can specify '-size -300,100' to ]"
echo " [search a region of size 400 that is centered 100 bp upstream of the peak center (useful if doing motif finding on ]"
echo " [putative TSS regions). If you have variable length regions, use the option '-size given' and HOMER will use the ]"
echo " [exact regions that were used as input.]"
echo " -l <string> [name of motifs that must be included in the final output]"
echo " [if multiple, separate them by a comma]"
echo " -M <file> [file containing TF motifs which must be included in the output]"
echo " -h [help]"
echo
exit 0
}
#### parse options ####
while getopts i:o:g:p:c:t:u:v:w:x:y:z:n:b:d:l:M:h ARG; do
case "$ARG" in
i) INDIR=$OPTARG;;
o) OUTDIR=$OPTARG;;
g) GENOME=$OPTARG;;
p) PROCESSOR=$OPTARG;;
c) NMOTIFS=$OPTARG;;
t) LENGTH=$OPTARG;;
u) REDUCE_THRESHOLD=$OPTARG;;
v) INFO=$OPTARG;;
w) MIN_T=$OPTARG;;
x) PVALUE=$OPTARG;;
y) MIN_P_T=$OPTARG;;
z) B=$OPTARG;;
n) S=$OPTARG;;
b) BKG_FILE=$OPTARG;;
d) SIZE=$OPTARG;;
l) MUST_INCLUDE_MOTIF=$OPTARG;;
M) MUST_INCLUDE_MOTIF_FILE=$OPTARG;;
h) HELP=1;;
esac
done
## usage, if necessary file and directories are given/exist
if [ -z "$INDIR" -o -z "$OUTDIR" -o "$HELP" ]; then
usage
fi
###################
#helperfunction
function wait_for_jobs_to_finish {
for job in `jobs -p`
do
echo $job
wait $job
done
echo $1
}
###############
<<"COMMENT1"
COMMENT1
echo -n "Create directory structure... "
if [ ! -d "$OUTDIR" ]; then
mkdir -p $OUTDIR
fi
if [ ! -d "$OUTDIR/motifs_logo" ]; then
mkdir -p $OUTDIR/motifs_logo
fi
echo "done"
echo -n "Populating files based on input genome, $GENOME (`date`).. "
if [ "$GENOME" == "mm9" ]; then
GENOME_FILE="/home/pundhir/project/genome_annotations/mouse.mm9.genome"
REPEAT_FILE="/home/pundhir/project/genome_annotations/mouse.mm9.simpleRepeat.gz"
GENOME_MOTIF="mm9r"
elif [ "$GENOME" == "hg19" ]; then
GENOME_FILE="/home/pundhir/project/genome_annotations/human.hg19.genome"
REPEAT_FILE="/home/pundhir/project/genome_annotations/human.hg19.simpleRepeat.gz"
GENOME_MOTIF="hg19r"
elif [ "$GENOME" == "mm10" ]; then
GENOME_FILE="/home/pundhir/project/genome_annotations/mouse.mm10.genome"
REPEAT_FILE="/home/pundhir/project/genome_annotations/mouse.mm10.simpleRepeat.gz"
GENOME_MOTIF="mm10r"
elif [ "$GENOME" == "hg38" ]; then
GENOME_FILE="/home/pundhir/project/genome_annotations/human.hg38.genome"
REPEAT_FILE="/home/pundhir/project/genome_annotations/human.hg38.simpleRepeat.gz"
GENOME_MOTIF="hg38r"
elif [ "$GENOME" == "danRer7" ]; then
GENOME_FILE="/home/pundhir/project/genome_annotations/zebrafish.danRer7.genome"
REPEAT_FILE="/home/pundhir/project/genome_annotations/zebrafish.danRer7.simpleRepeat.gz"
GENOME_MOTIF="danRer7r"
else
echo "Presently the program only support analysis for mm9, hg19 or danRer7"
echo
usage
fi
echo done
## parse input directories
oIFS=$IFS
IFS=","
INDIRS=($INDIR)
IFS=$oIFS
<<"COMMENT"
COMMENT
for CLASS in DENOVO KNOWN; do
echo -n "Parse through each input directory to create a list of all $CLASS motifs found at any motifAna run (`date`).. "
for (( i=0; i<${#INDIRS[@]}; i++)); do
ID=$(echo ${INDIRS[$i]} | perl -ane '$_=~s/\/$//g; $_=~s/^.*\///g; print $_;')
head -n $NMOTIFS ${INDIRS[$i]}/"$CLASS"_MOTIFS.TXT | perl -ane 'print "'$ID'\t$F[0]\t$F[1]\t$F[2]\t$F[3]\t$F[4]\t$F[5]\n";'
#done > $OUTDIR/"$CLASS"_MOTIFS.TXT
done | perl -ane '$id=$F[1]; $id=~s/^.*BestGuess\://g; $id=~s/\s+.*//g; $score=$id; $score=~s/^.*\(//g; $score=~s/\)//g; $id=~s/\(.*//g; if(defined($seen{$id}) && $score > $seen{$id}{'score'}) { $seen{$id}{'score'}=$score; $seen{$id}{'entry'}=$_; } elsif(!defined($seen{$id})) { $seen{$id}{'score'}=$score; $seen{$id}{'entry'}=$_; } END { foreach(keys(%seen)) { print "$seen{$_}{'entry'}"; } }' > $OUTDIR/"$CLASS"_MOTIFS.TXT
for MOTIF_FILE in $(cut -f 6 $OUTDIR/"$CLASS"_MOTIFS.TXT); do
cat $MOTIF_FILE
done > $OUTDIR/"$CLASS"_MOTIFS.MOTIFS
echo "done"
## check, if any "must include" motifs are provided
if [ ! -z "$MUST_INCLUDE_MOTIF" -a -s "$MUST_INCLUDE_MOTIF_FILE" ]; then
zless $MUST_INCLUDE_MOTIF_FILE >> $OUTDIR/"$CLASS"_MOTIFS.MOTIFS
fi
echo -n "Parse through each input directory to compute enrichment of motifs (`date`).. "
ADDITIONAL_ARGUMENT=""
echo -n "Determine all TFBS motifs enriched in the input regions... "
#if [ ! -z "$INFO" ]; then
# ADDITIONAL_ARGUMENT="$ADDITIONAL_ARGUMENT -info $INFO"
#fi
#if [ ! -z "$MIN_T" ]; then
# ADDITIONAL_ARGUMENT="$ADDITIONAL_ARGUMENT -minT $MIN_T"
#fi
for (( i=0; i<${#INDIRS[@]}; i++)); do
ID=$(echo ${INDIRS[$i]} | perl -ane '$_=~s/\/$//g; $_=~s/^.*\///g; print $_;')
if [ ! -d "${INDIRS[$i]}/$CLASS'_motifDynAna'" ]; then
mkdir -p ${INDIRS[$i]}/$CLASS"_motifDynAna"
fi
if [ -z "$BKG_FILE" ]; then
#findMotifsGenome.pl ${INDIRS[$i]}/REGIONS_INTEREST.bed $GENOME_MOTIF ${INDIRS[$i]}/$CLASS"_motifDynAna" -mcheck $OUTDIR/"$CLASS"_MOTIFS.MOTIFS -mknown $OUTDIR/"$CLASS"_MOTIFS.MOTIFS -p $PROCESSOR -nomotif -len $LENGTH -reduceThresh $REDUCE_THRESHOLD -S $S $ADDITIONAL_ARGUMENT -size $SIZE
findMotifsGenome.pl ${INDIRS[$i]}/REGIONS_INTEREST.bed $GENOME_MOTIF ${INDIRS[$i]}/$CLASS"_motifDynAna" -mcheck $OUTDIR/"$CLASS"_MOTIFS.MOTIFS -mknown $OUTDIR/"$CLASS"_MOTIFS.MOTIFS -p $PROCESSOR -nomotif -len $LENGTH -S $S $ADDITIONAL_ARGUMENT -size $SIZE -seqlogo
else
#findMotifsGenome.pl ${INDIRS[$i]}/REGIONS_INTEREST.bed $GENOME_MOTIF ${INDIRS[$i]}/$CLASS"_motifDynAna" -mcheck $OUTDIR/"$CLASS"_MOTIFS.MOTIFS -mknown $OUTDIR/"$CLASS"_MOTIFS.MOTIFS -p $PROCESSOR -nomotif -len $LENGTH -reduceThresh $REDUCE_THRESHOLD -S $S $ADDITIONAL_ARGUMENT -bg $BKG_FILE -size $SIZE
findMotifsGenome.pl ${INDIRS[$i]}/REGIONS_INTEREST.bed $GENOME_MOTIF ${INDIRS[$i]}/$CLASS"_motifDynAna" -mcheck $OUTDIR/"$CLASS"_MOTIFS.MOTIFS -mknown $OUTDIR/"$CLASS"_MOTIFS.MOTIFS -p $PROCESSOR -nomotif -len $LENGTH -S $S $ADDITIONAL_ARGUMENT -bg $BKG_FILE -size $SIZE -seqlogo
fi
if [ -z "$BKG_FILE" ]; then
TARGET_TOTAL=$(grep "^>" ${INDIRS[$i]}/target.fa | wc -l)
BKG_TOTAL=$(grep "^>" ${INDIRS[$i]}/background.fa | wc -l)
else
TARGET_TOTAL=$(cat ${INDIRS[$i]}/REGIONS_INTEREST.bed | wc -l)
BKG_TOTAL=$(cat $BKG_FILE | wc -l)
fi
if [ -f "${INDIRS[$i]}/knownResults.html" ]; then
#grep -v Name ${INDIRS[$i]}/$CLASS"_motifDynAna"/knownResults.txt | perl -ane '$pval=sprintf("%0.20f", $F[2]); $my_pval=sprintf("%0.20f", '$PVALUE'); $F[6]=~s/\%//g; $F[8]=~s/\%//g; if($pval < $my_pval && $F[6]>='$MIN_P_T' && $F[8]<='$B') { chomp($_); $_=~s/\%//g; print "'$ID'\t$_\t'$TARGET_TOTAL'\t'$BKG_TOTAL'\n"; }'
## we want to include all motifs irrespective of their p-values
grep -v Name ${INDIRS[$i]}/$CLASS"_motifDynAna"/knownResults.txt | perl -ane 'chomp($_); $_=~s/\%//g; print "'$ID'\t$_\t'$TARGET_TOTAL'\t'$BKG_TOTAL'\n";'
fi
done | sort -k 1,1 -k 2,2 > $OUTDIR/"$CLASS"_MOTIF_ENRICHMENT_DYNAMICS.TXT
echo "done"
echo -n "Plot motif dynamics (`date`).. "
if [ ! -z "$MUST_INCLUDE_MOTIF" -a -s "$MUST_INCLUDE_MOTIF_FILE" ]; then
motifDynAna.R -i $OUTDIR/"$CLASS"_MOTIF_ENRICHMENT_DYNAMICS.TXT -o $OUTDIR/"$CLASS"_motif_dynamics.pdf -l $MUST_INCLUDE_MOTIF
else
motifDynAna.R -i $OUTDIR/"$CLASS"_MOTIF_ENRICHMENT_DYNAMICS.TXT -o $OUTDIR/"$CLASS"_motif_dynamics.pdf
fi
echo "done"
echo -n "Plot logo for significant motifs (`date`).. "
for motifInfo in $(cut -f 1,2,6 $OUTDIR/"$CLASS"_MOTIFS.TXT | perl -ane '$F[0]=~s/[^a-zA-Z0-9]+/_/g; $F[1]=~s/\(.*//g; $F[1]=~s/[^a-zA-Z0-9]+/_/g; $key="$F[0]_$F[1]"; if(!$seen{$key}) { print "$key#$F[2]\n"; $seen{$key}=1; }'); do
logoFile=$(echo $motifInfo | sed -E 's/\#.*//g');
motifFile=$(echo $motifInfo | sed -E 's/^.*\#//g');
motif2Logo.pl $motifFile -pdf -o $OUTDIR/motifs_logo/$logoFile;
done
if [ ! -z "$MUST_INCLUDE_MOTIF" -a -s "$MUST_INCLUDE_MOTIF_FILE" ]; then
oIFS=$IFS
IFS=","
MUST_INCLUDE_MOTIFS=($MUST_INCLUDE_MOTIF)
IFS=$oIFS
for (( i=0; i<${#MUST_INCLUDE_MOTIFS[@]}; i++)); do
extractMotifs -i $OUTDIR/"$CLASS"_MOTIFS.MOTIFS -j ${MUST_INCLUDE_MOTIFS[$i]} | motif2Logo.pl - -pdf -o $OUTDIR/motifs_logo/${MUST_INCLUDE_MOTIFS[$i]}
done
fi
echo "done"
done
## OLD
<<"COMMENT"
echo -n "Parse through each input directory to compute enrichment of motifs (`date`).. "
for (( i=0; i<${#INDIRS[@]}; i++)); do
ID=$(echo ${INDIRS[$i]} | perl -ane '$_=~s/\/$//g; $_=~s/^.*\///g; print $_;')
if [ -f "${INDIRS[$i]}/target.fa" -a -f "${INDIRS[$i]}/background.fa" ]; then
homer2 find -i ${INDIRS[$i]}/target.fa -m $OUTDIR/DENOVO_MOTIFS.MOTIFS 2> /dev/null | perl -ane '$id="$F[0]_$F[3]"; print "$id\n";' | sort | uniq > $OUTDIR/$ID.target.fa.homer
homer2 find -i ${INDIRS[$i]}/background.fa -m $OUTDIR/DENOVO_MOTIFS.MOTIFS 2> /dev/null | perl -ane '$id="$F[0]_$F[3]"; print "$id\n";' | sort | uniq > $OUTDIR/$ID.background.fa.homer
TARGET_TOTAL=$(grep "^>" ${INDIRS[$i]}/target.fa | wc -l)
BKG_TOTAL=$(grep "^>" ${INDIRS[$i]}/background.fa | wc -l)
for MOTIF in $(grep "^>" $OUTDIR/DENOVO_MOTIFS.MOTIFS | cut -f 2); do
TARGET_COUNT=$(grep $MOTIF $OUTDIR/$ID.target.fa.homer | wc -l)
BKG_COUNT=$(grep $MOTIF $OUTDIR/$ID.background.fa.homer | wc -l)
echo "$MOTIF $TARGET_COUNT $TARGET_TOTAL $BKG_COUNT $BKG_TOTAL"
done | perl -ane 'if($F[1]==0) { $F[1]=1; } if($F[3]==0) { $F[3]=1; } $obs=$F[1]/$F[2]; $exp=$F[3]/$F[4]; $log_odd_score=sprintf("%0.2f", log($obs/$exp)/log(2)); $obs=sprintf("%0.2f", $obs); $exp=sprintf("%0.2f", $exp); print "$F[0]\t$obs\t$exp\t$log_odd_score\t'$ID'\t$F[1]\t$F[2]\t$F[3]\t$F[4]\n";'
rm $OUTDIR/$ID.target.fa.homer
rm $OUTDIR/$ID.background.fa.homer
fi
done > $OUTDIR/MOTIF_ENRICHMENT_DYNAMICS.TXT
echo "done"
COMMENT