Ordinary diffmaps

efxtools/diffmaps/diffmap.py

@@ -0,0 +1,104 @@
+#!/usr/bin/env python
+"""
+Make an ordinary difference map between columns of data.
+
+The naming convention chosen for inputs is `on` and `off`, such
+that the generated differrence map will be `|F_on| - |F_off|`.
+"""
+
+import argparse
+import numpy as np
+import reciprocalspaceship as rs
+
+from efxtools.diffmaps.weights import compute_weights
+from efxtools.utils.io import subset_to_FSigF
+
+
+def parse_arguments():
+    """Parse commandline arguments"""
+    parser = argparse.ArgumentParser(
+        formatter_class=argparse.RawTextHelpFormatter, description=__doc__
+    )
+
+    # Required arguments
+    parser.add_argument(
+        "-on",
+        "--onmtz",
+        nargs=3,
+        metavar=("mtz", "data_col", "sig_col"),
+        required=True,
+        help=("MTZ to be used as `on` data. Specified as (filename, F, SigF)"),
+    )
+    parser.add_argument(
+        "-off",
+        "--offmtz",
+        nargs=3,
+        metavar=("mtz", "data_col", "sig_col"),
+        required=True,
+        help=("MTZ to be used as `off` data. Specified as (filename, F, SigF)"),
+    )
+    parser.add_argument(
+        "-r",
+        "--refmtz",
+        nargs=2,
+        metavar=("ref", "phi_col"),
+        required=True,
+        help=(
+            "MTZ containing isomorphous phases to be used. "
+            "Specified as (filename, Phi)."
+        ),
+    )
+
+    # Optional arguments
+    parser.add_argument(
+        "-a",
+        "--alpha",
+        type=float,
+        default=0.0,
+        help="alpha value for computing difference map weights (default=0.0)",
+    )
+    parser.add_argument(
+        "-o", "--outfile", default="diffmap.mtz", help="Output MTZ filename"
+    )
+
+    return parser.parse_args()
+
+
+def main():
+
+    # Parse commandline arguments
+    args = parse_arguments()
+    refmtz, phi_col = args.refmtz
+
+    # Read MTZ files
+    onmtz = subset_to_FSigF(*args.onmtz, {args.onmtz[1]: "F", args.onmtz[2]: "SigF"})
+    offmtz = subset_to_FSigF(
+        *args.offmtz, {args.offmtz[1]: "F", args.offmtz[2]: "SigF"}
+    )
+
+    ref = rs.read_mtz(refmtz)
+    ref.rename(columns={phi_col: "Phi"}, inplace=True)
+    ref = ref.loc[:, ["Phi"]]
+    if not isinstance(ref["Phi"].dtype, rs.PhaseDtype):
+        raise ValueError(
+            f"{args.Phi} is not a phases column in {args.mtz2}. Try again."
+        )
+
+    diff = onmtz.merge(offmtz, on=["H", "K", "L"], suffixes=("_on", "_off"))
+    diff["DF"] = diff["F_on"] - diff["F_off"]
+    diff["SigDF"] = np.sqrt((diff["SigF_on"] ** 2) + (diff["SigF_off"] ** 2))
+
+    # Compute weights
+    diff["W"] = compute_weights(diff["DF"], diff["SigDF"], alpha=args.alpha)
+    diff["W"] = diff["W"].astype("Weight")
+
+    # Join with phases and write map
+    common = diff.index.intersection(ref.index).sort_values()
+    diff = diff.loc[common]
+    diff["Phi"] = ref.loc[common, "Phi"]
+    diff.infer_mtz_dtypes(inplace=True)
+    diff.write_mtz(args.outfile)
+
+
+if __name__ == "__main__":
+    main()

efxtools/diffmaps/internaldiffmap.py

@@ -9,6 +9,7 @@
 import gemmi
 
 from efxtools.diffmaps.weights import compute_weights
+from efxtools.utils.io import subset_to_FSigF
 
 
 def parse_arguments():
