plot gene/transcript read coverage in R

plot_coverage.R

library("wiggleplotr")
library("dplyr")
library("GenomicRanges")
library("GenomicFeatures")
library("EnsDb.Mmusculus.v79")
library("ensembldb")

setwd(<working_dir>)
list.files()

# Load ensembl data base for mouse
EnsDb.Mmusculus.v79

# Plot all transcripts for a specified gene
plotTranscriptsFromEnsembldb(EnsDb.Mmusculus.v79, 
                             gene_names = "Plagl1") 


# Track data format:
# 1. sample_id
# 2. condition
# 3. scaling_factor [1 means no scaling]
# 4. bigWig (paths to bigWig files)
# 5. track_id
# 6. colour_group

# Load sample layout file that normally comes with the project (sample_ids, conditions, etc)
sampleInfo <- read.table(<sample_layout>)

# Create a vector holding paths to bigwig file
bigWigs <- paste(sampleInfo$SampleID, ".bw", sep = "")
# create track data frame
track_data <- data.frame(sample_id = sampleInfo$SampleID,
                         condition = sampleInfo$Group,
                         scaling_factor = rep(1, dim(sampleInfo)[1]),
                         bigWig = bigWigs,
                         track_id=sampleInfo$SampleID,
                         colour_group=sampleInfo$Group) 
track_data

# Plot coverage for all samples and all transcript for the select gene
plotCoverageFromEnsembldb(EnsDb.Mmusculus.v79,
                          "Plagl1",
                          track_data = track_data,
                          heights = c(2,1),
                          fill_palette = getGenotypePalette())

# Limit track data to 3 representative samples
select <- c(<s1>, <s2>, <s3>)
track_data <- track_data[which(track_data$sample_id %in% select),]
track_data

# Plot coverage for select samples and one transcript that appears to be expressed
plotCoverageFromEnsembldb(EnsDb.Mmusculus.v79,
                          "Plagl1",
                          transcript_ids = "ENSMUST00000121646",
                          track_data = track_data,
                          heights = c(2,1),
                          fill_palette = getGenotypePalette())