Fingerprinting sequences using gaps

[andrewyatz]

Following on from the previous discussions today I spent some time writing up a bit of code to find gaps in a FASTA file. It's in Python so you don't need to worry about Perl. I ran this over chromosome 22 from UCSC and Ensembl. I also verified this against the [AGP file from UCSC](wget https://hgdownload.cse.ucsc.edu/goldenpath/hg38/bigZips/hg38.agp.gz).

$ curl -s https://hgdownload.cse.ucsc.edu/goldenpath/hg38/bigZips/hg38.agp.gz  | gzip -dc | grep --regexp="\t[NU]\t" | grep "chr22\t"
chr22   1   10000   1   N   10000   telomere    no  na
chr22   10001   10510000    2   N   10500000    short_arm   no  na
chr22   10784644    10834643    5   N   50000   contig  no  na
chr22   10874573    10924572    7   N   50000   contig  no  na
chr22   10966725    11016724    9   N   50000   contig  no  na
chr22   11068988    11118987    12  N   50000   contig  no  na
chr22   11160922    11210921    14  N   50000   contig  no  na
chr22   11378057    11428056    16  N   50000   contig  no  na
chr22   11497338    11547337    19  N   50000   contig  no  na
chr22   11631289    11681288    23  N   50000   contig  no  na
chr22   11724630    11774629    25  N   50000   contig  no  na
chr22   11977556    12027555    27  N   50000   contig  no  na
chr22   12225589    12275588    29  N   50000   contig  no  na
chr22   12438691    12488690    31  N   50000   contig  no  na
chr22   12641731    12691730    33  N   50000   contig  no  na
chr22   12726205    12776204    35  N   50000   contig  no  na
chr22   12818138    12868137    37  N   50000   contig  no  na
chr22   12904789    12954788    39  N   50000   contig  no  na
chr22   12977326    12977425    41  N   100 contig  no  na
chr22   12986172    12994027    43  N   7856    scaffold    yes paired-ends
chr22   13011654    13014130    45  N   2477    scaffold    yes paired-ends
chr22   13021323    13021422    47  N   100 contig  no  na
chr22   13109445    13109544    49  N   100 contig  no  na
chr22   13163678    13163777    51  N   100 contig  no  na
chr22   13227313    13227412    53  N   100 contig  no  na
chr22   13248083    13248182    55  N   100 contig  no  na
chr22   13254853    13254952    57  N   100 contig  no  na
chr22   13258198    13258297    59  N   100 contig  no  na
chr22   13280859    13280958    61  N   100 contig  no  na
chr22   13285144    13285243    63  N   100 contig  no  na
chr22   14419455    14419554    65  N   100 contig  no  na
chr22   14419895    14419994    67  N   100 contig  no  na
chr22   14420335    14420434    69  N   100 contig  no  na
chr22   14421633    14421732    71  N   100 contig  no  na
chr22   15054319    15154318    73  N   100000  contig  no  na
chr22   16279673    16302843    99  N   23171   scaffold    yes unspecified
chr22   16304297    16305427    101 N   1131    scaffold    yes paired-ends
chr22   16307049    16307605    103 N   557 scaffold    yes paired-ends
chr22   16310303    16310402    105 N   100 scaffold    yes paired-ends
chr22   16313517    16314010    107 N   494 scaffold    yes paired-ends
chr22   18239130    18339129    142 N   100000  contig  no  na
chr22   18433514    18483513    144 N   50000   contig  no  na
chr22   18659565    18709564    146 N   50000   contig  no  na
chr22   49973866    49975365    1119    N   1500    contig  no  na
chr22   50808469    50818468    1153    N   10000   telomere    no  na

I then ran the gap finder code on two files. One from UCSC which had to be pulled out using their twoBitToFa tool and the second from Ensembl's FTP site.

Both representations of chromosome 22 produced a gap fingerprint identical to the AGP file above.

The biggest issue with this method is if there are no gaps in the underlying sequence, then this cannot be used to define a meaningful fingerprint. As said on the call the rise of full length assemblies and improvements in sequencing methods will mean there is limited value here but there is an argument that says this is useful.

I would suggest that if we take on this idea that it's an optional extension

[sveinugu]

Thanks, @andrewyatz, for your investigations! I see you base the algorithm on a cutoff (10 bases) to determine whether a gap is to be considered consecutive. Although this produces the same result as the AGP file in your test case, at least theoretically it might fail in other cases. For such an extension, why couldn't one just make directly use of the AGP file instead of basing the calculations on the FASTA sequence?

I agree that this should be considered an optional extension. I raised the idea earlier without getting much feedback, if I recall correctly, so I assumed this would be out of scope of the standard. But I do think it really fits well together with other arrays we have been considering, such as alphabet and topology.

[andrewyatz]

The reason I was trying to investigate this method was because if we can base another metric on the sequence content rather than another out of band file or data source there would be an algorithm capable of doing this for any sequence.

As an update I just redid this. Worryingly yes the gaps did start to differ. This is a shame. It seems though either my code is broken badly (very possible) or the edits in each of these sequences is too much to recreate the gaps found in the AGP file