Genotyping chip data lift-over to reference genome build GRCh38/hg38 V.2

docker pull seqinfoteam/genotype-liftover-imputation-protocols:<tag>

echo PATH=$PATH:/path/to/plink/:/path/to/R/:/path/to/bcftools/ \
    >> $HOME/.bashrc
source $HOME/.bashrc

echo export BCFTOOLS_PLUGINS=/path/to/bcftools/plugins >> $HOME/.bashrc
source $HOME/.bashrc

install.packages('data.table', type = 'source', 
repos = 'http://Rdatatable.github.io/data.table')

for CHR in {1..23}; do
    bcftools norm -m -any ref_data_chr${CHR}.vcf.gz \
    -Oz -o ref_data_split_multiallelic_chr${CHR}.vcf.gz
done

# Generate a chromosome renaming file
for CHR in {1..23} X ; do 
    echo ${CHR} chr${CHR}
done >> chr_names.txt

# Multiple processing commands piped together
for CHR in {1..22} X; do
    bcftools annotate --rename-chrs chr_names.txt \
        ref_data_split_multiallelic_chr${CHR}.vcf.gz \
        -Oz -o ref_data_chr${CHR}.vcf.gz
done

for CHR in {1..23}; do
    # Drop genotypes
    bcftools view -G ref_data_AF_chr${CHR}.vcf.gz \
    -Oz -o ref_data_drop_genotypes_chr${CHR}.vcf.gz
    # Generate index
    bcftools index -t \
        ref_data_drop_genotypes_chr${CHR}.vcf.gz
done

# List the files
echo ref_data_drop_genotypes_chr{1..23}.vcf.gz | \
tr ' ' '\n' > files_for_concatenation.txt

# Concatenate each chr into a single VCF
bcftools concat -f files_for_concatenation.txt \
-Oz -o ref_data.vcf.gz

wget http://www.well.ox.ac.uk/~wrayner/strand/update_build.sh
wget http://www.well.ox.ac.uk/~wrayner/strand/your_chip_platform_specific_chip_strandfile-b38-strand.zip

FASTA=/path/to/your_reference_genome.fasta
BCFTOOLS_PLUGINS=/path/to/bcftools_installation/plugins

DATASET=your_dataset_prefix

# Convert to VCF utilizing fasta file
plink2 \
    --bfile ${DATASET} \
    --recode vcf id-paste=iid bgz \
    --ref-from-fa \
    --fa ${FASTA} \
    --output-chr chrM \
    --out ${DATASET}

BCFTOOLS_PLUGINS=/path/to/bcftools_installation/plugins

DATASET=your_dataset_prefix
VCF=your_dataset_prefix.vcf.gz

# Extract frequencies with BCFtools
# Requires environment variable BCFTOOLS_PLUGINS
export BCFTOOLS_PLUGINS=${BCFTOOLS_PLUGINS}
bcftools +fill-tags ${VCF}\
     -Oz -o ${DATASET}_AF.vcf.gz -- -t AF
bcftools query \
    -f '%CHROM\t%CHROM\_%POS\_%REF\_%ALT\t%REF\t%ALT\t%INFO/AF\n' \
    ${DATASET}_AF.vcf.gz | \
sed '1iCHR\tSNP\tREF\tALT\tAF' > ${DATASET}.frq

#!/bin/env Rscript --no-save

# Required packages
library(data.table) # For fast fread()

# Input variables
args <- commandArgs(TRUE)
indataset <- args[1]
infile <- args[2]
ref_dataset <- args[3]
refname <- args[4]
af_diff_limit <- as.numeric(args[5])

# Read in the frequency files
chip <- fread(infile, header = T)
ref_data <- fread(ref_dataset, header = T)

# Take an intersection of the reference and chip data 
# based on SNP column (in format CHR_POS_REF_ALT)
isec <- merge(ref_data, chip, by = "SNP")

# Exclude if AF is not within the range
exclude <- !abs(isec$AF.x - isec$AF.y) < af_diff_limit

# Non-ref data variants
nonref <- chip[!(chip$SNP) %in% (isec$SNP)]

# Save the plot as jpg
png(paste0(indataset, "_", refname, "_AF.png"), 
    width = 600, height = 600)
# Plot first all and then excludable variants
plot(isec$AF.x, isec$AF.y, col=1, pch=20,
    main=paste0(indataset, " vs. ", refname, " AF"),
    xlab="Reference data AF",
    ylab="Chip data AF")
points(isec[exclude]$AF.x, isec[exclude]$AF.y, 
    col=2, pch=20)
# Draw a legend
legend("topleft", legend=c(
    paste0("Concordant AF, n = ", nrow(isec[!exclude])),
    paste0("High AF difference, n = ", nrow(isec[exclude])),
    paste0("Non-ref variants, n= ", nrow(nonref))),
    col=c("black", "red", "white"), pch=20, cex=1.2)
dev.off()

