Calling for one specific SNP from multiple sequencing runs

Question

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.

Answer 1

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