Snakemake Report

Jobs 1

Input files

/net/eichler/vol26/projects/chm13_t2t/nobackups/Assembly_analysis/SEDEF/chm13_v1.1_plus38Y_masked.fasta.fai
results/annotated_windows.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/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/annotated_windows.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/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 0x7f396afbd5e0>
at 0x7f396afbd670>

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}_{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 4

Input files

at 0x7f396afc55e0>

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 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}.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 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}_{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 0x7f396afbdee0>

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}.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 1

Input files

results/long_and_filtered_snv_over_windows.bed.gz

Output files

results/figures/cumulative_divergence.svg
results/figures/cumulative_divergence_per_hap.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 <- snakemake@input[[1]]
outfile <- snakemake@output[[1]]
outfile2 <- snakemake@output[[2]]

df <- read_in_snv_windows(infile)


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

pal <- COLORS[unique(df$region)]
df$region <- factor(df$region, levels = names(COLORS))

#
# make plot
#
fakeadd <- 0.005
plot.df <- df %>%
    filter(region %in% c("Unique", "SD", "chrX")) %>%
    data.table()
plot.df[per_div == 0]$per_div <- fakeadd

p <- ggplot() +
    stat_ecdf(
        data = plot.df,
        aes(per_div, color = region),
        size = 1.5, alpha = 0.75
    ) +
    scale_x_log10(
        limits = c(fakeadd, 20),
        breaks = c(fakeadd, 0.01, 0.1, 1, 10),
        labels = c("0.00", "0.01", "0.10", "1.00", "10.0")
    ) +
    annotation_logticks(sides = "b") +
    scale_fill_manual(values = pal) +
    scale_color_manual(values = pal) +
    xlab(glue("% divergence of 10 kbp windows (1 kbp slide)")) +
    ylab("Cumulative fraction of windows") +
    ggtitle("Divergence of 10 kbp windows aligned to T2T-CHM13 v1.1",
        subtitle = "(Minumum 1 Mbp alignment, SD windows are at least 95% SD)"
    ) +
    theme_cowplot() +
    theme(legend.position = "top", legend.title = element_blank()) +
    guides(fill = guide_legend(ncol = length(pal) / 2))
ggsave(outfile, width = 12, height = 8, plot = p)


#
# all haplotypes
#
plot.df2 <- plot.df %>%
    filter(region %in% c("SD", "Unique")) %>%
    data.table()

p2 <- p +
    stat_ecdf(
        data = plot.df2,
        aes(per_div, group = paste0(hap, region), color = region),
        alpha = 0.5,
        size = 0.1,
        linetype = "dashed"
    )
ggsave(outfile2, width = 12, height = 8, plot = p2)

Jobs 1

Input files

