Merge remote-tracking branch 'origin/master' into release

README.md

@@ -1,29 +1,23 @@
-# UCSF PyEM
-UCSF PyEM is a collection of Python modules and command-line utilities for electron microscopy of biological samples.
+# UCSF pyem
+UCSF pyem is a collection of Python modules and command-line utilities for electron microscopy of biological samples.
+
+Documentation for the programs can be found in their usage text, comments in code, and in the Wiki of this repository.
 
 The entire collection is licensed under the terms of the GNU Public License, version 3 (GPLv3).
 
-Copyright information is listed within each individual file. Current copyright holders include:
- * Eugene Palovcak (UCSF)
- * Daniel Asarnow (UCSF)
+# How to cite
 
-Documentation for the programs can be found in their usage text, comments in code, and in the Wiki of this repository.
+Please cite the DOI for the UCSF pyem code repository.
 
-## Programs
- 1. `projection_subtraction.py` - Perform projection subtraction using per-particle FRC normalization.
- + `recenter.py` - Recenter particles on the center-of-mass of corresponding 2D class-averages.
- + `angdist.py` - Graph angular distributions on polar scatter plots. Supports particle subset selection.
- + `pyem/star.py` - Alter .star files. Supports dropping arbitrary fields, Euler angles, etc.
- + `project.py` - Project a map according to .star file entries (angles, CTF, etc.).
- + `csparc2star.py` - Convert Cryosparc metadata files to Relion .star format.
+# Installation
 
