Stacks

Description

According to the page of Stacks, Stacks is a software pipeline for building loci from short read sequences, such as those generated on the Illumina platform. It was developed to work with restriction enzyme-based data, such as RAD-seq, for genetic map construction and population genomics and phylogeography studies.

Available Versions

  • stacks/2.4 (default)

Loading the Module

# Load Stacks
module load stacks/2.4

# Verify installation
ustacks -h

Basic Pipeline with Stacks

submit_stacks_basic.sh
#!/bin/bash
#SBATCH -J stacks_basic
#SBATCH -N 1
#SBATCH -c 8
#SBATCH -t 48:00:00
#SBATCH --mem=32G

export INPUT="samples.txt"
export OUTPUT="stacks_output/"

module load stacks/2.4

mkdir -p stacks_output

# Process each sample with ustacks
i=1
while read sample; do
    job-nanny ustacks -f ${sample}.fq.gz -o stacks_output \
                    -i $i -m 3 -M 2 -p $SLURM_CPUS_PER_TASK
    let "i+=1"
done < samples.txt

# Build catalog
job-nanny cstacks -P stacks_output -M popmap.txt \
                 -n 1 -p $SLURM_CPUS_PER_TASK

# Run sstacks
job-nanny sstacks -P stacks_output -M popmap.txt \
                 -p $SLURM_CPUS_PER_TASK

# Populations
job-nanny populations -P stacks_output -M popmap.txt \
                     -r 0.8 -p 1 --min-maf 0.05 \
                     --write-single-snp --vcf --genepop --structure

Script Processing

submit_stacks_script.sh
#!/bin/bash
#SBATCH -J stacks_script
#SBATCH -N 1
#SBATCH -c 16
#SBATCH -t 72:00:00
#SBATCH --mem=64G
#SBATCH --output=stacks_%j.out

export INPUT="ustacks.sh *.fq.gz"
export OUTPUT="*.tsv.gz"
export LARGE_FILES="true"

module load stacks/2.4
job-nanny ./ustacks.sh
ustacks.sh
#!/bin/bash

# Sample list
samples=(
    "90415"
    "90414"
    "90426"
    "90427"
    "90422"
    "90420"
)

# Working directory
work_dir=./

# Process each sample
i=1
for sample in "${samples[@]}"; do
    echo "Processing sample $sample (index $i)"

    ustacks -p 100 -t gzfastq -m 3 -M 4 -i $i \
            -f ${work_dir}/${sample}.fq.gz \
            -o ${work_dir}

    let "i+=1"
done

echo "All samples processed"

Population Map File (popmap)

popmap.txt
sample1    population_A
sample2    population_A
sample3    population_A
sample4    population_B
sample5    population_B
sample6    population_B
sample7    population_C
sample8    population_C
sample9    population_C

Job Array for Multiple Samples

submit_stacks_array.sh
#!/bin/bash
#SBATCH -J stacks_array
#SBATCH --array=1-9
#SBATCH -N 1
#SBATCH -c 8
#SBATCH -t 24:00:00
#SBATCH --mem=16G

SAMPLES=(
    "sample1"
    "sample2"
    "sample3"
    "sample4"
    "sample5"
    "sample6"
    "sample7"
    "sample8"
    "sample9"
)

SAMPLE=${SAMPLES[$SLURM_ARRAY_TASK_ID-1]}
INDEX=$SLURM_ARRAY_TASK_ID
export INPUT="${SAMPLE}.fq.gz"
export OUTPUT="stacks_${SAMPLE}/"

module load stacks/2.4

mkdir -p stacks_${SAMPLE}

job-nanny ustacks -f ${SAMPLE}.fq.gz -o stacks_${SAMPLE} \
                -i $INDEX -m 3 -M 2 -p $SLURM_CPUS_PER_TASK

Population Analysis

submit_populations.sh
#!/bin/bash
#SBATCH -J populations
#SBATCH -N 1
#SBATCH -c 8
#SBATCH -t 12:00:00
#SBATCH --mem=16G

export INPUT="stacks_output/"
export OUTPUT="populations_output/"

module load stacks/2.4

job-nanny populations -P stacks_output -M popmap.txt \
                     -r 0.8 -p 2 --min-maf 0.05 \
                     --write-single-snp \
                     --vcf --genepop --structure \
                     --fstats --phylip --phylip-var \
                     -t $SLURM_CPUS_PER_TASK

SNP Filtering

submit_filter_snps.sh
#!/bin/bash
#SBATCH -J filter_snps
#SBATCH -N 1
#SBATCH -n 1
#SBATCH -t 04:00:00
#SBATCH --mem=8G

export INPUT="populations.snps.vcf"
export OUTPUT="filtered_snps/"

module load stacks/2.4
module load vcftools/0.1.16

mkdir -p filtered_snps

# Filter SNPs by quality
vcftools --vcf populations.snps.vcf \
         --minQ 30 \
         --max-missing 0.8 \
         --maf 0.05 \
         --hwe 0.001 \
         --recode --recode-INFO-all \
         --out filtered_snps/filtered

# Calculate statistics
vcftools --vcf filtered_snps/filtered.recode.vcf \
         --site-mean-depth \
         --out filtered_snps/depth_stats

vcftools --vcf filtered_snps/filtered.recode.vcf \
         --het \
         --out filtered_snps/heterozygosity

Results Analysis

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

module load R/4.0.2

# Generate statistics with R
cat > analyze_stacks.R << 'EOF'
library(ggplot2)
library(dplyr)

# Read population summary
if (file.exists("populations.log")) {
    log_data <- readLines("populations.log")

    # Extract number of SNPs
    snps_line <- grep("Kept", log_data, value=TRUE)
    cat("SNPs retained:", snps_line, "\n")
}

# Read heterozygosity data
if (file.exists("filtered_snps/heterozygosity.het")) {
    het <- read.table("filtered_snps/heterozygosity.het", header=TRUE)

    p <- ggplot(het, aes(x=HET_RATE)) +
         geom_histogram(bins=30, fill="steelblue", color="black") +
         theme_minimal() +
         labs(title="Heterozygosity distribution",
              x="Heterozygosity rate", y="Frequency")

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

# PCA if VCF exists
if (file.exists("filtered_snps/filtered.recode.vcf")) {
    system("plink --vcf filtered_snps/filtered.recode.vcf --pca --out pca")

    pca <- read.table("pca.eigenvec", header=FALSE)
    colnames(pca)[1:2] <- c("FID", "IID")

    p <- ggplot(pca, aes(x=V3, y=V4)) +
         geom_point() +
         theme_minimal() +
         labs(title="PCA of SNPs",
              x="PC1", y="PC2")

    ggsave("pca_plot.png", p, width=8, height=6)
}
EOF

Rscript analyze_stacks.R

References

See also