convert2pcawgGenesetFormatRandomized.pl

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

###############################################################################
###############################################################################
### Converts file to PCAWG-5_GeneSetsFormat.                                ###
###############################################################################
###############################################################################

use strict;
use warnings;
use Getopt::Long;
no autovivification;
use List::Util qw[min max sum];

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

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

&GetOptions(
    "help!"       => \$help,           # print this help
    "type:s"      => \$type,           # if not defined first column will be assumed as type otherwise e.g. Uterus-AdenoCA
    "isDominant!" => \$isDominant,     # input file described lost of interactions due to dominant transfer switch
) or die "\nTry \"$0 -h\" for a complete list of options\n\n";

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

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

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

$| = 1;

# zcat $databases/pcawg/histo/pcawg.rnaseq.extended.metadata.aliquot_id.tsv.gz | tail -n +2 | perl -ne 'chomp($_); @a=split(/\t/); $h{$a[56]} .= "$a[1]:" if($h{$a[56]}!~/$a[1]/); END{print "my \%project = ("; foreach $k(sort keys %h){$h{$k}=~s/\:$//; print "\"$k\", \"$h{$k}\",\n"}}'

##############################################################################
### 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 = "zcat missingInteractionsInDTU_filt_int_anno_mutations.tsv.gz | $0 -type PanCan -min 0.5 -cod\n";

    print "\n\tUsage: " .  $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 %lostIntWithMut = ();