# Store the high AF difference variants
output <- rbind(c("chip variants", nrow(chip)),
                c("intersection", nrow(isec)),
                c("intersection/chip variants", 
                     (nrow(isec))/nrow(chip)),
                c("non-ref variants", 
                      nrow(chip) - nrow(isec)),
                c("high AF difference", 
                      nrow(isec[exclude])))

write.table(output, 
    paste0(indataset, "_", refname, "_comparison.txt"),
    quote=F, row.names=F, col.names=F, sep = "\t")

FASTA=/path/to/your_reference_genome.fasta
DATASET=your_dataset_prefix

# Align to reference fasta file and 
# keep only RSIDs from the multiallelic sites
bcftools norm -f ${FASTA} -c ws ${DATASET}.vcf.gz -Ou | \
bcftools view -m 3 -Ou | \
bcftools query -f '%ID\n' \
> ${DATASET}_nonreference_alleles.rsid

plink \
    --bfile ${DATASET} \
    --flip ${DATASET}_nonreference_alleles.rsid \
    --make-bed \
    --out ${DATASET}_nonreference_flipped

#!/bin/env Rscript --no-save

# Required packages
library(data.table)

# Input variables
args <- commandArgs(TRUE)
indataset <- args[1]
infile <- args[2]
ref_dataset <- args[3]

# Read in the frequency files
chip <- fread(infile, header = T)
ref_data <- fread(ref_dataset, header =T)

# Take an intersection of the reference and chip data
# based on SNP column (in format CHR_POS_REF_ALT)
isec <- merge(ref_data, chip, by = "SNP")

# Exclude if AF is not within the range
exclude <- !abs(isec$AF.x - isec$AF.y) < 0.1
discrepant <- isec[exclude]

# Add column with a flipped AF value for the chip data
discrepant$AF_flip <- 1-discrepant$AF.y

# Test if the flipped AF matches to the reference data
flippable <- abs(discrepant$AF.x - discrepant$AF_flip) 
    < 0.1

# Keep only ambiguous alleles
ambiguous <- (discrepant$REF.y %in% c("A", "T") 
             & discrepant$ALT.y %in% c("A", "T")) |
         (discrepant$REF.y %in% c("C", "G") 
             & discrepant$ALT.y %in% c("C", "G"))

# Generate output in format CHR \t POS
output <- strsplit(discrepant[flippable & ambiguous]$SNP, "_")
output <- do.call(rbind, output)[,1:2]

# Store the output
write.table(output, paste0(indataset, "_flippable_AF.txt"),
quote = F, row.names = F, col.names = F, sep = "\t")

DATASET=your_dataset_prefix

Rscript --no-save /path/to/find_ambiguous.R \
    ${DATASET} \
    ${DATASET}.frq \
    /path/to/ref_data.frq

bcftools view -T ${DATASET}_flippable_AF.txt \
    ${DATASET}.vcf.gz | \
bcftools query -f '%ID\n' \
    > ${DATASET}_ambiguous_for_flipping.rsid

plink \
    --bfile ${DATASET} \
    --flip ${DATASET}_ambiguous_for_flipping.rsid \
    --make-bed \
    --out ${DATASET}_ambiguous_flipped

FASTA=/path/to/your_reference_genome.fasta
BCFTOOLS_PLUGINS=/path/to/bcftools_installation/plugins

DATASET=your_dataset_prefix
REFERENCE_VCF=ref_data.vcf.gz

# Export environmental variable BCFTOOLS_PLUGINS
export BCFTOOLS_PLUGINS=${BCFTOOLS_PLUGINS}

# Use fixref with reference data RSIDs
bcftools +fixref ${DATASET}.vcf.gz \
        -Oz -o ${DATASET}_fixref.vcf.gz \
        -- -f ${FASTA} \
        -i ${REFERENCE_VCF}

# Sort the file 
bcftools sort ${DATASET}_fixref.vcf.gz \
    -Oz -o ${DATASET}_sorted.vcf.gz
mv ${DATASET}_sorted.vcf.gz ${DATASET}_fixref.vcf.gz

FASTA=/path/to/your_reference_genome.fasta
BCFTOOLS_PLUGINS=/path/to/bcftools_installation/plugins

DATASET=your_dataset_prefix

# Export environmental variable BCFTOOLS_PLUGINS
export BCFTOOLS_PLUGINS=${BCFTOOLS_PLUGINS}

# Fixref force flip
bcftools +fixref ${DATASET}.vcf.gz \
        -Oz -o ${DATASET}_fixref_force.vcf.gz \
        -- -f ${FASTA} \
        -m flip

FASTA=/path/to/your_reference_genome.fasta

DATASET=your_dataset_prefix

# Align to reference genome
bcftools norm -f ${FASTA} -c ws ${DATASET}.vcf.gz -Ou | \
bcftools view -m 2 -M 2 -Ou | \
bcftools norm -d none -Oz -o ${DATASET}_aligned.vcf.gz