Chapter 4 SNP calling and filtering

We assessed the quality of raw reads using multiqc (Ewels et al. 2016) and trimmed them with trimmomatic (Bolger et al. 2014). We kept only pair-end reads without adaptors and a phred score above 15 in a sliding window of 4. Seventy percent of trimmed reads mapped off-targets using bwa (Li & Durbin 2009). We thus mapped trimmed reads on the hybrid reference built for the sequence capture experiment using bwa (Li & Durbin 2009), picard (Broad Institute 2018), samtools (Li et al. 2009) and bedtools (Quinlan & Hall 2010). We called variants for each individual using HaplotypeCaller, aggregated variants using GenomicsDBImport and jointly-genotyped individuals using GenotypeGVCFs all in GATK4 software (Auwera et al. 2013). We filtered biallelic SNPs with a quality above 30, a quality by depth above 2, a Fisher strand bias below 60 and a strand odds ratio above 3 using GATK4 (Auwera et al. 2013). Finally, we filtered individuals and SNPs for missing data with a maximum of 95% and 15% of missing data per individual and SNP, respectively, using plink2 (Chen et al. 2019). We obtained 454,262 biallelic SNPs over 385 individuals without outgroups, that we used for population genetic analyses. Since low-frequency alleles and linkage disequilibrium will bias the number of fixed loci and increase the number of false-positives in genomic scans for outliers (Foll & Gaggiotti 2008), we built a second dataset for quantitative genomics and genomic scans, filtering variants with a minor allele frequency above 5% (18 individuals) and with linkage disequilibrium \(r^2<0.99\). We further removed admixed individuals (see population genetic analyses for criteria) and retained 70,737 biallelic SNPs over 372 individuals.

4.1 Quality Check

We received demultiplexed libraries from sequencing. We checked sequences quality combining already produced fastqc and compared them with originally furnished (i) baits, (ii) targets, and (iii) references:

  1. Multi Quality Check: we used multiqc to combined fastqc inputs for every library (1002 for forward and reverse individuals) and check sequences, counts, quality and GC content
  2. Trimming: we trimmed sequences removing bad quality and adaptors sequences
  3. Targets mapping: we mapped 10 libraries on targets to check proportion of off-targets sequences
  4. Reference mapping: we mapped 10 libraries on hybrid reference to check proportion of off-reference sequences, and assess the need for de novo assembly of captured sequences (in case of a high proportion of off-reference sequences)

4.1.1 Multi Quality Check

We used multiqc to combine fastqc inputs for every library (1002 for forward and reverse individuals) and checked sequences, counts, quality and GC content.

cd ~/Documents/BIOGECO/PhD/data/Eschweilera_Paracou/Sequences/quality
multiqc fastqc
mkdir multiqc
mv multiqc_data/ multiqc_report.html L1.fastqc.tar.gz L2.fastqc.tar.gz multiqc

4.1.1.1 Counts

We had a big heterogeneity of sample representativity (215 000 fold), but 85% of samples had more than 66 6667 sequences (ca 1M targets / 150 bp * 10X). Moreover, duplicated sequences were obviously more present in over-represented individuals, probably more linked to PCR biases than ro sequencing issues.

Sequence counts.

Figure 4.1: Sequence counts.

4.1.1.2 Quality

Sequence quality was very good as the Phred score is above 25 for every base on all positions across all sequences.

Phred score.

Figure 4.2: Phred score.

4.1.1.3 GC content

The mean GC content was 41.5 and only a few sequences had non expected global GC content or GC content across the sequence.

GC content across sequences.

Figure 4.3: GC content across sequences.

GC content within sequences.

Figure 4.4: GC content within sequences.

4.1.2 Trimming

We listed all libraries in text files and trimmed all libraries with trimmomatic in pair end (PE) into paired and unpaired compressed fastq files (fq.gz). We trimmed the adaptor (ILLUMINACLIP) of our protocol (TruSeq3-PE) with a seed mismatch of 2 (mismatched count allowed), a threshold for clipping palindrome of 30 (authorized match for ligated adapters), a threshold for simple clip of 10 (match between adapter and sequence), a minimum adaptor length of 2, and keeping both reads each time (keepBothReads). We trimmed sequences on phred score with a minimum of 15 in a sliding window of 4 (SLIDINGWINDOW:4:15) without trimming the beginning (LEADING:X) or the end (TRAILING:X). Without surprise due to the high quality check of sequencing, trimming resulted in 99.91% of paired trimmed reads compared to raw reads (4.5). Thus the main issue of our dataset was more the representativity of sequences more than their quality.

data.frame(libraries = list.files(file.path(path, "Sequences", "raw"))) %>% 
  mutate(libraries = gsub("_R[12].fastq.gz", "", libraries)) %>% 
  unique() %>% 
  write_tsv(path = file.path(path, "Sequences", "libraries.txt"), col_names = F)
