addBindingSiteInfo2gffFile.pl

#!/usr/bin/perl -w
# Author: Abdullah Kahraman
# Date: 15.09.2016

###############################################################################
###############################################################################
### adds DIMA binding site information to a gff file.                       ###  
###############################################################################
###############################################################################

use strict;
use warnings;
use Getopt::Long;
no autovivification;

my (
    # variable for parameters which are read in from commandline
    $help,
    $gffFile,
    $limit,
    $domainFile,
    $stringDimaFile,
    $genomeDir,
    $verbose,

    #global variables
    $gen2aaMap,
    $aa2genMap,
   );

##############################################################################
### read all needed parameters from commandline ##############################

&GetOptions(
    "help!"            => \$help,           # print this help
    "gffFile=s"        => \$gffFile,        # e.g. Homo_sapiens.GRCh37.75.gff
    "limit:s"          => \$limit,          # Limit annotations to lines having this string, e.g. exonic_part
    "domainFile=s"     => \$domainFile,     # e.g. missingIsoformDomains.tsv.gz
    "stringDimaFile=s" => \$stringDimaFile, # e.g. stringDimaInteraction.tsv.gz
    "genomeDir=s"      => \$genomeDir,      # e.g. genome2aa_v75
    "verbose!"         => \$verbose,        # print additional information on computation progress
) or die "\nTry \"$0 -h\" for a complete list of options\n\n";

##############################################################################

# help
if ($help) {printHelp(); exit}

##############################################################################
### SETTINGS #################################################################
##############################################################################

$| = 1;

##############################################################################
### SUBROUTINES ##############################################################
############################################################################## 

###############################################################################
sub printHelp {
###############################################################################
    # prints a help about the using and parameters of this scripts 
    # (execute if user types commandline parameter -h)
    # param:  no paramaters
    # return: no return value

    my (
	$usage,
	$sourceCode,
	@rows,
	$row,
	$option,
	$scriptInfo,
	$example,
       );

    $usage = "$0\n";


    print "\nUsage: " .  $usage . "\n";

    print "Valid options are:\n\n";
    open(MYSELF, "$0") or
      die "Cannot read source code file $0: $!\n";
    $sourceCode .= join "", <MYSELF>;
    close MYSELF;
    $sourceCode =~ s/^.+?\&GetOptions\(\n//s;
    $sourceCode =~ s/\n\).+$//s;
    @rows = split /\n/, $sourceCode;
    foreach $row (@rows){
        $option = $row;
	$option =~ s/\s+\"//g;
	$option =~ s/\"\s.+\#/\t\#/g;
	$option =~ s/=./\t<value> [required]/;
	$option =~ s/:./\t<value> [optional]/;
	$option =~ s/!/\t<non value> [optional]/;

	$row =~ s/^.*//;
	print "\t";
	printf("%-1s%-30s%-30s\n", "-",$option,$row);

    } # end of foreach $row (@rows)
    print "\n";
    print "Options may be abreviated, e.g. -h for --help\n\n";

    $example  = "$0";
}

