ga4gh / ga4gh/refget

Fingerprinting sequences using gaps

Open
#29 3 comments 0 reactions 0 assignees View on GitHub
refget-seqcol schema-term
Dominant language
No language data
Stars
17
Forks
10
Avg merge
11m
Merged PRs (30d)
1

Description

Following on from the previous discussions today I spent some time writing up a [bit of code to find gaps in a FASTA file](https://gist.github.com/andrewyatz/652020eb5af7e731c4b2e0db2d42e555). 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](https://hgdownload.cse.ucsc.edu/goldenpath/hg38/bigZips/p13/hg38.p13.2bit) and the second from [Ensembl's FTP site](http://ftp.ensembl.org/pub/current_fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.chromosome.22.fa.gz).

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

Contributor guide

Open the contributing guide

Assessment

This issue has not been assessed yet.

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.