Software recommendation: find DNA sequence distribution over entire transcript

Question

I would like to create a density/histogram of the distribution of a particular DNA sequence over the entire transcript using R and/or command line tools. From here, I would like to use the coordinates of the bins to map the intron-exon diagram below the plot. In this way, I can tell the distribution of various DNA sequences across the entire transcript visually and see differences among alternatively spliced isoforms.

I have already found the number of occurrences of the sequences I am interested in using the seqinr and biostrings packages. I have been using biomaRt to get the exon start and stop sequences as well as the full gene sequence. Between all this, I feel I have all the tools I need but I am missing how to put it all together.

The issue I am having is that biomaRt returns multiple exon sequences with identical start positions but different end positions for a total of 37 exons sequences. Meanwhile, Genatlas and NCBI are returning anywhere from 19 to 21 exons for Grin1, my gene of interest. How can I put all of these tools together to get the plot I need? Is there a better approach?

Finally, if there is an R package or command line tool that already has this function built in, that would be ideal.

Example code for one transcript. This code is actually from an earlier project where I used multiple transcripts but it works fine with one as well. I masked the sequences I'm searching for as "XXXX" at the request of my project manager.

library(magrittr)
library(biomaRt)
library(org.Mm.eg.db)
library(biostrings)
library(seqinr)

### Select only the longest sequences
getMaxLengthSequences <- function(x) {
  seqMax <- split(x, factor(x$mgi_symbol))
  seqMax <- lapply(seqMax, function(x) max(x$length)) %>% 
    unlist() %>% 
    as.data.frame()
  seqMax$mgi_symbol <- row.names(seqMax)
  row.names(seqMax) <- NULL
  names(seqMax) <- c("length", "mgi_symbol")
  seqUnique <- merge(seqMax, x)
  seqUnique <- seqUnique[order(seqUnique$mgi_symbol), ]

  return(seqUnique)
}

### Get sequence for gene of interest
getSequences <- function(genelist, seqType) {
  ensembl <- useMart("ensembl", dataset = "mmusculus_gene_ensembl")

  sequ <- biomaRt::getSequence(id = genelist,
                               type = "mgi_symbol",
                               seqType = seqType,
                               mart = ensembl)

  ### Remove unavailable sequences
  sequ <- sequ[!grepl("Sequence unavailable", sequ[ ,seqType], fixed = TRUE), ]

  ### Find the length of each sequence then pick the max length sequence if there are multiple sequences returned from BioMart
  sequ <- sequ[order(sequ$mgi_symbol), ]
  sequ$length <- sapply(sequ[ ,seqType], nchar)
  sequ <- getMaxLengthSequences(sequ)

  return(sequ)
}

### Get motif counts of interest (substrings are masked)
getmotifs <- function(x) {

  seq <- x %>% 
    tolower() %>% 
    s2c()

  output <- c("XXXX" = count(seq, 4)[c("xxxx", "xxxx")] %>% sum(),
              "XXXX" = count(seq, 4)[c("xxxx", "xxxx")] %>% sum(),
              "XXX" = count(seq, 3)["xxx"] %>% sum())

  return(output)
}

### Using the above functions: Get sequence and then find motifs
seq <- getSequences("Grin1")
targetMotifs <- matrix(sapply(seq[ ,"coding"], getmotifs), nrow = nrow(targets), ncol = 3, byrow = TRUE)

targetMotifTable <- data.frame(Length = seq$length)
rownames(targetMotifTable) <- make.names(seq$mgi_symbol, unique = TRUE)
targetMotifTable <- cbind(targetMotifTable, targetMotifs) %>% as.data.frame()
names(targetMotifTable) <- c("Length", "XXXX", "XXXX", "XXX")

### Use biomaRt to get exons positions
ensembl <- useMart("ensembl", dataset = "mmusculus_gene_ensembl")

gb <- getBM(attributes = c('ensembl_exon_id', "exon_chrom_start","exon_chrom_end","gene_exon"), 
            filters = "mgi_symbol", 
            values="Grin1", 
            mart = ensembl, 
            bmHeader = TRUE)

Answer 1

If I understand the question correctly, you'd like to plot the positions of the matches to you motif along with a gene model that shows the positions of introns and exons for the different transcripts. This can be accomplished fairly easily with ggbio:

library(EnsDb.Mmusculus.v79)
library(ggbio)
library(biomaRt)
library(stringr)
library(dplyr)

# plot gene model with ggbio
gene_model <- autoplot(EnsDb.Mmusculus.v79, ~ symbol == "Grin1")

ensembl <- useMart("ensembl", dataset = "mmusculus_gene_ensembl")

# You need the full length gene sequence to match the gene model
seq <- biomaRt::getSequence(id = "Grin1",
                         type = "mgi_symbol",
                         seqType = "gene_exon_intron",
                         mart = ensembl)

# retrieve position of gene in genome to match up with gene model
gene_start=biomaRt::getBM(attributes=c("start_position"), filters = 'external_gene_name', values = 'Grin1', mart=ensembl)

# use stringR to find positions of exact matches
motif <- "ggcc"
motif_positions <- as.data.frame(str_locate_all(seq$gene_exon_intron, toupper(motif))) %>%
  mutate(start=start+gene_start[[1]]-1)

# make plot similar to Llopis's answer 
motif_plot <- ggplot(motif_positions, aes(x=start)) +
    geom_point(y=1) +
    scale_y_continuous(limits=c(0,1))

# combine plots
tracks(motif= motifs, gene_model=gene_model)

This gives all motif matches for the gene sequence, not just the transcripts. It would be quite straightforward to filter out positions that overlap introns

Answer 2

I assume you have the sequence in seq¹:

# Set the input
seq <- getSequences("GRIN1")
motif <- "ggcc"
# Convert to lower case
seq2 <- tolower(seqinr::s2c(as.character(seq)))
istarts <- seq(from = 1 , to = length(seq2), by = 1)
# Create words of the same size as the motif
oligos <- seq2[istarts]
for (i in 2:nchar(motif)) {
        oligos <- paste(seq2, seq2[istarts + i - 1], sep = "")
}
# Remove exceeding words (words with NA, might need tweaking to motif size
oligos <- oligos[-c(length(oligos):(length(oligos)-nchar(motif)+2))]
# Find motif in the sequence
hits <- grep(motif, oligos, ignore.case = TRUE)
#Prepare the plot
positions <- ifelse(1:length(seq2) %in% hits, 1, 0)
plot(1:length(seq2), positions, col = c("white", "black")[position+1], pch = 19)

The plot is not the density, but you can choose how to plot it now, as you might be interested in windows sizes of 5 nucleotides or 20.

From here you can iterate through multiple sequences. If you want to find several motifs you can take advantage of the oligos object if they are of the same size.

There is the seqPattern package in Bioconductor for more advanced pattern matching. But surprisingly it doesn't have this simple motif search. I think it would be easier to do this with perl and pyhton than in R.

---

1: In the code posted the function getSequence returns an error.