A 'bcftools' script for:

• Extracting SNP data from GBS data in vcf file format
• Filtering out raw SNPs to a usable set of SNPs

I. Extracting the headers:

bcftools view -h  input_file.vcf.gz | awk '/CHROM/{print}' | cut -f 1-5,10-  > output_file_header.txt

h, header; 'awk' stands for "Aho, Weinberger, Kernighan (authors)" which is a printing and row-based filtering function, ‘awk’ prints the output; ‘cut’ selects columns

II. Extracting the SNPs from the GBS data:

bcftools query -f '%CHROM:%POS\t%CHROM\t%POS\t%REF\t%ALT[\t%GT]\n' input_file.vcf.gz | \
  awk '{gsub ("0/0","A");gsub ("0/1","H");gsub ("1/1","B");gsub ("./.","NA");print}'  > output_file_without_filter.txt

'f', format; '%CHROM', 'CHROM' column; '\t' tab separated; '\n' new line; gsub is a function to substitute (sub) something from everywhere in the file (global)- global substitution (gsub); '0', reference allele; '1', alternative allele; 'A', Parent 1 (reference parent in the alignment) and 'B', Parent 2 (this is an example of bi-parental breeding population)

Count the row numbers to keep tract of filtering outcomes later on:

awk 'END { print NR }'  output_file _auto_without_filter.txt

Total number of rows: XXX,XXX

III. Filtering out the extracted SNPs:

• There are several filtering options
• Filtering can be done separately using an individual option to see the outcome of the filter
• Or, filtering can be done using all the options together using '&' in the command line

Filter 1A: Filtering against read depth (DP):

bcftools query -i'FMT/DP>10' -f '%CHROM:%POS\t%CHROM\t%POS\t%REF\t%ALT[\t%GT]\n'  
input_file.vcf.gz |\
  awk '{gsub ("0/0","A");gsub ("0/1","H");gsub ("1/1","B");gsub ("./.","NA");print}'  > output_file_filter1A_DP.txt

'i', include; 'FMT', Format; 'DP', read depth; 'GQ', Genotype quality

Count the row numbers:

awk 'END { print NR }'  output_file_filter1A_DP.txt

Total number of rows: XXX

Filter 1B: Filtering against genotype quality (GQ):

bcftools query -i'FMT/GQ>20' -f '%CHROM:%POS\t%CHROM\t%POS\t%REF\t%ALT[\t%GT]\n'  
input_file.vcf.gz |\
  awk '{gsub ("0/0","A");gsub ("0/1","H");gsub ("1/1","B");gsub ("./.","NA");print}'  > output_file_filter1B_GQ.txt

Count the row numbers:

awk 'END { print NR }'  output_file_filter1B_GQ.txt

Total number of rows: XXX

Filter 1C: Filtering against minor allele frequency (MAF):

bcftools query -i'MAF>=0.05' -f '%CHROM:%POS\t%CHROM\t%POS\t%REF\t%ALT[\t%GT]\n'  
input_file.vcf.gz |\
awk '{gsub ("0/0","A");gsub ("0/1","H");gsub ("1/1","B");gsub ("./.","NA");print}'  > output_file_filter1C_MAF.txt

Count the row numbers:

awk 'END { print NR }'  output_file_filter1C_MAF.txt

Total number of rows: XXX

Filter 1D: Filtering against missing genotypes (F_MISSING):

bcftools query -i'F_MISSING<=0.1' -f '%CHROM:%POS\t%CHROM\t%POS\t%REF\t%ALT[\t%GT]\n'  
input_file.vcf.gz |\
  awk '{gsub ("0/0","A");gsub ("0/1","H");gsub ("1/1","B");gsub ("./.","NA");print}'  > output_file_filter1D_F_MISSING.txt

Count the row numbers:

awk 'END { print NR }'  output_file_filter1D_F_MISSING.txt

Total number of rows: XXX

Filter 1E: Filtering against DP and GQ:

bcftools query -i'FMT/DP>10 & FMT/GQ>20' -f '%CHROM:%POS\t%CHROM\t%POS\t%REF\t%ALT[\t%GT]\n'  
input_file.vcf.gz |\
  awk '{gsub ("0/0","A");gsub ("0/1","H");gsub ("1/1","B");gsub ("./.","NA");print}'  > output_file_filter1E_DP_GQ.txt

'&' enforces filters within a sample; '&&' enforces filters but not necessarily within a sample

Count the row numbers:

awk 'END { print NR }'  output_file_filter1E_DP_GQ.txt

Total number of rows: XXX

Filter 1F: Filtering against DP, GQ, MAF and F_MISSING:

bcftools query -i'FMT/DP>10 & FMT/GQ>20 & MAF>=0.05 & F_MISSING<=0.1' -f '%CHROM:%POS\t%CHROM\t%POS\t%REF\t%ALT[\t%GT]\n'  
input_file.vcf.gz |\
  awk '{gsub ("0/0","A");gsub ("0/1","H");gsub ("1/1","B");gsub ("./.","NA");print}'  > output_file_filter1F_DP_GQ_MAF_F_MISSING.txt

Count the row numbers:

awk 'END { print NR }'  output_file_filter1F_DP_GQ_MAF_F_MISSING.txt

Total number of rows: XXX

Filter 2A: Remove the uncalled, monomorphic, and heterozygous sites from Parent 1 i.e keep only A:

Manually find the column number for Parent 1 in the 'header' text file and confirm:

awk '{print $8}' output_file_header.txt

Return: Parent 1

So, Parent 1 column number is confirmed to be $8; use this in the filtering command below

Filtering:

awk '{ if($8 == "A") print $0;}' output_file_filter1F_DP_GQ_MAF_F_MISSING.txt > output_file_filter2A_clean_Parent 1.txt

Count the row numbers:

awk 'END { print NR }'  output_file_filter2A_clean_Parent 1.txt

Total number of rows: XXX

Filter 3A: Remove the uncalled, monomorphic, and heterozygous sites from Parent 2 i.e keep only B:

Find the column number for Parent 2:

awk '{print $103}' output_file_header.txt

Returns: Parent 2

So, Parent 2 column number is $103

Filtering:

awk '{ if($103 == "B") print $0;}'  output_file_filter2A_clean_Parent 1.txt  > output_file_filter3A_clean_Parent 2.txt

Count the row numbers:

awk 'END { print NR }'  output_file_filter3A_clean_Parent 2.txt

Total number of rows: XXX