<jupyter_start><jupyter_text>Load the [StructuralVariantAnnotation](https://bioconductor.org/packages/devel/bioc/vignettes/StructuralVariantAnnotation/inst/doc/vignettes.html) package<jupyter_code>suppressPackageStartupMessages(require(StructuralVariantAnnotation))
suppressPackageStartupMessages(require(dplyr))
suppressPackageStartupMessages(require(ggplot2))<jupyter_output><empty_output><jupyter_text>[SV type inference](https://github.com/PapenfussLab/gridss/blob/7b1fedfed32af9e03ed5c6863d368a821a4c699f/example/simple-event-annotation.RL9)<jupyter_code>infer_svtype <- function(gr)
{
  gr$svtype <-
    ifelse(
      seqnames(gr) != seqnames(partner(gr)),
      "TRA", # Using TRA instead of ITX or BP
      ifelse(
        gr$insLen >= abs(gr$svLen) * 0.7,
        "INS",
        ifelse(
          strand(gr) == strand(partner(gr)),
          "INV",
          ifelse(xor(
            start(gr) < start(partner(gr)), strand(gr) == "-"
          ), "DEL",
          "DUP")
        )
      )
    )
  return(gr)
}
load_bed <- function(bed_file)
{
    bed_regions <- rtracklayer::import(bed_file)
    # set NCBI seqlevels
    seqlevelsStyle(bed_regions) <- "NCBI"
    return(bed_regions)
}
load_bedpe <- function(bedpe_file, filter_regions)
{
  bedpe_gr <- pairs2breakpointgr(rtracklayer::import(bedpe_file))
  bedpe_gr <- filter_regions(bedpe_gr, load_bed(filter_regions), mode='remove')
  return(bedpe_gr)
}
load_vcf <- function(vcf_file, svtype, caller, filter_regions)
{
  
    # Load VCF file
    vcf_gr <-
      VariantAnnotation::readVcf(vcf_file)
    
    # set NCBI seqlevels
    seqlevelsStyle(vcf_gr) <- 'NCBI'
    
    if(caller=='survivor')
    {
        # update info header
        info(header(vcf_gr)) <- rbind(info(header(vcf_gr)),
                                      data.frame(Number = '4', Type = 'String', Description = 'DELLY CT'))
        # update info
        info(vcf_gr) <- cbind(info(vcf_gr), 
                              data.frame(CT=factor(rep('3to5',nrow(info(vcf_gr))), 
                                                   levels=c('5to5', '3to3', '3to5', '5to3'))))
        # TRA
        idx <- which(info(vcf_gr)$SVTYPE=='TRA')
        info(vcf_gr)$CT[idx[seq(1, length(idx), by=2)]]  <- '3to3'
        info(vcf_gr)$CT[idx[seq(2, length(idx), by=2)]]  <- '5to5'
        
        # INV
        idx <- which(info(vcf_gr)[['SVTYPE']]=='INV')
        info(vcf_gr)$CT[idx[seq(1, length(idx), by=2)]]  <- '3to5'
        info(vcf_gr)$CT[idx[seq(2, length(idx), by=2)]]  <- '5to3'
        
        # other SVTYPEs
        idx <- which(!info(vcf_gr)[['SVTYPE']]%in%c('INV','TRA'))
        info(vcf_gr)$CT[idx]  <- '3to5'
        
        # SURVIVOR simSV assigns LowQual to all artificial SVs
        vcf_gr <- vcf_gr[rowRanges(vcf_gr)$FILTER%in%c("LowQual")]
        
        if(svtype == 'INS')
        {
            info(vcf_gr)$END <- end(ranges(rowRanges(vcf_gr)))
        }

    }else{
        # Keep only SVs that passed the filtering (PASS or .)
        vcf_gr <- vcf_gr[rowRanges(vcf_gr)$FILTER%in%c("PASS",".")]
    }

    if (caller == 'lumpy')
    {
      # Read evidence support as a proxy for QUAL
      support <- unlist(info(vcf_gr)$SU)
      fixed(vcf_gr)$QUAL <- support
    } else if (caller == 'delly')
    {
      # Split-read support plus Paired-end read support as a proxy for QUAL
      sr_support <- info(vcf_gr)$SR
      sr_support[is.na(vcf_gr)] <- 0
      fixed(vcf_gr)$QUAL <-
        sr_support + info(vcf_gr)$PE
    }
    
    vcf_gr <- breakpointRanges(vcf_gr)
    vcf_gr <- infer_svtype(vcf_gr)
    
    # Select only one SV type
    vcf_gr <- vcf_gr[which(vcf_gr$svtype == svtype)]

    # Select SVs >= 50 bp
    if(!svtype %in% c('TRA','INS'))
      {
        vcf_gr <- vcf_gr[abs(vcf_gr$svLen) >= 50]
      }
    
    #Filter regions
    vcf_gr <- filter_regions(vcf_gr, load_bed(filter_regions), mode='remove')
    return(vcf_gr)
}
filter_regions <- function(regions_to_filter, ref_regions, mode='remove')
{
    print(length(regions_to_filter))
  if (mode == 'keep')
  {
    result <- regions_to_filter[overlapsAny(regions_to_filter, ref_regions) &
         overlapsAny(partner(regions_to_filter), ref_regions), ]
  } else if (mode == 'remove'){
    result <- regions_to_filter[!(
      overlapsAny(regions_to_filter, ref_regions) |
        overlapsAny(partner(regions_to_filter), ref_regions)
    ), ]
  }
    print(length(result))
    return(result)
}<jupyter_output><empty_output><jupyter_text>Load the truth set<jupyter_code>truth_set <- load_vcf(truth_set_file, svtype, 'survivor', regions_for_filtering)
# truth_set <- truth_set[truth_set$sourceId==svtype]
length(truth_set)<jupyter_output>Warning message:
“info fields with no header: CHR2,END,SVLEN,SVMETHOD,SVTYPE,AF”
Warning message:
“info fields with no header: CHR2,END,SVLEN,SVMETHOD,SVTYPE,AF,CT”
Warning message:
“info fields with no header: CHR2,END,SVLEN,SVMETHOD,SVTYPE,AF,CT”
Warning message:
“info fields with no header: CHR2,END,SVLEN,SVMETHOD,SVTYPE,AF,CT”
Warning message:
“info fields with no header: CHR2,END,SVLEN,SVMETHOD,SVTYPE,AF,CT”
Warning message:
“info fields with no header: CHR2,END,SVLEN,SVMETHOD,SVTYPE,AF,CT”
Warning message:
“info fields with no header: CHR2,END,SVLEN,SVMETHOD,SVTYPE,AF,CT”
Warning message in .breakpointRanges(x, ...):
“Found 1000 duplicate row names (duplicates renamed).”<jupyter_text>Load the SV callsets<jupyter_code>sv_regions <- list()
sv_regions[['CNN']] <- load_bedpe(callsets[['CNN']], regions_for_filtering)
sv_regions[['GRIDSS']] <- load_vcf(callsets[['GRIDSS']], svtype, 'gridss', regions_for_filtering)
sv_regions[['Manta']] <- load_vcf(callsets[['Manta']], svtype, 'manta', regions_for_filtering)
sv_regions[['Lumpy']] <- load_vcf(callsets[['Lumpy']], svtype, 'lumpy', regions_for_filtering)
sv_regions[['DELLY']] <- load_vcf(callsets[['DELLY']], svtype, 'delly', regions_for_filtering)
# exclude callsets with zero calls
sv_regions <- sv_regions[sapply(sv_regions, length)!=0]
# rename NA. column of CNN into svLen
names(mcols(sv_regions[['CNN']]))[names(mcols(sv_regions[['CNN']]))=='NA.'] <- 'svLen'<jupyter_output><empty_output><jupyter_text>Add SV caller name<jupyter_code>for (c in names(sv_regions))
{
  sv_regions[[c]]$caller <- c
}
seqlengths(truth_set) <- c("12" = 2000000, "22" = 2000000)
for (c in names(sv_regions))
{
seqlengths(sv_regions[[c]]) <- c("12" = 2000000, "22" = 2000000)
}
hits <- findBreakpointOverlaps(
      sv_regions[['CNN']],
      sv_regions[['CNN']],
      # read pair based callers make imprecise calls.
      # A margin around the call position is required when matching with the truth set
      maxgap = 100,
      # Since we added a maxgap, we also need to restrict the mismatch between the
      # size of the events. We don't want to match a 100bp deletion with a
      # 5bp duplication. This will happen if we have a 100bp margin but don't also
      # require an approximate size match as well
      sizemargin = 0.25,
      ignore.strand = TRUE,
      # We also don't want to match a 20bp deletion with a 20bp deletion 80bp away
      # by restricting the margin based on the size of the event, we can make sure
      # that simple events actually do overlap
      restrictMarginToSizeMultiple = 0.5
    )
data <- as.data.frame(hits)
getfirst <- function(x){x[1]}
get_best_qual <- function(x){x[which.max(sv_regions[['CNN']]$QUAL[x])]}
aggdata <- aggregate(data$queryHits, by=list(data$subjectHits), FUN=get_best_qual)
sv_regions[['CNN']] <- sv_regions[['CNN']][unique(aggdata$x)]
sv_regions[['CNN']] <- sv_regions[['CNN']][sv_regions[['CNN']]$partner %in% names(sv_regions[['CNN']])]<jupyter_output><empty_output><jupyter_text>Compute overlap<jupyter_code>for (c in names(sv_regions))
{
  if(svtype == 'INS'){
      
  sv_regions[[c]]$truth_matches <-
    countBreakpointOverlaps(
    sv_regions[[c]],
    truth_set,
    # using a smaller margin for insertions, insertion location should be precise
    maxgap = 5,
    # sizemargin cannot be used for insertions
    # sizemargin = 0.25,
    ignore.strand = TRUE,
    restrictMarginToSizeMultiple = 0.5,
    # countOnlyBest cannot be used for insertions
    # countOnlyBest = TRUE
    )
      
  }else{
      
  sv_regions[[c]]$truth_matches <-
    countBreakpointOverlaps(
      sv_regions[[c]],
      truth_set,
      # read pair based callers make imprecise calls.
      # A margin around the call position is required when matching with the truth set
      maxgap = 200,
      # Since we added a maxgap, we also need to restrict the mismatch between the
      # size of the events. We don't want to match a 100bp deletion with a
      # 5bp duplication. This will happen if we have a 100bp margin but don't also
      # require an approximate size match as well
      sizemargin = 0.25,
      ignore.strand = TRUE,
      # We also don't want to match a 20bp deletion with a 20bp deletion 80bp away
      # by restricting the margin based on the size of the event, we can make sure
      # that simple events actually do overlap
      restrictMarginToSizeMultiple = 0.5,
      # Some callers make duplicate calls and will sometimes report a variant multiple
      # times with slightly different bounds. countOnlyBest prevents these being
      # double-counted as multiple true positives.
      countOnlyBest = TRUE
    )
    }
}
    
sv_regions_unlisted <- unlist(GRangesList(sv_regions))<jupyter_output><empty_output><jupyter_text>Plotting Precision and Recall as in [StructuralVariantAnnotation vignette](https://bioconductor.org/packages/devel/bioc/vignettes/StructuralVariantAnnotation/inst/doc/vignettes.html)<jupyter_code>main.title <- paste("SURVIVOR test data", mode, svtype, sep='-')
names(main.title) <- c(sample)

ggplot(
  as.data.frame(sv_regions_unlisted) %>%
    dplyr::select(QUAL, caller, truth_matches) %>%
    dplyr::group_by(caller, QUAL) %>%
    dplyr::summarise(calls = n(),
                     tp = sum(truth_matches > 0)) %>%
    dplyr::group_by(caller) %>%
    dplyr::arrange(dplyr::desc(QUAL)) %>%
    dplyr::mutate(
      cum_tp = cumsum(tp),
      cum_n = cumsum(calls),
      cum_fp = cum_n - cum_tp,
      precision = cum_tp / cum_n,
      recall = cum_tp / length(truth_set)
    )
) +
  aes(x = recall,
      y = precision,
      colour = caller) +
  geom_point() +
  geom_line() +
  scale_y_continuous(labels = scales::percent) +
  scale_x_continuous(
    labels = scales::percent,
    sec.axis = sec_axis( ~ (.) * length(truth_set), name = "true positives")
  ) +
  labs(title = main.title[sample])
filename <- paste("SURVIVOR", mode, svtype, sep='_')
ggsave(file=paste(filename, 'png', sep='.'))
make_percent <- function(x){
  signif(x*100,digits = 4)
}<jupyter_output><empty_output><jupyter_text>Summarize results<jupyter_code>res.df <- as.data.frame(sv_regions_unlisted) %>%
    dplyr::select(caller, truth_matches) %>%
    dplyr::group_by(caller) %>%
    dplyr::summarise(calls = n(),
                   TP = sum(truth_matches > 0)) %>%
    dplyr::group_by(caller) %>%
    dplyr::mutate(
    FP = calls - TP,
    precision = signif(TP / calls, digits = 4),
    recall = signif(TP / length(truth_set), digits = 4)
    )
res.df$F1_score = with(res.df, 2 * (precision * recall) / (precision + recall))
res.df$precision <- make_percent(res.df$precision)
res.df$recall <- make_percent(res.df$recall)
res.df$F1_score <- make_percent(res.df$F1_score)
res.df
write.csv(file=paste(filename, 'csv', sep='.'), res.df, quote=FALSE, row.names=FALSE)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1379220327184128, 'fraction_numerical': 0.014792899408284023, 'mean_word_length': 3.624949698189135, 'pattern_counts': {'":': 0, '<': 101, '<?xml version=': 0, '>': 48, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '237306', 'license': 'Apache-2.0', 'n_tokens_mistral': 4367, 'n_tokens_neox': 4028, 'n_words': 832, 'path': 'notebooks/benchmark/Benchmark.ipynb', 'repo_name': 'renskir/sv-channels', 'size': '11492'}