sos run GWAS_QC.ipynb -h

[global]
# the output directory for generated files
parameter: cwd = path("output")
# A string to identify your analysis run
parameter: name = ""
# PLINK binary files (either BED/BIM/FAM or PGEN/PVAR/PSAM format)
parameter: genoFile = paths
# The path to the file that contains the list of samples to remove (format FID, IID)
parameter: remove_samples = path('.')
# The path to the file that contains the list of samples to keep (format FID, IID)
parameter: keep_samples = path('.')
# The path to the file that contains the list of variants to keep
parameter: keep_variants = path('.')
# The path to the file that contains the list of variants to exclude
parameter: exclude_variants = path('.')
# Kinship coefficient threshold for related individuals
# (e.g first degree above 0.25, second degree above 0.125, third degree above 0.0625)
parameter: kinship = 0.0625
# For cluster jobs, number commands to run per job
parameter: job_size = 1
# Wall clock time expected
parameter: walltime = "5h"
# Memory expected
parameter: mem = "16G"
# Number of threads
parameter: numThreads = 20
# Software container option
parameter: container = ""
parameter: entrypoint= ""
# use this function to edit memory string for PLINK input
from sos.utils import expand_size
cwd = path(f"{cwd:a}")

# Determine if the file is in PLINK1 (BED/BIM/FAM) or PLINK2 (PGEN/PVAR/PSAM) format
def determine_plink_format(file_path):
    """
    Determine the PLINK file format based on file extensions and companion files.
    
    Args:
        file_path (str or Path): Path to the input file
    
    Returns:
        str: 'plink1' or 'plink2'
    """
    # Convert to string if it's a Path object
    file_path = str(file_path)
    
    # Check direct file extensions
    if file_path.endswith('.bed'):
        return 'plink1'
    elif file_path.endswith('.pgen'):
        return 'plink2'
    
    # If the file doesn't have a standard extension, try to infer format
    try:
        # Remove the file extension if present
        base_path = file_path.rsplit('.', 1)[0] if '.' in file_path else file_path
        
        # Check for PLINK1 companion files
        plink1_companion_files = [
            f"{base_path}.bim",
            f"{base_path}.fam"
        ]
        
        # Check for PLINK2 companion files
        plink2_companion_files = [
            f"{base_path}.pvar",
            f"{base_path}.psam"
        ]
        
        # Check PLINK1 format
        if all(os.path.exists(f) for f in plink1_companion_files):
            return 'plink1'
        
        # Check PLINK2 format
        if all(os.path.exists(f) for f in plink2_companion_files):
            return 'plink2'
    
    except Exception as e:
        print(f"Error determining PLINK format: {e}")
    
    # Default to PLINK1 if can't determine
    return 'plink1'


# Get the appropriate PLINK command based on the input file format
def get_plink_command_prefix(file_path):
    format_type = determine_plink_format(file_path)
    if format_type == 'plink1':
        return "--bfile"
    else:  # plink2
        return "--pfile"
        
# Generate the appropriate file extension based on the requested format
def get_output_extension(output_format, is_prune=False):
    if output_format == 'plink1':
        return '.bed' if not is_prune else '.prune.bed'
    else:  # plink2
        return '.pgen' if not is_prune else '.prune.pgen'
        
# Choose the make-bed or make-pgen command based on desired output format
def get_make_command(output_format):
    if output_format == 'plink1':
        return '--make-bed'
    else:  # plink2
        return '--make-pgen'

# Inference of relationships in the sample to identify closely related individuals
[king_1]
# PLINK binary file
parameter: kin_maf = 0.01
input: genoFile
output: f'{cwd}/{_input:bn}{("."+name) if name else ""}.kin0'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'

plink_command = get_plink_command_prefix(genoFile)

