[Giffredo]

Ok, I think it is a good idea..

Anyway thank dpryan for all your suggests. I appreciated them so much

[dpryan]

You might consider starting another thread. I'm unfamiliar with varscan's output and you're more likely to find someone who is if you start a thread mentioning varscan in the subject.

[Giffredo]

before starting to make a script, I tried VarScan..

Reading input from Pileup_merge.pileup (VARSCAN)
chrM 1 G 4 3 G:3:1:21:1:0:3:0
chrM 2 G 4 4 G:4:1:24:1:0:4:0
chrM 3 A 4 4 A:4:1:24:1:0:4:0
chrM 4 T 4 3 T:3:1:23:1:0:3:0
chrM 5 C 2 1 C:1:1:24:1:0:1:0
chrM 6 C 4 3 C:3:1:23:1:0:3:0

mpileup file:
chrM 1 G 4 ^*,^7,^7,^*, 5529
chrM 2 G 4 ,,,, 7;8=
chrM 3 A 4 ,,,, 7:7=
chrM 4 T 4 ,,,, 465=
chrM 5 C 2 ,, 19
chrM 6 C 4 ,,,, 465=

A part that VARSCAN give me error after 300 positions (I ll think after about this error), I don t see the link between these 2 output..

[dpryan]

Yeah, that'd be pretty straight-forward to do in most languages. I should note that there is a C API for this whole process if you prefer that.

[Giffredo]

ok I understood the output that mpileup give me.

Now I think that should exist a "script" to translate the output.. I mean: if I have ....,. in the output it means I have 6 reads in total in that X position in which 5 reads are forward and 1 is reverse. I think that for a normal py script, to do this kind of translation could be easy..

[dpryan]

You'll have a call for every read covering that site, note the depth column. While Phred scores are just log10 values, they're usually encoded in ASCII by first adding 33 to the Phred score. The information provided in a pileup won't be useful for splicing, but it could help with finding RNA editing and other modifications (at least if they're observable in sequencing). Whether you want to use the pileup or just the VCF is sort of up to you.

[Giffredo]

"," is a base call, it just means "the same as the reference, on the reverse orientation". The base quality scores are in the last column, they're phred encoded. "@", for example, is 31.

OK.. but I have still some doubts: I know "," is a base call and for this reason I don t understand why I have more than one "," for only one base in one position...

..and what is Phred code? I know the ASCII code but Phred for me is only a number derived from a logarithmic operation...

My goal is measure the editing sites position, splicing, post transcription modifications in general from my mutant sample mRNA. And I want to make a statistical work on these variations.

[dpryan]

You need to read the manual a bit more. "," is a base call, it just means "the same as the reference, on the reverse orientation". The base quality scores are in the last column, they're phred encoded. "@", for example, is 31.

Regarding mpileup vs. pileup, the output is the same. The only difference is that mpileup can deal with multiple files at once (it just tacks on an extra 3 columns per file).

Perhaps it would help if you mentioned what your actual goal is.

[Giffredo]

I think it is not...
I expected something like this but with the right symbol in the right positions in order to have the possibility to translate and understand them!
where are the read bases? and the quality value? why so many ',' if the position is one?
and..
what mean @ it is not in the legend.. for me there is characters mismatch because of file wrong conversion...

Comes on the question: using pmileup is possible to reach other types of results different from these i reached so far?

[dpryan]

That text file is exactly the output that the manual is describing, so I'm not sure what else you were expecting.

[Giffredo]

chrM 136 A 6 ,,,,,, 896774
chrM 137 A 6 ,,,,,, ?=@===
chrM 138 A 6 ,,,,,, ?=<=8<
chrM 139 C 6 ,,,,,, 887801
chrM 140 C 6 ,,,,,, @>==9;
chrM 141 C 6 ,,,,,, ><==9;
chrM 142 T 6 ,,,,,, 747701

the manual say: each line represents a genomic position, consisting of chromosome name, coordinate, reference base, read bases, read qualities and alignment mapping qualities. Information on match, mismatch, indel, strand, mapping quality and start and end of a read are all encoded at the read base column.

it is not like expected... this txt is not useful at all.

[Giffredo]

Ok thanks!!! I will try.
About & I know.. I put it because of force of habit.. ...

[dpryan]

.. then if I use..

samtools mpileup -f .fas .sorted.bam | tee teeoutput.txt&

...I ll reach the table that I pine for!!

Well, it depends on what sort of table you're after. If you just want the pileup and not variant calls then:

samtools mpileup -f .fas .sorted.bam > output.txt

would seem to do what you want. There'd be no need to pipe things to tee in that case. Also, there's no reason to always end commands with "&", particularly if they'll just be printing stuff to the screen.

[Giffredo]

So, are you telling me that the output from

samtools mpileup -f .fas .sorted.bam

is a .txt?

.. then if I use..

samtools mpileup -f .fas .sorted.bam | tee teeoutput.txt&

...I ll reach the table that I pine for!!

[dpryan]

I just replied on biostars, but this won't work since bcftools is expecting BCF, and you're giving it mpileup, which is text.