import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
import sys 
import argparse
import pandas as pd
import numpy as np
import os

### This for purely explorative purposes. The labelling is roughly done on geographic coordinates, but


def visualize_pcs(data, outfile, pc_a=1, pc_b=2):
    """
    Visualize PCA done with plink
    :param data: df, DataFrame with eigenvectors and geographic coordinates
    :param outfile: str, where to save Figure to
    :param pc_a: int, first PC to plot
    :param pc_b: int, second PC to plot
    """
    np.random.seed(42)
    colors = np.where(
        (data.Longitude >= -16.7) & (data.Longitude <= 17) & (data.Latitude >= 2.5) & (data.Latitude <= 20), 'green', 0)
    colors = np.where((data.Longitude > 17) & (data.Latitude >= -5.0) & (data.Latitude <= 20), 'blue', colors)
    colors = np.where((data.Latitude <= -5), 'orange', colors)
    colors = np.where((data.Latitude > 15) & (data.Longitude >= 35), 'pink', colors)
    colors = np.where((data.Latitude > 27) & (data.Longitude < 15), 'red', colors)
    color_dict = {'Middle East': 'pink',
                  "West Africa": "green",
                  "East Africa": "blue",
                  "South Africa": "orange",
                  "Europe": "red"}
    data['colors'] = colors
    fig, ax = plt.subplots()
    ax.scatter(data.loc[:, f'PC{pc_a}'], data.loc[:, f'PC{pc_b}'].values, c=data.colors, alpha=0.5)
    ax.set_xlabel(f'PC{pc_a}')
    ax.set_ylabel(f'PC{pc_b}')
    legend_elements = [Line2D([0], [0], color=col, marker='o', ls='', label=pop) for pop, col in color_dict.items()]
    ax.legend(handles=legend_elements, loc='upper center', ncol=4, bbox_to_anchor=(0.5, -.14))
    fig.savefig(outfile, bbox_inches='tight')
    plt.close()


def main(argv):
    parser = argparse.ArgumentParser(description='Visualize PCA performed with plink2')
    parser.add_argument('--eigenvec', help='.eigenvec file that was generated by plink. ')
    parser.add_argument('--coords', help='file containing geographic coordinates of sample locations.')
    parser.add_argument('--plot_dir', help='Directory where to save plots')
    args = parser.parse_args()
    if args.plot_dir.endswith('/'):
        plot_dir = args.plot_dir
    else:
        plot_dir = args.plot_dir + '/'
    if not os.path.isdir(plot_dir):
        os.makedirs(plot_dir)
    eigenvectors = pd.read_csv(args.eigenvec, sep='\t')
    eigenvectors.set_index('IID', inplace=True)
    coordinates = pd.read_csv(args.coords, sep='\t', header=None, names=['FID', "IID", "Latitude", "Longitude"])
    # join
    data = eigenvectors.join(coordinates.set_index('IID'))
    outfile = plot_dir + args.eigenvec.split('/')[-1].split('.eigenvec')[0] + "_pca.png"
    visualize_pcs(data, outfile)


if __name__ == '__main__':
    main(sys.argv[1:])

import numpy as np
import pkg_resources
from sklearn.impute import SimpleImputer
from pandas_plink import read_plink
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from feems.utils import prepare_graph_inputs
from feems import SpatialGraph, Viz
from feems.cross_validation import run_cv
import sys
import argparse
import pandas as pd

# change matplotlib fonts
plt.rcParams["font.family"] = "Arial"
plt.rcParams["font.sans-serif"] = "Arial"


def create_spatial_graph(data_path, prefix, coords_file):
    """
    Create SpatialGraph object
    :param data_path: str, data directory
    :param prefix: str, prefix for plink files in data directory
    :param coords_file: str, path to file with sampling coordinates
    :return: SpatialGraph
    """
    data_path_feems = pkg_resources.resource_filename("feems", "data/")
    # read the genotype data and mean impute missing data
    (bim, fam, G) = read_plink("{}/{}".format(data_path, prefix))
    imp = SimpleImputer(missing_values=np.nan, strategy="mean")
    genotypes = imp.fit_transform((np.array(G)).T)

    # setup graph
    coords = pd.read_csv(coords_file, sep='\t', header=None,
                         names=['fid_c', 'iid_c', 'latitude', 'longitude'])
    coords = fam.set_index('iid').join(coords.set_index('iid_c')).loc[:, ["longitude", "latitude"]].values
    grid_path = "{}/grid_100.shp".format(data_path_feems)  # path to discrete global grid

    # graph input files
    _, edges, grid, _ = prepare_graph_inputs(coord=coords,
                                             ggrid=grid_path,
                                             translated=False,
                                             buffer=0,)

    # # construct spatial graph object
    sp_graph = SpatialGraph(genotypes, coords, grid, edges, scale_snps=True)
    return sp_graph