bash: container=container, expand= "${ }", stderr = f'{_output}.stderr', stdout = f'{_output}.stdout', entrypoint=entrypoint
    plink2 \
      ${plink_command} ${_input:n} \
      --make-king-table \
      --king-table-filter ${kinship} \
      ${('--keep %s' % keep_samples) if keep_samples.is_file() else ""} \
      ${('--remove %s' % remove_samples) if remove_samples.is_file() else ""} \
      --min-af ${kin_maf} \
      --max-af ${1-kin_maf} \
      --out ${_output:n} \
      --threads ${numThreads} \
      --memory ${int(expand_size(mem) * 0.9)/1e6} 
    
bash: expand= "${ }", stderr = f'{_output}.stderr', stdout = f'{_output}.stdout', container = container, entrypoint=entrypoint
    i="${_output}"
    output_size=$(ls -lh $i | cut -f 5 -d ' ')
    output_rows=$(cat $i | wc -l | cut -f 1 -d ' ')
    output_column=$(cat $i | head -1 | wc -w)
    output_preview=$(cat $i | grep -v "##" | head | cut -f 1,2,3,4,5,6)
    
    printf "output_info: %s\noutput_size: %s\noutput_rows: %s\noutput_column: %s\noutput_preview:\n%s\n" \
        "$i" "$output_size" "$output_rows" "$output_column" "$output_preview"

# Select a list of unrelated individual with an attempt to maximize the unrelated individuals selected from the data 
[king_2: shared = "related_id" ]
related_id = [x.strip() for x in open(_input).readlines() if not x.startswith("#")]
output: f'{_input:n}.related_id'
with open(_output, 'a'):
    pass
