Figures

if (!requireNamespace("BiocManager", quietly=TRUE)) {
    install.packages("BiocManager")
}
#BiocManager::install("svaNUMT")
BiocManager::install("PapenfussLab/svaNUMT")
BiocManager::install("plyranges")
BiocManager::install("VariantAnnotation")
BiocManager::install("StructuralVariantAnnotation")
BiocManager::install("TxDb.Hsapiens.UCSC.hg19.knownGene")
BiocManager::install("BSgenome.Hsapiens.UCSC.hg19")
BiocManager::install("S4Vectors")
install.packages("readr")
install.packages("circlize")
install.packages("ggplot2")

library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(dplyr)
library(circlize)
library(svaRetro)
library(svaNUMT)
library(plyranges)
library(readr)
library(ggplot2)
wkdir <- getwd()

Circos plot

BCLAF1 (Fig. 2b)

Retrieve breakpoints supporting the BCLAF1 RT event.

RT_vcf <- readVcf(system.file("extdata", "diploidSV.vcf", package = "svaRetro"))
RT_gr <- StructuralVariantAnnotation::breakpointRanges(RT_vcf, nominalPosition=TRUE)
RT <- rtDetect(RT_gr, TxDb.Hsapiens.UCSC.hg19.knownGene, maxgap=50, minscore=0.3)
rt_chr_prefix <- c(RT$BCLAF1$junctions, RT$BCLAF1$insSite)
seqlevelsStyle(rt_chr_prefix) <- "UCSC"
pairs <- breakpointgr2pairs(rt_chr_prefix)
#pairs

Creating a dataframe with BCLAF1 transcript information.

seqlevelsStyle(TxDb.Hsapiens.UCSC.hg19.knownGene) <- "UCSC"
BCLAF1_txs <- transcriptsByOverlaps(TxDb.Hsapiens.UCSC.hg19.knownGene, GRanges("chr6:136578001-136610989"), columns=c("tx_name", "exon_id"))
BCLAF1_exons <- exonsBy(TxDb.Hsapiens.UCSC.hg19.knownGene, by="tx", use.names=TRUE)[BCLAF1_txs$tx_name]
BCLAF1_df <- bind_rows(mapply(function(x,y) mutate(as.data.frame(x), transcript = y), BCLAF1_exons, names(BCLAF1_exons), SIMPLIFY = FALSE)) %>% 
    mutate(gene = "chr6", exon = exon_rank) %>% 
    dplyr::select(gene, start, end, transcript, exon)

Generating the BCLAF1 circos plot.

#pdf(paste0(wkdir,"/circosRT_BCLAF1_txs.pdf"), width = 5, height = 5)
n = max(tapply(BCLAF1_df$transcript, BCLAF1_df$gene, function(x) length(unique(x))))
circos.genomicInitialize(bind_rows(BCLAF1_df, 
                                   tibble(gene = c("chr5", "chr6"), 
                                          start = c(110270000,136577501),
                                          end = c(110290000,136611489), 
                                          transcript = c("q22.1","q23.3"),
                                          exon = c(0,0))), 
                         tickLabelsStartFromZero=FALSE)

circos.genomicTrack(BCLAF1_df, ylim = c(0.2, n + 0.2), 
    panel.fun = function(region, value, ...) {
        all_tx = unique(value$transcript)
        for(i in seq_along(all_tx)) {
            l = value$transcript == all_tx[i]
            # for each transcript
            current_tx_start = min(region[l, 1])
            current_tx_end = max(region[l, 2])
            circos.lines(c(current_tx_start, current_tx_end), 
                c(n - i + 1, n - i + 1), col = "#324e64") #txs lines
            
            circos.genomicRect(region[l, , drop = FALSE], ytop = n - i + 1 + 0.3, 
                ybottom = n - i + 1 - 0.3, col = "#f3d197", border = NA) #txs exons
        }
}, bg.col = c("#aacfdd", "#f0855f"), bg.border = NA, track.height = 0.2) #chrs
circos.genomicLink(as.data.frame(S4Vectors::first(pairs)), 
                   as.data.frame(S4Vectors::second(pairs)),
                   col=c(rep("#2b4b9b", 4), "#e52822", "#2b4b9b", "#e52822")) #bps

circos.clear()
dev.off()

## null device 
##           1

Simulated NUMT (Fig. 3b)

numtS <- read_table(paste0(wkdir, "/../vsDinumt/numtS.txt"), col_names = FALSE)

## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   X1 = col_double(),
##   X2 = col_character(),
##   X3 = col_double(),
##   X4 = col_double(),
##   X5 = col_character(),
##   X6 = col_double(),
##   X7 = col_character()
## )

colnames(numtS) <- c("bin", "seqnames", "start", "end", "name", "score", "strand")
numtS <- GRanges(numtS) %>% `seqlevelsStyle<-`("NCBI")

library(BSgenome.Hsapiens.UCSC.hg19)

## Loading required package: BSgenome

