Answers.Body: <p>Short answer for short reads: you can't. Not reliably, at least.</p>
  <p>There's a known error mode with aligning short reads to reference genomes where it's quite good at mis-calling reads within repetitive regions. See, for example, <a href="https://doi.org/10.1093/gigascience/giz104" rel="nofollow noreferrer">this paper</a>, which compares Illumina vs Nanopore for detecting bacterial strains.</p>
  <p>If you know that a region is likely to contain a segmental duplication and/or deletion (e.g. due to unexpected changes in read coverage), then it'd be a good idea to double-check any discovered variants (or unexpected non-variants) using a different technology that does not have the same coverage issues. I recommend a long-read sequencing technology that can handle regions at least as large as the segmental duplication of interest (e.g. PacBio up to ~15kb, Nanopore for anything longer), so that you have a high chance of mapping reads uniquely.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-07-05T08:20:12.827Z
Answers.Text: I agree that a longer read technique would help but unfortunately that isn't an option for me as I am processing the data generated by others. I need to improve the results of using short reads and don't have the option of using long. But I don't see why you say "you can't". It turns out I most certainly can, and the approach described in the paper you linked to in the question (thanks, I should have included that) works. The complication is dealing with different repeat numbers and ploidy, but the general approach certainly seems to work so far and is giving more accurate results.
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-07-05T04:46:53.230Z
Answers.LastActivityDate: 2021-07-05T04:46:53.230Z
Answers.UserDisplayName: gringer
Body: <p>This is a follow-up of <a href="https://bioinformatics.stackexchange.com/q/8975/298">my other question</a>. I have been having trouble calling variants in the human SMN1 and SMN2 genes, because the human genome has a large segmental duplication there and these two genes (and their surrounding sequences) are essentially identical :</p>
  <p><a href="https://i.sstatic.net/0UiCX.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/0UiCX.png" alt="UCSC browser screenshot" /></a>
  <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"><sub>Click here to see the original image on the UCSC browser.</sub></a></p>
  <p>Since variants in the SMN1 gene have been linked to disease, it is important for me to be able to call them in my NGS data. However, because of this segmental duplication, aligning reads to this region correctly is very hard since the presence of the duplication makes each read align almost perfectly to separate locations in the genome, leading to low alignment scores. The low scores will, in turn, block variant callers from calling in those regions.</p>
  <p>I thought that the only real solution for this would be to mask the SMN2 gene and its associated duplication in my reference genome which would allow the aligner to align all reads to the region of SMN1. I would then miss any real variants in SMN2, but since the primary gene of interest for  spinal muscular atrophy is SMN1 and not SMN2, I would be OK with that.</p>
  <p>My questions are:</p>
  <ol>
  <li><p>Is this a reasonable approach? I was told that <a href="https://bioinformatics.stackexchange.com/questions/8975/how-can-i-locate-duplicated-regions-in-a-sequence#comment12265_8976">&quot;generally, masking segdups is not recommended&quot;</a> by someone who most certainly knows what he's talking about. So is this not a good approach for my stated goal?</p>
  </li>
  <li><p>If not, then what else can I do that will let me call variants in the SMN1 gene?</p>
  </li>
  </ol>
  <p>I would like something similar to the approach described <a href="https://genomebiology.biomedcentral.com/articles/10.1186/s13059-019-1707-2" rel="nofollow noreferrer">here</a>: mask one of the regions to rescue mutations in the other.</p>
  
Tags: <phylogenetics><variant-calling><reference-genome>
Title: How should I deal with segmental duplications when aligning NGS reads to a reference genome?
ThreadText: # How should I deal with segmental duplications when aligning NGS reads to a reference genome?
  
  (Asked by: username_0 on 2019-07-15)
  
  This is a follow-up of [my other question](https://bioinformatics.stackexchange.com/q/8975/298). I have been having trouble calling variants in the human SMN1 and SMN2 genes, because the human genome has a large segmental duplication there and these two genes (and their surrounding sequences) are essentially identical :
  
  [![UCSC browser screenshot](https://i.sstatic.net/0UiCX.png)](https://i.sstatic.net/0UiCX.png)
  [Click here to see the original image on the UCSC 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)
  
  Since variants in the SMN1 gene have been linked to disease, it is important for me to be able to call them in my NGS data. However, because of this segmental duplication, aligning reads to this region correctly is very hard since the presence of the duplication makes each read align almost perfectly to separate locations in the genome, leading to low alignment scores. The low scores will, in turn, block variant callers from calling in those regions.
  
  I thought that the only real solution for this would be to mask the SMN2 gene and its associated duplication in my reference genome which would allow the aligner to align all reads to the region of SMN1. I would then miss any real variants in SMN2, but since the primary gene of interest for spinal muscular atrophy is SMN1 and not SMN2, I would be OK with that.
  
  My questions are:
  
  1. Is this a reasonable approach? I was told that ["generally, masking segdups is not recommended"](https://bioinformatics.stackexchange.com/questions/8975/how-can-i-locate-duplicated-regions-in-a-sequence#comment12265_8976) by someone who most certainly knows what he's talking about. So is this not a good approach for my stated goal?
  2. If not, then what else can I do that will let me call variants in the SMN1 gene?
  
  I would like something similar to the approach described [here](https://genomebiology.biomedcentral.com/articles/10.1186/s13059-019-1707-2): mask one of the regions to rescue mutations in the other.
  
  **Tags:** phylogenetics, variant-calling, reference-genome
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2021-07-05. Score: 1
  
  Short answer for short reads: you can't. Not reliably, at least.
  
  There's a known error mode with aligning short reads to reference genomes where it's quite good at mis-calling reads within repetitive regions. See, for example, [this paper](https://doi.org/10.1093/gigascience/giz104), which compares Illumina vs Nanopore for detecting bacterial strains.
  
  If you know that a region is likely to contain a segmental duplication and/or deletion (e.g. due to unexpected changes in read coverage), then it'd be a good idea to double-check any discovered variants (or unexpected non-variants) using a different technology that does not have the same coverage issues. I recommend a long-read sequencing technology that can handle regions at least as large as the segmental duplication of interest (e.g. PacBio up to ~15kb, Nanopore for anything longer), so that you have a high chance of mapping reads uniquely.
  
  > **Comments:**
  > **username_0**: I agree that a longer read technique would help but unfortunately that isn't an option for me as I am processing the data generated by others. I need to improve the results of using short reads and don't have the option of using long. But I don't see why you say "you can't". It turns out I most certainly can, and the approach described in the paper you linked to in the question (thanks, I should have included that) works. The complication is dealing with different repeat numbers and ploidy, but the general approach certainly seems to work so far and is giving more accurate results.  
  
  ---
  
  