Answers.Body: <p>I just recently discovered the ampliCan R package from the Valen lab after asking this question. I highly recommend it. It does all of what I was expecting &amp; more.</p>
  
  <p>To install: <a href="https://bioconductor.org/packages/release/bioc/html/amplican.html" rel="nofollow noreferrer">https://bioconductor.org/packages/release/bioc/html/amplican.html</a></p>
  
  <p>To run (&amp; other explanations): <a href="https://rdrr.io/bioc/amplican/f/vignettes/amplicanOverview.Rmd" rel="nofollow noreferrer">https://rdrr.io/bioc/amplican/f/vignettes/amplicanOverview.Rmd</a></p>
  
  <p>This is the original paper &amp; Do not forget to credit the people!
  <a href="https://genome.cshlp.org/content/early/2019/03/08/gr.244293.118.abstract" rel="nofollow noreferrer">https://genome.cshlp.org/content/early/2019/03/08/gr.244293.118.abstract</a></p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-09-18T15:35:52.823Z
Answers.LastActivityDate: 2019-09-18T15:35:52.823Z
Answers.UserDisplayName: francoiskroll
Answers.Body: <p>With 150bp reads, you might be able to blastx all the unique sequences you've got against the target proteins.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-03-27T18:29:38.503Z
Answers.LastActivityDate: 2019-03-27T18:29:38.503Z
Answers.UserDisplayName: swbarnes2
Body: <p>I have Illumina MiSeq paired-end reads from 150bp amplicons mapped to my reference genome (> 1000X coverage). These reads have indels that may or may not induce frameshifts. If the indel induces a frameshift (i.e. indel is not a multiple of 3bp), I want to know if that leads to a stop codon later in the read. Is there any package/method that would allow me to do this?</p>
  
  <p>Few extra details</p>
  
  <ul>
  <li>I might be wrong but I think this is different than gene prediction as I would want the tool to scan each read in the correct frame; not try different frames to detect ORFs etc.</li>
  <li>I expect mosaicism, so I do not want to compute a single assembly/consensus and predict stop codons in it.</li>
  </ul>
  
  <p>Here is how I imagine such a script/tool would work:</p>
  
  <ol>
  <li>Take the position of the first translated codon of the exon I am interested in (from Ensembl for instance?). <em>Eg. gene abc/exon1 starts at chr1:1000.</em></li>
  </ol>
  
  <p>Then for each read:</p>
  
  <ol start="2">
  <li>Look at which position the read is mapped to my reference. <em>Eg. read 1 is mapped to chr1:1050.</em></li>
  <li>Compute the position of the first full codon in the correct frame in my read. <em>Eg. start codon at chr1:1000 but my read starts at chr1:1050, the first full translated codon in my read should start at chr1:1052 (I think?).</em></li>
  <li>Split all my read into 3-mers starting at that position and build something like a 3-mer/position dictionnary. </li>
  <li>Look for stop codons and report the positions.</li>
  </ol>
  
  <p>Does it exist?</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-03-27T12:57:12.353Z
Comments.Text: Could you give us a bit more context? Is this targeted sequencing of a eukaryote? Or whole genome bacterial sequencing? So as input you have a reference genome (FASTA), annotated/predicted genes (GTF/GFF3), and your NGS reads (FASTQ/BAM)? And you'd like to predict stop codons based on the read mappings? Do you have reason to believe the annotated stop codons are incorrect? Are you looking for nonsense mutations?
Comments.UserDisplayName: Daniel Standage
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-03-27T14:08:01.900Z
Comments.Text: @DanielStandage It is targeted sequencing of small (150bp) amplicons from zebrafish genomic DNA. I am sequencing after inducing mutations, i.e. my reads contain various indels. This is why I need to do this for each read in my alignment. I am looking for nonsense mutations, but not only. Eg. there is a 2bp deletion early in the read, I want to know if there will be a stop codon downstream as a result when the RNA gets translated.
Comments.UserDisplayName: francoiskroll
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-03-27T16:10:20.340Z
Comments.Text: Let us [continue this discussion in chat](https://chat.stackexchange.com/rooms/91609/discussion-between-terdon-and-francoiskroll).
Comments.UserDisplayName: terdon
Tags: <phylogenetics><sequence-annotation><illumina>
Title: How to predict stop codons in Illumina reads?
ThreadText: # How to predict stop codons in Illumina reads?
  
  (Asked by: username_0 on 2019-03-27)
  
  I have Illumina MiSeq paired-end reads from 150bp amplicons mapped to my reference genome (> 1000X coverage). These reads have indels that may or may not induce frameshifts. If the indel induces a frameshift (i.e. indel is not a multiple of 3bp), I want to know if that leads to a stop codon later in the read. Is there any package/method that would allow me to do this?
  
  Few extra details
  
  * I might be wrong but I think this is different than gene prediction as I would want the tool to scan each read in the correct frame; not try different frames to detect ORFs etc.
  * I expect mosaicism, so I do not want to compute a single assembly/consensus and predict stop codons in it.
  
  Here is how I imagine such a script/tool would work:
  
  1. Take the position of the first translated codon of the exon I am interested in (from Ensembl for instance?). *Eg. gene abc/exon1 starts at chr1:1000.*
  
  Then for each read:
  
  2. Look at which position the read is mapped to my reference. *Eg. read 1 is mapped to chr1:1050.*
  3. Compute the position of the first full codon in the correct frame in my read. *Eg. start codon at chr1:1000 but my read starts at chr1:1050, the first full translated codon in my read should start at chr1:1052 (I think?).*
  4. Split all my read into 3-mers starting at that position and build something like a 3-mer/position dictionnary.
  5. Look for stop codons and report the positions.
  
  Does it exist?
  
  **Tags:** phylogenetics, sequence-annotation, illumina
  
  > **Comments:**
  > **username_1**: Could you give us a bit more context? Is this targeted sequencing of a eukaryote? Or whole genome bacterial sequencing? So as input you have a reference genome (FASTA), annotated/predicted genes (GTF/GFF3), and your NGS reads (FASTQ/BAM)? And you'd like to predict stop codons based on the read mappings? Do you have reason to believe the annotated stop codons are incorrect? Are you looking for nonsense mutations?  
  > **username_0**: @username_1 It is targeted sequencing of small (150bp) amplicons from zebrafish genomic DNA. I am sequencing after inducing mutations, i.e. my reads contain various indels. This is why I need to do this for each read in my alignment. I am looking for nonsense mutations, but not only. Eg. there is a 2bp deletion early in the read, I want to know if there will be a stop codon downstream as a result when the RNA gets translated.  
  > **username_2**: Let us [continue this discussion in chat](https://chat.stackexchange.com/rooms/91609/discussion-between-username_2-and-username_0).  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_0 on 2019-09-18. Score: 0 | Accepted answer
  
  I just recently discovered the ampliCan R package from the Valen lab after asking this question. I highly recommend it. It does all of what I was expecting & more.
  
  To install: <https://bioconductor.org/packages/release/bioc/html/amplican.html>
  
  To run (& other explanations): <https://rdrr.io/bioc/amplican/f/vignettes/amplicanOverview.Rmd>
  
  This is the original paper & Do not forget to credit the people!
  <https://genome.cshlp.org/content/early/2019/03/08/gr.244293.118.abstract>
  
  
  ---
  
  ### Answer by username_3 on 2019-03-27. Score: 1
  
  With 150bp reads, you might be able to blastx all the unique sequences you've got against the target proteins.
  
  
  ---
  
  