genomeMT <- BSgenome.Hsapiens.UCSC.hg19$chrMT

#library(circlize)
#proj_path="/Users/dong.rn/Documents/svaPkgs/simReads/simNUMT/"
all_numt.gr <- breakpointRanges(readVcf(paste0(wkdir, "/../simReads/simNUMT/chr1_numt_pe_HS25.sv.vcf")), nominalPosition=TRUE)
all_numt <- numtDetect(all_numt.gr, numtS, genomeMT, 
                       max_ins_dist = 10000, maxgap_numtS=100, 
                       min_len=50, min.Align=0.8)

numt_chr_prefix <- c(all_numt$MT$NU$`1`[[7]], all_numt$MT$MT$`1`[[7]])
GenomeInfoDb::seqlevelsStyle(numt_chr_prefix) <- "UCSC"
numt_pairs <- breakpointgr2pairs(numt_chr_prefix)
numt_pairs

cytoband_df = system.file(package = "circlize", "extdata", "cytoBand.txt")
#start and end must be specified as 'numeric' instead of 'integer'!
cytoband_df = read.table(cytoband_df, colClasses = c("character", "numeric",
    "numeric", "character", "character"), sep = "\t")
colnames(cytoband_df) <- c('seqnames', 'start', 'end', 'region', 'name')

#pdf(paste0(wkdir,"/circosNUMT_sim7.pdf"), width = 5, height = 5)
circos.initializeWithIdeogram(filter(cytoband_df, seqnames=="chr1") %>% 
                                  rbind(data.frame(seqnames="chrM", start=1, end=16571, 
                                                   region=NA, name=NA)),
                              sector.width = c(0.5, 0.5))
circos.genomicLink(as.data.frame(S4Vectors::first(numt_pairs)), as.data.frame(S4Vectors::second(numt_pairs)))

circos.clear()
dev.off()

## null device 
##           1

Bar plot

The statistics below are acquired from analysis done in sim_reads_rt.Rmd.

sim_reads_results <- tibble(category=rep(c('detected', 'undetected bps', 
                      'unmapped insSeqs', 'misaligned reads',
                      'detected', 'undetected bps', 
                      'single insSite', 'unmapped insSeqs'), each = 2),
       count=c(470, 443, 23, 45, 5, 11, 2, 1, #RT
               352, 238, 14, 75, 88, 118, 46, 69), #NUMT
       caller=rep(c('GRIDSS', 'Manta'), 8),
       group=rep(c('RT', 'NUMT'), each=8))

stacked bars (Fig. 4)

sim_reads_results %>% tidyr::unite("set", group:caller, remove = FALSE) %>% 
    ggplot(., aes(set, count, fill=reorder(category, count))) +
    geom_bar(stat="identity", width=0.5) +
    #scale_fill_manual(values = c('#58595b','#7d489c','#ece953','#49a942','#2372b9')) + 
    #scale_fill_manual(values = c('#8dd3c7','#ffffb3','#bebada','#fb8072','#80b1d3')) +
    scale_fill_brewer(palette="Pastel1") +
    geom_text(aes(label=count), color="black", size=3, position = position_stack(vjust = 0.5)) +
    coord_flip() +
    labs(x="", y="", fill="category") +
    theme_minimal()

#ggsave(filename = paste0(wkdir, "/simRTnNUMT.pdf"), width = 9, height = 3)
#ggsave("~/Documents/PhD project/SVEnsemble/output/seminar/simRTnNUMT.png", width = 9, height = 4.5)

barplots

simulated RTs:

sim_reads_results %>% filter(group=="RT") %>% 
    mutate(category=stringr::str_wrap(category, width = 11)) %>% 
    ggplot(., aes(category, count, fill=caller)) + 
    geom_bar(stat="identity", position=position_dodge()) +
    geom_text(aes(label=count), vjust=-0.1, color="black",
            position = position_dodge(0.9), size=3.5) +
    scale_fill_brewer(palette="Pastel1") + 
    xlab("category of simulated RTs") +
    #coord_flip() +
    ylim(0, 500) +
    #facet_grid(.~group) +
    theme_minimal()

#ggsave(filename = paste0(wkdir, "/simRT.pdf"), width = 5, height = 5)

simulated NUMTs:

sim_reads_results %>% filter(group=="NUMT") %>% 
    mutate(category=stringr::str_wrap(category, width = 11)) %>% 
    ggplot(., aes(category, count, fill=caller)) + 
    geom_bar(stat="identity", position=position_dodge()) +
    geom_text(aes(label=count), vjust=-0.1, color="black",
            position = position_dodge(0.9), size=3.5) +
    scale_fill_brewer(palette="Pastel1") + 
    xlab("category of simulated NUMTs") +
    #coord_flip() +
    ylim(0, 500) +
    #facet_grid(.~group) +
    theme_minimal()

#ggsave(filename = paste0(wkdir, "/simNUMT.pdf"), width = 5, height = 5)