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

CC BY-SA 4.0
2021-05-03T13:15:20.363Z
2021-05-03T13:15:20.363Z
Kamil S Jaron
<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>

CC BY-SA 4.0
2021-05-05T13:06:20.420Z
2021-05-05T13:06:20.420Z
gringer
<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>

<variant-calling><bash><covid-19><sars-cov-2>
Calling for one specific SNP from multiple sequencing runs
# 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, 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.


---

