calculate total number of unique positions in target regions

[rama]

Hi there,

I have a bed file containing target regions and I want to calculate the total number of bases covered by these regions(total length of target regions) and also the number of unique base positions covered(since some of the target regions are over lapping).

thank you in advance for your kind help.
Rama

[vivek_]

For unique bases:

cat <File.bed> | awk '{for(i=$1;i<$2;i++) print i}' | sort | uniq | wc -l

For all bases:

cat <File.bed> | awk '{for(i=$1;i<$2;i++) print i}' | wc -l

[rama]

Hi Vivek,

Thanks for your reply.

I am not familiar with awk, I tried your commands but both gave 0 as result.
here I have copied first few lines from my bed file.

chr1 43814962 43815050 . MPL
chr1 115256504 115256554 . NRAS
chr1 115258717 115258765 . NRAS
chr2 29432657 29432711 . ALK
chr2 29443688 29443741 . ALK
chr2 209113105 209113153 . IDH1

don't know if I should change any of the arguments with the command.



thanks again for you

[vivek_]

Hi,

Try these:

For unique bases:

cat yourFile.bed | awk '{for(i=$2;i<$3;i++) print $1"\t"i}' | sort | uniq | wc -l

For all bases:

cat yourFile.bed | awk '{for(i=$2;i<$3;i++) print $1"\t"i}' | wc -l

[rama]

Thank you Vivek, it worked like a charm.

just one more request, how do I group by chromosome to see the number of bases per chromosome covered by the target regions.

[vivek_]

This should work

cat yourFile.bed | awk '{for(i=$2;i<$3;i++) print $1"\t"i}' | sort | uniq | awk '{print $1}' | sort | uniq -c

[rama]

Thank you vivek.

[tedtoal]

vivek_'s solution works, but it is very slow for large BED files. Here's a short Perl program to do it, along with bedtools:

    RGN_SIZE='
        $chr = "";
        $sum = 0;
        $sumAll = 0;
        $start = 0;
        while ($ln = <STDIN>){
            chomp $ln;
            @A = split("\t", $ln);
            if ($A[0] ne $chr){
                if ($chr ne ""){
                    print("$chr\t$sum\n");
                    $sumAll += $sum;
                }
                $chr = $A[0];
                $sum = 0;
                $start = 0; # 0-based, next position after last one counted.
            }
            if ($A[1] > $start){ $start = $A[1]; }
            if ($A[2] > $start){ $sum += $A[2] - $A[1]; $start = $A[2]; }
        }
        print("$chr\t$sum\n");
        $sumAll += $sum;
        print("Total: $sumAll\n");
        '

    bedtools sort -i <bedfile_name> | perl -e "$RGN_SIZE"