-## Library modules
- 1. `pyem/mrc.py` - Simple, standalone MRC I/O functions.
- + `pyem/star.py` - Parse and write .star files. Uses pandas.DataFrame as a backend.
+To install UCSF pyem, please follow the
+[instructions](https://github.com/asarnow/pyem/wiki/Install-pyem-with-Miniconda) in the project wiki.
 
-## Other files
- 1. `activate` - Place in `EMAN2/bin` to turn EMAN2 into a Python virtual environment.
+# Exporting from cryoSPARC to Relion
 
+The most popular feature of UCSF pyem is exporting particle metadata from cryoSPARC.
+Detailed [instructions](https://github.com/asarnow/pyem/wiki/Export-from-cryoSPARC-v2) can be found in the project wiki.
 
-(C) 2016 Daniel Asarnow  
+(C) 2016-2019 Daniel Asarnow  
 University of California, San Francisco

csparc2star.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python2.7
+#!/usr/bin/env python
 # Copyright (C) 2016 Daniel Asarnow
 # University of California, San Francisco
 #
@@ -41,9 +41,8 @@ def main(args):
         try:
             df = metadata.parse_cryosparc_2_cs(cs, passthroughs=args.input[1:], minphic=args.minphic, boxsize=args.boxsize, swapxy=args.swapxy)
         except (KeyError, ValueError) as e:
-            log.error(e.message)
-            log.error("A passthrough file may be required (check inside the cryoSPARC 2+ job directory)")
-            log.debug(e, exc_info=True)
+            log.error(e, exc_info=True)
+            log.error("Required fields could not be mapped. Are you using the right input file(s)?")
             return 1
     else:
         log.debug("Detected CryoSPARC 0.6.5 .csv file")
@@ -58,10 +57,8 @@ def main(args):
 
     if args.copy_micrograph_coordinates is not None:
         coord_star = pd.concat(
-            (star.parse_star(inp, keep_index=False) for inp in
+            (star.parse_star(inp, keep_index=False, augment=True) for inp in
              glob(args.copy_micrograph_coordinates)), join="inner")
-        star.augment_star_ucsf(coord_star)
-        star.augment_star_ucsf(df)
         key = star.merge_key(df, coord_star)
         log.debug("Coordinates merge key: %s" % key)
         if args.cached or key == star.Relion.IMAGE_NAME:
@@ -79,7 +76,7 @@ def main(args):
         df = star.transform_star(df, r, inplace=True)
 
     # Write Relion .star file with correct headers.
-    star.write_star(args.output, df)
+    star.write_star(args.output, df, resort_records=True)
     log.info("Output fields: %s" % ", ".join(df.columns))
     return 0
 

ctf2star.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python3
+#!/usr/bin/env python
 # Copyright (C) 2019 Daniel Asarnow
 # University of California, San Francisco
 #

map.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python2.7
+#!/usr/bin/env python
 # Copyright (C) 2017 Daniel Asarnow
 # University of California, San Francisco
 #
@@ -30,20 +30,15 @@
 from pyem import vop
 from scipy.ndimage import affine_transform
 from scipy.ndimage import shift
-from scipy.ndimage import zoom
+import warnings
+warnings.filterwarnings('ignore', '.*output shape of zoom.*')
 
 
 def main(args):
     log = logging.getLogger(__name__)
-    log.setLevel(logging.INFO)
     hdlr = logging.StreamHandler(sys.stdout)
-    if args.quiet:
-        hdlr.setLevel(logging.ERROR)
-    elif args.verbose:
-        hdlr.setLevel(logging.INFO)
-    else:
-        hdlr.setLevel(logging.WARN)
     log.addHandler(hdlr)
+    log.setLevel(logging.getLevelName(args.loglevel.upper()))
 
     data, hdr = read(args.input, inc_header=True)
     final = None
@@ -80,11 +75,32 @@ def main(args):
         args.apix = hdr["xlen"] / hdr["nx"]
         log.info("Using computed pixel size of %f Angstroms" % args.apix)
 
-    if args.target and args.matrix:
+    if args.apix_out is not None:
+        if args.scale is not None:
+            log.warn("--apix-out supersedes --scale")
+        args.scale = args.apix / args.apix_out
+    elif args.scale is not None:
+            args.apix_out = args.apix / args.scale
+    elif args.boxsize is not None:
+        args.scale = box[0] / np.double(args.boxsize)
+
+    if args.apix_out is None:
+        args.apix_out = args.apix
+
+    if args.boxsize is None:
+        if args.scale is None:
+            args.boxsize = box[0]
+            args.scale = 1
+        else:
+            args.boxsize = np.int(box[0] * args.scale)
+
+    log.info("Volume will be scaled by %f to size %d @ %f A/px" % (args.scale, args.boxsize, args.apix_out))
+
+    if args.target and args.transform:
         log.warn("Target pose transformation will be applied after explicit matrix")
-    if args.euler is not None and (args.target is not None or args.matrix is not None):
+    if args.euler is not None and (args.target is not None or args.transform is not None):
         log.warn("Euler transformation will be applied after target pose transformation")
-    if args.translate is not None and (args.euler is not None or args.target is not None or args.matrix is not None):
+    if args.translate is not None and (args.euler is not None or args.target is not None or args.transform is not None):
         log.warn("Translation will be applied after other transformations")
 
     if args.origin is not None:
@@ -98,17 +114,26 @@ def main(args):
         args.origin = center
         log.info("Origin set to box center, %s" % (args.origin * args.apix))
 
-    if not (args.target is None and args.euler is None and args.matrix is None and args.boxsize is None) \
+    if not (args.target is None and args.euler is None and args.transform is None and args.boxsize is None) \
             and vop.ismask(data) and args.spline_order != 0:
         log.warn("Input looks like a mask, --spline-order 0 (nearest neighbor) is recommended")
 
-    if args.matrix is not None:
+    if args.transform is not None:
         try:
-            r = np.array(json.loads(args.matrix))
+            args.transform = np.array(json.loads(args.transform))
         except:
-            log.error("Matrix format is incorrect")
+            log.error("Transformation matrix must be in JSON/Numpy format")
             return 1
-        data = vop.resample_volume(data, r=r, t=None, ori=None, order=args.spline_order)
+        r = args.transform[:, :3]
+        if args.transform.shape[1] == 4:
+            t = args.transform[:, -1] / args.apix
+            t = r.dot(args.origin) + t - args.origin
+            t = -r.T.dot(t)
+        else:
+            t = 0
+        log.debug("Final rotation: %s" % str(r).replace("\n", "\n" + " " * 16))
+        log.debug("Final translation: %s (%f px)" % (str(t), np.linalg.norm(t)))
+        data = vop.resample_volume(data, r=r, t=t, ori=None, order=args.spline_order, invert=args.invert)
 
     if args.target is not None:
         try:
@@ -122,7 +147,9 @@ def main(args):
         r = vec2rot(args.target)
         t = np.linalg.norm(args.target)
         log.info("Euler angles are %s deg and shift is %f px" % (np.rad2deg(rot2euler(r)), t))
-        data = vop.resample_volume(data, r=r, t=args.target, ori=ori, order=args.spline_order, invert=args.target_invert)
+        log.debug("Final rotation: %s" % str(r).replace("\n", "\n" + " " * 16))
+        log.debug("Final translation: %s (%f px)" % (str(t), np.linalg.norm(t)))
+        data = vop.resample_volume(data, r=r, t=args.target, ori=ori, order=args.spline_order, invert=args.invert)
 
     if args.euler is not None:
         try:
@@ -144,19 +171,17 @@ def main(args):
         args.translate -= args.origin
         data = shift(data, -args.translate, order=args.spline_order)
 
-    if args.boxsize is not None:
-        args.boxsize = np.double(args.boxsize)
-        data = zoom(data, args.boxsize / box, order=args.spline_order)
-        args.apix = args.apix * box[0] / args.boxsize
-
     if final is None:
         final = data
 
     if args.final_mask is not None:
         final_mask = read(args.final_mask)
         final *= final_mask
 
-    write(args.output, final, psz=args.apix)
+    if args.scale != 1 or args.boxsize != box[0]:
+        final = vop.resample_volume(final, scale=args.scale, output_shape=args.boxsize, order=args.spline_order)
+
+    write(args.output, final, psz=args.apix_out)
     return 0
 
 
@@ -167,7 +192,7 @@ def main(args):
     parser.add_argument("input", help="Input volume (MRC file)")
     parser.add_argument("output", help="Output volume (MRC file)")
     parser.add_argument("--apix", "--angpix", "-a", help="Pixel size in Angstroms", type=float)
-    parser.add_argument("--mask", help="Final mask to apply after any operations", dest="final_mask")
+    parser.add_argument("--mask", help="Final mask (applied before scaling)", dest="final_mask")
     parser.add_argument("--transpose", help="Swap volume axes order", metavar="a1,a2,a3")
     parser.add_argument("--normalize", "-n", help="Convert map densities to Z-scores", action="store_true")
     parser.add_argument("--reference", "-r", help="Normalization reference volume (MRC file)")
@@ -176,15 +201,17 @@ def main(args):
     parser.add_argument("--pfac", help="Padding factor for 3D FFT", type=int, default=2)
     parser.add_argument("--origin", help="Origin coordinates in Angstroms (volume center by default)", metavar="x,y,z")
     parser.add_argument("--target", help="Target pose (view axis and origin) coordinates in Angstroms", metavar="x,y,z")
-    parser.add_argument("--target-invert", help="Undo target pose transformation", action="store_true")
+    parser.add_argument("--invert", help="Invert the transformation", action="store_true")
+    parser.add_argument("--target-invert", action="store_true", dest="invert", help=argparse.SUPPRESS)
     parser.add_argument("--euler", help="Euler angles in degrees (Relion conventions)", metavar="phi,theta,psi")
     parser.add_argument("--translate", help="Translation coordinates in Angstroms", metavar="x,y,z")
-    parser.add_argument("--matrix",
+    parser.add_argument("--transform",
                         help="Transformation matrix (3x3 or 3x4 with translation in Angstroms) in Numpy/json format")
-    parser.add_argument("--boxsize", help="Set the output box dimensions, accounting for pixel size", type=int)
+    parser.add_argument("--boxsize", help="Output box size (pads/crops with --scale or --apix-out, otherwise scales)", type=int)
+    parser.add_argument("--scale", help="Scale factor for output pixel size", type=float)
+    parser.add_argument("--apix-out", help="Pixel size in output (similar to --scale)", type=float)
     parser.add_argument("--spline-order",
                         help="Order of spline interpolation (0 for nearest, 1 for trilinear, default is cubic)",
                         type=int, default=3, choices=np.arange(6))
-    parser.add_argument("--quiet", "-q", help="Print errors only", action="store_true")
-    parser.add_argument("--verbose", "-v", help="Print info messages", action="store_true")
+    parser.add_argument("--loglevel", "-l", type=str, default="WARNING", help="Logging level and debug output")
     sys.exit(main(parser.parse_args()))

mask.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python2.7
+#!/usr/bin/env python
 # Copyright (C) 2017 Daniel Asarnow
 # University of California, San Francisco
 #
@@ -40,7 +40,7 @@ def main(args):
         base_map = read(args.base_map, inc_header=False)
         base_mask = binarize_volume(base_map, args.threshold, minvol=args.minvol, fill=args.fill)
         total_width = args.extend + args.edge_width
-        excl_mask = binary_dilate(mask, args.extend+args.edge_width, strel=args.relion)
+        excl_mask = binary_dilate(mask, total_width, strel=args.relion)
         base_mask = binary_dilate(base_mask, args.extend, strel=args.relion)
         base_mask = base_mask &~ excl_mask
         if args.overlap > 0:
@@ -53,7 +53,7 @@ def main(args):
         se = binary_sphere(args.extend, False)
         mask = binary_closing(mask, structure=se, iterations=1)
     final = mask.astype(np.single)
-    if args.edge_width is not None:
+    if args.edge_width != 0:
         dt = distance_transform_edt(~mask)  # Compute *outward* distance transform of mask.
         idx = (dt <= args.edge_width) & (dt > 0)  # Identify edge points by distance from mask.
         x = np.arange(1, args.edge_width + 1)  # Domain of the edge profile.
@@ -84,7 +84,7 @@ def main(args):
     parser.add_argument("--extend", "-e", help="Structuring element size for dilating initial mask",
                         type=int, default=0)
     parser.add_argument("--edge-width", "-w", help="Width for soft edge",
-                        type=int)
+                        type=int, default=0)
     parser.add_argument("--edge-profile", "-p", help="Soft edge profile type",
                         choices=["sinusoid"],
                         default="sinusoid")

par2star.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python2.7
+#!/usr/bin/env python
 # Copyright (C) 2017 Daniel Asarnow
 # University of California, San Francisco
 #
@@ -47,7 +47,7 @@ def main(args):
     df = df.sort_values(by="C")
 
     df = metadata.par2star(df, data_path=args.stack, apix=args.apix, cs=args.cs,
-                           ac=args.ac, kv=args.voltage, invert_eulers=args.relion)
+                           ac=args.ac, kv=args.voltage, invert_eulers=args.invert_eulers)
 
     if args.cls is not None:
         df = star.select_classes(df, args.cls)
@@ -68,7 +68,8 @@ def main(args):
     parser.add_argument("--voltage", "--kv", "-v", help="Acceleration voltage (kV)", type=float)
     parser.add_argument("--min-occ", help="Minimum occupancy for inclusion in output", type=float)
     parser.add_argument("--class", "-c", help="Classes to preserve in output", action="append", dest="cls")
-    parser.add_argument("--relion", help="Use if .par file came from importing .star in cisTEM", action="store_false")
+    parser.add_argument("--relion", help=argparse.SUPPRESS, action="store_true")
+    parser.add_argument("--invert-eulers", help="Invert Euler angles (generally unnecessary)", action="store_true")
     parser.add_argument("--loglevel", "-l", type=str, default="WARNING", help="Logging level and debug output")
     sys.exit(main(parser.parse_args()))
 

pose.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python2.7
+#!/usr/bin/env python
 # -*- coding: utf-8 -*-
 # Copyright (C) 2018 Daniel Asarnow
 # University of California, San Francisco

project.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python2.7
+#!/usr/bin/env python
 # Copyright (C) 2016 Daniel Asarnow, Eugene Palovcak
 # University of California, San Francisco
 #

projection_subtraction.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python2.7
+#!/usr/bin/env python
 # -*- coding: utf-8 -*-
 # Copyright (C) 2015-2018 Daniel Asarnow, Eugene Palovcak
 # University of California, San Francisco

pyem/algo/algo.py

@@ -1,4 +1,3 @@
-#!/usr/bin/env python2.7
 # Copyright (C) 2018 Daniel Asarnow
 # University of California, San Francisco
 #

pyem/algo/algo_numba.py

@@ -1,4 +1,3 @@
-#!/usr/bin/env python2.7
 # Copyright (C) 2018 Daniel Asarnow
 # University of California, San Francisco
 #

pyem/ctf.py

@@ -1,4 +1,3 @@
-#!/usr/bin/env python2.7
 # -*- coding: utf-8 -*-
 # Copyright (C) 2018 Daniel Asarnow, Eugene Palovcak
 # University of California, San Francisco