Stacks

Descrição

De acordo com a página do Stacks, Stacks é um pipeline de software para construir loci a partir de sequências de leitura curta, como as geradas na plataforma Illumina. Foi desenvolvido para trabalhar com dados baseados em enzimas de restrição, como RAD-seq, para construção de mapas genéticos e estudos de genômica populacional e filogeografia.

Versões Disponíveis

  • stacks/2.4 (default)

Carregando o Módulo

# Carregar Stacks
module load stacks/2.4

# Verificar instalação
ustacks -h

Pipeline Básico com 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

# Processar cada amostra com ustacks
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
done < samples.txt

# Construir catálogo
job-nanny cstacks -P stacks_output -M popmap.txt \
                 -n 1 -p $SLURM_CPUS_PER_TASK

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

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

Processamento com Script

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

# Lista de amostras
samples=(
    "90415"
    "90414"
    "90426"
    "90427"
    "90422"
    "90420"
)

# Diretório de trabalho
work_dir=./

# Processar cada amostra
i=1
for sample in "${samples[@]}"; do
    echo "Processando amostra $sample (índice $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 "Todas as amostras processadas"

Arquivo de População (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 para Múltiplas Amostras

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

Análise de Populações

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

Filtragem de SNPs

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

# Filtrar SNPs por qualidade
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

# Calcular estatísticas
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

Análise de Resultados

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

# Gerar estatísticas com R
cat > analyze_stacks.R << 'EOF'
library(ggplot2)
library(dplyr)

# Ler sumário de populações
if (file.exists("populations.log")) {
    log_data <- readLines("populations.log")

    # Extrair número de SNPs
    snps_line <- grep("Kept", log_data, value=TRUE)
    cat("SNPs retidos:", snps_line, "\n")
}

# Ler dados de heterozigosidade
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="Distribuição de heterozigosidade",
              x="Taxa de heterozigosidade", y="Frequência")

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

# PCA se houver VCF
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 dos SNPs",
              x="PC1", y="PC2")

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

Rscript analyze_stacks.R

Referências

Ver também