read_tsv(file.path(path, "Sequences", "libraries.txt"), col_names = "Library") %>% 
  sample_n(10) %>% 
  write_tsv(path = file.path(path, "Sequences", "libraries_mapping.txt"), col_names = F)
#!/bin/bash
#SBATCH --time=36:00:00
#SBATCH -J trimming
#SBATCH -o trimming_output.out
#SBATCH -e trimming_error.out
#SBATCH --mem=20G
#SBATCH --cpus-per-task=1
#SBATCH --mail-type=BEGIN,END,FAIL
module load bioinfo/Trimmomatic-0.36
for library in $(cat libraries.txt)
do
  java -jar $TRIM_HOME/trimmomatic.jar PE \
    raw/"$library"_R1.fastq.gz raw/"$library"_R2.fastq.gz \
    trimmed/paired/"$library"_R1_paired.fq.gz trimmed/unpaired/"$library"_R1_unpaired.fq.gz \
    trimmed/paired/"$library"_R2_paired.fq.gz trimmed/unpaired/"$library"_R2_unpaired.fq.gz \
    ILLUMINACLIP:TruSeq3-PE.fa:2:30:10:2:keepBothReads \
    SLIDINGWINDOW:4:15
done
cat trimmed/paired_stat.txt
for file in $(ls trimmed/paired)
do
  zcat trimmed/paired/$file | echo $file"   "$((`wc -l`/4)) >> trimmed/paired_stat.txt
done
cat trimmed/unpaired_stat.txt
for file in $(ls trimmed/unpaired)
do
  zcat trimmed/unpaired/$file | echo $file"   "$((`wc -l`/4)) >> trimmed/unpaired_stat.txt
done
Trimming results.

Figure 4.5: Trimming results.

4.1.3 Targets mapping

We mapped every library on the hybrid reference to check off-reference sequences, and assess the need for de novo assembly, in case many sequences would not map on the reference. Globally we had a low coverage of the reference (median of 19%, Fig. 4.6) but reads were 79% to 88% on-reference (Fig. 4.1) ! Finally, we had a median of 4Mb covered with 10X on reference, which is 4 times what we designed in probes. Consequently, we won’t need de novo assembly and will proceed to read mapping for every library on the built reference, already partly annotated.

#!/bin/bash
#SBATCH --time=36:00:00
#SBATCH -J targetsMapping
#SBATCH -o targetsMapping_output.out
#SBATCH -e targetsMapping_error.out
#SBATCH --mem=20G
#SBATCH --cpus-per-task=1
#SBATCH --mail-type=BEGIN,END,FAIL

module purge
module load bioinfo/bwa-0.7.15
module load bioinfo/picard-2.14.1
module load bioinfo/samtools-1.4
module load bioinfo/bedtools-2.26.0

targets=../Baits/files-Symphonia/target-sequences.fas
bwa index $targets
for library in $(cat libraries_mapping.txt)
do
  rg="@RG\tID:${library}\tSM:${library}\tPL:HiSeq4K"
  bwa mem -M -R "${rg}" -t 16 $targets trimmed/paired/"$library"_R1_paired.fq.gz trimmed/paired/"$library"_R2_paired.fq.gz > targetsMapping/sam/"${library}.sam"
  java -Xmx4g -jar $PICARD SortSam I=targetsMapping/sam/"${library}.sam" O=targetsMapping/bam/"${library}".bam SORT_ORDER=coordinate
  java -Xmx4g -jar $PICARD BuildBamIndex I=targetsMapping/bam/"${library}".bam O=targetsMapping/bam/"${filename}".bai
  samtools index targetsMapping/bam/"${library}".bam
  bedtools bamtobed -i =targetsMapping/bam/"${library}".bam > targetsMapping/bed/"${library}".bed
  bedtools merge -i targetsMapping/bed/"${library}".bed > targetsMapping/merged_bed/"${library}".bed
