Hi!

Wondering if I can pick your brains about potential faults in my pipeline. I am trying to find variants (SNPs, indels) in baboon whole genome samples using bcftools.

My plan is to filter based on a depth of greater than 5 and quality of greater than 30. Further, I'd like to filter for allele depth of greater than 2, with a ratio of >.2 (since these animals are diploid I'd expect a heterozygous individual to have alleles at a ratio of .5).

My questions are: 1) Is this an appropriate filter? I am struggling to identify best practices so I was wondering if any one had a standard they would recommend. Is there anything I need to be careful of?

2) I'm confused by multiallelic base calling. I use bcftools norm -m to split multiallelic sites. But I'm confused how that works if the individual in my pipeline is heterozygous at a given locus, with both alleles not matching the reference.

For example: if the individual is heterozygous for a T, C, at a locus, then my file would look like:

chr1 100 . A T,C

Once I split multiallelic sites, it would look like:

chr1 100 . A T

chr1 100 . A C

So then each would be evaluated through my filters separately, correct? And if they didn't pass the filters, would bcftools consensus default to calling the reference?

3) I'm getting an error: The site CM018188.1:6899541 overlaps with another variant, skipping... when I run bcftools consensus. Why would that happen if I have already phased? Should I not split multiallelic sites?

My code is posted below for reference.

# --------------------------------------------------------------
# -----------------------------------------------------------------------------------------------
# Variant calling:
# -----------------------------------------------------------------------------------------------
# Run variant caller
bcftools mpileup \
  -f ${BASE}/ref/papioanubis.fasta \
  -R ${BASE}/bed_files/genes.bed \
  -Ou \
  ${BAMFILE}.rg.bam | \
bcftools call \
  -m \
  --threads 4 \
  -Oz \
  -o ${RESULTS}/${SAMPLE}_calls.vcf.gz

# -----------------------------------------------------------------------------------------------
# Normalize:
# -----------------------------------------------------------------------------------------------
# Normalize VCF to fix overlapping variants and left-align indels
bcftools norm \
  -f ${BASE}/ref/papioanubis.fasta \
  -m -any \
  ${RESULTS}/${SAMPLE}_calls.vcf.gz \
  -Oz -o ${RESULTS}/${SAMPLE}_calls_norm.vcf.gz

# -----------------------------------------------------------------------------------------------
# Filter:
# -----------------------------------------------------------------------------------------------
# FILTER low quality variants
bcftools filter \
  -i 'INFO/DP>5 && QUAL>30' \
  ${RESULTS}/${SAMPLE}_calls_norm.vcf.gz \
  -Oz -o  ${RESULTS}/${SAMPLE}_calls_norm_filtered_pre.vcf.gz

# FILTER low ALLELE depth by filtering out reads with heterozygosity less than .2
# The VAF filter does not exist here so I chose instead to calculate my own
# However, this gets a little tricky when we have two potential alleles that don't match the reference
# require both alt alleles to have depth >= 2
# and allow ratio check to be skipped when ref depth is 0
bcftools filter \
  -i '(AD[0:0] > 0 && AD[0:1] >= 2 && AD[0:1]/(AD[0:0]+AD[0:1]) >= 0.2) || 
      (AD[0:0] = 0 && AD[0:1] >= 2)' \
  ${RESULTS}/${SAMPLE}_calls_norm.vcf.gz \
  -Oz -o ${RESULTS}/${SAMPLE}_calls_filtered_norm.vcf.gz


# -----------------------------------------------------------------------------------------------
# File checks:
# -----------------------------------------------------------------------------------------------
# Check if the file exists and is not empty
if [ ! -s ${RESULTS}/${SAMPLE}_calls_filtered_norm.vcf.gz ]; then
    echo "ERROR: output file is empty for ${SAMPLE}"
    exit 1
fi

# Check if the output is a valid VCF
if ! bcftools view ${RESULTS}/${SAMPLE}_calls_filtered_norm.vcf.gz > /dev/null 2>&1; then
    echo "ERROR: output is not a valid VCF for ${SAMPLE}"
    exit 1
fi

bcftools index -f ${RESULTS}/${SAMPLE}_calls_filtered_norm.vcf.gz

# -----------------------------------------------------------------------------------------------
# Phasing:
# -----------------------------------------------------------------------------------------------
# Phase variants
whatshap phase \
  --sample ${SAMPLE} \
  -o ${RESULTS}/${SAMPLE}_filtered_calls_phased.vcf.gz \
  --reference ${BASE}/ref/papioanubis.fasta \
  ${RESULTS}/${SAMPLE}_calls_filtered_norm.vcf.gz \
  ${BAMFILE}.rg.bam

bcftools index -f ${RESULTS}/${SAMPLE}_filtered_calls_phased.vcf.gz

# -----------------------------------------------------------------------------------------------
# Get consensus:
# -----------------------------------------------------------------------------------------------

# Need to look at low-coverage bases
# I am masking any regions with coverage less than 5 reads
samtools depth -a -b ${BASE}/bed_files/genes.bed ${BAMFILE}.rg.bam \
    | awk '$3 < 5 {print $1"\t"$2-1"\t"$2}' \
    > ${BASE}/qc/${SAMPLE}_low_coverage.bed

# Get the consensus sequence for each chromosome
# Masking any low-coverage area
# Add a chain, which will maintain numbering from the reference sequence
# POTENTIAL ISSUE WITH SKIPPING VARIANTS AT MULTIALLELIC SITES
bcftools consensus \
  --sample ${SAMPLE} \
  -f ${BASE}/ref/papioanubis.fasta \
  -H 1 \
  -m ${BASE}/qc/${SAMPLE}_low_coverage.bed \
  --chain ${BASE}/qc/${SAMPLE}_hap1.chain \
  ${RESULTS}/${SAMPLE}_filtered_calls_phased.vcf.gz \
  > ${RESULTS}/${SAMPLE}_hap1.fasta

bcftools consensus \
  --sample ${SAMPLE} \
  -f ${BASE}/ref/papioanubis.fasta \
  -H 2 \
  -m ${BASE}/qc/${SAMPLE}_low_coverage.bed \
  --chain ${BASE}/qc/${SAMPLE}_hap2.chain \
  ${RESULTS}/${SAMPLE}_filtered_calls_phased.vcf.gz \
  > ${RESULTS}/${SAMPLE}_hap2.fasta