Snakemake Report

Rule properties
Input files
Jobs	1
results/annotation/annotation_sizes.tbl
Output files
results/annotation/annotation_sizes_wide.tbl

Jobs 1

Input files

Output files

results/annotation/annotation_sizes.tbl

Conda software stack

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

Code

        printf "hap	anno	size
" > {output}
        cat {input.combos} | sort >> {output}

Jobs 846

Input files

results/annotation/combos/{sm}_{h}_intersection_{anno1}_and_{anno2}.bed.gz

Output files

temp/regions/{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}_intersection_{anno1}_and_{anno2}.tbl

Conda software stack

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

Code

        printf "{wildcards.sm}_{wildcards.h}	" > {output}
        if [ {wildcards.anno1} == {wildcards.anno2} ]; then
            printf "{wildcards.anno1}_size	" >> {output}
        else
            printf "{wildcards.anno1}_{wildcards.anno2}_size	" >> {output}
        fi 
        awk '{{sum += $3 - $2}}END{{print sum}}' <(gunzip -c {input}) >> {output}

Jobs 846

Input files

results/syntenic_and_callable/{sm}_{h}.bed.gz
temp/anno/{anno1}.bed.gz
temp/anno/{anno2}.bed.gz

Output files

results/annotation/combos/{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}_intersection_{anno1}_and_{anno2}.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

1 2	`bedtools intersect -sorted -a {input.bed} -b {input.anno1} \| bedtools intersect -sorted -a - -b {input.anno2} \| bedtools sort -i - \| bedtools merge -i - \| gzip -c > {output}`

Jobs 94

Input files

at 0x7f1f7c4ddee0>
at 0x7f1f7c4dde50>

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>=2.30
seqtk
vcftools
bcftools
htslib
csvtk

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 5

Input files

at 0x7f1f7c4718b0>

Output files

temp/anno/{anno}.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

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

Jobs 1

Input files

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

Output files

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

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

Rule properties
Input files
Jobs	1
temp/annotated_windows.bed.gz
Output files
results/windows.bed.gz

Jobs 1

Input files

temp/haplotype_coverage.bed.gz
temp/anno/Sat.bed.gz
temp/anno/SD.bed.gz
temp/anno/Unique.bed.gz
temp/anno/TRF.bed.gz
temp/anno/RM.bed.gz

Output files

