Assembling unmapped reads reveals hidden variation in South Asian genomes

Domain: neuroscience
Published: 2026-04-10
Editorial outcome: major_revision

Abstract
Conventional genome mapping-based approaches systematically overlook genetic variation, particularly in regions that substantially differ from the reference. To explore this hidden variation, here we examine unmapped and poorly mapped reads from the genomes of 640 human individuals from South Asian populations in the 1000 Genomes Project and the Simons Genome Diversity Project. We assemble tens of megabases of non-redundant sequence in tens of thousands of large contigs, a significant portion of which is present in both South Asian and other populations. We demonstrate that much of this sequence is not discovered by traditional variant discovery approaches even when using complete genomes and pangenomes. Across 20,000 placed contigs, we find 8,215 intersections with 106 protein coding genes and over 15,000 placements within 1 kbp of a known GWAS hit. We use long-read data from a subset of samples to validate the majority of their assembled sequences, align RNA-seq data to identify hundreds of unplaced contigs with transcriptional potential, and query existing nucleotide databases to infer the origins of the remaining unplaced sequences. Our results highlight the limitations of even the most complete reference genomes and provide a model for understanding the distribution of hidden variation in any human population.

Introduction
Genome sequencing studies have not uniformly spanned the spectrum of human diversity. As of 2019, >80% of participants in genomic studies were of predominantly European ancestries, with poor representation from other geographic regions 1,2. While recent sequencing and assembly efforts are now addressing this gap 3–5, the persisting biases hinder fundamental biological research and clinical interpretation 6,7.

One such underrepresented region is South Asia, home to more than 1.9 billion people 8, with another 35-40 million people in its diaspora 9–11. Population genomic research focused on South Asia has uncovered deep population structure, runs of homozygosity 12 attributed to founder effects and endogamy, as well as intricate histories of both ancient and recent admixture, all of which are influenced by societal and cultural factors 13–17. Although these studies are invaluable, South Asians compose only ~2% of individuals in genome-wide association studies (GWAS), and even within this broader group, certain ethnicities are over- or underrepresented 1 (Supplementary Figure 1). Prior efforts to catalog South Asian diversity are largely out-dated, based on genotyping arrays or low-coverage sequencing data, or performed using historical versions of the reference genome18,19. While there is an abundance of low-coverage sequencing data from individuals of South Asian ancestries, curated, high-coverage, open-access datasets from South Asian samples are lacking. Two notable exceptions are the 1000 Genomes Project (1KGP) and Simons Genome Diversity Project (SGDP), which offer high-coverage whole-genome sequencing data from 640 South Asian individuals across 24 populations (Figure 1a-b).

Even with high-quality data, researchers face additional biases that impede genetic studies. Most conventional genomic analyses involve alignment to a reference genome, but a single reference genome may obscure variation in regions poorly represented by the reference. Even for the complete T2T-CHM13 reference genome, insertions relative to the reference and highly divergent haplotypes may fail to align 4,5,20. This reference bias was highlighted in the African Pangenome project (APG) 21, which found hundreds of megabases of sequence present in African individuals that were absent from GRCh38. Similar efforts with other populations have similarly identified tens of megabases of sequence that are missing in widely used reference genomes 22–24. Even current pangenomes cannot fully overcome this challenge, as human populations possess an excess of ultra-rare variation that will not be adequately captured by any modest-sized panel of references 3,25.

In this work, we cataloged genetic diversity within South Asian populations and tested the limitations of reference genomes and pangenomes (Figure 1, Supplementary Figure 2). To this end, we re-analyzed existing short read sequencing data from 640 South Asian individuals across 24 populations in the 1KGP 26,27 and SGDP 28 datasets and assembled sequences that are present in these samples but absent from reference genomes (Figure 1c). This process involved assembling unmapped reads into larger contigs, which has been demonstrated to reveal variant sequences in both human and non-human samples 21,29–33. We then validated thousands of large contigs assembled from these unmapped reads through alignment to long-read based assemblies 34, effectively anchoring approximately 20,000 contigs on T2T-CHM13 and facilitating interpretation. Many of these placements overlapped with biologically relevant regions, including thousands of intersections with protein coding genes or in proximity to GWAS hits. In contrast, conventional pipelines failed to detect and place most of these large contigs and ignored most of the unplaced assembled sequence. Through alignment of existing RNA-seq data, we also uncovered several contigs with transcriptional evidence, underscoring their functional and phenotypic potential. Together, our research highlights the limitations of both single and pangenome references and presents a generalizable framework and open-source software pipelines for recovering hidden variation from existing sequencing data.

Results
Alignment and assembly of data from the 1KGP SAS set
We first aligned reads from the 601 1KGP SAS individuals against the T2T-CHM13 v2.0 4 and GRCh38 35 linear references. Alignment to the recent T2T-CHM13 reference improved the alignment rate (+0.7-1.0%) compared to GRCh38 and increased the number of reads concordantly aligned to a single locus (+1.0%) (Figure 2a). This improvement was expected, since T2T-CHM13 adds hundreds of millions of bases of resolved sequence. However, despite this improvement, 1.3-1.7% of all reads per individual remain unaligned against T2T-CHM13. These reads represent an opportunity for discovery, as they may contain variant sequences these individuals possess relative to the chosen reference, as well as reads that fail to align due to sequencing errors or other technical artifacts. In addition to extracting fully unaligned reads, we also examined reads with poor alignment quality scores (MAPQ < 20), which constituted 20-30x the number of fully unaligned reads. Results from the alignment and assembly of these reads are provided in Supplementary Note 1.

Unaligned reads were then assembled into larger contigs whose origins could be more easily assessed. Using the MEGAHIT genome assembler 36 on reads that failed to align against T2T-CHM13, we initially assembled ~1 Mbp of sequence per individual. However, much of this sequence was composed of short contigs, potentially reflecting contamination or assembly errors. Using an approach inspired by the African Pangenome study 21, this set was filtered to exclude contigs <1 kbp in length. This reduced the amount of assembled sequence to approximately 550 kbp per individual. These contigs were then screened using BLAST 37 and Centrifuge 38 against databases of known bacterial and microbial sequences, removing contigs with any matches to these contaminants. This step removed few contigs (< 0.5%), but increased confidence that the remaining sequences were of human origin. Of the contaminant contigs that were removed, the majority matched to adapter sequences or common human contaminants such as herpes virus.

After filtering, we assembled an average of ~550 kbp of sequence per SAS individual, totaling approximately 350 Mbp across the entire 1KGP SAS set (Figure 2b-c). While three outlier individuals were apparent, self-alignments of their contigs revealed nearly-duplicated assembled sequence, thereby resolving these cases (Figure 2b). Across all samples the size of these contigs varied greatly, with an average of 70 contigs >5 kbp and 18 contigs >10 kbp per individual. The average N50 for these filtered assemblies was 1.95 kbp, and the average size of the largest contig assembled per individual was 15-18 kbp, with the largest contigs in the entire set spanning over 50 kbp. The analogous approach using GRCh38 assembled an average of 2 Mbp of sequence per individual, significantly exceeding the above results from T2T-CHM13 (two-tailed Welch’s t-test, t = -134.45, p-value < 1 x 10-308), and underscoring the greater completeness of the T2T reference.

Application of RepeatMasker 39 showed that 4.39% of the assembled sequence consisted of short interspersed nuclear elements (SINEs) (1.73% ALUs, 2.66% MIRs), 17.89% long interspersed nuclear elements (LINEs) (13.85% LINE1, 3.89% LINE2, 0.11% L3/CR1), 8.54% long tandem repeat elements (LTRs) (3.40 ERV-L, 3.52% ERVL-MaLRs, 1.04% ERV Class I, 0.04 ERV Class II) and 2.53% DNA elements, together constituting a total of 33.40% interspersed repeats (Supplementary Table 1). An additional 2.2% of sequence consisted of small RNA (0.03%), satellites (0.64%), simple repeats (1.21%), and low complexity DNA (0.32%). In comparison, the contigs assembled from GRCh38 unmapped reads contained substantially more SINEs (8.5%), satellites (35.5%) and simple repeats (17%), but a lower proportional composition (but greater absolute composition) of LINEs (7.5%), LTR elements (3.56%) and DNA elements (0.9%).

We further evaluated the impact of the complete T2T-CHM13 reference by performing “cross alignment”, where the contigs assembled from T2T-CHM13 unaligned reads were aligned against those assembled from GRCh38 unaligned reads (Supplementary Figure 3). We observed that the T2T-CHM13 contigs were largely a subset of the GRCh38 contigs, with 70% having a high identity, high mapping quality alignment against a GRCh38 contig. Conversely, 75% GRCh38 contigs did not align well to the T2T-CHM13 set, but 70% aligned well to the T2T-CHM13 reference itself, reflecting its completeness and mirroring results from the African Pangenome study 21. Our results suggest that use of GRCh38 obscures genetic variation, particularly within erroneously assembled or unresolved regions that were later resolved in T2T-CHM13. For these reasons, unless otherwise noted, we focused our subsequent analyses on the contigs assembled from T2T-CHM13-unaligned reads.

Alignment and assembly of data from the SGDP SAS set
We repeated the alignment and assembly steps with data from 39 SAS individuals across 19 populations from the SGDP dataset. For 37 out of the 39 individuals, assembly of their unaligned reads against T2T-CHM13 (0.9 - 1.7% of the total input reads) and filtering to >1 kbp non-contaminant contigs resulted in 1.4 - 1.8 Mbp of assembled sequence per individual. Each assembly contained ~500 contigs with average N50s of 3.2 kbp, and the largest assembled sequences reached ~40 kbp (Supplementary Figure 4). We attribute the larger span of non-contaminate contigs to the higher sequencing depth of SGDP (32-66x, average coverage of 43x) compared to 1KGP (~30x).

Two outlier individuals from the Kusunda group exhibited markedly different alignment and assembly results. Specifically, ~10% fewer reads from these samples aligned to the reference and the unaligned reads assembled into >200x as much sequence as other SGDP samples (Supplementary Figure 5). BLAST and Centrifuge did not categorize these sequences as known human sequence or contaminants, and most quality metrics were in line with other SGDP SAS samples. Full details of these two outliers are provided in Supplementary Figure 6, but we hypothesize that the excess sequences originate from a mixture of unknown microbial contaminants and sequencing artifacts, and we therefore excluded these two samples from subsequent analyses.

Pangenomic improvements in alignment and assembly
We repeated the assembly process with reads that failed to align to the HPRC Minigraph and Minigraph-Cactus v1.0 pangenome references 3. These HPRC pangenome references were constructed by augmenting the T2T-CHM13 reference with the genome sequences from 47 individuals, including one individual from the South Asian continental group 3. Against the Minigraph reference, alignment rates rose by 0.3-1.1% per individual, which in some cases halved the number of unaligned reads compared to results using T2T-CHM13. The expanded Minigraph-Cactus reference further improved the alignment rates, though this improvement was marginal (<0.1%). For the SAS individual also present in the Minigraph and Minigraph-Cactus reference, we observed near perfect alignment across the genome, with the overall alignment rate rising from 98.7% against T2T-CHM13 to 99.8% against the Minigraph reference.

Despite the sizable declines in the number of unaligned reads, the reductions in the amounts of assembled sequence were not as dramatic (Figure 2d-f). Using reads that failed to align to the Minigraph pangenome reference, we assembled 450-500 kbp of sequence per 1KGP SAS individual in >1 kbp non-contaminant contigs, a ~10% decrease in sequence compared to the same individuals using T2T-CHM13 (Welch’s t-test, two-tailed p-value of 2.03 x 10-13; Figure 2e-f). This proportionally modest decline suggests that many of the newly aligned reads would have assembled into short contigs discarded during our filtering step. The number of >1 kbp non-contaminant contigs declined even less (Figure 2d), suggesting parts of previously assembled contigs had been resolved. The average N50 of the filtered assemblies was 1.9 kbp and the largest assembled contigs per individual were 15-18 kbp (maximum ~50 kbp). The augmented Minigraph-Cactus reference resulted in further modest declines in the amount of assembled sequence (<5%) and number of contigs that composed this sequence (<3.5%).

Comparison to non-SAS populations
To contextualize these results, we compared to analogous results obtained from samples from non-SAS populations. Specifically, we repeated the same alignment and assembly steps using 10 randomly selected individuals, 5 XXA and 5 XY, from each of the 21 non-SAS populations in the 1KGP dataset. Across the two linear references, the results with the non-SAS set were similar to those from the SAS set, with approximately three times as much sequence assembled per individual from reads that failed to align to GRCh38 versus reads that failed to align to T2T-CHM13. However, individuals from the SAS set produced on average 10-15% less assembled sequence than individuals from the seven African (AFR) populations (two-tailed Welch's t-test; t = -36.8; two-tailed p-value = 1.103 x 10-46), and on average 5% more sequence than the individuals from the non-AFR, non-SAS populations ((two-tailed Welch's t-test; t = 6.22; two-tailed p-value = 5.39 x 10-9) (Supplementary Figure 7, Supplementary Table 2). Performing population level comparisons, where the sequence assembled across 601 SAS individuals are compared just to the 10 individuals from each non-SAS population, we see all seven AFR populations have a two-tailed p-value < 1.76 x 10-7, and the non-SAS populations CHS, CDX, PEL, GBR, FIN and CEU have a two-tailed p-value < 0.05 (Supplementary Table 2).

It is worth noting that while T2T-CHM13 exhibits the greatest genetic similarity to 1KGP samples from the European (EUR) continental group, the amount of sequence assembled for samples from EUR populations is nearly the same as any other non-AFR population. This underscores the fact that human genetic diversity tends to be either globally common (present in all populations at relatively similar frequencies) or rare (geographically restricted, but only because it is rare and is also rare within a given population) 40. More broadly, these results are consistent with previous knowledge that 1) genetic diversity in African populations exceeds that in non-African populations as a result of the out-of-Africa bottleneck 41,42 and 2) there is an abundance of rare variation in human genomes due to rapid population growth 43.

Overall, these results highlight the improvements offered by pangenomes, which recover common variation that previously failed to map. However, substantial amounts of sequence can be assembled from reads that fail to align even to these pangenome references, highlighting the need for alternative approaches beyond reference-based alignment.

Placement of assembled contigs
Placing assembled contigs back into the reference genome may be possible through direct alignment of the contigs to the reference genome or using data from mapped mate pair reads 21. Using the direct alignment approach, approximately 40% of contigs aligned at least partially to the T2T-CHM13 reference, declining to ~8% upon required alignment lengths of at least 500 bp. This equated to an average of ~350 contigs that remained unaligned per individual, and we hypothesized that these may contain any non-reference sequence. In contrast, the mate pair approach placed fewer contigs than direct alignment (3,952 / 1.9%), but addressed cases where contigs diverged from the reference sequence. Combining these two methods, 19,691 of 199,564 contigs (9.8%) >1 kbp in length were placed against T2T-CHM13, with enrichments noted on chromosomes 6, 9, and Y (Supplementary Figure 8), potentially driven by their higher repeat content (e.g., the HLA locus on chr6). Of the placed contigs, and restricting to the 21 SAS individuals with available long-read data, >85% could be aligned with high identity against a personalized long read assemblies. While extending across all chromosomes, these placements were not evenly distributed throughout the genome, with enrichments noted on chromosomes 6, 9, and Y (Supplementary Figure 8), potentially driven by their higher repeat content (e.g., the HLA locus on chr6).