my %int = ();
my @genes = ();
my @l = ();
while(my $l = <>) {
    chomp($l);
    if($l =~ /^#/) {
	next;
    }
    push(@l, $l);

    my @a = split(/\t/, $l);
    push(@genes, "$a[26]\t$a[27]");
}
foreach my $l (@l) {
    my ($tissue, $ensg, $enspCanon, $enstCanon, $enstIso,
	$fdr,
	$dtuCanon, $dtuIso,
	$qValueCanon, $qValueIso,
	$medCanonExpressionInNormal, $medCanonExpressionInCancer, $medIsoformExpressionInNormal, $medIsoformExpressionInCancer,
	$relNumberOfMissedInteractionsOfIsoformInPerc, $relNumberOfMissedInteractionsInNormalInGeneInPerc, $relNumberOfMissedInteractionsInCancerInGeneInPerc,
	$pfam1, $domain1, $region1, $pfam2, $domain2, $enspCanon2, $stringScore, $dimaScore,
	$ensg2, $geneName1, $geneName2,
	$stringIntN1, $stringDensityScore1, $stringRank1, $stringIntN2, $stringDensityScore2, $stringRank2,
	$enspCanonCosmic, $enspCanonCosmicDist, $enspCanonCosmicMinScore,
	$enspCanon2cosmic, $enspCanon2cosmicDist, $enspCanon2cosminMinScore,
	$enspRandom, $enspRandomCosmic, $enspRandomCosmicDist, $enspRandomCosmicMinScore, $randomDist, $randomMinScore,
	$keggCancerPathway1, $keggCancerPathway2,
	$keggPathway1, $keggPathway2,
	$mutationInfos, $eqtlInfos) = split(/\t/, $l);

    $type = $tissue if(!defined $type);

    next if($stringScore eq "-");

    if(defined $isDominant) {
	$relNumberOfMissedInteractionsInCancerInGeneInPerc = $relNumberOfMissedInteractionsInNormalInGeneInPerc/$relNumberOfMissedInteractionsOfIsoformInPerc;
	$qValueCanon = 0;
    }
    my @mutationInfos = split(/,/, $mutationInfos);
    my @eqtlInfos = split(/,/, $eqtlInfos);
    my $sampleBarcode = "-";
    foreach my $mutationInfo (@mutationInfos) {
	next if($mutationInfo eq "-");
	my ($regionType, $chr, $region, $varType, $mRegion, $tumorSampleBarcode) = split(/\:/, $mutationInfo);

	my $gene1 = $geneName1;
	my $gene2 = $geneName2;
	if(($regionType !~ /pc\.cds/ and $regionType !~ /pc\.ss/)) {
	    $randomGenes = $genes[int(rand($#genes + 1))];
	    ($gene1, $gene2) = split(/\t/, $randomGenes);
	}
	$lostIntWithMut{$type}->{$gene1}->{$gene2}->{$regionType}->{$relNumberOfMissedInteractionsInCancerInGeneInPerc}->{$qValueCanon}->{$tissue} = $_;
    }
}

print "#\n";
print "# Result file in PCAWG-5 Gene Sets format listing sets of genes with associated somatic mutations that in addition have overexpressed or\n";
if(defined $isDominant) {
    print "# a switch in the most dominant transcript (MDT) ";
} else {
    print "# differentially used transcripts (DTU) ";
}
print "that partially lack protein binding sites (from DIMA DB v3.0 and 3did-08/2016) found in the canonical isoform.\n";
print "# Data Type:             Functional region of associated mutations (see https://www.synapse.org/#!Synapse:syn5259890).\n";
print "# Sample set:            Cohort.\n";
print "# Gene1:Gene2:...:GeneN: Mutated gene with MDT switch and its protein-protein interaction partners that are unable to bind the mutated gene.\n";
print "# Gene Set Name:         Mutated gene with MDT switch.\n"; 
print "# Score:                 Fraction of interaction partners that are unable to bind mutated gene with MDT switch [values in %].\n";
if(defined $isDominant) {
    print "# FDR:                   -.\n";
} else {
    print "# FDR:                   Mean adjusted p-value of differentially used canonical isoform in entire gene set.\n";
}
print "# Comments:              GTEx normal tissue/PCAWG project code where mutated gene with MDI switch were found.\n";
print "#\n";
print "# For questions, please contact abdullah.kahraman\@gmail.com\n";
print "#\n";

print "#Data_Type\tSample_set\tGene1:Gene2:...:GeneN\tGene_Set_Name\tScore\tFDR\tComments\n";
foreach my $sampleId (sort keys %lostIntWithMut) {
    foreach my $gene1 (sort keys %{$lostIntWithMut{$sampleId}}) {
	my $dataTypes = "";
	my $geneSet = "";
	my $dtu = "";
	my $meanIntPerc = 0;
	my $nInt = 0;
	my $meanFdr = 0;
	my $nFdr = 0;
	my $cancerType = "";
	foreach my $gene2 (sort keys %{$lostIntWithMut{$sampleId}->{$gene1}}) {
	    foreach my $dataType (sort keys %{$lostIntWithMut{$sampleId}->{$gene1}->{$gene2}}) {
		$dtu .= "$gene1:" if($dtu !~ /\b$gene1\b/);
		$geneSet .= "$gene1:" if($geneSet !~ /\b$gene1\b/); 
		$geneSet .= "$gene2:" if($geneSet !~ /\b$gene2\b/); 
		$dataTypes .= "$dataType:" if($dataTypes !~ /\b$dataType\b/);
			    
		foreach my $relMissedInteractionInPerc (sort keys %{$lostIntWithMut{$sampleId}->{$gene1}->{$gene2}->{$dataType}}) {
		    $meanIntPerc += $relMissedInteractionInPerc;
		    $nInt++;
		    foreach my $fdr (sort keys %{$lostIntWithMut{$sampleId}->{$gene1}->{$gene2}->{$dataType}->{$relMissedInteractionInPerc}}) {
			$meanFdr += $fdr;
			$nFdr++;
			foreach my $type (sort keys %{$lostIntWithMut{$sampleId}->{$gene1}->{$gene2}->{$dataType}->{$relMissedInteractionInPerc}->{$fdr}}) {
			    $cancerType .= "$type:" if($cancerType !~ /\b$type\b/);
			}
		    }
		}
	    }
	}
	$dataTypes =~ s/\:$//;
	$geneSet =~ s/\:$//;
	$dtu =~ s/\:$//;
	$cancerType =~ s/\:$//;
	printf("%s\t%s\t%s\t%s\t%.1f\t%s\t%s\n", $dataTypes, $type, $geneSet, $dtu, ($meanIntPerc/$nInt)*100, (defined $isDominant ? "-" : sprintf("%.3e", $meanFdr/$nFdr)), $cancerType) if($dataTypes ne "");
    }
}