temp/annotated_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=$(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>=2.30
seqtk
vcftools
bcftools
htslib
csvtk

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

Output files

temp/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

1 2	`bedtools makewindows -g {input.fai} -w {params.window_size} -s {params.step_size} \| bedtools sort -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>=2.30
seqtk
vcftools
bcftools
htslib
csvtk

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>=2.30
seqtk
vcftools
bcftools
htslib
csvtk

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 0x7f1f7c488e50>

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>=2.30
seqtk
vcftools
bcftools
htslib
csvtk

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/all_snv_exploded.bed.gz
results/annotation/annotation_sizes.tbl

Output files

results/tables/snv_per_kbp.tbl
results/tables/snv_per_kbp_wide.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")
threads <- 8
threads <- snakemake@threads
outfile <- snakemake@output[[1]]
outfile2 <- snakemake@output[[2]]
outfile3 <- snakemake@output[[3]]

snv_f <- "results/all_snv_exploded.bed.gz"
snv_f <- snakemake@input[[1]]
df <- fread(snv_f, showProgress = TRUE, nThread = threads) %>%
    mutate(hap = paste0(SAMPLE, gsub("h", "_", HAP))) %>%
    data.table()

if (F) {
    #
    # add paired annotation columns
    #
    anno_cols <- sort(names(df)[grepl("anno_", names(df))])
    for (anno in anno_cols) {
        type <- gsub("anno_", "", anno)
        if (type != "SD") {
            new_col <- paste0("anno_SD+", type)
            print(new_col)
            df[[new_col]] <- 0
            df[(anno_SD == 1 & df[[anno]] == 1), new_col] <- 1
        }
    }
    # remove columns that wont be annotated
    df <- df[rowSums(df[, ..anno_cols]) > 0]
    # get new longer pairs
    paired_anno_cols <- sort(names(df)[grepl("anno_", names(df))])
}
#
# make regions definitions
#
paired_anno_cols <- c("anno_SD", "anno_Unique")
keep_cols <- c("hap", "ID", paired_anno_cols)
df <- df[anno_SD > 0 | anno_Unique > 0] %>%
    select(all_of(keep_cols)) %>%
    pivot_longer(all_of(paired_anno_cols),
        names_to = "region",
        names_prefix = "anno_"
    ) %>%
    filter(value > 0) %>%
    data.table()

#
# read in region sizes
#
sizes_f <- "results/annotation/annotation_sizes.tbl"
sizes_f <- snakemake@input[[2]]
region_sizes <- fread(sizes_f) %>%
    mutate(region = gsub("_", "+", gsub("_size", "", anno))) %>%
    data.table()

out_df <- df %>%
    merge(region_sizes, by = c("hap", "region")) %>%
    # , allow.cartesian = TRUE) %>%
    group_by(hap, region) %>%
    summarize(
        "# SNVs" = n(),
        Mbp = unique(size) / 1e6,
        "# SNVs per 10 kbp" = 1e4 * n() / unique(size),
    )

out_df_wide <- out_df %>%
    pivot_wider(
        names_from = region,
        values_from = c("# SNVs", "Mbp", "# SNVs per 10 kbp"),
        names_sort = TRUE,
        names_glue = "{region} {.value}",
    ) %>%
    select(order(colnames(.))) %>%
    data.table()

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

fileConn <- file(paste0(outfile3, "/index.html"))
x <- kable(out_df,
    format = "html",
    align = "l",
    digits = 2,
    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/Unique.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}

Jobs 1

Input files

results/tables/snv_per_kbp.tbl
config/hprc_year1_sample_metadata.txt

Output files

results/figures/violin.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
r-knitr

Code

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

infile <- "results/tables/snv_per_kbp.tbl"
infile2 <- ".test/hprc_year1_sample_metadata.txt"

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

pop <- fread(infile2, fill = TRUE)

df <- fread(infile) %>%
    mutate(Sample = gsub("_(1|2)", "", hap)) %>%
    merge(pop, by = "Sample", all.x = T)

p <- df %>%
    ggplot(
        aes(x = region, y = `# SNVs per 10 kbp`, color = region, fill = region)
    ) +
    geom_violin(alpha = 0.5) +
    geom_jitter(width = 0.2) +
    facet_row(~Superpopulation) +
    scale_fill_manual(values = c(COLOR1, COLOR2)) +
    scale_color_manual(values = c(COLOR1, COLOR2)) +
    theme_cowplot() +
    theme(legend.position = "none") +
    xlab("Genomic region")

ggsave(outfile,
    plot = p,
    width = 9, height = 4,
)

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

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 <- p2
print(height)
ggsave(snakemake@output[[1]],
    plot = fig,
    width = 16, height = height,
    units = "in",
    limitsize = FALSE
)

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 94

Input files

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

Output files

temp/add_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}_snvs_haplotype_coverage.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.windows} | head -n 1 || :)
        HEADER="${{HEADER}}	num_snv	hap"
        echo $HEADER

        gunzip -c {input.windows}             | grep "{wildcards.sm}_{wildcards.h}"             | bedtools coverage                 -a -                 -b {input.snv}                 -counts -sorted             | sed "s/$/\t{wildcards.sm}_{wildcards.h}/"             | sed "1s/^/${{HEADER}}\n/"             | csvtk cut -C "$" -tT -f -haps             > {output}

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

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 = "(Minimum 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)

Workflow

Figures

Tables

Statistics