Snakemake Report

Jobs 1

Input files

/net/eichler/vol26/projects/chm13_t2t/nobackups/Assembly_analysis/SEDEF/chm13_v1.1_plus38Y_masked.fasta.fai
results/snv_count_annotated_haplotype_coverage.bed.gz

Output files

results/uncallable_regions.bed.gz

Conda software stack

snakemake
snakefmt
black
snakedeploy
pandas
minimap2>=2.20
samtools>=1.11
bedtools
seqtk
vcftools
bcftools
htslib

Code

1 2	`bedtools complement -i <(bedtools merge -i {input.windows}) -g <(cat {input.fai} \| sort -k 1,1 -k2,2n ) \| gzip -c > {output}`

Jobs 1

Input files

temp/snvs_haplotype_coverage.bed.gz
temp/anno/Sat.bed.gz
temp/anno/SD.bed.gz
temp/anno/TRF.bed.gz
temp/anno/RM.bed.gz

Output files

results/snv_count_annotated_haplotype_coverage.bed.gz

Conda software stack

snakemake
snakefmt
black
snakedeploy
pandas
minimap2>=2.20
samtools>=1.11
bedtools
seqtk
vcftools
bcftools
htslib

Code

        HEADER=$(gunzip -c {input.windows} | head -n 1 || :)
        HEADER="${{HEADER}}	{params.annotation_names}"
        echo $HEADER

        bedtools annotate -i {input.windows}                 -files {input.annotation_files}             | bedtools sort -i -             | sed "1s/^/${{HEADER}}\n/"             | gzip -c > {output}

Jobs 1

Input files

temp/haplotype_coverage.bed.gz
temp/syntenic_and_callable/snv_HG002_1.bed.gz
temp/syntenic_and_callable/snv_HG002_2.bed.gz
temp/syntenic_and_callable/snv_HG00438_1.bed.gz
temp/syntenic_and_callable/snv_HG00438_2.bed.gz
temp/syntenic_and_callable/snv_HG005_1.bed.gz
temp/syntenic_and_callable/snv_HG005_2.bed.gz
temp/syntenic_and_callable/snv_HG00621_1.bed.gz
temp/syntenic_and_callable/snv_HG00621_2.bed.gz
temp/syntenic_and_callable/snv_HG00673_1.bed.gz
temp/syntenic_and_callable/snv_HG00673_2.bed.gz
temp/syntenic_and_callable/snv_HG00733_1.bed.gz
temp/syntenic_and_callable/snv_HG00733_2.bed.gz
temp/syntenic_and_callable/snv_HG00735_1.bed.gz
temp/syntenic_and_callable/snv_HG00735_2.bed.gz
temp/syntenic_and_callable/snv_HG00741_1.bed.gz
temp/syntenic_and_callable/snv_HG00741_2.bed.gz
temp/syntenic_and_callable/snv_HG01071_1.bed.gz
temp/syntenic_and_callable/snv_HG01071_2.bed.gz
temp/syntenic_and_callable/snv_HG01106_1.bed.gz
temp/syntenic_and_callable/snv_HG01106_2.bed.gz
temp/syntenic_and_callable/snv_HG01109_1.bed.gz
temp/syntenic_and_callable/snv_HG01109_2.bed.gz
temp/syntenic_and_callable/snv_HG01123_1.bed.gz
temp/syntenic_and_callable/snv_HG01123_2.bed.gz
temp/syntenic_and_callable/snv_HG01175_1.bed.gz
temp/syntenic_and_callable/snv_HG01175_2.bed.gz
temp/syntenic_and_callable/snv_HG01243_1.bed.gz
temp/syntenic_and_callable/snv_HG01243_2.bed.gz
temp/syntenic_and_callable/snv_HG01258_1.bed.gz
temp/syntenic_and_callable/snv_HG01258_2.bed.gz
temp/syntenic_and_callable/snv_HG01358_1.bed.gz
temp/syntenic_and_callable/snv_HG01358_2.bed.gz
temp/syntenic_and_callable/snv_HG01361_1.bed.gz
temp/syntenic_and_callable/snv_HG01361_2.bed.gz
temp/syntenic_and_callable/snv_HG01891_1.bed.gz
temp/syntenic_and_callable/snv_HG01891_2.bed.gz
temp/syntenic_and_callable/snv_HG01928_1.bed.gz
temp/syntenic_and_callable/snv_HG01928_2.bed.gz
temp/syntenic_and_callable/snv_HG01952_1.bed.gz
temp/syntenic_and_callable/snv_HG01952_2.bed.gz
temp/syntenic_and_callable/snv_HG01978_1.bed.gz
temp/syntenic_and_callable/snv_HG01978_2.bed.gz
temp/syntenic_and_callable/snv_HG02055_1.bed.gz
temp/syntenic_and_callable/snv_HG02055_2.bed.gz
temp/syntenic_and_callable/snv_HG02080_1.bed.gz
temp/syntenic_and_callable/snv_HG02080_2.bed.gz
temp/syntenic_and_callable/snv_HG02109_1.bed.gz
temp/syntenic_and_callable/snv_HG02109_2.bed.gz
temp/syntenic_and_callable/snv_HG02145_1.bed.gz
temp/syntenic_and_callable/snv_HG02145_2.bed.gz
temp/syntenic_and_callable/snv_HG02148_1.bed.gz
temp/syntenic_and_callable/snv_HG02148_2.bed.gz
temp/syntenic_and_callable/snv_HG02257_1.bed.gz
temp/syntenic_and_callable/snv_HG02257_2.bed.gz
temp/syntenic_and_callable/snv_HG02486_1.bed.gz
temp/syntenic_and_callable/snv_HG02486_2.bed.gz
temp/syntenic_and_callable/snv_HG02559_1.bed.gz
temp/syntenic_and_callable/snv_HG02559_2.bed.gz
temp/syntenic_and_callable/snv_HG02572_1.bed.gz
temp/syntenic_and_callable/snv_HG02572_2.bed.gz
temp/syntenic_and_callable/snv_HG02622_1.bed.gz
temp/syntenic_and_callable/snv_HG02622_2.bed.gz
temp/syntenic_and_callable/snv_HG02630_1.bed.gz
temp/syntenic_and_callable/snv_HG02630_2.bed.gz
temp/syntenic_and_callable/snv_HG02717_1.bed.gz
temp/syntenic_and_callable/snv_HG02717_2.bed.gz
temp/syntenic_and_callable/snv_HG02723_1.bed.gz
temp/syntenic_and_callable/snv_HG02723_2.bed.gz
temp/syntenic_and_callable/snv_HG02818_1.bed.gz
temp/syntenic_and_callable/snv_HG02818_2.bed.gz
temp/syntenic_and_callable/snv_HG02886_1.bed.gz
temp/syntenic_and_callable/snv_HG02886_2.bed.gz
temp/syntenic_and_callable/snv_HG03098_1.bed.gz
temp/syntenic_and_callable/snv_HG03098_2.bed.gz
temp/syntenic_and_callable/snv_HG03453_1.bed.gz
temp/syntenic_and_callable/snv_HG03453_2.bed.gz
temp/syntenic_and_callable/snv_HG03486_1.bed.gz
temp/syntenic_and_callable/snv_HG03486_2.bed.gz
temp/syntenic_and_callable/snv_HG03492_1.bed.gz
temp/syntenic_and_callable/snv_HG03492_2.bed.gz
temp/syntenic_and_callable/snv_HG03516_1.bed.gz
temp/syntenic_and_callable/snv_HG03516_2.bed.gz
temp/syntenic_and_callable/snv_HG03540_1.bed.gz
temp/syntenic_and_callable/snv_HG03540_2.bed.gz
temp/syntenic_and_callable/snv_HG03579_1.bed.gz
temp/syntenic_and_callable/snv_HG03579_2.bed.gz
temp/syntenic_and_callable/snv_NA18906_1.bed.gz
temp/syntenic_and_callable/snv_NA18906_2.bed.gz
temp/syntenic_and_callable/snv_NA19240_1.bed.gz
temp/syntenic_and_callable/snv_NA19240_2.bed.gz
temp/syntenic_and_callable/snv_NA20129_1.bed.gz
temp/syntenic_and_callable/snv_NA20129_2.bed.gz
temp/syntenic_and_callable/snv_NA21309_1.bed.gz
temp/syntenic_and_callable/snv_NA21309_2.bed.gz

