3
$\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$
  • $\begingroup$ These are very well annotated sequences though. It's the human genome and well known genes. Gene prediction won't give me anything new, all I really want is to identify the extent of this specific duplication. I know it extends beyond the limits of the SMN1/2 genes, I just need to know how far so I can mask one of the two and use the resulting genome as a reference for aligning short reads. $\endgroup$ – terdon Jul 12 at 11:18
  • $\begingroup$ @KamilSJaron I've never used those. Could you post an answer explaining how they would help? $\endgroup$ – terdon Jul 12 at 11:57
  • $\begingroup$ @KamilSJaron here you go: ncbi.nlm.nih.gov/nuccore/NC_000005.9 $\endgroup$ – terdon Jul 12 at 12:23
  • $\begingroup$ Did you try to plot a dotplot (in both axis the sequence and the point is black if both positions do have the same nucleotide) ? This way you can identify similar regions (those that follow a up-left to bottom-right diagonal pattern) and find the positions to mask. $\endgroup$ – llrs Jul 12 at 13:43
  • $\begingroup$ @llrs that's exactly what I'm looking for. I've seen the kind of plot you're talking about, for instance in the output of blast2seq but I don't know how I would use that when comparing the same sequence against itself since all positions are the same by definition. I need to find segmental duplications within the same sequence. I'm now looking at the segmental duplication track the ucsc browser has which was suggested by Chris. $\endgroup$ – terdon Jul 12 at 14:02
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.

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$
  • $\begingroup$ This looks perfect, thanks! I'll give it a go and come back and accept after trying. Might not be before Monday though :) $\endgroup$ – terdon Jul 12 at 15:38
  • $\begingroup$ Yeah. 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. $\endgroup$ – Kamil S Jaron Jul 12 at 16:06
  • 1
    $\begingroup$ UCSC has a segdup track. Multiple tools can do whole-genome homology search such as mashmap or minimap2 (command line here). For close homology, they are 10-100 times faster than older tools. $\endgroup$ – user172818 Jul 12 at 23:33
  • 1
    $\begingroup$ @user172818 yes, the segdup track is what I've been looking at after Chris Rands mentioned it in a (now deleted) comment. That looks promising. As does minimap2, thanks! By the way, do you agree that the best (only?) way to get bwa to align to such troublesome regions is to mask one of the two segmental duplications? Or am I barking up the wrong tree here? $\endgroup$ – terdon Jul 13 at 13:19
  • $\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 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.