Intersection with annotated elements
We intersected the placed contigs with the “JHU RefSeqv110 + Liftoff v5.1” gene annotations from the T2T Consortium using bedtools 44. Across the SAS set, we identified 52,727 intersections with known transcripts, 12,021 with known genes, 5,050 with known exons, and 3,015 with known CDS regions (Figure 3a), with substantial heterogeneity noted across the genome (Figure 3b; Supplementary Figure 9). These intersections affected a total of 1,189 unique transcripts, 172 genes (106 protein-coding), 33 exons and 6 CDS regions (Figure 3c). While the majority of the 106 protein coding genes intersected with contigs from one or few individuals, some reflected more common variation (Supplementary Figure 10). Qualitatively similar results were observed when applying the same approach to samples from non-SAS populations (Figure 3d). Gene intersections in the SAS set included 11 genes where mutations have been related to eye disorders (AGBL1, CPAMD8, CRB1, HMCN1, KIF14, LIN9, MYO7A, PPEF1, TTC39C, USH2A and WDR17), as well as genes involved thyroid function (TSHR), facial dysmorphia (GPC3, CPA6, OPHN1, OR51A4 and SNP1), ciliary function (CC2D2A, DRC1, NME8, TTC39C), and cancers (DLEC1, DMBT1, LIN9, PTPRO, SEM1, SYK and ZNF292 (Supplementary Data 1; Supplementary Data 2-3). The observed count of 12,021 gene intersections was substantially lower than the number of intersections obtained through random placement of contigs against the reference genome (14,500 intersections, averaged over three runs), consistent with the action of purifying selection against common structural variation in genic versus non-genic regions.

Linkage disequilibrium within the shared contig set
The presence/absence of each shared contig across all 640 SAS individuals allowed us to measure linkage disequilibrium (LD) between pairs of contigs and identify “blocks” of placed contigs that likely belong to the same larger sequence. Specifically, we computed the squared Pearson correlation coefficient (r2) between all pairs of 13,875 shared contigs and found 644,578 pairs of contigs with r2 > 0.9 (0.33% of all pairs). Of these, 194,428 (0.1% of all pairs, 30% of pairs with r2 > 0.9) are between contigs on the same chromosome and 56,935 are between contigs within 10 kbp of each other. We also observed 4,142 pairs where a placed contig is in strong LD with an unplaced contig (r2 > 0.9 across the 640 individuals). In these cases, the placed contig can act as an “anchor”, allowing us to narrow down the potential location of the unplaced contig within the genome. Within these anchor contigs, we identified 3,334 intersections with 24 genes across 13 chromosomes (Figure 3e-f), as well as 58 intersections with known GWAS hits.

Analysis of unplaced sequence
Despite placing approximately 20,000 contigs from the SAS set against T2T-CHM13, we assembled over 180,000 contigs that could not be placed. By querying this set of unplaced contigs using BLAST, we found that >90% had a high similarity match with non-reference human sequence or reference and non-reference primate sequences (Supplementary Data 4), supporting their non-contaminant human origin. We also observed a small fraction of contigs with high similarity to non-primate mammalian sequences, potentially reflecting highly conserved loci.

Analyzing pangenome placements
The limited availability of pangenome tools precluded extension of the placement approach to the graph pangenomes. As an alternative, we performed conventional alignment to the linear references that comprise the pangenome, identifying 8.5% of the assembled contigs (16,854/199,564) with an alignment length >500 bp, thus slightly exceeding the number of contig alignments to T2T-CHM13 alone (15,739).

Of the contigs that failed to be placed against T2T-CHM13, 905 of 179,873 (0.5%) aligned to the HPRC Minigraph pangenome. These alignments occur throughout the genome, with no specific chromosome or region especially enriched for these placements. These results highlight the power of pangenome approaches while emphasizing the need for centralized annotations and auxiliary data to facilitate downstream interpretation.

Validation with long reads
The initial data release from the 1KGP ONT effort 45 contains long read data from 21 individuals present in our SAS set, which can be used to help validate the assembled contigs from these individuals. We constructed simple assemblies for these 21 individuals with their long reads using the Flye assembler 46. The assemblies contain 2.81 to 2.84 Gb of sequence, with N50s ranging from 31.9 to 49.2 Mb and mean coverage of 29 to 47x. Although these assemblies are far from telomere-to-telomere completeness, we observed that 40-50% of these previously unaligned reads aligned well to the personalized long read assemblies of the respective samples. The “rescue” of these unaligned reads highlights the limitations of a single reference genome and confirms that much of the missing sequence is variation from loci that are easily assembled from long-reads.

For the same 21 samples, we also aligned the contigs we previously assembled from reads that failed to align to T2T-CHM13 to the corresponding long-read assemblies, finding that 85% of contigs aligned with high identity (Figure 4a). These included 92.7% (331/357) placed contigs that intersect with genes (Figure 4b). We hypothesize that the remaining 15% of contigs originate from more complex regions that may require tailored assembly approaches or additional sequencing with complementary platforms.

We repeat this long read validation using a second dataset from Schloissnig et al. 47, and perform the same analysis on 20 SAS individuals from this set. The Flye assemblies for these individuals are comparable to those from the 1KGP set (2.65 to 2.82 Gb of sequence, N50s of 25.6 to 45.7 Mb, and mean coverage of 14 to 33x), and initial alignment sees 40-50% of previously unaligned reads now align to the personalized long read assembly. Contig alignment is slightly lower with this data (80-87%), but in line with the fact this dataset is substantially lower coverage than the 1KGP ONT set.

Shared sequence analysis
Summing over all SAS samples, we assembled 410 Mbp of sequence in > 1 kbp, non-contaminant contigs from short reads that failed to align to T2T-CHM13. However, we anticipated that much of this sequence would be redundant, as multiple individuals may share an allele (identical by descent, or more rarely, by recurrent mutation) that is divergent from the reference (Supplementary Figure 11). To address such cases, we performed an all-vs-all alignment of the assembled contigs using Minimap2 48, and considered contigs with >90% alignment to have originated from a shared allele. In total, this procedure collapsed the 199,564 contigs totalling 410 Mbp of sequence into 13,875 shared contigs and 2,137 unique (i.e., “non-redundant”) contigs totalling 47 Mbp of sequence. The allele frequency spectrum demonstrates that a large fraction of the contigs are rare (present in <10% of the combined SAS set), with the majority of these present in just 1 or 2 individuals, consistent with the known abundance of rare variation in human populations. (Figure 4c, Supplementary Figure 12-13). We note that the presence of trios distorts the frequency spectrum compared to expectations for an unrelated sample, but the practical effect is minimal in practice (e.g., enriching for doubletons versus singletons). Further analysis of the shared sequences between members of the trios can be found in the Supplementary Notes (“Patterns of inheritance and similarity across contigs assembled in 85 1KGP trios”) and Supplementary Methods.

We were also interested to quantify the amount of sequence from these SAS samples that is present in genomes of individuals from non-SAS populations. To quantify this proportion, we aligned read data from the selected 210 non-SAS individuals against the non-redundant SAS contigs. We found that, on average, 3.8 - 4.5% of reads from each individual align to this set of contigs, with individuals of all 21 non-SAS populations exhibiting roughly the same alignment rate. This observation further highlights the fact that most variation is rare, and that the variation that is common tends to be shared across all populations.

Comparison to traditional variant finding approaches
To evaluate the strengths and limitations of our "rescued read assembly" approach, we compared our results to those obtained from conventional short-read-based SV callers, starting with Manta 49 and LUMPY 50 and focusing on insertions (given that these are also the focus of our assembly strategy). Manta identified an average of 2,300 insertions per individual, but only an average of nine insertions per individual were longer than 500 bp (192 in total, with a maximum of 16 found in HG03898), and no insertions were longer than 1 kbp. In contrast, our approach assembled an average of 299 (SD = 9.3) contigs longer than 1 kbp per individual, an average of 36 of which could be placed against T2T-CHM13. The placed contigs overlapped some of the smaller and many of the larger (>500 bp) insertions called by Manta but also included hundreds of placed sequences that were not detected by Manta (Supplementary Figure 14-15), 85% of which were validated by long reads. LUMPY is largely focused on deletion calling and is limited in its ability to detect small insertions, but is able to detect large insertions through identification of breakends 50. Across the 21 SAS individuals, LUMPY detected none of these large insertions, despite our approach identifying and placing several large contigs in these samples, further validated with long reads.

We also compared the sequence assembled by our approach to the sequences found by PopIns2 51, which detects large non-reference sequences by assembling reads that do not have a high quality reference alignment. Across a random selection of 10 1KGP individuals (one XX and one XY from each of the five SAS 1KGP populations), PopIns2 assembled an average of 6,264 contigs >1 kbp per individual from reads that fail to align to T2T-CHM13, in broad accordance with our results (Supplementary Notes, Supplementary Figure 16). Comparing the placement of these contigs against the contigs assembled and placed by our approach, we found that of the 267 contigs placed against these 10 individuals, 210 overlapped a placed Popins2 contig. However, our approach assembled more large contigs than PopIns2 (average of 510 contigs >5 kbp and 95 contigs > 10 kbp by our approach vs 432 contigs >5 kbp and 72 contigs >10 kbp by PopIns2) and also places a larger fraction of such contigs (10-12% of our contigs vs 8% of PopIns2 contigs) based on alignment, LD, and mapped mate pair reads.

RNA-seq analysis of hidden variation
We used RNA-seq data from the MAGE project 52, obtained from lymphoblastoid cell lines from 731 globally diverse individuals from 1KGP, to evaluate the functional potential of the previously assembled contigs. We specifically focused on the 140 MAGE samples from SAS populations that overlap with our study. We aligned RNA-seq reads from each individual to their contigs that were assembled from reads that failed to align to T2T-CHM13. Per individual, ~800-1,000 reads aligned to their respective contigs, with the majority (>80%) spread across few contigs. While this represents a small fraction of the input RNA-seq dataset (<0.1%), this is not unexpected, as the total length of assembled contigs per individual is <0.02% the size of the human genome. Across the 140 individuals, 90% of contigs with >100 RNA-seq alignments had BLAST hits to non-reference human and primate sequences, some of which is annotated as functional (Supplementary Data 5-6). These include hits to annotated non-reference sequences from chromosomes 1, 2, 3, 7, 9, 15, X and Y.

We also aligned the RNA-Seq data from all 140 individuals against a combined set of filtered contigs from all 601 1KGP SAS individuals and identified 200 contigs with >100 read alignments and 19 contigs with >1,000 read alignments (Figure 4e-f, Supplementary Figure 17-19), 83% of which had BLAST hits to non-reference human and primate sequences (Supplementary Data 5-6). Across both individual-specific and combined RNA-seq alignments, we found that the majority (>80%) of highly aligned-to contigs appear in just one or a few individuals, with few classified as widely shared sequences (Figure 4d).

Intersecting the placement information for these contigs with the RNA-seq alignments, we find that the majority of highly aligned-to contigs could be placed via alignment to T2T-CHM13 (140/200 of contigs with >100 alignments, 12/19 contigs with >10 alignments), with most of the placements tracing to the Y chromosome. Aligning these contigs against the HPRC Minigraph pangenome reference modestly increased the rate of alignment (144/200 of contigs with >100 alignments, 13/19 contigs with >10 alignments), with few additional contig alignments to chromosomes 19 and Y.

Intersection of assembled contigs with GWAS hits
We next explored the phenotypic relevance of the assembled contigs by comparing the placed sequences against existing catalogs of trait-associated loci. We intersected the set of placed contigs against the GWAS v1.0 catalog (lifted over to T2T-CHM13 coordinates 4,5) and identified 2,812 cases where a placed contig overlapped an annotated GWAS hit. These hits were largely concentrated on chromosome 6 (1,995), Y (601), and 5 (208), with small counts noted on other autosomes (Figure 5b). Loosening our search criteria, we identified >6,000 contigs within 100 bp of 36 unique GWAS hits and >15,000 contigs within 1 kbp of 112 unique GWAS hits (Figure 5a-b), the vast majority of which were validated with long-read data (<100 bp distance: 215/254, 84.6%; <1 kbp distance: 506/587, 86.2%; Supplementary Figure 20).

Finally, we investigated if any of the placed contigs overlapped with or were near trait associated loci from the UK Biobank (UKB) 53,54 and East London Genes & Health (ELGH) study 55 (Figure 5c). The former contains data from ~400,000 individuals of all ancestries, while the latter contains data from ~37,000 SAS individuals in the UK. Across the 21 biomarker traits from the Pan-UK Biobank 56 and 42 biomarker traits from the the Genes and Health (ELGH) 55 dataset, we identified only 5 placed contigs within 10 kbp of significant locus at a relaxed p-value threshold of 1 × 10-6. Expanding to a 1 Mbp window surrounding GWAS hits, we identified 8-14 traits close to placed contigs as we varied the p-value threshold from 1 × 10-8 to 1 × 10-6 respectively (Supplementary Data 5). However, we did not find any placements within 1 Mbp of a GWAS hit within the smaller UK Biobank “CSA” (Central/South Asian) ancestry group (n = ~7,000) until we lowered our p-value threshold to 1 × 10-5. This is not unexpected - especially as ~90% of our assembled contigs remain unplaced - and does not rule out the functional potential of these sequences. Understanding the functional implications of these complex loci will require future association studies based on genotype data encompassing these technically challenging loci 57.

Methods
Detailed guides for several of these sections, containing the exact syntax used, can be found on the project GitHub page, inside the “Guides” subdirectory.

Ethics Statement
This study uses publicly accessible, entirely anonymized human sequencing and GWAS data, and does not generate any new data from new human participants. As such, ethical approval was not required.

Data
We used read data from the 1000 Genomes project (1KGP) 27 and the Simons Genome Diversity Project (SGDP) 28, totaling 640 individuals from 24 South Asian populations across the two sets. The read data from these sets is high coverage (>30x), high quality, short paired reads. The 601 SAS individuals present in the 1KGP set are from five South Asian populations. Three are in the South Asian diaspora - Gujarati in Houston (GIH), Sri Lanka Tamil in the UK (STU) and Indian Telegu in the UK (ITU) - and two are in South Asia itself - Punjabi in Lahore (PJL) and Bengali in Bangladesh (BEB). The 601 individuals are largely evenly distributed between the five populations, with 109-144 individuals from each. In the 1KGP set, there a re 282 XX individuals and 319 XY individuals. The SGDP set contributes 39 individuals, spread across 19 different groups within India, Pakistan, Nepal and Bangladesh. Most of these are in pairs from each group, with a single individual from one population (Khonda-Dora from India) and four from another (Punjabi from Pakistan) the only exceptions. In the SGDP set, there are 8 XX individuals and 31 XY individuals. All sex and demographic information is from the 1KGP and SGDP datasets, both of which state the information was self-reported by participants during the data collection process.

These reads are aligned against two linear references and two pangenome builds. The two linear references are GRCh38 35 and T2T-CHM13 v2.0 4,20. We used the Minigraph v1.0 and Minigraph-Cactus v1.0 builds of the recently released Human Pangenome Reference Consortium (HPRC) draft pangenome references. Both pangenome references are constructed over the T2T-CHM13 reference and include variation present in 47 selected individuals from the 1KGP set, one of which is a SAS individual present in our set. The Minigraph reference contains only SVs larger than >50bp, while the Minigraph-Cactus reference is further augmented to include base-level alignments for other variants.

We performed validation of our assembled contigs using long read data from the 1KGP ONT effort 34, which is available for 21 SAS individuals in our 1KGP set, and using long read data for 20 SAS individuals from the Schloissnig et al. cohort 47. Our functional analysis was augmented using RNA-Seq data from the MAGE dataset 52, which provides RNA-seq data for 140 SAS individuals. For our GWAS and Biobank analysis, we used the standard T2T-CHM13 GWAS catalog and the phenotype manifests from the Pan-UK Biobank 56 and the East London Genes and Health effort 55, focusing on the quantitative biomarker traits in each.

More details about the data characteristics of the read data and its alignments to the linear references can be found in Supplementary Figures 21-24.

Read alignment and extraction
Reads from the 640 SAS individuals in the 1KGP and SGDP sets were aligned to GRCh38 and T2T-CHM13 using bowtie2 v2.4.1 66, and the unaligned and poorly aligned (quality value < 20) reads were extracted using samtools v1.14 67. Unaligned reads of two forms are extracted from these alignments: both reads where the read and its pair are both unaligned, as well as single unaligned reads with a mapped mate read.

Assembling unaligned reads
We tested three assemblers on a selected group of 10 SAS individuals (two from each of the 1KGP populations): SPAdes v3.15.3 68, MaSuRCA v4.0.9 69 and MEGAHIT v1.2.9 36. There was some variation in the number of short contigs (SPADes produces far more than the other two), but the number of large non-contaminant contigs varies by no more than 5% between assemblers. All assemblers were tested using their default parameters, with the four sets of unaligned reads extracted in the previous step passed in as input. The biggest difference between these assemblers is the variation in runtime. We opted for the fastest of these three assemblers - MEGAHIT - in an effort to speed up analysis and reduce costs.

Filtering and processing assembled sequence
As in the APG analysis, we first discarded all contigs less than 1 kbp in size. This was done by filtering the assembly FASTA file for sequences longer than the target length, with all other sequences being discarded.

We also screened these sequences for contaminants using two tools: Centrifuge v1.0.3 38 and BLAST 37. For Centrifuge, all assembled sequences were compared against the “Bacteria, Archaea, Viruses, Human” database, and sequences categorized as non-human are removed. With BLAST, we used the “refseq_genomic” database, and removed all sequences that have a high scoring match to a non-human sequence. The contaminant sequences can be identified by parsing the resulting outputs, and these sequences were then cut from the FASTA file. After the filtering stage, we were left with a list of >1 kbp, non-contaminant contigs from each of the SAS individuals for use in the later stages.

Repeat content in the assembled contigs and reference genomes is detected, masked and annotated using RepeatMasker v4.1.2-p1 39. Using the assembled contigs as input, RepeatMasker detects and masks SINEs, LINEs, LTR and DNA elements, small RNA and satellites, and depending on the use of the “nolow” parameter, simple repeats and low complexity elements.

Assembly of non-SAS sequences
We selected 10 individuals - 5 XX and 5 XY - from each of the 21 non-SAS 1KGP populations.This was done by iterating through a list of participants sorted by participant ID, and choosing the first 5 XX and 5 XY from each population. As done with the SAS set, these individuals’ reads were then aligned against T2T-CHM13 using bowtie2.4.1, the unaligned and poorly aligned reads were extracted, the reads were assembled using MEGAHIT with default parameters, and filtered for size and contaminants.

Reference vs reference comparison
Another analysis we performed to highlight the improvements made in T2T-CHM13 over the previous GRCh38 reference is “cross-alignment”. Contigs assembled from unaligned reads against T2T-CHM13 were aligned against those assembled from unaligned reads against GRCh38, and vice versa. We also aligned the unaligned contigs from both references against both references, to demonstrate how T2T-CHM13 is more effective at “rescuing” these contigs. These alignments were done using Minimap2 v2.22 48, and the number of contigs that align or fail to align from and against each reference is noted.

Categorizing sequence as shared or unique
We evaluated how much of the sequence in the T2T-CHM13 unaligned contigs was unique to the individual it is present in and how much was present in multiple individuals using all-vs-all alignments of the contigs to themselves.

The all-vs-all alignment was performed using Minimap2 v2.22 48 with the strictest “Long assembly to reference mapping” default parameters (asm5), and the resulting output file was used to collapse the shared sequence. To do this, we repeatedly iterate through the output of the all-vs-all alignment and consider contigs with alignment lengths to each other greater than 90% of their own length to be the same sequence. We selected the larger of these sequences to be the representative “super-contig”, and we noted the other sequence as being equivalent to or is derived from this one. We repeated this process until a full iteration through the set occurs without any more sequences being collapsed down, leaving a finalized list of “super-contigs”, a list of which contigs have collapsed into each super-contig, and a list of sequences that have undergone no collapsing and are therefore unique. Combining the list of super-contigs with all non-collapsed sequences gave us our finalized, “non-redundant” set of SAS contigs.

This approach acts as a floor on the amount of shared sequence, as it only merges contigs that are almost entirely identical. However, this approach does not consider contigs that share small subsequences or overlap only slightly at the ends. These situations would be best captured through all-vs-all comparisons and collapsing with no threshold for minimum alignment length. Alternatively, using a similarity metric such as Mash distance 70 could quickly estimate the similarity between contigs and decide which ones to merge. Such an approach would account for the case where two contigs have smaller regions of similarity dispersed throughout their length. For Mash distance calculations we used Mash v2.3 70 with default parameters (1000 hashes, sketch and compute overlap in one step) on the Cartesian product of the set of individuals and themselves.

Assessing non-SAS contribution to SAS contigs
For each of the 210 non-SAS individuals, consisting of 5 XX and 5 XY individuals from the 21 non-SAS 1KGP populations, we aligned their reads using bowtie2 v2.4.1 to the finalized non-redundant set of SAS contigs. We then parsed the alignment output logs to extract their alignment rates, and then further analyzed these counts to generate sub- and super-population summaries.

Placing assembled sequence
In order to place the contigs against the reference, we followed two approaches: direct alignment of the contigs using Minimap2, or using mate-pair reads that have been mapped to the reference.

For direct alignment, we aligned the contigs against the chosen reference, and reported the locations of all contigs with alignments over >500bp. This alignment and filtering stage was performed using Minimap2 and then a simple post-alignment filtering based on the length of the alignment, which is stored in one column of the output .paf file.

The mate-pair placement method is the approach taken in the construction of the African Pangenome 21. During our read extraction step, we explicitly extracted unaligned reads that had mapped mates, and once we knew which contig an unaligned read was assembled into, we could use the mapped mate to anchor the entire contig.

Intersections with and analysis of annotated elements
We used the “JHU RefSeqv110 + Liftoff v5.1” annotations generated by the T2T consortium and intersected this with the list of placements generated in the previous step. We computed these intersections using bedtools v2.30.0 44, using the “-wb” parameter to output the details of the intersection annotated element. The bedtools output could then be further analyzed using Python scripts to generate counts for the number of placements per gene, the number of times a gene is intersected with, and other such counts.

We obtained details on the gene types, status, summaries and expressions from their NCBI pages. This can be done using the NCBI E-utilities or simply by scraping the corresponding webpage for each affected gene - we primarily used the latter, due to technical issues with the E-utilities tool. The vast majority of the summary and expression data come from RefSeq and “HPA RNA-seq normal tissues” project 71.

GWAS intersections with placed contigs
We used the same steps as the annotated region intersections, except this time using the “GWAS v1.0” annotations generated by the T2T consortium. We used “bedtools intersect” with the same configuration to generate the list of all direct overlaps between GWAS elements and placed contigs.

In order to account for placements that are close to but not strictly overlapping the GWAS elements, we implemented a Python script to compare the coordinates of placed contigs against those of known GWAS elements, and compute the distance between all placed contigs in a chromosome and all its GWAS sites. This was done by computing the distance between all GWAS sites and the start and end points of each contig, thus allowing us to both identify all overlaps but also pre-compute all distances between all pairs of contigs and GWAS sites. We could then alter the threshold we want to test for, and count the number of interactions with a distance lower than this threshold.

Analysis of unplaced sequence
We queried the unplaced sequences against the BLAST “nt” database to ascertain their potential function, and extracted the top 50 highest scoring hits for each contig. We did this by parsing the BLAST output to extract the lines containing the top matches after each query sequence, before tallying up the number of sequences included in those “top hits”, and outputting these to a table.

Linkage disequilibrium of shared sequence
From the 13,875 sequences present in multiple individuals, we generated a 640 x 13,875 binary grid representing the presence (1) or absence (0) of each of the shared sequences (or those that have been collapsed down into them) in each of the 640 SAS individuals. Each column in the grid represents the presence or absence of a given contig in each of the 640 individuals in the SAS set, and two contigs that are in LD with each other will have highly similar columns.

We computed a pairwise r2 between each of the 13,875 contigs, and then filtered out pairs that had a score higher than our chosen threshold. We then used our placement data to identify pairs of contigs placed within the same chromosome or within a threshold distance, and to identify pairs where one contig is placed and the other is not.

This approach “links” pairs of contigs with high degrees of association; it may be possible to extend this to tie together and place larger groups of contigs or incorporate other data into the association 33.

Pangenome analysis
We tested a range of pangenome aligners, including GraphAligner 72, VG 73, vg giraffe 74, Minigraph 75 and Minigraph-Cactus 76. We found the vg suites of tools to be easiest to use, but the outputs across the aligners are largely consistent. Regardless of tool, we either forced the output of the aligners to be in the GAM format, or used existing conversion tools to create GAM files. These GAM files can then be converted to the standard SAM/BAM/CRAM format for read extraction. Once in the SAM format we continued to filter and extract reads exactly as we did against the linear references, and repeat the same assembly and filtering steps for these contigs.

Placement was also more complicated due to there being multiple parallel variant sequences, differences in coordinate systems and indexing between aligners (or even different versions of the same aligner), and the increased computational complexity of now performing alignment against a branching sequence graph. We obtained basic coordinates and counts from parsing the various alignment files generated by the different aligners, all of which are run with the finalized list of contigs and the pangenome reference as input and default parameters.

Once placed, there are few resources (GWAS catalogs, annotations, identified significant loci) that have been fully translated into the coordinates and formats used in these pangenome builds. Therefore, we did not replicate this analysis against the pangenome references.

Long read validation
We used long read data for 21 SAS individuals from the 1KGP ONT and 20 SAS individuals from Schloissnig et al. in an effort to validate our contigs 34. We first construct basic long read assemblies for each of the 21 individuals. This was computed using the Flye assembler 46, with default Flye v2.8.1 parameters for “raw-ONT” reads, and with just the long read sets as input.

The previously unaligned short reads were then aligned against these assemblies using bowtie2, and the assembled unaligned read contigs can be aligned using Minimap2. We used the same tool versions (bowtie2 v2.4.1, and Minimap2 v2.22) and parameters for each of these as the read alignment and contig mapping stages in the analysis of the short read linear reference data.

For short reads, we filtered the resulting SAM alignment file to obtain a list of the reads that now align as well as a list of the reads that remain unaligned. Similarly, the output of the Minimap2 alignment of contigs to this individual-specific long read reference was filtered based on the alignment quality and length.

Long read validation of contig placements
We overlapped the list of contigs with high quality alignments in each of the 21 individuals against the list of placed contigs and their intersections with gene and GWAS elements, giving us a list of long read validated intersections. This was done by iterating through the list of placements that intersect an annotated element or GWAS site, and checking if the placed contig had an alignment to the long read assembly.

RNA-Seq analysis
We used STAR v2.7.10 77 to align the RNA-Seq data obtained from the MAGE dataset 52 from each of the 140 individuals against the short read contigs assembled from their unaligned reads against T2T-CHM13, and selected all contigs with large numbers of RNA-Seq alignments.

We tested both single and two-pass alignments, and found little to no changes between these two approaches. Alignment was also tested without splice junctions or annotations, and with a list of known human annotations; once again there was no substantial difference. The remaining numerical and memory parameters (including “genomeSAindexNbases” and “limitGenomeGenerateRAM”) were left as default, or amended to STAR’s recommendation.

Once we had identified the most aligned-to contigs, we used the BLAST “nt” database to query for these sequences and report all top hits. As we did for the analysis of unplaced sequence, we parsed the BLAST output for each query sequence and output the top 50 highest scoring hits for each contig.

We also performed an exhaustive all-vs-all alignment of all the RNA-Seq reads across the 140 individuals against all the contigs assembled across the 640 SAS individuals. The parameters for this alignment were slightly adjusted, mainly so that STAR can handle the sharp increase in the number of reads (notably, the “genomeSAindexNbases” and “limitGenomeGenerateRAM” are adjusted as per STAR’s requirements). Similar to the previous analysis, we extracted the most aligned to contigs (>200 RNA-Seq alignments per contig), and ran them through the same BLAST query process. The top 50 BLAST hits (if there are that many) for each contig were then extracted.

Biobank analysis
We used the UK Biobank Phenotype Manifest from the Pan-UK Biobank effort 56 to obtain our chosen set of 21 quantitative biomarker traits. Using the provided data for each trait, we were able to extract a list of all chosen loci in the genome and their potential association with each of the traits. The manifest represents this association as the negative log of the p-value, and this p-value is computed for both the entire UK Biobank set and for each of the individual ancestry groups within the Biobank. For our analysis, we used the p-values across the “CSA” ancestry group and the entire UK Biobank set. The Genes and Health dataset 55 is formatted in a similar way, and contains the associations between all chosen loci in the genome and 42 quantitative biomarker traits.

For each trait, we then “swept” across all potential loci, and identified all loci that have significance above our chosen p-value threshold. For all significant loci within 50 kbp of each other we selected the highest scoring loci (i.e. the one with the most significant p-value) as the representative “peak” for the entire group. We computed the pairwise distance between all selected peaks for a given trait and all placed contigs across the SAS set, and identified all peak/placed contig pairs that are within 1 Mbp of each other. This process is repeated for all 21 biomarker traits in the UK Biobank, and for the 42 quantitative biomarker traits in the Genes and Health dataset.

Comparison to existing variant discovery tools
We compared the assembled contigs from 21 SAS 1KGP individuals with LR data to the variant calls made in those individuals using Manta 49 and Lumpy 50. Both tools were run with default parameters using the output of the alignment stage. Manta uses the .cram file from the alignment stage directly, while there is an initial preprocessing step to extract the split and discordant reads for Lumpy.

Both tools output a vcf file that detailed the structural variants they find, which we parsed to determine the location, type (“INS” for insertions, “BND” for breakends of large insertions in Lumpy) and size of the SVs. We then compared this against the list of contigs assembled from these individuals, the contig placements in this set, and which of the contigs were validated with long reads.

For PopIns2 51, we selected 10 samples from the 1KGP SAS set, and followed the steps outlined in the PopIns2 Github repository using default parameters, starting from the alignments of these individuals’ reads against T2T-CHM13. We then parsed the resulting assembly and vcf files to compare against the assembled and placed sequence from our approach.

Plotting and Statistical Analyses
The analyses in this study were computed using a range of established bioinformatics tools and custom scripts, and plotting was performed using existing Python libraries, such as sklearn 78 or matplotlib 79. Each existing tool used in this study has been mentioned in the relevant methods section.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Das and colleagues present a significant reanalysis of unmapped sequence reads from 640 individuals from both the\nThousand Genomes Project and the Simons Diversity Project to assess the ability to identify \"lost\" sequences failing to map\nto the Human Reference Genome. The authors validate these primarily short read sequences with long-read sequences\nfrom multiple individuals, RNA-seq data and BLAST analyses to characterize the extent of novel variation found in South\nAsian populations. The manuscript is exquisitely written and presented and will provide needed data on population variation\nand further expand upon the need for the human pangenome to include many more ethnicities. My only real issue with the\nmanuscript is the apparent neglect in comparing their analysis pipeline with other works that have previously utilized\nunmapped reads from species including zebrafish (Faber-Hammond and Brown 2016, [Zebrafish]), dog (Holden et al.,\n2018), human (Chekanov et al., 2010; Genovese et al., 2013; Liu et al, 2014; Faber-Hammond and Brown, 2016[Human\nGenetics]) and other. While the current project expanded upon the results of these researchers, using long-read sequences,\nthe T2T reference genome, and RNA-seq data to provide more in depth assessment of unaligned contigs, the finding of\nsignificant sequence data not being mapped to reference genomes is consistent.\nThe methodology, detail, interpretation, and conclusions made by the authors are all exceptional.\nI have no specific revisions as the manuscript was exceptionally well written and complete.", "reviewer_2": "The authors explored genetic variation often overlooked by conventional genome mapping methods by analyzing\nunmapped and poorly mapped reads from 640 South Asian individuals—601 from the 1000 Genomes Project and 39 from\nthe Simons Genome Diversity Project. They assembled tens of megabases of novel, non-redundant sequences and\ncompared them across global populations, revealing substantial hidden variation that remains undetected even with\ncomplete genome and pangenome references. Their analysis uncovered thousands of contigs intersecting protein-coding\ngenes and GWAS loci, and these findings were further supported using long-read sequencing and RNA-seq data. This work\nhighlights an innovative strategy for uncovering overlooked genetic diversity. Notably, the authors took an innovative\napproach to better characterize South Asian genetic variation and made all code publicly available, which is an important\ncontribution to transparency and reproducibility. However, some observed patterns in the data remain unexplained, and\nadditional analyses could further enhance the impact and interpretability of the study.\nMajors:\n1- The authors used all South Asian individuals in the most recent 1000G high-coverage release “601 1KGP SAS\nindividuals”. However, in this release, 85 trios (mother-father-children) from South Asia are included. The relatedness\nbetween these individuals can have an impact on different measures of the study: e.g. Linkage Disequilibrium or number of\nshared unmapped reads, as they will share more haplotypes than unrelated individuals. It’s something important to consider\nand discuss, as a quality check or as a measure of heritability of unmapped reads for example.\n2- In the PCA of unaligned reads, individuals cluster by dataset (SGDP vs. 1000G) rather than by population, which is\nunexpected given known patterns of genetic diversity and suggests the presence of technical batch effects. This raises\nconcerns about the direct comparability of the two datasets and may influence the interpretation of population-specific\nsignals. The authors mentioned the potential impact of sequencing depth on differences between datasets “We attribute the\nlarger span of non-contaminate contigs due to the higher sequencing depth of SGDP compared to 1KGP” but no specific\nvalues are provided in the main text or Methods. Including these details would help assess the extent to which sequencing\ndepth differences might explain the observed clustering, and help discussing their potential impact on the study’s\nconclusions.\n3-a. Still related to the PCA of unaligned reads, I noticed that the 1000 Genomes individuals appear to separate into two\ndistinct clusters, which do not correspond to population structure. This raises the possibility that other factors—such as the\nsex of the individuals—may be driving this separation, especially given the reported enrichment of assembled contigs on\nchromosome Y. Investigating this further could provide valuable insight and strengthen the study’s findings.\n-b. Additionally, the journal encourages to follow to the SAGER guidelines\n(https://researchintegrityjournal.biomedcentral.com/articles/10.1186/s41073-016-0007-6), which recommend conducting\nanalyses that can reveal sex-related differences. In this context, it would also be helpful to report the number of male and\nfemale individuals included in the study.\n4- The study focuses on South Asian populations, yet it does not engage with key aspects of South Asian evolutionary\nhistory—such as the region’s high genetic diversity between sub-populations, high patterns of homozygosity, or population\nstructure—which could provide important context for interpreting the findings. For example, do individuals who share the\nsame unmapped reads are also more closely related (e.g. share more Identical-by-Descent segments)? Exploring this could\nhelp clarify whether the shared sequences reflect deep population structure, recent relatedness, etc.\nAdditionally, while GenomeAsia 100K Project is briefly mentioned, other major regional efforts such as GenomeAsia Phase\n2 dataset (Wall et al. 2023) or the ongoing Genome India project (Bhattacharyya et al. 2025) are not cited. Including and\nsituating the work within the broader landscape of South Asian genomics would enrich the relevance and impact of the\nstudy.\n5- The authors state: “We found that, on average, 3.8–4.5% of reads from each individual align to this set of contigs, with\nindividuals of all 21 non-SAS populations exhibiting roughly the same alignment rate.” This is an intriguing result, and it\nwould be valuable to contextualize it—are these alignment rates consistent with other “non-redundant” set of population-\nspecific contigs? Also, is this a genome-wide pattern (or is it linked to one chromosome/one genomic region)?\n6- The authors report that “individuals in the SAS set had on average 10–15% less assembled sequence than the\nindividuals in the seven African (AFR) populations, and on average 5% more sequence than the individuals in the non-AFR,\nnon-SAS cohorts.” While these differences are potentially important, it would be helpful to know whether they are statistically\nsignificant. More broadly, the manuscript would benefit from the inclusion of formal statistical tests throughout to support key\nclaims and enhance the rigor of the findings.\nMinors:\n1-The figures are not cited in order within the text—for example, Figure 2b is referenced after Figure 2f. Additionally, some\nfigures such as Figure 1 (a and b), Figure 2a, and Figure 3d do not appear to be cited at all. Notably, Figure 2a compares\nalignment rates for a method that is only described in the Methods section; this figure should either be removed or properly\ncited in the main text alongside the corresponding method. Improving the order and completeness of figure citations will\nenhance the clarity and flow of the manuscript.\n2-Figure 1A would be clearer and more informative if it included the sample sizes for each group.\n3-Figures S6, S8, S9, S13, S18, S19, and S26 are difficult to read due to low resolution and/or small font size. Improving the\nclarity and readability of these figures would greatly enhance the quality and accessibility of the supplementary material.\n4-Page 6 line 16, the acronym “APG” is used without being defined.", "reviewer_3": "The authors present analysis of assembled contigs from South Asian individuals. The extra sequences were found to\noverlap coding sequences and GWAS hits. The methodology is straight forward and fairly complete, and the multiple well\nexecuted vignettes give functional context for what missing sequences were assembled. The LD placement of contigs is a\nclever trick (similar to work by Giulio Genovese https://pmc.ncbi.nlm.nih.gov/articles/PMC3971454/ - may be worth noting\nhere). The low overlap with Manta insertions is expected.\nSome gaps in analysis remain:\nComments.\nSignificant. Unless I’m missing something, the assemblies are not released with the manuscript. Because this is analysis of\npublic, open consent data, the data being released for review is a condition for acceptance. These should be uploaded to\nzenodo, and it will require an additional round of reviews as the reviewers examine the contents of the assembled contigs.\nSignificant. The number of contigs that are not annotated as having homology (placement) to the reference is high. Similar\nexperience comparing inserted sequences back to the genome in HPRC analysis showed that the majority of inserted\nsequences are repetitive elements. Running blast against NR is a great first step, but more analysis should be done\n(partially relating to the next comment).\nSignificant. The repeat content of the new contigs should be annotated. One would expect the majority of the new\nsequences to be derived from mobile elements and tandem repeat expansions.\nMinor. The GWAS hits should be compared to the number of hits that one would have randomly shuffling the aligned\ncontigs, taking into account repeat content of the contigs (e.g., similar repeat content of the locations of the shuffled\nalignments as original alignments).\nMinor. This is more of a commentary than a comment, but the data that were analyzed were all public. So why focus only on\nSouth Asian individuals and just 10 individuals from other populations? The assembly pipeline is clearly robust. The context\nof 1.8B individuals is important, but having additional populations (more than 10 samples) to gauge how much sequence is\nmissing from these individuals versus from other populations, for context.\nMinor. Lines 3-6, 10-14, page 15 offer qualitative descriptions of results (“most of…”). This should be changed to offer\nspecific percentages.\nMinor. The assembled sequences with RNA-seq alignments should be searched in the HPRC genomes.\nMinor. The shared sequence analysis should state the 5% divergence allowed, based on the commands in the supplement.\nThis is important for the context of how divergent sequences can be and considered shared.\nMinor. Page 11, gene intersections. If these are intronic, they are probably of low functional impact. If there is some\nadditional context to believe that eye-related disorders are an important pathology in SA populations, it should be\nmentioned.\nMinor. Page 14. If there are any complete exons that are mapped to the novel contigs, it would increase the impact of the\nmanuscript for them to be included in a main figure showing the context of the exon.", "reviewer_4": "The proposed manuscript by Arun Das et al. uses de novo assembly of unmapped short-read sequences to explore hidden\nvariation in South Asian populations not represented in the current reference genome(s). The authors report results from\napproximately 20,000 placed contigs that they were able to anchor to the reference genome and show how a subset of these\ncontigs could be validated using existing long-read sequencing data. Since non-European ancestries are underrepresented\nin genomic studies, this manuscript fills a gap in the understanding of genetic variation in South Asian populations that could\npotentially benefit future population-wide studies. Apart from some shortcomings in the analysis of these contigs, the short-\nread assembly approach is well-described and sound. The manuscript could be improved by showing how these unplaced\nand placed contigs can be used in future studies using augmented references, for instance, and genotyping of the contigs\nwould be useful to demonstrate the relevance of this hidden variation in GWAS studies.\nDetailed comments:\n* The proposed method does not appear to generate genotypes for all placed contigs based on reference and alternative\nalignments. This would be useful for some of the proposed downstream analyses. For example, the GWAS analysis is\ncurrently based solely on proximity to GWAS hits. Are any of these contigs (or variants within them) in high LD to GWAS\nvariants? Similarly, the current comparison with PopIns2 at the contig level could be extended to include genotypes. For\nbenchmarking, the SAS trios that are part of the extended 1000 Genomes cohort (Byrska-Bishop et al.) should be used to\ne\nto evaluate Mendelian inheritance errors of shared contigs.\n* It would be useful to broadly classify the genome sequence present in these contigs. Are they predominantly tandem\nrepeats like VNTRs, to what extent do the assembled contigs contain mobile element sequences, and what percentage is\nlargely unique sequence? In particular, L1s and SVAs tend to be larger than 3 kbp, and these common mobile element\ninsertions (MEIs) have been relatively well identified by previous alignment-based studies, even when specific insertions are\nmissing from the reference.\n* HPRC Release 1 and HGSVC collected fewer than 10 samples from SAS populations, and although HPRC Release 2 will\nincrease this number to 36, the authors rightly point out that larger studies are needed to capture rare South Asian variation.\nHowever, parallel, medium-coverage long-read projects such as Gustafson et al. and Schloissnig et al., with approximately\nonly validate the majority of the reported contigs but also place a larger proportion of the contigs in the reference genome,\nwhich could ease the use of these contigs in future studies (see next comment).\n* The concluding sentence of the manuscript, that the analysis presented in the manuscript provides \"a model for expanding\nbeyond a single reference genome\" should be better supported by the results. How do the authors envision the use of this\nresource in future South Asian population studies? How should these placed and unplaced contigs be integrated into\npangenome graphs, or how should linear reference genomes be extended with these additional contigs? Does this improve\nvariant detection, and are additional variants on these contigs reliably detected?\n* I found it surprising that only 14% of the assembled contigs are rare (<10% of individuals), as this frequency distribution\ndoes not seem to correspond to the observations in 1000 Genomes, but these were for variants, not contigs. The manuscript\nstill reports greater African diversity on the contig-level, so it would be useful to discuss why there are rather few rare contigs\n- is the assembly approach perhaps underpowered for rare variation?\n* Many of the cited bioRxiv preprints have already been published, sometimes years ago, such as Aganezov et al. (Science,\n2022), Nurk et al. (Science, 2022), Taylor er al. (Nature, 2024) or Gustafson et al. (Genome Res, 2024).\n* Minor comments: Please clarify the term \"non-LOC\" gene. Figure S32 was illegible in my printout. Perhaps the\ngrayscale could be adjusted?"}, "author_response": "To the editor,\nThank you and the reviewers for your helpful comments. Please find attached our revised\nmanuscript along with our point-by-point response to the reviewer queries. In the revised\nmanuscript we have highlighted our major changes in the text (for sections that have moved\nwithin the manuscript but have not had any content changed, we have highlighted the title) and\nin the point-by-point we include the reviewer comments in black and our responses in blue.\nText added to the manuscript specifically to address a particular comment appears in our\nresponse, and is in blue and italicized.\nThe major changes were as follows:\n1. Investigated the repeat content of the assembled contigs using RepeatMasker.\n2. Aligned previously unplaced sequence and sequences with large numbers of RNA-seq\nalignments to the HPRC pangenome, in an effort to provide more context for the\nsequence.\n3. Added details and discussion addressing the distribution of shared and unique\nsequences, and the impact of related trios within the set.\n4. Further investigated the potential technical artifacts behind the results observed in PCA\nclustering.\n5. Added a second long read sequencing dataset to further validate the assembled contigs\nfrom 20 more individuals.\n6. Updated Figure 1A, and restructured the sub-plots in Figure 3, 4 and 5 to improve\nreadability and flow.\n7. Released the assembled contigs and auxiliary on Zenodo, in addition to their prior\navailability on AnVIL via the project GitHub.\nWe also made several edits throughout the text and figures to clarify particular points of\nuncertainty or confusion, to explain certain experimental choices, to add citations to existing\nrelated work, and to improve the readability of the manuscript. These edits are not highlighted,\nbut are visible in the “track changes” version of the manuscript and supplement that are\nincluded in this submission (the two _TRACKED.docx files).\nThe GitHub page (https://github.com/arun96/SouthAsianGenomeDiversity) has also been\nupdated accordingly, and updated versions of all figures and supplementary tables can be\naccessed from there (in addition to being uploaded with this submission). Data is now available\non Zenodo at https://zenodo.org/records/17419004.\nThank you for your consideration,\nDr. Arun Das\nREVIEWER COMMENTS\nReviewer #1 (Remarks to the Author):\nDas and colleagues present a significant reanalysis of unmapped sequence reads from 640\nindividuals from both the Thousand Genomes Project and the Simons Diversity Project to\nassess the ability to identify \"lost\" sequences failing to map to the Human Reference Genome.\nThe authors validate these primarily short read sequences with long-read sequences from\nmultiple individuals, RNA-seq data and BLAST analyses to characterize the extent of novel\nvariation found in South Asian populations. The manuscript is exquisitely written and presented\nand will provide needed data on population variation and further expand upon the need for the\nhuman pangenome to include many more ethnicities. My only real issue with the manuscript is\nthe apparent neglect in comparing their analysis pipeline with other works that have previously\nutilized unmapped reads from species including zebrafish (Faber-Hammond and Brown 2016,\n[Zebrafish]), dog (Holden et al., 2018), human (Chekanov et al., 2010; Genovese et al., 2013;\nLiu et al, 2014; Faber-Hammond and Brown, 2016[Human Genetics]) and other. While the\ncurrent project expanded upon the results of these researchers, using long-read sequences, the\nT2T reference genome, and RNA-seq data to provide more in depth assessment of unaligned\ncontigs, the finding of significant sequence data not being mapped to reference genomes is\nconsistent.\nThe methodology, detail, interpretation, and conclusions made by the authors are all\nexceptional.\nI have no specific revisions as the manuscript was exceptionally well written and complete.\nThank you for your comments! And thank you for highlighting prior work that has utilized\nunmapped reads to assemble and discover variation. We have updated the introduction to\ninclude a reference to past human and non-human efforts before the African Pangenome\nProject, which was one of our major inspirations.\nPage 2, Line 40: “This process involved assembling unmapped reads into larger contigs, which\nhas been demonstrated to reveal novel and variant sequences in both human and non-human\nsamples [21,29–33].”\nReviewer #1 (Remarks on code availability):\nThe github site has a very detailed readme file with an explanation of the pipeline and how to\nutilize this. Any researcher versed in basic bioinformatic techniques should be able to utilize this\npipeline for their own analysis.\nThank you!\nReviewer #2 (Remarks to the Author):\nThe authors explored genetic variation often overlooked by conventional genome mapping\nmethods by analyzing unmapped and poorly mapped reads from 640 South Asian individuals—\n601 from the 1000 Genomes Project and 39 from the Simons Genome Diversity Project. They\nassembled tens of megabases of novel, non-redundant sequences and compared them across\nglobal populations, revealing substantial hidden variation that remains undetected even with\ncomplete genome and pangenome references. Their analysis uncovered thousands of contigs\nintersecting protein-coding genes and GWAS loci, and these findings were further supported\nusing long-read sequencing and RNA-seq data. This work highlights an innovative strategy for\nuncovering overlooked genetic diversity. Notably, the authors took an innovative approach to\nbetter characterize South Asian genetic variation and made all code publicly available, which is\nan important contribution to transparency and reproducibility. However, some observed patterns\nin the data remain unexplained, and additional analyses could further enhance the impact and\ninterpretability of the study.\nThank you for your comments!\nMajor:\n1- The authors used all South Asian individuals in the most recent 1000G high-coverage release\n“601 1KGP SAS individuals”. However, in this release, 85 trios (mother-father-children) from\nSouth Asia are included. The relatedness between these individuals can have an impact on\ndifferent measures of the study: e.g. Linkage Disequilibrium or number of shared unmapped\nreads, as they will share more haplotypes than unrelated individuals. It’s something important to\nconsider and discuss, as a quality check or as a measure of heritability of unmapped reads for\nexample.\nThe majority of our analyses are not affected by relatedness. However, as noted by the\nreviewers, there are a few sections where trio relatedness may influence our results:\n1. The allele frequency spectrum will be slightly distorted compared to expectations for an\nunrelated sample, for example with an enrichment of doubletons (shared between one parent\nand offspring) versus singletons. The allele frequency spectrum is simply presented as a\nqualitative observation, and this modest distortion is now noted in the main text.\nPage 14, Line 20: “The allele frequency spectrum demonstrates that a large fraction of the\ncontigs are rare (present in <10% of the combined SAS set), with the majority of these present\nin just 1 or 2 individuals, consistent with the known abundance of rare variation in human\npopulations. (Figure 4c, Supplementary Figure 12-13). We note that the presence of trios\ndistorts the frequency spectrum compared to expectations for an unrelated sample, but the\npractical effect is minimal in practice (e.g., enriching for doubletons versus singletons).”\n2. In relation to the point above, variants that occur in multiple samples due to relatedness will\nbe double-counted as overlapping with genes or GWAS hits. This is addressed by our analyses\nthat collapse results to unique genes and GWAS hits across the entire sample.\nPage 11, Line 6: “These intersections affected a total of 1,189 unique transcripts, 172 genes\n(106 protein-coding), 33 exons and 6 CDS regions (Figure 3c).”\nPage 18, Line 10: “Loosening our search criteria, we identified >6,000 contigs within 100 bp of\n36 unique GWAS hits and >15,000 contigs within 1 kbp of 112 unique GWAS hits.”\n3. In LD based placement, the presence of related individuals in the set will inflate LD genome-\nwide, but crucially is not expected to influence local peaks of extreme LD that inform our\nplacements. This is because we detect association between contigs based on their presence\nand absence across all 640 individuals in the set (computed as the Pearson correlation\ncoefficient (r2) between two 640 x 1 arrays, where each index in the array is presence/absence\nin each individual). The addition of a related parent/child will only affect the presence or\nabsence of a contig in at most one other individual, and does not affect the remaining 600+\nindividuals across which the contigs may or may not be associated. Therefore, we do not expect\nthis to affect the majority of contigs. The only contigs that may be substantially affected are\n“doubletons” shared between one parent and a child; these doubletons may see slightly\nelevated levels of association through the presence of both individuals in the set.\nFurthermore, many of the hard to place contigs are rare variants, occurring in just one person\n(and impossible to place using LD). The relatedness of a parent/child will only add a single\noccurrence of this contig, while not affecting its occurrence in, and thus association across, the\nremainder of the set.\n2- In the PCA of unaligned reads, individuals cluster by dataset (SGDP vs. 1000G) rather than\nby population, which is unexpected given known patterns of genetic diversity and suggests the\npresence of technical batch effects. This raises concerns about the direct comparability of the\ntwo datasets and may influence the interpretation of population-specific signals. The authors\nmentioned the potential impact of sequencing depth on differences between datasets “We\nattribute the larger span of non-contaminate contigs due to the higher sequencing depth of\nSGDP compared to 1KGP” but no specific values are provided in the main text or Methods.\nIncluding these details would help assess the extent to which sequencing depth differences\nmight explain the observed clustering, and help discussing their potential impact on the study’s\nconclusions.\nWe have attempted to address these issues in the following ways:\n1. Added coverage details: We have amended the manuscript to include the exact coverage\ndepths of the 1KGP and SGDP sets in the “Alignment and assembly of data from the SGDP\nSAS set” section, to better highlight that difference.\nPage 6, Line 42: “We attribute the larger span of non-contaminate contigs to the higher\nsequencing depth of SGDP (32-66x, average coverage of 43x) compared to 1KGP (~30x).”\n2. Sequencing artifacts: As discussed in the manuscript, there are two individuals from SGDP\nwith obvious sequencing artifacts, further highlighting that there are known and recognized\nbatch effects between the two sets. These individuals are excluded from our results in an effort\nto minimize their impact on the study.\n3. Added constraints during clustering: In Supplementary Figure 29, we included a range of\nPCA plots with a number of different filtering criteria, including one where we include contigs if\nthey are present in both the 1KGP and SGDP set (in an effort to reduce batch effects).\nHowever, we still observe clustering based on the dataset and not by population.\n4. Accounting for coverage: We performed downsampling on the SGDP reads to reduce their\ncoverage to the level of the 1KGP set; however, we still observed clustering based on dataset\nand not by population. Similarly, constraining to only high coverage contigs did not remove the\neffect.\nIn short, despite our various tests and filtering operations, we did not attribute the PCA clusters\nto a single biological or technical variable, and instead conclude that a more subtle combination\nof multiple factors (data characteristics, artifacts, underlying population structure, relatedness,\netc.) must account for these patterns. This is evidenced by the relatively small proportion of the\nvariance explained by the first PC (<15%), dropping immediately to <3% for PC2, which are the\nmain PCs on which the batches are separating.\nHowever, we do not believe that these artifacts distort our primary conclusions for a few\nreasons. We remain confident that the sequences being assembled are true, non-contaminant\nhuman sequences (as evidenced by our filtering and auxiliary analyses), and that the methods\nused to group and place them remain robust. Our shared sequence analysis evaluates the\ndistribution of contigs across one or more individuals regardless of population, and none of our\nanalysis or results are dependent on a clear division between the various subpopulations.\nFinally, due to the complicated history and underlying substructure of sub-populations in South\nAsia, a clear split was unlikely to occur, and the high degree of similarity between these\npopulations means that even small technical artifacts will be captured on the first principal\ncomponents.\n3-a. Still related to the PCA of unaligned reads, I noticed that the 1000 Genomes individuals\nappear to separate into two distinct clusters, which do not correspond to population structure.\nThis raises the possibility that other factors—such as the sex of the individuals—may be driving\nthis separation, especially given the reported enrichment of assembled contigs on chromosome\nY. Investigating this further could provide valuable insight and strengthen the study’s findings.\n-b. Additionally, the journal encourages to follow to the SAGER guidelines\n(https://researchintegrityjournal.biomedcentral.com/articles/10.1186/s41073-016-0007-6), which\nrecommend conducting analyses that can reveal sex-related differences. In this context, it would\nalso be helpful to report the number of male and female individuals included in the study.\nThank you for your comments. We have amended the manuscript to include a PCA plot colored\nby karyotype (XX/XY), now found in Supplementary Figure 30. As seen in that figure, the\nclustering is not driven by the sex of the individuals. This is also confirmed by inspecting the\nPCA loadings, which trace to markers across the genome and are not concentrated on the sex\nchromosomes. Thus, the same conclusions regarding separation of SGDP and 1KGP samples\nthat we detailed in the previous response (i.e., subtle batch effects) also apply within the 1KGP\ndataset. Nevertheless, as also noted in the previous response, we are confident that the major\nconclusions of our study do not hinge on the absence of such technical variation and are\nvalidated through orthogonal means (e.g., based on independent long-read sequence data).\nAs suggested by the reviewer, we now report the number of male and female individuals in the\nstudy in the “Data” section in the Methods.\nPage 23 Line 10 and 14: “In the 1KGP set, there are 282 XX individuals and 319 XY\nindividuals… In the SGDP set, there are 8 XX individuals and 31 XY individuals.”\n4- The study focuses on South Asian populations, yet it does not engage with key aspects of\nSouth Asian evolutionary history—such as the region’s high genetic diversity between sub-\npopulations, high patterns of homozygosity, or population structure—which could provide\nimportant context for interpreting the findings. For example, do individuals who share the same\nunmapped reads are also more closely related (e.g. share more Identical-by-Descent\nsegments)? Exploring this could help clarify whether the shared sequences reflect deep\npopulation structure, recent relatedness, etc.\nThis is a good suggestion. We have added a section to the introduction of the manuscript\nhighlighting the unique challenges and research opportunities regarding South Asian\nevolutionary history.\nPage 2, Line 9: “Population genomic research focused on South Asia has uncovered deep\npopulation structure, runs of homozygosity [12] attributed to founder effects and endogamy, as\nwell as intricate histories of both ancient and recent admixture, all of which are influenced by\nsocietal and cultural factors [12, 13–17]. Although these studies are invaluable, South Asians\ncompose only ~2% of individuals in genome-wide association studies (GWAS), and even within\nthis broader group, certain ethnicities are over- or underrepresented [1] (Supplementary Figure\n1).”\nTo investigate sharedness in the set, we plotted the number of individuals each shared contig\nwas present in (Supplementary Figure 13). We see the expected number of singleton and\ndoubleton variants, consistent with what we know about variation in human populations and the\nnumber of trios in the set. We also see a large number of contigs present in larger subsets of\nindividuals, highlighting the high amount of sequence shared between non-trios across the\nentire set. This is in line with past work discussing the distribution of variation in South Asia.\nWe also compared the similarity in shared contigs between individuals to the lengths of their\nshared IBD segments (from the IGSR Data Portal https://www.internationalgenome.org/data-\nportal/). We did not find a clear association between these two metrics, with a Pearson\ncorrelation coefficient ~0.05 (plot included below). This comparison only takes into account long\n(>1kbp) segments of sequence that have already undergone collapsing if they were >90%\nsimilar; in reality many, much shorter segments may be shared, and contigs that were >90%\nsimilar but not identical would be considered to be different sequences.\nIt would be interesting to see this type of analysis attempted when working with a biobank-scale\ndataset (such as that used by Williams et al, 2025), using short sequences or even reads, or\nwith a different measure of similarity between each pair of individuals.\nAdditionally, while GenomeAsia 100K Project is briefly mentioned, other major regional efforts\nsuch as GenomeAsia Phase 2 dataset (Wall et al. 2023) or the ongoing Genome India project\n(Bhattacharyya et al. 2025) are not cited. Including and situating the work within the broader\nlandscape of South Asian genomics would enrich the relevance and impact of the study.\nThis is an important point. In our work, we focused on open-access data resources, but our\nanalytical approach can be generalized to any relevant dataset (including the two mentioned in\nthe comment). This is a research goal that we are eager to enable and pursue.\nBhattacharyya et al is already cited in the discussion, but we have added a citation to Wall et al.\n2023 into the updated discussion.\nPage 21, Line 29: “The generation of more high quality DNA sequencing and genome\nassemblies, both within South Asia [12,56,57] and globally, will further expand the amount of\nvariant sequence, while auxiliary data modalities will help validate these sequences and\ncharacterize their functional potential [56]”\n5- The authors state: “We found that, on average, 3.8–4.5% of reads from each individual align\nto this set of contigs, with individuals of all 21 non-SAS populations exhibiting roughly the same\nalignment rate.” This is an intriguing result, and it would be valuable to contextualize it—are\nthese alignment rates consistent with other “non-redundant” set of population-specific contigs?\nAlso, is this a genome-wide pattern (or is it linked to one chromosome/one genomic region)?\nWe do not generate a “non-redundant” set of population specific contigs for the non-SAS\ngroups, as assembling contigs specific to other continental groups was outside the scope of our\nstudy. All workflows used to generate the contigs and the non-redundant set are available for\npublic use, so future studies can repeat this with other, non-SAS populations.\nHowever, given the distribution of genetic variation across human populations (either globally\ncommon or rare even within populations), we anticipate this result would hold across any\ncontinental group. The one exception would be the African populations, where there is a subset\nof common variation within Africa that was not sampled during the out-of-Africa bottleneck. This\nis partially seen in the African Pangenome effort (Sherman et al, 2019) where alignment rates\nagainst non-reference global sequences improve when using variation found in African\nindividuals, as common variation in non-AFR populations is also present within Africa.\nFor the genome-wide distribution of these alignments, we did not observe a concentration in any\nparticular genomic region. However, as only a small fraction of the set is placed, we cannot yet\ndetermine how the reads aligning to unplaced contigs would be distributed in the genome—an\nintriguing question for future work\n6- The authors report that “individuals in the SAS set had on average 10–15% less assembled\nsequence than the individuals in the seven African (AFR) populations, and on average 5% more\nsequence than the individuals in the non-AFR, non-SAS cohorts.” While these differences are\npotentially important, it would be helpful to know whether they are statistically significant. More\nbroadly, the manuscript would benefit from the inclusion of formal statistical tests throughout to\nsupport key claims and enhance the rigor of the findings.\nThank you for your comment. As suggested, we now use t-tests to compare the amount of\nsequence assembled per individual between the different pairs of continental groups:\nPage 9, line 39: “However, individuals from the SAS set produced on average 10-15% less\nassembled sequence than individuals from the seven African (AFR) populations (two-tailed\nWelch's t-test; t = -36.8; two-tailed p-value = 1.103 x 10-46), and on average 5% more\nsequence than the individuals from the non-AFR, non-SAS populations ((two-tailed Welch's t-\ntest; t = 6.22; two-tailed p-value = 5.39 x 10-9) (Supplementary Figure 7, Supplementary Table\n9). Performing population level comparisons, where the sequence assembled across 601 SAS\nindividuals are compared just to the 10 individuals from each non-SAS population, we see all\nseven AFR populations have a two-tailed p-value < 1.76 x 10-7, and the non-SAS populations\nCHS, CDX, PEL, GBR, FIN and CEU have a two-tailed p-value < 0.05 (Supplementary Table\n9).”\nSecond, we compared the difference in post-filtering, >1 Kb sequence assembled from CHM13\nand GRCh38 unmapped reads.\nPage 6, Line 8: “The analogous approach using GRCh38 assembled an average of 2 Mbp of\nsequence per individual, significantly exceeding the above results from T2T-CHM13 (two-tailed\nWelch’s t-test, t = -134.45, p-value < 1 x 10-308), and underscoring the greater completeness of\nthe T2T reference.”\nAnd third, we compared the difference in post-filtering, >1 Kb sequence assembled from\nCHM13 and Minigraph unmapped reads.\nPage 9, Line 20: “Using reads that failed to align to the Minigraph pangenome reference, we\nassembled 450-500 kbp of sequence per 1KGP SAS individual in >1 kbp non-contaminant\ncontigs, a ~10% decrease in sequence compared to the same individuals using T2T-CHM13\n(Welch’s t-test, two-tailed p-value of 2.03 x 10-13; Figure 2e-f).”\nWe agree that the application of formal statistical tests has improved the rigor and\ninterpretability of our quantitative statements.\nMinor:\n1-The figures are not cited in order within the text—for example, Figure 2b is referenced after\nFigure 2f. Additionally, some figures such as Figure 1 (a and b), Figure 2a, and Figure 3d do not\nappear to be cited at all. Notably, Figure 2a compares alignment rates for a method that is only\ndescribed in the Methods section; this figure should either be removed or properly cited in the\nmain text alongside the corresponding method. Improving the order and completeness of figure\ncitations will enhance the clarity and flow of the manuscript.\nThank you. We have corrected this by adding references to Figures 1A and 1B into the\nintroduction, added a reference to 2A in the first paragraph of the results section, and a\nreference to 3D in the “Intersection with annotated elements” section on page 11.\nPage 11, Line 9: “Qualitatively similar results were observed when applying the same approach\nto samples from non-SAS populations (Figure 3d)”\nWe have also made a number of changes to Figures 3-5 to ensure the order is kept when\nreferring to their subfigures, and changed the order of some supplementary figures to try and\nmaintain the order in which they are referenced.\n2-Figure 1A would be clearer and more informative if it included the sample sizes for each\ngroup.\nWe have edited the figure to enlarge the circles and add the sample sizes for each of the 1KGP\ngroups, added a bar chart breaking down the exact counts from each group, and added details\nto the caption for Figure 1A.\n3-Figures S6, S8, S9, S13, S18, S19, and S26 are difficult to read due to low resolution and/or\nsmall font size. Improving the clarity and readability of these figures would greatly enhance the\nquality and accessibility of the supplementary material.\nWe apologize for the difficulty in reading these figures. These figures become readable when\nfull-sized, which was not an option during the submission of the current version of the\nmanuscript. This is especially true for Figures S8 (now S11), S13 (now S10), S18, S19, and\nS26, where we are displaying hundreds of datapoints,as well as for figures where there are\nmultiple panels (S6 and S9 (now S12)). For figures with a corresponding table (S13 (now S10),\n18, 19, 26 and 27), we have added a reference to the table in the caption.\nFull size, high resolution versions of all figures will be submitted with the final version of the\nmanuscript, and they can also all be accessed through the project GitHub page.\n4-Page 6 line 16, the acronym “APG” is used without being defined.\nThank you for catching this. We have expanded this to “African Pangenome”, and added the\nacronym “APG” alongside the first appearance of the project in the introduction.\nPage 2, Line 28: “This reference bias was highlighted in the African Pangenome project\n(APG)..”\nPage 6, Line 31: “...reflecting its completeness and mirroring results from the African\nPangenome study…”\nReviewer #2 (Remarks on code availability):\nThe GitHub repository provided by the authors is exceptionally thorough and well-maintained. It\nstands out as one of the most complete and user-friendly repositories I’ve encountered. The\nclarity of the documentation, the reproducibility of the workflows, and the accessibility of the\ncode all reflect a strong commitment to transparency and open science. This resource will likely\nbe valuable not only for replicating the study but also for other researchers working in related\nareas.\nThank you!\nReviewer #3 (Remarks to the Author):\nThe authors present analysis of assembled contigs from South Asian individuals. The extra\nsequences were found to overlap coding sequences and GWAS hits. The methodology is\nstraight forward and fairly complete, and the multiple well executed vignettes give functional\ncontext for what missing sequences were assembled. The LD placement of contigs is a\nclever trick (similar to work by Giulio Genovese https://pmc.ncbi.nlm.nih.gov/articles/PMC3971454/ -\nmay be worth noting here). The low overlap with Manta insertions is expected.\nSome gaps in analysis remain:\nThank you for your comments. We have reviewed the work by Genovese et al., which we agree\nis relevant, and have added a brief citation to this paper in the methods section.\nPage 27, Line 23: “This approach “links” pairs of contigs with high degrees of association; it may\nbe possible to extend this to tie together and place larger groups of contigs or incorporate other\ndata into the association [33]”.\nSignificant.\nUnless I’m missing something, the assemblies are not released with the manuscript. Because\nthis is analysis of public, open consent data, the data being released for review is a condition for\nacceptance. These should be uploaded to zenodo, and it will require an additional round of\nreviews as the reviewers examine the contents of the assembled contigs.\nThe final assembled contigs, as well as several intermediate steps, were all made available on\nAnVIL when the manuscript was submitted. We have now also uploaded the final contig\nassemblies to Zenodo (doi.org/10.5281/zenodo.17419004, see Data Availability statement for\nfull details).\nThe number of contigs that are not annotated as having homology (placement) to the reference\nis high. Similar experience comparing inserted sequences back to the genome in HPRC\nanalysis showed that the majority of inserted sequences are repetitive elements. Running blast\nagainst NR is a great first step, but more analysis should be done (partially relating to the next\ncomment).\nIn the original version of the manuscript, we had already attempted to align all contigs\n(regardless of placement against CHM13) against the HPRC Minigraph pangenome, and saw\nvery similar results to those observed against CHM13. We have now repeated this process,\ncombining the outputs of several pangenome aligners and using more tailored parameter\nsettings.\nCompared to CHM13, we are now able to place ~1,000 more contigs with alignment only\nagainst the HPRC pangenomes. ~900 of these were originally unplaced, and this accounts for\n0.5% of the previously unplaced contigs. These placements were all across the HPRC\npangenome build, with no significant enrichment in any region of the genome. Due to the lack of\ncentralized annotations and equivalent pangenome tools, we cannot replicate the linear\nreference analysis for these placed contigs, but we hope to be able to with future pangenome\nreleases.\nWe have added these counts to the “Analyzing pangenome placements” section of the\nmanuscript, as a follow on to the BLAST analysis reported in the previous “Analysis of unplaced\nsequence” section.\nPage 13, Line 31: “The limited availability of pangenome tools precluded extension of the\nplacement approach to the graph pangenomes. As an alternative, we performed conventional\nalignment to the linear references that comprise the pangenome, identifying 8.5% of the\nassembled contigs (16,854/199,564) with an alignment length >500 bp, thus slightly exceeding\nthe number of contig alignments to T2T-CHM13 alone (15,739).\nOf the contigs that failed to be placed against T2T-CHM13, 905 of 179,873 (0.5%) aligned to the\nHPRC Minigraph pangenome. These alignments occur throughout the genome, with no specific\nchromosome or region especially enriched for these placements. These results highlight the\npower of pangenome approaches while emphasizing the need for centralized annotations and\nauxiliary data to facilitate downstream interpretation.”\nWe also analyze the repeat content of these sequences as you suggested (see below).\nThe repeat content of the new contigs should be annotated. One would expect the majority of\nthe new sequences to be derived from mobile elements and tandem repeat expansions.\nWe have run RepeatMasker on the final assembled contigs, and added the key results to the\nmanuscript.These details, as well as details on the syntax used to run RepeatMasker, have\nbeen added to the main text, methods and as a supplementary note. The masked sequence\nand outputs from RepeatMasker are available via the GitHub page for the project.\nMain Text: Page 6, Line 13: “Application of RepeatMasker [39] showed that 4.39% of the\nassembled sequence consisted of short interspersed nuclear elements (SINEs) (1.73% ALUs,\n2.66% MIRs), 17.89% long interspersed nuclear elements (LINEs) (13.85% LINE1, 3.89%\nLINE2, 0.11% L3/CR1), 8.54% long tandem repeat elements (LTRs) (3.40 ERV-L, 3.52% ERVL-\nMaLRs, 1.04% ERV Class I, 0.04 ERV Class II) and 2.53% DNA elements, together constituting\na total of 33.40% interspersed repeats (Supplementary Table 8). An additional 2.2% of\nsequence consisted of small RNA (0.03%), satellites (0.64%), simple repeats (1.21%), and low\ncomplexity DNA (0.32%). In comparison, the contigs assembled from GRCh38 unmapped reads\ncontained substantially more SINEs (8.5%), satellites (35.5%) and simple repeats (17%), but a\nlower proportional composition (but greater absolute composition) of LINEs (7.5%), LTR\nelements (3.56%) and DNA elements (0.9%).”\nMethods: Page 24, Line 27: “Repeat content in the assembled contigs and reference genomes\nis detected, masked and annotated using RepeatMasker v4.1.2-p1 33. Using the assembled\ncontigs as input, RepeatMasker detects and masks SINEs, LINEs, LTR and DNA elements,\nsmall RNA and satellites, and depending on the use of the “nolow” parameter, simple repeats\nand low complexity elements (Supplementary Note 11).”\nMinor.\nThe GWAS hits should be compared to the number of hits that one would have randomly\nshuffling the aligned contigs, taking into account repeat content of the contigs (e.g., similar\nrepeat content of the locations of the shuffled alignments as original alignments).\nOur original intention was not to make a claim about enrichment, but simply to use GWAS data\nand other annotations to interpret the functional and phenotypic potential of newly recovered\nsequences.\nNevertheless, comparing the overlaps of placed contigs with genes against the number of\noverlaps seen when randomly placing contigs, we see that the random placements result in\nmore gene intersections (~14.5k across 3 runs) than the observed counts (12,021). This is as\nexpected, given that genes are generally highly conserved regions under purifying selection. We\nalso observe similar results when counting the number of GWAS sites within 100 or 1,000\nbases of a placed contig. We have now integrated this result into the manuscript as follows:\nPage 11, Line 16: “The observed count of 12,021 gene intersections was substantially lower\nthan the number of intersections obtained through random placement of contigs against the\nreference genome (14,500 intersections, averaged over three runs), consistent with the action\nof purifying selection against common structural variation in genic versus non-genic regions.”\nThis is more of a commentary than a comment, but the data that were analyzed were all public.\nSo why focus only on South Asian individuals and just 10 individuals from other populations?\nThe assembly pipeline is clearly robust. The context of 1.8B individuals is important, but having\nadditional populations (more than 10 samples) to gauge how much sequence is missing from\nthese individuals versus from other populations, for context.\nWhile there is certainly merit in expanding this approach to individuals from all over the world,\nthere are several reasons why we have focused on South Asia specifically.\n1. South Asian individuals are among the least well represented individuals in modern genomics\nrelative to their census population size.\n2. While the majority of our work is population agnostic, expanding GWAS analysis to globally\ndiverse populations can be complicated, as populations differ in both allele frequencies and\nstructure of linkage disequilibrium, not to mention the potential for gene-by-environment (GxE)\nand gene-by-gene (GxG, i.e., epistatic) interactions.\n3. From a practical perspective, there are substantial costs associated with cloud-based\nanalyses at this scale ($5-10 per sample). We have focused our efforts on 640 SAS samples to\nkeep these costs manageable, but we are interested in scaling up the approach conditional on\nfuture funding.\n4. Building on the last point, our pipeline and methodology is now available and usable by the\nresearch community, enabling application to other cohorts.\nIn summary, our goal was to improve genomic resources for a large understudied population,\nwhile creating resources and analysis pipelines that can be used for other populations in the\nfuture.\nLines 3-6, 10-14, page 15 offer qualitative descriptions of results (“most of…”). This should be\nchanged to offer specific percentages.\nWe have updated the comparisons in that section to include exact percentages instead of\nqualitative descriptions.\nPage 16, Line 5: “We used RNA-seq data from the MAGE project 49, obtained from\nlymphoblastoid cell lines from 731 globally diverse individuals from 1KGP, to evaluate the\nfunctional potential of the previously assembled contigs. We specifically focused on the 140\nMAGE samples from SAS populations that overlap with our study. We aligned RNA-seq reads\nfrom each individual to their contigs that were assembled from reads that failed to align to T2T-\nCHM13. Per individual, ~800-1,000 reads aligned to their respective contigs, with the majority\n(>80%) spread across few contigs. While this represents a small fraction of the input RNA-seq\ndataset (<0.1%), this is not unexpected, as the total length of assembled contigs per individual\nis <0.02% the size of the human genome. Across the 140 individuals, 90% of contigs with >100\nRNA-seq alignments had BLAST hits to non-reference human and primate sequences, some of\nwhich is annotated as functional (Supplementary Table 3-4). These include hits to annotated\nnon-reference sequences from chromosomes 1, 2, 3, 7, 9, 15, X and Y.\nWe also aligned the RNA-Seq data from all 140 individuals against a combined set of filtered\ncontigs from all 601 1KGP SAS individuals and identified 200 contigs with >100 read alignments\nand 19 contigs with >1,000 read alignments (Figure 4e-f, Supplementary Figure 17-19), 83% of\nwhich had BLAST hits to non-reference human and primate sequences (Supplementary Tables\n3-4). Across both individual-specific and combined RNA-seq alignments, we found that the\nmajority (>80%) of highly aligned-to contigs appear in just one or a few individuals, with few\nclassified as widely shared sequences (Figure 4d).”\nThe assembled sequences with RNA-seq alignments should be searched in the HPRC\ngenomes.\nUsing the alignments to HPRC mentioned in a previous comment, we see that only a few (4) of\nhigh-scoring (>100 aligned reads) contigs that were unplaced against CHM13 are now placed\nagainst the HPRC Minigraph build. We have added these details to the ‘RNA-seq Analysis’\nsection of the manuscript.\nPage 16, Line 26: “Intersecting the placement information for these contigs with the RNA-seq\nalignments, we find that the majority of highly aligned-to contigs could be placed via alignment\nto T2T-CHM13 (140/200 of contigs with >100 alignments, 12/19 contigs with >10 alignments),\nwith most of the placements tracing to the Y chromosome. Aligning these contigs against the\nHPRC Minigraph pangenome reference modestly increased the rate of alignment (144/200 of\ncontigs with >100 alignments, 13/19 contigs with >10 alignments), with few additional contig\nalignments to chromosomes 19 and Y.”\nThe shared sequence analysis should state the 5% divergence allowed, based on the\ncommands in the supplement. This is important for the context of how divergent sequences can\nbe and considered shared.\nThank you for pointing this out. We have updated the “Shared sequence analysis” section to\ninclude a mention of this, and specifically included the fact we used the “asm5” preset for\nMinimap2 in the corresponding Methods section.\nMain text: Page 14, Line 35: “To address such cases, we performed an all-vs-all alignment of\nthe assembled contigs using Minimap2 [45], and considered contigs with >90% alignment to\nhave originated from a shared allele.”\nMethods: Page 25, Line 14: “The all-vs-all alignment was performed using Minimap2 v2.22 [41]\nwith the strictest “Long assembly to reference mapping” default parameters (asm5), and the\nresulting output file was used to collapse the shared sequence.”\nPage 11, gene intersections. If these are intronic, they are probably of low functional impact. If\nthere is some additional context to believe that eye-related disorders are an important pathology\nin SA populations, it should be mentioned.\nIn this work, we do not make any claims about particular disorders being enriched or more likely\nin SA populations, and instead just report the number of gene intersections and the most\ncommonly intersected genes. We also report the number of exons intersected on the same\npage in the previous section, and full details of all gene and exon intersections are provided in\nthe supplementary tables.\nThe 1KGP dataset does not contain any phenotype data, so future analyses in other datasets\nwould be needed to confirm if these intersections resulted in any clinical outcomes.\nPage 14. If there are any complete exons that are mapped to the novel contigs, it would\nincrease the impact of the manuscript for them to be included in a main figure showing the\ncontext of the exon.\nWe do not explicitly look for any complete exons mapped to the novel contigs. However,\nthrough our alignment of RNA-seq reads back to the assembled novel contigs, we do observe a\nnumber of alignments where reads appear to align to a section of the contig in a manner similar\nto alignments to an exon. Two of the clearest examples for this are in GM21112 (contig\nk141_235) and in HG03717 (contig k141_768), where there is large pileup of RNA-seq\nalignments to a specific location in the contig (both of which are unplaced against the CHM13\nreference). However, it is hard to be definitive that these or other examples are indeed complete\nexons, although the alignment data support our central claim that some of these sequences are\ntranscribed.\nWe have added the IGV screenshots of these two contigs and their RNA-seq alignments to the\nsupplement, in Supplementary Figure 33.\nReviewer #4 (Remarks to the Author):\nThe proposed manuscript by Arun Das et al. uses de novo assembly of unmapped short-read\nsequences to explore hidden variation in South Asian populations not represented in the current\nreference genome(s). The authors report results from approximately 20,000 placed contigs that\nthey were able to anchor to the reference genome and show how a subset of these contigs\ncould be validated using existing long-read sequencing data. Since non-European ancestries\nare underrepresented in genomic studies, this manuscript fills a gap in the understanding of\ngenetic variation in South Asian populations that could potentially benefit future population-wide\nstudies. Apart from some shortcomings in the analysis of these contigs, the short-read assembly\napproach is well-described and sound. The manuscript could be improved by showing how\nthese unplaced and placed contigs can be used in future studies using augmented references,\nfor instance, and genotyping of the contigs would be useful to demonstrate the relevance of this\nhidden variation in GWAS studies.\nThank you for your comments!\nDetailed comments\nThe proposed method does not appear to generate genotypes for all placed contigs based on\nreference and alternative alignments. This would be useful for some of the proposed\ndownstream analyses. For example, the GWAS analysis is currently based solely on proximity\nto GWAS hits. Are any of these contigs (or variants within them) in high LD to GWAS variants?\nSimilarly, the current comparison with PopIns2 at the contig level could be extended to include\ngenotypes. For benchmarking, the SAS trios that are part of the extended 1000 Genomes\ncohort (Byrska-Bishop et al.) should be used to evaluate Mendelian inheritance errors of shared\ncontigs.\nThank you for your comments.\nWhile we do not generate an explicit genotype matrix, we do generate a presence/absence\nmatrix for each of the ~14K collapsed contigs against each of the 640 SAS individuals. We use\nthis matrix to perform LD computations with known variants and other sequences, and functions\nas a “collapsed” genotype matrix.\nGiven the data we have, there are several limitations that prevent us from generating full\ngenotypes. First, it is non-trivial to determine the zygosity of these variants, as read coverage\nand other indicators of zygosity are variable across the contigs. Second, most variants are not\npresent in homozygous alt form (unless the reference genome itself carries the rare variant).\nThird, the short read lengths severely restrict our ability to phase the contigs. Therefore, we\nmostly find and report heterozygous, unphased variants.\nIt would be useful to broadly classify the genome sequence present in these contigs. Are they\npredominantly tandem repeats like VNTRs, to what extent do the assembled contigs contain\nmobile element sequences, and what percentage is largely unique sequence? In particular, L1s\nand SVAs tend to be larger than 3 kbp, and these common mobile element insertions (MEIs)\nhave been relatively well identified by previous alignment-based studies, even when specific\ninsertions are missing from the reference.\nWe have run RepeatMasker on the final assembled contigs, and added the key results to the\nmanuscript (also see R3 above).\nThese details, as well as details on the syntax used to run RepeatMasker, have been added to\nthe main text, methods and as a supplementary note. The masked sequence and outputs from\nRepeatMasker will be made available via the GitHub page for the project.\nMain Text: Page 6, Line 13: “Application of RepeatMasker [39] showed that 4.39% of the\nassembled sequence consisted of short interspersed nuclear elements (SINEs) (1.73% ALUs,\n2.66% MIRs), 17.89% long interspersed nuclear elements (LINEs) (13.85% LINE1, 3.89%\nLINE2, 0.11% L3/CR1), 8.54% long tandem repeat elements (LTRs) (3.40 ERV-L, 3.52% ERVL-\nMaLRs, 1.04% ERV Class I, 0.04 ERV Class II) and 2.53% DNA elements, together constituting\na total of 33.40% interspersed repeats (Supplementary Table 8). An additional 2.2% of\nsequence consisted of small RNA (0.03%), satellites (0.64%), simple repeats (1.21%), and low\ncomplexity DNA (0.32%). In comparison, the contigs assembled from GRCh38 unmapped reads\ncontained substantially more SINEs (8.5%), satellites (35.5%) and simple repeats (17%), but a\nlower proportional composition (but greater absolute composition) of LINEs (7.5%), LTR\nelements (3.56%) and DNA elements (0.9%).”\nMethods: Page 24, Line 27: “Repeat content in the assembled contigs and reference genomes\nis detected, masked and annotated using RepeatMasker v4.1.2-p1 33. Using the assembled\ncontigs as input, RepeatMasker detects and masks SINEs, LINEs, LTR and DNA elements,\nsmall RNA and satellites, and depending on the use of the “nolow” parameter, simple repeats\nand low complexity elements (Supplementary Note 11).”\nHPRC Release 1 and HGSVC collected fewer than 10 samples from SAS populations, and\nalthough HPRC Release 2 will increase this number to 36, the authors rightly point out that\nlarger studies are needed to capture rare South Asian variation. However, parallel, medium-\ncoverage long-read projects such as Gustafson et al. and Schloissnig et al., with approximately\nfrom both studies could not only validate the majority of the reported contigs but also place a\nlarger proportion of the contigs in the reference genome, which could ease the use of these\ncontigs in future studies (see next comment).\nThank you for pointing out those two studies. Note, we were already using the data from\nGustafson et al. to perform long read validation of our assembly contigs.\nOur method is currently focused on using unmapped short reads to recover previously missing\nvariation, so we have opted to not use the long reads to aid placement. Rather, they remain a\ntool to validate these sequences. Future studies can build on our current approach to\nincorporate unaligned (or, as is more likely to be the case, poorly aligned) long reads into\nassembled novel and variant contigs.\nTo that end, we have repeated the analysis performed with the Gustafson et al. samples with 20\nsamples from the Schloissnig et al. set, and we found similar results. Specifically, ~40-50% of\nthe previously unaligned reads align to the basic Flye assembly, and >80% of all assembled\ncontigs have high identity alignments to their personalized assembly. We have added these\ndetails to the manuscript (in the “Validation with long reads section”) and added details of the\nSchloissnig et al. data to the relevant methods section.\nPage 14, Line 22: “We repeat this long read validation using a second dataset from Schloissnig\net al. [40], and perform the same analysis on 20 SAS individuals from this set. The Flye\nassemblies for these individuals are comparable to those from the 1KGP set (2.65 - 2.82 Gb of\nsequence, N50s of 25.6 - 45.7 Mb, and mean coverage of 14-33x), and initial alignment sees\n40-50% of previously unaligned reads now align to the personalized long read assembly. Contig\nalignment is slightly lower with this data (80-87%), but in line with the fact this dataset is\nsubstantially lower coverage than the 1KGP ONT set.”\nThe concluding sentence of the manuscript, that the analysis presented in the manuscript\nprovides \"a model for expanding beyond a single reference genome\" should be better supported\nby the results. How do the authors envision the use of this resource in future South Asian\npopulation studies? How should these placed and unplaced contigs be integrated into\npangenome graphs, or how should linear reference genomes be extended with these additional\ncontigs? Does this improve variant detection, and are additional variants on these contigs\nreliably detected?\nWe have added to the discussion to address these points, and explain how these sequences\nmay be used and the resulting benefits of including them.\nPage 21, Line 37: “The variant sequences discovered using our approach can be used in\nmultiple downstream applications. For example, the sequences could be stored as alternate\nsequence panels or incorporated into a tiered reference querying approach [62], recovering\nreads that fail to align against a single reference. Alternatively, the placed contigs could be used\nto augment pangenome builds, with the placements serving as the coordinates in the\npangenome from which these sequences branch. The addition of novel sequences to the\npangenome has been demonstrated to improve both variant detection and reduce variant\ndiscovery errors, motivating this application [3].”\nI found it surprising that only 14% of the assembled contigs are rare (<10% of individuals), as\nthis frequency distribution does not seem to correspond to the observations in 1000 Genomes,\nbut these were for variants, not contigs. The manuscript still reports greater African diversity on\nthe contig-level, so it would be useful to discuss why there are rather few rare contigs - is the\nassembly approach perhaps underpowered for rare variation?\nRevisiting the section being discussed, we believe that we incorrectly reported the 14% fraction;\nthis was a count from another experiment that was erroneously left in this section. As seen in\nFigure 4e (now 4c), substantially more than 14% of the assembled contigs are rare. We have\nupdated this section, and we apologize for the confusion.\nPage 14, Line 39: “The allele frequency spectrum demonstrates that a large fraction of the\ncontigs are rare (present in <10% of the combined SAS set), with the majority of these present\nin just 1 or 2 individuals, consistent with the known abundance of rare variation in human\npopulations. (Figure 4c, Supplementary Figure 12-13). We note that the presence of trios\ndistorts the frequency spectrum compared to expectations for an unrelated sample, but the\npractical effect is minimal in practice (e.g., enriching for doubletons versus singletons).”\nHowever, you are correct in saying that the frequency distribution still does not fully correspond\nto what we expect from a population genetics perspective. We believe that a combination of\ndifferent factors are behind this observation.\n1. Some of the results reported in our paper are skewed by our decision to “bin” the data, both\nfor plotting and when reporting our results. Varying the bin sizes and groupings can distort the\ndistribution; if the bin size is too small (e.g., reporting the exact number of contigs present in 3%\nof the population), the overall pattern becomes hard to see, whereas larger bin sizes (e.g. \nreporting all contigs between 0-10%, or other large ranges) can exaggerate certain trends or\nsquash smaller signals. To address this, we use raw counts in addition to percentages when\nreporting shared sequence in the manuscript, and all individual level sequences (plus the\ngenotype matrix) are available on Zenodo as of this submission.\n2. As mentioned in your comment, we are not working with variants of all lengths, but rather with\nonly larger contigs (>1Kb sequences). This may further distort the counts, as we have\nsubstantially fewer variants than if we also included small variants.\n3. There are several trios (85) in the 1KGP dataset, which will affect the number of singleton\nvariants detected in the data. We have addressed this in the manuscript, and the effect it may\nhave on the counts.\nPage 15, Line 3: “We note that the presence of trios distorts the frequency spectrum compared\nto expectations for an unrelated sample, but the practical effect is minimal in practice (e.g.,\nenriching for doubletons versus singletons).”\n4. There is also a “squashing” effect that happens through our approach to detecting shared\nsequence. We combine contigs that are 90% similar into the same sequence, which results in\nus losing singleton variants that are highly similar but still distinct. Removing this tolerance and\nrequiring an exact match will cause the number of singleton variants to rise.\n5. As shown in our analysis of the SGDP outliers (and in the PCA analysis), technical artifacts\ndo exist in this data, and this too may affect the distribution of variants.\nMany of the cited bioRxiv preprints have already been published, sometimes years ago, such as\nAganezov et al. (Science, 2022), Nurk et al. (Science, 2022), Taylor er al. (Nature, 2024) or\nGustafson et al. (Genome Res, 2024).\nWe have updated the citations to refer to the published versions of Nurk et al., Aganezov et al.,\nTaylor et al., Gustafson et al. and Stolyarova et al. We apologize for the earlier errors.\nMinor comments\nPlease clarify the term \"non-LOC\" gene. Figure S32 was illegible in my printout. Perhaps the\ngrayscale could be adjusted?\nThis refers to how genes are named and annotated by Gencode v35. When a published gene\nsymbol (e.g., BRCA2) is not available, and orthologs have not yet been determined, Genecode\nwill provide a symbol that is constructed as 'LOC' + the GeneID. We use the term “non-LOC” to\nindicate how we have removed these unnamed genes from our analyses, so that we primarily\nrestrict our analysis to well-annotated genes from Gencode v35.\nWe apologize for the issues with Figure S32. The goal of the figure was to highlight the few, but\nsignificant, associations present in the all-vs-all comparison. We hope that the full size figure will\nmake this clearer to see.\nReviewer #4 (Remarks on code availability):\nThe code base consists largely of WDL files with commands, which are also listed in the\nappendix. This mostly looks good to me, but I haven't tried these WDL workflows.\nThank you for your comment."}, {"round": 2, "reviewer_comments": {"reviewer_2": "I would like to thank the authors for carefully considering my previous comments, responding thoughtfully, and even\nperforming additional analyses. However, one important aspect still seems to be missing: the representativeness and global\npopulation context of the variation found in these unmapped reads.\nIn the discussion, the authors note that they compared these newly assembled sequences to sequences from non-SAS\npopulations to understand how the variation is distributed within and across populations. Yet, as they also note in their\nresponses to my earlier comments, there is no clear correlation between contig sharing and individual relatedness (e.g.,\nindividuals sharing more IBD do not tend to share more unmapped reads) or population structure (e.g., among 1000\nGenomes Project populations such as those from Sri Lanka, Bangladesh, and Pakistan). As the authors acknowledge, these\npatterns suggest that much of the observed variation may be driven by batch effects or other technical artifacts: consistent\nwith their finding that approximately 14% of the presence/absence variation in shared contigs can be explained by dataset\ndifferences.\nWhile this is not affecting the primary conclusions, it represents an important limitation that could be further quantified. For\nexample, assessing the heritability of unmapped reads between parents and children, as trios are available in the dataset,\ncould help evaluate how effectively current sequencing and analytical methods capture this variation. There will always be a\ntechnical component, but such an analysis could help determine whether the unmapped reads are being captured\nconsistently enough, or whether the detection remains too incomplete or too random to reveal meaningful population\nstructure.\nFinally, this limitation should be explicitly stated: as currently analyzed, these unmapped reads cannot yet be reliably used\nfor population genetics studies, since their relationship to population structure remains unclear. Consequently, they should\nalso not be used for GWAS or similar analyses at this stage.", "reviewer_3": "I thank the authors for posting the assembled sequences in this revision. My concerns have been addressed.", "reviewer_4": "The revised manuscript by Das et al. now includes an additional analysis of repeat content of the assembled contigs. This\nshows that contigs assembled from GRCh38 unmapped reads contain a substantially higher proportion of satellite repeats\ncompared to T2T-CHM13 (as expected) and a large proportion of mobile elements. An additional long-read study was also\nincorporated to verify the assembly results, which strengthens confidence in the reported contig set. The data generated in\nthis work have now also been deposited at Zenodo. All issues identified in my initial review have been addressed\nsatisfactorily. While some of my earlier comments could have been explored in greater depth (e.g., Mendelian errors in trios),\nI recognize that the previously mentioned lack of contig genotypes is a limitation for such analyses. I have no further major\ncomments, aside from a few minor errors and typographical issues noted during review:\n* Figure 3: Panel c is missing in the Figure caption\n* Figure 4: Inset in (a) is … should be (e), I think, as panel (a) has no inset.\n* ST8 – Repeat analysis: “ALIUs” should be ALUs"}, "author_response": "REVIEWER COMMENTS\nReviewer #2 (Remarks to the Author):\nI would like to thank the authors for carefully considering my previous comments, responding\nthoughtfully, and even performing additional analyses. However, one important aspect still\nseems to be missing: the representativeness and global population context of the variation\nfound in these unmapped reads.\nIn the discussion, the authors note that they compared these newly assembled sequences to\nsequences from non-SAS populations to understand how the variation is distributed within and\nacross populations. Yet, as they also note in their responses to my earlier comments, there is no\nclear correlation between contig sharing and individual relatedness (e.g., individuals sharing\nmore IBD do not tend to share more unmapped reads) or population structure (e.g., among\nAs the authors acknowledge, these patterns suggest that much of the observed variation may\nbe driven by batch effects or other technical artifacts: consistent with their finding that\napproximately 14% of the presence/absence variation in shared contigs can be explained by\ndataset differences.\nWhile this is not affecting the primary conclusions, it represents an important limitation that could\nbe further quantified. For example, assessing the heritability of unmapped reads between\nparents and children, as trios are available in the dataset, could help evaluate how effectively\ncurrent sequencing and analytical methods capture this variation. There will always be a\ntechnical component, but such an analysis could help determine whether the unmapped reads\nare being captured consistently enough, or whether the detection remains too incomplete or too\nrandom to reveal meaningful population structure.\nFinally, this limitation should be explicitly stated: as currently analyzed, these unmapped reads\ncannot yet be reliably used for population genetics studies, since their relationship to population\nstructure remains unclear. Consequently, they should also not be used for GWAS or similar\nanalyses at this stage.\nThank you for your comments.\nTo address your first concern and better understand the limitations of this approach, we\nperformed the following trio analysis.\nAcross the 85 1KGP trios in the set, we counted the number of contigs in the child that were\npresent in either or both parents, and counted the number of contigs present in each parent that\nwere also found in the child. We found that on average 92% of the child contigs were also found\nin one or both of their parents (distribution plotted below). This percentage provides strong\nsupport for what is expected for a child to inherit from its parents, with slight distortions caused\nby noise added during assembly, the selection of only large contigs, the collapsing of similar\ncontigs, and genuine de novo contigs, all of which will lower the fraction below 100% (see below\nfor a discussion of these factors).\nWe also found that on average 74% of the contigs found in parents were also found in their\nchildren. These numbers indicate that our approach, and the assembled sequences, do capture\nthe expected patterns of Mendelian inheritance (distributions plotted below for male and female\nparents), especially when considering the additional noise added through variance during\nassembly, selection of only large contigs, and the grouping of similar contigs that are not\nidentical.\nThe level of inheritance from the parent to the child is also inflated by the number of contigs\nshared between parents. To remove this effect, we also looked at how many contigs are present\nonly in one parent, that are inherited by the child. Doing this, we see that the fraction of contigs\ninherited drops to 57% on average. This fraction is in line with expectations for predominately\nheterozygous contigs in the parent, with slight distortions caused by noise added during\nassembly or during the collapsing of similar contigs (in addition to expected biological variance).\nWe also looked at the sharedness between the parents in each pair, as similarity here will\nfurther increase the number of parent contigs found in the child. We found that, on average,\n~50% of the contigs in each parent was also found in the other parent. This high degree of\nsimilarity will inflate the fraction of each parent’s contigs found in their child, as both parents are\ncontributing the same contigs. Note many of these will occur from the reference genome itself\ncarrying a rare deletion so that both parents will share a common “insertion”.\nThere are a few considerations to be made that could be distorting all these numbers. First,\nthere is noise and variance introduced during the assembly process, caused by both genuine\nsequencing variability and different choices made by the assembler (e.g. a small shift if read\nplacement could influence if a region can be assembled). Second, as noted in our previous\nresponses, we focus our analysis on large contigs (≥1 kbp), which can create thresholding\nissues when a child’s or parent’s contig is slightly shorter than this limit and thus not considered\nby this analysis. Finally, these numbers are taken from our final set of collapsed unique and\nshared contigs which uses a 90% similarity threshold to consider two contigs to be equivalent\nand collapsed together.\nTo evaluate the impact of this collapsing, we align the contigs from the individuals in 5 trios to\nthe other individuals in the same trio, and require 99% identity alignments between contigs to\ncall them the same contig. When using this threshold, we find that 83% of the contigs in the\nchild are found in one or both of their parents, and 56% of the contigs present in a parent are\nfound in their child. The similarity between parents also drops using this metric: on average, the\ntwo parents shared less than 20% of their contigs. The primary reason behind this drop is due to\na stricter threshold emphasizing the variation that occurs during the assembly process, with\ncontigs with similar but not nearly-identical subsequences no longer meeting the threshold for\nbeing called identical.\nThis analysis, and the corresponding plots, has been added to the supplementary section of the\nmanuscript, and a reference to this has been added to the main manuscript. The full size\nversion of these figures have also been added to the GitHub repository, along with all existing\nfigures in the paper.\nMain manuscript, “Shared Sequence Analysis” section, Page 15 Line 5:\n“Further analysis of the shared sequences between members of the trios can be found in the\nSupplementary Results (“Patterns of inheritance and similarity across contigs assembled in 85\n1KGP trios”) and Supplementary Methods.”\nSupplementary Results, “Patterns of inheritance and similarity across contigs assembled in 85\n1KGP trios” section, Page 2 Line 26:\n“Across the 85 1KGP trios in the set, we looked at the number of contigs in the child that were\npresent in either or both parents, and calculated the number of contigs present in each parent\nthat were also found in the child. We found that on average 92% of the child contigs were also\nfound in one or both of their parents (Supplementary Figure 34). This number is consistent with\nwhat we expect a child to inherit from its parents, with slight distortions caused by noise added\nduring assembly, the selection of only large contigs, or during the collapsing of similar contigs (in\naddition to genuine de novo contigs which will lower the fraction below 100%). For example, our\nchosen threshold of 1 Kbp for minimum contig length may result in a 999 base long contig from\nthe parent being discarded, while the corresponding 1,001 base contig remains in the child,\nlowering the fraction of inherited contigs.\nWe also found that on average 74% of the contigs found in parents were also found in their\nchildren (Supplementary Figure 35). This level of inheritance from parent to child is inflated by\nthe number of contigs shared between parents. To remove this effect, we also looked at how\nmany contigs are present only in one parent, that are inherited by the child. Doing this, we see\nthat the fraction of contigs inherited drops to 57% on average (Supplementary Figure 36). This\nnumber is consistent with predominately heterozygous contigs in the parents, with slight\ndistortions caused by noise added during assembly or during the collapsing of similar contigs (in\naddition to expected biological variance).\nTo quantify the extent to which interparent similarity affects both these counts, we also\nevaluated the sharedness between the parents in each pair. We found that, on average, 50% of\nthe contigs in each parent was found in the other parent. This high degree of similarity will\ninflate the fraction of each parent’s contigs found in their child, as both parents are contributing\nthe same contigs.\nAcross all these counts (child contigs in the parents, parent contigs in the child, interparent\nsimilarity), there are a few key factors that distort the observed counts from what we expect.\nFirst, there is noise and variance introduced during the assembly process, caused by both\ngenuine sequencing variability and due to different random choices made by the assembler.\nSecond, these are only large contigs, and do not represent the whole range of inherited\nvariation. Finally, these numbers are taken from our final set of collapsed unique and shared\ncontigs; during this collapsing, contigs that are 90% similar are called the same contig, resulting\nin contigs that are not identical being considered to be the same.\nTo evaluate the impact of the final factor, we align the contigs from individuals in 5 trios to\nthe other individuals in the same trio, and require 99% identity alignments between contigs to\ncall them the same contig. When we do this, we see that 83% of the contigs in the child are found in\none or both of their parents (down from 92%), and 56% of the contigs present in a parent are\nfound in their child (down from 74%). The similarity between parents also drops using this\nmetric: on average, the two parents shared less than 20% of their contigs.”\nSupplementary Methods, “Trio Analysis” section, Page 4 Line 20:\n“Trio contig analysis is carried out by comparing the collapsed contigs present in each member\nof the trio against the other two members. To calculate the fraction of contigs a child inherits\nfrom its parents, each of the child’s contigs are checked against the combined set of their parent\ncontigs. For the fraction of contigs a parent passes to a child, we perform the reverse and\ncompare a parent’s contigs against those found in the child. To remove the impact interparent\nsimilarity plays in this inheritance, we first compare the two sets of parent contigs, remove all\ncontigs shared by both parents, and then for each parent compare their private contigs against\nthe child set.\nFor our 99% identity analysis, we re-align the contigs from all members of a trio against each\nother, repeated across 5 trios, using the Minimap2 v2.22 setting specified in the main method\nsection (“Categorizing sequence as shared or unique”). We then count the number of contigs in\neach pairwise alignment that have identity >99%, and report that count over the total number of\ncontigs as the fraction of matched contigs.”\nWe have also added the following text to the manuscript, in the Discussion, to explicitly state the\nlimitations of this analysis:\nMain manuscript, Page 22, Line 1:\n“The approach presented here, and the contigs assembled through this approach, serve as a\nmethod to identify interesting candidate regions where variation absent in the reference genome\nmay occur. Given the lower sensitivity and specificity compared, for example, to conventional\nalignment based genotyping of SNVs, results from any application of these unmapped-read\nvariants for population genetic scans, GWAS, or similar analyses would need to be validated\nwith more targeted approaches. Nevertheless, the trio-based transmission analysis and other\nconfirmatory analyses in the study support the validity of most of these variants, and thus\nsupport their use for downstream exploration and discovery.”\nReviewer #3 (Remarks to the Author):\nI thank the authors for posting the assembled sequences in this revision. My concerns have\nbeen addressed.\nThank you for your comment.\nReviewer #4 (Remarks to the Author):\nThe revised manuscript by Das et al. now includes an additional analysis of repeat content of\nthe assembled contigs. This shows that contigs assembled from GRCh38 unmapped reads\ncontain a substantially higher proportion of satellite repeats compared to T2T-CHM13 (as\nexpected) and a large proportion of mobile elements. An additional long-read study was also\nincorporated to verify the assembly results, which strengthens confidence in the reported contig\nset. The data generated in this work have now also been deposited at Zenodo. All issues\nidentified in my initial review have been addressed satisfactorily. While some of my earlier\ncomments could have been explored in greater depth (e.g., Mendelian errors in trios), I\nrecognize that the previously mentioned lack of contig genotypes is a limitation for such\nanalyses. I have no further major comments, aside from a few minor errors and typographical\nissues noted during review:\n* Figure 3: Panel c is missing in the Figure caption\n* Figure 4: Inset in (a) is … should be (e), I think, as panel (a) has no inset.\n* ST8 – Repeat analysis: “ALIUs” should be ALUs\nThank you for your comments.\nWe have addressed the minor errors and typographical issues you mentioned. We have made\nthe following changes:\n● Added “(c) Counts of the number of unique gene elements in T2T-CHM13 v2.0\nintersecting placed contigs.” to the Figure 3 caption.\n● Amended the caption to Figure 4 to be “Inset in (e) is the” instead of (a)\n● Corrected “ALIUs” to “ALUs” in ST8. This will also be reflected in the Zenodo upload.\nReviewer #4 (Remarks on code availability):\nThe GitHub repository looked fine to me already in the first review.\nThank you!"}, {"round": 3, "reviewer_comments": {"reviewer_2": "I thank the authors for thoughtfully addressing my comments and for performing additional analyses.\nMy concerns have been addressed."}, "author_response": ""}]

Decision letter