done
touch readsMappingStat.txt
for file in $(ls bam/*.bam)
do
  samtools flagstat $file | echo $file"   "$(grep "mapped (") >> readsMappingStat.txt
done
Reads alignment coverage on targets. Distribution has been cut at 2000X.

Figure 4.6: Reads alignment coverage on targets. Distribution has been cut at 2000X.

Table 4.1: Reads mapped on targets statistics.
Library Reads mapped Percentage of reads mapped
P7-3-2806 358925 28.22
BCI-SG14 950 27.15
BCI-SG47 16677 27.18
P11-2-240 1064 19.25
P14-2-2842 607526 23.85
P2-2-675 499249 28.01
P4-2-2657 784026 26.77
P5-3-2202 722215 30.28
P6-3-2800 474588 20.19
P6-4-2867 1210288 19.31
P7-3-2806 358925 28.22

4.1.4 Reference mapping

We mapped every library on hybrid reference to check off-reference sequences, and assess de novo usefulness. Globally we had a low coverage of the reference (median of 19%, 4.7) but reads were 79% to 88% on-reference (4.2) ! Finally, we had a median of 4Mb covered with 10X on reference, which is 4 times what we designed in probes. Consequently, we won’t need de novo assembly and will proceed to read mapping for every library on the built reference, already partly annotated.

#!/bin/bash
#SBATCH --time=36:00:00
#SBATCH -J referenceMapping
#SBATCH -o treferenceMapping_output.out
#SBATCH -e referenceMapping_error.out
#SBATCH --mem=20G
#SBATCH --cpus-per-task=1
#SBATCH --mail-type=BEGIN,END,FAIL

module purge
module load bioinfo/bwa-0.7.15
module load bioinfo/picard-2.14.1
module load bioinfo/samtools-1.4
module load bioinfo/bedtools-2.26.0

cat ../../Symphonia_Genomic/neutral_selection/merged.fasta > referenceMapping/reference.fasta
reference=referenceMapping/reference.fasta 
bwa index $reference
for library in $(cat libraries_mapping.txt)
do
  rg="@RG\tID:${library}\tSM:${library}\tPL:HiSeq4K"
  bwa mem -M -R "${rg}" -t 16 $reference trimmed/paired/"$library"_R1_paired.fq.gz trimmed/paired/"$library"_R2_paired.fq.gz > referenceMapping/sam/"${library}.sam"
  java -Xmx4g -jar $PICARD SortSam I=referenceMapping/sam/"${library}.sam" O=referenceMapping/bam/"${library}".bam SORT_ORDER=coordinate
  java -Xmx4g -jar $PICARD BuildBamIndex I=referenceMapping/bam/"${library}".bam O=referenceMapping/bam/"${filename}".bai
  samtools index targetsMapping/bam/"${library}".bam
  bedtools bamtobed -i =referenceMapping/bam/"${library}".bam > referenceMapping/bed/"${library}".bed
  bedtools merge -i referenceMapping/bed/"${library}".bed > referenceMapping/merged_bed/"${library}".bed
done
touch readsMappingStat.txt
for file in $(ls bam/*.bam)
do
  samtools flagstat $file | echo $file"   "$(grep "mapped (") >> readsMappingStat.txt
done
Reads alignment coverage on reference. Distribution has been cut at 2000X.

Figure 4.7: Reads alignment coverage on reference. Distribution has been cut at 2000X.

Table 4.2: Reads mapped on reference statistics.
Library Reads mapped Percentage of reads mapped
BCI-SG14 3232 85.28
BCI-SG47 57142 85.57
P11-2-240 4669 78.74
P14-2-2842 2384919 85.85
P2-2-675 1684774 86.47
P4-2-2657 2717886 85.82
P5-3-2202 2276779 87.75
P6-3-2800 2161522 84.74
P6-4-2867 5707686 83.93
P7-3-2806 1153783 83.91

4.2 Mapping

We proceeded to read mapping for every library on the built reference, already partly annotated:

  1. Repeats merging: 41 libraries were repeated, we merged their fastq files before mapping to increase their information before variant calling
  2. Reads mapping: we mapped every library in pair end with bwa mem on the hybrid reference from Ivan Scotti and Sanna Olsson used to build the targets
  3. Reference sequences: we built bedtools for every alignment in order to list sequences with matches in the reference to be used to reduce the explored reference area in variant calling

4.2.1 Repeats merging

41 libraries were repeated, we merged their FASTQ to increase their information before variant calling. Merging repeats confirmed the presence of all 430 individuals at the end of the alignment (402 from Paracou, 20 from herbariums, and 8 from BCI, Itubera and La Selva).

#!/bin/bash
#SBATCH --time=36:00:00
#SBATCH -J compression
#SBATCH -o compression_output.out
#SBATCH -e compression_error.out
#SBATCH --mem=4G
#SBATCH --cpus-per-task=1
#SBATCH --mail-type=BEGIN,END,FAIL

folder=trimmed.paired.joined/
name=symcapture.trimmed.paired.joined

module purge
mkdir trimmed.paired
for file in $(ls paired/*)
do
  mv $file trimmed.paired/$(basename $(echo $file | sed -e 's/_[[:alpha:]]*-[[:alpha:]]*-[[:alnum:]]*//'))
done
rm -r paired
cp ../libraries.txt ./
cat libraries.txt | sed -e 's/_[[:alpha:]]*-[[:alpha:]]*-[[:alnum:]]*_L00[56]//' | sort | uniq | sed -e 's/-b//' | sort | uniq > libraries.uniq.txt
mkdir trimmed.paired.joined
for ind in $(cat libraries.uniq.txt)
do
  echo $ind
  cat trimmed.paired/$ind*_R1_paired.fq.gz > trimmed.paired.joined/"$ind"_R1_paired.fq.gz
  cat trimmed.paired/$ind*_R2_paired.fq.gz > trimmed.paired.joined/"$ind"_R2_paired.fq.gz
done
rm -r trimmed.paired
tar -zcvf $name.tar.gz $folder

4.2.2 Reads mapping

We mapped every library in pair end with bwa mem on the hybrid reference from Ivan Scotti and Sanna Olsson used to build the targets (32 alignments with 2 processes on 64 cores of 1 node of the genologin computer cluster at Genotoul, Toulouse). We had globally a good mapping with more than 80% of the reads mapped for 98% of the libraries (Figure 4.8).

#!/bin/bash
#SBATCH --time=48:00:00
#SBATCH -J mapping
#SBATCH -o mapping_output.out
#SBATCH -e mapping_error.out
#SBATCH --mem=160G
#SBATCH --cpus-per-task=64
#SBATCH --mail-type=BEGIN,END,FAIL

module purge
module load bioinfo/bwa-0.7.15
module load bioinfo/picard-2.14.1
module load bioinfo/samtools-1.4
module load bioinfo/bedtools-2.26.0

task(){
    echo MAPPING "$1"
    bwa mem -M -R "@RG\tID:$1\tSM:$1\tPL:HiSeq4K" \
        -t 2 \
        reference/reference.fasta \
        trimming/trimmed.paired.joined/"$1"_R1_paired.fq.gz \
        trimming/trimmed.paired.joined/"$1"_R2_paired.fq.gz \
        > mapping/sam/"$1".sam
    rm trimming/trimmed.paired.joined/"$1"_R1_paired.fq.gz
    rm trimming/trimmed.paired.joined/"$1"_R2_paired.fq.gz
    java -Xmx4g -jar $PICARD SortSam \
        I=mapping/sam/"$1".sam \
        O=mapping/bam2/"$1".bam
        SORT_ORDER=coordinate
    rm mapping/sam/"$1".sam
    samtools index mapping/bam2/"$1".bam
}

mkdir mapping/sam
mkdir mapping/bam2
N=32
(
for library in $(cat unmapped.txt)
do 
   ((i=i%N)); ((i++==0)) && wait
   task "$library" &
done
)
rm -r mapping/sam
# folders
mkdir mappingStat
touch readsMappingStat.txt

# test
file=$(ls bam*/*.bam | head -n 1)
module load bioinfo/samtools-1.4 ; samtools flagstat $file | echo $file"   "$(grep "mapped (") >> readsMappingStat.txt
cat readsMappingStat.txt
rm readsMappingStat.txt
touch readsMappingStat.txt

# sarray
for file in $(ls bam*/*.bam); do echo 'module load bioinfo/samtools-1.4 ; samtools flagstat '$file' | echo '$file'"   "$(grep "mapped (") >> readsMappingStat.txt'; done > mappingStat.sh
sarray -J mappingStat -o mappingStat/%j.out -e mappingStat/%j.err -t 1:00:00 --mail-type=BEGIN,END,FAIL  mappingStat.sh

# clean
rm -r mappingStat
Mapping result

Figure 4.8: Mapping result

4.2.3 Reference sequences

We built bedtools for every alignment in order to list sequences with matches in the reference to be used to reduce the explored reference areas in variant calling. 99.98% of reference sequences had at least one library matching (Fig. 4.9). Consequently we used all sequences from the reference in the variant calling. Some reference sections were underrepresented in our libraries but they have been removed at the SNP filtering stage..

# folders
mkdir bed
mkdir bed.out

# test
file=$(ls bam*/*.bam | head -n 1)
module load bioinfo/bedtools-2.26.0 ; bedtools bamtobed -i $file > bed/$(basename "${file%.*}").bed
rm bed/*

# sarray
for file in $(ls bam*/*.bam); do echo 'module load bioinfo/bedtools-2.26.0 ; file='$file' ; bedtools bamtobed -i $file > bed/$(basename "${file%.*}").bed'; done > bed.sh
sarray -J bed -o bed.out/%j.out -e bed.out/%j.err -t 1:00:00 --mail-type=BEGIN,END,FAIL  bed.sh

# clean
rm -r bed.out
  
# statistics
mkdir bed.out
touch referenceMappedStats.txt
for file in $(ls bed/*.bed); do echo "cut $file -f1 | sort | uniq | awk -v file="$(basename "${file%.*}")" '{print \$1, file}' >> referenceMappedStats.txt" ; done > bed.sh
sarray -J bed -o bed.out/%j.out -e bed.out/%j.err -t 1:00:00 --mail-type=BEGIN,END,FAIL  bed.sh
rm -r bed.out
Sequences from reference alignment with reads from all libraries.

Figure 4.9: Sequences from reference alignment with reads from all libraries.

4.3 Variant call

We used GATK as it has apparently similar performance to other variant callers (Supernat et al. 2018) and was more known by Myriam. For that we used the following pipeline:

  1. Variant calling Run the HaplotypeCaller on each sample’s BAM files to create single-sample gVCFs using the .g.vcf extension for the output file.
  2. Data aggregation Aggregate the GVCF files and feed in one GVCF with GenomicsDBImport to be genotyped
  3. Joint genotyping Run GenotypeGVCFs on all of them together to create the raw SNP and indel VCFs that are usually emitted by the callers.

4.3.1 Variant calling

Run the HaplotypeCaller on each sample’s BAM files to create single-sample gVCFs using the .g.vcf extension for the output file.
We used sarray which is much more powerful than sbatch in parallel computing.

# folders
mkdir variantCalling/gvcf5

# test
file=$(ls mapping/bam5/*.bam | head -n 1)
srun --mem=20G --pty bash
module load bioinfo/gatk-4.1.2.0 ; gatk --java-options "-Xmx20G" HaplotypeCaller -R reference/reference.fasta -I $file -O variantCalling/gvcf4/$(basename "${file%.*}").g.vcf.gz  -ERC GVCF
exit
rm variantCalling/gvcf4/*

# sarray
for file in $(ls mapping/bam5/*.bam); do echo "module load bioinfo/gatk-4.1.2.0 ; gatk --java-options \"-Xmx20G\" HaplotypeCaller -R reference/reference.fasta -I $file -O variantCalling/gvcf5/$(basename "${file%.*}").g.vcf.gz  -ERC GVCF"; done > haplo5.sh
mkdir haplo5
sarray -J haplo5 -o haplo5/%j.out -e haplo5/%j.err -t 48:00:00 --mem=20G --mail-type=BEGIN,END,FAIL  haplo5.sh

# clean
rm -r haplo5
rm -r tmp

4.3.2 Data aggregation

We aggregated the GVCF files and fed them into one GVCF database with GenomicsDBImport to be genotyped. Beware, GATK 4.0.0.0 does not deal with multiple intervals when using GenomicsDBImport, so we used GATK 4.1.2.0. We divided the step into several interval files of a maximum of 1000 sequences computed in parallel to speed up the operation. NB, we tested the pipeline with 3 individual haplotypes and 10 intervals of 100 sequences run in parallel; and it took 24 minutes. Consequently with 10 fold more sequences per interval we may increase to 4H, and the effect of 10 fold more individual haplotypes is hard to assess.

Due to a memory overload on the cluster I boosted the sarray to 24G per node beside limiting gatk java session to 20G, still the overload is strange as if gatk was opening a parallel session of 20G java. We should not decrease batch size as a batch of 50 individuals means that we will use 9 batches ! If memory issues persist we may decrease the interval length (currently 1000 sequences) and increase the number of jobs in sarray. We may even decrease intervals to 100 sequences resulting in more than 800 jobs and run them by batch of 100 if they are really more efficient.

Running a first batch of 10 samples on 100 sequences took 35 minutes. Thus 432 samples should take ca 1 day and 40 minutes. This is the first run of DB among 8, so in total DB build should take 8 days ! But on the other hand joint genotyping might be launched on each run as soon as they finish. So we might clean the first vcf and obtain a preview of population genetics structure with the first 100 sequences.

# Sample map
touch sample_map.txt
for file in $(ls gvcf*/*.g.vcf.gz)
do 
   echo -e $(basename "${file%.*}")"\t"$file >> sample_map.txt
done

# seq lists
mkdir reference.sequences.lists
cut ../reference/reference.fasta.fai -f1 > reference.sequences.lists/reference.sequences.list
cd reference.sequences.lists
split -l 100 -d reference.sequences.list reference.sequences_ --additional-suffix=.list
rm reference.sequences.list
ls | wc -l
cd ..

# folders
mkdir tmp
mkdir symcaptureDB

# test
srun --mem=24G --pty bash
file=$(ls reference.sequences.lists/ | head -n 1)
module load bioinfo/gatk-4.1.2.0 ; gatk --java-options "-Xmx20g -Xms20g" GenomicsDBImport --genomicsdb-workspace-path symcaptureDB/"${file%.*}".DB -L reference.sequences.lists/$file --sample-name-map sample_map.txt --batch-size 50 --tmp-dir=tmp
exit
rm -r symcaptureDB/*
rm tmp/*

# sarray

for file in $(ls reference.sequences.lists/); do echo "module load bioinfo/gatk-4.1.2.0 ; gatk --java-options \"-Xmx20g -Xms20g\" GenomicsDBImport --genomicsdb-workspace-path symcaptureDB/\"${file%.*}\".DB -L reference.sequences.lists/$file --sample-name-map sample_map.txt --batch-size 10 --consolidate"; done > combine.sh
split -l 207 -d combine.sh combine_ --additional-suffix=.sh
rm combine.sh
mkdir combine
sarray -J combine -o combine/%j.out -e combine/%j.err -t 48:00:00 --mem=40G --mail-type=BEGIN,END,FAIL combine.sh
rm combine_00.sh
rm -r combine
rm -r tmp

# clean
rm -r combine
rm -r tmp

4.3.3 Joint genotyping

We joint-genotyped individuals with GenotypeGVCFs on all of them together to create the raw SNP and indel VCFs that are usually emitted by the callers. We divided the step into several intervals with a maximum of 1000 sequences computed in parallel to speed up the operation (similarly to the previous step). NB, we tested the pipeline with 3 individual haplotypes and 10 intervals of 100 sequences run in parallel; and it took 6 minutes. Consequently with 10 fold more sequences per interval we may increase to 1H, and the effect of 10 fold more individual haplotypes is hard to assess. Then we merged genotypes of all intervals with GatherVcfs from Picard in one raw VCF to be filtered.

# folders
mkdir tmp
mkdir symcapture.vcf.gz

# test
file=$(ls reference.sequences.lists/ | head -n 1)
srun --mem=20G --pty bash
module load bioinfo/gatk-4.1.2.0 ; gatk --java-options "-Xmx20g" GenotypeGVCFs -R ../reference/reference.fasta -L reference.sequences.lists/$file -V gendb://symcaptureDB/$file -O symcapture.vcf.gz/"${file%.*}".vcf.gz
exit
rm tmp/*
rm symcapture.vcf.gz/*

# sarray

for file in $(ls reference.sequences.lists); do echo "module load bioinfo/gatk-4.1.2.0 ; gatk --java-options \"-Xmx20g\" GenotypeGVCFs -R ../reference/reference.fasta -L reference.sequences.lists/$file -V gendb://symcaptureDB/${file%.*}.DB -O symcapture.vcf.gz/${file%.*}.vcf.gz"; done > genotype.sh
mkdir genotype
sarray -J genotype -o genotype.array/%j.out -e genotype.array/%j.err -t 48:00:00 --mem=20G --mail-type=BEGIN,END,FAIL genotype.array.sh

# clean
rm -r genotype.array
rm -r tmp

# merge
echo -e '#!/bin/bash\n#SBATCH --time=48:00:00\n#SBATCH -J gather\n#SBATCH -o gather.out\n#SBATCH -e gather.err\n#SBATCH --mem=20G\n#SBATCH --cpus-per-task=1\n#SBATCH --mail-type=BEGIN,END,FAIL\nmodule load bioinfo/picard-2.14.1\njava -Xmx20g -jar $PICARD GatherVcfs \' > gather.sh
for file in $(ls symcapture.vcf.gz/*.gz)
do
    echo -e '\tI='$file' \' >> gather.sh
done
echo -e '\tO=symcapture.all.raw.vcf.gz\n' >> gather.sh

4.4 Variant filtering

We filtered the previously produced raw vcf with several steps:

  • Gather raw vcf files which resulted in 26 813 513 variants over 432 individuals
  • Biallelic raw vcf filtering which resulted in 19 242 294 variants over 432 individuals
  • SNP biallelic vcf filtering which resulted in 17 521 879 variants over 432 individuals
  • Filters biallelic snps which resulted in 15 531 866 variants over 432 individuals
  • Missing filtered biallelic snp vcf filtering which resulted in 454 262 variants over 406 individuals
  • Paracou filtered & non missing biallelic snp vcf filtering which resulted in 454 262 variants over 385 individuals

4.4.1 Gather

We first gathered all raw vcf files. Individuals genotyping lost 42 reference scaffolds (over 878, 4%) in few batches of individuals, which blocked the functioning of gatk CombineVariants. We thus removed the variants associated with these 42 reference scaffolds. We obtained 26 813 513 variants.

mkdir out
mkdir missing_ind
for file in $(ls symcapture.filtered.vcf/*.vcf.gz) ; do file=$(basename $file) ; file=${file%.*} ; echo "module load bioinfo/tabix-0.2.5 ; module load bioinfo/vcftools-0.1.15 ; vcftools --gzvcf symcapture.filtered.vcf/$file.gz --missing-indv -c > missing_ind/$file.missing.txt"; done > missingInd.sh
sarray -J missingInd -o out/%j.missingInd.out -e out/%j.missingInd.err -t 1:00:00 --mail-type=BEGIN,END,FAIL missingInd.sh
for file in $(ls missing_ind/*.missing.txt) ; do awk '{{if (NR!=1) print FILENAME"\t"$0}}' $file ; done > missingInd.txt
rm - r out
rm -r missing_ind
echo -e '#!/bin/bash\n#SBATCH --time=48:00:00\n#SBATCH -J gather\n#SBATCH -o gather.out\n#SBATCH -e gather.err\n#SBATCH --mem=20G\n#SBATCH --cpus-per-task=1\n#SBATCH --mail-type=BEGIN,END,FAIL\nmodule load bioinfo/picard-2.14.1\njava -Xmx20g -jar $PICARD GatherVcfs \' > gather.sh
for file in $(cat nonmissing.list)
do
    echo -e '\tI=symcapture.raw.vcf/'$(basename $file)' \' >> gather.sh
done
echo -e '\tO=symcapture.all.raw.vcf.gz\n' >> gather.sh
zcat symcapture.all.raw.vcf.gz | grep "#contig" | wc -l

4.4.2 Biallelic

We then used bcftools to limit data to biallelic variants (--max-alleles 2), resulting in 19 242 294 biallelic variants.

srun --mem=80G --cpus-per-task=8 --pty bash
module load bioinfo/bcftools-1.8
bcftools view --max-alleles 2 --threads 8 symcapture.all.raw.vcf.gz | bgzip -c --threads 8 > symcapture.all.biallelic.vcf.gz
bcftools stats  --threads 8 symcapture.all.biallelic.vcf.gz

4.4.3 SNP

We then used gatk to limit data to biallelic snps, resulting in 17 521 879 biallelic snps.

module load bioinfo/gatk-4.1.2.0
gatk IndexFeatureFile \
  -F symcapture.all.biallelic.vcf.gz
gatk SelectVariants \
  -V symcapture.all.biallelic.vcf.gz \
  -select-type SNP \
  -O symcapture.all.biallelic.snp.vcf.gz
gatk IndexFeatureFile \
  -F symcapture.all.biallelic.snp.vcf.gz
module load bioinfo/bcftools-1.8
bcftools stats  --threads 8 symcapture.all.biallelic.snp.vcf.gz

4.4.4 Filters

We filtered the biallelic snp vcf with following filters (name, filter, description), resulting in 15 531 866 filtered biallelic snps,
using next histograms to set and test parameters values :

  • Quality (QUAL) QUAL < 30: represents the likelihood of the site to be homozygous across all samples, we filter out variants having a low quality score (4.10)
  • Quality depth (QD) QD < 2: filter out variants with low variant confidence (4.10)
  • Fisher strand bias (FS) FS > 60: filter out variants based on Phred-scaled p-value using Fisher’s exact test to detect strand bias (4.10)
  • Strand odd ratio (SOR) SOR < 3: filter out variants based on Phred-scaled p-value used to detect strand bias (4.10)
vcftools --gzvcf reference.sequences_00.vcf.gz --missing-indv -c
vcftools --gzvcf reference.sequences_00.vcf.gz --missing-site -c > missing.txt
vcftools --gzvcf reference.sequences_00.vcf.gz --site-quality -c > QUAL.txt
vcftools --gzvcf reference.sequences_00.vcf.gz \
  --get-INFO AC \
  --get-INFO AF \
  --get-INFO QD \
  --get-INFO FS \
  --get-INFO SOR \
  -c > INFO.txt
Quality, quality by depth, Fisher strand and strand odds ratios for biallelic SNPs.

Figure 4.10: Quality, quality by depth, Fisher strand and strand odds ratios for biallelic SNPs. 

module load bioinfo/gatk-4.1.2.0
gatk VariantFiltration \
  -V symcapture.all.biallelic.snp.vcf.gz \
  --filter-expression "QUAL < 30.0 || QD < 2.0 || FS > 60.0 || SOR > 3.0" \
  --filter-name "FAIL" \
  -O symcapture.all.biallelic.snp.intermediate.vcf.gz
gatk SelectVariants \
  -V symcapture.all.biallelic.snp.intermediate.vcf.gz \
  --exclude-filtered \
  -O symcapture.all.biallelic.snp.filtered.vcf.gz
module load bioinfo/bcftools-1.8
bcftools stats --threads 8 symcapture.all.biallelic.snp.filtered.vcf.gz
gatk IndexFeatureFile \
  -F symcapture.all.biallelic.snp.filtered.vcf.gz

4.4.5 Missing data

Missing data filtering is a bit more tricky because missing data of SNPs and individuals are related, e.g. removing individuals with a lot of missing data results in the decrease of SNPs. Ideally, we wanted to keep all individuals, but this would result in a lot of SNP loss because of least represented individuals. So we needed to choose a threshold for missing data for individuals --mind and SNPs --geno.

module load bioinfo/plink_high_contig_20190905
module load bioinfo/plink2_high_contig_20190905
mkdir filtered
plink2 --threads 8 --memory 80000 \
  --vcf symcapture.all.biallelic.snp.filtered.vcf.gz \
  --allow-extra-chr \
  --make-bed --out filtered/symcapture.all.biallelic.snp.filtered
cd filtered
plink --threads 8 --memory 80000 \
  --bfile symcapture.all.biallelic.snp.filtered \
  --allow-extra-chr --missing --het --freq --pca --freqx \
  --out symcapture.all.biallelic.snp.filtered
Missing data statistics for filtered biallelic SNPs before missing data filtering per individual.

Figure 4.11: Missing data statistics for filtered biallelic SNPs before missing data filtering per individual.

Missing data statistics for filtered biallelic SNPs before missing data filtering per SNP.

Figure 4.12: Missing data statistics for filtered biallelic SNPs before missing data filtering per SNP.

4.4.6 Normal filter

With a maximum of 95% of missing data per individual --mind 0.95 and a maximum of 15% of missing data per SNP -geno 0.15, we obtained 454 262 biallelic filtered snps for 406 individuals.

module load bioinfo/plink_high_contig_20190905
module load bioinfo/plink2_high_contig_20190905
mkdir nonmissing
plink2 --threads 8 --memory 80000 \
  --bfile filtered/symcapture.all.biallelic.snp.filtered \
  --allow-extra-chr \
  --mind 0.95 --geno 0.15 \
  --make-bed --out nonmissing/symcapture.all.biallelic.snp.filtered.nonmissing
cd nonmissing
plink --threads 8 --memory 80000 \
  --bfile symcapture.all.biallelic.snp.filtered.nonmissing \
  --allow-extra-chr --missing --het --freqx --pca \
  --out symcapture.all.biallelic.snp.filtered.nonmissing
Missing data statistics for filtered biallelic SNPs after missing data filtering (95% for individuals and 15% for SNPs) per individual.

Figure 4.13: Missing data statistics for filtered biallelic SNPs after missing data filtering (95% for individuals and 15% for SNPs) per individual.

Missing data statistics for filtered biallelic SNPs after missing data filtering (95% for individuals and 15% for SNPs) per SNP.

Figure 4.14: Missing data statistics for filtered biallelic SNPs after missing data filtering (95% for individuals and 15% for SNPs) per SNP.

Heterozygosity statistics for filtered biallelic SNPs after missing data filtering (95% for individuals and 15% for SNPs) per SNP.

Figure 4.15: Heterozygosity statistics for filtered biallelic SNPs after missing data filtering (95% for individuals and 15% for SNPs) per SNP.

4.4.7 Paracou

Finally, we subseted the filtered and biallelic SNPs to individuals from Paracou only, resulting in 385 remaining individuals (17 lost !).

module load bioinfo/plink2-v2.0_alpha2
mkdir paracou
plink2 \
  --bfile nonmissing/symcapture.all.biallelic.snp.filtered.nonmissing \
  --allow-extra-chr \
  --keep Paracou.fam \
  --make-bed --out paracou/symcapture.all.biallelic.snp.filtered.nonmissing.paracou
  # --recode vcf-iid

References

Auwera, G.A., Carneiro, M.O., Hartl, C., Poplin, R., Angel, G. del, Levy-Moonshine, A., Jordan, T., Shakir, K., Roazen, D., Thibault, J., Banks, E., Garimella, K.V., Altshuler, D., Gabriel, S. & DePristo, M.A. (2013). From FastQ Data to High-Confidence Variant Calls: The Genome Analysis Toolkit Best Practices Pipeline. Current Protocols in Bioinformatics, 43, 483–492. Retrieved from https://onlinelibrary.wiley.com/doi/abs/10.1002/0471250953.bi1110s43

Bolger, A.M., Lohse, M. & Usadel, B. (2014). Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics, 30, 2114–2120. Retrieved from https://academic.oup.com/bioinformatics/article-lookup/doi/10.1093/bioinformatics/btu170

Broad Institute. (2018). Picard Tools. Retrieved from http://broadinstitute.github.io/picard/

Chen, Z.L., Meng, J.M., Cao, Y., Yin, J.L., Fang, R.Q., Fan, S.B., Liu, C., Zeng, W.F., Ding, Y.H., Tan, D., Wu, L., Zhou, W.J., Chi, H., Sun, R.X., Dong, M.Q. & He, S.M. (2019). A high-speed search engine pLink 2 with systematic evaluation for proteome-scale identification of cross-linked peptides. Nature Communications, 10. Retrieved from http://dx.doi.org/10.1038/s41467-019-11337-z

Ewels, P., Magnusson, M., Lundin, S. & Käller, M. (2016). MultiQC: Summarize analysis results for multiple tools and samples in a single report. Bioinformatics, 32, 3047–3048.

Foll, M. & Gaggiotti, O. (2008). A genome-scan method to identify selected loci appropriate for both dominant and codominant markers: A Bayesian perspective. Genetics, 180, 977–993. Retrieved from http://www.ncbi.nlm.nih.gov/pubmed/17246615 http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=PMC1201091 http://www.ncbi.nlm.nih.gov/pubmed/18780740 http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=PMC2567396

Li, H. & Durbin, R. (2009). Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics, 25, 1754–1760.

Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., Marth, G., Abecasis, G. & Durbin, R. (2009). The Sequence Alignment/Map format and SAMtools. Bioinformatics, 25, 2078–2079.

Quinlan, A.R. & Hall, I.M. (2010). BEDTools: A flexible suite of utilities for comparing genomic features. Bioinformatics, 26, 841–842.

Supernat, A., Vidarsson, O.V., Steen, V.M. & Stokowy, T. (2018). Comparison of three variant callers for human whole genome sequencing. Scientific Reports, 8, 17851. Retrieved from http://www.nature.com/articles/s41598-018-36177-7