Output files

temp/snvs_haplotype_coverage.bed.gz

Conda software stack

snakemake
snakefmt
black
snakedeploy
pandas
minimap2>=2.20
samtools>=1.11
bedtools
seqtk
vcftools
bcftools
htslib

Code

1 2	`bedtools coverage -a {input.windows} -b {input.snv} -header -counts -sorted \| sed "1s/$/\tnum_snv/" \| gzip -c > {output}`

Jobs 1

Input files

temp/windows.bed.gz
results/syntenic_and_callable/HG002_1.bed.gz
results/syntenic_and_callable/HG002_2.bed.gz
results/syntenic_and_callable/HG00438_1.bed.gz
results/syntenic_and_callable/HG00438_2.bed.gz
results/syntenic_and_callable/HG005_1.bed.gz
results/syntenic_and_callable/HG005_2.bed.gz
results/syntenic_and_callable/HG00621_1.bed.gz
results/syntenic_and_callable/HG00621_2.bed.gz
results/syntenic_and_callable/HG00673_1.bed.gz
results/syntenic_and_callable/HG00673_2.bed.gz
results/syntenic_and_callable/HG00733_1.bed.gz
results/syntenic_and_callable/HG00733_2.bed.gz
results/syntenic_and_callable/HG00735_1.bed.gz
results/syntenic_and_callable/HG00735_2.bed.gz
results/syntenic_and_callable/HG00741_1.bed.gz
results/syntenic_and_callable/HG00741_2.bed.gz
results/syntenic_and_callable/HG01071_1.bed.gz
results/syntenic_and_callable/HG01071_2.bed.gz
results/syntenic_and_callable/HG01106_1.bed.gz
results/syntenic_and_callable/HG01106_2.bed.gz
results/syntenic_and_callable/HG01109_1.bed.gz
results/syntenic_and_callable/HG01109_2.bed.gz
results/syntenic_and_callable/HG01123_1.bed.gz
results/syntenic_and_callable/HG01123_2.bed.gz
results/syntenic_and_callable/HG01175_1.bed.gz
results/syntenic_and_callable/HG01175_2.bed.gz
results/syntenic_and_callable/HG01243_1.bed.gz
results/syntenic_and_callable/HG01243_2.bed.gz
results/syntenic_and_callable/HG01258_1.bed.gz
results/syntenic_and_callable/HG01258_2.bed.gz
results/syntenic_and_callable/HG01358_1.bed.gz
results/syntenic_and_callable/HG01358_2.bed.gz
results/syntenic_and_callable/HG01361_1.bed.gz
results/syntenic_and_callable/HG01361_2.bed.gz
results/syntenic_and_callable/HG01891_1.bed.gz
results/syntenic_and_callable/HG01891_2.bed.gz
results/syntenic_and_callable/HG01928_1.bed.gz
results/syntenic_and_callable/HG01928_2.bed.gz
results/syntenic_and_callable/HG01952_1.bed.gz
results/syntenic_and_callable/HG01952_2.bed.gz
results/syntenic_and_callable/HG01978_1.bed.gz
results/syntenic_and_callable/HG01978_2.bed.gz
results/syntenic_and_callable/HG02055_1.bed.gz
results/syntenic_and_callable/HG02055_2.bed.gz
results/syntenic_and_callable/HG02080_1.bed.gz
results/syntenic_and_callable/HG02080_2.bed.gz
results/syntenic_and_callable/HG02109_1.bed.gz
results/syntenic_and_callable/HG02109_2.bed.gz
results/syntenic_and_callable/HG02145_1.bed.gz
results/syntenic_and_callable/HG02145_2.bed.gz
results/syntenic_and_callable/HG02148_1.bed.gz
results/syntenic_and_callable/HG02148_2.bed.gz
results/syntenic_and_callable/HG02257_1.bed.gz
results/syntenic_and_callable/HG02257_2.bed.gz
results/syntenic_and_callable/HG02486_1.bed.gz
results/syntenic_and_callable/HG02486_2.bed.gz
results/syntenic_and_callable/HG02559_1.bed.gz
results/syntenic_and_callable/HG02559_2.bed.gz
results/syntenic_and_callable/HG02572_1.bed.gz
results/syntenic_and_callable/HG02572_2.bed.gz
results/syntenic_and_callable/HG02622_1.bed.gz
results/syntenic_and_callable/HG02622_2.bed.gz
results/syntenic_and_callable/HG02630_1.bed.gz
results/syntenic_and_callable/HG02630_2.bed.gz
results/syntenic_and_callable/HG02717_1.bed.gz
results/syntenic_and_callable/HG02717_2.bed.gz
results/syntenic_and_callable/HG02723_1.bed.gz
results/syntenic_and_callable/HG02723_2.bed.gz
results/syntenic_and_callable/HG02818_1.bed.gz
results/syntenic_and_callable/HG02818_2.bed.gz
results/syntenic_and_callable/HG02886_1.bed.gz
results/syntenic_and_callable/HG02886_2.bed.gz
results/syntenic_and_callable/HG03098_1.bed.gz
results/syntenic_and_callable/HG03098_2.bed.gz
results/syntenic_and_callable/HG03453_1.bed.gz
results/syntenic_and_callable/HG03453_2.bed.gz
results/syntenic_and_callable/HG03486_1.bed.gz
results/syntenic_and_callable/HG03486_2.bed.gz
results/syntenic_and_callable/HG03492_1.bed.gz
results/syntenic_and_callable/HG03492_2.bed.gz
results/syntenic_and_callable/HG03516_1.bed.gz
results/syntenic_and_callable/HG03516_2.bed.gz
results/syntenic_and_callable/HG03540_1.bed.gz
results/syntenic_and_callable/HG03540_2.bed.gz
results/syntenic_and_callable/HG03579_1.bed.gz
results/syntenic_and_callable/HG03579_2.bed.gz
results/syntenic_and_callable/NA18906_1.bed.gz
results/syntenic_and_callable/NA18906_2.bed.gz
results/syntenic_and_callable/NA19240_1.bed.gz
results/syntenic_and_callable/NA19240_2.bed.gz
results/syntenic_and_callable/NA20129_1.bed.gz
results/syntenic_and_callable/NA20129_2.bed.gz
results/syntenic_and_callable/NA21309_1.bed.gz
results/syntenic_and_callable/NA21309_2.bed.gz

