Answers.Body: <p>It could be an idea to fragment the long reads into small sequences, like simulating Illumina reads of 150 bp, and then map these small sequences against the original long reads and extract regions with a high coverage? </p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-02T14:20:35.313Z
Answers.LastActivityDate: 2017-06-02T14:20:35.313Z
Answers.UserDisplayName: Sergio Arredondo
Answers.Body: <p>I've sorted the visualisation out. Here are three alternative representations of repetitive structures for the same sequence:</p>
  <p><a href="https://i.sstatic.net/OIcNm.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/OIcNm.png" alt="Repeat patterns represented as a dotplot" /></a>
  <a href="https://i.sstatic.net/QfV3u.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/QfV3u.png" alt="Repeat patterns represented as a linear &quot;ripple&quot; plot" /></a>
  <a href="https://i.sstatic.net/L158Z.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/L158Z.png" alt="Repeat patterns represented as a semicircular &quot;sunset&quot; plot" /></a></p>
  <p>These were generated using the same R script, callable from the command line:</p>
  <pre><code> <span class="math-container">$ repaver.r -style dotplot -prefix dotplot MAOA_chrX_43645715-43756016.fa
   $</span> repaver.r -style profile -prefix profile MAOA_chrX_43645715-43756016.fa
   $ repaver.r -style semicircular -prefix semi MAOA_chrX_43645715-43756016.fa
  </code></pre>
  <p>More details about this are in the <a href="http://dx.doi.org/10.7490/f1000research.1116032.1" rel="nofollow noreferrer">presentation</a> I gave at Queenstown Research Week, 2018. I also wrote a chapter in a peer-reviewed eBook <a href="https://www.humboldtaustralia.org.au/wp-content/uploads/2018/09/Our-Changing-World-in-the-South-Pacific.pdf" rel="nofollow noreferrer">here</a>.</p>
  <p>This is fast enough that I can run it on the nanopore C. elegans genome in about half an hour, producing these plots for each contig. I don't quite have a method to iterate through this plot and pick out the dominant repeat length at each location, but that's a [relatively] simple extension of what I've already done.</p>
  <p>With a lot of optimisations for speed and memory consumption, I've now been able to run it on the full human genome. It takes a couple of days on my desktop computer (64GB RAM + SSD swap) to categorise the 100bp repeat structure of the assembled T2T/CHM13-v1.0 chromosomes.</p>
  <p>Code <a href="https://gitlab.com/gringer/bioinfscripts/-/blob/master/repaver.r" rel="nofollow noreferrer">here</a>.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-09-04T11:09:00.900Z
Answers.LastActivityDate: 2020-10-04T10:15:41.157Z
Answers.LastEditDate: 2020-10-04T10:15:41.157Z
Answers.UserDisplayName: gringer
Body: <p>After discovering a few difficulties with genome assembly, I've taken an interest in finding and categorising repetitive DNA sequences, such as this one from <em>Nippostrongylus brasiliensis</em> [each base is colour-coded as A: green; C: blue; G: yellow; T: red]:</p>
  
  <p><a href="https://i.sstatic.net/QSk09.png" rel="noreferrer"><img src="https://i.sstatic.net/QSk09.png" alt="Repeat sequence represented in a rectangular fashion"></a>
  <a href="https://i.sstatic.net/RjVLd.png" rel="noreferrer"><img src="https://i.sstatic.net/RjVLd.png" alt="Repeat sequence represented in a circular fashion"></a></p>
  
  <p>[FASTA file associated with this sequence can be found <a href="http://www.gringene.org/data/longRepeat_ch192_r13894.fa" rel="noreferrer">here</a>]</p>
  
  <p>These sequences with large repeat unit sizes are only detectable (and assembleable) using long reads (e.g. PacBio, nanopore) because any subsequence smaller than the unit length will not be able to distinguish between sequencing error and hitting a different location within the repeat structure. I have been tracking these sequences down in a bulk fashion by two methods:</p>
  
  <ol>
  <li>Running an all-vs-all mapping, and looking for sequences that map to themselves lots of times</li>
  <li>Carrying out a compression of the sequence (e.g. bzip2), and finding sequences that have a compression rate that is substantially higher than normal</li>
  </ol>
  
  <p>After I've found a suspicious sequence, I then want to be able to categorise the repeat (e.g. major repeat length, number of tandem repeats, repetitive sequence). This is where I'm getting stuck.</p>
  
  <p>For doing a "look, shiny" demonstration, I currently have a very manual process of getting these sequences into a format that I can visualise. My process is as follows:</p>
  
  <ol>
  <li>Use LAST to produce a dot plot of self-mapping for the mapping</li>
  <li>Visually identify the repetitive region, and extract out the region from the sequence</li>
  <li>Use a combination of <code>fold -w &lt;width&gt;</code> and <code>less -S</code> to visually inspect the sequence with various potential repeat unit widths to find the most likely repeat unit size</li>
  <li>Display the sequence in a rectangular and circular fashion using <a href="https://github.com/gringer/bioinfscripts/blob/master/seqmat.r" rel="noreferrer">my own script</a>, wrapping at the repeat unit length</li>
  </ol>
  
  <p>But that process is by no means feasible when I've got thousands of potential repetitive sequences to fish through.</p>
  
  <p>Is there any better way to do this? Given an arbitrary DNA sequence of length >10kb, how can I (in an automated fashion) find both the location of the repeat region, and also the unit length (bearing in mind that there might be multiple repeat structures, with unit lengths from 30bp to 10kb)?</p>
  
  <p>An example sequence can be found <a href="http://www.gringene.org/data/fd2.fa" rel="noreferrer">here</a>, which has a ~21kb repeat region with ~171bp repeat units about 1/3 of the way into the sequence.</p>
  
  <h2>A Kmer-based Analysis</h2>
  
  <p>I've now seen human sequences with repetitive regions <a href="https://twitter.com/gringene_bio/status/919190635306876928" rel="noreferrer">in excess of 10kb</a> (i.e. out of the range of most linked-read applications). My <a href="https://github.com/gringer/bioinfscripts/blob/master/fastx-rlength.pl" rel="noreferrer">current idea</a> is centred around creating hash tables of short sequences (currently 13-mers) and tracking their location:</p>
  
  <ol>
  <li>Process the sequence, storing the location of each kmer</li>
  <li>For each kmer, find out how many times it appears in the sequence</li>
  <li>For repeated kmers, find out how much of a gap there is between the next time that kmer occurs</li>
  <li>Report the median and modal gap length of repeated kmers, and statistics associated with their frequency in the sequence</li>
  </ol>
  
  <p>Some local repetitive regions may be lost in the statistics with this approach, it's hard to tell if there are multiple repetitive regions within a single sequence, and if the repeat units are themselves slightly repetitive (enough that a kmer is duplicated within a repeat unit), then the algorithm will under-report the repetitiveness (see step 3).</p>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-02T12:13:17.187Z