def cross_validate(sp_graph):
    """
    Perform cross validation to find optimal lambda
    :param sp_graph: SpatialGraph
    :return: float, lambda
    """
    # reverse the order of lambdas and alphas for warmstart
    lamb_grid = np.geomspace(1e-6, 1e2, 20)[::-1]

    # run cross-validation
    cv_err = run_cv(sp_graph, lamb_grid, n_folds=sp_graph.n_observed_nodes, factr=1e10)

    # average over folds
    mean_cv_err = np.mean(cv_err, axis=0)

    # argmin of cv error
    lamb_cv = float(lamb_grid[np.argmin(mean_cv_err)])
    return lamb_cv


def fit_final_graph(sp_graph, lamb_cv, prefix, output_path):
    """
    Fit final SpatialGraph using optimal lambda inferred from cross-validation
    :param sp_graph: SpatialGraph
    :param lamb_cv: float, lambda
    :param output_path: str, directory where to save output figure named feems_plot.pdf
    """
    sp_graph.fit(lamb_cv)
    projection = ccrs.PlateCarree()
    # Plot the FEEMS result
    fig = plt.figure(figsize=(9, 16), dpi=600)
    ax = fig.add_subplot(1, 1, 1, projection=projection)
    ax.set_extent([-25, 60, -40, 50], crs=ccrs.PlateCarree())
    v = Viz(ax, sp_graph, projection=projection, edge_width=.5,
            edge_alpha=1, edge_zorder=100, sample_pt_size=20,
            obs_node_size=7.5, sample_pt_color="black",
            cbar_font_size=5, cbar_ticklabelsize=4)
    v.draw_map()
    v.draw_edges(use_weights=True)
    v.draw_obs_nodes(use_ids=False)
    v.draw_edge_colorbar()
    fig.savefig(output_path + prefix + '_feems_plot.pdf', bbox_inches='tight', dpi=600)


def main(argv):
    parser = argparse.ArgumentParser()
    parser.add_argument('-d', '--data_path', help='Data directory. default=./data/', default='./data/')
    parser.add_argument('-p', '--prefix', help='Prefix to plink files')
    parser.add_argument('-c', '--coords', help='Path to file coordinates of samples. '
                                               'default=data/population_geo_coords.tab',
                        default='data/population_geo_coords.tab')
    parser.add_argument('-o', '--output_path', help='Output directory. Will save Figure called feems_plot.pdf '
                                                    'in this direcotry. default=./results/', default='./results/')
    args = parser.parse_args()
    data_path = args.data_path
    prefix = args.prefix
    coords = args.coords
    output_path = args.output_path
    if not data_path.endswith('/'):
        data_path += '/'
    if not output_path.endswith('/'):
        output_path += '/'
    sp_graph = create_spatial_graph(data_path, prefix, coords)
    lamb_cv = cross_validate(sp_graph)
    fit_final_graph(sp_graph, lamb_cv, prefix, output_path)


if __name__ == '__main__':
    main(sys.argv[1:])

import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import glob
import sys
import argparse
import seaborn as sns
import numpy.ma as ma
from pykrige.ok import OrdinaryKriging
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.io.shapereader as shpreader
import shapely.geometry as sgeom
from shapely.ops import unary_union
from shapely.prepared import prep
from itertools import permutations


def is_land(x, y, land):
    """
    Function to determine if a given point is on land or not
    :param x: float, longitude
    :param y: float, latitude
    :param land:
    :return: boolean
    """
    return land.contains(sgeom.Point(x, y))