temp/add_snv/HG002_1_snvs_haplotype_coverage.bed
temp/add_snv/HG002_2_snvs_haplotype_coverage.bed
temp/add_snv/HG00438_1_snvs_haplotype_coverage.bed
temp/add_snv/HG00438_2_snvs_haplotype_coverage.bed
temp/add_snv/HG005_1_snvs_haplotype_coverage.bed
temp/add_snv/HG005_2_snvs_haplotype_coverage.bed
temp/add_snv/HG00621_1_snvs_haplotype_coverage.bed
temp/add_snv/HG00621_2_snvs_haplotype_coverage.bed
temp/add_snv/HG00673_1_snvs_haplotype_coverage.bed
temp/add_snv/HG00673_2_snvs_haplotype_coverage.bed
temp/add_snv/HG00733_1_snvs_haplotype_coverage.bed
temp/add_snv/HG00733_2_snvs_haplotype_coverage.bed
temp/add_snv/HG00735_1_snvs_haplotype_coverage.bed
temp/add_snv/HG00735_2_snvs_haplotype_coverage.bed
temp/add_snv/HG00741_1_snvs_haplotype_coverage.bed
temp/add_snv/HG00741_2_snvs_haplotype_coverage.bed
temp/add_snv/HG01071_1_snvs_haplotype_coverage.bed
temp/add_snv/HG01071_2_snvs_haplotype_coverage.bed
temp/add_snv/HG01106_1_snvs_haplotype_coverage.bed
temp/add_snv/HG01106_2_snvs_haplotype_coverage.bed
temp/add_snv/HG01109_1_snvs_haplotype_coverage.bed
temp/add_snv/HG01109_2_snvs_haplotype_coverage.bed
temp/add_snv/HG01123_1_snvs_haplotype_coverage.bed
temp/add_snv/HG01123_2_snvs_haplotype_coverage.bed
temp/add_snv/HG01175_1_snvs_haplotype_coverage.bed
temp/add_snv/HG01175_2_snvs_haplotype_coverage.bed
temp/add_snv/HG01243_1_snvs_haplotype_coverage.bed
temp/add_snv/HG01243_2_snvs_haplotype_coverage.bed
temp/add_snv/HG01258_1_snvs_haplotype_coverage.bed
temp/add_snv/HG01258_2_snvs_haplotype_coverage.bed
temp/add_snv/HG01358_1_snvs_haplotype_coverage.bed
temp/add_snv/HG01358_2_snvs_haplotype_coverage.bed
temp/add_snv/HG01361_1_snvs_haplotype_coverage.bed
temp/add_snv/HG01361_2_snvs_haplotype_coverage.bed
temp/add_snv/HG01891_1_snvs_haplotype_coverage.bed
temp/add_snv/HG01891_2_snvs_haplotype_coverage.bed
temp/add_snv/HG01928_1_snvs_haplotype_coverage.bed
temp/add_snv/HG01928_2_snvs_haplotype_coverage.bed
temp/add_snv/HG01952_1_snvs_haplotype_coverage.bed
temp/add_snv/HG01952_2_snvs_haplotype_coverage.bed
temp/add_snv/HG01978_1_snvs_haplotype_coverage.bed
temp/add_snv/HG01978_2_snvs_haplotype_coverage.bed
temp/add_snv/HG02055_1_snvs_haplotype_coverage.bed
temp/add_snv/HG02055_2_snvs_haplotype_coverage.bed
temp/add_snv/HG02080_1_snvs_haplotype_coverage.bed
temp/add_snv/HG02080_2_snvs_haplotype_coverage.bed
temp/add_snv/HG02109_1_snvs_haplotype_coverage.bed
temp/add_snv/HG02109_2_snvs_haplotype_coverage.bed
temp/add_snv/HG02145_1_snvs_haplotype_coverage.bed
temp/add_snv/HG02145_2_snvs_haplotype_coverage.bed
temp/add_snv/HG02148_1_snvs_haplotype_coverage.bed
temp/add_snv/HG02148_2_snvs_haplotype_coverage.bed
temp/add_snv/HG02257_1_snvs_haplotype_coverage.bed
temp/add_snv/HG02257_2_snvs_haplotype_coverage.bed
temp/add_snv/HG02486_1_snvs_haplotype_coverage.bed
temp/add_snv/HG02486_2_snvs_haplotype_coverage.bed
temp/add_snv/HG02559_1_snvs_haplotype_coverage.bed
temp/add_snv/HG02559_2_snvs_haplotype_coverage.bed
temp/add_snv/HG02572_1_snvs_haplotype_coverage.bed
temp/add_snv/HG02572_2_snvs_haplotype_coverage.bed
temp/add_snv/HG02622_1_snvs_haplotype_coverage.bed
temp/add_snv/HG02622_2_snvs_haplotype_coverage.bed
temp/add_snv/HG02630_1_snvs_haplotype_coverage.bed
temp/add_snv/HG02630_2_snvs_haplotype_coverage.bed
temp/add_snv/HG02717_1_snvs_haplotype_coverage.bed
temp/add_snv/HG02717_2_snvs_haplotype_coverage.bed
temp/add_snv/HG02723_1_snvs_haplotype_coverage.bed
temp/add_snv/HG02723_2_snvs_haplotype_coverage.bed
temp/add_snv/HG02818_1_snvs_haplotype_coverage.bed
temp/add_snv/HG02818_2_snvs_haplotype_coverage.bed
temp/add_snv/HG02886_1_snvs_haplotype_coverage.bed
temp/add_snv/HG02886_2_snvs_haplotype_coverage.bed
temp/add_snv/HG03098_1_snvs_haplotype_coverage.bed
temp/add_snv/HG03098_2_snvs_haplotype_coverage.bed
temp/add_snv/HG03453_1_snvs_haplotype_coverage.bed
temp/add_snv/HG03453_2_snvs_haplotype_coverage.bed
temp/add_snv/HG03486_1_snvs_haplotype_coverage.bed
temp/add_snv/HG03486_2_snvs_haplotype_coverage.bed
temp/add_snv/HG03492_1_snvs_haplotype_coverage.bed
temp/add_snv/HG03492_2_snvs_haplotype_coverage.bed
temp/add_snv/HG03516_1_snvs_haplotype_coverage.bed
temp/add_snv/HG03516_2_snvs_haplotype_coverage.bed
temp/add_snv/HG03540_1_snvs_haplotype_coverage.bed
temp/add_snv/HG03540_2_snvs_haplotype_coverage.bed
temp/add_snv/HG03579_1_snvs_haplotype_coverage.bed
temp/add_snv/HG03579_2_snvs_haplotype_coverage.bed
temp/add_snv/NA18906_1_snvs_haplotype_coverage.bed
temp/add_snv/NA18906_2_snvs_haplotype_coverage.bed
temp/add_snv/NA19240_1_snvs_haplotype_coverage.bed
temp/add_snv/NA19240_2_snvs_haplotype_coverage.bed
temp/add_snv/NA20129_1_snvs_haplotype_coverage.bed
temp/add_snv/NA20129_2_snvs_haplotype_coverage.bed
temp/add_snv/NA21309_1_snvs_haplotype_coverage.bed
temp/add_snv/NA21309_2_snvs_haplotype_coverage.bed

