def METHOD_NAME():
    input_data = read_h5mu(par["input"])
    modality_data = input_data.mod[par["modality"]]
    var = modality_data.var
    layer = modality_data.X if not par['layer'] else modality_data.layers[par['layer']]
    if not issparse(layer):
        raise NotImplementedError("Expected layer to be in sparse format.")
    if isspmatrix_coo(layer):
        layer = csr_matrix(layer)
    layer.eliminate_zeros()
    # var statistics
    num_nonzero_obs = layer.getnnz(axis=0)
    obs_mean, _  = mean_variance_axis(layer, axis=0)
    pct_dropout = (1 - num_nonzero_obs / layer.shape[0]) * 100
    total_counts_obs = np.ravel(layer.sum(axis=0))
    # obs statistics
    num_nonzero_vars = layer.getnnz(axis=1)
    total_counts_var = np.ravel(layer.sum(axis=1))
    top_metrics = {}
    if par["top_n_vars"]:
        par["top_n_vars"] = sorted(par["top_n_vars"])
        distributions = get_top_from_csr_matrix(layer, par["top_n_vars"])
        top_metrics = {distribution_size: distribution * 100
                       for distribution_size, distribution 
                       in zip(par["top_n_vars"], distributions.T)}
    total_expr_qc = {}
    pct_expr_qc = {}
    if par["var_qc_metrics"]:
        for qc_metric in par["var_qc_metrics"]:
            if not qc_metric in var:
                raise ValueError(f"Value for --var_qc_metrics, {qc_metric} "
                                 f"not found in .var for modality {par['modality']}")
            qc_column = var[qc_metric]
            if qc_column.isna().any():
                if par["var_qc_metrics_fill_na_value"] is None:
                    raise ValueError(f"The .var column '{qc_metric}', selected by '--var_qc_metrics', contains NA values. "
                                     "It is ambiguous whether or not to include these values in the static calulation. "
                                     "You can explicitly map the NA values to 'False' or 'True using '--var_qc_metrics_fill_na_value'")
                else:
                    qc_column = qc_column.fillna(par['var_qc_metrics_fill_na_value'], inplace=False)
            qc_column = qc_column.values
            if not set(np.unique(qc_column)) == set((True, False)):
                raise ValueError(f"Column {qc_metric} in .var for modality {par['modality']} "
                                 f"must only contain boolean values")
            total_expr_qc[qc_metric] = np.ravel(layer[:, qc_column].sum(axis=1))
            pct_expr_qc[qc_metric] =  total_expr_qc[qc_metric] / total_counts_var * 100
    # Write all of the calculated statistics
    modality_data.var = modality_data.var.assign(
        **{"pct_dropout": pct_dropout,
           "num_nonzero_obs": num_nonzero_obs,
           "obs_mean": obs_mean,
           "total_counts": total_counts_obs})
    modality_data.obs = modality_data.obs.assign(
        **({"num_nonzero_vars": num_nonzero_vars,
            "total_counts": total_counts_var} | \
           {f"pct_{qc_metric}": col for qc_metric, col in pct_expr_qc.items()} | \
           {f"total_counts_{qc_metrix}": col for qc_metrix, col in total_expr_qc.items()}) | \
           {f"pct_of_counts_in_top_{n_top}_vars": col for n_top, col in top_metrics.items()})
    input_data.write(par["output"], compression=par["output_compression"])