Whole genome sequencing data of multiple individuals of Pakistani descent

Here we report whole genome sequencing of four individuals (H3, H4, H5, and H6) from a family of Pakistani descent. Whole genome sequencing yielded 1084.92, 894.73, 1068.62, and 1005.77 million mapped reads corresponding to 162.73, 134.21, 160.29, and 150.86 Gb sequence data and 52.49x, 43.29x, 51.70x, and 48.66x average coverage for H3, H4, H5, and H6, respectively. We identified 3,529,659, 3,478,495, 3,407,895, and 3,426,862 variants in the genomes of H3, H4, H5, and H6, respectively, including 1,668,024 variants common in the four genomes. Further, we identified 42,422, 39,824, 28,599, and 35,206 novel variants in the genomes of H3, H4, H5, and H6, respectively. A major fraction of the variants identified in the four genomes reside within the intergenic regions of the genome. Single nucleotide polymorphism (SNP) genotype based comparative analysis with ethnic populations of 1000 Genomes database linked the ancestry of all four genomes with the South Asian populations, which was further supported by mitochondria based haplogroup analysis. In conclusion, we report whole genome sequencing of four individuals of Pakistani descent.


Background & Summary
The completion of Human Genome Project ignited several large scale efforts to characterize variations in the human genome, which led to a comprehensive catalog of the common variants including single-nucleotide polymorphisms (SNPs) and insertions/deletions (indels), across the entire human genome 1,2 . Population-based genome reference datasets played an important role in elucidation of rare variants in specific populations 3,4 . So far, comprehensive genome reference datasets have been reported for African, Japanese, Korean, and Chinese populations [5][6][7][8] .
Advancements in next-generation sequencing technologies have impelled the development of a comprehensive catalog of genetic variants from different ethnic populations [9][10][11][12][13][14][15] . The 1000 Genomes Project reports human genetic variation profiles from 26 ethnic populations, including one Pakistani (Punjabi), two Indian (Gujarati and Telugu), one Bangladeshi (Bengali), and one Sri Lankan (Tamil) population-all descendants of the Indian subcontinent 15 .
Additionally, independent groups have recently published two Indian and two Pakistani genomes with an overall 25-30× sequencing coverage [16][17][18][19] . Recently, the GenomeAsia 100 K project reported genomes of 1,739 individuals, including 113 individuals of Pakistani origin (https://browser.genomeasia100k.org). We previously reported the whole genome sequencing of two Pakistani individuals 20 . Here, we report whole genome sequencing of four individuals of Pakistani descent.
Library preparation and next-generation sequencing. Whole genome sequencing was performed using the Illumina HiSeq X10 (Illumina, San Diego, CA, USA). Briefly, 1.0-2.0 µg of fragmented gDNA was used to prepare paired-end libraries with the TruSeq DNA PCR-Free Library Preparation Kit for four samples (H3, H4, H5, and H6) according to the manufacturer's instructions (Illumina Inc., San Diego, CA). All four libraries were sequenced using Illumina HiSeq X10 in paired-end fashion (2 × 150 bp; Illumina Inc.). The base calls were assigned through Illumina Real-Time Analysis software (Ver. 1.17.20) and binary base call (BCL) files were converted to flat-file format (qseq.txt) using Illumina BCL Converter software (Ver. 1.9.4).
Bioinformatics analysis. Paired-end raw reads were aligned to the human reference genome (GRCh38. p13) using Burrows-Wheeler Aligner-MEM (BWA-MEM; Ver. 0.7.17-r1188) without ALT-aware mode 21 . The quality of the read alignments was examined using CollectAlignmentSummaryMetrics from Picard Tools (Ver. 2.19.0; http://broadinstitute.github.io/picard). The duplicate reads were removed from the mapped reads using MarkDuplicates from Picard Tools. The variants including SNPs and indels were called using the Genome Analysis Tool Kit (GATK; Ver. 4.0) best-practices 12,22 . Briefly, the recalibration of base qualities of input reads was performed using GATK tools (BaseRecalibrator and ApplyBQSR). Subsequently, the SNPs, indels, and genotype of variants were identified using multiple tools i.e. HaplotypeCaller (in GVCF mode), GenotypeGVCFs, and VCFtools (Ver. 0.1.15) 23 . Alignment metrics were generated using CollectAlignmentSummaryMetrics and CollectInsertSizeMetrics from Picard Tools. Genome-wide read coverage was generated using Bedtools (Ver. Variant filtering and annotation. The variants identified through the GATK tool kit were further screened using the high-confidence regions characterized by Genome in a Bottle (GIAB) database 25 . The variants aligned within the large segmental duplication regions of the human genome were discarded while variants mapped to the high-confidence regions of GIAB were used in downstream analyses including Venn diagram generation using VennPainter 26 . Note: An allele (variant) with a minimum of 40% of the total reads mapped to reference allele is considered authentic. The filtered variants were annotated using clinEff (Ver. 1.0 h; http://www.dnaminer. com/clineff.html), a professional version of SNPEff 27 , designed for the prediction of functional effects of variants.
Variant calling. The CNVnator (Ver. 0.4.1) algorithm was used for the identification of copy number variations (CNVs) with a bin size of 1,000 and 10,0000 28 . The GIAB filtered variants (SNPs) were imported into the CNV analysis pipeline for plotting the B-allele frequency (BAF) along the read depths for all deletion and duplication events. ancestry prediction. The ancestral roots of H3, H4, H5, and H6 were examined using the algorithms of Peddy (Ver. 0.3.5) 29 . The study utilized the high-performance computational capabilities of the Biowulf Linux cluster at the National Institutes of Health, Bethesda, MD. PCA plots were created using SNPs genotype information obtained from VCF (variant call format) files (from whole genome sequencing data of H3, H4, H5, and H6) and comparing it with combined ethnic populations from the 1000 Genomes dataset.
In parallel, ancestral roots of H3, H4, H5, and H6 were examined through a comparative analysis with genomes of five different ethnic populations within the 1000 Genomes database. We randomly selected 96 samples from African, Ad Mixed American, East Asian, European, and South Asian populations for comparative analysis by the bcftools-isec algorithm. These variants from 1000 Genomes database and four genomes in VCF format were converted to BCF using bcftools (Ver. 1.8). The BCF files were converted to PLINK format using PLINK (Ver. 1.90b6.18) and PLINK algorithms were used to filter the variants to generate a list of markers in approximate linkage equilibrium for PCA analysis.

Data Records
The next-generation whole genome sequencing raw reads of H3, H4, H5, and H6 have been deposited in the NCBI Sequence Read Archive (SRA) with the accession number PRJNA596295 33 . The chromosomal distribution of the variants identified in H3, H4, H5, and H6 genomes is available at figshare 34 .

technical Validation
The next-generaton whole genome sequencing generated 1344.74, 1110.55, 1200.77, and 1142.35 million total reads for H3, H4, H5, and H6, respectively (Table 1) 33 . Quality control (QC) examination of the sequencing reads revealed that >99% of the sequencing data yielded a PHRED score of 30 or above (PHRED score of 30 represents the probability of 0.001 that the base call is wrong). Subsequent to QC examination and the removal of PCR duplicates (~10-18% of reads were marked duplicates and subsequently removed in downstream analysis), www.nature.com/scientificdata www.nature.com/scientificdata/ the majority of the reads (>99% of reads with a PHRED score ≥ 30) mapped to reference human genome (GRCh38.p13; Table 1 The evaluation of sequencing reads revealed that a significant fraction of the genomes of H3, H4, H5, and H6 exhibited 30-60x read coverage ( Fig. 1 and Table 2). Importantly, 5-6% of the genomes of H3, H4, H5, and H6 were not captured, representing 0x read coverage while approximately, 1%, 2%, and 5% of four genomes exhibited 1-10x, 10-20x, and 20-30x read coverage, respectively ( Fig. 1 and Table 2). A minor fraction i.e. <1% of the genomes of H3, H4, H5, and H6 exhibited 80-100x read coverage ( Fig. 1 and Table 2).  Table 1. Summary of the next-generation whole genome sequencing data.  www.nature.com/scientificdata www.nature.com/scientificdata/ Sequence analysis of the genome of H3 revealed a total of 3,529,659 variants including 3,035,369 SNPs and 494,290 indels. The SNPs were annotated against dbSNP (Ver. 150) that identified 7,553 novel variants (0.21% of the total variants) in the genome of H3 34 . A total of 494,290 indels including 34,869 novel indels (7.05% of the total indels) were identified in the H3 genome 34 .
Sequence analysis of the genome of H4 identifed 3,478,495 total variants including 2,996,403 SNPs and 482,092 indels while annotation of the SNPs identified 6,631 novel SNPs (0.19% of the total variants) in the genome of H4 34 . A total of 482,092 indels including 33,193 novel indels (6.88% of the total indels) were identified in the genome of H4 34   www.nature.com/scientificdata www.nature.com/scientificdata/ variants including 2,972,863 SNPs and 453,999 indels while annotation of the SNPs identified 6,703 novel SNPs (0.19% of the total variants) in the genome H6 34 . A total of 453,999 indels including 28,503 novel indels (6.28% of the total indels) were identified in the genome of H6 34 . www.nature.com/scientificdata www.nature.com/scientificdata/ Importantly, we identified a total of 1,668,024 variants including 1,666,232 variants reported previously and 1,792 novel SNPs common in the four genomes (Fig. 2a-c). Altogether, the variants common in the four genomes constitute nearly half of the total variants identified in each genome.
We examined the putative effect of the variants based on their location in the genome (exon, intron, etc.), functional impact (high, moderate, and low), and classification (synonymous vs. non-synonymous), etc. The analysis revealed that intergenic regions harbor the majority of SNPs consistent with the GIAB high-confidence variants. Furthermore, in contrast to intergenic variants, fewer variants were identified in the exons, splice site, and untranslated regions (UTRs) of the genome. Furthermore, >3 K, >82 K, and >28 K variants present in all four genomes were predicted to exhibit a putative high, moderate, and low impact, respectively.
We used CNVnator, an algorithm to characterize copy number variations (CNVs), to examine structural variants in the genomes of H3, H4, H5, and H6. The analysis identified a total of 4,269 copy number variation regions (CNVRs) common in four genomes, covering 305.95 Mb (9.53%) of the reference human genome (GRCh38.p13).
Although H3, H4, H5, and H6 belong to the Punjabi ethnic group of Pakistani population suggesting a close ancestral relationship with South Asian populations, we sought of genomic evidence to confirm the ancestral roots of the four genomes. We compared the SNP genotypes of H3, H4, H5, and H6 to the combined population of the 1000 Genomes project by the Peddy algorithm. The analysis localized the all four genomes within South Asian populations in principal component 1 and 3 (PC1 and PC3) (Fig. 3a-d; arrows pointing to samples shown as red circles in PCA plots) and on the edge of the South Asian populations in principal component 2 (PC2) towards the European populations (Fig. 3a-d). The localization of H3, H4, H5, and H6 in PC2 suggests some ancestral link with European populations.
In parallel, we performed an additional comparative analysis of the four genomes with the genomes of five different ethnic populations in the 1000 Genomes database. The analysis identified >94% overlap of variants in the genomes of H3, H4, H5, and H6 with South Asian populations (Table 3) with a small number of variants in the genomes of H3 (157,110), H4 (166,633), H5 (159,515), and H6 (163,635) genomes not present in South Asian populations (Table 3). We identified > 92% overlap of variants in the genomes of H3, H4, H5, and H6 with both European and Ad Mixed American populations (Table 3). Likewise, we identified > 88% and > 90% overlap of variants in the genomes of H3, H4, H5, and H6 with East Asian and African populations, respectively (Table 3). These data support the notion that H3, H4, H5, and H6 have a close ancestral relationship with South Asian populations (Fig. 4a-c).
www.nature.com/scientificdata www.nature.com/scientificdata/ To further confirm the results of SNP genotype based analysis, we performed mitochondria and Y chromosome based haplogroup analysis. The mitochondria genome analysis revealed M35b haplogroup in the H3 genome and M6 haplogroup in H4, H5, and H6 genomes. Both mitochondrial haplogroups (M35b and M6) have been mainly identified in South Asian populations 35,36 . The Y chromosome analysis identified G1a1b2a haplogroup in H3 and H5 genomes, suggesting a Middle Eastern origin. Taken together, the mitochondria haplogroup based analyses support the results of the SNP genotype based analysis and strengthen the notion that H3, H4, H5, and H6 have a close ancestral relationship with South Asian populations.
In conclusion, we have completed next-generation based whole genome sequencing of four individuals from a family of Pakistani descent. Importantly, nearly 1% of the total variants identified in each of the four genomes are novel and have not been reported previously. To the best of our knowledge, this is the first report of whole genome sequencing of four individuals from a family.