Merge branch 'development':

Major overhaul of the pipeline according to the changes and improvements
made since the bioRxiv preprint.

.gitignore

@@ -1,14 +1,16 @@
 __pycache__
 .snakemake
 .ipynb_checkpoints
+.sourcetrail
 data
 tools
 sandbox
 jupyter
 *.log
 *.aux
-publications/astronaut-telomeres-paper
-publications/methods-paper/figures/*-.png
-publications/methods-paper/figures/*/*-.png
-publications/methods-paper/figures/*.svg
-publications/methods-paper/figures/*/*.svg
+assets/paper/figures/*.png
+assets/paper/figures/*/*.png
+assets/paper/figures/*.svg
+assets/paper/figures/*/*.svg
+*.odt
+Snakefile

tools/generate-hg38ext.py → assets/generate-hg38ext.py

@@ -7,16 +7,31 @@
 from contextlib import contextmanager
 from binascii import hexlify
 from gzip import open as gzopen
-from tqdm import tqdm
 from itertools import chain
 from textwrap import fill
+from zlib import decompress
+from base64 import decodebytes
+
+try:
+    from tqdm import tqdm
+except ModuleNotFoundError:
+    def tqdm(it, desc=None, *args, **kwargs):
+        if desc is not None:
+            print(desc, end="...\n", file=stderr)
+        return it
 
 
 USAGE = """usage:
 {0} --local hg38.fasta stong2014.fasta
-    generate hg38ext.fa from local files
+    generate hg38ext.fa from local files and output to stdout
 {0} --remote
-    download appropriate assemblies and generate hg38ext.fa
+    download appropriate assemblies, generate hg38ext.fa, and output to stdout
+{0} --ecx
+    output the edgeCase indeX (hg38ext.fa.ecx) to stdout
+
+NOTE! This tool writes uncompressed data (FASTA or ECX) to stdout.
+You should pipe it into a file, for example:
+{0} --remote > hg38ext.fa
 """
 
 NCBI_FTP_DIR = "ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/405"
@@ -49,6 +64,39 @@
 COMPLEMENTS = dict(zip(ALPHABET, reversed(ALPHABET)))
 COMPLEMENT_PATTERN = compile(r'|'.join(COMPLEMENTS.keys()))
 