def build_mask(x, y, land):
    """
    Create a mask of points that are on land and should be plotted
    :param x: np.array, longitude values
    :param y: np.array, latitude values
    :param land:
    :return: np.array, boolean
    """
    # mesh
    xv, yv = np.meshgrid(x, y, indexing='xy')
    # initialize mask
    mask = np.zeros(xv.shape, dtype=bool)
    for i in range(x.shape[0]):
        for j in range(y.shape[0]):
            # check if point is land
            mask[j, i] = is_land(xv[j, i], yv[j, i], land)
    return mask


def sort_individuals(df, K):
    """
    Sort individuals by modal ancestry components and it's intensity
    :param df: pd.DataFrame, data frame containing ancestry porportions
    :param K: int, K used to run admixture
    :return: pd.DataFrame, sorted df
    """
    # get mean ancestry per component for each population
    mean_ancestry_pops = df.groupby('population').mean()
    # identify modal component for each population, i.e., the component with the highest ancestry proportion
    modal_clusters_pops = mean_ancestry_pops.idxmax(axis=1).sort_values()
    individual_orders = []
    # iterate over all k
    for i in range(K):
        # get all populations for which the current component is model
        pops = modal_clusters_pops[modal_clusters_pops == i].index.values.tolist()
        # sort the population by they mean ancestry proportions of the modal component
        pops_sorted = mean_ancestry_pops.loc[pops, i].sort_values(ascending=False).index.values
        inds = []
        # sort the individuals within a population by ancestry proportion
        for pop in pops_sorted:
            inds.extend(df.loc[df.population == pop, i].sort_values(ascending=False).index.values.tolist())
        individual_orders.extend(inds)

    df = df.loc[individual_orders, :]
    return df


def get_best_k_and_plot_best_k(qfiles, Kvals, fam, log_dir, colors, prefix):
    """
    Find best K across different ADMIXTURE runs using different K and plot admixture plot of best K separately
    :param qfiles: list, list of paths to aligned .Q files
    :param Kvals: list-like, K values that were used
    :param fam: df, pd.DataFrame of plink .fam file
    :param log_dir: str, path to directory with log files of ADMIXTURE runs
    :param colors: list, colors used for plotting
    :param prefix: str, prefix of output figure "_admixture_plot_best_K.pdf" will be appended
    :return: array-like, array-like; order of individuals, cross-validation error
    """
    min_k = min(Kvals)
    # get CV error of ADMIXTURE runs
    cv_error = np.zeros(len(Kvals))
    logfiles = glob.glob(log_dir + '*.log')
    # read CV error
    for log in logfiles:
        k = int(log.split('.')[-2])
        with open(log, 'r') as fp:
            for line in fp:
                if line.startswith("CV error"):
                    cv_error[k - min_k] = float(line.split(':')[1])
                    break
            fp.close()
    # get best k
    best_k = Kvals[np.argmin(cv_error)]
    df = pd.read_csv(qfiles[np.argmin(cv_error)], sep=' ', header=None, index_col=None)
    df = df.iloc[:, 6:]
    df.columns = [i for i in range(best_k)]
    df['population'] = fam.fid.values
    # sort individuals
    df_sorted = sort_individuals(df, best_k)
    # get order of idnviduals
    plot_best_k(df_sorted, best_k, min(cv_error), colors, prefix)
    return cv_error


def sort_individuals_specified_K(qfile, K, fam):
    """
    Sort individuals based on estimated ancestry proportions given a K
    :param qfile: str, path to aligned Q file
    :param K: int, K value
    :param fam: df, pd.DataFrame of plink .fam file
    :return: array-like,  order of individuals

    """
    df = pd.read_csv(qfile, sep=' ', header=None, index_col=None)
    df = df.iloc[:, 6:]
    df.columns = [i for i in range(K)]
    df['population'] = fam.fid.values
    # sort individuals
    df_sorted = sort_individuals(df, K)
    # get order of idnviduals
    individual_order = df_sorted.index.values
    return individual_order


