convert_vntrseek_resource_2_vcf

#!/usr/bin/perl

use warnings;
use strict;
use POSIX qw(ceil floor);
use File::Basename;
use File::Path;
use Getopt::Long;
use Pod::Usage;

=head1 NAME

convert_vntrseek_resource_2_vcf

=head1 SYNOPSIS

convert_vntrseek_resource_2_vcf [options] <filename>

    -v  verbose
    -d  debug
    -o  output contig directory.

    usage: convert_vntrseek_resource_2_vcf a.txt  -o output.bed
    
    Converts VNTRseek resource file to a BED file..

=head1 DESCRIPTION

=cut

#option variables
my $help;
my $outputVCFFile;


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

if(!GetOptions('h'=>\$help)
    || scalar(@ARGV)!=1)
{
    if ($help)
    {
        pod2usage(-verbose => 2);
    }
    else
    {
        pod2usage(1);
    }
}

my $inputFile = $ARGV[0];

open(IN, "$inputFile") || die "Cannot open $inputFile";


my $vcfHdr = <<HDR;
##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##source=TandemRepeatFinder 4.07b/VNTRSeek 1.08;
##INFO=<ID=MOTIF,Number=1,Type=String,Description="Canonical repeat motif">
##INFO=<ID=RU,Number=1,Type=String,Description="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=SCORE,Number=1,Type=Float,Description="Score of repeat tract">
##ALT=<ID=VNTR,Description="Variable Number of 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;

#Repeatid,FirstIndex,LastIndex,CopyNumber,FastaHeader,FlankingLeft1000,Pattern,ArraySequence,FlankingRight1000,Conserved
#175344010,
#1133237,
#1133267,
#4.428571,
#chr1,
#TGTATCCATGAGGCTCACCCCCTGTGCAGCTAGGATGGAACTGCACCCCCTCCTCAGAGCAGCCCCCTCAGGCCAGAGCTTCAGACCCTCAGAGCTTCAGTGGACAAACAGGGGAGCAAGGCTATAGGAAGGAAGCAATCTCCCAGAGCCATCTGTCTGAATGAGGACCAGGGCCTCTGCCTCCTGCTGGCCCAGTGGGCTCGGGGTGAGGAGGGACTGACTGGCGTCTGGGGCAGGTATGGCTGCTGGAGATGAATTCTAACCCAGCCCTGCACACCAACTGCGAGGTCCTGAAGGAGGTCATCCCAGGTGTGGTCATCGAGACCCTGGGTGAGCCTCCAAGCCCCCACCCCACACCCCCACAACCTGCCCAGTCTTGAATGCCTACCTGCCCTCTTCCCAAGTTCAGACCCCCACAGATGGGGCTCTTCCACCTCCCGAGGGCCTGTGTGGCAGCGCCGCcccctgccccaaccttccccacacctcagctcacctgggtcccacccctcacaccctccccaccccctccagccccagcctctgccccccactcaccctctctcccccacccGCACCAGACCTGGTGCTCGAGACCTTCCGGAAGAGCCTGCGCGGCCAGAAGATGTTGCCTCTGCTGTCCCAGCGCCGCTTCGTGCTCCTGCACAACGGTGAGGCCGACCCGCGGCCGCACCTGGGGGGCTCGTGCAGCCTCCGCCGCTGGCCGCCCCTGCCCACCCGCCAGGCCAAGTCCTCCGGGCCACCCATGCCGCATGCCCCAGACCAGCCGGGCGCCCGCAGGCCTGCGCCACCTCCCTTGGTGCCGCAGCGTCCCCGGCCACCCGGCCCCGACCTGGACAGCGCCCACGATGGGGAGCCCCAGGCCCCGGGCACGGAGCAGTCGGGCACAGGCAACAGGCACCCGGCGCAAGAGCCTTCCCCGGGGACAGCCAAGGAGGAACGCGAGGAGCCTGAGAACGCGAGGCCCTAGGGGCAGCCA,
#CCCGCGC,
#CCCGCGCCCAGCGCCCCGCGCCCCGCGCCCC,
#AGCCGTGCTGCCTGCCCTCAGGGACCTATAAAGCCCACTTTGCTACAAACACAGTCTCTGCAGCAGAGTCCTGCCACCTCAGCGCCCCTGAGCCCATGCCCACCGGGGTCACCCCCACCACACCAGCACTGGCCTGAAAGATGCTCCCCTGGGACTTCCCCAGGCCCCCACTAGCCCCCAGTCCATGGCTGGGCCCTGTCCCTGGGGAGACTGTGCTGAGGGAGGGGGTCTCCGCGGCTGGGCCCAGTCCTGGGGGGCTGTTCTGAGGGAGGGGGTTCCCCACGGCTGGTCTATACTCCCTCCCGGGCTGCCCCATCCTGGGCTCCCGCCCACCACCCTGTGCCCCTCAGCACCTGGGCCCTCCTGCCGCCCTGTGACCCCACGGCCCCCATGGCCTGGCCGGTGATCCACCATTGCCTGTCATCCCCAGGAAGGGTCTGAACCCCTCTGGAGCTCCAAAGGAGGAGCCCCTTGAGCCCAGGATTGCTGCCGTGGGCACCCACAGAGCTCGCCCCAGCCTGCCTGGCACAGAACCTTCCTGGGCCTTTAAGCCCCAAGCCTGGCGCTCCAGCTCCCTGAGCCCCCAACCTTCCTGCAGTGGCCCCAGGACTTAGCCAGGGGTGAGGGGCACAGGGCAGCTAGCCCCAACAAGAGGGTCCCACAGAAGACATGCTAAGCCCACACTGGCCAGGCCCTGGTCAGCCCTGTGTAGCCGCTGATGGTTCTGTTCCCCATGGGTCTGGGGGTTTGCGAGGAGAGCTCTGGCGTCCCGGGGCCTCTGAGGCGTCCACACTACCCACCCCCTCACTCATCTCCCACCCAGAGAGCTAAAGTGTATGGGCCTCTTCCCAATGGCTGGCTCAAGGTGTCTGTGCCCCCCAAGTCGAGACCTGAGGATCGGGGGGCCCCGCCCCCAGGCCCACCCCCCACATGCCCTGCCTCCGACTCCACACGGTCACGGCCAAGTCGCCGGTCCCTCCTCCTGCCCTGCTCCGCAGCC,
#0.967742

while (<IN>)
{
    if ($.!=1)
    {
        chomp;
        
        my ($repeatID, $beg1, $end1, $refCopyNo, $chrom, $lflank, $repeatUnit, $repeatTract, $rflank, $score) = split(",");
        
        $chrom =~ s/chr//;
        
        my $motif = motif2cannonical($repeatUnit);
        my $repeatLength = length($repeatTract);
        
        print "$chrom\t$beg1\t.\t$repeatTract\t<VNTR>\t.\tPASS\tMOTIF=$motif;RU=$repeatUnit;RL=$repeatLength;REF=$refCopyNo;SCORE=$score\n";  
    }
}

close(IN);

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';
    }
}