+HG38EXT_ECX_GZ = b'\n'.join([
+b'eJyVWEtvIzcMPk/+Ro9FAFFPsos9F8X2WmD3ZLhGdhMkcba2u2j/fUk9JvOwOG6AUMSY/CRqOB8p'
+b'/fT1Zf/t/Mvr/vy82x8Pj2+nj95Q/PD17fT8EYHsh8tpf7i03yA69Hc/PRwvp3+H03H/+jB8fzvL'
+b'/88wHB5Pb69v5zd++Lp/Ou7K7zLB8P30xOrhZX8+D+e3v0+Hh+Hl6fg8/PmyPzy/PJ0vd0bcYQDD'
+b'f1lCeZCFLGkIw2SVw+M3h8P9cH+XTezS0RahOWYTt3R0RWiO2cQvHX0RmmM2CUvHUITmmE3iELNj'
+b'bI6xCM0xm6TljKkIzTGb4HJGLEJzzCa0nJGK0ByzCZhVBpgq1RwoRrBcLzSpOhejVQqBrVJ1Lkac'
+b'Rnlq8S8K1OdFqhDFyHcgfJUqRDHi1EoNIk0gQpUqRDGKq12IVarOxSit9j9VqToXI1zNjFWqztnI'
+b'muXM1lSpOuecszAEU+YuI9SnRaoUUmaR1Alt9eEdwlapQmTTz8vgPxehemb4L0vPL0WoniXjBuuR'
+b'fPQM1DS3JF3XQfCFX623gC5YapozS/btIYRKv4TWhRCoadEsabiHECsPk7HGM8SoxSUf9xAqEUIg'
+b'5wIlN2p+SZM9hMqI4INJGF0ctbTkyx5CpUZw6JxP/DE0DZfE2UFwjSHBuYRJXmfT3JpBeyCNKcHx'
+b'+mPMIEVzaybtgTTG5Plt5FygptU1zBi1B9I4E8C74B1vcdNozak9kMaaYCKiS5ldikZrVu2BNN4E'
+b'w7kJgNQ0Mmte7YE0akRnffIemmLXBNmDaASJxkbej9CUuKbJHkQxoiFg5JfKfF2VVJ8XqUJQY1rP'
+b'6QkxNaUx8IRvOxDeNE6NkYjQNcWvGbcH0Zg1GDToefKq0JpxexC2Mi5TjTNIYdTiknl7CK5wbEgW'
+b'kucEqEpc8m/Pv2Qfv4NAiGWgCfPmLpvJ+/T0cHl83R93r/vL4bH5Px0f/uGW/M6HwrLAyVCkm3Cv'
+b'jsFriLUZBTQpS5w0o5veSetIb1w/tubUurz+8uU3Tr4Ng7Q+9TaMYFrLyuUrSz9h3q2dCK1LBeNB'
+b'QLxtfFufbwHYjd71xijc2MWGFLKMM87dXId/bwBtPjKU4jY2gJsAYezjYibNCHOC2gSIlUMScQ0W'
+b'GRurFHnjTqR5Z4GTHuO9syhYbgsLxx4DrA3oRs1PeowbsWjabeDYbeCs29Cx+ERnJh2Hk1xtWlhl'
+b'voIyL/chUtOSuZJ8CtCs5HtiamtauJJ+CpCbFluDQE2z82K7CeSnpc5TbMq81G3ChLHccaFNzI5V'
+b'oSspufXqY5wWPnnhVQmzwre5pjQpfkkqV1XirPhtwuB7AbRceWDU7KQA3hhZ2U+7+/V3TNxzhh+w'
+b'279cBk+B6Vykn6WBXCnxFyzXSwWK8+3+LrX2kD9SbmdFgu7G87JbPSHvPv1mkzQBdXLk75NPKHkI'
+b'03ZvZSiwbroaK7D2vXnkz1P4uGl2inbFHcQ9Z/IfxzJZ8rIlpwPvsWdgkW7W+a13xAuIrz/3TtRr'
+b'NyduYXDfLw8vO7gPxnzawY5fDe8fJ+qK+8IghhmgvWNrOKp8/ROv4/B2lmNXRaqruA6UkiA1NoYo'
+b'xwGIYBT/EkkUP25olUhmzKtGQtdxUkBMWeIET40HJVXRVEuucS5Lr/iXeCRXkZMolYXsrgU0vtjt'
+b'kNB2oPjAVeXVFOtEJXmDzQE8OKEQGXSQHBrK94JSH5V3teBcNbTQgUp8dBOBM0A9MEmjekai1XfU'
+b'gyhhBXHl3fhrvRT5lj3JH7NvBl027k4LELuglijKJybD9B6mskwnyJxb9fBgPRMwJLKjpuGUSFGu'
+b'Xs1gb410drzQIuWy0gWNSPw684ATUDVSklSjaim3TYGDjKOWFJwcKQk384bQrZHODgD9SFslJN+H'
+b'dk7OuHkIE2g9XslBqubgxJbknFc1p+CUeKWMkLTzNyfx/Miivtyk4MqlEZTBznD1eHMmYluCczE3'
+b'FKOGKlQJWQoOhwT+5pA7jcH2DoDcufbmcWQSl8o8xBsakOvTcP+bJ+p0OWDZmaMuI/2fafJuSQQC'
+b'L9XkShw/Ojs2r1P6FjkNmk/bIAgyhHWp6e+Kz9itQUKud8ZJj1iUK1VrHbnLEFxs8OZMmZ9i9bij'
+b'BsxHrFCG+fF4I+qUkSdXgXwa5aVVBVSoGnXMEDiEW4Oe3QbpIZOCmuQ95yFMUPV4IScn1GUAgk/S'
+b'z48aKUg1XMoIMMRbw51dX6nhcrr1UflXWwY3Qd0IN2ckVAfI5wbiHrhppCCVcMFmBOmJb87p+dWO'
+b'HnFQgJEMl+I8pBnwRtA5IUdjMtai874pQYWqUUuV/A80XjlL', b''])
+
 
 def revcomp(sequence):
     """Reverse-complement a sequence"""