Output files

temp/haplotype_coverage.bed.gz

Conda software stack

snakemake
snakefmt
black
snakedeploy
pandas
minimap2>=2.20
samtools>=1.11
bedtools
seqtk
vcftools
bcftools
htslib

Code

        NUM_COL=$(gunzip -c {input.windows} | head -n 1 | awk '{{print NF}}' || :)
        NAME_COL=$((NUM_COL+1))
        BOOL_COL=$((NUM_COL+2))
        echo $NAME_COL $BOOL_COL

        bedtools intersect -f 0.95 -C             -sorted -sortout             -a {input.windows} -b {input.beds}             -names {params.names}         | awk -v bool=$BOOL_COL '$bool != 0'         | bedtools merge -i -             -d {params.overlap}             -c $NAME_COL,$BOOL_COL -o distinct,sum         | sed '1s/^/#chr\tstart\tend\thaps\thap_count\n/'         | gzip -c > {output}

Jobs 1

Input files

/net/eichler/vol26/projects/chm13_t2t/nobackups/Assembly_analysis/SEDEF/chm13_v1.1_plus38Y_masked.fasta.fai
/net/eichler/vol26/projects/chm13_t2t/nobackups/Assembly_analysis/SEDEF/chm13_v1.1_plus38Y.cen.bed
/net/eichler/vol26/projects/chm13_t2t/nobackups/Assembly_analysis/SEDEF/chm13_v1.1_plus38Y.SDs.bed
/net/eichler/vol26/projects/chm13_t2t/nobackups/Assembly_analysis/Masked/chm13_v1.1_plus38Y_trf.bed
/net/eichler/vol26/projects/chm13_t2t/nobackups/Assembly_analysis/Masked/chm13_v1.1_plus38Y_repeatmasker.out.bed

