Convert pileup to bed

[gmarco]

Hello,

I'm currently using samtools mpileup like this:

samtools mpileup -l target_regions.bed file1.bam file2.bam file3.bam > mpileup-all

The output format is pileup here's a short example:

chr1	6384820	N	0	*	*	1	G	E	0	*	*	0	*	*	0	*	*	1	G	A0	*	*
chr1	6384821	N	0	*	*	1	G	H	0	*	*	0	*	*	0	*	*	1	G	H0	*	*
chr1	6384822	N	0	*	*	1	C	I	0	*	*	0	*	*	0	*	*	1	C	G0	*	*
chr1	6384825	N	0	*	*	1	C	I	0	*	*	0	*	*	0	*	*	1	C	G0	*	*
chr1	6384826	N	0	*	*	1	C	I	0	*	*	0	*	*	0	*	*	1	C	G0	*	*
chr1	6384827	N	0	*	*	1	C	I	0	*	*	0	*	*	0	*	*	1	C	G0	*	*
chr1	6384828	N	0	*	*	1	C	I	0	*	*	0	*	*	0	*	*	1	C	G0	*	*

Is there any way to convert his pileup format to bed. For this short example the desired output would be:

chr1	6384820	6384822
chr1	6384825	6384828

Thanks.

[AlexReynolds]

Perhaps a basic awk script would work:

$ awk 'BEGIN{\
first_start=$2;\
old_start=first_start;\
}\
{\
chr=$1;\
current_start=$2;\
if (current_start > (old_start + 1)) {\
print chr"\t"first_start"\t"old_start;\
first_start=current_start;\
old_start=first_start;\
}\
else {\
old_start=current_start;\
}\
}\
END {\
print chr"\t"first_start"\t"old_start;\
}' pileup.txt

[dgscofield]

A bed for your example would be

chr1	6384819	6384822
chr1	6384824	6384828

because bed is 0-based, right-open.

You could use my intervalBed script from here in a little pipeline:

samtools mpileup ... \
| cut -f1,2 \
| intervalBed header=0 val_col=0 trackname="pileupIntervals" > pileup.bed

The repository as well as the intervalBed script (written in awk) include more documentation.

[AlexReynolds]

I think dgscofield is right — while BAM is 0-based, the output from SAM is still 1-based. Easy to screw up. Subtracting one from the start coordinate in my awk script should fix the indexing.