@@ -112,7 +160,7 @@ def parser_iterator(filename, to_revcomp, desc="Parsing", entry_filter=lambda e:
 
 
 def generate_hg38ext(hg38, stong2014):
-    """Generate hg38ext from the hg38 and stong2014 FASTA files"""
+    """Generate hg38ext from the hg38 and stong2014 FASTA files, write to stdout"""
     subhaps = {mask.format(STONG_INFIX) for mask in STONG2014_SUBHAP_MASKS}
     to_revcomp = {mask.format(STONG_INFIX) for mask in TO_REVCOMP_MASKS}
     hg38_iterator = parser_iterator(hg38, to_revcomp, desc="Parsing reference")
@@ -125,6 +173,11 @@ def generate_hg38ext(hg38, stong2014):
     return 0
 
 
+def output_ecx():
+    """Write the hg38ext ECX to stdout"""
+    print(decompress(decodebytes(HG38EXT_ECX_GZ)).decode().rstrip("\n"))
+
+
 if __name__ == "__main__":
     # interpret command-line arguments and dispatch to subroutines:
     if (len(argv) == 2) and (argv[1] == "--remote"):
@@ -133,6 +186,8 @@ def generate_hg38ext(hg38, stong2014):
             returncode = generate_hg38ext(hg38, stong2014)
     elif (len(argv) == 4) and (argv[1] == "--local"):
         returncode = generate_hg38ext(hg38=argv[2], stong2014=argv[3])
+    elif (len(argv) == 2) and (argv[1] == "--ecx"):
+        returncode = output_ecx()
     else:
         print(USAGE.format(__file__).rstrip(), file=stderr)
         returncode = 1

assets/hg38ext.fa.ecx

@@ -1,74 +1,74 @@
-#flags:ucsc_mask_anchor=4096;fork=8192;tract_anchor=16384
-#entry	rname	pos	pos+1	chromosome	ucsc_rname	flag	prime	class	source	link	blacklist
-0	chr1	10000	10001	chr1	chr1	4096	5	ucsc_mask_anchor	hg38	-	-
-1	chr2	10000	10001	chr2	chr2	4096	5	ucsc_mask_anchor	hg38	-	-
-2	chr3	10000	10001	chr3	chr3	4096	5	ucsc_mask_anchor	hg38	-	-
-3	chr4	10000	10001	chr4	chr4	4096	5	ucsc_mask_anchor	hg38	-	-
-4	chr5	10000	10001	chr5	chr5	4096	5	ucsc_mask_anchor	hg38	-	-
-5	chr6	60000	60001	chr6	chr6	4096	5	ucsc_mask_anchor	hg38	-	-
-6	chr7	10000	10001	chr7	chr7	4096	5	ucsc_mask_anchor	hg38	-	-
-7	chr8	60000	60001	chr8	chr8	4096	5	ucsc_mask_anchor	hg38	-	-
-8	chr9	10000	10001	chr9	chr9	4096	5	ucsc_mask_anchor	hg38	-	-
-9	chr10	10000	10001	chr10	chr10	4096	5	ucsc_mask_anchor	hg38	-	-
-10	chr11	60000	60001	chr11	chr11	4096	5	ucsc_mask_anchor	hg38	-	-
-11	chr12	10000	10001	chr12	chr12	4096	5	ucsc_mask_anchor	hg38	-	-
-12	chr13	16000000	16000001	chr13	chr13	4096	5	ucsc_mask_anchor	hg38	-	-
-13	chr14	16000000	16000001	chr14	chr14	4096	5	ucsc_mask_anchor	hg38	-	-
-14	chr15	17000000	17000001	chr15	chr15	4096	5	ucsc_mask_anchor	hg38	-	-
-15	chr16	10000	10001	chr16	chr16	4096	5	ucsc_mask_anchor	hg38	-	-
-16	chr17	60000	60001	chr17	chr17	4096	5	ucsc_mask_anchor	hg38	-	-
-17	chr18	10000	10001	chr18	chr18	4096	5	ucsc_mask_anchor	hg38	-	-
-18	chr20	60000	60001	chr20	chr20	4096	5	ucsc_mask_anchor	hg38	-	-
-19	chr21	5010000	5010001	chr21	chr21	4096	5	ucsc_mask_anchor	hg38	-	-
-20	chr22	10510000	10510001	chr22	chr22	4096	5	ucsc_mask_anchor	hg38	-	-
-21	chrX	10000	10001	chrX	chrX	4096	5	ucsc_mask_anchor	hg38	-	-
-22	chrY	10000	10001	chrY	chrY	4096	5	ucsc_mask_anchor	hg38	-	-
-23	chr1	248946422	248946423	chr1	chr1	4096	3	ucsc_mask_anchor	hg38	-	-
-24	chr2	242183529	242183530	chr2	chr2	4096	3	ucsc_mask_anchor	hg38	-	-
-25	chr3	198235559	198235560	chr3	chr3	4096	3	ucsc_mask_anchor	hg38	-	-
-26	chr4	190204555	190204556	chr4	chr4	4096	3	ucsc_mask_anchor	hg38	-	-
-27	chr7	159335973	159335974	chr7	chr7	4096	3	ucsc_mask_anchor	hg38	-	-
-28	chr8	145078636	145078637	chr8	chr8	4096	3	ucsc_mask_anchor	hg38	-	-
-29	chr9	138334717	138334718	chr9	chr9	4096	3	ucsc_mask_anchor	hg38	-	-
-30	chr10	133787422	133787423	chr10	chr10	4096	3	ucsc_mask_anchor	hg38	-	-
-31	chr11	135076622	135076623	chr11	chr11	4096	3	ucsc_mask_anchor	hg38	-	-
-32	chr12	133265309	133265310	chr12	chr12	4096	3	ucsc_mask_anchor	hg38	-	-
-33	chr13	114354328	114354329	chr13	chr13	4096	3	ucsc_mask_anchor	hg38	-	-
-34	chr14	106883718	106883719	chr14	chr14	4096	3	ucsc_mask_anchor	hg38	-	-
-35	chr15	101981189	101981190	chr15	chr15	4096	3	ucsc_mask_anchor	hg38	-	-
-36	chr17	83247441	83247442	chr17	chr17	4096	3	ucsc_mask_anchor	hg38	-	-
-37	chr18	80263285	80263286	chr18	chr18	4096	3	ucsc_mask_anchor	hg38	-	-
-38	chr19	58607616	58607617	chr19	chr19	4096	3	ucsc_mask_anchor	hg38	-	-
-39	chr20	64334167	64334168	chr20	chr20	4096	3	ucsc_mask_anchor	hg38	-	-
-40	chr21	46699983	46699984	chr21	chr21	4096	3	ucsc_mask_anchor	hg38	-	-
-41	chr22	50808468	50808469	chr22	chr22	4096	3	ucsc_mask_anchor	hg38	-	-
-42	chrX	156030895	156030896	chrX	chrX	4096	3	ucsc_mask_anchor	hg38	-	-
-43	chrY	57217415	57217416	chrY	chrY	4096	3	ucsc_mask_anchor	hg38	-	-
-44	chr1	585988	585989	chr1	chr1	16384	5	riethman_match	hg38	-	dup?
+#flags:mask_anchor=4096;fork=8192;tract_anchor=16384
+#entry	rname	pos	pos+1	chromosome	main_rname	flag	prime	class	source	link	blacklist
+0	chr1	10000	10001	chr1	chr1	4096	5	mask_anchor	hg38	-	-
+1	chr2	10000	10001	chr2	chr2	4096	5	mask_anchor	hg38	-	-
+2	chr3	10000	10001	chr3	chr3	4096	5	mask_anchor	hg38	-	-
+3	chr4	10000	10001	chr4	chr4	4096	5	mask_anchor	hg38	-	-
+4	chr5	10000	10001	chr5	chr5	4096	5	mask_anchor	hg38	-	-
+5	chr6	60000	60001	chr6	chr6	4096	5	mask_anchor	hg38	-	-
+6	chr7	10000	10001	chr7	chr7	4096	5	mask_anchor	hg38	-	-
+7	chr8	60000	60001	chr8	chr8	4096	5	mask_anchor	hg38	-	-
+8	chr9	10000	10001	chr9	chr9	4096	5	mask_anchor	hg38	-	-
+9	chr10	10000	10001	chr10	chr10	4096	5	mask_anchor	hg38	-	-
+10	chr11	60000	60001	chr11	chr11	4096	5	mask_anchor	hg38	-	-
+11	chr12	10000	10001	chr12	chr12	4096	5	mask_anchor	hg38	-	-
+12	chr13	16000000	16000001	chr13	chr13	4096	5	mask_anchor	hg38	-	-
+13	chr14	16000000	16000001	chr14	chr14	4096	5	mask_anchor	hg38	-	-
+14	chr15	17000000	17000001	chr15	chr15	4096	5	mask_anchor	hg38	-	-
+15	chr16	10000	10001	chr16	chr16	4096	5	mask_anchor	hg38	-	-
+16	chr17	60000	60001	chr17	chr17	4096	5	mask_anchor	hg38	-	-
+17	chr18	10000	10001	chr18	chr18	4096	5	mask_anchor	hg38	-	-
+18	chr20	60000	60001	chr20	chr20	4096	5	mask_anchor	hg38	-	-
+19	chr21	5010000	5010001	chr21	chr21	4096	5	mask_anchor	hg38	-	-
+20	chr22	10510000	10510001	chr22	chr22	4096	5	mask_anchor	hg38	-	-
+21	chrX	10000	10001	chrX	chrX	4096	5	mask_anchor	hg38	-	-
+22	chrY	10000	10001	chrY	chrY	4096	5	mask_anchor	hg38	-	-
+23	chr1	248946422	248946423	chr1	chr1	4096	3	mask_anchor	hg38	-	-
+24	chr2	242183529	242183530	chr2	chr2	4096	3	mask_anchor	hg38	-	-
+25	chr3	198235559	198235560	chr3	chr3	4096	3	mask_anchor	hg38	-	-
+26	chr4	190204555	190204556	chr4	chr4	4096	3	mask_anchor	hg38	-	-
+27	chr7	159335973	159335974	chr7	chr7	4096	3	mask_anchor	hg38	-	-
+28	chr8	145078636	145078637	chr8	chr8	4096	3	mask_anchor	hg38	-	-
+29	chr9	138334717	138334718	chr9	chr9	4096	3	mask_anchor	hg38	-	-
+30	chr10	133787422	133787423	chr10	chr10	4096	3	mask_anchor	hg38	-	-
+31	chr11	135076622	135076623	chr11	chr11	4096	3	mask_anchor	hg38	-	-
+32	chr12	133265309	133265310	chr12	chr12	4096	3	mask_anchor	hg38	-	-
+33	chr13	114354328	114354329	chr13	chr13	4096	3	mask_anchor	hg38	-	-
+34	chr14	106883718	106883719	chr14	chr14	4096	3	mask_anchor	hg38	-	-
+35	chr15	101981189	101981190	chr15	chr15	4096	3	mask_anchor	hg38	-	-
+36	chr17	83247441	83247442	chr17	chr17	4096	3	mask_anchor	hg38	-	-
+37	chr18	80263285	80263286	chr18	chr18	4096	3	mask_anchor	hg38	-	-
+38	chr19	58607616	58607617	chr19	chr19	4096	3	mask_anchor	hg38	-	-
+39	chr20	64334167	64334168	chr20	chr20	4096	3	mask_anchor	hg38	-	-
+40	chr21	46699983	46699984	chr21	chr21	4096	3	mask_anchor	hg38	-	-
+41	chr22	50808468	50808469	chr22	chr22	4096	3	mask_anchor	hg38	-	-
+42	chrX	156030895	156030896	chrX	chrX	4096	3	mask_anchor	hg38	-	-
+43	chrY	57217415	57217416	chrY	chrY	4096	3	mask_anchor	hg38	-	-
+44	chr1	585988	585989	chr1	chr1	16384	5	riethman_match	hg38	-	inexact
 45	chr2	10262	10263	chr2	chr2	16384	5	riethman_match	hg38	-	-
 46	chr5	11807	11808	chr5	chr5	16384	5	riethman_match	hg38	-	-
-47	chr6	60000	60001	chr6	chr6	16384	5	riethman_match	hg38	-	-
-48	chr7	10232	10233	chr7	chr7	16384	5	riethman_match	hg38	-	-
-49	chr8	60000	60001	chr8	chr8	16384	5	riethman_match	hg38	-	-
+47	chr6	60000	60001	chr6	chr6	16384	5	riethman_match	hg38	-	inexact
+48	chr7	10232	10233	chr7	chr7	16384	5	riethman_match	hg38	-	inexact
+49	chr8	60000	60001	chr8	chr8	16384	5	riethman_match	hg38	-	inexact
 50	chr9	10353	10354	chr9	chr9	16384	5	riethman_match	hg38	-	-
 51	chr10	10419	10420	chr10	chr10	16384	5	riethman_match	hg38	-	-
-52	chr11	60000	60001	chr11	chr11	16384	5	riethman_match	hg38	-	-
+52	chr11	60000	60001	chr11	chr11	16384	5	riethman_match	hg38	-	inexact
 53	chr12	10575	10576	chr12	chr12	16384	5	riethman_match	hg38	-	-
 54	chr16	10027	10028	chr16	chr16	16384	5	riethman_match	hg38	-	-
 55	chr18	10615	10616	chr18	chr18	16384	5	riethman_match	hg38	-	-
-56	chr20	79359	79360	chr20	chr20	16384	5	riethman_match	hg38	-	dup?
-57	chr3	198235558	198235559	chr3	chr3	16384	3	riethman_match	hg38	-	-
-58	chr4	190122583	190122584	chr4	chr4	16384	3	riethman_match	hg38	-	dup?
+56	chr20	79359	79360	chr20	chr20	16384	5	riethman_match	hg38	-	inexact
+57	chr3	198235558	198235559	chr3	chr3	16384	3	riethman_match	hg38	-	inexact
+58	chr4	190122583	190122584	chr4	chr4	16384	3	riethman_match	hg38	-	inexact
 59	chr7	159335873	159335874	chr7	chr7	16384	3	riethman_match	hg38	-	-
 60	chr8	145073354	145073355	chr8	chr8	16384	3	riethman_match	hg38	-	-
 61	chr11	135076569	135076570	chr11	chr11	16384	3	riethman_match	hg38	-	-
 62	chr12	133264944	133264945	chr12	chr12	16384	3	riethman_match	hg38	-	-
 63	chr15	101980819	101980820	chr15	chr15	16384	3	riethman_match	hg38	-	-
 64	chr19	58607496	58607497	chr19	chr19	16384	3	riethman_match	hg38	-	-
-65	chr20	64286708	64286709	chr20	chr20	16384	3	riethman_match	hg38	-	dup?
+65	chr20	64286708	64286709	chr20	chr20	16384	3	riethman_match	hg38	-	inexact
 66	chr21	46699874	46699875	chr21	chr21	16384	3	riethman_match	hg38	-	-
 67	chr22	50807895	50807896	chr22	chr22	16384	3	riethman_match	hg38	-	-
-68	chrX	156029891	156029892	chrX	chrX	16384	3	riethman_match	hg38	-	-
+68	chrX	156029891	156029892	chrX	chrX	16384	3	riethman_match	hg38	-	inexact
 69	chr12_GL877875v1_alt	49533	49534	chr12	chr12	8192	5	fork	hg38	70	-
 70	chr12	55530	55531	chr12	chr12	8192	5	fork	hg38	69	-
 71	chr14_KI270846v1_alt	825824	825825	chr14	chr14_KI270846v1_alt	8192	3	fork	hg38	72	-
@@ -93,7 +93,7 @@
 90	2qtel_1-500K_1_12_12_rc	499999	500000	chr2	chr2	16384	3	tel_fork	riethman2014	-	-
 91	2qtel_1-500K_1_12_12_rc	468937	468938	chr2	chr2	8192	3	tel_fork	riethman2014	92	-
 92	chr2	242152486	242152487	chr2	chr2	8192	3	tel_fork	hg38	91	-
-93	9qtel_1-500K_1_12_12_rc	499999	500000	chr9	chr9	16384	3	tel_fork	riethman2014	-	-
+93	9qtel_1-500K_1_12_12_rc	499999	500000	chr9	chr9	16384	3	tel_fork	riethman2014	-	inexact
 94	9qtel_1-500K_1_12_12_rc	433984	433985	chr9	chr9	8192	3	tel_fork	riethman2014	95	-
 95	chr9	138192962	138192963	chr9	chr9	8192	3	tel_fork	hg38	94	-
 96	10qtel_1-500K_1_12_12_rc	499999	500000	chr10	chr10	16384	3	tel_fork	riethman2014	-	-

publications/methods-paper/COPYRIGHT → assets/paper/COPYRIGHT

@@ -6,5 +6,5 @@ unless otherwise noted, are protected by U.S. and International copyright laws.
 Reproduction and distribution, with or without modification, of these files
 without a written permission of the authors is prohibited.
 
-© 2020 Kirill Grigorev, Jonathan Foox, Christopher E. Mason
+© 2021 Kirill Grigorev, Jonathan Foox, Christopher E. Mason
 Institute for Computational Biomedicine, Weill Cornell Medicine

assets/paper/Snakefile

@@ -0,0 +1,82 @@
+from pandas import read_fwf
+from io import StringIO
+
+
+HG38EXT_ECX = "data/references/hg38/hg38ext.fa.ecx"
+DATA_DIR = "data/datasets/2021"
+MIN_MAP_OVERLAP = 500
+MIN_SUBTELOMERE_OVERLAP = 3000
+MIN_TELOMERE_OVERLAP = 3000
+MAX_READ_LENGTH = 100000
+TARGET = "tract_anchor"
+N_MOTIFS_TO_PLOT = 3
+MIN_CHROM_COVERAGE = 25
+SMALLEST_P_VALUE = 5e-324
+
+DATASETS = read_fwf(StringIO(str.strip("""
+group           subject  dataset      priority
+NA12878         HG001    11kb         1
+AshkenazimTrio  HG002    10kb         1
+AshkenazimTrio  HG002    15kb         1
+AshkenazimTrio  HG002    15kb_20kb    1
+AshkenazimTrio  HG003    15kb         1
+AshkenazimTrio  HG003    15kb_20kb    2
+AshkenazimTrio  HG004    15kb         1
+AshkenazimTrio  HG004    15kb_21kb    1
+ChineseTrio     HG005    11kb         1
+ChineseTrio     HG006    15kb_20kb    1
+ChineseTrio     HG006    hifi_google  1
+ChineseTrio     HG007    15kb_20kb    1
+ChineseTrio     HG007    hifi_google  1
+""")))
+
+"""Fully cannibalized datasets:
+AshkenazimTrio  HG004    15kb_20kb
+"""
+
+DATASETS["subject_pacbio_path"] = DATASETS.apply(
+    lambda row: "{}/PacBio/{}/{}".format(DATA_DIR, *row[:2]), axis=1,
+)
+DATASETS["dataset_pacbio_path"] = DATASETS.apply(
+    lambda row: "{}/PacBio/{}/{}/{}".format(DATA_DIR, *row[:3]), axis=1,
+)
+
+wildcard_constraints:
+    group="[^/]+", subject="[^/]+", dataset="[^/]+", name="[^/]+", kind="[^/]+",
+    arm="[pq]_arm",
+
+
+def get_sam_flags(arm, target=None):
+    if target:
+        if arm == "p_arm":
+            return "-f '{}' -F is_q".format(target)
+        elif arm == "q_arm":
+            return "-f is_q -f '{}'".format(target)
+        else:
+            raise ValueError("arm", arm)
+    else:
+        if arm == "p_arm":
+            return "-F is_q"
+        elif arm == "q_arm":
+            return "-f is_q"
+        else:
+            raise ValueError("arm", arm)
+
+
+include: "assets/paper/snakefiles/longread-motifs.snake"
+include: "assets/paper/snakefiles/shortread-support.snake"
+include: "assets/paper/snakefiles/shortread-motifs.snake"
+include: "assets/paper/snakefiles/bonferroni.snake"
+include: "assets/paper/snakefiles/densityplots.snake"
+include: "assets/paper/snakefiles/kmerscanner-all.snake"
+include: "assets/paper/snakefiles/levenshtein.snake"
+
+
+rule all:
+    input:
+        rules.densityplot_all.input,
+        rules.telbam_support_all.input,
+        rules.kmerscanner_all_motifs_all_subjects.input,
+        rules.kmerscanner_for_haploplots.input,
+        rules.levenshtein_all.input,
+        rules.repeatfinder_all_shortread.input,

assets/paper/figures/Figure_2.tex

@@ -0,0 +1,15 @@
+\documentclass{article}
+\usepackage[paperheight=840pt,paperwidth=700pt,margin=0in]{geometry}
+\usepackage[sfdefault]{roboto}
+\usepackage{graphicx}
+\usepackage{tikz}
+\usepackage[absolute,overlay]{textpos}
+    \setlength{\TPHorizModule}{1in}
+    \setlength{\TPVertModule}{1in}
+\begin{document}
+\begin{textblock}{13}(0.300,0.000)\includegraphics{Figure_2/densityplot-p_arm.pdf}\end{textblock}
+\begin{textblock}{13}(-0.05,0.200)\LARGE{(A)}\end{textblock}
+\begin{textblock}{13}(0.300,5.800)\includegraphics{Figure_2/densityplot-q_arm.pdf}\end{textblock}
+\begin{textblock}{13}(-0.05,6.000)\LARGE{(B)}\end{textblock}
+\begin{textblock}{13}(7.000,1.100)\includegraphics[width=2.30in,keepaspectratio]{Figure_2/densityplot-legend.pdf}\end{textblock}
+\end{document}

assets/paper/figures/Figure_4.tex

@@ -0,0 +1,20 @@
+\documentclass{article}
+\usepackage[paperheight=14.724in,paperwidth=12.108in,margin=0in]{geometry}
+\usepackage[sfdefault]{roboto}
+\usepackage{graphicx}
+\usepackage{tikz}
+\usepackage[absolute,overlay]{textpos}
+    \setlength{\TPHorizModule}{1in}
+    \setlength{\TPVertModule}{1in}
+\begin{document}
+\begin{textblock}{13}(-0.200,0.000)\includegraphics{Figure_4/chr2.pdf}\end{textblock}
+\begin{textblock}{13}(-0.200,1.727)\includegraphics{Figure_4/3ptel_1-500K_1_12_12.pdf}\end{textblock}
+\begin{textblock}{13}(-0.200,3.288)\includegraphics{Figure_4/4ptel_1-500K_1_12_12.pdf}\end{textblock}
+\begin{textblock}{13}(-0.200,7.862)\includegraphics{Figure_4/chr5.pdf}\end{textblock}
+\begin{textblock}{13}(-0.200,9.006)\includegraphics{Figure_4/chr9.pdf}\end{textblock}
+\begin{textblock}{13}(-0.200,10.067)\includegraphics{Figure_4/chr12.pdf}\end{textblock}
+\begin{textblock}{13}(-0.200,11.322)\includegraphics{Figure_4/17ptel_1_500K_1_12_12.pdf}\end{textblock}
+\begin{textblock}{13}(7.750,0)
+\includegraphics[width=4.000in,keepaspectratio]{Figure_4/legend.pdf}
+\end{textblock}
+\end{document}