Answers.Body: <p>If you know approximately where are the duplicated regions, you can cut them out from the genome and align the two block using tools like <a href="https://github.com/garviz/MUMmer" rel="nofollow noreferrer">MUMmer</a> or <a href="http://darlinglab.org/mauve/" rel="nofollow noreferrer">Mauve</a>. Although Mauve is not scalable solution (You can not scan genomes for duplicates with it) it's very handy for aligning small(ish) blocks. It can also align multiple sequences if you would like to check some tri/multiplications.</p>
  <p>I downloaded <a href="https://www.ncbi.nlm.nih.gov/nuccore/NC_000005.9" rel="nofollow noreferrer">the sequence of ch5</a> you linked and extracted big regions around the genes (300bp blocks). The first block is from 70100000 - 70400000 (with SMN1 sequence at 115768 - 150839) and the second from 69200000 - 69500000 (with SMN2 sequence at 140350 - 175383).</p>
  <p>I used progressive mauve to align the two (click on File -&gt; align with ProgressiveMauve -&gt; add Sequence ... -&gt; Align).</p>
  <p><a href="https://i.sstatic.net/RAT8G.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/RAT8G.png" alt="enter image description here" /></a></p>
  <p>I found that the duplicated region is almost 200k long and maps between 44690 - 225520 of the first block and 19760 - 200940 of the second block, which we still need translate to the original genome coordinates, but I believe I can leave this to you :-).</p>
  <p>P.S. you can also somehow upload annotations to Mauve, but I was lazy, so I &quot;hand drawn&quot; the arrow where I knew the genes are.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-07-13T13:47:35.447Z