Output files

temp/windows.bed.gz

Conda software stack

snakemake
snakefmt
black
snakedeploy
pandas
minimap2>=2.20
samtools>=1.11
bedtools
seqtk
vcftools
bcftools
htslib

Code

1 2	`bedtools makewindows -g {input.fai} -w {params.window_size} -s {params.step_size} \| bedtools sort -i - \| gzip -c > {output}`

Jobs 94

Input files

at 0x7fc971881280>
at 0x7fc971881310>

Output files

results/syntenic_and_callable/{sm,HG002|HG00438|HG005|HG00621|HG00673|HG00733|HG00735|HG00741|HG01071|HG01106|HG01109|HG01123|HG01175|HG01243|HG01258|HG01358|HG01361|HG01891|HG01928|HG01952|HG01978|HG02055|HG02080|HG02109|HG02145|HG02148|HG02257|HG02486|HG02559|HG02572|HG02622|HG02630|HG02717|HG02723|HG02818|HG02886|HG03098|HG03453|HG03486|HG03492|HG03516|HG03540|HG03579|NA18906|NA19240|NA20129|NA21309|all}_{h}.bed.gz

Conda software stack

snakemake
snakefmt
black
snakedeploy
pandas
minimap2>=2.20
samtools>=1.11
bedtools
seqtk
vcftools
bcftools
htslib

Code

1
2

        bedtools intersect -a {input.callable}             -b <( bedtools sort -i {input.aln} | awk '$3 - $2 >= {params.min_syntenic_size}' )             | bedtools sort -i -             | bedtools merge -i -             | sed 's/$/\t{wildcards.sm}_{wildcards.h}/g'             | gzip -c > {output}

Jobs 94

Input files

temp/syntenic_and_callable/snv_explode_{sm}.bed.gz
results/syntenic_and_callable/{sm}_{h}.bed.gz

Output files

temp/syntenic_and_callable/snv_{sm,HG002|HG00438|HG005|HG00621|HG00673|HG00733|HG00735|HG00741|HG01071|HG01106|HG01109|HG01123|HG01175|HG01243|HG01258|HG01358|HG01361|HG01891|HG01928|HG01952|HG01978|HG02055|HG02080|HG02109|HG02145|HG02148|HG02257|HG02486|HG02559|HG02572|HG02622|HG02630|HG02717|HG02723|HG02818|HG02886|HG03098|HG03453|HG03486|HG03492|HG03516|HG03540|HG03579|NA18906|NA19240|NA20129|NA21309|all}_{h}.bed.gz

Conda software stack

snakemake
snakefmt
black
snakedeploy
pandas
minimap2>=2.20
samtools>=1.11
bedtools
seqtk
vcftools
bcftools
htslib

Code

1 2	`gunzip -c {input.snv} \| awk -F$'\t' -v pat="h{wildcards.h}\|HAP" '$15 ~ pat' \| bedtools intersect -u -a - -b <(gunzip -c {input.callable}) -header \| bedtools sort -header -i - \| gzip -c > {output}`

Jobs 47

Input files

at 0x7fc971881b80>

Output files

temp/syntenic_and_callable/snv_explode_{sm,HG002|HG00438|HG005|HG00621|HG00673|HG00733|HG00735|HG00741|HG01071|HG01106|HG01109|HG01123|HG01175|HG01243|HG01258|HG01358|HG01361|HG01891|HG01928|HG01952|HG01978|HG02055|HG02080|HG02109|HG02145|HG02148|HG02257|HG02486|HG02559|HG02572|HG02622|HG02630|HG02717|HG02723|HG02818|HG02886|HG03098|HG03453|HG03486|HG03492|HG03516|HG03540|HG03579|NA18906|NA19240|NA20129|NA21309|all}.bed.gz

Conda software stack

snakemake
snakefmt
black
snakedeploy
pandas
minimap2>=2.20
samtools>=1.11
bedtools
seqtk
vcftools
bcftools
htslib

Code

import pandas as pd

