VCFtools

Description

According to the documentation of VCFtools, 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

# Load VCFtools
module load vcftools/0.1.16

# Verify installation
vcftools --version

Basic Filtering

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

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

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
keep_popA.txt
sample1
sample2
sample3
sample4
sample5

Statistics Calculation

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

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

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

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

See also