.. _vcftools:

========
VCFtools
========

.. contents:: In this section:
   :local:
   :depth: 2

Description
===========

According to the `documentation of VCFtools <https://vcftools.github.io/>`_, VCFtools is a set of tools for working with VCF (Variant Call Format) files, used to manipulate and analyze genetic variants. It allows filtering, comparing, and summarizing large-scale genetic variation data.

Available Versions
==================

* vcftools/0.1.16 (default)

Loading the Module
==================

.. code-block:: bash

   # Load VCFtools
   module load vcftools/0.1.16

   # Verify installation
   vcftools --version

Basic Filtering
===============

.. code-block:: bash
   :caption: submit_vcftools_filter.sh

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

   export INPUT="variants.vcf"
   export OUTPUT="filtered.vcf"

   module load vcftools/0.1.16

   job-nanny vcftools --vcf variants.vcf \
                     --minQ 30 \
                     --max-missing 0.8 \
                     --maf 0.05 \
                     --hwe 0.001 \
                     --recode --recode-INFO-all \
                     --out filtered

Depth Filtering
===============

.. code-block:: bash
   :caption: submit_vcftools_depth.sh

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

   export INPUT="variants.vcf"
   export OUTPUT="depth_filtered.vcf"

   module load vcftools/0.1.16

   job-nanny vcftools --vcf variants.vcf \
                     --min-meanDP 10 \
                     --max-meanDP 100 \
                     --minDP 5 \
                     --recode --recode-INFO-all \
                     --out depth_filtered

Population Filtering
====================

.. code-block:: bash
   :caption: submit_vcftools_pop.sh

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

   export INPUT="variants.vcf keep_popA.txt"
   export OUTPUT="popA.vcf"

   module load vcftools/0.1.16

   job-nanny vcftools --vcf variants.vcf \
                     --keep keep_popA.txt \
                     --recode --recode-INFO-all \
                     --out popA

.. code-block:: text
   :caption: keep_popA.txt

   sample1
   sample2
   sample3
   sample4
   sample5

Statistics Calculation
======================

.. code-block:: bash
   :caption: submit_vcftools_stats.sh

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

   export INPUT="filtered.vcf"
   export OUTPUT="stats/"

   module load vcftools/0.1.16

   mkdir -p stats

   # Allele frequencies
   job-nanny vcftools --vcf filtered.vcf --freq --out stats/freq

   # Site mean depth
   job-nanny vcftools --vcf filtered.vcf --site-mean-depth --out stats/depth

   # Site quality
   job-nanny vcftools --vcf filtered.vcf --site-quality --out stats/quality

   # Individual heterozygosity
   job-nanny vcftools --vcf filtered.vcf --het --out stats/het

   # Pairwise distance
   job-nanny vcftools --vcf filtered.vcf --hardy --out stats/hardy

   # Fst statistics
   job-nanny vcftools --vcf filtered.vcf --weir-fst-pop popA.txt \
                     --weir-fst-pop popB.txt --out stats/fst_popA_popB

Format Conversion
=================

.. code-block:: bash
   :caption: submit_vcftools_convert.sh

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

   export INPUT="filtered.vcf"
   export OUTPUT="converted/"

   module load vcftools/0.1.16

   mkdir -p converted

   # PLINK
   job-nanny vcftools --vcf filtered.vcf --plink --out converted/plink

   # PED (pedigree format)
   job-nanny vcftools --vcf filtered.vcf --plink-tped --out converted/ped

   # Loci
   job-nanny vcftools --vcf filtered.vcf --ldhat --out converted/ldhat

   # BeaST
   job-nanny vcftools --vcf filtered.vcf --BEAGLE-PL --out converted/beagle

   # Impute
   job-nanny vcftools --vcf filtered.vcf --IMPUTE --out converted/impute

Job Array for Multiple Chromosomes
==================================