df = pd.read_csv(snakemake.input.bed, sep="\t", low_memory=False)
df = df.astype(str)
df = (
    df.apply(lambda col: col.str.split(";").explode())
    .reset_index()
    .reindex(df.columns, axis=1)
)
# fix the genotypes for expanded
df.loc[(df["HAP"] == "h1") & (df["GT"] == "1|1"), "GT"] = "1|."
df.loc[(df["HAP"] == "h2") & (df["GT"] == "1|1"), "GT"] = ".|1"

df["SAMPLE"] = snakemake.wildcards.sm
df.to_csv(snakemake.output.snv, sep="\t", compression="gzip", index=False)

Jobs 4

Input files

at 0x7fc971881c10>

Output files

temp/anno/{anno}.bed.gz

Conda software stack

snakemake
snakefmt
black
snakedeploy
pandas
minimap2>=2.20
samtools>=1.11
bedtools
seqtk
vcftools
bcftools
htslib

Code

1 2	`bedtools sort -i {input.annotation_file} \| bedtools merge -i - \| gzip -c > {output}`

Jobs 1

Input files

temp/syntenic_and_callable/snv_HG002_1.bed.gz
temp/syntenic_and_callable/snv_HG002_2.bed.gz
temp/syntenic_and_callable/snv_HG00438_1.bed.gz
temp/syntenic_and_callable/snv_HG00438_2.bed.gz
temp/syntenic_and_callable/snv_HG005_1.bed.gz
temp/syntenic_and_callable/snv_HG005_2.bed.gz
temp/syntenic_and_callable/snv_HG00621_1.bed.gz
temp/syntenic_and_callable/snv_HG00621_2.bed.gz
temp/syntenic_and_callable/snv_HG00673_1.bed.gz
temp/syntenic_and_callable/snv_HG00673_2.bed.gz
temp/syntenic_and_callable/snv_HG00733_1.bed.gz
temp/syntenic_and_callable/snv_HG00733_2.bed.gz
temp/syntenic_and_callable/snv_HG00735_1.bed.gz
temp/syntenic_and_callable/snv_HG00735_2.bed.gz
temp/syntenic_and_callable/snv_HG00741_1.bed.gz
temp/syntenic_and_callable/snv_HG00741_2.bed.gz
temp/syntenic_and_callable/snv_HG01071_1.bed.gz
temp/syntenic_and_callable/snv_HG01071_2.bed.gz
temp/syntenic_and_callable/snv_HG01106_1.bed.gz
temp/syntenic_and_callable/snv_HG01106_2.bed.gz
temp/syntenic_and_callable/snv_HG01109_1.bed.gz
temp/syntenic_and_callable/snv_HG01109_2.bed.gz
temp/syntenic_and_callable/snv_HG01123_1.bed.gz
temp/syntenic_and_callable/snv_HG01123_2.bed.gz
temp/syntenic_and_callable/snv_HG01175_1.bed.gz
temp/syntenic_and_callable/snv_HG01175_2.bed.gz
temp/syntenic_and_callable/snv_HG01243_1.bed.gz
temp/syntenic_and_callable/snv_HG01243_2.bed.gz
temp/syntenic_and_callable/snv_HG01258_1.bed.gz
temp/syntenic_and_callable/snv_HG01258_2.bed.gz
temp/syntenic_and_callable/snv_HG01358_1.bed.gz
temp/syntenic_and_callable/snv_HG01358_2.bed.gz
temp/syntenic_and_callable/snv_HG01361_1.bed.gz
temp/syntenic_and_callable/snv_HG01361_2.bed.gz
temp/syntenic_and_callable/snv_HG01891_1.bed.gz
temp/syntenic_and_callable/snv_HG01891_2.bed.gz
temp/syntenic_and_callable/snv_HG01928_1.bed.gz
temp/syntenic_and_callable/snv_HG01928_2.bed.gz
temp/syntenic_and_callable/snv_HG01952_1.bed.gz
temp/syntenic_and_callable/snv_HG01952_2.bed.gz
temp/syntenic_and_callable/snv_HG01978_1.bed.gz
temp/syntenic_and_callable/snv_HG01978_2.bed.gz
temp/syntenic_and_callable/snv_HG02055_1.bed.gz
temp/syntenic_and_callable/snv_HG02055_2.bed.gz
temp/syntenic_and_callable/snv_HG02080_1.bed.gz
temp/syntenic_and_callable/snv_HG02080_2.bed.gz
temp/syntenic_and_callable/snv_HG02109_1.bed.gz
temp/syntenic_and_callable/snv_HG02109_2.bed.gz
temp/syntenic_and_callable/snv_HG02145_1.bed.gz
temp/syntenic_and_callable/snv_HG02145_2.bed.gz
temp/syntenic_and_callable/snv_HG02148_1.bed.gz
temp/syntenic_and_callable/snv_HG02148_2.bed.gz
temp/syntenic_and_callable/snv_HG02257_1.bed.gz
temp/syntenic_and_callable/snv_HG02257_2.bed.gz
temp/syntenic_and_callable/snv_HG02486_1.bed.gz
temp/syntenic_and_callable/snv_HG02486_2.bed.gz
temp/syntenic_and_callable/snv_HG02559_1.bed.gz
temp/syntenic_and_callable/snv_HG02559_2.bed.gz
temp/syntenic_and_callable/snv_HG02572_1.bed.gz
temp/syntenic_and_callable/snv_HG02572_2.bed.gz
temp/syntenic_and_callable/snv_HG02622_1.bed.gz
temp/syntenic_and_callable/snv_HG02622_2.bed.gz
temp/syntenic_and_callable/snv_HG02630_1.bed.gz
temp/syntenic_and_callable/snv_HG02630_2.bed.gz
temp/syntenic_and_callable/snv_HG02717_1.bed.gz
temp/syntenic_and_callable/snv_HG02717_2.bed.gz
temp/syntenic_and_callable/snv_HG02723_1.bed.gz
temp/syntenic_and_callable/snv_HG02723_2.bed.gz
temp/syntenic_and_callable/snv_HG02818_1.bed.gz
temp/syntenic_and_callable/snv_HG02818_2.bed.gz
temp/syntenic_and_callable/snv_HG02886_1.bed.gz
temp/syntenic_and_callable/snv_HG02886_2.bed.gz
temp/syntenic_and_callable/snv_HG03098_1.bed.gz
temp/syntenic_and_callable/snv_HG03098_2.bed.gz
temp/syntenic_and_callable/snv_HG03453_1.bed.gz
temp/syntenic_and_callable/snv_HG03453_2.bed.gz
temp/syntenic_and_callable/snv_HG03486_1.bed.gz
temp/syntenic_and_callable/snv_HG03486_2.bed.gz
temp/syntenic_and_callable/snv_HG03492_1.bed.gz
temp/syntenic_and_callable/snv_HG03492_2.bed.gz
temp/syntenic_and_callable/snv_HG03516_1.bed.gz
temp/syntenic_and_callable/snv_HG03516_2.bed.gz
temp/syntenic_and_callable/snv_HG03540_1.bed.gz
temp/syntenic_and_callable/snv_HG03540_2.bed.gz
temp/syntenic_and_callable/snv_HG03579_1.bed.gz
temp/syntenic_and_callable/snv_HG03579_2.bed.gz
temp/syntenic_and_callable/snv_NA18906_1.bed.gz
temp/syntenic_and_callable/snv_NA18906_2.bed.gz
temp/syntenic_and_callable/snv_NA19240_1.bed.gz
temp/syntenic_and_callable/snv_NA19240_2.bed.gz
temp/syntenic_and_callable/snv_NA20129_1.bed.gz
temp/syntenic_and_callable/snv_NA20129_2.bed.gz
temp/syntenic_and_callable/snv_NA21309_1.bed.gz
temp/syntenic_and_callable/snv_NA21309_2.bed.gz

