Workflow Steps and Code Snippets

import bacchus.hic as bch
import bacchus.io as bcio
import bacchus.plot as bcp
import cooler
import copy
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import numpy as np
import os
from os.path import join
import pandas as pd
from scipy import ndimage
import seaborn as sns

# Import snakemake parameters.
annotation_file = snakemake.input.annotation
mat_wt_file = str(snakemake.input.cool_wt)
mat_rf_file = str(snakemake.input.cool_rf)
rna_wt_file = snakemake.input.rna_wt
rna_rf_file = snakemake.input.rna_rf
cov_wt_file = snakemake.input.cov_wt
cov_rf_file = snakemake.input.cov_rf
gc_file = snakemake.input.gc
epod_file = snakemake.input.EPODs
res = int(snakemake.params.res)
cmap = snakemake.params.cmap
outdir = str(snakemake.params.outdir)
width = snakemake.params.width
pos = snakemake.params.positions

# Create outdir if necessary.
os.makedirs(outdir, exist_ok=True)


def import_annotation_gff(annotation_file):
    """Function to create a table of the gene positions from the gff file."""
    annotation = pd.DataFrame(
        columns=["type", "start", "end", "strand", "name"]
    )
    with open(annotation_file, "r") as file:
        for line in file:
            # Header.
            if line.startswith("#"):
                continue
            # Stop at the fasta sequences.
            elif line.startswith(">"):
                break
            else:
                line = line.split("\t")
                if line[2] in ["gene", "tRNA", "rRNA"]:
                    if line[2] == "gene":
                        name = line[8].split("Name=")[-1].split(";")[0]
                        # Extract gene position.
                        annot = {
                            "type": line[2],
                            "start": int(line[3]),
                            "end": int(line[4]),
                            "strand": line[6],
                            "gene_name": name,
                        }
                        annotation = annotation.append(annot, ignore_index=True)
    return annotation


# Import files.
annotation = import_annotation_gff(annotation_file)
mat_wt = cooler.Cooler(f"{mat_wt_file}::/resolutions/{res}").matrix(
    balance=True, sparse=False
)[:]
mat_wt[np.isnan(mat_wt)] = 0
mat_rf = cooler.Cooler(f"{mat_rf_file}::/resolutions/{res}").matrix(
    balance=True, sparse=False
)[:]
mat_rf[np.isnan(mat_rf)] = 0
rna_wt, _ = bcio.extract_big_wig(rna_wt_file, binning=100)
rna_rf, _ = bcio.extract_big_wig(rna_rf_file, binning=100)
cov_wt, _ = bcio.extract_big_wig(cov_wt_file, binning=1000)
cov_rf, _ = bcio.extract_big_wig(cov_wt_file, binning=1000)
gc_content = pd.read_csv(gc_file, sep="\t", header=None).iloc[:, 1]
epod = pd.read_csv(epod_file, sep="\t", header=None)

# Make the rotation of the matrix to plot them
mat_wt_rot = copy.copy(mat_wt)
mat_wt_rot = bch.interpolate_white_lines(mat_wt_rot)
mat_wt_rot[np.isnan(mat_wt_rot)] = 0
mat_wt_rot = ndimage.rotate(mat_wt_rot, 45, reshape=True)
mat_rf_rot = copy.copy(mat_rf)
mat_rf_rot = bch.interpolate_white_lines(mat_rf_rot)
mat_rf_rot[np.isnan(mat_rf_rot)] = 0
mat_rf_rot = ndimage.rotate(mat_rf_rot, 45, reshape=True)