def plot_krigin(qfile, K, path_to_coords, vmin, fam, colors, prefix):
    """
    Interpolate admixture proportion in space using Kriging method
    :param qfile: str, path to qfile corresponding to K
    :param K: int, K used for admixture run
    :param path_to_coords:
    :param vmin: float, minimum admixture proportion required for Kriging method
    :param fam: df, pd.DataFrame of plink .fam file
    :param colors: list, colors used for plotting
    :param prefix: str, prefix for figure
    """
    # The parameters are tailored to our specific case and may not give meaningful results in other cases
    if K > 4:
        raise AttributeError("Only K<= 4 is supported")
    land_shp_fname = shpreader.natural_earth(resolution='50m',
                                             category='physical', name='land')

    land_geom = unary_union(list(shpreader.Reader(land_shp_fname).geometries()))
    # create land object
    land = prep(land_geom)

    # read data
    coords = pd.read_csv(path_to_coords, sep='\t', header=None, names=['fid', 'lat', 'long'], usecols=[0, 2, 3])
    coords.drop_duplicates('fid', inplace=True)
    df = pd.read_csv(qfile, sep=' ', header=None, index_col=None)
    df = df.iloc[:, 6:]
    df.columns = [i for i in range(K)]
    df['population'] = fam.fid.values
    df = df.join(coords.set_index('fid'), on='population')

    # longitude and latidude coords of interest
    long_coords = np.arange(-25, 60, 0.15, dtype=float)
    lat_coords = np.arange(-40, 50, 0.15, dtype=float)
    # color maps
    cmaps = [
        sns.color_palette(f"blend:white,{colors[0]}", as_cmap=True),
        sns.color_palette(f"blend:white,{colors[1]}", as_cmap=True),
        sns.color_palette(f"blend:white,{colors[2]}", as_cmap=True),
        sns.color_palette(f"blend:white,{colors[3]}", as_cmap=True),
    ]

    projection = ccrs.PlateCarree()

    fig = plt.figure(figsize=(9, 16))
    ax = fig.add_subplot(1, 1, 1, projection=projection)
    ax.set_extent([-25, 60, -40, 50], crs=ccrs.PlateCarree())

    # run kriging for each K
    perms = list(set(permutations([i for i in range(K - 1)])))
    perms = [perm + (K-1,) for perm in perms]
    for perm in perms:
        for i in perm:
            cmap = cmaps[i]
            proportions = df.loc[:, i].values
            ok = OrdinaryKriging(df.long.values,
                                 df.lat.values,
                                 proportions, variogram_model='linear',
                                 coordinates_type='geographic', nlags=2)
            z, ss = ok.execute('grid', long_coords, lat_coords)
            # get mask
            mask = build_mask(long_coords, lat_coords, land)
            # normalize data
            data = (z.data - z.data.min()) / (z.data.max() - z.data.min())
            # update mask
            mask[data < vmin] = False
            masked_data = ma.masked_array(data, ~mask)
            # plot
            ax.imshow(np.flip(masked_data, 0), extent=[-25, 60, -40, 50], cmap=cmap,
                      interpolation='none', vmax=1, vmin=vmin, alpha=1 / len(perms))
    # formatting
    ax.add_feature(cfeature.COASTLINE, color='black')
    ax.add_feature(cfeature.BORDERS, linestyle=':', color='black')
    # plot populations
    ax.scatter(coords.long, coords.lat, color='black', s=8)
    ax.axis('off')
    fig.savefig(f"{prefix}_admixture_kriging_K{K}.pdf", bbox_inches='tight', dpi=600)