Output files

results/all_snv_exploded.bed.gz

Conda software stack

snakemake
snakefmt
black
snakedeploy
pandas
minimap2>=2.20
samtools>=1.11
bedtools
seqtk
vcftools
bcftools
htslib

Code

1 2	`gunzip -c {input.snv} \| bedtools sort -i - -header \| gzip -c > {output}`

Jobs 47

Input files

temp/syntenic_and_callable/snv_{sm}_1.bed.gz
temp/syntenic_and_callable/snv_{sm}_2.bed.gz

Output files

results/syntenic_and_callable/snv_{sm,HG002|HG00438|HG005|HG00621|HG00673|HG00733|HG00735|HG00741|HG01071|HG01106|HG01109|HG01123|HG01175|HG01243|HG01258|HG01358|HG01361|HG01891|HG01928|HG01952|HG01978|HG02055|HG02080|HG02109|HG02145|HG02148|HG02257|HG02486|HG02559|HG02572|HG02622|HG02630|HG02717|HG02723|HG02818|HG02886|HG03098|HG03453|HG03486|HG03492|HG03516|HG03540|HG03579|NA18906|NA19240|NA20129|NA21309|all}.bed.gz

Conda software stack

snakemake
snakefmt
black
snakedeploy
pandas
minimap2>=2.20
samtools>=1.11
bedtools
seqtk
vcftools
bcftools
htslib

Code

import pandas as pd

li = []
for filename in snakemake.input.beds:
    tdf = pd.read_csv(filename, index_col=None, sep="\t", low_memory=False)
    li.append(tdf)

df = pd.concat(li, axis=0, ignore_index=True)
df = df.astype(str)

id_cols = ["ID"]
if "SAMPLE" in df.columns:
    id_cols = ["ID", "SAMPLE"]

implode = (
    df.sort_values(by=["#CHROM", "POS", "ID", "HAP"])
    .groupby(id_cols)
    .agg(lambda x: ";".join(x.unique()))
    .reset_index()
)
implode.loc[(implode["GT"] == "1|.;.|1"), "GT"] = "1|1"

implode = implode.loc[:, df.columns]

implode.to_csv(snakemake.output.snv, sep="\t", compression="gzip", index=False)

Jobs 1

Input files

results/snv_count_annotated_haplotype_coverage.bed.gz

Output files

results/figures/cumulative_divergence.svg

Conda software stack

r::r>=4.0
r-r.utils
r-base
r-essentials
r-tidyverse
r-scales
r-cowplot
r-plotly
r-glue
r::r-rcolorbrewer
r-data.table
r-ggrepel
r-argparse
r-ggnewscale
r-svglite
bioconductor-genomicranges
bioconductor-karyoploter
r-ggforce