Output files

results/long_and_filtered_snv_over_windows.bed.gz

Conda software stack

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

Code

        HEADER=$(head -n 1 {input.snv[1]} || :)
        echo $HEADER

        sort -m -k 1,1 -k2,2n {input.snv}             | grep -v "^#"             | sed "1s/^/${{HEADER}}\n/"             | gzip -c             > {output}

Jobs 1

Input files

results/long_and_filtered_snv_over_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

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)"
    )

# read in long windows
called_regions <- read_in_snv_windows(snakemake@input[["windows"]])
called_regions$region <- factor(called_regions$region, levels = names(COLORS))
p2 <- called_regions %>%
    group_by(hap, region) %>%
    do(get_num_bp(.)) %>%
    ggplot(aes(
        x = V1 / 1e6,
        y = hap,
        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) +
    scale_fill_manual(values = COLORS) +
    facet_wrap(. ~ region, scales = "free") +
    ylab("") +
    xlab("Mbp") +
    theme_cowplot() +
    theme(legend.position = "top")

height <- length(unique(called_regions$hap)) / 2
# fig <- cowplot::plot_grid(p, p2)
fig <- p2
print(height)
ggsave(snakemake@output[[1]],
    plot = fig,
    width = 16, height = height,
    units = "in",
    limitsize = FALSE
)

Jobs 1

Input files

results/all_snv_exploded.bed.gz
results/long_and_filtered_snv_over_windows.bed.gz

Output files

results/tables/snv_per_kbp.tbl
results/tables

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
r-knitr

Code

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

infile <- "results/all_snv_exploded.bed.gz"
infile2 <- "results/long_and_filtered_snv_over_windows.bed.gz"

infile <- snakemake@input[[1]]
infile2 <- snakemake@input[[2]]
outfile <- snakemake@output[[1]]
outfile2 <- snakemake@output[[2]]

df <- fread(infile)
long_windows <- read_in_snv_windows(infile2)

dply_reduce <- function(df, groupvars) {
    tmp <- as.data.frame(reduce(toGRanges(as.data.frame(df))))
    tmp$region_size <- sum(tmp$end - tmp$start)
    tmp
}


hap_annotated_windows <- long_windows %>%
    filter(region %in% c("SD", "Unique")) %>%
    group_by(region, hap) %>%
    group_modify(dply_reduce) %>%
    select(-region, region) %>%
    select(-hap, hap) %>%
    data.table()


hap_region_overlaps <- function(df, groupvars) {
    gr <- toGRanges(as.data.frame(df))
    df.hap <- hap_annotated_windows[hap == groupvars$hap]
    gr2 <- toGRanges(df.hap)
    o <- findOverlaps(gr, gr2)
    cbind(
        df[queryHits(o), ],
        df.hap[subjectHits(o), c("region", "region_size")]
    )
}

df_snv_annotated <- df %>%
    mutate(hap = paste0(SAMPLE, gsub("h", "_", HAP))) %>%
    group_by(hap) %>%
    group_modify(hap_region_overlaps) %>%
    select(-hap, hap) %>%
    data.table()

out_df <- df_snv_annotated %>%
    group_by(hap, region) %>%
    summarize(
        "# SNVs" = n(),
        Mbp = unique(region_size) / 1e6,
        "# SNVs per kbp" = 1e3 * n() / unique(region_size),
    ) %>%
    pivot_wider(
        names_from = region,
        values_from = c("# SNVs", "Mbp", "# SNVs per kbp"),
        names_sort = TRUE,
        names_glue = "{region} {.value}",
    ) %>%
    data.table()

fwrite(out_df, outfile, sep = "\t", quote = FALSE, row.names = FALSE)

fileConn <- file(paste0(outfile2, "/index.html"))
x <- kable(out_df,
    format = "html",
    caption = "Average divergence statistics"
)
write(x, fileConn)
close(fileConn)

Jobs 1

Input files

temp/snv/anno_HG002_1.bed
temp/snv/anno_HG002_2.bed
temp/snv/anno_HG00438_1.bed
temp/snv/anno_HG00438_2.bed
temp/snv/anno_HG005_1.bed
temp/snv/anno_HG005_2.bed
temp/snv/anno_HG00621_1.bed
temp/snv/anno_HG00621_2.bed
temp/snv/anno_HG00673_1.bed
temp/snv/anno_HG00673_2.bed
temp/snv/anno_HG00733_1.bed
temp/snv/anno_HG00733_2.bed
temp/snv/anno_HG00735_1.bed
temp/snv/anno_HG00735_2.bed
temp/snv/anno_HG00741_1.bed
temp/snv/anno_HG00741_2.bed
temp/snv/anno_HG01071_1.bed
temp/snv/anno_HG01071_2.bed
temp/snv/anno_HG01106_1.bed
temp/snv/anno_HG01106_2.bed
temp/snv/anno_HG01109_1.bed
temp/snv/anno_HG01109_2.bed
temp/snv/anno_HG01123_1.bed
temp/snv/anno_HG01123_2.bed
temp/snv/anno_HG01175_1.bed
temp/snv/anno_HG01175_2.bed
temp/snv/anno_HG01243_1.bed
temp/snv/anno_HG01243_2.bed
temp/snv/anno_HG01258_1.bed
temp/snv/anno_HG01258_2.bed
temp/snv/anno_HG01358_1.bed
temp/snv/anno_HG01358_2.bed
temp/snv/anno_HG01361_1.bed
temp/snv/anno_HG01361_2.bed
temp/snv/anno_HG01891_1.bed
temp/snv/anno_HG01891_2.bed
temp/snv/anno_HG01928_1.bed
temp/snv/anno_HG01928_2.bed
temp/snv/anno_HG01952_1.bed
temp/snv/anno_HG01952_2.bed
temp/snv/anno_HG01978_1.bed
temp/snv/anno_HG01978_2.bed
temp/snv/anno_HG02055_1.bed
temp/snv/anno_HG02055_2.bed
temp/snv/anno_HG02080_1.bed
temp/snv/anno_HG02080_2.bed
temp/snv/anno_HG02109_1.bed
temp/snv/anno_HG02109_2.bed
temp/snv/anno_HG02145_1.bed
temp/snv/anno_HG02145_2.bed
temp/snv/anno_HG02148_1.bed
temp/snv/anno_HG02148_2.bed
temp/snv/anno_HG02257_1.bed
temp/snv/anno_HG02257_2.bed
temp/snv/anno_HG02486_1.bed
temp/snv/anno_HG02486_2.bed
temp/snv/anno_HG02559_1.bed
temp/snv/anno_HG02559_2.bed
temp/snv/anno_HG02572_1.bed
temp/snv/anno_HG02572_2.bed
temp/snv/anno_HG02622_1.bed
temp/snv/anno_HG02622_2.bed
temp/snv/anno_HG02630_1.bed
temp/snv/anno_HG02630_2.bed
temp/snv/anno_HG02717_1.bed
temp/snv/anno_HG02717_2.bed
temp/snv/anno_HG02723_1.bed
temp/snv/anno_HG02723_2.bed
temp/snv/anno_HG02818_1.bed
temp/snv/anno_HG02818_2.bed
temp/snv/anno_HG02886_1.bed
temp/snv/anno_HG02886_2.bed
temp/snv/anno_HG03098_1.bed
temp/snv/anno_HG03098_2.bed
temp/snv/anno_HG03453_1.bed
temp/snv/anno_HG03453_2.bed
temp/snv/anno_HG03486_1.bed
temp/snv/anno_HG03486_2.bed
temp/snv/anno_HG03492_1.bed
temp/snv/anno_HG03492_2.bed
temp/snv/anno_HG03516_1.bed
temp/snv/anno_HG03516_2.bed
temp/snv/anno_HG03540_1.bed
temp/snv/anno_HG03540_2.bed
temp/snv/anno_HG03579_1.bed
temp/snv/anno_HG03579_2.bed
temp/snv/anno_NA18906_1.bed
temp/snv/anno_NA18906_2.bed
temp/snv/anno_NA19240_1.bed
temp/snv/anno_NA19240_2.bed
temp/snv/anno_NA20129_1.bed
temp/snv/anno_NA20129_2.bed
temp/snv/anno_NA21309_1.bed
temp/snv/anno_NA21309_2.bed

Output files

results/all_snv_exploded.bed.gz

Conda software stack

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

Code

        HEADER=$(head -n 1 {input.snv[1]} || :)
        echo $HEADER

        sort -m -k 1,1 -k2,2n {input.snv}             | grep -v "^#"             | sed "1s/^/${{HEADER}}\n/"             | gzip -c             > {output}

Jobs 94

Input files

temp/syntenic_and_callable/snv_{sm}_{h}.bed.gz
temp/anno/Sat.bed.gz
temp/anno/SD.bed.gz
temp/anno/TRF.bed.gz
temp/anno/RM.bed.gz

Output files

temp/snv/anno_{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}_{h}.bed

Conda software stack

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

Code

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

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

Workflow

Figures

Tables

Statistics