count A/T/C/G/N in a aligned site

[johnsequence]

Hello,

I am wondering which tools I should use to count the A/T/C/G/N number/percentage in a given aligned site after alignment. Thanks in advance.

John

[johnsequence]

it looks like no one knows this and if there is not a tool over there, I would like to write a simple one for this basic function and am wondering if any journal like to publish it

[kopi-o]

Have you looked at samtools pileup?

[gringer]

perl's tr can be used for this, if you use it after something like mpileup. It's also possible there's a SAMtool to do it as well.

$valueA = tr/Aa//;$valueC = tr/Cc//;$valueG = tr/Gg//;$valueT = tr/Tt//;
print("A: $valueA; C: $valueC; G: $valueG; T: $valueT\n");

[johnsequence]

Thanks for response. I tried samtools pileup, but the output was a large file, even including a long sequence like....a..t....I only need simple output atcg and their percentage. I might not use it correctly, but I already tried twice and I think perl is great! Thanks.

[Emilie]

Hi John,

You can use mpileup for a given site or position:
samtools mpileup -r 1:18321-18323 bam_file.bam

You will obtain something like that:
1 18321 N 14 TTTTTTTttttttt BHEHEHHICF#HHB
1 18322 N 14 TTTTTTTttttttt DHIHEIIHE@#HGB
1 18323 N 14 CCCCCCCccccccc ,HIHBIIICF#EGB

Then you can cut the 5th column and count the bases as gringer proposed.

Emilie

[johnsequence]

Emilie

Thank you for your post that give me a sense of the output file from mpileup. In this case, if I have one million sites like 1:18321, I just need to count a million line. It should be fast. Good idea!

John