VariantCalling.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(data.table)\n",
    "library(reshape2)\n",
    "library(plyr)\n",
    "library(ggplot2)\n",
    "library(dplyr)\n",
    "library(stringr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_dir <- \"round1_plots\"\n",
    "if (!dir.exists(output_dir)) {\n",
    "    dir.create(output_dir)\n",
    "  }\n",
    "\n",
    "dir <- \"round1\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "file_paths"
    ]
   },
   "outputs": [],
   "source": [
    "counts_horizontal <- paste(dir, \"/cromwell-executions/PlotAssemblies/a0ca5fcf-0de9-45d6-a869-d9dfe02aab42/call-make_plots/inputs/781347057/write_lines_ebd5d08c92e59058d17818c89d633bc4.tmp\", sep=\"/\")\n",
    "str_typeCount <- paste(dir, \"/cromwell-executions/PlotAssemblies/a0ca5fcf-0de9-45d6-a869-d9dfe02aab42/call-make_plots/inputs/781347057/write_lines_befdd0a3f211c9c7a606cb786c2ef863.tmp\", sep=\"/\")\n",
    "nonRep_typeCount <- paste(dir, \"/cromwell-executions/PlotAssemblies/a0ca5fcf-0de9-45d6-a869-d9dfe02aab42/call-make_plots/inputs/781347057/write_lines_f545feba9f416c5a8b503021d0ccf29e.tmp\", sep=\"/\")\n",
    "nonRep_lengths <- paste(dir, \"/cromwell-executions/PlotAssemblies/a0ca5fcf-0de9-45d6-a869-d9dfe02aab42/call-make_plots/inputs/781347057/write_lines_b02c8c823386764afc19702f118d72b0.tmp\", sep=\"/\")\n",
    "sv_str_VennInput <- paste(dir, \"/cromwell-executions/PlotAssemblies/a0ca5fcf-0de9-45d6-a869-d9dfe02aab42/call-make_plots/inputs/781347057/write_lines_8fd3b25c2db5902381597cca500ae5c7.tmp\", sep=\"/\")\n",
    "str_lengths <- paste(dir, \"/cromwell-executions/PlotAssemblies/a0ca5fcf-0de9-45d6-a869-d9dfe02aab42/call-make_plots/inputs/781347057/write_lines_8c91babb438271a4a5ed61658f765e07.tmp\", sep=\"/\")\n",
    "happy_extended <- paste(dir, \"/cromwell-executions/PlotAssemblies/a0ca5fcf-0de9-45d6-a869-d9dfe02aab42/call-make_plots/inputs/781347057/write_lines_51d18ba14618c9201abe7dd4c349537f.tmp\", sep=\"/\")\n",
    "genomecov_noGaps <- paste(dir, \"/cromwell-executions/PlotAssemblies/a0ca5fcf-0de9-45d6-a869-d9dfe02aab42/call-make_plots/inputs/781347057/write_lines_30955bb5b8b013be558f1085616415b0.tmp\", sep=\"/\")\n",
    "segDup_typeCount <- paste(dir, \"/cromwell-executions/PlotAssemblies/a0ca5fcf-0de9-45d6-a869-d9dfe02aab42/call-make_plots/inputs/781347057/write_lines_afac60d03a7738fd00046629cb726fdf.tmp\", sep=\"/\")\n",
    "indel_lengths_horizontal <- paste(dir, \"/cromwell-executions/PlotAssemblies/a0ca5fcf-0de9-45d6-a869-d9dfe02aab42/call-make_plots/inputs/781347057/write_lines_3ac88ce248d9a14dcb45ec84e0b76397.tmp\", sep=\"/\")\n",
    "happy_het_counts_horizontal <- paste(dir, \"/cromwell-executions/PlotAssemblies/a0ca5fcf-0de9-45d6-a869-d9dfe02aab42/call-make_plots/inputs/781347057/write_lines_ce8f7658892f3bf96781581c6b143878.tmp\", sep=\"/\")\n",
    "segDup_lengths <- paste(dir, \"/cromwell-executions/PlotAssemblies/a0ca5fcf-0de9-45d6-a869-d9dfe02aab42/call-make_plots/inputs/781347057/write_lines_7e4f72479bf2a79fb395e8c3acb2f461.tmp\", sep=\"/\")\n",
    "sv_segDup_VennInput <- paste(dir, \"/cromwell-executions/PlotAssemblies/a0ca5fcf-0de9-45d6-a869-d9dfe02aab42/call-make_plots/inputs/781347057/write_lines_c42dc6aeb8c1cd59cac11f6b46f19019.tmp\", sep=\"/\")\n",
    "indel_counts <- paste(dir, \"/cromwell-executions/PlotAssemblies/a0ca5fcf-0de9-45d6-a869-d9dfe02aab42/call-make_plots/inputs/781347057/write_lines_9c59639608aa522cd2dad0c069615972.tmp\", sep=\"/\")\n",
    "sv_counts_horizontal <- paste(dir, \"/cromwell-executions/PlotAssemblies/a0ca5fcf-0de9-45d6-a869-d9dfe02aab42/call-make_plots/inputs/781347057/write_lines_b59c4353065cb5e37e21711f924e0469.tmp\", sep=\"/\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "populations <- read.table(file=pop, header=FALSE, sep=\"\\t\", col.names=c(\"Sample\", \"Population\", \"Superpopulation\", \"Sex\"))\n",
    "small_variants_files <- read.table(file=happy_extended, header=FALSE, col.names=c(\"file\"))\n",
    "indels_files <- read.table(file=indel_counts, header=FALSE, col.names=c(\"file\"))\n",
    "sv_files <- read.table(file=sv_counts_horizontal, header=FALSE, col.names=c(\"file\"))\n",
    "genome_cov_files <- read.table(file=genomecov_noGaps, header=FALSE, col.names=c(\"file\"))\n",
    "nonRep_lengths_files <- read.table(file=nonRep_lengths, header=FALSE, col.names=c(\"file\"))\n",
    "str_lengths_files <- read.table(file=str_lengths, header=FALSE, col.names=c(\"file\"))\n",
    "segDup_lengths_files <- read.table(file=segDup_lengths, header=FALSE, col.names=c(\"file\"))\n",
    "counts_files <- read.table(file=counts_horizontal, header=FALSE, col.names=c(\"file\"))\n",
    "nonRep_types_files <- read.table(file=nonRep_typeCount, header=FALSE, col.names=c(\"file\"))\n",
    "str_types_files <- read.table(file=str_typeCount, header=FALSE, col.names=c(\"file\"))\n",
    "segDup_types_files <- read.table(file=segDup_typeCount, header=FALSE, col.names=c(\"file\"))\n",
    "het_fates_files <- read.table(file=happy_het_counts_horizontal, header=FALSE, col.names=c(\"file\"))\n",
    "str_venn_files <- read.table(file=sv_str_VennInput, header=FALSE, col.names=c(\"file\"))\n",
    "segDup_venn_files <- read.table(file=sv_segDup_VennInput, header=FALSE, col.names=c(\"file\"))\n",
    "cigar_indel_lengths_file <- read.table(file=indel_lengths_horizontal, header=FALSE, col.names=c(\"file\"))\n",
    "assembly_file <- read.table(file=assembly_tsv, header=FALSE, col.names=c(\"Name\", \"Fasta1\", \"Fasta2\", \"Truth\", \"ID\", \"DataType\"), sep=\"\\t\")\n",
    "assemblies <- assembly_file[,c(\"Name\", \"ID\", \"DataType\")]\n",
    "assemblies$Squashed <- rep(FALSE, length(assemblies$Name))\n",
    "assemblies$Squashed[grepl(\"empty.fa\", assembly_file$Fasta2)] <- TRUE\n",
    "assemblies$ID <- as.character(assemblies$ID)\n",
    "assemblies$ID[assemblies$Squashed] <- paste(assemblies$ID[assemblies$Squashed], \"*\", sep=\"\")\n",
    "assemblies$ID <- factor(assemblies$ID, levels=str_sort(unique(assemblies$ID), numeric=TRUE))\n",
    "\n",
    "tsv_final <- assemblies[,c(\"Name\", \"ID\", \"DataType\")]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "cigar_indel_lengths <- data.frame()\n",
    "for (in_file in cigar_indel_lengths_file$file) {\n",
    "    assembly <- gsub(\".indel_lengths.horizontal.txt\", \"\", basename(in_file))\n",
    "    df <- read.table(file=in_file, header=TRUE, sep=\"\\t\")\n",
    "    df$Assembly <- assembly\n",
    "    cigar_indel_lengths <- rbind(cigar_indel_lengths, df)\n",
    "}\n",
    "cigar_indel_lengths$Metric = factor(cigar_indel_lengths$Metric, levels=c(\"1 bp\", \"2 bp\", \"3-10 bp\", \"11-50 bp\", \"51-100 bp\", \"101-1000 bp\", \"1-10 kb\", \"10 kb+\"))\n",
    "\n",
    "cigar_indel_lengths <- merge(cigar_indel_lengths, assemblies, by.x=\"Assembly\", by.y=\"Name\")\n",
    "ggplot(cigar_indel_lengths, aes(x=ID, y=Value, fill=DataType)) +\n",
    "    geom_bar(stat=\"identity\", position=\"dodge\") +\n",
    "    facet_grid(Metric~Type, scales=\"free\") +\n",
    "    theme_bw() +\n",
    "    theme(axis.text.x = element_text(angle=45, hjust=1, size=7)) +\n",
    "    xlab(\"ID (*=Squashed assembly)\") +\n",
    "    ggtitle(\"Counts of indel variants detected in cigar string by size\")\n",
    "ggsave(paste(output_dir, \"indel_cigar_counts.png\", sep=\"/\"))\n",
    "\n",
    "short_deletions <- cigar_indel_lengths %>%\n",
    "    filter(Type==\"del\") %>% filter(Metric == \"1 bp\" | Metric == \"2 bp\" | Metric == \"3-10 bp\" | Metric == \"11-50 bp\") %>%\n",
    "    group_by(ID) %>% summarize(sum=sum(Value))\n",
    "tsv_final <- merge(tsv_final, short_deletions %>%\n",
    "                   select(ID, sum) %>%\n",
    "                   rename(deletions=sum), by=\"ID\")\n",
    "short_insertions <- cigar_indel_lengths %>%\n",
    "    filter(Type==\"ins\") %>% filter(Metric == \"1 bp\" | Metric == \"2 bp\" | Metric == \"3-10 bp\" | Metric == \"11-50 bp\") %>%\n",
    "    group_by(ID) %>% summarize(sum=sum(Value))\n",
    "tsv_final <- merge(tsv_final, short_insertions %>%\n",
    "                   select(ID, sum) %>%\n",
    "                   rename(insertions=sum), by=\"ID\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "counts <- data.frame()\n",
    "for (in_file in counts_files$file) {\n",
    "    assembly <- gsub(\".counts.horizontal.txt\", \"\", basename(in_file))\n",
    "    df <- read.table(file=in_file, header=FALSE, sep=\"\\t\", col.names=c(\"Metric\", \"Value\"))\n",
    "    df$Assembly <- assembly\n",
    "    counts <- rbind(counts, df)\n",
    "}\n",
    "\n",
    "counts <- merge(counts, assemblies, by.x=\"Assembly\", by.y=\"Name\")\n",
    "ggplot(counts %>% filter(str_detect(Metric, \"^Cigar\")) %>% filter(!str_detect(Metric, \"bp\")), aes(x=ID, y=Value, fill=DataType)) +\n",
    "    geom_bar(stat=\"identity\", position=\"dodge\") +\n",
    "    facet_grid(Metric~., scales=\"free\") +\n",
    "    theme_bw() +\n",
    "    theme(axis.text.x = element_text(angle=45, hjust=1, size=7)) +\n",
    "    xlab(\"ID (*=Squashed assembly)\") +\n",
    "    ggtitle(\"Counts of variants detected in cigar string\")\n",
    "ggsave(paste(output_dir, \"cigar_counts.png\", sep=\"/\"))\n",
    "\n",
    "ggplot(counts %>% filter(str_detect(Metric, \"^Split\")), aes(x=ID, y=Value, fill=DataType)) +\n",
    "    geom_bar(stat=\"identity\", position=\"dodge\") +\n",
    "    facet_grid(Metric~., scales=\"free\") +\n",
    "    theme_bw() +\n",
    "    theme(axis.text.x = element_text(angle=45, hjust=1, size=7)) +\n",
    "    xlab(\"ID (*=Squashed assembly)\") +\n",
    "    ggtitle(\"Counts of variants detected from split alignments\")\n",
    "ggsave(paste(output_dir, \"split_counts.png\", sep=\"/\"))\n",
    "\n",
    "ggplot(counts %>% filter(str_detect(Metric, \"^total\")) %>% filter(!str_detect(Metric, \"our\")), aes(x=ID, y=Value, fill=DataType)) +\n",
    "    geom_bar(stat=\"identity\", position=\"dodge\") +\n",
    "    facet_wrap(\"Metric\", scales=\"free\") +\n",
    "    theme_bw() +\n",
    "    theme(axis.text.x = element_text(angle=45, hjust=1, size=7)) +\n",
    "    xlab(\"ID (*=Squashed assembly)\") +\n",
    "    ggtitle(\"Total 50bp+ SV calls combined cigar and split alignments\")\n",
    "ggsave(paste(output_dir, \"total_sv_counts.png\", sep=\"/\"))\n",
    "\n",
    "tsv_final <- merge(tsv_final, counts %>% filter(Metric==\"Cigar substitutions\") %>%\n",
    "                   select(ID, Value) %>%\n",
    "                   rename(snv_count=Value), by=\"ID\")\n",
    "tsv_final <- merge(tsv_final, counts %>% filter(Metric==\"total combined calls\") %>%\n",
    "                   select(ID, Value) %>%\n",
    "                   rename(sv_count=Value), by=\"ID\")\n",
    "tsv_final <- merge(tsv_final, counts %>% filter(Metric==\"total nonRep combined calls\") %>%\n",
    "                   select(ID, Value) %>%\n",
    "                   rename(sv_nonrep_count=Value), by=\"ID\")\n",
    "tsv_final <- merge(tsv_final, counts %>% filter(Metric==\"total str combined calls\") %>%\n",
    "                   select(ID, Value) %>%\n",
    "                   rename(sv_str_count=Value), by=\"ID\")\n",
    "tsv_final <- merge(tsv_final, counts %>% filter(Metric==\"total segDup combined calls\") %>%\n",
    "                   select(ID, Value) %>%\n",
    "                   rename(sv_segdup_count=Value), by=\"ID\")\n",
    "head(tsv_final)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "small_variants <- data.frame()\n",
    "for (in_file in small_variants_files$file) {\n",
    "    assembly <- gsub(\".happy.extended.csv\", \"\", basename(in_file))\n",
    "    df <- read.table(file=in_file, header=TRUE, sep=\",\")\n",
    "    df$Assembly <- assembly\n",
    "    small_variants <- rbind(small_variants, df)\n",
    "}\n",
    "small_variants$Subset <- recode(small_variants$Subset, \"*\"=\"ALL\")\n",
    "small_variants$Type <- factor(small_variants$Type, levels=c(\"SNP\", \"INDEL\"))\n",
    "small_variants <- merge(small_variants, assemblies, by.x=\"Assembly\", by.y=\"Name\")\n",
    "small_variants$variant_discordance <- (small_variants$TRUTH.FN + small_variants$QUERY.FP)/(small_variants$TRUTH.TP + small_variants$TRUTH.FN + small_variants$QUERY.FP)\n",
    "small_variants$detection_sensitivity <- (small_variants$TRUTH.TP + small_variants$FP.gt)/small_variants$TRUTH.TOTAL\n",
    "small_variants$detection_specificity <- (small_variants$QUERY.TP + small_variants$FP.gt)/small_variants$QUERY.TOTAL\n",
    "\n",
    "ggplot(small_variants %>% filter(Subset == \"ALL\" & Subtype==\"*\") %>% filter(Filter == \"PASS\" | Filter==\"*\"), aes(x=ID, y=METRIC.Recall, fill=DataType)) +\n",
    "    geom_bar(stat=\"identity\", position=\"dodge\") +\n",
    "    facet_grid(Type~.) +\n",
    "    theme_bw() +\n",
    "    theme(axis.text.x = element_text(angle=45, hjust=1, size=7)) +\n",
    "    xlab(\"ID (*=Squashed assembly)\") +\n",
    "    ggtitle(\"Hap.py sensitivity to truthset\")\n",
    "ggsave(paste(output_dir, \"happy.sensitivity.png\", sep=\"/\"))\n",
    "\n",
    "ggplot(small_variants %>% filter(Subset == \"TS_contained\" & Subtype==\"*\") %>% filter(Filter == \"PASS\" | Filter==\"*\"), aes(x=ID, y=METRIC.Precision, fill=DataType)) +\n",
    "    geom_bar(stat=\"identity\", position=\"dodge\") +\n",
    "    facet_grid(Type~.) +\n",
    "    theme_bw() +\n",
    "    theme(axis.text.x = element_text(angle=45, hjust=1, size=7)) +\n",
    "    xlab(\"ID (*=Squashed assembly)\") +\n",
    "    ggtitle(\"Hap.py specificity in confident region\")\n",
    "ggsave(paste(output_dir, \"happy.specificity.png\", sep=\"/\"))\n",
    "\n",
    "ggplot(small_variants %>% filter(Subset == \"ALL\" & Subtype==\"*\") %>% filter(Filter == \"PASS\" | Filter==\"*\"), aes(x=ID, y=detection_sensitivity, fill=DataType)) +\n",
    "    geom_bar(stat=\"identity\", position=\"dodge\") +\n",
    "    facet_grid(Type~.) +\n",
    "    theme_bw() +\n",
    "    theme(axis.text.x = element_text(angle=45, hjust=1, size=7)) +\n",
    "    xlab(\"ID (*=Squashed assembly)\") +\n",
    "    ggtitle(\"Hap.py sensitivity to truthset (detection only, no gt)\")\n",
    "ggsave(paste(output_dir, \"happy.sensitivity_detection.png\", sep=\"/\"))\n",
    "\n",
    "ggplot(small_variants %>% filter(Subset == \"TS_contained\" & Subtype==\"*\") %>% filter(Filter == \"PASS\" | Filter==\"*\"), aes(x=ID, y=detection_specificity, fill=DataType)) +\n",
    "    geom_bar(stat=\"identity\", position=\"dodge\") +\n",
    "    facet_grid(Type~.) +\n",
    "    theme_bw() +\n",
    "    theme(axis.text.x = element_text(angle=45, hjust=1, size=7)) +\n",
    "    xlab(\"ID (*=Squashed assembly)\") +\n",
    "    ggtitle(\"Hap.py specificity in confident region (detection only, no gt)\")\n",
    "ggsave(paste(output_dir, \"happy.specificity_detection.png\", sep=\"/\"))\n",
    "\n",
    "ggplot(small_variants %>% filter(Subset == \"TS_contained\" & Subtype==\"*\") %>% filter(Filter == \"PASS\" | Filter==\"*\"), aes(x=ID, y=variant_discordance, fill=DataType)) +\n",
    "    geom_bar(stat=\"identity\", position=\"dodge\") +\n",
    "    facet_grid(Type~.) +\n",
    "    theme_bw() +\n",
    "    theme(axis.text.x = element_text(angle=45, hjust=1, size=7)) +\n",
    "    xlab(\"ID (*=Squashed assembly)\") +\n",
    "    ggtitle(\"Variant discordance (TRUTH.FN + QUERY.FP)/(TRUTH.TP + TRUTH.FN + QUERY.FP)\")\n",
    "ggsave(paste(output_dir, \"happy.variant_discordance.png\", sep=\"/\"))\n",
    "\n",
    "tsv_final <- merge(tsv_final, small_variants[,c(\"ID\", \"detection_sensitivity\", \"detection_specificity\", \"Subset\", \"Subtype\", \"Filter\", \"Type\")] %>% \n",
    "                   filter(Subset == \"TS_contained\" & Subtype==\"*\") %>% \n",
    "                   filter(Filter == \"PASS\" | Filter==\"*\") %>%\n",
    "                   filter(Type == \"SNP\") %>%\n",
    "                   select(ID, detection_sensitivity, detection_specificity) %>%\n",
    "                   rename(snv_sensitivity = detection_sensitivity) %>%\n",
    "                   rename(snv_specificity = detection_specificity), by=\"ID\")\n",
    "                   \n",
    "tsv_final <- merge(tsv_final, small_variants[,c(\"ID\", \"detection_sensitivity\", \"detection_specificity\", \"Subset\", \"Subtype\", \"Filter\", \"Type\")] %>% \n",
    "                   filter(Subset == \"TS_contained\" & Subtype==\"*\") %>% \n",
    "                   filter(Filter == \"PASS\" | Filter==\"*\") %>%\n",
    "                   filter(Type == \"INDEL\") %>%\n",
    "                   select(ID, detection_sensitivity, detection_specificity) %>%\n",
    "                   rename(indel_sensitivity = detection_sensitivity) %>%\n",
    "                   rename(indel_specificity = detection_specificity), by=\"ID\")\n",
    "head(tsv_final)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "indels <- data.frame()\n",
    "for (in_file in indels_files$file) {\n",
    "    assembly <- gsub(\".indels.indel.counts.horizontal.txt\", \"\", basename(in_file))\n",
    "    df <- read.table(file=in_file, header=FALSE, sep=\"\\t\", col.names=c(\"Metric\", \"Value\"))\n",
    "    df$Assembly <- assembly\n",
    "    indels <- rbind(indels, df)\n",
    "}\n",
    "indels$slop <- 50\n",
    "indels$slop[grepl(\"10 bp\", indels$Metric)] <- 10\n",
    "indels$slop[grepl(\"1 bp\", indels$Metric)] <- 1\n",
    "truth_name <- \"GiaB_HG002\"\n",
    "indels$Metric <- recode(indels$Metric,\n",
    "        \"total GiaB_HG002 het indels\"=\"totalTruthHet\",\n",
    "        \"pairtopair for GiaB_HG002 indels vs our calls, 10 bp slop\"=\"truthAllMatched\",\n",
    "        \"pairtopair for GiaB_HG002 indels vs our calls, 50 bp slop\"=\"truthAllMatched\",\n",
    "        \"pairtopair for GiaB_HG002 indels vs our calls, 1 bp slop\"=\"truthAllMatched\",\n",
    "        \"pairtopair for GiaB_HG002 het to homalt, 50 bp slop\"=\"truthHetToHomalt\",\n",
    "        \"pairtopair for GiaB_HG002 het to homalt, 10 bp slop\"=\"truthHetToHomalt\",\n",
    "        \"pairtopair for GiaB_HG002 het to homalt, 1 bp slop\"=\"truthHetToHomalt\",\n",
    "        \"pairtopair for GiaB_HG002 het indels vs our calls, 50 bp slop\"=\"truthHetMatched\",\n",
    "        \"pairtopair for GiaB_HG002 het indels vs our calls, 10 bp slop\"=\"truthHetMatched\",\n",
    "        \"pairtopair for GiaB_HG002 het indels vs our calls, 1 bp slop\"=\"truthHetMatched\",\n",
    "        \"pairtopair for GiaB_HG002 het to het, 50 bp slop\"=\"truthHetToHet\",\n",
    "        \"pairtopair for GiaB_HG002 het to het, 10 bp slop\"=\"truthHetToHet\",\n",
    "        \"pairtopair for GiaB_HG002 het to het, 1 bp slop\"=\"truthHetToHet\")\n",
    "for (assembly in unique(indels$Assembly)) {\n",
    "    row <- indels[indels$Assembly==assembly & indels$Metric==\"totalTruthHet\",]\n",
    "    row$slop <- 10\n",
    "    indels <- rbind(indels, row)\n",
    "    row$slop <- 1\n",
    "    indels <- rbind(indels, row)\n",
    "}\n",
    "indels2 <- dcast(indels, Assembly+slop~Metric, value.var=\"Value\")\n",
    "indels2$truthHetToHomalt <- indels2$truthHetMatched - indels2$truthHetToHet\n",
    "indels2$truthHetToHomref <- indels2$totalTruthHet - indels2$truthHetMatched\n",
    "indels <- melt(indels2, id.vars=c(\"Assembly\", \"slop\"), variable.name=\"Metric\", value.name=\"Value\")\n",
    "\n",
    "indels <- merge(indels, assemblies, by.x=\"Assembly\", by.y=\"Name\")\n",
    "ggplot(indels %>% filter(Metric == \"truthHetToHomalt\" | Metric == \"truthHetToHet\" | Metric == \"truthHetToHomref\"), aes(x=ID, y=Value, fill=Metric)) +\n",
    "    geom_bar(stat=\"identity\", position=\"stack\") +\n",
    "    facet_grid(slop~Squashed, scales=\"free\") +\n",
    "    theme_bw() +\n",
    "    theme(axis.text.x = element_text(angle=45, hjust=1, size=5)) +\n",
    "    xlab(\"ID (*=Squashed assembly)\") +\n",
    "    ggtitle(\"Fates of heterozygous indels from truthset\")\n",
    "ggsave(paste(output_dir, \"truthHetIndelFates.png\", sep=\"/\"), width=8, height=6, units=\"in\")\n",
    "\n",
    "ggplot(indels %>% filter(Metric == \"truthHetToHomalt\" | Metric == \"truthHetToHomref\"), aes(x=ID, y=Value, fill=Metric)) +\n",
    "    geom_bar(stat=\"identity\", position=\"stack\") +\n",
    "    facet_grid(slop~Squashed, scales=\"free\") +\n",
    "    theme_bw() +\n",
    "    theme(axis.text.x = element_text(angle=45, hjust=1, size=6)) +\n",
    "    xlab(\"ID (*=Squashed assembly)\") +\n",
    "    ggtitle(\"Fates of heterozygous indels from truthset\")\n",
    "ggsave(paste(output_dir, \"truthHetIndelFatesErrorsOnly.png\", sep=\"/\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "sv <- data.frame()\n",
    "for (in_file in sv_files$file) {\n",
    "    assembly <- gsub(\".sv.counts.horizontal.txt\", \"\", basename(in_file))\n",
    "    df <- read.table(file=in_file, header=FALSE, sep=\"\\t\", col.names=c(\"Metric\", \"Value\"))\n",
    "    df$Assembly <- assembly\n",
    "    sv <- rbind(sv, df)\n",
    "}\n",
    "\n",
    "sv$Region <- \"ALL\"\n",
    "sv$Region[grepl(\"str\", sv$Metric)] <- \"STR\"\n",
    "sv$Region[grepl(\"segDup\", sv$Metric)] <- \"SEGDUP\"\n",
    "sv$Region[grepl(\"nonRep\", sv$Metric)] <- \"NONREP\"\n",
    "\n",
    "sv$Filter <- \"ALL\"\n",
    "sv$Filter[grepl(\"PASS\", sv$Metric)] <- \"PASS\"\n",
    "\n",
    "sv$Metric <- gsub(\"str \", \"\", sv$Metric)\n",
    "sv$Metric <- gsub(\"nonRep \", \"\", sv$Metric)\n",
    "sv$Metric <- gsub(\"segDup \", \"\", sv$Metric)\n",
    "sv$Metric <- gsub(\"str\", \"\", sv$Metric)\n",
    "sv$Metric <- gsub(\"nonRep\", \"\", sv$Metric)\n",
    "sv$Metric <- gsub(\"segDup\", \"\", sv$Metric)\n",
    "\n",
    "sv$Metric <- gsub(\"PASS \", \"\", sv$Metric)\n",
    "sv$Metric <- gsub(\"PASS\", \"\", sv$Metric)\n",
    "sv$Metric <- gsub(\"for \", \"\", sv$Metric)\n",
    "sv$Metric <- gsub(\", 50 bp slop\", \"\", sv$Metric)\n",
    "sv$Metric <- gsub(\" $\", \"\", sv$Metric)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "sv$Overlap <- \"pairtopair 50 bp slop\"\n",
    "sv$Overlap[grepl(\"50\", sv$Metric)] <- \"reciprocal 50% overlap\"\n",
    "sv$Overlap[grepl(\"10\", sv$Metric)] <- \"reciprocal 10% overlap\"\n",
    "sv$Metric <- gsub(\"%\", \"\", sv$Metric)\n",
    "sv$Metric <- gsub(\"50 \", \"\", sv$Metric)\n",
    "sv$Metric <- gsub(\"10 \", \"\", sv$Metric)\n",
    "\n",
    "sv$Region <- recode(sv$Region,\n",
    "                    NONREP=paste(\"NONREP n=\", unique(sv$Value[sv$Metric==\"total GiaB_HG002 calls\" & sv$Region==\"NONREP\" & sv$Filter==\"PASS\"]), sep=\"\"),\n",
    "                    SEGDUP=paste(\"SEGDUP n=\", unique(sv$Value[sv$Metric==\"total GiaB_HG002 calls\" & sv$Region==\"SEGDUP\" & sv$Filter==\"PASS\"]), sep=\"\"),\n",
    "                    STR=paste(\"STR n=\", unique(sv$Value[sv$Metric==\"total GiaB_HG002 calls\" & sv$Region==\"STR\" & sv$Filter==\"PASS\"]), sep=\"\"))\n",
    "\n",
    "sv$Overlap[sv$Metric==\"total GiaB_HG002 bed\"] <- \"reciprocal 50% overlap\"\n",
    "rows <- sv[sv$Metric==\"total GiaB_HG002 bed\",]\n",
    "rows$Overlap <- \"reciprocal 10% overlap\"\n",
    "sv <- rbind(sv, rows)\n",
    "\n",
    "rows <- sv[sv$Metric==\"SVs in confident region intersect\",]\n",
    "rows$Filter <- \"PASS\"\n",
    "sv <- rbind(sv, rows)\n",
    "\n",
    "rows <- sv[sv$Metric==\"SVs in confident region both ends\",]\n",
    "rows$Filter <- \"PASS\"\n",
    "sv <- rbind(sv, rows)\n",
    "\n",
    "rows <- sv[sv$Metric==\"SVs in confident region coverage gt .9\",]\n",
    "rows$Filter <- \"PASS\"\n",
    "sv <- rbind(sv, rows)\n",
    "\n",
    "sv2 <- dcast(sv, Assembly+Region+Filter+Overlap~Metric, value.var=\"Value\")\n",
    "sv2$\"pairtopair GiaB_HG002 vs combined calls\"[sv2$Region==\"ALL\" & sv2$Filter==\"PASS\" & sv2$Overlap==\"pairtopair 50 bp slop\"] <-\n",
    "    sv2$\"pairtopair GiaB_HG002 vs combined calls\"[sv2$Region==\"NONREP n=5652\" & sv2$Filter==\"PASS\" & sv2$Overlap==\"pairtopair 50 bp slop\"] +\n",
    "    sv2$\"pairtopair GiaB_HG002 vs combined calls\"[sv2$Region==\"SEGDUP n=48\" & sv2$Filter==\"PASS\" & sv2$Overlap==\"pairtopair 50 bp slop\"] +\n",
    "    sv2$\"pairtopair GiaB_HG002 vs combined calls\"[sv2$Region==\"STR n=6982\" & sv2$Filter==\"PASS\" & sv2$Overlap==\"pairtopair 50 bp slop\"]\n",
    "sv2$sensitivity <- sv2$\"pairtopair GiaB_HG002 vs combined calls\"/sv2$\"total GiaB_HG002 calls\"\n",
    "sv2$sensitivity[is.na(sv2$sensitivity)] <- sv2$\"Reciprocal overlap GiaB_HG002\"[is.na(sv2$sensitivity)]/sv2$\"total GiaB_HG002 bed\"[is.na(sv2$sensitivity)]\n",
    "sv2$specificity <- sv2$\"pairtopair GiaB_HG002 vs combined calls\"/sv2$\"SVs in confident region both ends\"\n",
    "sv2$count <- sv2$\"pairtopair GiaB_HG002 vs combined calls\"\n",
    "sv2$count[is.na(sv2$count)] <- sv2$\"Reciprocal overlap GiaB_HG002\"[is.na(sv2$count)]\n",
    "\n",
    "sv2 <- merge(sv2, assemblies, by.x=\"Assembly\", by.y=\"Name\")\n",
    "ggplot(sv2 %>% filter(Region != \"ALL\") %>% filter(Filter == \"PASS\"), aes(x=ID, y=count, fill=DataType)) +\n",
    "    geom_bar(stat=\"identity\", position=\"dodge\") +\n",
    "    facet_grid(Region~Overlap, scales=\"free\") +\n",
    "    theme_bw() +\n",
    "    theme(axis.text.x = element_text(angle=45, hjust=1, size=6)) +\n",
    "    xlab(\"ID (*=Squashed assembly)\") +\n",
    "    ggtitle(\"SV sensitivity counts\")\n",
    "ggsave(paste(output_dir, \"sv.sensitivityCounts.png\", sep=\"/\"), width=10, height=10, units=\"in\")\n",
    "\n",
    "ggplot(sv2 %>% filter(Region != \"ALL\") %>% filter(Filter == \"PASS\"), aes(x=ID, y=sensitivity, fill=DataType)) +\n",
    "    geom_bar(stat=\"identity\", position=\"dodge\") +\n",
    "    facet_grid(Region~Overlap, scales=\"free\") +\n",
    "    theme_bw() +\n",
    "    theme(axis.text.x = element_text(angle=45, hjust=1, size=6)) +\n",
    "    xlab(\"ID (*=Squashed assembly)\") +\n",
    "    ggtitle(\"SV sensitivity\")\n",
    "ggsave(paste(output_dir, \"sv.sensitivity.png\", sep=\"/\"), width=10, height=10, units=\"in\")\n",
    "\n",
    "tsv_final <- merge(tsv_final, sv2 %>% filter(Region == \"ALL\") %>% filter(Filter == \"PASS\") %>%\n",
    "                   filter(Overlap == \"pairtopair 50 bp slop\") %>%\n",
    "                   filter(Region == \"ALL\") %>%\n",
    "                   select(ID, sensitivity, specificity) %>%\n",
    "                   rename(sv_sensitivity=sensitivity) %>%\n",
    "                   rename(sv_specificity=specificity), by=\"ID\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "spec <- sv2 %>% filter(Region == \"ALL\" & Filter==\"PASS\")\n",
    "spec <- spec %>% select(ID, DataType, \"SVs in confident region both ends\", \"SVs in confident region intersect\", \"SVs in confident region coverage gt .9\", \"pairtopair GiaB_HG002 vs combined calls\")\n",
    "spec <- melt(spec, id.vars=c(\"ID\", \"pairtopair GiaB_HG002 vs combined calls\", \"DataType\"))\n",
    "spec$specificity <- spec$\"pairtopair GiaB_HG002 vs combined calls\"/spec$value\n",
    "spec <- spec %>% rename(specificity_method=variable)\n",
    "spec$specificity_method <- recode(spec$specificity_method, \"SVs in confident region both ends\"=\"both ends\")\n",
    "spec$specificity_method <- recode(spec$specificity_method, \"SVs in confident region intersect\"=\"intersect 100%\")\n",
    "spec$specificity_method <- recode(spec$specificity_method, \"SVs in confident region coverage gt .9\"=\"covered >= 90%\")\n",
    "head(spec)\n",
    "ggplot(spec, aes(x=ID, y=specificity, fill=DataType)) +\n",
    "    geom_bar(stat=\"identity\", position=\"dodge\") +\n",
    "    facet_grid(.~specificity_method, scales=\"free\") +\n",
    "    theme_bw() +\n",
    "    theme(axis.text.x = element_text(angle=45, hjust=1, size=6)) +\n",
    "    xlab(\"ID (*=Squashed assembly)\") +\n",
    "    ggtitle(\"SV specificity\")\n",
    "ggsave(paste(output_dir, \"sv.specificity.png\", sep=\"/\"), width=10, height=10, units=\"in\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "genome_cov <- data.frame()\n",
    "cov_summary <- data.frame()\n",
    "for (in_file in genome_cov_files$file) {\n",
    "    assembly <- gsub(\".genomecov.noGaps.hist\", \"\", basename(in_file))\n",
    "    cov <- read.table(file=in_file, header=FALSE, sep=\"\\t\", col.names=c(\"chrom\", \"cov\", \"bp\"))\n",
    "    cov <- cov %>% filter(chrom %in% c(\"chr1\", \"chr2\", \"chr3\", \"chr4\", \"chr5\", \"chr6\", \"chr7\", \"chr8\", \"chr9\", \"chr10\", \"chr11\", \"chr12\", \"chr13\", \"chr14\", \"chr15\", \"chr16\", \"chr17\", \"chr18\", \"chr19\", \"chr20\", \"chr21\", \"chr22\", \"chrX\", \"chrY\"))\n",
    "    gt9 <- cov %>% filter(cov > 9)\n",
    "    lt10 <- cov %>% filter(cov < 10)\n",
    "    gt9_zoomed <- gt9 %>% group_by(chrom) %>% summarize(bp=sum(bp))\n",
    "    gt9_zoomed$cov <- 10\n",
    "    cov_zoomed <- rbind(lt10, gt9_zoomed)\n",
    "    cov_zoomed$Assembly <- assembly\n",
    "    genome_cov <- rbind(genome_cov, cov_zoomed)\n",
    "    ggplot(cov_zoomed, aes(x=cov, y=bp)) +\n",
    "        geom_bar(stat=\"identity\", position=\"dodge\") +\n",
    "        facet_wrap(\"chrom\", scales=\"free\") +\n",
    "        scale_x_discrete(name=\"Coverage\") +\n",
    "        theme_bw() +\n",
    "        ggtitle(assembly)\n",
    "    ggsave(paste(output_dir, paste(assembly, \"no_gaps_coverage.png\", sep=\".\"), sep=\"/\"))\n",
    "    covDip <- cov %>% filter(!grepl(\"X\", chrom)) %>% filter(!grepl(\"Y\", chrom))\n",
    "    gt9 <- covDip %>% filter(cov > 9)\n",
    "    lt10 <- covDip %>% filter(cov < 10)\n",
    "    gt9_zoomed <- gt9 %>% group_by(chrom) %>% summarize(bp=sum(bp))\n",
    "    gt9_zoomed$cov <- 10\n",
    "        covDip_zoomed <- rbind(lt10, gt9_zoomed)\n",
    "    summary_bp <- covDip_zoomed %>% group_by(cov) %>% summarize(sum=as.numeric(sum(as.numeric(bp))))\n",
    "    total_bp <- sum(summary_bp$sum)\n",
    "    summary_bp$percent <- summary_bp$sum/total_bp\n",
    "    summary_bp$Assembly <- assembly\n",
    "    cov_summary <- rbind(cov_summary, summary_bp)\n",
    "\n",
    "    ggplot(summary_bp, aes(x=cov, y=percent)) +\n",
    "        geom_bar(stat=\"identity\", position=\"dodge\") +\n",
    "        scale_x_discrete(name=\"Coverage\") +\n",
    "        theme_bw() +\n",
    "        ggtitle(assembly)\n",
    "    ggsave(paste(output_dir, paste(assembly, \"no_gaps_coverage_summary.png\", sep=\".\"), sep=\"/\"))\n",
    "}\n",
    "genome_cov <- merge(genome_cov, assemblies, by.x=\"Assembly\", by.y=\"Name\")\n",
    "ggplot(genome_cov, aes(x=cov, y=bp)) +\n",
    "    geom_boxplot() +\n",
    "    facet_grid(Squashed~chrom, scales=\"free\") +\n",
    "    scale_x_discrete(name=\"Coverage\") +\n",
    "    theme_bw() +\n",
    "    ggtitle(\"Coverage by chromosome, excludes ref gaps\")\n",
    "ggsave(paste(output_dir, \"no_gaps_coverage.png\", sep=\"/\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "cov_summary <- merge(cov_summary, assemblies, by.x=\"Assembly\", by.y=\"Name\")\n",
    "cov_summary$cov <- as.character(cov_summary$cov)\n",
    "cov_summary$cov[cov_summary$cov==\"10\"] <- \"10+\"\n",
    "cov_summary$cov <- factor(cov_summary$cov, levels=c(\"0\", \"1\", \"2\", \"3\", \"4\", \"5\", \"6\", \"7\", \"8\", \"9\", \"10+\"))\n",
    "cov_summary$Squashed_text <- \"Squashed\"\n",
    "cov_summary$Squashed_text[!cov_summary$Squashed] <- \"Not squashed\"\n",
    "ggplot(cov_summary, aes(x=cov, y=percent)) +\n",
    "    geom_boxplot() +\n",
    "    facet_grid(.~Squashed_text, scales=\"free\") +\n",
    "    scale_x_discrete(name=\"Coverage\") +\n",
    "    theme_bw() +\n",
    "    ggtitle(\"Coverage for chr1-22, excludes ref gaps\")\n",
    "ggsave(paste(output_dir, \"no_gaps_coverage_summary.png\", sep=\"/\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "for (in_file in nonRep_lengths_files$file) {\n",
    "    assembly <- gsub(\".ours.nonRep.lengths.txt\", \"\", basename(in_file))\n",
    "    sv_nonrep <- read.table(in_file, header=FALSE, col.names=c(\"type\", \"length\"))\n",
    "    sv_nonrep$region <- \"nonRep\"\n",
    "    sv_str <- read.table(gsub(\"nonRep\", \"str\", in_file), header=FALSE, col.names=c(\"type\", \"length\"))\n",
    "    sv_str$region <- \"str\"\n",
    "    sv_segdup <- read.table(gsub(\"nonRep\", \"segDup\", in_file), header=FALSE, col.names=c(\"type\", \"length\"))\n",
    "    sv_segdup$region <- \"segDup\"\n",
    "    sv <- rbind(sv_nonrep, sv_str, sv_segdup) %>% filter(type!=\"GAP\")\n",
    "    sv$length_zoom <- sv$length\n",
    "    sv$length_zoom[sv$length_zoom < -10000000] <- -1050.5\n",
    "    sv$length_zoom[sv$length_zoom < -1000000] <- -1040.5\n",
    "    sv$length_zoom[sv$length_zoom < -100000] <- -1030.5\n",
    "    sv$length_zoom[sv$length_zoom < -10000] <- -1020.5\n",
    "    sv$length_zoom[sv$length_zoom < -1000 & sv$length_zoom != -1050.5 & sv$length_zoom != -1040.5 & sv$length_zoom != -1030.5 & sv$length_zoom != -1020.5] <- -1010.5\n",
    "    sv$length_zoom[sv$length_zoom > 1000 & sv$length_zoom < 10000] <- 1010.5\n",
    "    sv$length_zoom[sv$length_zoom > 10000 & sv$length_zoom < 100000] <- 1020.5\n",
    "    sv$length_zoom[sv$length_zoom > 100000 & sv$length_zoom < 1000000] <- 1030.5\n",
    "    sv$length_zoom[sv$length_zoom > 1000000 & sv$length_zoom < 10000000] <- 1040.5\n",
    "    sv$length_zoom[sv$length_zoom > 10000000] <- 1050.5\n",
    "    sv$length_zoom[sv$length_zoom==-1050.5] <- -1050\n",
    "    sv$length_zoom[sv$length_zoom==-1040.5] <- -1040\n",
    "    sv$length_zoom[sv$length_zoom==-1030.5] <- -1030\n",
    "    sv$length_zoom[sv$length_zoom==-1020.5] <- -1020\n",
    "    sv$length_zoom[sv$length_zoom==-1010.5] <- -1010\n",
    "    sv$length_zoom[sv$length_zoom==1010.5] <- 1010\n",
    "    sv$length_zoom[sv$length_zoom==1020.5] <- 1020\n",
    "    sv$length_zoom[sv$length_zoom==1030.5] <- 1030\n",
    "    sv$length_zoom[sv$length_zoom==1040.5] <- 1040\n",
    "    sv$length_zoom[sv$length_zoom==1050.5] <- 1050\n",
    "    ggplot(sv, aes(x=length_zoom)) +\n",
    "        geom_histogram(binwidth=10) +\n",
    "        facet_grid(region~type, scales=\"free\") +\n",
    "        theme_bw() +\n",
    "        theme(axis.text.x = element_text(angle = 45, hjust = 1, size=7)) +\n",
    "        ggtitle(assembly)\n",
    "    ggsave(file=paste(output_dir, paste(assembly, \"sv_lengths.png\", sep=\".\"), sep=\"/\"))\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "svCounts <- data.frame()\n",
    "for (in_file in nonRep_types_files$file) {\n",
    "    assembly <- gsub(\".ours.nonRep.typeCount.txt\", \"\", basename(in_file))\n",
    "    for (type in c(\"nonRep\", \"str\", \"segDup\")) {\n",
    "        df <- read.table(gsub(\"nonRep\", type, in_file), header=FALSE, sep=\" \", col.names=c(\"Value\", \"Metric\"))\n",
    "        df$Assembly <- assembly\n",
    "        df$Type <- type\n",
    "        svCounts <- rbind(svCounts, df)\n",
    "    }\n",
    "}\n",
    "svCounts$Sample <- gsub(\"\\\\..*\", \"\", svCounts$Assembly)\n",
    "svCounts <- merge(svCounts, populations, by=\"Sample\")\n",
    "svCounts <- merge(svCounts, assemblies, by.x=\"Assembly\", by.y=\"Name\")\n",
    "\n",
    "ggplot(svCounts %>% filter(Type==\"nonRep\") %>% filter(!str_detect(Metric, \"^G\")), aes(x=ID, y=Value, fill=DataType)) +\n",
    "    geom_bar(stat=\"identity\", position=\"dodge\") +\n",
    "    facet_wrap(\"Metric\", scales=\"free\") +\n",
    "    theme_bw() +\n",
    "    theme(axis.text.x = element_text(angle=45, hjust=1, size=6)) +\n",
    "    ggtitle(\"Non-repetitive SV counts by type\")\n",
    "    ggsave(paste(output_dir, \"nonrepSvCountsSquashed.png\", sep=\"/\"), width=12, height=8, units=\"in\")\n",
    "\n",
    "ggplot(svCounts %>% filter(Type==\"str\") %>% filter(!str_detect(Metric, \"^G\")), aes(x=ID, y=Value, fill=DataType)) +\n",
    "    geom_bar(stat=\"identity\", position=\"dodge\") +\n",
    "    facet_wrap(\"Metric\", scales=\"free\") +\n",
    "    theme_bw() +\n",
    "    theme(axis.text.x = element_text(angle=45, hjust=1, size=6)) +\n",
    "    ggtitle(\"STR SV counts by type\")\n",
    "ggsave(paste(output_dir, \"strSvCountsSquashed.png\", sep=\"/\"), width=12, height=8, units=\"in\")\n",
    "\n",
    "ggplot(svCounts %>% filter(Type==\"segDup\") %>% filter(!str_detect(Metric, \"^G\")), aes(x=ID, y=Value, fill=DataType)) +\n",
    "    geom_bar(stat=\"identity\", position=\"dodge\") +\n",
    "    facet_wrap(\"Metric\", scales=\"free\") +\n",
    "    theme_bw() +\n",
    "    theme(axis.text.x = element_text(angle=45, hjust=1, size=6)) +\n",
    "    ggtitle(\"SegDup SV counts by type\")\n",
    "ggsave(paste(output_dir, \"segDupSvCountsSquashed.png\", sep=\"/\"), width=12, height=8, units=\"in\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "a <- svCounts %>% filter(!str_detect(Metric, \"^G\")) %>% group_by(ID, Metric) %>% summarize(count=sum(Value))\n",
    "b <- dcast(a, ID~Metric)\n",
    "head(a)\n",
    "head(b)\n",
    "tsv_final <- merge(tsv_final, b, by=\"ID\")\n",
    "head(tsv_final)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "het <- data.frame()\n",
    "for (in_file in het_fates_files$file) {\n",
    "    assembly <- gsub(\".happy.het.counts.horizontal.txt\", \"\", basename(in_file))\n",
    "    df <- read.table(in_file, header=FALSE, col.names=c(\"Type\", \"Origin\", \"Fate\", \"Value\"))\n",
    "    df$Assembly <- assembly\n",
    "    het <- rbind(het, df)\n",
    "}\n",
    "het$Type <- factor(het$Type, levels=c(\"SNP\", \"INDEL\"))\n",
    "het <- merge(het, assemblies, by.x=\"Assembly\", by.y=\"Name\")\n",
    "\n",
    "ggplot(het, aes(x=ID, y=Value, fill=Fate)) +\n",
    "    geom_bar(stat=\"identity\", position=\"stack\") +\n",
    "    facet_grid(Type~Origin, scales=\"free\") +\n",
    "    theme_bw() +\n",
    "    theme(axis.text.x = element_text(angle=45, hjust=1, size=7)) +\n",
    "    ggtitle(\"Fate of variants from truthset\")\n",
    "ggsave(paste(output_dir, \"het_fates.png\", sep=\"/\"))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "het_totals <- het %>% group_by(ID, Origin, Type) %>% summarize(total=sum(Value))\n",
    "het_totals <- dcast(het_totals, ID~Origin+Type) %>%\n",
    "                rename(het_snp_total=HET_SNP) %>% rename(het_indel_total=HET_INDEL) %>%\n",
    "                rename(hom_snp_total=HOMALT_SNP) %>% rename(hom_indel_total=HOMALT_INDEL)\n",
    "het_tp <- het %>% filter(Fate==\"HET\" | Fate==\"HOMALT\") %>% group_by(ID, Origin, Type) %>% summarize(total=sum(Value))\n",
    "het_tp <- dcast(het_tp, ID~Origin+Type) %>%\n",
    "                rename(het_snp_tp=HET_SNP) %>% rename(het_indel_tp=HET_INDEL) %>%\n",
    "                rename(hom_snp_tp=HOMALT_SNP) %>% rename(hom_indel_tp=HOMALT_INDEL)\n",
    "\n",
    "head(het_totals)\n",
    "head(het_tp)\n",
    "tsv_final <- merge(tsv_final, het_totals, by=\"ID\")\n",
    "tsv_final <- merge(tsv_final, het_tp, by=\"ID\")\n",
    "tsv_final$het_snp_sensitivity <- tsv_final$het_snp_tp/tsv_final$het_snp_total\n",
    "tsv_final$het_indel_sensitivity <- tsv_final$het_indel_tp/tsv_final$het_indel_total\n",
    "tsv_final$hom_snp_sensitivity <- tsv_final$hom_snp_tp/tsv_final$hom_snp_total\n",
    "tsv_final$hom_indel_sensitivity <- tsv_final$hom_indel_tp/tsv_final$hom_indel_total"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "rep <- data.frame()\n",
    "for (in_file in str_venn_files$file) {\n",
    "    assembly <- gsub(\".sv.str.VennInput.txt\", \"\", basename(in_file))\n",
    "    df <- read.table(in_file, header=FALSE, row.names=c(\"A\", \"B\", \"C\", \"AandB\", \"BandC\", \"AandC\", \"AandBandC\"))\n",
    "    df2 <- transpose(df)\n",
    "    colnames(df2) <- rownames(df)\n",
    "    df2$Assembly <- assembly\n",
    "    df2$Region <- \"str\"\n",
    "    rep <- rbind(rep, df2)\n",
    "    df <- read.table(gsub(\"str\", \"segDup\", in_file), header=FALSE, row.names=c(\"A\", \"B\", \"C\", \"AandB\", \"BandC\", \"AandC\", \"AandBandC\"))\n",
    "    df2 <- transpose(df)\n",
    "    colnames(df2) <- rownames(df)\n",
    "    df2$Assembly <- assembly\n",
    "    df2$Region <- \"segDup\"\n",
    "    rep <- rbind(rep, df2)\n",
    "}\n",
    "rep$totalGiaBPass <- rep$A\n",
    "rep$totalGiaB <- rep$B\n",
    "rep$totalOurs <- rep$C\n",
    "rep$sensitivity <- rep$BandC/rep$totalGiaB\n",
    "rep$sensitivityPass <- rep$AandC/rep$totalGiaBPass\n",
    "rep$count <- rep$BandC\n",
    "rep$countPass <- rep$AandC\n",
    "rep <- merge(rep, assemblies, by.x=\"Assembly\", by.y=\"Name\")\n",
    "rep2 <- melt(rep, id.vars=c(\"ID\", \"Region\", \"DataType\"))\n",
    "rep2$value <- as.numeric(rep2$value)\n",
    "rep2$Region <- recode(rep2$Region, segDup=paste(\"segDup n=\", unique(rep$totalGiaBPass[rep$Region==\"segDup\"]), sep=\"\"),\n",
    "                      str=paste(\"str n=\", unique(rep$totalGiaBPass[rep$Region==\"str\"]), sep=\"\"))\n",
    "ggplot(rep2 %>% filter(variable==\"sensitivityPass\"), aes(x=ID, y=value, fill=DataType)) +\n",
    "    geom_bar(stat=\"identity\", position=\"dodge\") +\n",
    "    facet_grid(Region~., scales=\"free\") +\n",
    "    theme_bw() +\n",
    "    theme(axis.text.x = element_text(angle=45, hjust=1, size=7)) +\n",
    "    ggtitle(\"Sensitivity to repetitive element detection\")\n",
    "ggsave(paste(output_dir, \"truthset.repetitive.sensitivity.png\", sep=\"/\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "summary(tsv_final)\n",
    "write.table(tsv_final %>%\n",
    "            select(ID, Name, snv_count, deletions, insertions, sv_count, snv_sensitivity, indel_sensitivity, snv_specificity, indel_specificity, sv_sensitivity, sv_specificity, het_snp_sensitivity, het_indel_sensitivity, hom_snp_sensitivity, hom_indel_sensitivity, DEL, INS, DUP, INV, TRANS, sv_nonrep_count, sv_str_count, sv_segdup_count) %>%\n",
    "            arrange(ID),\n",
    "            file=paste(output_dir, \"stats.tsv\", sep=\"/\"), row.names=FALSE, quote=FALSE, sep=\"\\t\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}