Code

source("workflow/scripts/setup.R")

infile <- "results/snv_count_annotated_haplotype_coverage.bed.gz"
infile <- snakemake@input[[1]]
outfile <- snakemake@output[[1]]
df <- fread(infile)

annotations <- colnames(df)[7:ncol(df)]
df$region <- "Other"
for (anno in annotations) {
    df[df$region == "Other" & df[[anno]] > 0.90]$region <- anno
}
df[SD < 0.2 & Sat < 0.2 && region == "Other"]$region <- "Unique"

chrX <- copy(df[df[["#chr"]] == "chrX"])
chrX$region <- "chrX"
df <- rbind(df, chrX)

pal <- COLORS[unique(df$region)]
df$region <- factor(df$region)

df <- df %>%
    mutate(snv_per_kbp = 1000 * num_snv / (hap_count * (end - start))) %>%
    filter(snv_per_kbp > 0) %>%
    data.table()



fakeadd <- 0.1
p <- ggplot() +
    stat_ecdf(
        data = df,
        aes(snv_per_kbp + fakeadd, color = region),
        size = 1.5, alpha = 0.75
    ) +
    scale_x_log10(
        limits = c(fakeadd, NA),
        breaks = c(fakeadd, 0.1, 1, 10),
        labels = c("0", "0.10", "1.00", "10.00")
    ) +
    annotation_logticks(sides = "b") +
    scale_fill_manual(values = pal) +
    scale_color_manual(values = pal) +
    xlab(glue("% divergence of 10 kbp windows (1 kbp slide) aligned to T2T CHM13")) +
    ylab("Cumulative fraction of 5kbp windows") +
    theme_cowplot() +
    theme(legend.position = "top", legend.title = element_blank()) +
    guides(fill = guide_legend(ncol = length(pal) / 2))
p
ggsave(outfile, width = 8, height = 5, plot = p)

Jobs 1

Input files

results/syntenic_and_callable/HG002_1.bed.gz
results/syntenic_and_callable/HG002_2.bed.gz
results/syntenic_and_callable/HG00438_1.bed.gz
results/syntenic_and_callable/HG00438_2.bed.gz
results/syntenic_and_callable/HG005_1.bed.gz
results/syntenic_and_callable/HG005_2.bed.gz
results/syntenic_and_callable/HG00621_1.bed.gz
results/syntenic_and_callable/HG00621_2.bed.gz
results/syntenic_and_callable/HG00673_1.bed.gz
results/syntenic_and_callable/HG00673_2.bed.gz
results/syntenic_and_callable/HG00733_1.bed.gz
results/syntenic_and_callable/HG00733_2.bed.gz
results/syntenic_and_callable/HG00735_1.bed.gz
results/syntenic_and_callable/HG00735_2.bed.gz
results/syntenic_and_callable/HG00741_1.bed.gz
results/syntenic_and_callable/HG00741_2.bed.gz
results/syntenic_and_callable/HG01071_1.bed.gz
results/syntenic_and_callable/HG01071_2.bed.gz
results/syntenic_and_callable/HG01106_1.bed.gz
results/syntenic_and_callable/HG01106_2.bed.gz
results/syntenic_and_callable/HG01109_1.bed.gz
results/syntenic_and_callable/HG01109_2.bed.gz
results/syntenic_and_callable/HG01123_1.bed.gz
results/syntenic_and_callable/HG01123_2.bed.gz
results/syntenic_and_callable/HG01175_1.bed.gz
results/syntenic_and_callable/HG01175_2.bed.gz
results/syntenic_and_callable/HG01243_1.bed.gz
results/syntenic_and_callable/HG01243_2.bed.gz
results/syntenic_and_callable/HG01258_1.bed.gz
results/syntenic_and_callable/HG01258_2.bed.gz
results/syntenic_and_callable/HG01358_1.bed.gz
results/syntenic_and_callable/HG01358_2.bed.gz
results/syntenic_and_callable/HG01361_1.bed.gz
results/syntenic_and_callable/HG01361_2.bed.gz
results/syntenic_and_callable/HG01891_1.bed.gz
results/syntenic_and_callable/HG01891_2.bed.gz
results/syntenic_and_callable/HG01928_1.bed.gz
results/syntenic_and_callable/HG01928_2.bed.gz
results/syntenic_and_callable/HG01952_1.bed.gz
results/syntenic_and_callable/HG01952_2.bed.gz
results/syntenic_and_callable/HG01978_1.bed.gz
results/syntenic_and_callable/HG01978_2.bed.gz
results/syntenic_and_callable/HG02055_1.bed.gz
results/syntenic_and_callable/HG02055_2.bed.gz
results/syntenic_and_callable/HG02080_1.bed.gz
results/syntenic_and_callable/HG02080_2.bed.gz
results/syntenic_and_callable/HG02109_1.bed.gz
results/syntenic_and_callable/HG02109_2.bed.gz
results/syntenic_and_callable/HG02145_1.bed.gz
results/syntenic_and_callable/HG02145_2.bed.gz
results/syntenic_and_callable/HG02148_1.bed.gz
results/syntenic_and_callable/HG02148_2.bed.gz
results/syntenic_and_callable/HG02257_1.bed.gz
results/syntenic_and_callable/HG02257_2.bed.gz
results/syntenic_and_callable/HG02486_1.bed.gz
results/syntenic_and_callable/HG02486_2.bed.gz
results/syntenic_and_callable/HG02559_1.bed.gz
results/syntenic_and_callable/HG02559_2.bed.gz
results/syntenic_and_callable/HG02572_1.bed.gz
results/syntenic_and_callable/HG02572_2.bed.gz
results/syntenic_and_callable/HG02622_1.bed.gz
results/syntenic_and_callable/HG02622_2.bed.gz
results/syntenic_and_callable/HG02630_1.bed.gz
results/syntenic_and_callable/HG02630_2.bed.gz
results/syntenic_and_callable/HG02717_1.bed.gz
results/syntenic_and_callable/HG02717_2.bed.gz
results/syntenic_and_callable/HG02723_1.bed.gz
results/syntenic_and_callable/HG02723_2.bed.gz
results/syntenic_and_callable/HG02818_1.bed.gz
results/syntenic_and_callable/HG02818_2.bed.gz
results/syntenic_and_callable/HG02886_1.bed.gz
results/syntenic_and_callable/HG02886_2.bed.gz
results/syntenic_and_callable/HG03098_1.bed.gz
results/syntenic_and_callable/HG03098_2.bed.gz
results/syntenic_and_callable/HG03453_1.bed.gz
results/syntenic_and_callable/HG03453_2.bed.gz
results/syntenic_and_callable/HG03486_1.bed.gz
results/syntenic_and_callable/HG03486_2.bed.gz
results/syntenic_and_callable/HG03492_1.bed.gz
results/syntenic_and_callable/HG03492_2.bed.gz
results/syntenic_and_callable/HG03516_1.bed.gz
results/syntenic_and_callable/HG03516_2.bed.gz
results/syntenic_and_callable/HG03540_1.bed.gz
results/syntenic_and_callable/HG03540_2.bed.gz
results/syntenic_and_callable/HG03579_1.bed.gz
results/syntenic_and_callable/HG03579_2.bed.gz
results/syntenic_and_callable/NA18906_1.bed.gz
results/syntenic_and_callable/NA18906_2.bed.gz
results/syntenic_and_callable/NA19240_1.bed.gz
results/syntenic_and_callable/NA19240_2.bed.gz
results/syntenic_and_callable/NA20129_1.bed.gz
results/syntenic_and_callable/NA20129_2.bed.gz
results/syntenic_and_callable/NA21309_1.bed.gz
results/syntenic_and_callable/NA21309_2.bed.gz
results/snv_count_annotated_haplotype_coverage.bed.gz