Answers.Text: @KamilSJaron I can't get this to work. It keeps running out of memory even on my laptop which has 32G, of which 20 were available. I even tried aligning 2 files of only 100bp each and it crashed. I managed to get it to run on a server with 41G of RAM available, but I can't get the graphical output that way. I see you ran it on macOS, so it seems like a bug of the Linux version :(
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-07-12T15:21:26.860Z
Answers.LastActivityDate: 2021-07-05T04:00:47.227Z
Answers.LastEditDate: 2021-07-05T04:00:47.227Z
Answers.UserDisplayName: Kamil S Jaron
Answers.Body: <p><a href="https://github.com/Edinburgh-Genome-Foundry/geneblocks" rel="nofollow noreferrer">GeneBlocks</a> is a Python library for comparing DNA sequences. It can be used to find common blocks (identical regions) in a group of DNA sequences. Using a modified version of your example (I duplicated a small region):</p>
  <pre class="lang-py prettyprint-override"><code>import geneblocks
  
  sequence = geneblocks.load_record('sequence.gb', name='seq')
  sequence2 = geneblocks.load_record('sequence.gb', name='seq2')
  blocks = geneblocks.CommonBlocks.from_sequences([sequence, sequence2])
  
  blocks.plot_common_blocks()
  </code></pre>
  <p><a href="https://i.sstatic.net/WtTtT.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/WtTtT.png" alt="GeneBlocks plot" /></a></p>
  <p>Note that it wasn't designed to plot only one sequence, so as a workaround the sequence is given twice, with different names. As an alternative, you can use one sequence and plot with <a href="https://github.com/Edinburgh-Genome-Foundry/DnaFeaturesViewer" rel="nofollow noreferrer">DNA Features Viewer</a>:</p>
  <pre class="lang-py prettyprint-override"><code>import dna_features_viewer
  
  blocks = geneblocks.CommonBlocks.from_sequences([sequence])
  seq_to_plot = blocks.sequences_with_annotated_blocks()[&quot;seq&quot;]
  graphic_record = dna_features_viewer.BiopythonTranslator().translate_record(seq_to_plot)
  ax, _ = graphic_record.plot(figure_width=10, strand_in_label_threshold=7)
  </code></pre>
  <p><a href="https://i.sstatic.net/CnSv3.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/CnSv3.png" alt="DNA Features Viewer plot" /></a></p>
  <p>See the <a href="https://edinburgh-genome-foundry.github.io/Geneblocks/" rel="nofollow noreferrer">documentation</a> for more details. A web interface is also provided at EGF CUBA: <a href="https://cuba.genomefoundry.org/find-common-blocks" rel="nofollow noreferrer">find common blocks</a>.</p>
  <p><em>Disclaimer: I'm the current maintainer of GeneBlocks</em></p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-11-26T18:26:48.493Z
Answers.LastActivityDate: 2020-11-26T18:45:15.327Z
Answers.LastEditDate: 2020-11-26T18:45:15.327Z
Answers.UserDisplayName: Peter
Body: <p>I am facing an issue when trying to align short reads against a region in human chr5. The two Sensory Motor Neuron genes, (SMN1 and SMN2) are almost 100% identical and this causes the aligner to fail to align reads correctly since each read matches perfectly to two separate regions of the reference, leading to very low alignment scores.</p>
  <p>If I BLAT the <a href="https://www.ncbi.nlm.nih.gov/nuccore/NG_008691.1?report=fasta" rel="nofollow noreferrer">SMN1 gene's full genomic sequence</a> against the human genome (hg19) on UCSC, I get (only showing the 1st two lines):</p>
  <pre><code>QUERY       SCORE START   END QSIZE IDENTITY  CHROM  STRAND  START       END   SPAN
  --------------------------------------------------------------------------------------------------------------
  NG_008691.1 35072     1 35072 35072   100.0%  chr5   +    70215768  70250839  35072
  NG_008691.1 34279     1 35072 35072    99.8%  chr5   +    69340350  69375383  35034
  </code></pre>
  <p>And <a href="http://genome.ucsc.edu/cgi-bin/hgTracks?db=hg19&amp;lastVirtModeType=default&amp;lastVirtModeExtraState=&amp;virtModeType=default&amp;virtMode=0&amp;nonVirtPosition=&amp;position=chr5%3A69311399%2D70350568&amp;hgsid=739135289_3EijLmAz1FaCi8vZojAHjpMDEwU0" rel="nofollow noreferrer">on the browser</a>:</p>
  <p><a href="https://i.sstatic.net/4osMP.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/4osMP.png" alt="UCSC genome browser BLAT results" /></a></p>
  <p>To get around this, I want to mask the SMN2 gene. However, the region of similarity extends beyond the annotated boundaries of the gene, so I would like to identify the full extent of this duplication.</p>
  <p>So, given a large nucleotide sequence in fasta format, how can I compare it to itself and identify identical or near identical regions?</p>
  
Tags: <fasta><genome><sequence-analysis><sequence-homology><duplications>
Title: How can I locate duplicated regions in a sequence?
ThreadText: # How can I locate duplicated regions in a sequence?
  
  (Asked by: username_0 on 2019-07-12)
  
  I am facing an issue when trying to align short reads against a region in human chr5. The two Sensory Motor Neuron genes, (SMN1 and SMN2) are almost 100% identical and this causes the aligner to fail to align reads correctly since each read matches perfectly to two separate regions of the reference, leading to very low alignment scores.
  
  If I BLAT the [SMN1 gene's full genomic sequence](https://www.ncbi.nlm.nih.gov/nuccore/NG_008691.1?report=fasta) against the human genome (hg19) on UCSC, I get (only showing the 1st two lines):
  
  ```
  QUERY       SCORE START   END QSIZE IDENTITY  CHROM  STRAND  START       END   SPAN
  --------------------------------------------------------------------------------------------------------------
  NG_008691.1 35072     1 35072 35072   100.0%  chr5   +    70215768  70250839  35072
  NG_008691.1 34279     1 35072 35072    99.8%  chr5   +    69340350  69375383  35034
  
  ```
  
  And [on the browser](http://genome.ucsc.edu/cgi-bin/hgTracks?db=hg19&lastVirtModeType=default&lastVirtModeExtraState=&virtModeType=default&virtMode=0&nonVirtPosition=&position=chr5%3A69311399%2D70350568&hgsid=739135289_3EijLmAz1FaCi8vZojAHjpMDEwU0):
  
  [![UCSC genome browser BLAT results](https://i.sstatic.net/4osMP.png)](https://i.sstatic.net/4osMP.png)
  
  To get around this, I want to mask the SMN2 gene. However, the region of similarity extends beyond the annotated boundaries of the gene, so I would like to identify the full extent of this duplication.
  
  So, given a large nucleotide sequence in fasta format, how can I compare it to itself and identify identical or near identical regions?
  
  **Tags:** fasta, genome, sequence-analysis, sequence-homology, duplications
  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_1 on 2019-07-12. Score: 2
  
  If you know approximately where are the duplicated regions, you can cut them out from the genome and align the two block using tools like [MUMmer](https://github.com/garviz/MUMmer) or [Mauve](http://darlinglab.org/mauve/). Although Mauve is not scalable solution (You can not scan genomes for duplicates with it) it's very handy for aligning small(ish) blocks. It can also align multiple sequences if you would like to check some tri/multiplications.
  
  I downloaded [the sequence of ch5](https://www.ncbi.nlm.nih.gov/nuccore/NC_000005.9) you linked and extracted big regions around the genes (300bp blocks). The first block is from 70100000 - 70400000 (with SMN1 sequence at 115768 - 150839) and the second from 69200000 - 69500000 (with SMN2 sequence at 140350 - 175383).
  
  I used progressive mauve to align the two (click on File -> align with ProgressiveMauve -> add Sequence ... -> Align).
  
  [![enter image description here](https://i.sstatic.net/RAT8G.png)](https://i.sstatic.net/RAT8G.png)
  
  I found that the duplicated region is almost 200k long and maps between 44690 - 225520 of the first block and 19760 - 200940 of the second block, which we still need translate to the original genome coordinates, but I believe I can leave this to you :-).
  
  P.S. you can also somehow upload annotations to Mauve, but I was lazy, so I "hand drawn" the arrow where I knew the genes are.
  
  > **Comments:**
  > **username_0**: @username_1 I can't get this to work. It keeps running out of memory even on my laptop which has 32G, of which 20 were available. I even tried aligning 2 files of only 100bp each and it crashed. I managed to get it to run on a server with 41G of RAM available, but I can't get the graphical output that way. I see you ran it on macOS, so it seems like a bug of the Linux version :(  
  
  ---
  
  ### Answer by username_2 on 2020-11-26. Score: 2
  
  [GeneBlocks](https://github.com/Edinburgh-Genome-Foundry/geneblocks) is a Python library for comparing DNA sequences. It can be used to find common blocks (identical regions) in a group of DNA sequences. Using a modified version of your example (I duplicated a small region):
  
  ```py
  import geneblocks
  
  sequence = geneblocks.load_record('sequence.gb', name='seq')
  sequence2 = geneblocks.load_record('sequence.gb', name='seq2')
  blocks = geneblocks.CommonBlocks.from_sequences([sequence, sequence2])
  
  blocks.plot_common_blocks()
  
  ```
  
  [![GeneBlocks plot](https://i.sstatic.net/WtTtT.png)](https://i.sstatic.net/WtTtT.png)
  
  Note that it wasn't designed to plot only one sequence, so as a workaround the sequence is given twice, with different names. As an alternative, you can use one sequence and plot with [DNA Features Viewer](https://github.com/Edinburgh-Genome-Foundry/DnaFeaturesViewer):
  
  ```py
  import dna_features_viewer
  
  blocks = geneblocks.CommonBlocks.from_sequences([sequence])
  seq_to_plot = blocks.sequences_with_annotated_blocks()["seq"]
  graphic_record = dna_features_viewer.BiopythonTranslator().translate_record(seq_to_plot)
  ax, _ = graphic_record.plot(figure_width=10, strand_in_label_threshold=7)
  
  ```
  
  [![DNA Features Viewer plot](https://i.sstatic.net/CnSv3.png)](https://i.sstatic.net/CnSv3.png)
  
  See the [documentation](https://edinburgh-genome-foundry.github.io/Geneblocks/) for more details. A web interface is also provided at EGF CUBA: [find common blocks](https://cuba.genomefoundry.org/find-common-blocks).
  
  *Disclaimer: I'm the current maintainer of GeneBlocks*
  
  
  ---
  
  