chip.pl

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

###############################################################################
###############################################################################
### bla bla bla bla bla bla bla bla bla bla bla bla bla bla bla bla         ###
###############################################################################
###############################################################################

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

my (
    # variable for parameters which are read in from commandline
    $help,
    $addVUS,
    );

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

&GetOptions(
    "help!"   => \$help,     # print this help
    "vus!"    => \$addVUS,   # adds VUS to statistics
) or die "\nTry \"$0 -h\" for a complete list of options\n\n";

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

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

##############################################################################
### SETTINGS #################################################################
##############################################################################
my $R = Statistics::R->new();

##############################################################################
### 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";
}

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


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

my %genes = ();
my %muts = ();
my %groups = ();
my %groupIDs = ();
my %groupsN = ();

my $header = <>;
while(my $l = <>) {
    chomp($l);
    my ($groupID, $uszID, $fmiID, $dob, $gender, $disease, $specimen, $mtpStatus,
	$sv, $cna, $rearr, $svVUS, $cnaVUS, $rearrVUS,
	$tmbStatus, $tmb, $msStatus, $msi, $loh,
	$dna) = split(/\t/, $l);

    $sv .= $svVUS if(defined $addVUS);
    $sv =~ s/,$//;
    $groupID =~ s/ /_/g;
    $groupID =~ s/-/_/g;
    $groupID = "none" if($groupID eq "");
    $groupIDs{$groupID} = 1;
    if(exists $groupsN{$groupID}) {
	$groupsN{$groupID}++;
    } else {
	$groupsN{$groupID} = 1;
    }

    my @sv = split(/,/, $sv);
    foreach my $mut (@sv) {
	my $gene = $mut;
	$gene =~ s/ .*//;
	my $pos = $mut;
	$pos =~ s/.* (.*) .*/$1/;
	my $freq = $mut;
	$freq =~ s/.* //;
	$freq =~ s/[\(\)\%]//g;
	my $isSomatic = $freq > 40 ? "no" : "yes";

	# a gene with the same frequency of a mutation should not exist
	if(exists $muts{$isSomatic}->{$gene}->{$freq}->{$pos}) {
	    $muts{$isSomatic}->{$gene}->{$freq}->{$pos}++;
	} else {   # Add muts to results
	    $muts{$isSomatic}->{$gene}->{$freq}->{$pos} = 1;
	}
	
	# Add gene counts to results
	if(exists $genes{$isSomatic}->{$gene}) {
	    $genes{$isSomatic}->{$gene}++;
	} else {
	    $genes{$isSomatic}->{$gene} = 1;
	}

	# Add cohort code to results
	if(exists $groups{$isSomatic}->{$gene}->{$groupID}) {
	    $groups{$isSomatic}->{$gene}->{$groupID}++;
	} else {
	    $groups{$isSomatic}->{$gene}->{$groupID} = 1;
	}
	
	# Create stats for all muts
	$isSomatic = "no+yes";
	$muts{$isSomatic}->{$gene}->{$freq}->{$pos} = 1;
    	if(exists $genes{$isSomatic}->{$gene}) {
	    $genes{$isSomatic}->{$gene}++;
	} else {
	    $genes{$isSomatic}->{$gene} = 1;
	}
	if(exists $groups{$isSomatic}->{$gene}->{$groupID}) {
	    $groups{$isSomatic}->{$gene}->{$groupID}++;
	} else {
	    $groups{$isSomatic}->{$gene}->{$groupID} = 1;
	}
    }
}

$R->startR();
# write out header
print "#IsSomatic(MAF<=40%)\tGene\tN\tMutations";
foreach my $groupID (sort {$a cmp $b} keys %groupIDs) {
    print "\t".$groupID."_count_".$groupsN{$groupID};
}
print "\tChiSqTestPvalue(#IPS vs #no_IPS)\n";

# no do data crunching
foreach my $isSomatic (sort {$b cmp $a} keys %genes) {
    foreach my $gene (sort{$genes{$isSomatic}->{$b}<=>$genes{$isSomatic}->{$a} || $a cmp $b} keys %{$genes{$isSomatic}}) {

	printf("%s\t%s\t%i\t", $isSomatic, $gene, $genes{$isSomatic}->{$gene});

	my $col = "";
	foreach my $freq (sort {$b <=> $a} keys %{$muts{$isSomatic}->{$gene}}) {
	    foreach my $pos (sort keys %{$muts{$isSomatic}->{$gene}->{$freq}}) {
		$col .= sprintf("%s(%s),", $pos, "$freq%");
	    }
	}
	$col =~ s/,$//;
	print $col;
	foreach my $groupID (sort {$a cmp $b} keys %groupIDs) {
	    if(exists $groups{$isSomatic}->{$gene}->{$groupID}) {
		my $count = $groups{$isSomatic}->{$gene}->{$groupID};
		print "\t$count";
		$R->send(qq`$groupID = $count`);
	    } else {
		$R->send(qq`$groupID = 0`);
		print "\t-";
	    }
	}
	if(exists $groups{$isSomatic}->{$gene}->{"IPS"} or exists $groups{$isSomatic}->{$gene}->{"no_IPS"}) {
	    $R->send(qq`test = binom.test(IPS, IPS+no_IPS)`);
	    $R->send(q`cat(test$p.value)`);
	    printf("\t%.2e", $R->read());
	} else {
	    print "\t-";
	}
	print "\n";
    }
}