adding coverage

[nmukherjee]

Is there a way to add a read coverage track above (or below) the transcripts? Either from a bigwig or bam file?

[egustavsson]

I am not sure if this a solution you are after but you can use e.g. plot_grid() to plot coverage above or below the ggtranscript plot. I would just make sure to use the same x-axis ccordinates for both plots. Here is an example function to plot coverage from a bigwig:

  TranscriptCoveragePlot <-
  function(seqnames, start, end, strand, gene_id, gtf, coverage) {

    # Filter long read gtf/gff for gene of interest
    gtf_filtered <- gtf[gtf$gene_id == gene_id]

    # loci used to filter data
    locus_subset <- GRanges(seqnames = seqnames, ranges = IRanges(start = start, end = end), strand = strand)

    # coverage data
    coverage_data <- rtracklayer::import.bw(coverage) %>% subsetByOverlaps(., locus_subset) %>% as.data.frame()

    # Plot transcripts
    exons <- data.frame(gtf_filtered) %>% dplyr::filter(type == "exon")
    introns <- exons %>% to_intron(group_var = "transcript_id")
    CDS <- data.frame(gtf_filtered) %>% dplyr::filter(type == "CDS")

    transcript_plot <-
      exons %>%
      ggplot(aes(
        xstart = start, 
        xend = end, 
        y = transcript_id)) +
      geom_range(fill = "white",
                 height = 0.25) +
      geom_range(data = CDS) +
      geom_intron(
        data = introns,
        arrow.min.intron.length = 500,
        arrow = grid::arrow(ends = "first", length = grid::unit(0.1, "inches"))
        ) +
      labs(y = "Transcript name", x = "") +
      xlim(start(locus_subset), end(locus_subset))

    # coverage data
    coverage_plot <-
      coverage_data %>%
      ggplot(aes(
        xmin = start,
        xmax = end,
        ymin = 0,
        ymax = score
        )) +
      geom_rect(show.legend = F, alpha = 0.8) +
      xlim(start(locus_subset), end(locus_subset))

    # Final plot
    transcript_coverage_plot <-
      plot_grid(
        transcript_plot,
        coverage_plot,
        ncol = 1,
        align = "hv",
        rel_heights = c(1, 3.5),
        axis = "lr"
        )

    return(transcript_coverage_plot)}

[dzhang32]

Thanks @egustavsson for the response! Closing this as resolved, happy to reopen if not.