def plot_region(
    M1,
    M1_rot,
    M2,
    M2_rot,
    rna1,
    rna2,
    cov1,
    cov2,
    annotation,
    binning,
    width,
    zoom_ini,
    outfile,
    split=False,
):
    rna_binning = 100
    pal = sns.color_palette("Paired")

    # Defined values to specify the borders of the matrices and tables
    zoom = [zoom // binning for zoom in zoom_ini]
    width = width / 1000
    row1 = len(M1_rot) // 2 - int(np.sqrt(2) * width)
    row2 = len(M1_rot) // 2 + int(np.sqrt(2) * width)
    col1 = int(zoom[0] * np.sqrt(2))
    col2 = int(zoom[1] * np.sqrt(2))
    annotation_zoom = annotation[
        np.logical_and(
            annotation.start > zoom_ini[0], annotation.end < zoom_ini[1]
        )
    ]
    ymax = np.nanmax(
        np.concatenate(
            (
                rna1[zoom_ini[0] // rna_binning : zoom_ini[1] // rna_binning],
                rna2[zoom_ini[0] // rna_binning : zoom_ini[1] // rna_binning],
            )
        )
    )
    max_cov = np.nanmax(cov1) * 1.02
    min_cov = np.nanmin(cov1) * 0.98

    # Define panels depending on split or not.
    if split and ymax > 1500:
        a = 4
        fig, ax = plt.subplots(
            7,
            2,
            figsize=(16, 12),
            gridspec_kw={"height_ratios": [1, 30, 1, 6, 12, 3, 3]},
            sharex=False,
        )

        # Parameters to merge panel 2 and 3.
        d = 0.02
        D = 0.04
        for j in range(2):
            ax[a - 1, j].set_xlim(zoom_ini[0] // 1000, zoom_ini[1] // 1000)
            ax[a - 1, j].tick_params(axis="both", labelsize=15)
            ax[a - 1, j].spines["bottom"].set_visible(False)
            ax[a - 1, j].spines["top"].set_visible(False)
            ax[a - 1, j].spines["right"].set_visible(False)
            ax[a - 1, j].set_ylim(ymax - 525, ymax)
            ax[a - 1, j].get_xaxis().set_visible(False)

            # Add the small cut between the two panels
            kwargs = dict(
                transform=ax[a - 1, j].transAxes, color="k", clip_on=False
            )
            ax[a - 1, j].plot((-d, +d), (-D, +D), **kwargs)  # top-left diagonal
            kwargs.update(
                transform=ax[a, j].transAxes
            )  # switch to the bottom axes
            ax[a, j].plot(
                (-d, +d), (1 - d, 1 + d), **kwargs
            )  # bottom-left diagonal

        # Plot the RNA on the split panel
        ax[a - 1, 0].fill_between(
            np.arange(
                zoom_ini[0] // 1000, zoom_ini[1] // 1000, rna_binning / 1000
            ),
            rna1[zoom_ini[0] // rna_binning : zoom_ini[1] // rna_binning],
            color="black",
        )
        ax[a - 1, 1].fill_between(
            np.arange(
                zoom_ini[0] // 1000, zoom_ini[1] // 1000, rna_binning / 1000
            ),
            rna2[zoom_ini[0] // rna_binning : zoom_ini[1] // rna_binning],
            color="black",
        )

    else:
        a = 3
        fig, ax = plt.subplots(
            6,
            2,
            figsize=(16, 12),
            gridspec_kw={"height_ratios": [1, 30, 1, 15, 3, 3]},
            sharex=False,
        )

    # Plot expressed genes - blue are forward - red are reversed
    for j in range(2):
        ax[0, j].set_xlim(zoom_ini[0], zoom_ini[1])
        ax[0, j].get_xaxis().set_visible(False)
        ax[0, j].get_yaxis().set_visible(False)
        ax[2, j].set_xlim(zoom_ini[0], zoom_ini[1])
        ax[2, j].get_xaxis().set_visible(False)
        ax[2, j].get_yaxis().set_visible(False)
    pos = 1
    for i in range(len(annotation_zoom)):
        # Extract annotation information
        annot = annotation_zoom.iloc[
            i,
        ]
        strand = annot.strand
        start = annot.start // rna_binning
        end = annot.end // rna_binning
        name = annot.gene_name
        # Defined color depending on the strand
        if strand == "+":
            color = pal.as_hex()[1]
        else:
            color = pal.as_hex()[5]
        # Print it only if it transcribed (10% most transcribed genes)
        if np.mean(rna1[start:end]) >= 120.7628998139441:
            pos = pos * -1
            for j in range(2):
                ax[0, j].add_patch(
                    patches.Rectangle(
                        (start * 100, 0),
                        (end - start) * 100,
                        1,
                        edgecolor=color,
                        facecolor=color,
                        fill=True,
                    )
                )
                ax[0, j].text(
                    x=(start + ((end - start) / 2)) * 100,
                    y=pos + 0.5,
                    s=name,
                    rotation=90 * pos,
                    wrap=True,
                )
    color = pal.as_hex()[3]
    for i in epod.index:
        start = epod.loc[i, 3]
        end = epod.loc[i, 4]
        if (
            (start > zoom_ini[0])
            and (start < zoom_ini[1])
            or (end > zoom_ini[0])
            and (end < zoom_ini[1])
        ):
            start = max(start, zoom_ini[0]) // rna_binning
            end = min(end, zoom_ini[1]) // rna_binning
            for j in range(2):
                ax[2, j].add_patch(
                    patches.Rectangle(
                        (start * 100, 0),
                        (end - start) * 100,
                        1,
                        edgecolor=color,
                        facecolor=color,
                        fill=True,
                    )
                )
    #                 ax[2, j].text(
    #                     x=(start + ((end - start) / 2)) * 100,
    #                     y=pos + 0.5,
    #                     s="EPOD",
    #                     rotation=90 * pos,
    #                     wrap=True,
    #                 )

    # Plot the matrices
    ax[1, 0].set_ylabel("Genomic distance (kb)", fontsize=15)
    for j in range(2):
        ax[1, j].get_xaxis().set_visible(False)
        ax[1, j].tick_params(axis="both", labelsize=15)

    # Plot the matrice 1
    im = ax[1, 0].imshow(
        M1_rot[row1:row2, col1:col2],
        cmap=cmap,
        interpolation="none",
        vmin=0,
        vmax=np.nanpercentile(M1[zoom[0] : zoom[1], zoom[0] : zoom[1]], 97),
        extent=(col1 / np.sqrt(2), col2 / np.sqrt(2), width, -width),
    )

    # Plot the matrice 2
    ax[1, 1].imshow(
        M2_rot[row1:row2, col1:col2],
        cmap=cmap,
        interpolation="none",
        vmin=0,
        vmax=np.nanpercentile(M1[zoom[0] : zoom[1], zoom[0] : zoom[1]], 97),
        extent=(col1 / np.sqrt(2), col2 / np.sqrt(2), width, -width),
    )

    # RNAseq legends and plot at the bottom panel.
    ax[a, 0].set_ylabel("RNA count (CPM)", fontsize=15)
    for j in range(2):
        #         ax[a, j].set_xlabel("Genomic coordinates (kb)", size=15)
        ax[a, j].get_xaxis().set_visible(False)
        ax[a, j].tick_params(axis="both", labelsize=15)
        ax[a, j].set_xlim(zoom_ini[0] // 1000, zoom_ini[1] // 1000)
        ax[a, j].spines["top"].set_visible(False)
        ax[a, j].spines["right"].set_visible(False)
        if split:
            if ymax > 1500:
                ax[a, j].set_ylim(0, 1050)
            else:
                ax[a, j].set_ylim(0, 1500)
        else:
            ax[a, j].set_ylim(0, ymax)

    # Plot RNA at the bottom panel
    ax[a, 0].fill_between(
        np.arange(zoom_ini[0] // 1000, zoom_ini[1] // 1000, rna_binning / 1000),
        rna1[zoom_ini[0] // rna_binning : zoom_ini[1] // rna_binning],
        color="black",
    )
    ax[a, 1].fill_between(
        np.arange(zoom_ini[0] // 1000, zoom_ini[1] // 1000, rna_binning / 1000),
        rna2[zoom_ini[0] // rna_binning : zoom_ini[1] // rna_binning],
        color="black",
    )

    # GC content
    gc_binning = 100
    for j in range(2):
        ax[a + 1, j].plot(
            np.arange(
                zoom_ini[0] // 1000, zoom_ini[1] // 1000, gc_binning / 1000
            ),
            gc_content[
                (zoom_ini[0] + 250)
                // gc_binning : (zoom_ini[1] + 250)
                // gc_binning
            ],
            color="black",
        )
        ax[a + 1, j].set_xlim(zoom_ini[0] // 1000, zoom_ini[1] // 1000)
        ax[a + 1, j].tick_params(axis="both", labelsize=15)
        ax[a + 1, j].get_xaxis().set_visible(False)
    ax[a + 1, 0].set_ylabel("GC", fontsize=15)

    # Coverage
    cov_binning = 1000
    ax[a + 2, 0].plot(
        np.arange(zoom_ini[0] // 1000, zoom_ini[1] // 1000, cov_binning / 1000),
        cov1[
            (zoom_ini[0] + 250)
            // cov_binning : (zoom_ini[1] + 250)
            // cov_binning
        ],
        color="black",
    )
    ax[a + 2, 1].plot(
        np.arange(zoom_ini[0] // 1000, zoom_ini[1] // 1000, cov_binning / 1000),
        cov2[
            (zoom_ini[0] + 250)
            // cov_binning : (zoom_ini[1] + 250)
            // cov_binning
        ],
        color="black",
    )
    for j in range(2):
        ax[a + 2, j].set_xlim(zoom_ini[0] // 1000, zoom_ini[1] // 1000)
        ax[a + 2, j].set_ylim(min_cov, max_cov)
        ax[a + 2, j].set_xlabel("Genomic coordinates (kb)", size=15)
        ax[a + 2, j].tick_params(axis="both", labelsize=15)
    ax[a + 2, 0].set_ylabel("HiC\ncoverage\n(CPM)", fontsize=15)

    # Colorbar
    cbar = fig.colorbar(
        im, ax=ax.ravel().tolist(), shrink=0.3, anchor=(1.2, 0.75)
    )
    cbar.ax.tick_params(labelsize=15)

    # Save the fig and adjust margins
    if split and ymax > 1500:
        plt.subplots_adjust(wspace=0.2, hspace=0.1)
    else:
        plt.subplots_adjust(wspace=0.2, hspace=0.1)
    plt.savefig(outfile, bbox_inches="tight", dpi=200)


for split in [True, False]:
    if split:
        outfile = join(
            outdir,
            f"region_{pos[0]}_{pos[1]}_split.pdf",
        )
    else:
        outfile = join(
            outdir,
            f"region_{pos[0]}_{pos[1]}.pdf",
        )
    position = [int(p) * 1000 for p in pos]
    plot_region(
        mat_wt,
        mat_wt_rot,
        mat_rf,
        mat_rf_rot,
        rna_wt,
        rna_rf,
        cov_wt,
        cov_rf,
        annotation,
        res,
        width,
        position,
        outfile,
        split=split,
    )

shell:
    "snakemake --nolock --report {output}"

shell:
    "snakemake --nolock --report {output}"

shell:
    "snakemake --nolock --report {output}"

shell:
    "snakemake --nolock --report {output}"

#Snakemake interface
args = list()
if (exists('snakemake')) {
    args = list(estimatedRelFile = snakemake@input$estRel, 
                trueRelFile = snakemake@input$trueRel)

} else {
    warning("Not running with snakemake. Using test arguments!!!")
    args = list(estimatedRelFile = '~/tmp/tribes/TFEur/FF-EUR-15-30-2-mut_BiSnp_EurAF:0.01_LD_PH_GRM-allchr_IBD.csv', 
                trueRelFile = '~/tmp/tribes/TFEur/FF-EUR-15-30-2-mut_relations.txt')
}
print(args)

project_root <- here::here()
renv::load(project_root)

# import libraries
library(magrittr)
library(optparse)
library(SingleCellExperiment)

# source helper functions
source(file.path(project_root, "scripts", "utils", "integration-helpers.R"))

# Set up optparse options
option_list <- list(
  make_option(
    opt_str = c("-l", "--library_file"),
    type = "character",
    default = file.path(project_root, "sample-info", "hca-processed-libraries.tsv"),
    help = "path to file listing all libraries that are to be converted"
  ),
  make_option(
    opt_str = c("-g", "--grouping_var"),
    type = "character",
    default = "project_name",
    help = "Column name present in the library metadata file to use for grouping SCE objects
    and merging."
  ),
  make_option(
    opt_str = c("--groups_to_integrate"),
    type = "character",
    default = "All",
    help = "Groups present in `grouping_var` column of metadata file to create merged SCE objects for.
      Default is 'All' which produces a merged object for each group in the metadata file. 
      Specify groups by using a vector, e.g. c('group1','group2')"
  ),
  make_option(
    opt_str = c("--add_celltype"),
    action = "store_true",
    default = FALSE,
    help = "Boolean indicating whether or not celltype data should be added to the 
     individual SCE object prior to merging."
  ),
  make_option(
    opt_str = c("--celltype_info"),
    type = "character",
    default = file.path(project_root, "sample-info", "hca-celltype-info.tsv"),
    help = "Path to file containing cell type information for each SCE object. 
      Must contain columns named `library_biomaterial_id`, `celltype`, and `barcode`."
  ),
  make_option(
    opt_str = c("--batch_column"),
    type = "character",
    default = "batch",
    help = "The name of the column in colData that indicates the batches for each cell,
      typically this corresponds to the library id. Default is 'batch'."
  ),
  make_option(
    opt_str = c("--random_merge"),
    default = FALSE,
    action = "store_true",
    help = "Used to indicate whether or not to merge SCE objects in a random order. Default is FALSE."
  ),
  make_option(
    opt_str = c("--subset_hvg"),
    default = FALSE,
    action = "store_true",
    help = "Indicates whether or not to subset the merged SCE object by highly variable genes.
      If --subset_hvg is used, the merged SCE object will only contain genes
      identified as highly variable genes."
  ),
  make_option(
    opt_str = c("--pca_use_all_genes"),
    default = FALSE,
    action = "store_true",
    help = "Indicates whether or not to use the all genes as input to performing
      principal component analysis. Otherwise only highly variable genes are used
      as input."
  ),
  make_option(
    opt_str = c("-n", "--num_hvg"),
    dest = "num_genes",
    type = "integer",
    default = 5000,
    help = "Number of highly variable genes to select."
  ),
  make_option(
    opt_str = c("--merged_sce_dir"),
    type = "character",
    default = file.path(project_root, "results", "human_cell_atlas", "merged_sce"),
    help = "path to folder where all merged SCE objects files will be saved as RDS files"
  ),
  make_option(
    opt_str = c("--seed"),
    type = "integer",
    default = NULL,
    help = "random seed to set prior to merging"
  )
)

# Setup ------------------------------------------------------------------------

# Parse options
opt <- parse_args(OptionParser(option_list = option_list))

set.seed(opt$seed)

# check that num genes provided is an integer
if(!is.integer(opt$num_genes)){
  stop("--num_hvg must be an integer.")
}

# checks that provided metadata files exist
if(!file.exists(opt$library_file)){
  stop("--library_file provided does not exist.")
}

# read in library metadata and grab unfiltered sce file paths
library_metadata_df <- readr::read_tsv(opt$library_file)

# check that cell type file exists if using add_celltype option 
if(opt$add_celltype){
  if(!file.exists(opt$celltype_info)){
    stop("--celltype_info file provided does not exist.")
  }
  celltype_info_df <- readr::read_tsv(opt$celltype_info)
}

# check that grouping variable is present
if(!opt$grouping_var %in% colnames(library_metadata_df)){
  stop("Must provide a grouping_var that is a column in the library metadata file.")
}

# define groups to integrate
groups <- library_metadata_df %>%
  dplyr::pull(opt$grouping_var) %>%
  unique()

# check that groups to integrate isn't set to All 
if(length(opt$groups_to_integrate) == 1 && (opt$groups_to_integrate == "All")){
  groups_to_integrate <- groups
} else {
  # if All is not used then unlist groups, using space, needed to parse a list from snakemake
  groups_to_integrate <- unlist(stringr::str_split(opt$groups_to_integrate, " "))

  # check that specified groups are present in grouping_var column 
  if(!any(groups_to_integrate %in% groups)){
    stop("Provided `--groups_to_integrate` must also be present in the `--grouping_var` column of 
         the library metadata file.")
  }
}

# subset metadata file to only contain groups to integrate 
library_metadata_df <- library_metadata_df %>%
  dplyr::filter(.data[[opt$grouping_var]] %in% groups_to_integrate)

# grab library ids 
library_ids <- library_metadata_df %>%
  dplyr::pull(library_biomaterial_id)

# setup output directory
create_dir(opt$merged_sce_dir)

# Identify SCE files -----------------------------------------------------------

# grab all unique directories corresponding to projects considered for integration
input_sce_dirs <- library_metadata_df %>%
  dplyr::pull(integration_input_dir) %>%
  unique()

# find SCE files that match library ID, and throw an error if any are missing.
sce_files <- find_input_sce_files(library_ids, input_sce_dirs)

# Merge by group ---------------------------------------------------------------

# get the library IDs from the SCE file names so that we can name the SCEs in the correct order
library_ids_sce_order <- stringr::str_extract(sce_files, 
                                              pattern = paste(library_ids, collapse = "|"))

sce_file_df <- data.frame(sce_files = sce_files,
                          library_id= library_ids_sce_order) %>%
  dplyr::left_join(library_metadata_df, by = c("library_id" = "library_biomaterial_id")) %>%
  dplyr::select(library_id, opt$grouping_var, sce_files)

# group dataframe by the grouping variable
grouped_sce_file_df <- split(sce_file_df, sce_file_df[,opt$grouping_var])

# create a list of SCE lists that is named by the grouping variable with
# each individual inner SCE list named by the library IDs
create_grouped_sce_list <- function(sce_info_dataframe,
                                    celltype_info_df = NULL,
                                    add_celltype){

  library_sce_list = list()
  for (library_idx in 1:length(sce_info_dataframe$library_id)){

    # read sce list for each library
    sce <- readr::read_rds(sce_info_dataframe$sce_files[library_idx])
    library_name <- sce_info_dataframe$library_id[library_idx]

    if(add_celltype){
      # if celltype info is provided add to sce object 
      if(!is.null(celltype_info_df)){
        # check that library has cell type information
        if(library_name %in% unique(celltype_info_df$library_biomaterial_id)){

          # filter celltype info to only have info for specified library 
          filtered_celltype_info <- celltype_info_df %>%
            dplyr::filter(library_biomaterial_id == library_name) %>%
            dplyr::select(barcode, celltype)

          # add celltype info 
          sce <- add_celltype_info(sce_object = sce, 
                                   celltype_info_df = filtered_celltype_info) 

          # add flag indicating that cell type information is available 
          metadata(sce)$celltype_info_available <- TRUE 
        }

        # only add celltype column/metadata if add celltype is yes, but no celltype data is available 
      } else {
        # note that no cell type information is available
        colData(sce)$celltype <- NA
        metadata(sce)$celltype_info_available <- FALSE
      } 
    }

    # create a list for each group named by the library IDs
    library_sce_list[[library_name]] <- sce

  }

  return(library_sce_list)

}

# create grouped sce list with/without celltype addition
if(opt$add_celltype){
  grouped_sce_list <- grouped_sce_file_df %>%
    purrr::map(create_grouped_sce_list, celltype_info_df, add_celltype = opt$add_celltype) 
} else {
  grouped_sce_list <- grouped_sce_file_df %>%
    purrr::map(create_grouped_sce_list, add_celltype = opt$add_celltype)
}

# create a list of merged SCE objects by group
#  In this default usage, a batch column named `batch` will get created
merged_sce_list <- grouped_sce_list %>%
  purrr::map(combine_sce_objects, 
             preserve_rowdata_columns = c("Gene", "gene_names", "ensembl_ids"),
             random_merge = opt$random_merge)

# HVG and dim reduction --------------------------------------------------------

# apply HVG calculation to list of merged SCEs
# object will only be subset to HVG if subset_hvg is true
merged_sce_list <- merged_sce_list %>%
  purrr::map(~ set_var_genes(.x,
                             num_genes = opt$num_genes,
                             subset_hvg = opt$subset_hvg,
                             batch_column = opt$batch_column))

# add PCA and UMAP
# if --pca_use_all_genes is used, use all genes, otherwise only HVG are used
if(opt$pca_use_all_genes){
  merged_sce_list <- merged_sce_list %>%
    purrr::map( ~ perform_dim_reduction(.x,
                                        var_genes = rownames(.x),
                                        pca_type = "multi"))
} else {
  merged_sce_list <- merged_sce_list %>%
    purrr::map( ~ perform_dim_reduction(.x,
                                        var_genes = metadata(.x)$variable_genes,
                                        pca_type = "multi"))
}


# Write RDS --------------------------------------------------------------------

# create paths to merged SCE files
# named with the name of the sce list which corresponds to the grouping variable, not library ID
merged_sce_files <- file.path(opt$merged_sce_dir,
                              paste0(names(merged_sce_list),
                                     "_merged_sce.rds"))

# export files
purrr::walk2(merged_sce_list, merged_sce_files, readr::write_rds)

shell:  """
    mkdir -p '{params.external_software_dir}'
    cd '{params.external_software_dir}'

    rm -rf Flye
    git clone https://github.com/sebschmi/Flye
    cd Flye
    git checkout 38921327d6c5e57a59e71a7181995f2f0c04be75

    mv bin/flye bin/flye.disabled # rename such that snakemake does not delete it
    """

shell:  """
    cd '{params.flye_directory}'

    export CXX=x86_64-conda-linux-gnu-g++
    export CC=x86_64-conda-linux-gnu-gcc
    # export INCLUDES=-I/usr/include/ # Somehow this is not seen by minimap's Makefile, so we had to change it in our custom version of Flye
    # The following also doesn't seem to work when building minimap, so again we had to modify minimap's Makefile
    # export LD_LIBRARY_PATH=$CONDA_PREFIX/lib:${{LD_LIBRARY_PATH:=''}} # Redirect library path to include conda libraries
    # make # This does not create the python script anymore

    /usr/bin/env python3 setup.py install

    mv bin/flye.disabled bin/flye # was renamed such that snakemake does not delete it
    """

shell:  """
    mkdir -p '{params.external_software_dir}'
    cd '{params.external_software_dir}'

    rm -rf rust-mdbg
    git clone https://github.com/sebschmi/rust-mdbg
    cd rust-mdbg
    git checkout 4ff0122a8c63210820ba0341fa7365d6ac216612

    cargo fetch

    # rename such that snakemake does not delete them
    mv utils/magic_simplify utils/magic_simplify.original
    mv utils/multik utils/multik.original
    """