3
$\begingroup$

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.

$\endgroup$
0

2 Answers 2

4
$\begingroup$

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.

$\endgroup$
0
2
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.