##############################################################################
sub readGffFile {
    my %gff = ();
    open (F1, $gffFile) or die "ERROR: Failed to open gffFile $gffFile: $!\n\n";
    my $n = 0;
    while(my $l = <F1>) {
	next if($l =~ /^#/);
	chomp($l);
	my ($seqname, $source, $feature, $start, $end, $score, $strand, $frame, $attribute)  = split(/\t/, $l);
	my $chr = $seqname;
	$gff{$n++}->{$chr}->{$start}->{$end} = $l;
    }
    close(F1);

    return \%gff;
}

##############################################################################
sub readDomainFile {
    my %domain = ();
    open (F2, "zcat $domainFile|") or die "ERROR: Failed to open domainFile $domainFile: $!\n\n";
    while(my $l = <F2>) {
	chomp($l);
	my ($stringIdI, $seqLengthI, $stringIdC, $seqLengthC, $stringIdG, $domainInfo, $seqLengthD,
	    $numberOfCanonicalDomains, $numberOfMissingDomains, $ratio, $intN, $intR) = split(/\t/, $l);
	$stringIdI =~ s/^.*\.//;
	$stringIdG =~ s/^.*\.//;
	$domain{$stringIdG}->{$domainInfo} = $l;
    }
    close(F2);

    return \%domain;
}
##############################################################################
sub readStringDimaFile {
    my %stringDima = ();
    open (F3, "zcat $stringDimaFile|") or die "ERROR: Failed to open the stringDima File $stringDimaFile: $!\n\n";
    while(my $l = <F3>) {
	chomp($l);
	my ($ensp1, $seqLength1, $ensg1, $domainInfo1, $domainLength1, $numberOfCanonicalDomains1,
	    $ensp2, $seqLength2, $ensg2, $domainInfo2, $domainLength2, $numberOfCanonicalDomains2,
	    $stringScore, $dimaScore) = split(/\t/, $l);

	$stringDima{$domainInfo1}->{$ensg1}->{$ensg2} = $domainInfo2;
	$stringDima{$domainInfo2}->{$ensg2}->{$ensg1} = $domainInfo1;
    }
    close(F3);

    return \%stringDima;
}

##############################################################################
sub read_genome2aaMap {

    (my $file) = @_;
    
    print STDERR "Reading in $file\n" if($verbose);
    open(F4, "zcat $file |") or die "\nERROR: Failed to unzip genome2aaMap file $file: $!\n\n";
    my $map = ();
    my ($ensp, $ensg, $enst, $geneNames);
    while(my $l = <F4>) {
	chomp($l);
	next if($l =~ /^#/);
	if ($l =~ /^>(.*)/) {
	    ($ensp, $ensg, $enst, $geneNames) = split(/\s+/, $1);
	    next;
	}
	(my $chr, my $genomePos, my $strand, my $codonNo, my $nt, my $aa, my $aaNo) = split(/\t/, $l);
	$gen2aaMap->{$chr}->{$genomePos}->{$ensg} = $l;
	$aa2genMap->{$ensg}->{$aaNo}->{$codonNo} = $nt;
    }
    close(F4);
}

##############################################################################
### END OF SUBROUTINES########################################################
############################################################################## 


############
### MAIN ###
############

if(!defined $gffFile or !defined $domainFile or !defined $genomeDir) {
    print STDERR "\n\tERROR: Missing options. Try $0 -help to get the list of all required options\n\n";
    exit 0;
}
print STDERR "Reading in StringDima interaction file ... " if($verbose);
my $interaction = &readStringDimaFile();
print STDERR "done\n" if($verbose);
print STDERR "Reading in Domain file ... " if($verbose);
my $domain = &readDomainFile();
print STDERR "done\n" if($verbose);
print STDERR "Reading in files in genome directory ...\n" if($verbose);
&read_genome2aaMap("$genomeDir/1_genome2aaMap.txt.gz");
&read_genome2aaMap("$genomeDir/2_genome2aaMap.txt.gz");
&read_genome2aaMap("$genomeDir/3_genome2aaMap.txt.gz");
&read_genome2aaMap("$genomeDir/4_genome2aaMap.txt.gz");
&read_genome2aaMap("$genomeDir/5_genome2aaMap.txt.gz");
&read_genome2aaMap("$genomeDir/6_genome2aaMap.txt.gz");
&read_genome2aaMap("$genomeDir/7_genome2aaMap.txt.gz");
&read_genome2aaMap("$genomeDir/8_genome2aaMap.txt.gz");
&read_genome2aaMap("$genomeDir/9_genome2aaMap.txt.gz");
&read_genome2aaMap("$genomeDir/10_genome2aaMap.txt.gz");
&read_genome2aaMap("$genomeDir/11_genome2aaMap.txt.gz");
&read_genome2aaMap("$genomeDir/12_genome2aaMap.txt.gz");
&read_genome2aaMap("$genomeDir/13_genome2aaMap.txt.gz");
&read_genome2aaMap("$genomeDir/14_genome2aaMap.txt.gz");
&read_genome2aaMap("$genomeDir/15_genome2aaMap.txt.gz");
&read_genome2aaMap("$genomeDir/16_genome2aaMap.txt.gz");
&read_genome2aaMap("$genomeDir/17_genome2aaMap.txt.gz");
&read_genome2aaMap("$genomeDir/18_genome2aaMap.txt.gz");
&read_genome2aaMap("$genomeDir/19_genome2aaMap.txt.gz");
&read_genome2aaMap("$genomeDir/20_genome2aaMap.txt.gz");
&read_genome2aaMap("$genomeDir/21_genome2aaMap.txt.gz");
&read_genome2aaMap("$genomeDir/22_genome2aaMap.txt.gz");
&read_genome2aaMap("$genomeDir/X_genome2aaMap.txt.gz");
&read_genome2aaMap("$genomeDir/Y_genome2aaMap.txt.gz");
print STDERR "Done\n" if($verbose);
print STDERR "Starting with mapping ...\n" if($verbose);

open (F1, $gffFile) or die "ERROR: Failed to open gffFile $gffFile: $!\n\n";
while(my $l = <F1>) {
    next if($l =~ /^#/);
    chomp($l);
    my ($seqname, $source, $feature, $start, $end, $score, $strand, $frame, $attribute)  = split(/\t/, $l);
    my $chr = $seqname;

    print $l;
    if(defined $limit and $l !~ /$limit/) {
	print "\n";
	next;
    } else {
	my $domainOut = "; pfamDomain \"";
	my $intOut = "; stringIteractionPartner \"";
	for(my $i = $start; $i <= $end; $i++) {
	    if(exists $gen2aaMap->{$chr}->{$i}) {
		foreach my $ensg (sort keys %{$gen2aaMap->{$chr}->{$i}}) {
		    if(exists $domain->{$ensg}) {
			foreach my $domainInfo (sort keys %{$domain->{$ensg}}) {
			    $domainOut .= "$domainInfo+" if($domainOut !~ /\b$domainInfo\b/);
			    if(exists $interaction->{$domainInfo}->{$ensg}) {
				foreach my $ensg2 (sort keys %{$interaction->{$domainInfo}->{$ensg}}) {
				    $intOut .= "$ensg2+" if($intOut !~ /\b$ensg2\b/);
				}
			    }
			}
		    }
		}
	    }
	}
	$domainOut =~ s/\+$//g;
	$intOut =~ s/\+$//g;
	print "$domainOut\"$intOut\"\n";
    }
}
close(F1);
print STDERR "ALL DONE\n" if($verbose);