done_if(len(related_id) == 0, msg = f"No related individuals detected from {_input}.")
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'
R:  container=container, expand= "${ }", stderr = f'{_output}.stderr', stdout = f'{_output}.stdout', entrypoint=entrypoint

    # By Rui Dong and Derek Lamb, Columbia Neurology
    # This function enhances the standard plinkQC::relatednessFilter approach through:
    #
    # 1. Graph-based preprocessing:
    #    - For large relatedness networks, breaks down dense clusters of related individuals
    #    - Removes highly connected nodes (individuals with many relatives) first
    #    - Prevents memory issues when processing large components (>20 individuals)
    #    - Reduces the input size for the more intensive plinkQC filtering
    #
    # 2. Additional checks:
    #    - Verifies that all related pairs are properly handled
    #    - Runs up to 20 iterations to ensure complete filtering
    #    - Catches edge cases that might be missed in a single pass
    #
    # This implementation is particularly helpful for datasets with large family 
    # clusters or complex relatedness structures where the standard approach
    # might struggle.    

    filter_relatedness <- function(
      relatedness,                                  # Path to the relatedness file (required)
      relatednessTh = 0.0625,                       # Threshold for relatedness
      analysis_type = "maximize_unrelated",         # Type of analysis: "maximize_unrelated" or "maximize_cases"
      output_prefix = "output_prefix",              # Base name for output file
      otherCriterion = NULL,                        # Optional additional criterion for filtering
      otherCriterionTh = NULL,                      # Threshold for the additional criterion
      otherCriterionThDirection = "ge",             # Direction of the threshold
      relatednessIID1 = "IID1",                     # Column name for first individual ID
      relatednessIID2 = "IID2",                     # Column name for second individual ID
      relatednessFID1 = NULL,                       # Column name for first family ID
      relatednessFID2 = NULL,                       # Column name for second family ID
      relatednessRelatedness = "PI_HAT",            # Column name for relatedness value
      pheno_file = NULL,                            # File for phenotype information
      pheno_col = "pheno",                          # Column name for phenotype
      otherCriterionIID = "IID",                    # Column name for additional criterion ID
      otherCriterionMeasure = NULL,                 # Column name for additional criterion measure
      verbose = FALSE                               # Whether to print verbose output
    ) {
      # Load required packages
      suppressMessages({
        library(tidyverse)
        library(igraph)
        library(plinkQC)
        library(data.table)
      })
      
      # Validate inputs
      if (is.null(relatedness)) {
        stop("Relatedness file path is required")
      }
      
      if (analysis_type == "maximize_cases" && is.null(pheno_file)) {
        stop("Must provide phenotype file when analysis_type is 'maximize_cases'")
      }
      
      # Read the relatedness file
      full_kin0 <- as.data.frame(relatedness)
      
      # Initial pruning of highly related individuals
      working_graph <- full_kin0 |>
        filter(!!sym(relatednessRelatedness) >= relatednessTh) |>
        select(!!sym(relatednessIID1), !!sym(relatednessIID2)) |>
        graph_from_data_frame(directed = FALSE)
      
      working_comp <- components(working_graph)
      
      # Set parameters
      graph_size_th = 20
      reduce_fraction = 0.05
      high_related_indiv <- c()
      
      # Reduce components to max 20
      while (max(working_comp$csize > graph_size_th)) {
        
        if (verbose) {
          print(paste0("The largest component graph has ", max(working_comp$csize), 
                       " individuals. Removing the most connected ", 
                       scales::percent(reduce_fraction), " of them."))
        }
        
        # Prepare for loop
        large_comp_ids <- which(working_comp$csize > graph_size_th)
        nodes_to_remove <- c()
        
        for (comp_id in large_comp_ids) {
          # Get nodes/degrees for component
          comp_nodes <- V(working_graph)[working_comp$membership == comp_id]
          comp_degrees <- degree(working_graph, v = comp_nodes)
          
          # Remove highly related individuals
          num_to_remove <- ceiling(length(comp_nodes) * reduce_fraction)
          high_degree_nodes <- names(sort(comp_degrees, decreasing = TRUE))[1:num_to_remove]
          
          # Append to removal list
          nodes_to_remove <- c(nodes_to_remove, high_degree_nodes)
        }
        
        # Update for next loop
        high_related_indiv <- c(high_related_indiv, nodes_to_remove)
        working_graph <- delete_vertices(working_graph, nodes_to_remove)
        working_comp <- components(working_graph)
      }
      
      kin0 <- full_kin0 |>
        filter(!(!!sym(relatednessIID1) %in% high_related_indiv) & 
               !(!!sym(relatednessIID2) %in% high_related_indiv))
      
      # Process based on analysis type
      if (analysis_type == "maximize_unrelated") {
        # Run the plinkQC relatedness Filter
        rel <- plinkQC::relatednessFilter(
          relatedness = kin0,
          otherCriterion = otherCriterion,
          relatednessTh = relatednessTh,
          relatednessIID1 = relatednessIID1,
          relatednessIID2 = relatednessIID2,
          otherCriterionTh = otherCriterionTh,
          otherCriterionThDirection = otherCriterionThDirection,
          relatednessFID1 = relatednessFID1,
          relatednessFID2 = relatednessFID2,
          relatednessRelatedness = relatednessRelatedness,
          otherCriterionIID = otherCriterionIID,
          otherCriterionMeasure = otherCriterionMeasure,
          verbose = verbose
        )$failIDs
        
        all_exclude <- rel$IID
        
      } else if (analysis_type == "maximize_cases") {
        # Load phenotype information
        related_individuals <- unique(c(kin0 |> pull(relatednessIID1), kin0 |> pull(relatednessIID2)))
        
        df_pheno <- fread(pheno_file) |>
          drop_na(!!sym(pheno_col)) |>
          filter(IID %in% related_individuals)
        
        # Match kinship and phenotype information
        related_pheno_individuals <- df_pheno |> 
          pull(IID)
        
        related_cases <- df_pheno |> 
          filter(!!sym(pheno_col) == 1) |> 
          pull(IID)
        
        related_controls <- df_pheno |> 
          filter(!!sym(pheno_col) == 0) |> 
          pull(IID)
        
        kin0 <- kin0 |>
          filter(!!sym(relatednessIID1) %in% related_pheno_individuals & 
                 !!sym(relatednessIID2) %in% related_pheno_individuals)
        
        # Optimize cases
        df_case_kin <- kin0 |>
          filter(!!sym(relatednessIID1) %in% related_cases & 
                 !!sym(relatednessIID2) %in% related_cases)
        
        rel_cases <- plinkQC::relatednessFilter(
          relatedness = df_case_kin,
          otherCriterion = otherCriterion,
          relatednessTh = relatednessTh,
          relatednessIID1 = relatednessIID1,
          relatednessIID2 = relatednessIID2,
          otherCriterionTh = otherCriterionTh,
          otherCriterionThDirection = otherCriterionThDirection,
          relatednessFID1 = relatednessFID1,
          relatednessFID2 = relatednessFID2,
          relatednessRelatedness = relatednessRelatedness,
          otherCriterionIID = otherCriterionIID,
          otherCriterionMeasure = otherCriterionMeasure,
          verbose = verbose
        )$failIDs
        
        # Remove controls related to cases
        related_cases_keep <- setdiff(related_cases, rel_cases$IID)
        
        related_controls_exclude <- c()
        
        for (i in 1:nrow(kin0)) {
          iid1 <- kin0 |> pull(!!sym(relatednessIID1)) |> nth(i) 
          iid2 <- kin0 |> pull(!!sym(relatednessIID2)) |> nth(i)
          
          if (iid1 %in% related_cases_keep & iid2 %in% related_controls) {
            related_controls_exclude <- c(related_controls_exclude, iid2)
            
          } else if(iid2 %in% related_cases_keep & iid1 %in% related_controls) {
            related_controls_exclude <- c(related_controls_exclude, iid1)
            
          }
        }
        
        # Subset controls to only unrelated controls
        related_controls_keep <- setdiff(related_controls, related_controls_exclude)
        
        df_control_kin <- kin0 |>
          filter(!!sym(relatednessIID1) %in% related_controls_keep & 
                 !!sym(relatednessIID2) %in% related_controls_keep)
        
        rel_controls <- plinkQC::relatednessFilter(
          relatedness = df_control_kin,
          otherCriterion = otherCriterion,
          relatednessTh = relatednessTh,
          relatednessIID1 = relatednessIID1,
          relatednessIID2 = relatednessIID2,
          otherCriterionTh = otherCriterionTh,
          otherCriterionThDirection = otherCriterionThDirection,
          relatednessFID1 = relatednessFID1,
          relatednessFID2 = relatednessFID2,
          relatednessRelatedness = relatednessRelatedness,
          otherCriterionIID = otherCriterionIID,
          otherCriterionMeasure = otherCriterionMeasure,
          verbose = verbose
        )$failIDs
        
        # Handle leftover individuals
        all_exclude <- c(rel_cases$IID, related_controls_exclude, rel_controls$IID)
        
      } else {
        # Directly filter on kinship threshold
        rel <- kin0 %>% dplyr::filter(!!sym(relatednessRelatedness) >= relatednessTh)
        
        all_exclude <- unique(c(rel |> pull(relatednessIID1), rel |> pull(relatednessIID2)))
      }
      
      # Check for leftovers
      df_related <- kin0 |>
        filter(!(!!sym(relatednessIID1) %in% all_exclude) & 
               !(!!sym(relatednessIID2) %in% all_exclude)) |>
        filter(!!sym(relatednessRelatedness) > relatednessTh)
      print("df_related")
      print(df_related)
      if (nrow(df_related) > 0) {
        if (verbose) {
          print(paste0("Warning: After first plinkQC pass, there are still ", 
                      nrow(df_related), " related subjects in dataset."))
        }
        
        # Iterate until all subjects are removed
        iter <- 0
        while (nrow(df_related) > 0 & iter < 20) {
          # Rerun plinkQC
          additional_exclude <- plinkQC::relatednessFilter(
            relatedness = df_related,
            otherCriterion = otherCriterion,
            relatednessTh = relatednessTh,
            relatednessIID1 = relatednessIID1,
            relatednessIID2 = relatednessIID2,
            otherCriterionTh = otherCriterionTh,
            otherCriterionThDirection = otherCriterionThDirection,
            relatednessFID1 = relatednessFID1,
            relatednessFID2 = relatednessFID2,
            relatednessRelatedness = relatednessRelatedness,
            otherCriterionIID = otherCriterionIID,
            otherCriterionMeasure = otherCriterionMeasure,
            verbose = verbose
          )$failIDs
          
          # Add to exclude list
          all_exclude <- c(all_exclude, additional_exclude$IID)
          
          # Update df_related and iter
          df_related <- kin0 |>
            filter(!(!!sym(relatednessIID1) %in% all_exclude) & 
                   !(!!sym(relatednessIID2) %in% all_exclude)) |>
            filter(!!sym(relatednessRelatedness) > relatednessTh)
          iter <- iter + 1
        }
        
        # Print status
        if (verbose) {
          if (nrow(df_related) == 0) {
            print(paste0("All related subjects successfully removed after ", iter, " iterations."))
          } else if (nrow(df_related) > 0) {
            stop(paste0("Error: After 20 plinkQC passes, there are still ", 
                       nrow(df_related), " subjects in dataset."))
          }
        }
      }
      
      # Add in high-related individuals
      all_exclude <- c(all_exclude, high_related_indiv)
      
      # List of all excluded subjects
      dat <- data.frame(IID = all_exclude, FID = as.character(all_exclude)) # Because IID and FID are the same
      
      # Output the related ids
      output_related <- paste0(output_prefix, "_kinship_cutoff_", relatednessTh, ".related_id")
      write.table(dat, output_related, quote = FALSE, row.names = FALSE, col.names = FALSE)
      
      if (verbose) {
        cat("\n[INFO] There are", nrow(dat), "excluded individuals identified using a kinship threshold of", relatednessTh, "\n")
        cat("[OUTPUT] Excluded individual IDs have been written to the file:\n")
        cat("      ", normalizePath(output_related), "\n")
      }
      
      # Return the excluded IDs
      return(dat)
    }
    
  
    # main code
    suppressMessages({
      library(tidyverse)
      library(data.table)
    })
    
    # Read input and set column names to match your original code
    kin0 <- read.table(${_input:r}, header=FALSE, stringsAsFactors=FALSE)
    colnames(kin0) <- c("FID1","ID1","FID2","ID2","NSNP","HETHET","IBS0","KINSHIP")
    
    # Always use the filter_relatedness function with maximize_unrelated
    result <- filter_relatedness(
      relatedness = kin0,  
      relatednessTh = ${kinship},
      analysis_type = "maximize_unrelated",
      relatednessIID1 = "ID1",
      relatednessIID2 = "ID2",
      relatednessFID1 = "FID1", 
      relatednessFID2 = "FID2",
      relatednessRelatedness = "KINSHIP",
      verbose = FALSE
    )
    
    # Format output to match your original code
    tmp1 <- kin0[,1:2]
    tmp2 <- kin0[,3:4]
    colnames(tmp1) = colnames(tmp2) = c("FID", "ID")
    lookup <- dplyr::distinct(rbind(tmp1,tmp2))
    dat <- lookup[which(lookup[,2] %in% result$IID),]
    # Write output in the exact format of the original script
    cat("There are", nrow(dat), "related individuals using a kinship threshold of ${kinship}\n")
    write.table(dat, ${_output:r}, quote=FALSE, row.names=FALSE, col.names=FALSE)
    