def plot_best_k(df, K, cv, colors, prefix):
    """
    Create horizontal bar plot of ADMIXTURE results
    :param df: pd.DataFrame, sorted DataFrame with ancestry proportions
    :param K: int, K used for admixture run
    :param colors: list, list of colors to use
    :param cv: float, CV error of ADMIXTURE run
    :param prefix: str, prefix for figure
    """
    fig, ax = plt.subplots(figsize=(2, 8))
    populations = df.population.values
    y_coords = [0]
    y_ticks = []
    y_labels = []
    prev_pop = populations[0]
    pop_coord_start = 0
    cumulative_padding = 0
    padding = 0.15
    for i, pop in enumerate(populations[1:], 1):
        if prev_pop != pop:
            # get ytick positions and labels
            y_ticks.append((pop_coord_start + i + cumulative_padding - 1) / 2)
            y_labels.append(prev_pop)
            prev_pop = pop
            # add white space between populations
            cumulative_padding += padding
            pop_coord_start = i + cumulative_padding
        # determine y coords --> add white space between populations
        y_coords.append(i + cumulative_padding)
    y_ticks.append((pop_coord_start + i + cumulative_padding) / 2)
    y_labels.append(prev_pop)
    prev_k = []
    # do plotting
    for i in range(K):
        if i > 0:
            ax.barh(y_coords, df.loc[:, i], left=df.loc[:, prev_k].sum(axis=1), height=1, color=colors[i])
        else:
            ax.barh(y_coords, df.loc[:, i], height=1, color=colors[i])
        prev_k.append(i)
    # formatting
    ax.set_yticks(y_ticks)
    ax.set_yticklabels([x.replace('Sabue', 'Chabu').replace('Xun', '!Xun').replace('_ERR_', "_Errachidia_")\
                        .replace('Juhoansi', "Ju|'Hoan").replace('SEBantu', 'BantuSAfrica')\
                        .replace('SWBantu', 'Herero').replace('Haihom', 'Hai||om').replace("_TIZ_", "_Tiznit_")\
                        .replace('_Sen_', '_Sened_').replace('_Chen_', '_Chenini_').replace('_', " ")\
                        .replace("BER", "Berber")
                        for x in y_labels], fontsize=5)
    ax.set_xticks([])
    ax.set_xticklabels([])
    ax.invert_yaxis()
    ax.set_xlabel(f"K={K}\nCV={round(cv, 4)}", fontsize=8)
    ax.set_ylim([len(populations) + cumulative_padding - 0.5, -0.5])
    ax.set_xlim([0, 1])
    fig.savefig(f"{prefix}_admixture_plot_best_K.pdf", bbox_inches='tight', dpi=600)


def plot_all_k(qfiles, Kvals, cv_error, fam, individual_order, colors, prefix):
    """
    Plot admixture proportions across all K
    :param qfiles: list, list of paths to aligned .Q files
    :param Kvals: list-like, K values that were used
    :param cv_error: array-like, cross-validation errors of individual runs
    :param fam: pd.DataFrame, df of plink .fam file
    :param individual_order: list, order of individudals in which to plot them to guarantee consistency
    :param colors: list, colors used for plotting
    :param prefix: str, prefix of output figure "_admixture_plots.pdf" will be appended
    """
    fig, ax = plt.subplots(1, len(qfiles), figsize=(8, 11))
    plt.subplots_adjust(wspace=0.15)
    min_k = min(Kvals)
    for qfile, K, cv in zip(qfiles, Kvals, cv_error):
        df = pd.read_csv(qfile, sep=' ', header=None, index_col=None)
        df = df.iloc[:, 6:]
        df.columns = [i for i in range(K)]
        df['population'] = fam.fid.values
        df_sorted = df.loc[individual_order]
        plot_admixture_proportions(df_sorted, K, ax[K - min_k], cv, colors, min_k)
    fig.savefig(prefix + "_admixture_plots.pdf", bbox_inches='tight', dpi=600)
    plt.close()


def plot_admixture_proportions(df, K, ax, cv, colors, min_k):
    """
    Helper to plot admixture proportions of for specific value of K
    :param df: pd.DataFrame, sorted DataFrame with ancestry proportions
    :param K: int, K used for admixture run
    :param ax: ax object to plot on
    :param cv: float, CV error of ADMIXTURE run
    :param colors: list, list of colors to use
    :param min_k: int, minimum k that was used in any run
    :return: ax
    """
    populations = df.population.values
    y_coords = [0]
    y_ticks = []
    y_labels = []
    prev_pop = populations[0]
    pop_coord_start = 0
    cumulative_padding = 0
    padding = 0.15
    for i, pop in enumerate(populations[1:], 1):
        if prev_pop != pop:
            # get ytick positions and labels
            y_ticks.append((pop_coord_start + i + cumulative_padding - 1) / 2)
            y_labels.append(prev_pop)
            prev_pop = pop
            # add white space between populations
            cumulative_padding += padding
            pop_coord_start = i + cumulative_padding
        # determine y coords --> add white space between populations
        y_coords.append(i + cumulative_padding)
    y_ticks.append((pop_coord_start + i + cumulative_padding) / 2)
    y_labels.append(prev_pop)
    prev_k = []
    # do plotting
    for i in range(K):
        if i > 0:
            ax.barh(y_coords, df.loc[:, i], left=df.loc[:, prev_k].sum(axis=1), height=1, color=colors[i])
        else:
            ax.barh(y_coords, df.loc[:, i], height=1, color=colors[i])
        prev_k.append(i)
    # formatting
    if K > min_k:
        ax.set_yticks([])
    else:
        ax.set_yticks(y_ticks)
        ax.set_yticklabels([x.replace('Sabue', 'Chabu').replace('Xun', '!Xun').replace('_ERR_', "_Errachidia_")\
                           .replace('Juhoansi', "Ju|'Hoan").replace('SEBantu', 'BantuSAfrica')\
                           .replace('SWBantu', 'Herero').replace('Haihom', 'Hai||om').replace("_TIZ_", "_Tiznit_")\
                           .replace('_Sen_', '_Sened_').replace('_Chen_', '_Chenini_').replace('_', " ")\
                           .replace("BER", "Berber")
                            for x in y_labels], fontsize=5)
    ax.set_xticks([])
    ax.set_xticklabels([])
    ax.invert_yaxis()
    ax.set_xlabel(f"K={K}\nCV={round(cv, 4)}", fontsize=5)
    ax.set_ylim([len(populations) + cumulative_padding - 0.5, -0.5])
    ax.set_xlim([0, 1])
    return ax