.. code-block:: bash
   :caption: submit_vcftools_array.sh

   #!/bin/bash
   #SBATCH -J vcf_array
   #SBATCH --array=1-22
   #SBATCH -N 1
   #SBATCH -n 1
   #SBATCH -t 06:00:00
   #SBATCH --mem=8G

   CHR=$SLURM_ARRAY_TASK_ID
   export INPUT="chr${CHR}.vcf"
   export OUTPUT="chr${CHR}_filtered/"

   module load vcftools/0.1.16

   mkdir -p chr${CHR}_filtered

   job-nanny vcftools --vcf chr${CHR}.vcf \
                     --minQ 30 \
                     --max-missing 0.8 \
                     --maf 0.05 \
                     --recode --recode-INFO-all \
                     --out chr${CHR}_filtered/chr${CHR}

PCA Analysis
============

.. code-block:: bash
   :caption: submit_vcf_pca.sh

   #!/bin/bash
   #SBATCH -J vcf_pca
   #SBATCH -N 1
   #SBATCH -n 1
   #SBATCH -t 08:00:00
   #SBATCH --mem=16G

   export INPUT="filtered.vcf"
   export OUTPUT="pca_results/"

   module load vcftools/0.1.16
   module load R/4.0.2

   mkdir -p pca_results

   # Extract genotype matrix from VCF
   vcftools --vcf filtered.vcf --012 --out pca_results/genotypes

   # PCA with R
   cat > pca_analysis.R << 'EOF'
   library(ggplot2)

   # Load data
   genotypes <- as.matrix(read.table("pca_results/genotypes.012", header=FALSE))
   samples <- read.table("pca_results/genotypes.012.indv", header=FALSE)$V1

   # Replace missing values (-1) with column mean
   for (i in 1:ncol(genotypes)) {
       col <- genotypes[,i]
       col_mean <- mean(col[col != -1], na.rm=TRUE)
       col[col == -1] <- col_mean
       genotypes[,i] <- col
   }

   # PCA
   pca_result <- prcomp(genotypes, scale.=TRUE, center=TRUE)

   # Explained variance
   var_explained <- summary(pca_result)$importance[2, ] * 100

   # Save coordinates
   write.table(pca_result$x[,1:10], "pca_results/pca_coords.txt", 
               row.names=samples, col.names=paste0("PC",1:10), quote=FALSE)

   # Save explained variance
   write.table(var_explained[1:10], "pca_results/variance_explained.txt", 
               row.names=FALSE, col.names=FALSE)

   # Plot PC1 vs PC2
   pca_df <- data.frame(pca_result$x[,1:2])
   colnames(pca_df) <- c("PC1", "PC2")
   pca_df$Sample <- samples

   p <- ggplot(pca_df, aes(x=PC1, y=PC2)) +
        geom_point(size=2, alpha=0.7) +
        theme_minimal() +
        labs(x = paste0("PC1 (", round(var_explained[1], 1), "%)"),
             y = paste0("PC2 (", round(var_explained[2], 1), "%)"),
             title = "PCA of SNPs") +
        theme(plot.title = element_text(hjust=0.5))

   ggsave("pca_results/pca_plot.png", p, width=8, height=6)

   # Scree plot
   scree_df <- data.frame(PC = 1:10, Var = var_explained[1:10])

   p2 <- ggplot(scree_df, aes(x=PC, y=Var)) +
         geom_bar(stat="identity", fill="steelblue") +
         theme_minimal() +
         labs(x = "Principal Component", y = "Explained Variance (%)",
              title = "Scree Plot") +
         theme(plot.title = element_text(hjust=0.5))

   ggsave("pca_results/scree_plot.png", p2, width=8, height=6)

   print("PCA completed successfully!")
   EOF

   Rscript pca_analysis.R

References
==========

- Documentation: https://vcftools.github.io/
- Manual: https://vcftools.github.io/man_latest.html
- GitHub: https://github.com/vcftools/vcftools

.. seealso::

   - :ref:`stacks` - RAD-seq pipeline
   - :ref:`running_simulations` - How to submit jobs