bash: expand= "${ }", stderr = f'{_output}.stderr', stdout = f'{_output}.stdout', container = container, entrypoint=entrypoint
    i="${_output}"
    output_size=$(ls -lh $i | cut -f 5 -d ' ')
    output_rows=$(cat $i | wc -l | cut -f 1 -d ' ')
    output_column=$(cat $i | head -1 | wc -w)
    output_preview=$(cat $i | grep -v "##" | head | cut -f 1,2,3,4,5,6)
    
    printf "output_info: %s\noutput_size: %s\noutput_rows: %s\noutput_column: %s\noutput_preview:\n%s\n" \
        "$i" "$output_size" "$output_rows" "$output_column" "$output_preview"

# Split genotype data into related and unrelated samples, if related individuals are detected
[king_3]
depends: sos_variable("related_id")
input: output_from(2), genoFile
plink_command = get_plink_command_prefix(_input[1])
output_format = determine_plink_format(_input[1])
make_command = get_make_command(output_format)
output_ext = get_output_extension(output_format)
output: unrelated_bed = f'{cwd}/{_input[0]:bn}.unrelated{output_ext}',
        related_bed = f'{cwd}/{_input[0]:bn}.related{output_ext}'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output[0]:bn}'
bash: expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout', container = container, entrypoint=entrypoint
    plink2 \
      ${plink_command} ${_input[1]:n} \
      --remove ${_input[0]} \
      ${('--keep %s' % keep_samples) if keep_samples.is_file() else ""} \
      ${make_command} \
      --out ${_output[0]:n} \
      --threads ${numThreads} \
      --memory ${int(expand_size(mem) * 0.9)/1e6} --new-id-max-allele-len 1000 --set-all-var-ids chr@:#_\$r_\$a 

    if [ ${len(related_id)} -ne 0 ] ; then
    plink2 \
      ${plink_command} ${_input[1]:n} \
      --keep ${_input[0]} \
      ${make_command} \
      --out ${_output[1]:n} \
      --threads ${numThreads} \
      --memory ${int(expand_size(mem) * 0.9)/1e6} --new-id-max-allele-len 1000 --set-all-var-ids chr@:#_\$r_\$a 
    else
       touch ${_output[1]}
    fi