Comments.Text: Would RepeatMasker help you here? Or does it only find known repetitive elements and regions of low complexity? I've never tried to use it to identify simple repeats.
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-02T12:21:14.977Z
Comments.Text: These repeats have high-complexity repeat units; I haven't found the regions being marked by tantan or dust, and suspect that RepeatMasker will be similar. The web service at least isn't able to find the repeat region in the test sequence that I've linked.
Comments.UserDisplayName: gringer
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-03T13:34:43.197Z
Comments.Text: I know about [this blog post](https://dazzlerblog.wordpress.com/2016/04/01/detecting-and-soft-masking-repeats/) from Gene Myers where he is using strategy of masking repeats in long reads to speed up genome assembly. I am not sure if it is exactly what you want, but I thought it might be useful.
Comments.UserDisplayName: Kamil S Jaron
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-03T13:42:59.153Z
Comments.Text: Thank you. TANmask (at the very end of the post) does look very much like the sort of algorithm I'm looking for. It will at least solve the boundary problem, if not the repeat unit length statistic. Can you write that up as an answer?
Comments.UserDisplayName: gringer
Tags: <nanopore><long-reads><repeat-elements>
Title: Finding the location and unit length of repetitive sequences within a long read
ThreadText: # Finding the location and unit length of repetitive sequences within a long read
  
  (Asked by: username_0 on 2017-06-02)
  
  After discovering a few difficulties with genome assembly, I've taken an interest in finding and categorising repetitive DNA sequences, such as this one from *Nippostrongylus brasiliensis* [each base is colour-coded as A: green; C: blue; G: yellow; T: red]:
  
  [![Repeat sequence represented in a rectangular fashion](https://i.sstatic.net/QSk09.png)](https://i.sstatic.net/QSk09.png)
  [![Repeat sequence represented in a circular fashion](https://i.sstatic.net/RjVLd.png)](https://i.sstatic.net/RjVLd.png)
  
  [FASTA file associated with this sequence can be found [here](http://www.gringene.org/data/longRepeat_ch192_r13894.fa)]
  
  These sequences with large repeat unit sizes are only detectable (and assembleable) using long reads (e.g. PacBio, nanopore) because any subsequence smaller than the unit length will not be able to distinguish between sequencing error and hitting a different location within the repeat structure. I have been tracking these sequences down in a bulk fashion by two methods:
  
  1. Running an all-vs-all mapping, and looking for sequences that map to themselves lots of times
  2. Carrying out a compression of the sequence (e.g. bzip2), and finding sequences that have a compression rate that is substantially higher than normal
  
  After I've found a suspicious sequence, I then want to be able to categorise the repeat (e.g. major repeat length, number of tandem repeats, repetitive sequence). This is where I'm getting stuck.
  
  For doing a "look, shiny" demonstration, I currently have a very manual process of getting these sequences into a format that I can visualise. My process is as follows:
  
  1. Use LAST to produce a dot plot of self-mapping for the mapping
  2. Visually identify the repetitive region, and extract out the region from the sequence
  3. Use a combination of `fold -w <width>` and `less -S` to visually inspect the sequence with various potential repeat unit widths to find the most likely repeat unit size
  4. Display the sequence in a rectangular and circular fashion using [my own script](https://github.com/username_0/bioinfscripts/blob/master/seqmat.r), wrapping at the repeat unit length
  
  But that process is by no means feasible when I've got thousands of potential repetitive sequences to fish through.
  
  Is there any better way to do this? Given an arbitrary DNA sequence of length >10kb, how can I (in an automated fashion) find both the location of the repeat region, and also the unit length (bearing in mind that there might be multiple repeat structures, with unit lengths from 30bp to 10kb)?
  
  An example sequence can be found [here](http://www.gringene.org/data/fd2.fa), which has a ~21kb repeat region with ~171bp repeat units about 1/3 of the way into the sequence.
  
  ## A Kmer-based Analysis
  
  I've now seen human sequences with repetitive regions [in excess of 10kb](https://twitter.com/gringene_bio/status/919190635306876928) (i.e. out of the range of most linked-read applications). My [current idea](https://github.com/username_0/bioinfscripts/blob/master/fastx-rlength.pl) is centred around creating hash tables of short sequences (currently 13-mers) and tracking their location:
  
  1. Process the sequence, storing the location of each kmer
  2. For each kmer, find out how many times it appears in the sequence
  3. For repeated kmers, find out how much of a gap there is between the next time that kmer occurs
  4. Report the median and modal gap length of repeated kmers, and statistics associated with their frequency in the sequence
  
  Some local repetitive regions may be lost in the statistics with this approach, it's hard to tell if there are multiple repetitive regions within a single sequence, and if the repeat units are themselves slightly repetitive (enough that a kmer is duplicated within a repeat unit), then the algorithm will under-report the repetitiveness (see step 3).
  
  **Tags:** nanopore, long-reads, repeat-elements
  
  > **Comments:**
  > **username_1**: Would RepeatMasker help you here? Or does it only find known repetitive elements and regions of low complexity? I've never tried to use it to identify simple repeats.  
  > **username_0**: These repeats have high-complexity repeat units; I haven't found the regions being marked by tantan or dust, and suspect that RepeatMasker will be similar. The web service at least isn't able to find the repeat region in the test sequence that I've linked.  
  > **username_2**: I know about [this blog post](https://dazzlerblog.wordpress.com/2016/04/01/detecting-and-soft-masking-repeats/) from Gene Myers where he is using strategy of masking repeats in long reads to speed up genome assembly. I am not sure if it is exactly what you want, but I thought it might be useful.  
  > **username_0**: Thank you. TANmask (at the very end of the post) does look very much like the sort of algorithm I'm looking for. It will at least solve the boundary problem, if not the repeat unit length statistic. Can you write that up as an answer?  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_0 on 2018-09-04. Score: 3 | Accepted answer
  
  I've sorted the visualisation out. Here are three alternative representations of repetitive structures for the same sequence:
  
  [![Repeat patterns represented as a dotplot](https://i.sstatic.net/OIcNm.png)](https://i.sstatic.net/OIcNm.png)
  [![Repeat patterns represented as a linear "ripple" plot](https://i.sstatic.net/QfV3u.png)](https://i.sstatic.net/QfV3u.png)
  [![Repeat patterns represented as a semicircular "sunset" plot](https://i.sstatic.net/L158Z.png)](https://i.sstatic.net/L158Z.png)
  
  These were generated using the same R script, callable from the command line:
  
  ```
   $ repaver.r -style dotplot -prefix dotplot MAOA_chrX_43645715-43756016.fa
   $ repaver.r -style profile -prefix profile MAOA_chrX_43645715-43756016.fa
   $ repaver.r -style semicircular -prefix semi MAOA_chrX_43645715-43756016.fa
  
  ```
  
  More details about this are in the [presentation](http://dx.doi.org/10.7490/f1000research.1116032.1) I gave at Queenstown Research Week, 2018. I also wrote a chapter in a peer-reviewed eBook [here](https://www.humboldtaustralia.org.au/wp-content/uploads/2018/09/Our-Changing-World-in-the-South-Pacific.pdf).
  
  This is fast enough that I can run it on the nanopore C. elegans genome in about half an hour, producing these plots for each contig. I don't quite have a method to iterate through this plot and pick out the dominant repeat length at each location, but that's a [relatively] simple extension of what I've already done.
  
  With a lot of optimisations for speed and memory consumption, I've now been able to run it on the full human genome. It takes a couple of days on my desktop computer (64GB RAM + SSD swap) to categorise the 100bp repeat structure of the assembled T2T/CHM13-v1.0 chromosomes.
  
  Code [here](https://gitlab.com/username_0/bioinfscripts/-/blob/master/repaver.r).
  
  
  ---
  
  ### Answer by username_3 on 2017-06-02. Score: 2
  
  It could be an idea to fragment the long reads into small sequences, like simulating Illumina reads of 150 bp, and then map these small sequences against the original long reads and extract regions with a high coverage?
  
  
  ---
  
  