How can I locate duplicated regions in a sequence?

Question

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 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:

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?

Answer 1

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):

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()

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:

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)

See the documentation for more details. A web interface is also provided at EGF CUBA: find common blocks.

Disclaimer: I'm the current maintainer of GeneBlocks

Answer 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 or 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 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).

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.