# Filter SNPs and select individuals 
[qc_no_prune, qc_1 (basic QC filters)]
# minimum MAF filter to use. 0 means do not apply this filter.
parameter: maf_filter = 0.0
# maximum MAF filter to use. 0 means do not apply this filter.
parameter: maf_max_filter = 0.0
# minimum MAC filter to use. 0 means do not apply this filter.
parameter: mac_filter = 0.0
# maximum MAC filter to use. 0 means do not apply this filter.
parameter: mac_max_filter = 0.0 
# Maximum missingess per-variant
parameter: geno_filter = 0.1
# Maximum missingness per-sample
parameter: mind_filter = 0.1
# HWE filter -- a very lenient one
parameter: hwe_filter = 1e-15
# Other PLINK arguments e.g snps_only, write-samples, etc
parameter: other_args = []
# Only output SNP and sample list, rather than the PLINK binary format of subset data
parameter: meta_only = False
# Remove duplicate variants
parameter: rm_dups = False
# Add option to process dosage
parameter: treat_dosage_missing = False

fail_if(not (keep_samples.is_file() or keep_samples == path('.')), msg = f'Cannot find ``{keep_samples}``')
fail_if(not (keep_variants.is_file() or keep_variants == path('.')), msg = f'Cannot find ``{keep_variants}``')
fail_if(not (remove_samples.is_file() or remove_samples == path('.')), msg = f'Cannot find ``{remove_samples}``')

