5
$\begingroup$

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:

UCSC genome browser BLAT results

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?

$\endgroup$
0
2
$\begingroup$

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

GeneBlocks plot

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)

DNA Features Viewer plot

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

$\endgroup$
2
$\begingroup$

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

enter image description here

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.

$\endgroup$
1
  • $\begingroup$ @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 :( $\endgroup$ – terdon Jul 13 '19 at 13:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.