Answers.Body: <p>Getting SNPs from whole-genome data is usually done via mapping sequencing reads on the reference (bwa-mem or bowtie are two popular mappers) and then using one of variant calling tools (like freebayes, samtools, or GATK). However, I would also look up if there is something virus-specific.</p>
  <p>According to <a href="https://www.biostars.org/p/93673/#93676" rel="nofollow noreferrer">this answer</a>, you could get variant calls of a specific region using samtools if you give it a bed file that contains the region of interest (three columns <code>chr    start    stop</code> format). For example using human reference</p>
  <pre><code>samtools mpileup -uDl regions.bed -f hg19.fa file1.sorted.bam file2.sorted.bam | bcftools view -bvcg - &gt; RAL_samtools.raw.bcf
  </code></pre>
  <p>You could also specify a bunch of files, not just two and it would do the variant calling on them all.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-05-03T13:15:20.363Z
Answers.LastActivityDate: 2021-05-03T13:15:20.363Z
Answers.UserDisplayName: Kamil S Jaron
Answers.Body: <p><em>Answer from @atpoint [transferred from comment]:</em></p>
  <p>Alternatively, you can use <code>samtools view</code> to only retain reads that overlap the region you are interested in and then feed these into the mpileup. IIRC there was something with that <code>-l</code> option in mpileup that caused trouble but I do not remember what it was. I think it still called all variants based on all reads and then simply filtered for the target regions in <code>-l</code>. The <code>view</code> approach would probably be much faster. Just keep that in mind as an alternative.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-05-05T13:06:20.420Z
Answers.LastActivityDate: 2021-05-05T13:06:20.420Z
Answers.UserDisplayName: gringer
Body: <p>Just learning the basics of bioinformatics and bash shell scripting.</p>
  <p>Is there a correct way to search for a specific allele with one or even multiple sequencing runs? I am searching for prevalence of the C2416T allele in the SARS-CoV-2 genome within the ORF1a RNA polymerase gene as a mini project.</p>
  <p>I downloaded a bunch of SARS-COV-2 sequencing runs, I think I need to convert them to bam file type than I am not sure how to call for one specific SNP, especially in a loop for 5-7 runs. Can this be done with bash?</p>
  <p>I was looking into variant calling from a bam file type using vcf and I can't find a way to call for one specific snp on multiple bam files at once.</p>
  <p>Ideally I would want to filter the variant out of the .bcf files or use grep since I sort of understand those the most and I’m not working with a large amount of data or a human genome.</p>
  
Tags: <variant-calling><bash><covid-19><sars-cov-2>
Title: Calling for one specific SNP from multiple sequencing runs
ThreadText: # Calling for one specific SNP from multiple sequencing runs
  
  (Asked by: username_0 on 2021-05-03)
  
  Just learning the basics of bioinformatics and bash shell scripting.
  
  Is there a correct way to search for a specific allele with one or even multiple sequencing runs? I am searching for prevalence of the C2416T allele in the SARS-CoV-2 genome within the ORF1a RNA polymerase gene as a mini project.
  
  I downloaded a bunch of SARS-COV-2 sequencing runs, I think I need to convert them to bam file type than I am not sure how to call for one specific SNP, especially in a loop for 5-7 runs. Can this be done with bash?
  
  I was looking into variant calling from a bam file type using vcf and I can't find a way to call for one specific snp on multiple bam files at once.
  
  Ideally I would want to filter the variant out of the .bcf files or use grep since I sort of understand those the most and I’m not working with a large amount of data or a human genome.
  
  **Tags:** variant-calling, bash, covid-19, sars-cov-2
  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_1 on 2021-05-03. Score: 4
  
  Getting SNPs from whole-genome data is usually done via mapping sequencing reads on the reference (bwa-mem or bowtie are two popular mappers) and then using one of variant calling tools (like freebayes, samtools, or GATK). However, I would also look up if there is something virus-specific.
  
  According to [this answer](https://www.biostars.org/p/93673/#93676), you could get variant calls of a specific region using samtools if you give it a bed file that contains the region of interest (three columns `chr start stop` format). For example using human reference
  
  ```
  samtools mpileup -uDl regions.bed -f hg19.fa file1.sorted.bam file2.sorted.bam | bcftools view -bvcg - > RAL_samtools.raw.bcf
  
  ```
  
  You could also specify a bunch of files, not just two and it would do the variant calling on them all.
  
  
  ---
  
  ### Answer by username_2 on 2021-05-05. Score: 2
  
  *Answer from @atpoint [transferred from comment]:*
  
  Alternatively, you can use `samtools view` to only retain reads that overlap the region you are interested in and then feed these into the mpileup. IIRC there was something with that `-l` option in mpileup that caused trouble but I do not remember what it was. I think it still called all variants based on all reads and then simply filtered for the target regions in `-l`. The `view` approach would probably be much faster. Just keep that in mind as an alternative.
  
  
  ---
  
  