Genome-wide haplotype counting in diatoms S. robusta and P. tricornutum - processed datasets

SNP calling: SNP calling on ILLUMINA short read sequencing was done using GATK HaplotypeCaller 3.7.0. In short, adapters and reads with the quality score below 20 were removed from ILLUMINA reads using BBduk2 with minlen=35 qtrim=rl trimq=20 hdist=1 tbo tpe options and custom adapter reference file. Next, the respective reads were aligned to S. robusta v1 assembly CAICTM010000001-CAICTM010004752 (European Nucleotide Archive) or P. tricornutum v2 assembly2 GCA_000150955.2 (European Nucleotide Archive) using Burrows-Wheeler Alignment Tool (BWA) algorithm BWA-MEM with -M option. Unmapped and multi-mapped reads were removed using SAMtools view with -h -F 4 -q 1 options. Aligned reads were then sorted using picard-tools 1.8.0 SortSam and duplicate reads were marked with MarkDuplicates and indexed with BuildBamIndex. Read base quality scores were adjusted by two round of recalibration. Here, SNPs and indels were called by GATK HaplotypeCaller and filtered with a set of hard filters using SelectVariants; QD < 2.0, FS > 60.0, MQ < 40.0, MQRankSum < -12.5, ReadPosRankSum < -8.0 for SNPs and QD < 2.0, FS > 200.0, ReadPosRankSum < -20.0 for indels. Recalibration table was generated with BaseRecalibrator and recalibrated reads were printed with PrintReads. After a second round of recalibration, germline SNPs were called using HaplotypeCaller with --genotyping_mode DISCOVERY. Next, reliable biallelic SNPs were selected using SelectVariants with --restrictAllelesTo BIALLELIC -selectType SNP and QD < 2.0, QUAL < 30.0, SOR > 3.0, FS > 60.2, MQ < 40.0, MQRankSum < -12.5, ReadPosRankSum < -8.0, AF > 0.2 and DP < 10 options. Repeat regions and low complexity DNA sequence in S. robusta and P. tricornutum were identified using RepeatModeler 1.0.9 and masked using RepeatMasker 4.0.5 and SNPs in these regions were removed from the dataset using BEDtools subtract algorithm. Finally, selected fields (CHROM, POS, REF, ALT) from the SNPs dataset were extracted from the vcf file to a table and split into independent files by contig/chromosome using awk. S. robusta PacBio reads processing: Circular Consensus Sequences (CCS) were obtained with smrtanalysis 2.3.0 (PacBio) with minFullPasses 0 option. CCS reads were then self-correct using canu1.4 with canu_correct genomeSize=136.0m errorRate=0.035 -pacbio-raw options and trimmed with canu_trim genomeSize=136.0m errorRate=0.035 -pacbio-corrected options. Corrected reads were mapped to the reference genome using BLASR with -sam -clipping soft options. The CIGAR string was corrected with samfixcigar, soft-clipped bases were removed with biostar84452 from jvarkit and uniquely mapped reads with mapping quality >20 were selected using SAMtools. Coverage was estimated using GATK 3.7.0 DepthOfCoverage. SAMtools view and awk were used to split the PacBio reads to separate files per contig. P. tricornutum MinION reads processing: MinION reads were self-corrected using canu version 1.4 with canu_correct genomeSize=30m errorRate=0.144 -nanopore-raw options. Reads were aligned to the genome using GraphMap with default settings and uniquely mapped reads were selected using SAMtools view. The CIGAR string was corrected with samfixcigar, soft-clipped bases were removed with biostar84452 from jvarkit. Coverage was estimated using GATK 3.7.0 DepthOfCoverage. SAMtools view and awk were used to split the PacBio reads to separate files per contig. Haplotype counting: The haplotype counting was done with a custom script based on bash, awk and Sam2Tsv from jvarkit (https://doi.org/10.5281/zenodo.4001752). In short, a record for every base in each processed PacBio/MinION read with the position, reference and the actual base was obtained using Sam2Tsv from jvarkit. Next, only positions of SNPs selected in ILLUMINA reads were retained. The record was divided into fixed windows of max 1kb from the first SNP and haplotypes for selected sites were written for each read separately. Reads containing indel or another base than reference and alternative base at selected SNP position or not covering the 1kb region were removed and the number of haplotypes and number of supporting reads for each haplotype was counted using awk. Loci with multiple haplotypes were selected from the record of the number of haplotypes per 1 kb with the number of supporting reads with following conditions: at least three haplotypes had to be supported each by at least 2 reads, the locus had to be at least 100 bp long and the coverage had to be below 100x in order to remove repeat regions that were not masked. Visualization of haplotype counting data was done using Circos and karyoploteR.