--bam-list

Provide a file containing sample names and normal/tumor BAM locations for each. Use the tab- delimited format [sample_name normal_bam tumor_bam] per line. This tool only needs sample_name, so all other columns can be skipped. The sample_name must be the same as the tumor sample names used in the MAF file (16th column, with the header Tumor_Sample_Barcode).

EOS ); }

sub _doc_authors { return <<EOS Nathan D. Dees, Ph.D. Qunyuan Zhang, Ph.D. EOS }

sub execute { # Parse input arguments my $self = shift; my $bam_list = $self->bam_list; my $maf_file = $self->maf_file; my $output_file = $self->output_file; my $gene_list = $self->gene_list; my $permutations = $self->permutations; my @all_sample_names; # names of all the samples, no matter if it's mutated or not my @genes_to_test; # the genes which will be tested for relationships my $skip_non_coding = $self->skip_non_coding; my $skip_silent = $self->skip_silent;

# Parse out the names of the samples which should match the names in the MAF file
my $sampleFh = IO::File->new( $bam_list ) or die "Couldn't open $bam_list. $!\n";
my $line_count = 0;
while( my $line = $sampleFh->getline )
{
    $line_count++;
    next if ( $line =~ m/^#/ );
    chomp( $line );
    my ( $sample ) = split( /\t/, $line );
    if ($sample) {
        push( @all_sample_names, $sample );
    } else {
        warn("could not parse sample name from line $line_count of --bam_list");
    }
}
$sampleFh->close;

# If user-specified, parse out the names of the genes that we are limiting our tests to
if( defined $gene_list )
{
    my $geneFh = IO::File->new( $gene_list ) or die "Couldn't open $gene_list. $!\n";
    while( my $line = $geneFh->getline )
    {
        next if ( $line =~ m/^#/ );
        chomp( $line );
        my ( $gene ) = split( /\t/, $line );
        push( @genes_to_test, $gene );
    }
    $geneFh->close;
}

# Create sample-gene matrix
my $matrix_file = $self->create_sample_gene_matrix($maf_file, \@all_sample_names, \@genes_to_test, $skip_non_coding, $skip_silent);

# Perform mutation-relation test using R
my $R_cmd = "R --slave --args < " . __FILE__ . ".R $matrix_file $permutations $output_file";
print "$R_cmd\n";
WIFEXITED(system $R_cmd) or croak "Couldn't run: $R_cmd ($?)";

return(1);
}

sub create_sample_gene_matrix { my $self = shift; my ( $maf_file, $all_sample_names_ref, $genes_to_test_ref, $skip_non_coding, $skip_silent ) = @_; my @all_sample_names = @{$all_sample_names_ref}; my @genes_to_test = @{$genes_to_test_ref};

# Create hash of mutations from the MAF file
my %mutations;
my %all_genes;

# Parse the MAF file
my $mafFh = IO::File->new( $maf_file ) or die "Couldn't open $maf_file. $!\n";
while( my $line = $mafFh->getline )
{
    next if( $line =~ m/^(#|Hugo_Symbol)/ );
    chomp $line;
    my @cols = split( /\t/, $line );
    my ( $gene, $mutation_class, $sample ) = ( $cols[0], $cols[8], $cols[15] );

    # If the mutation classification is odd, quit with error
    if( $mutation_class !~ m/^(Missense_Mutation|Nonsense_Mutation|Nonstop_Mutation|Splice_Site|Translation_Start_Site|Frame_Shift_Del|Frame_Shift_Ins|In_Frame_Del|In_Frame_Ins|Silent|Intron|RNA|3'Flank|3'UTR|5'Flank|5'UTR|IGR|Targeted_Region|De_novo_Start_InFrame|De_novo_Start_OutOfFrame)$/ )
    {
        print STDERR "Unrecognized Variant_Classification \"$mutation_class\" in MAF file for gene $gene\n";
        print STDERR "Please use TCGA MAF Specification v2.3.\n";
        return undef;
    }

    # If user wants, skip Silent mutations, or those in Introns, RNA, UTRs, Flanks, IGRs, or the ubiquitous Targeted_Region
    if(( $skip_non_coding && $mutation_class =~ m/^(Intron|RNA|3'Flank|3'UTR|5'Flank|5'UTR|IGR|Targeted_Region)$/ ) ||
        ( $skip_silent && $mutation_class =~ m/^Silent$/ ))
    {
        print "Skipping $mutation_class mutation in gene $gene.\n";
        next;
    }

    $all_genes{$gene}++;
    $mutations{$sample}{$gene}++;
}
$mafFh->close;

# If the user specified a gene list, then check for genes that are not in the MAF
if( scalar( @genes_to_test ) > 0 )
{
    for( my $i = 0; $i < scalar( @genes_to_test ); ++$i )
    {
        unless( defined $all_genes{$genes_to_test[$i]} )
        {
            print "Skipping ", $genes_to_test[$i], " from specified gene-list since it was not found in the MAF file\n";
            splice( @genes_to_test, $i, 1 );
        }
    }
}
else
{
    @genes_to_test = sort keys %all_genes;
}

# Write the input matrix to a file for use by the R code
my $matrix_file;
unless( $matrix_file = $self->mutation_matrix_file ) {
    $matrix_file = Genome::Sys->create_temp_file_path();
}

my $matrix_fh = new IO::File $matrix_file,"w";
# Print input matrix file header
my $header = join("\t","Sample",@genes_to_test);
$matrix_fh->print("$header\n");

# Print mutation relation input matrix
for my $sample (sort @all_sample_names) {
    $matrix_fh->print($sample);
    for my $gene (@genes_to_test) {
        if (exists $mutations{$sample}{$gene}) {
            $matrix_fh->print("\t1");
        }
        else {
            $matrix_fh->print("\t0");
        }
    }
    $matrix_fh->print("\n");
}

return $matrix_file;
}

1;