Finding some specific SNPs in aligned sam file

Question

[hasani.iut6]

Hi

I have used SNPlocs.Hsapiens.dbSNP144.GRCh38 package in Rstudio to access to the SNPs of chromosome one, furthermore I have aligned sam file that I am interested to find my selected SNPs in them. In other words, for each aligned read, I'm interested to find all of my SNPs in them and also the allele of that SNPs in the reads.

How can I do this?

Thanks all.

EDIT:

Suppose that I can make a VCF file from my specific SNPs and then by using BEDtools I can find overlap between my SNPs VCF file and SAM file of aligned reads but I want a better solution .

[Titus]

You can use a variant caller like Varscan or GATK

[hasani.iut6]

Using GATK cause to variant calling but I have a list of SNPs, that I am interested in find them in my aligned reads

[ATpoint]

Use your aligned reads to call variants on them. There are plenty of tools and tutorials out there. Spend some quality time on reading them. Also, dbSNP144 is ancient. dbSNP151 is the current one for hg38, and it contains notably more variants as dbSNP150, I think because liftovers from hg19 versions of the database were included. Also, it includes allele frequencies for each variant based on the TOPMED project, a large-scale sequencing endeavour by the NIH, some thirty-thousand samples sequenced. This will allow you to make better statements on the abundance of the SNPs in the population.

[hasani.iut6]

Thank you for your response and for saying about dbSNP151.

Actually, I am interested to find some specific SNPs in my aligned reads. I don't want to do variant calling in my aligned reads.

[WouterDeCoster]

Finding SNPs IS variant calling.

[hasani.iut6]

Yes, you are right and maybe I said my question wrongly. I want to find my specific SNPs in aligned reads.

[finswimmer]

Hello hasani.iut6 ,

could you please explain a little bit more about what is your goal?

Doe you want to know which SNPs are in your alignment file? Do you want to extract reads that support a given variant? ...?

fin swimmer

[hasani.iut6]

Thank you for your response fin.

Actually, I want to know which reads contain my which SNPs and I want to know which nucleotide is in my aligned reads in the position of SNPs.

[finswimmer]

Like this?

Answer 1

[Benn]

The best solution is to do variant calling first, like everyone is suggesting in the comments. There are pipelines defined to do this variant calling, so why reinvent the wheel? When you have all the variants called, you can annotate your variants, and using dbSNP you can find your specific SNPs. Use like people suggest here already in the comments something like GATK or VARSCAN followed by annovar.