Get All Alleles Of 1 Column In Bam File

Question

[Manon Delahaye]

Hi all, I have bam files from which I would like to extract 1 specific column (ie chr:position). I won't to have all alleles, even if they come from low quality reads. And even if the position is not a SNP (ie all alleles are the same or almost).

I tried to use "samtools view" giving the position but then the output is all reads that overlap the position so it's not exactly what I am looking for.

Does anyone has any idea how I could solve my problem ?

Thanks a lot. M.D

Answer 1

[Manon Delahaye]

I managed to do what I was hoping for.

samtools mpileup -Bf myReference.fasta -r chr9:myPos-myPos [?]

-B allows for all location (not only SNPs)

Then with some combination of grep and others I could count my alleles.

Thank you for your help ! M.D

Answer 2

[Sean Davis]

I think you are looking for:

samtools mpileup ....

There are a number of flags that you may want to tweak to get exactly the output you want.

Answer 3

[Manon Delahaye]

Hi, Thanks I will have a look. However, I thought that this command is for SNP calling. And I want the alleles at the position even if there is no variant at that position.

Do you think that is included in one of the option ?

I will have a better and closer look tomorrow.

Thanks. MD

[Rm]

samtools pileup -cf reference.fasta Input.bam

[Sean Davis]

Just a note that pileup is deprecated (or not even available) in newer samtools versions. Mpileup is the replacement.

Answer 4

[guns.guru]

Hi M.D

It's been a long time, since you posted this but I have the exact same problem i.e., to extract the allele at particular positions, whether they are SNPs or not.

How did you achieve this?

Thanks

[Devon Ryan]

Please start a new thread. Also, the answer is samtools mpileup, as the others mentioned above.

[guns.guru]

Thanks for looking into this dpryan! I did start a new thread here: Extract allele at particular positions from BAM files I did infer that "samtools mpileup" is the way to go... however, I couldn't make it work! I'll post more details on that thread. Hope you could take a look there again...