input: genoFile, group_by=1
plink_command = get_plink_command_prefix(_input)
output_format = determine_plink_format(_input)
make_command = get_make_command(output_format) if not meta_only else "--write-snplist --write-samples"
output_ext = get_output_extension(output_format) if not meta_only else ".snplist"
output: f'{cwd}/{_input:bn}{("." + name) if name else ""}.plink_qc{".extracted" if keep_variants.is_file() else ""}{output_ext}'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'
bash: container=container, expand= "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', entrypoint=entrypoint
    plink2 \
      ${plink_command} ${_input:n} \
      --allow-extra-chr \
      ${('--maf %s' % maf_filter) if maf_filter > 0 else ''} \
      ${('--max-maf %s' % maf_max_filter) if maf_max_filter > 0 else ''} \
      ${('--mac %s' % mac_filter) if mac_filter > 0 else ''} \
      ${('--max-mac %s' % mac_max_filter) if mac_max_filter > 0 else ''} \
      ${('--geno %s' % geno_filter) if geno_filter > 0 else ''} ${('dosage' if treat_dosage_missing else '')} \
      ${('--hwe %s' % hwe_filter) if hwe_filter > 0 else ''} \
      ${('--mind %s' % mind_filter) if mind_filter > 0 else ''} ${('dosage' if treat_dosage_missing else '')} \
      ${('--keep %s' % keep_samples) if keep_samples.is_file() else ""} \
      ${('--remove %s' % remove_samples) if remove_samples.is_file() else ""} \
      ${('--exclude %s' % exclude_variants) if exclude_variants.is_file() else ""} \
      ${('--extract %s' % keep_variants) if keep_variants.is_file() else ""} \
      ${make_command} \
      ${("") if not rm_dups else "--rm-dup force-first 'list'"} \
      ${paths(["--%s" % x for x in other_args]) if other_args else ""} \
      --out ${_output:n} \
      --threads ${numThreads} \
      --memory ${int(expand_size(mem) * 0.9)/1e6} --new-id-max-allele-len 1000 --set-all-var-ids chr@:#_\$r_\$a 
        
