.\" Automatically generated by Pod::Man 4.11 (Pod::Simple 3.35)
.\"
.\" Standard preamble:
.\" ========================================================================
.de Sp \" Vertical space (when we can't use .PP)
.if t .sp .5v
.if n .sp
..
.de Vb \" Begin verbatim text
.ft CW
.nf
.ne \\$1
..
.de Ve \" End verbatim text
.ft R
.fi
..
.\" Set up some character translations and predefined strings.  \*(-- will
.\" give an unbreakable dash, \*(PI will give pi, \*(L" will give a left
.\" double quote, and \*(R" will give a right double quote.  \*(C+ will
.\" give a nicer C++.  Capital omega is used to do unbreakable dashes and
.\" therefore won't be available.  \*(C` and \*(C' expand to `' in nroff,
.\" nothing in troff, for use with C<>.
.tr \(*W-
.ds C+ C\v'-.1v'\h'-1p'\s-2+\h'-1p'+\s0\v'.1v'\h'-1p'
.ie n \{\
.    ds -- \(*W-
.    ds PI pi
.    if (\n(.H=4u)&(1m=24u) .ds -- \(*W\h'-12u'\(*W\h'-12u'-\" diablo 10 pitch
.    if (\n(.H=4u)&(1m=20u) .ds -- \(*W\h'-12u'\(*W\h'-8u'-\"  diablo 12 pitch
.    ds L" ""
.    ds R" ""
.    ds C` ""
.    ds C' ""
'br\}
.el\{\
.    ds -- \|\(em\|
.    ds PI \(*p
.    ds L" ``
.    ds R" ''
.    ds C`
.    ds C'
'br\}
.\"
.\" Escape single quotes in literal strings from groff's Unicode transform.
.ie \n(.g .ds Aq \(aq
.el       .ds Aq '
.\"
.\" If the F register is >0, we'll generate index entries on stderr for
.\" titles (.TH), headers (.SH), subsections (.SS), items (.Ip), and index
.\" entries marked with X<> in POD.  Of course, you'll have to process the
.\" output yourself in some meaningful fashion.
.\"
.\" Avoid warning from groff about undefined register 'F'.
.de IX
..
.nr rF 0
.if \n(.g .if rF .nr rF 1
.if (\n(rF:(\n(.g==0)) \{\
.    if \nF \{\
.        de IX
.        tm Index:\\$1\t\\n%\t"\\$2"
..
.        if !\nF==2 \{\
.            nr % 0
.            nr F 2
.        \}
.    \}
.\}
.rr rF
.\"
.\" Accent mark definitions (@(#)ms.acc 1.5 88/02/08 SMI; from UCB 4.2).
.\" Fear.  Run.  Save yourself.  No user-serviceable parts.
.    \" fudge factors for nroff and troff
.if n \{\
.    ds #H 0
.    ds #V .8m
.    ds #F .3m
.    ds #[ \f1
.    ds #] \fP
.\}
.if t \{\
.    ds #H ((1u-(\\\\n(.fu%2u))*.13m)
.    ds #V .6m
.    ds #F 0
.    ds #[ \&
.    ds #] \&
.\}
.    \" simple accents for nroff and troff
.if n \{\
.    ds ' \&
.    ds ` \&
.    ds ^ \&
.    ds , \&
.    ds ~ ~
.    ds /
.\}
.if t \{\
.    ds ' \\k:\h'-(\\n(.wu*8/10-\*(#H)'\'\h"|\\n:u"
.    ds ` \\k:\h'-(\\n(.wu*8/10-\*(#H)'\`\h'|\\n:u'
.    ds ^ \\k:\h'-(\\n(.wu*10/11-\*(#H)'^\h'|\\n:u'
.    ds , \\k:\h'-(\\n(.wu*8/10)',\h'|\\n:u'
.    ds ~ \\k:\h'-(\\n(.wu-\*(#H-.1m)'~\h'|\\n:u'
.    ds / \\k:\h'-(\\n(.wu*8/10-\*(#H)'\z\(sl\h'|\\n:u'
.\}
.    \" troff and (daisy-wheel) nroff accents
.ds : \\k:\h'-(\\n(.wu*8/10-\*(#H+.1m+\*(#F)'\v'-\*(#V'\z.\h'.2m+\*(#F'.\h'|\\n:u'\v'\*(#V'
.ds 8 \h'\*(#H'\(*b\h'-\*(#H'
.ds o \\k:\h'-(\\n(.wu+\w'\(de'u-\*(#H)/2u'\v'-.3n'\*(#[\z\(de\v'.3n'\h'|\\n:u'\*(#]
.ds d- \h'\*(#H'\(pd\h'-\w'~'u'\v'-.25m'\f2\(hy\fP\v'.25m'\h'-\*(#H'
.ds D- D\\k:\h'-\w'D'u'\v'-.11m'\z\(hy\v'.11m'\h'|\\n:u'
.ds th \*(#[\v'.3m'\s+1I\s-1\v'-.3m'\h'-(\w'I'u*2/3)'\s-1o\s+1\*(#]
.ds Th \*(#[\s+2I\s-2\h'-\w'I'u*3/5'\v'-.3m'o\v'.3m'\*(#]
.ds ae a\h'-(\w'a'u*4/10)'e
.ds Ae A\h'-(\w'A'u*4/10)'E
.    \" corrections for vroff
.if v .ds ~ \\k:\h'-(\\n(.wu*9/10-\*(#H)'\s-2\u~\d\s+2\h'|\\n:u'
.if v .ds ^ \\k:\h'-(\\n(.wu*10/11-\*(#H)'\v'-.4m'^\v'.4m'\h'|\\n:u'
.    \" for low resolution devices (crt and lpr)
.if \n(.H>23 .if \n(.V>19 \
\{\
.    ds : e
.    ds 8 ss
.    ds o a
.    ds d- d\h'-1'\(ga
.    ds D- D\h'-1'\(hy
.    ds th \o'bp'
.    ds Th \o'LP'
.    ds ae ae
.    ds Ae AE
.\}
.rm #[ #] #H #V #F C
.\" ========================================================================
.\"
.IX Title "Genome::Model::Tools::Music::PathScan 3pm"
.TH Genome::Model::Tools::Music::PathScan 3pm "2020-11-06" "perl v5.30.3" "User Contributed Perl Documentation"
.\" For nroff, turn off justification.  Always turn off hyphenation; it makes
.\" way too many mistakes in technical documents.
.if n .ad l
.nh
.IP "\-\-pathway\-file" 4
.IX Item "--pathway-file"
.RS 4
.PD 0
.ie n .IP "This is a tab-delimited file prepared from a pathway database (such as \s-1KEGG\s0), with the columns: [path_id, path_name, class, gene_line, diseases, drugs, description] The latter three columns are optional (but are available on \s-1KEGG\s0). The gene_line contains the ""entrez_id:gene_name"" of all genes involved in this pathway, each separated by a ""|"" symbol." 8
.el .IP "This is a tab-delimited file prepared from a pathway database (such as \s-1KEGG\s0), with the columns: [path_id, path_name, class, gene_line, diseases, drugs, description] The latter three columns are optional (but are available on \s-1KEGG\s0). The gene_line contains the ``entrez_id:gene_name'' of all genes involved in this pathway, each separated by a ``|'' symbol." 8
.IX Item "This is a tab-delimited file prepared from a pathway database (such as KEGG), with the columns: [path_id, path_name, class, gene_line, diseases, drugs, description] The latter three columns are optional (but are available on KEGG). The gene_line contains the entrez_id:gene_name of all genes involved in this pathway, each separated by a | symbol."
.PD
For example, a line in the pathway-file would look like:
.Sp
.Vb 1
\&  hsa00061    Fatty acid biosynthesis    Lipid Metabolism    31:ACACA|32:ACACB|27349:MCAT|2194:FASN|54995:OXSM|55301:OLAH
.Ve
.Sp
Ensure that the gene names and entrez IDs used match those used in the \s-1MAF\s0 file. Entrez IDs are
not mandatory (use a 0 if Entrez \s-1ID\s0 unknown). But if a gene name in the \s-1MAF\s0 does not match any
gene name in this file, the entrez IDs are used to find a match (unless it's a 0).
.RE
.RS 4
.RE
.IP "\-\-gene\-covg\-dir" 4
.IX Item "--gene-covg-dir"
.RS 4
.PD 0
.ie n .IP "This is usually the gene_covgs subdirectory created when you run ""music bmr calc-covg"". It should contain files for each sample that report per-gene covered base counts." 8
.el .IP "This is usually the gene_covgs subdirectory created when you run ``music bmr calc-covg''. It should contain files for each sample that report per-gene covered base counts." 8
.IX Item "This is usually the gene_covgs subdirectory created when you run music bmr calc-covg. It should contain files for each sample that report per-gene covered base counts."
.RE
.RS 4
.RE
.IP "\-\-bam\-list" 4
.IX Item "--bam-list"
.RS 4
.IP "Provide a file containing sample names and normal/tumor \s-1BAM\s0 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 \s-1MAF\s0 file (16th column, with the header Tumor_Sample_Barcode)." 8
.IX Item "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)."
.RE
.RS 4
.RE
.IP "\-\-bmr" 4
.IX Item "--bmr"
.RS 4
.ie n .IP "The overall background mutation rate. This can be calculated using ""music bmr calc-bmr""." 8
.el .IP "The overall background mutation rate. This can be calculated using ``music bmr calc-bmr''." 8
.IX Item "The overall background mutation rate. This can be calculated using music bmr calc-bmr."
.RE
.RS 4
.RE
.IP "\-\-genes\-to\-ignore" 4
.IX Item "--genes-to-ignore"
.RS 4
.IP "A comma-delimited list of genes to ignore from the \s-1MAF\s0 file. This is useful when there are recurrently mutated genes like \s-1TP53\s0 which might mask the significance of other genes." 8
.IX Item "A comma-delimited list of genes to ignore from the MAF file. This is useful when there are recurrently mutated genes like TP53 which might mask the significance of other genes."
.RE
.RS 4
.RE
.PD
.PP
\&\s-1EOS\s0
  );
}
.PP
sub _doc_authors {
    return <<\s-1EOS\s0
 Michael Wendl, Ph.D.
\&\s-1EOS\s0
}
.PP
sub _doc_credits {
    return <<\s-1EOS\s0
This module uses reformatted copies of data from the Kyoto Encyclopedia of Genes and Genomes (\s-1KEGG\s0) database:
.PP
.Vb 3
\& * KEGG \- http://www.genome.jp/kegg/
\&EOS
\&}
.Ve
.PP
sub execute
{
    my \f(CW$self\fR = shift;
    my \f(CW$covg_dir\fR = \f(CW$self\fR\->gene_covg_dir;
    my \f(CW$bam_list\fR = \f(CW$self\fR\->bam_list;
    my \f(CW$pathway_file\fR = \f(CW$self\fR\->pathway_file;
    my \f(CW$maf_file\fR = \f(CW$self\fR\->maf_file;
    my \f(CW$output_file\fR = \f(CW$self\fR\->output_file;
    my \f(CW$bgd_mut_rate\fR = \f(CW$self\fR\->bmr;
    my \f(CW$genes_to_ignore\fR = \f(CW$self\fR\->genes_to_ignore;
    my \f(CW$min_mut_genes_per_path\fR = \f(CW$self\fR\->min_mut_genes_per_path;
    my \f(CW$skip_non_coding\fR = \f(CW$self\fR\->skip_non_coding;
    my \f(CW$skip_silent\fR = \f(CW$self\fR\->skip_silent;
.PP
.Vb 6
\&    # Check on all the input data before starting work
\&    print STDERR "MAF file not found or is empty: $maf_file\en" unless( \-s $maf_file );
\&    print STDERR "Directory with gene coverages not found: $covg_dir\en" unless( \-e $covg_dir );
\&    print STDERR "List of samples not found or is empty: $bam_list\en" unless( \-s $bam_list );
\&    print STDERR "Pathway info file not found or is empty: $pathway_file\en" unless( \-s $pathway_file );
\&    return undef unless( \-s $maf_file && \-e $covg_dir && \-s $bam_list && \-s $pathway_file );
\&
\&    # Build a hash to quickly lookup the genes whose mutations should be ignored
\&    my %ignored_genes = ();
\&    if( defined $genes_to_ignore )
\&    {
\&        %ignored_genes = map { $_ => 1 } split( /,/, $genes_to_ignore );
\&    }
\&
\&    # PathScan uses a helluva lot of hashes \- all your RAM are belong to it
\&    my %sample_gene_hash; # sample => array of genes (based on maf)
\&    my %gene_path_hash; # gene => array of pathways (based on path_file)
\&    my %path_hash; # pathway => all the information about the pathways in the database
\&    my %sample_path_hash; # sample => pathways (based on %sample_gene_hash and %gene_path_hash)
\&    my %path_sample_hits_hash; # path => sample => hits,mutated_genes
\&    my %gene_sample_cov_hash; # gene => sample => coverage
\&    my @all_sample_names; # names of all the samples, no matter if it\*(Aqs mutated or not
\&    my %id_gene_hash; # entrez id => gene (based on first two columns in MAF)
\&
\&    # 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\*(Aqt open $bam_list. $!\en";
\&    while( my $line = $sampleFh\->getline )
\&    {
\&        next if ( $line =~ m/^#/ );
\&        chomp( $line );
\&        my ( $sample ) = split( /\et/, $line );
\&        push( @all_sample_names, $sample );
\&    }
\&    $sampleFh\->close;
\&
\&    # Read coverage data calculated by the Music::Bmr::CalcCovg
\&    $covg_dir =~ s/(\e/)+$//; # Remove trailing forward slashes if any
\&    read_CoverageFiles( $covg_dir, \e@all_sample_names, \e%gene_sample_cov_hash );
\&
\&    #build gene => average_coverage hash for population test
\&    my %gene_cov_hash;
\&    foreach my $gene ( keys %gene_sample_cov_hash )
\&    {
\&        my $total_cov = 0;
\&        my $sample_num = scalar( @all_sample_names );
\&        $total_cov += $gene_sample_cov_hash{$gene}{$_} foreach( @all_sample_names );
\&        $gene_cov_hash{$gene} = int( $total_cov / $sample_num );
\&    }
\&
\&    #build %sample_gene_hash based on maf
\&    my $maf_fh = IO::File\->new( $maf_file );
\&    while( my $line = $maf_fh\->getline )
\&    {
\&        next if( $line =~ m/^(#|Hugo_Symbol)/ );
\&        chomp( $line );
\&        my @cols = split( /\et/, $line );
\&        my ( $gene, $entrez_id, $mutation_class, $tumor_sample ) = ( $cols[0], $cols[1], $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|Silent|In_Frame_Ins|Intron|RNA|3\*(AqFlank|3\*(AqUTR|5\*(AqFlank|5\*(AqUTR|IGR|Targeted_Region|De_novo_Start_InFrame|De_novo_Start_OutOfFrame)$/ )
\&        {
\&            print STDERR "Unrecognized Variant_Classification $mutation_class in MAF file.\en";
\&            print STDERR "Please use TCGA MAF Specification v2.3.\en";
\&            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\*(AqFlank|3\*(AqUTR|5\*(AqFlank|5\*(AqUTR|IGR|Targeted_Region)$/ ) ||
\&          ( $skip_silent && $mutation_class =~ m/^Silent$/ ))
\&        {
\&            print STDERR "Skipping $mutation_class mutation in gene $gene.\en";
\&            next;
\&        }
\&
\&        # Check that the user followed instructions and named each sample correctly
\&        unless( grep( /^$tumor_sample$/, @all_sample_names ))
\&        {
\&            print STDERR "Sample $tumor_sample in MAF file does not match any in $bam_list\en";
\&            return undef;
\&        }
\&
\&        next if( defined $ignored_genes{$gene} ); # Ignore variants in genes that user wants ignored
\&        $id_gene_hash{$entrez_id} = $gene unless( $entrez_id eq \*(Aq\*(Aq or $entrez_id == 0 or $entrez_id !~ m/^\ed+$/ );
\&        push( @{$sample_gene_hash{$tumor_sample}}, $gene ) unless( grep /^$gene$/, @{$sample_gene_hash{$tumor_sample}} );
\&    }
\&    $maf_fh\->close;
\&
\&    my $path_fh = IO::File\->new( $pathway_file );
\&    while( my $line = $path_fh\->getline )
\&    {
\&        chomp( $line );
\&        next if( $line =~ /^(#|ID)/ ); #Skip headers
\&
\&        my ( $path_id, $name, $class, $gene_line, $diseases, $drugs, $description ) = split( /\et/, $line );
\&        my @genes = split( /\e|/, $gene_line ); #Each gene is in the format "EntrezID:GeneSymbol"
\&        $diseases =~ s/\e|/, /g; #Change the separators to commas
\&        $drugs =~ s/\e|/, /g; #Change the separators to commas
\&        $path_hash{$path_id}{name} = $name unless( $name eq \*(Aq\*(Aq );
\&        $path_hash{$path_id}{class} = $class unless( $class eq \*(Aq\*(Aq );
\&        $path_hash{$path_id}{diseases} = $diseases unless( $diseases eq \*(Aq\*(Aq );
\&        $path_hash{$path_id}{drugs} = $drugs unless( $drugs eq \*(Aq\*(Aq );
\&        $path_hash{$path_id}{description} = $description unless( $description eq \*(Aq\*(Aq );
\&        @{$path_hash{$path_id}{gene}} = ();
\&
\&        foreach my $gene ( @genes )
\&        {
\&            my ( $entrez_id, $gene_symbol ) = split( /:/, $gene );
\&            unless( $entrez_id eq \*(Aq\*(Aq or $entrez_id == 0 or $entrez_id !~ m/^\ed+$/ )
\&            {
\&                # Use the gene name from the MAF file if the entrez ID matches
\&                $gene_symbol = $id_gene_hash{$entrez_id} if( defined $id_gene_hash{$entrez_id} );
\&            }
\&            push( @{$gene_path_hash{$gene_symbol}}, $path_id ) unless( grep /^$path_id$/, @{$gene_path_hash{$gene_symbol}} );
\&            unless( grep /^$gene_symbol$/, @{$path_hash{$path_id}{gene}} )
\&            {
\&                push( @{$path_hash{$path_id}{gene}}, $gene_symbol );
\&            }
\&        }
\&    }
\&    $path_fh\->close;
\&
\&    #build a sample => pathway hash
\&    foreach my $sample ( keys %sample_gene_hash )
\&    {
\&        foreach my $gene ( @{$sample_gene_hash{$sample}} )
\&        {
\&            if( defined $gene_path_hash{$gene} )
\&            {
\&                foreach my $pathway ( @{$gene_path_hash{$gene}} )
\&                {
\&                    push( @{$sample_path_hash{$sample}}, $pathway ) unless( grep /^$pathway$/, @{$sample_path_hash{$sample}} );
\&                }
\&            }
\&        }
\&    }
\&
\&    #build path_sample_hits_hash, for population test
\&    foreach my $sample ( keys %sample_path_hash )
\&    {
\&        foreach my $path ( @{$sample_path_hash{$sample}} )
\&        {
\&            my $hits = 0;
\&            my @mutated_genes = (); #Mutated genes in this sample belonging to this pathway
\&            my @mutated_genes_in_sample = @{$sample_gene_hash{$sample}};
\&            foreach my $gene ( @{$path_hash{$path}{gene}} )
\&            {
\&                if( grep /^$gene$/, @mutated_genes_in_sample ) #if this gene is mutated in this sample (in maf)
\&                {
\&                    $hits++;
\&                    push( @mutated_genes, $gene );
\&                }
\&            }
\&            if( $hits > 0 )
\&            {
\&                $path_sample_hits_hash{$path}{$sample}{hits} = $hits;
\&                $path_sample_hits_hash{$path}{$sample}{mutated_genes} = \e@mutated_genes;
\&            }
\&        }
\&    }
\&
\&    #Calculation of p value
\&    my %data; #For printing
\&    my @pvals;
\&    foreach my $path ( sort keys %path_hash )
\&    {
\&        my @pathway_genes = @{$path_hash{$path}{gene}};
\&        my @gene_sizes = ();
\&        foreach my $gene ( @pathway_genes )
\&        {
\&            if( defined $gene_cov_hash{$gene} )
\&            {
\&                my $avg_cov = int( $gene_cov_hash{$gene} );
\&                push( @gene_sizes, $avg_cov ) if( $avg_cov > 3 );
\&            }
\&        }
\&
\&        #If this pathway doesn\*(Aqt have any gene coverage, skip it
\&        next unless( scalar( @gene_sizes ) > 0 );
\&
\&        my @num_hits_per_sample; #store hits info for each patient
\&        my @mutated_samples = sort keys %{$path_sample_hits_hash{$path}};
\&
\&        foreach my $sample ( @all_sample_names )
\&        {
\&            my $hits = 0;
\&            #if this sample has mutation
\&            if( grep /^$sample$/, @mutated_samples )
\&            {
\&                $hits = $path_sample_hits_hash{$path}{$sample}{hits};
\&            }
\&            push( @num_hits_per_sample, $hits );
\&        }
\&
\&        #If this pathway doesn\*(Aqt have any mutated genes in any samples, skip it
\&        next unless( scalar( @num_hits_per_sample ) > 0 );
\&
\&        my $hits_ref = \e@num_hits_per_sample;
\&
\&        ########### MCW ADDED
\&        # FIND MAX NUMBER OF HITS IN A SAMPLE
\&        my $max_hits = 0;
\&        foreach my $hits_in_sample ( @num_hits_per_sample )
\&        {
\&            $max_hits = $hits_in_sample if( $hits_in_sample > $max_hits );
\&        }
\&        ########### MCW ADDED
\&
\&        my $pop_obj = Genome::Model::Tools::Music::PathScan::PopulationPathScan\->new( \e@gene_sizes );
\&        if( scalar( @gene_sizes ) >= 3 )
\&        {
\&            ########### MCW ADDED
\&            if( $max_hits > 15 )
\&            {
\&                $pop_obj\->assign( 5 );
\&            }
\&            else
\&            {
\&                $pop_obj\->assign( 3 );
\&            }
\&            ########### MCW ADDED
\&            #$pop_obj\->assign(3);
\&        }
\&        elsif( @gene_sizes == 2 )
\&        {
\&            $pop_obj\->assign( 2 );
\&        }
\&        else
\&        {
\&            $pop_obj\->assign( 1 );
\&        }
\&
\&        $pop_obj\->preprocess( $bgd_mut_rate, $hits_ref );  #mwendl\*(Aqs new fix
\&
\&        my $pval = $pop_obj\->population_pval_approx($hits_ref);
\&        $data{$pval}{$path}{samples} = \e@mutated_samples;
\&        $data{$pval}{$path}{hits} = $hits_ref;
\&        push( @pvals, $pval ); # For calculation of FDR
\&    }
\&
\&    # Calculate False Discovery Rates (Benjamini\-Hochberg FDR) for the p\-values
\&    my $pval_cnt = scalar( @pvals );
\&    my %fdr_hash;
\&    for( my $i = 0; $i < $pval_cnt; $i++ )
\&    {
\&        my $fdr = $pvals[$i] * $pval_cnt / ( $pval_cnt \- $i );
\&        $fdr = 1 if $fdr > 1;
\&        $fdr_hash{$pvals[$i]} = $fdr;
\&    }
\&
\&    # Print two output files, one more detailed than the other
\&    my $out_fh = IO::File\->new( $output_file, ">" );
\&    my $out_detailed_fh = IO::File\->new( "$output_file\e_detailed", ">" );
\&    $out_fh\->print( "Pathway\etName\etClass\etSamples_Affected\etTotal_Variations\etp\-value\etFDR\en" );
\&
\&    foreach my $pval ( sort { $a <=> $b } keys %data )
\&    {
\&        foreach my $path ( sort keys %{$data{$pval}} )
\&        {
\&            # Skip this pathway if it has fewer affected genes than the user wants
\&            my %mutated_gene_hash;
\&            my @samples = @{$data{$pval}{$path}{samples}};
\&            foreach my $sample ( @samples )
\&            {
\&                foreach my $gene ( @{$path_sample_hits_hash{$path}{$sample}{mutated_genes}} )
\&                {
\&                    $mutated_gene_hash{$gene}++;
\&                }
\&            }
\&            next unless ( scalar( keys %mutated_gene_hash ) >= $min_mut_genes_per_path );
\&
\&            # Print detailed output to a separate output file
\&            $out_detailed_fh\->print( "Pathway: $path\en" );
\&            $out_detailed_fh\->print( "Name: ", $path_hash{$path}{name}, "\en" ) if( defined $path_hash{$path}{name} );
\&            $out_detailed_fh\->print( "Class: ", $path_hash{$path}{class}, "\en" ) if( defined $path_hash{$path}{class} );
\&            $out_detailed_fh\->print( "Diseases: ", $path_hash{$path}{diseases}, "\en" ) if( defined $path_hash{$path}{diseases} );
\&            $out_detailed_fh\->print( "Drugs: ", $path_hash{$path}{drugs}, "\en" ) if( defined $path_hash{$path}{drugs} );
\&            $out_detailed_fh\->print( "P\-value: $pval\en", "FDR: ", $fdr_hash{$pval}, "\en" );
\&            $out_detailed_fh\->print( "Description: ", $path_hash{$path}{description}, "\en" );
\&
\&            my @hits = @{$data{$pval}{$path}{hits}};
\&            foreach my $sample ( @samples )
\&            {
\&                my @mutated_genes = @{$path_sample_hits_hash{$path}{$sample}{mutated_genes}};
\&                $out_detailed_fh\->print( "$sample:" );
\&                $out_detailed_fh\->print( join ",", @mutated_genes );
\&                $out_detailed_fh\->print( "\en" );
\&            }
\&            my ( $mutSampleCnt, $totalMutGenes ) = ( 0, 0 );
\&            $out_detailed_fh\->print( "Samples with mutations (#hits): " );
\&            for( my $i = 0; $i < scalar( @all_sample_names ); ++$i )
\&            {
\&                if( $hits[$i] > 0 )
\&                {
\&                    $out_detailed_fh\->print( "$all_sample_names[$i]($hits[$i]) " );
\&                    $mutSampleCnt++;
\&                    $totalMutGenes += $hits[$i];
\&                }
\&            }
\&            $out_detailed_fh\->print( "\en\en" );
\&
\&            # Print tabulated output to the main output file
\&            my ( $path_name, $path_class ) = ( "\-", "\-" );
\&            $path_name = $path_hash{$path}{name} if( defined $path_hash{$path}{name} );
\&            $path_class = $path_hash{$path}{class} if( defined $path_hash{$path}{class} );
\&            $out_fh\->print( "$path\et$path_name\et$path_class\et$mutSampleCnt\et$totalMutGenes\et",
\&                "$pval\et", $fdr_hash{$pval}, "\en" );
\&        }
\&    }
\&    $out_detailed_fh\->close;
\&    $out_fh\->close;
\&    return 1;
\&}
.Ve
.PP
# Reads files for each sample which are formatted as tab-separated lines each showing the number of
# bases with sufficient coverage in a gene.
sub read_CoverageFiles
{
    my ( \f(CW$covg_dir\fR, \f(CW$all_samples_ref\fR, \f(CW$gene_sample_cov_hash_ref\fR ) = ( \f(CW$_\fR[0], \f(CW$_\fR[1], \f(CW$_\fR[2] );
.PP
.Vb 9
\&    # Read per\-gene covered base counts for each sample
\&    foreach my $sample ( @{$all_samples_ref} )
\&    {
\&        # If the file doesn\*(Aqt exist, quit with error. The Music::Bmr::CalcCovg step is incomplete
\&        unless( \-s "$covg_dir/$sample.covg" )
\&        {
\&            print STDERR "Couldn\*(Aqt find $sample.covg in $covg_dir. (music bmr calc\-covg possibly incomplete)\en";
\&            exit 1;
\&        }
\&
\&        my $covgFh = IO::File\->new( "$covg_dir/$sample.covg" );
\&        while( my $line = $covgFh\->getline )
\&        {
\&            next if( $line =~ m/^#/ );
\&            my ( $gene, undef, $covd_bases ) = split( /\et/, $line );
\&            $gene_sample_cov_hash_ref\->{$gene}{$sample} = $covd_bases;
\&        }
\&        $covgFh\->close;
\&    }
\&}
.Ve
.PP
1;