Hi Guys,
I have a SAM file containing my alignment. What I would like to have is a table containing all positions on my chromosome with information about how many reads contain a certain base. The table below might be better to understand what I mean:
CHR | Position | Number of reads that contain an A | Number of reads that contain a T | Number of reads that contain a C | Number of reads that contain a G.
Do you know any tool that would summarise this for me? I have tried it with mpileup but it can't do it.
Many thanks in advance!
•
written
by
tomislav.ilicic •
Mpileup does exactly what you want. Can you tell us why it won't work? Perhaps a quick perl one liner could fix your problem with mpileup.
EDIT - pileup - not mpileup
EDIT -
samtools pileup dmel.sort.bam > dmel.pileup
The newest versions of samtools have retired pileup so that function doesn't exist any more (which has caused me to need to update a lot of my variant calling pipelines).
Perhaps this is part of the problem for the user?
Really? Hmm, could you perhaps tell me how 1) I can output all positions on the chromosome, not only the snps and 2) where in the output I can find find the number of reads not only covering the base (DP attribute), but also the proportions of each allele? I just can't find it.
This is a comments so it should be posted using the "comment button." I edited my original post. I should have said pileup not mpileup. This should work.
Sorry, I don't often use forms so I forget sometimes. and thank you for the answer!
See this http://en.wikipedia.org/wiki/Pileup_format
5th column can be used to get the allele specific depth. All the positions that have at least one read aligned or spanning it should be in the output file depending on the mapping quality threshold and base quality threshold you have chosen. The position with coverage zero wont be displayed in the output.
thank you for the help, very much appreciated
Hi,
Just FYI. I have run pileup as described above with the default parameters and I want to add that not specifying the reference fasta file (with -f) leads to not containing any positions that are marked to be mapped to the references ("." or ","). In addition it always outputs "N" as the Ref allele, which makes sense as no reference was specified. However, if I do specify the reference file, it shows me the mapped positions, but in total there is much less output than in the other file.
I am not sure how to proceed/which file to take?
Any comments on this?
Thanks!
Yes. If you don't give reference fasta then it should behave in the way you just discussed.
Log in to answer this question.
