FREALIGN/estimateCTF_CTFFIND4_plotReslim.py at master · leschzinerlab/FREALIGN · GitHub

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
#!/usr/bin/env python

from pylab import *
from matplotlib.colors import LogNorm
import numpy as np
import matplotlib.pyplot as plt
import optparse
from sys import *
import os,sys,re
import linecache

#=========================
def setupParserOptions():
        parser = optparse.OptionParser()
        parser.set_usage("%prog --starfile=<relion_star_file>")
        parser.add_option("--input",dest="input",type="string",metavar="FILE",
                help="Output ctf file from estimateCTF_CTFFIND4.py ('ctf_param.txt')")
        parser.add_option("--binsize",dest="bin",type="int",metavar="INT",default=5,
                help="Optional: bin size for histogram of euler angles (Default=1 Angstrom)")
        parser.add_option("-d", action="store_true",dest="debug",default=False,
                help="debug")
        options,args = parser.parse_args()

        if len(args) > 0:
                parser.error("Unknown commandline options: " +str(args))

        if len(sys.argv) < 3:
                parser.print_help()
                sys.exit()
        params={}
        for i in parser.option_list:
                if isinstance(i.dest,str):
                        params[i.dest] = getattr(options,i.dest)
        return params

#==============================
def checkConflicts(params):

        if not os.path.exists(params['input']):
            print 'Error: File %s does not exist' %(params['input'])
            sys.exit()

#===============================
def getRelionColumnIndex(star,rlnvariable):

    counter=50
    i=1

    while i<=50:

        line=linecache.getline(star,i)

        if len(line)>0:
            if len(line.split())>1:
                if line.split()[0] == rlnvariable:
                    return line.split()[1][1:]

        i=i+1

#===============================
def plot_single(star,colnum,debug,binsize):

    #open star file
    f1=open(star,'r')

    #remove temporary file
    tmp='tmp_relion_star.star'
    if os.path.exists(tmp):
        os.remove(tmp)

    #open for writing into new tmp file without header
    o1=open(tmp,'w')

    for line in f1:
        if l.split()[-1] == 'ResLimit':
            continue
        o1.write(line)

    o1.close()
    f1.close()

    usecolumn=int(colnum-1)

    eulers=np.loadtxt(tmp,usecols=[usecolumn])

    #Get number of particles
    tot=len(open(tmp,'r').readlines())

    #Create bins:
    bins=np.arange(3,30,params['bin'])

    #Plot histogram
    plt.hist(eulers,bins=bins)
    plt.title("Micrograph resolution histogram")
    plt.xlabel("Resolution of CTFFIND4 confidence limit (Angstroms)")
    plt.ylabel("Number of micrographs")
    plt.show()

#==============================
if __name__ == "__main__":

        params=setupParserOptions()

        #Check that file exists & relion euler designation is real
        checkConflicts(params)

        #Get column number for euler designation
        plot_single(params['input'],6,params['debug'],params['bin'])