bash: expand= "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container = container, entrypoint=entrypoint
    i="${_output}"
    output_size=$(ls -lh $i | cut -f 5 -d ' ')
    printf "output_info: %s\noutput_size: %s\n" "$i" "$output_size" >> ${_output:n}.stdout

# LD prunning and remove related individuals (both ind of a pair)
# Plink2 has multi-threaded calculation for LD prunning
[qc_2 (LD pruning)]
# Window size
parameter: window = 50
# Shift window every 10 snps
parameter: shift = 10
parameter: r2 = 0.1
stop_if(r2==0)
plink_command = get_plink_command_prefix(_input)
output_format = determine_plink_format(_input)
make_command = get_make_command(output_format)
output_ext = get_output_extension(output_format, is_prune=True)
output: bed=f'{cwd}/{_input:bn}{output_ext}', prune=f'{cwd}/{_input:bn}.prune.in'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output[0]:bn}'
bash: container=container, expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout', entrypoint=entrypoint
    plink2 \
    ${plink_command} ${_input:n} \
    --allow-extra-chr \
    --rm-dup force-first \
    --indep-pairwise ${window} ${shift} ${r2}  \
    --out ${_output["prune"]:nn} \
    --threads ${numThreads} \
    --memory ${int(expand_size(mem) * 0.9)/1e6}
   
    plink2 \
    ${plink_command} ${_input:n} \
    --allow-extra-chr \
    --extract ${_output['prune']} \
    ${make_command} \
    --out ${_output['bed']:n} \
    --threads ${numThreads} \
    --memory ${int(expand_size(mem) * 0.9)/1e6}
    
bash: expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout', container = container, entrypoint=entrypoint
    i="${_output[0]}"
    output_size=$(ls -lh $i | cut -f 5 -d ' ')
    printf "output_info: %s\noutput_size: %s\n" "$i" "$output_size" >> ${_output[0]:n}.stdout
    i="${_output[1]}"
    output_size=$(ls -lh $i | cut -f 5 -d ' ')
    output_rows=$(zcat $i | wc -l | cut -f 1 -d ' ')
    output_column=$(zcat $i | head -1 | wc -w)
    output_preview=$(cat $i | grep -v "##" | head | cut -f 1,2,3,4,5,6) 
    printf "output_info: %s\noutput_size: %s\noutput_rows: %s\noutput_column: %s\noutput_preview:\n%s\n" \
        "$i" "$output_size" "$output_rows" "$output_column" "$output_preview" >> ${_output[1]}.stdout

