convert_lobstr_trf_bed_to_vcf

#!/usr/bin/perl -w

use warnings;
use strict;
use POSIX;
use Getopt::Long;
use File::Path;
use File::Basename;
use Pod::Usage;

=head1 NAME

convert_lobstr_bed_to_vcf

=head1 SYNOPSIS

 convert_lobstr_bed_to_vcf [options]

=head1 DESCRIPTION

 Converts lobstr trf BED file to VCF file.

=cut

#option variables
my $help;
my $verbose;
my $debug;
my $refFASTAFile;

#initialize options
Getopt::Long::Configure ('bundling');

if(!GetOptions ('h'=>\$help, 'v'=>\$verbose, 'd'=>\$debug,
  'r:s'=>\$refFASTAFile)
  || !defined($refFASTAFile)
  || scalar(@ARGV)!=1)
{
    if ($help)
    {
        pod2usage(-verbose => 2);
    }
    else
    {
        pod2usage(1);
    }
}

my $lobstrTrfBEDFile = $ARGV[0];

open(BED, "$lobstrTrfBEDFile") || die "Cannot open  $lobstrTrfBEDFile";

#note that lobSTR BED file interval is not 0 based half open but instead is 1 based closed

#chr1   10001   10468   6   77.2    6   10001   10468   789 33  51  0   15  1.43    TAACCC  .
#chr1   54713   54817   4   25.8    4   54713   54817   149 0   26  0   72  0.91    TTTC    .
#chr1   66161   66630   2   262.0   2   66161   66630   339 49  0   0   49  1.08    TA  66205,66632,5,202;
#chr1   83792   84041   4   64.8    4   83792   84041   335 72  0   28  0   0.86    AAAG    .
#chr1   99000   99042   4   10.8    4   99000   99042   86  23  0   0   76  0.78    TTTA    .
#chr1   99047   99116   1   70.0    1   99047   99116   68  0   11  0   88  0.51    T   .

#chr1   : chromosome
#10001  : start position
#10468  : end position
#6      : period size
#77.2   : copy number
#789    : purity score
#33     : %tage A
#51     : %tage C
#0      : %tage G
#15     : %tage T
#1.43   :
#TAACCC : STR repeat unit
#.

#CCCAGCCCACCCATCCCATCCCCGCCCACCCATCCCATCCC   <STR>   .   PASS    RU=CCCAT;RL=43;REF=8.6;PSCORE=41;COMP=17,4,9,68;ENTROPY=1.35

my $vcfHdr = <<HDR;
##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##source=TandemRepeatFinder/lobSTR_3.0.2;
##INFO=<ID=RU,Number=1,Type=String,Description="Repeat motif">
##INFO=<ID=MOTIF,Number=1,Type=String,Description="Canonical repeat motif">
##INFO=<ID=RL,Number=1,Type=Integer,Description="Reference STR track length in bp">
##INFO=<ID=REF,Number=1,Type=Float,Description="Reference copy number">
##INFO=<ID=TRF_SCORE,Number=1,Type=Integer,Description="TRF Score for M/I/D as 2/-7/-7">
##INFO=<ID=COMP,Number=4,Type=Integer,Description="Composition(%) of bases in repeat tract">
##INFO=<ID=ENTROPY,Number=1,Type=Float,Description="Entropy measure of repeat tract">
##ALT=<ID=STR,Description="Short Tandem Repeat">
##assembly=/net/fantasia/home/atks/ref/genome/hs37d5.fa
##contig=<ID=1,length=249250621>
##contig=<ID=2,length=243199373>
##contig=<ID=3,length=198022430>
##contig=<ID=4,length=191154276>
##contig=<ID=5,length=180915260>
##contig=<ID=6,length=171115067>
##contig=<ID=7,length=159138663>
##contig=<ID=8,length=146364022>
##contig=<ID=9,length=141213431>
##contig=<ID=10,length=135534747>
##contig=<ID=11,length=135006516>
##contig=<ID=12,length=133851895>
##contig=<ID=13,length=115169878>
##contig=<ID=14,length=107349540>
##contig=<ID=15,length=102531392>
##contig=<ID=16,length=90354753>
##contig=<ID=17,length=81195210>
##contig=<ID=18,length=78077248>
##contig=<ID=19,length=59128983>
##contig=<ID=20,length=63025520>
##contig=<ID=21,length=48129895>
##contig=<ID=22,length=51304566>
##contig=<ID=X,length=155270560>
##contig=<ID=Y,length=59373566>
#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO
HDR

print $vcfHdr;

while(<BED>)
{
    s/\r?\n?$//;

    if (/^#/)
    {
        print "$_\n";
    }
    else
    {
        my ($chrom, $start, $end, $periodSize, $refCopyNo,
            $dummy1, $dummy2, $dummy3,
            $purityScore, $A, $C, $G, $T, $entropy, $repeatUnit, $overlap) = split("\t");

        $chrom =~ s/chr//;
        my $ref = `vt seq -i $chrom:$start-$end -r $refFASTAFile -q`;
        my $alt = "<STR>";

        my $refCopyNoBp = int($refCopyNo * length($repeatUnit) + 0.5);
        my $motif = motif2cannonical($repeatUnit);
        
        print "$chrom\t$start\t.\t$ref\t$alt\t.\tPASS\tRU=$repeatUnit;MOTIF=$motif;RL=$refCopyNoBp;REF=$refCopyNo;TRF_SCORE=$purityScore;COMP=$A,$G,$T,$C;ENTROPY=$entropy\n";
    }
}

close(BED);

sub motif2cannonical
{
    my $motif = shift;
    my $motifrc = reverseComplement($motif);

    my @motifs = ();

    for my $i (0..(length($motif)-1))
    {
        push(@motifs, shiftBase($motif, $i));
    }
    
    @motifs = sort {$a cmp $b} @motifs;
    
    return $motifs[0];
}

sub reverseComplement
{
    my $seq = shift;

    my $len = length($seq);
    my $rc = "";

    for my $i (0..($len-1))
    {
        $rc .= complement(substr($seq, ($len-1-$i),1));
    }

    return $rc;
}

sub shiftBase
{
    my ($seq, $i) = @_;
    my $len = length($seq);
    if ($i>=$len) {$i = $i%%$len;}
    my $shifted = substr($seq, $i, ($len-$i)) . substr($seq, 0, $i);

    return $shifted;
}

sub complement
{
    my $base = shift;

    if ($base eq 'A')
    {
        return 'T';
    }
    elsif ($base eq 'C')
    {
        return 'G';
    }
    elsif ($base eq 'G')
    {
        return 'C';
    }
    elsif ($base eq 'T')
    {
        return 'A';
    }
}