def main(argv):
    parser = argparse.ArgumentParser()
    parser.add_argument('-i', '--input_dir',
                        help='Input directory of "converted" files, i.e., .Q files aligned with CLUMPAK')
    parser.add_argument('-l', '--logfile_dir',
                        help='Directory containing log files from admixture runs to extract CV error.')
    parser.add_argument('-f', '--fam', help='Plink fam file to extract individual and family IDs')
    parser.add_argument('-c', '--coords', help='File with individual geographic coordinates')
    parser.add_argument('--vmin', type=float,
                        help='Only admixture proportions greater vmin are considered for Kriging method. default=0.5',
                        default=0.5)
    parser.add_argument('-k', '--K_kriging', type=int, help='K to use for plotting Kriging')
    parser.add_argument('-o', '--output_prefix', help='Output prefix')
    args = parser.parse_args()
    input_dir = args.input_dir
    log_dir = args.logfile_dir
    if not input_dir.endswith('/'):
        input_dir += '/'
    if not log_dir.endswith('/'):
        log_dir += '/'
    # get input q files
    aligned_q_files = glob.glob(input_dir + "*.converted")
    # read fam file
    fam = pd.read_csv(args.fam, sep='\t', names=['fid', 'iid'], usecols=[0, 1])
    # get k values
    Kvals = [int(f.split('.')[-3]) for f in aligned_q_files]
    # sort based on K
    aligned_q_files = np.array(aligned_q_files)[np.argsort(Kvals)]
    Kvals = np.sort(Kvals)
    max_k = max(Kvals)

    # initialize colors
    colors = ["purple", "orange", "blue", "green"]
    # pad colors
    additional_colors = []
    while max_k > len(colors) + len(additional_colors):
        additional_colors.append((np.random.random(), np.random.random(), np.random.random()))
        additional_colors = list(set(additional_colors))
    colors.extend(additional_colors)

    cv_error = get_best_k_and_plot_best_k(aligned_q_files, Kvals, fam, log_dir, colors, args.output_prefix)
    individual_order = sort_individuals_specified_K(aligned_q_files[np.where(Kvals == args.K_kriging)[0][0]],
                                                    args.K_kriging, fam)
    plot_all_k(aligned_q_files, Kvals, cv_error, fam, individual_order, colors, args.output_prefix)
    plot_krigin(aligned_q_files[np.where(Kvals == args.K_kriging)[0][0]], args.K_kriging, args.coords, args.vmin, fam,
                colors, args.output_prefix)


if __name__ == '__main__':
    main(sys.argv[1:])

shell:
    "wget -q {params.url}; unzip {params.zip}; rm {params.zip}; mv {output} {output}_tmp; "
    "unzip -d {output}_tmp {output}_tmp/26_03_2015_CLUMPAK.zip; mv {output}_tmp/26_03_2015_CLUMPAK/CLUMPAK .;"
    "rm -r {output}_tmp; chmod +x {output}/distruct/*; chmod +x {output}/CLUMPP/*; chmod +x {output}/mcl/*"