# This workflow extracts overlapping samples for genotype data with phenotype data, and output the filtered sample genotype list as well as sample phenotype list
[genotype_phenotype_sample_overlap]
# A genotype fam file
parameter: genoFile = path
# A phenotype file, can be bed.gz or tsv
parameter: phenoFile = path
# If this file is provided, a genotype/phenotype sample name match will be performed
# It must contain two column names: genotype_id, sample_id
parameter: sample_participant_lookup = path(".")
depends: executable('tabix'), executable('bgzip')
input: genoFile, phenoFile
output: f'{cwd:a}/{path(_input[1]):bn}.sample_overlap.txt', f'{cwd:a}/{path(_input[1]):bn}.sample_genotypes.txt'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'
R: expand = "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout', container = container, entrypoint=entrypoint
    # Load required libraries
    library(dplyr)
    library(readr)
    library(data.table)

    # Read data files; let read_delim auto-determine the delimiter
    input_file <- "${_input[0]}"
    
    # Determine if the input is a FAM file (PLINK1) or a PSAM file (PLINK2)
    if (grepl("\\.psam$", input_file)) {
        # For PSAM file (PLINK2 format)
        genoFam <- fread(${_input[0]:ar}, header=TRUE)
        # Standardize column names (PSAM typically uses '#IID' or '#FID IID')
        if ("#FID" %in% colnames(genoFam) && "IID" %in% colnames(genoFam)) {
            genoFam <- genoFam %>% select(`#FID`, IID)
            colnames(genoFam) <- c("V1", "V2")  # Make column names match PLINK1 format for consistency
        } else if ("#IID" %in% colnames(genoFam)) {
            genoFam <- genoFam %>% mutate(V1 = `#IID`) %>% select(V1, `#IID`)
            colnames(genoFam) <- c("V1", "V2")  # Make column names match PLINK1 format for consistency
        } else {
            # Try to handle other possible PSAM formats
            if (ncol(genoFam) >= 2) {
                genoFam <- genoFam[, 1:2]
                colnames(genoFam) <- c("V1", "V2")  # Make column names match PLINK1 format for consistency
            } else {
                stop("Cannot determine the structure of the PSAM file. Expected columns not found.")
            }
        }
    } else {
        # For FAM file (PLINK1 format) - assume standard format with no header
        genoFam <- fread(${_input[0]:ar}, header=FALSE)
    }
    
    phenoFile <- read_delim(${_input[1]:ar}, col_names=TRUE)
    if (${"TRUE" if sample_participant_lookup.is_file() else "FALSE"}) {
        sample_lookup <- fread(${sample_participant_lookup:ar}, header=TRUE)
        colnames(sample_lookup) <- c("sample_id", "genotype_id") # FIXME: This is for the old lookup table with columns c("sample_id", "participant_id")
        # rename phenotype file according to lookup file
        colnames(phenoFile)[-c(1:4)] <- phenoFile %>%
          colnames() %>%
          .[-c(1:4)] %>%
          match(sample_lookup$sample_id) %>%
          sample_lookup$genotype_id[.]
      output_file <- paste0(${phenoFile:nnr}, '.rename_sample.bed')
      # rename phenotype file
      phenoFile$start <- as.integer(phenoFile$start) 
      phenoFile$end <- as.integer(phenoFile$end) 
      fwrite(phenoFile, output_file, sep = '\t')
      # Compress the file using bgzip
      system(paste("bgzip -c", output_file, ">", paste0(output_file, ".gz")))
      # Create the index file using tabix
      system(paste("tabix -p bed", paste0(output_file, ".gz")))
      system(paste("rm", output_file))
    } else {
        sample_lookup <- cbind(genoFam[,2], genoFam[,2])
        colnames(sample_lookup) <- c("genotype_id", "sample_id")
    }
    sample_lookup <- sample_lookup %>%
    filter(
        genotype_id %in% genoFam$V2,
        sample_id %in% colnames(phenoFile)
    )
    
    genoFam %>%
    filter(
        V2 %in% sample_lookup$genotype_id,
    ) %>%
    select(V1, V2) %>%
    fwrite(${_output[1]:r}, col.names=FALSE, sep="\t")

    sample_lookup %>%
    fwrite(${_output[0]:r}, sep="\t")