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
|
#!/bin/bash
#
#
# Plot the results from ctffind using gnuplot
#
# Alexis Rohou, June 2014
#
# Copyright 2014 Howard Hughes Medical Institute. All rights reserved.
# Use is subject to Janelia Farm Research Campus Software Copyright 1.1
# license terms ( http://license.janelia.org/license/jfrc_copyright_1_1.html )
#
# Last modification: Tapu Shaik, April 2018
# Alexis Rohou, Dec 2018 (adding score to plot)
# Parse arguments
if [ $# -ne 1 ]
then
echo "Usage: `basename $0` output_from_ctffind_avrot.txt"
exit 65
fi
input_fn=$1
output_fn=${1%%.???}.pdf
#input_summary_fn=$(ls ${1%%_avrot.???}.${1##*.})
input_summary_fn="${1%%_avrot.???}.${1##*.}"
# Check whether gnuplot is available
if ! hash gnuplot 2>/dev/null; then
echo "Gnuplot was not found on your system. Please make sure it is installed."
exit -1
fi
# Check whether gnuplot >= 4.6 is available
gnuplot_version=$( gnuplot --version | awk '{print $2}' )
compare_result=$( echo "$gnuplot_version >= 4.6" | bc )
if [ $compare_result -le 0 ]; then
echo "This script requires gnuplot version >= 4.6, but you have version $gnuplot_version"
exit -1
fi
# Define a neat function
function transpose_ignore_comments {
gawk '
{
if ($1 ~ /^#/) {
print $0
} else {
for (i=1; i<=NF; i++) {
a[NR,i] = $i
}
}
}
NF>p { p = NF }
END {
for(j=1; j<=p; j++) {
str=a[1,j]
for(i=2; i<=NR; i++){
str=str" "a[i,j];
}
print str
}
}' $1
}
# CTFFind outputs data in lines, but gnuplot wants things in columns
transpose_ignore_comments $input_fn > /tmp/tmp.txt
# Let's grab useful values
pixel_size=$(gawk 'match($0,/Pixel size: ([0-9.]*)/,a) {print a[1]}' $input_fn)
number_of_micrographs=$(gawk 'match($0,/Number of micrographs: ([0-9]*)/,a) {print a[1]}' $input_fn)
mic_name=$(grep "Input file:" $input_fn | awk 'BEGIN {FS="Input file: "} {print $2}' | cut -d " " -f1 | sed 's/_/\\\_/g')
lines_per_micrograph=6
if [ -f $input_summary_fn ]; then
# Let's grab values from the summary file
i=0
while read -a myArray
do
if [[ ${myArray[0]} != \#* ]]; then
df_one[++i]=${myArray[1]}
df_two[i]=${myArray[2]}
angast[i]=${myArray[3]}
pshift[i]=${myArray[4]}
score[i]=${myArray[5]}
maxres[i]=${myArray[6]}
fi
done < $input_summary_fn
else
df_one[0]="0"
df_two[0]="0"
angast[0]="0"
pshift[0]="0"
score[0]="0"
maxres[0]="0"
fi
# Run Gnuplot
gnuplot > /dev/null 2>&1 <<EOF
#cat <<EOF > temp.txt
set border linewidth 1.5
if (strstrt(GPVAL_TERMINALS, 'pdfcairo') > 0) {
set terminal pdfcairo size 10.0cm,4.0cm enhanced font 'Arial,9'
set output '$output_fn'
} else {
set terminal postscript size 10,4 enhanced font 'Arial,14'
set output '${output_fn%%.pdf}.ps'
}
# color definitions
set style line 1 lc rgb '#0060ad' lt 1 lw 2 pt 7 # --- blue
set style line 2 lc rgb 'red' lt 1 lw 1 pt 7 fill transparent solid 0.5 # --- red
set style line 3 lc rgb 'orange' lt 3 lw 1 pt 7 # --- orange
set style line 4 lc rgb 'light-blue' lt 1 lw 2 pt 7 # --- light blue
set style line 5 lc rgb 'green' lt 1 lw 2 pt 7
set style line 6 lc rgb 'gray' lt 1 lw 2 pt 7
set xlabel 'Spatial frequency(1/Å)'
set ylabel 'Amplitude (or cross-correlation)'
set yrange [-0.1:1.1]
set key outside
defocus_1_values="${df_one[*]}"
defocus_2_values="${df_two[*]}"
angast_values="${angast[*]}"
pshift_values="${pshift[*]}"
score_values="${score[*]}"
maxres_values="${maxres[*]}"
do for [current_micrograph=1:$number_of_micrographs] {
def_1=sprintf('%.0f Å',word(defocus_1_values,current_micrograph)+0)
def_2=sprintf('%.0f Å',word(defocus_2_values,current_micrograph)+0)
angast=sprintf('%.1f °',word(angast_values,current_micrograph)+0)
pshift=sprintf('%.2f rad',word(pshift_values,current_micrograph)+0)
score=sprintf('%.3f',word(score_values,current_micrograph)+0)
maxres=sprintf('%.2f Å',word(maxres_values,current_micrograph)+0)
# set title 'Micrograph '.current_micrograph." of $number_of_micrographs\n{/*0.5 Defocus 1: ".def_1.' | Defocus 2: '.def_2.' | Azimuth: '.angast.' | Phase shift: '.pshift.' | Score: '.score.'}'
set title '$mic_name, '.current_micrograph." of $number_of_micrographs micrographs\n{/*0.5 Defocus 1: ".def_1.' | Defocus 2: '.def_2.' | Azimuth: '.angast.' | Phase shift: '.pshift.' | Score: '.score.' | MaxRes: '.maxres.'}'
plot '/tmp/tmp.txt' using (\$1):(column(4+(current_micrograph-1)*$lines_per_micrograph)) w lines ls 3 title 'CTF fit', \
'' using (\$1):(column(5+(current_micrograph-1)*$lines_per_micrograph)) w lines ls 1 title 'Quality of fit', \
'' using (\$1):(column(3+(current_micrograph-1)*$lines_per_micrograph)) w lines ls 5 title 'Amplitude spectrum'
}
EOF
# Convert the postscript file to a pdf file if needed
if [ -f ${output_fn%%.pdf}.ps ]; then
ps2pdf ${output_fn%%.pdf}.ps $output_fn
fi
# Fix the title of the PDF file
if hash pdftk 2>/dev/null; then
echo "InfoKey: Title" > /tmp/pdftemp.txt
echo "InfoValue: $1" >> /tmp/pdftemp.txt
mv ${output_fn} /tmp/tmp.pdf
pdftk /tmp/tmp.pdf update_info /tmp/pdftemp.txt output $output_fn
elif hash exiftool 2>/dev/null; then
exiftool -Title="$1" ${output_fn}
else
echo "Neither pdftk nor exiftool were found on your system. Please install one of them to get improved metadata in your output PDF file."
fi
# Make a PNG version
#convert -density 300 -flatten ${output_fn} ${output_fn%%.pdf}.png
# Let the user know where to find results
echo " "
echo "Generated plots of 1D profiles and fits: ${output_fn}"
echo " "
|
