Structure

In this section:

Description
Available Versions
Loading the Module
Job Submission
Job Array for Multiple K Values
Multiple Runs for Each K
Input File (mainparams)
Input File (extraparams)
Data File
Results Analysis
Q-matrix Plot
References

Description 

According to the page of Structure, the structure program is a free software package for using multi-locus genotypic data to investigate population structure. Its uses include inferring the presence of distinct populations, assigning individuals to populations, studying hybrid zones, identifying migrants and admixed individuals, and estimating population allele frequencies.

Available Versions 

structure/2.3.4 (default)

Loading the Module 

# Load Structure
module load structure/2.3.4

# Verify installation
structure -h

Job Submission 

submit_structure.sh

#!/bin/bash
#SBATCH -J structure
#SBATCH -N 1
#SBATCH -n 1
#SBATCH -t 24:00:00
#SBATCH --mem=8G

export INPUT="input"
export OUTPUT="output_K5"

module load structure/2.3.4

job-nanny structure -K 5 -o output_K5

Job Array for Multiple K Values 

submit_structure_array.sh

#!/bin/bash
#SBATCH -J structure_array
#SBATCH --array=1-10
#SBATCH -N 1
#SBATCH -n 1
#SBATCH -t 24:00:00
#SBATCH --mem=8G

K=$SLURM_ARRAY_TASK_ID
export INPUT="input"
export OUTPUT="output_K${K}/"

module load structure/2.3.4

mkdir -p output_K${K}
cd output_K${K}
cp ../input .

job-nanny structure -K $K -o output

Multiple Runs for Each K 

submit_structure_replicates.sh

#!/bin/bash
#SBATCH -J structure_reps
#SBATCH --array=1-50
#SBATCH -N 1
#SBATCH -n 1
#SBATCH -t 24:00:00
#SBATCH --mem=8G

# 5 K values × 10 replicates
K=$(( (SLURM_ARRAY_TASK_ID-1) / 10 + 1 ))
REP=$(( (SLURM_ARRAY_TASK_ID-1) % 10 + 1 ))

export INPUT="input"
export OUTPUT="K${K}_rep${REP}/"

module load structure/2.3.4

mkdir -p K${K}_rep${REP}
cd K${K}_rep${REP}
cp ../input .

# Different seeds for each replicate
SEED=$((12345 + SLURM_ARRAY_TASK_ID))

job-nanny structure -K $K -o output -D $SEED

Input File (mainparams)

mainparams

#define INFILE input
#define OUTFILE structure_out
#define NUMINDS 100
#define NUMLOCI 10
#define PLOIDY 2
#define MISSING -9
#define ONEROWPERIND 0
#define LABEL 1
#define POPDATA 1
#define POPFLAG 0
#define LOCDATA 0
#define PHENOTYPE 0
#define EXTRACOLS 0
#define MARKERNAMES 1
#define RECESSIVEALLELES 0
#define MAPDISTANCES 0

Input File (extraparams)

extraparams

#define BURNIN 10000
#define NUMREPS 20000
#define NOADMIX 0
#define LINKAGE 0
#define USEPOPINFO 0
#define LOCPRIOR 0
#define FREQSCORR 1
#define ONEFST 0
#define INFERALPHA 1
#define POPALPHAS 0
#define ALPHA 1.0
#define INFERLAMBDA 0
#define LAMBDA 1.0
#define COMPUTEPROB 1
#define PFROMPOPFLAGONLY 0
#define ANCESTDIST 0
#define STARTATPOPINFO 0
#define METROFREQ 10
#define UPDATEFREQ 1

Data File 

input

Format: Rows represent individuals, columns represent loci with encoded genotypes.

Results Analysis 

analyze_structure.sh

#!/bin/bash
#SBATCH -J analyze_structure
#SBATCH -N 1
#SBATCH -n 1
#SBATCH -t 02:00:00
#SBATCH --mem=4G

module load R/4.0.2

# Create analysis script
cat > analyze_structure.R << 'EOF'
library(ggplot2)
library(dplyr)

# Collect results for different K
results <- data.frame()

for (k in 1:10) {
    files <- list.files(pattern = paste0("K", k, "_rep.*/output_f"),
                       recursive = TRUE, full.names = TRUE)

    if (length(files) > 0) {
        # Extract likelihoods
        for (f in files) {
            lines <- readLines(f)
            like_line <- grep("Estimated Ln Prob of Data", lines, value=TRUE)
            if (length(like_line) > 0) {
                like <- as.numeric(gsub(".*= ", "", like_line))
                results <- rbind(results, data.frame(K=k, Likelihood=like))
            }
        }
    }
}

# Plot average likelihood by K
summary <- results %>%
    group_by(K) %>%
    summarise(mean_like = mean(Likelihood, na.rm=TRUE),
              sd_like = sd(Likelihood, na.rm=TRUE),
              n = n())

p <- ggplot(summary, aes(x=K, y=mean_like)) +
     geom_point(size=3) +
     geom_line() +
     geom_errorbar(aes(ymin=mean_like-sd_like, ymax=mean_like+sd_like), width=0.2) +
     theme_minimal() +
     labs(title="Likelihood by K", x="K", y="Mean likelihood")

ggsave("likelihood_plot.png", p, width=8, height=6)

# Evanno statistics (delta K)
if (nrow(summary) > 2) {
    summary$L_prime <- c(NA, diff(summary$mean_like))
    summary$L_doubleprime <- c(NA, diff(summary$L_prime), NA)
    summary$deltaK <- abs(summary$L_doubleprime) / summary$sd_like

    write.csv(summary, "evanno_stats.csv", row.names=FALSE)

    p2 <- ggplot(summary, aes(x=K, y=deltaK)) +
          geom_point(size=3) +
          geom_line() +
          theme_minimal() +
          labs(title="Delta K (Evanno)", x="K", y="Delta K")

    ggsave("deltaK_plot.png", p2, width=8, height=6)
}
EOF

Rscript analyze_structure.R

Q-matrix Plot 

plot_qmatrix.sh

#!/bin/bash
#SBATCH -J plot_qmatrix
#SBATCH -N 1
#SBATCH -n 1
#SBATCH -t 01:00:00
#SBATCH --mem=2G

module load R/4.0.2

cat > plot_q.R << 'EOF'
library(ggplot2)
library(reshape2)

# Read Structure results file (Q format)
read_structure_q <- function(file, K) {
    lines <- readLines(file)
    # Skip header
    data_lines <- lines[grep("^[0-9]", lines)]

    q_matrix <- matrix(0, nrow=length(data_lines), ncol=K)
    for (i in 1:length(data_lines)) {
        vals <- as.numeric(strsplit(data_lines[i], "\\s+")[[1]])
        q_matrix[i,] <- vals[length(vals) - (K-1):0]
    }

    return(q_matrix)
}

# Example for K=3
q <- read_structure_q("output_K3/output_f", 3)

df <- melt(q)
colnames(df) <- c("Individual", "Population", "Proportion")
df$Population <- factor(df$Population)

p <- ggplot(df, aes(x=Individual, y=Proportion, fill=Population)) +
     geom_bar(stat="identity", width=1) +
     theme_minimal() +
     labs(title="Population structure (K=3)",
          x="Individual", y="Proportion") +
     scale_fill_brewer(palette="Set1")

ggsave("structure_plot_K3.png", p, width=12, height=4)
EOF

Rscript plot_q.R