text: user1488: 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)
  
  
  user48: I assume you have the sequence in seq1:
  # 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.
  
  user581: 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
  
original_text: 1488: <p>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. </p>
  
  <p>I have already found the number of occurrences of the sequences I am interested in using the <code>seqinr</code> and <code>biostrings</code> packages. I have been using <code>biomaRt</code> 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. </p>
  
  <p>The issue I am having is that <code>biomaRt</code> returns multiple exon sequences with identical start positions but different end positions for a total of 37 exons sequences. Meanwhile, <a href="http://genatlas.medecine.univ-paris5.fr/fiche.php?symbol=grin1" rel="nofollow noreferrer">Genatlas</a> and <a href="https://www.ncbi.nlm.nih.gov/IEB/Research/Acembly/av.cgi?db=mouse&amp;c=Gene&amp;l=Grin1" rel="nofollow noreferrer">NCBI</a> 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?</p>
  
  <p>Finally, if there is an R package or command line tool that already has this function built in, that would be ideal.</p>
  
  <p>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.</p>
  
  <pre><code>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)
  </code></pre>
   48: <p>I assume you have the sequence in seq<sup>1</sup>:</p>
  
  <pre><code># Set the input
  seq &lt;- getSequences("GRIN1")
  motif &lt;- "ggcc"
  # Convert to lower case
  seq2 &lt;- tolower(seqinr::s2c(as.character(seq)))
  istarts &lt;- seq(from = 1 , to = length(seq2), by = 1)
  # Create words of the same size as the motif
  oligos &lt;- seq2[istarts]
  for (i in 2:nchar(motif)) {
          oligos &lt;- paste(seq2, seq2[istarts + i - 1], sep = "")
  }
  # Remove exceeding words (words with NA, might need tweaking to motif size
  oligos &lt;- oligos[-c(length(oligos):(length(oligos)-nchar(motif)+2))]
  # Find motif in the sequence
  hits &lt;- grep(motif, oligos, ignore.case = TRUE)
  #Prepare the plot
  positions &lt;- ifelse(1:length(seq2) %in% hits, 1, 0)
  plot(1:length(seq2), positions, col = c("white", "black")[position+1], pch = 19)
  </code></pre>
  
  <p>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. 
  <a href="https://i.stack.imgur.com/6CWmu.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/6CWmu.png" alt="positions where the motif is found in a sequence"></a></p>
  
  <p>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.</p>
  
  <p>There is the <a href="http://bioconductor.org/packages/seqPattern/" rel="nofollow noreferrer">seqPattern</a> 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.</p>
  
  <hr>
  
  <p>1: In the code posted the function getSequence returns an error.</p>
   581: <p>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:</p>
  <pre><code>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)
  </code></pre>
  <p><a href="https://i.stack.imgur.com/JHJ8x.png" rel="noreferrer"><img src="https://i.stack.imgur.com/JHJ8x.png" alt="Gene Model with Motif matchs" /></a></p>
  <p>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</p>
  