Software recommendation: find DNA sequence distribution over entire transcript 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 &lt;- function(x) {
  seqMax &lt;- split(x, factor(x$mgi_symbol))
  seqMax &lt;- lapply(seqMax, function(x) max(x$length)) %&gt;% 
    unlist() %&gt;% 
    as.data.frame()
  seqMax$mgi_symbol &lt;- row.names(seqMax)
  row.names(seqMax) &lt;- NULL
  names(seqMax) &lt;- c("length", "mgi_symbol")
  seqUnique &lt;- merge(seqMax, x)
  seqUnique &lt;- seqUnique[order(seqUnique$mgi_symbol), ]

  return(seqUnique)
}

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

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

  ### Remove unavailable sequences
  sequ &lt;- 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 &lt;- sequ[order(sequ$mgi_symbol), ]
  sequ$length &lt;- sapply(sequ[ ,seqType], nchar)
  sequ &lt;- getMaxLengthSequences(sequ)

  return(sequ)
}

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

  seq &lt;- x %&gt;% 
    tolower() %&gt;% 
    s2c()

  output &lt;- c("XXXX" = count(seq, 4)[c("xxxx", "xxxx")] %&gt;% sum(),
              "XXXX" = count(seq, 4)[c("xxxx", "xxxx")] %&gt;% sum(),
              "XXX" = count(seq, 3)["xxx"] %&gt;% sum())

  return(output)
}

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

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

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

gb &lt;- getBM(attributes = c('ensembl_exon_id', "exon_chrom_start","exon_chrom_end","gene_exon"), 
            filters = "mgi_symbol", 
            values="Grin1", 
            mart = ensembl, 
            bmHeader = TRUE)
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 &lt;- autoplot(EnsDb.Mmusculus.v79, ~ symbol == &quot;Grin1&quot;)

ensembl &lt;- useMart(&quot;ensembl&quot;, dataset = &quot;mmusculus_gene_ensembl&quot;)

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

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

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

# make plot similar to Llopis's answer 
motif_plot &lt;- 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