@@ -74,28 +75,19 @@ def main():
 
     # Parse commandline arguments
     args = parse_arguments()
-    inputmtz, f_col, sigf_col = args.inputmtz
     refmtz, phi_col = args.refmtz
 
     # Read MTZ files
-    mtz = rs.read_mtz(inputmtz)
+    mtz = subset_to_FSigF(
+        *args.inputmtz, {args.inputmtz[1]: "F", args.inputmtz[2]: "SigF"}
+    )
     ref = rs.read_mtz(refmtz)
 
     # Canonicalize column names
-    mtz.rename(columns={f_col: "F", sigf_col: "SigF"}, inplace=True)
-    mtz = mtz[["F", "SigF"]]
     ref.rename(columns={phi_col: "Phi"}, inplace=True)
     ref = ref[["Phi"]]
 
     # Error checking of datatypes
-    if not isinstance(mtz["F"].dtype, rs.StructureFactorAmplitudeDtype):
-        raise ValueError(
-            f"{args.F} is not a structure factor amplitude in {args.mtz1}. Try again."
-        )
-    if not isinstance(mtz["SigF"].dtype, rs.StandardDeviationDtype):
-        raise ValueError(
-            f"{args.SigF} is not a structure factor amplitude in {args.mtz1}. Try again."
-        )
     if not isinstance(ref["Phi"].dtype, rs.PhaseDtype):
         raise ValueError(
             f"{args.Phi} is not a phases column in {args.mtz2}. Try again."

efxtools/utils/io.py

@@ -0,0 +1,51 @@
+import reciprocalspaceship as rs
+
+
+def subset_to_FSigF(mtzpath, data_col, sig_col, column_names_dict={}):
+    """
+    Utility function for reading MTZ and returning DataSet with F and SigF.
+
+    Parameters
+    ----------
+    mtzpath : str, filename
+        Path to MTZ file to read
+    data_col : str, column name
+        Column name for data column. If Intensity is specified, it will be
+        French-Wilson'd.
+    sig_col : str, column name
+        Column name for sigma column. Must select for a StandardDeviationDtype.
+    column_names_dict : dictionary
+        If particular column names are desired for the output, this can be specified
+        as a dictionary that includes `data_col` and `sig_col` as keys and what
+        values they should map to.
+
+    Returns
+    -------
+    rs.DataSet
+    """
+    mtz = rs.read_mtz(mtzpath)
+
+    # Check dtypes
+    if not isinstance(
+        mtz[data_col].dtype, (rs.StructureFactorAmplitudeDtype, rs.IntensityDtype)
+    ):
+        raise ValueError(
+            f"{data_col} must specify an intensity or |F| column in {mtzpath}"
+        )
+    if not isinstance(mtz[sig_col].dtype, rs.StandardDeviationDtype):
+        raise ValueError(
+            f"{sig_col} must specify a standard deviation column in {mtzpath}"
+        )
+
+    # Run French-Wilson if intensities are provided
+    if isinstance(mtz[data_col].dtype, rs.IntensityDtype):
+        scaled = rs.algorithms.scale_merged_intensities(
+            mtz, data_col, sig_col, mean_intensity_method="anisotropic"
+        )
+        mtz = scaled.loc[:, ["FW-F", "FW-SIGF"]]
+        mtz.rename(columns={"FW-F": data_col, "FW-SIGF": sig_col}, inplace=True)
+    else:
+        mtz = mtz.loc[:, [data_col, sig_col]]
+
+    mtz.rename(columns=column_names_dict, inplace=True)
+    return mtz

setup.py

@@ -14,6 +14,7 @@
             "efxtools.scaleit=efxtools.scaleit.scaleit:main",
             "efxtools.internal_diffmap=efxtools.diffmaps.internaldiffmap:main",
             "efxtools.ccsym=efxtools.stats.ccsym:main",
+            "efxtools.diffmap=efxtools.diffmaps.diffmap:main",
         ]
     },
 )