Output files

results/figures/callable_space.svg

Conda software stack

r::r>=4.0
r-r.utils
r-base
r-essentials
r-tidyverse
r-scales
r-cowplot
r-plotly
r-glue
r::r-rcolorbrewer
r-data.table
r-ggrepel
r-argparse
r-ggnewscale
r-svglite
bioconductor-genomicranges
bioconductor-karyoploter
r-ggforce

Code

source("workflow/scripts/setup.R")


all.files <- Sys.glob("results/syntenic_and_callable/H*bed.gz")
all.files <- snakemake@input[["beds"]]
mylist <- lapply(all.files, fread)
df <- rbindlist(mylist, idcol = TRUE)
colnames(df) <- c("ID", "chr", "start", "end", "haplotype")

p <- df %>%
    group_by(haplotype) %>%
    summarize(Mbp = sum(end - start) / 1e6) %>%
    ggplot(aes(x = Mbp, y = haplotype, label = comma(Mbp))) +
    geom_bar(stat = "identity") +
    geom_vline(xintercept = 3100, color = RED, linetype = "dashed", size = 1) +
    geom_text(hjust = 0) +
    theme_cowplot() +
    ggtitle("Callable space for each haplotype assebmly",
        subtitle = "(Only regions with at least 1 Mbp of synteny considered)"
    )
p




# called_regions <- fread("results/snv_count_annotated_haplotype_coverage.bed.gz")
called_regions <- fread(snakemake@input[["windows"]])
p2 <- called_regions %>%
    mutate(region = case_when(
        SD > 0.95 & Sat < 0.7 ~ "SD",
        Sat > 0.7 ~ "Sat",
        TRUE ~ "Unique"
    )) %>%
    separate_rows(haps, sep = ",") %>%
    group_by(haps, region) %>%
    do(get_num_bp(.)) %>%
    ggplot(aes(
        x = V1 / 1e6,
        y = haps,
        label = comma(round(V1 / 1e6)), fill = region
    )) +
    geom_bar(stat = "identity") +
    geom_vline(xintercept = 3100, color = RED, linetype = "dashed", size = 1) +
    geom_text(hjust = 0) +
    facet_col(. ~ region, scales = "free") +
    ylab("") +
    xlab("Mbp") +
    theme_cowplot() +
    theme(legend.position = "top")

fig <- cowplot::plot_grid(p, p2)
fig
ggsave(snakemake@output[[1]], plot = fig, width = 16, height = 9)

Workflow

Figures

Statistics