Genome...c::Smg(3pm)

💖 Genome...c::Smg(3pm) 💖

okay so bestie, gather round 🫶 this is the man page for Genome...c::Smg(3pm), filed under User Contributed Perl Documentation, and it is genuinely giving "the documentation is literally the source code." it's a Perl module. it does the SMG thing. it is so much smarter than it looks, no cap.

we've got two documented flags to go over, and then we are reading the actual code together, because that's just how this page rolls. 💅

✨ OPTIONS ✨

• --bmr-modifier-file — bestie can hand over a BMR modifier for each and every gene in the ROI file, and it acts as a multiplier for the categorized background mutation rates, applied before those get tested against the gene's categorized mutation rates. what is it even for? correcting for regional or systematic bias in mutation rates across the genome that might be correlated to CpG deamination, or to DNA repair processes like transcription-coupled repair or mismatch repair 🧬. (mutation rates have also been associated with DNA replication timing, where higher mutation rates are seen in late replicating regions — yeah, we know 👀.) key deets: the same per-gene multiplier gets used on each mutation category of BMR, and any genes from the ROI file that aren't in the BMR modifier file get tested against unmodified overall BMRs per mutation category. also BMR modifiers of <=0 are not permitted, because that's just silly. periodt.
• --skip-low-mr-genes — genes with consistently lower MRs than the BMRs across mutation categories might show up in the results as an SMG (by CT or LRT). if such genes are not of interest to you, they may be assigned a p-value of 1. bonus: this should also speed things up 🏃‍♀️💨. caveat tho, and it matters: genes with higher Indel or Truncation rates than the background will not be skipped, even if the gene's overall MR is lower than the BMR. and if bmr-modifiers are applied, this step uses the modified BMRs instead. iykyk.

🦋 THE CODE 🦋

the page literally opens mid-thought — we get the tail end of some heredoc and then it's straight into the good stuff. iconic, honestly.

EOS
  ); }

🫶 the authors, credit where it's due

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

👑 execute — the main character

the comments in here are doing a lot of the heavy lifting, so i'm just gonna let them speak. note the return 1; at the bottom of the sub, and the lone 1; holding down the end of the file like a little period at the end of a sentence ✨

sub execute {
    my $self = shift;
    my $gene_mr_file = $self->gene_mr_file;
    my $output_file = $self->output_file;
    my $output_file_detailed = $output_file . "_detailed";
    my $max_fdr = $self->max_fdr;
    my $skip_low_mr_genes = $self->skip_low_mr_genes;
    my $bmr_modifier_file = $self->bmr_modifier_file;
    my $processors = $self->processors;

    # Check on all the input data before starting work
    print STDERR "Gene mutation rate file not found or is empty: $gene_mr_file\n" unless( -s $gene_mr_file );
    print STDERR "BMR modifier file not found or is empty: $bmr_modifier_file\n" unless( !defined $bmr_modifier_file || -s $bmr_modifier_file );
    return undef unless( -s $gene_mr_file && ( !defined $bmr_modifier_file || -s $bmr_modifier_file ));

    # If BMR modifiers were provided, then load them, and create another gene_mr_file with modified BMRs
    if( defined $bmr_modifier_file ) {
        my $inBmrModFh = IO::File->new( $bmr_modifier_file ) or die "Couldn't open $bmr_modifier_file. $!\n";
        my %bmr_modifier = ();
        while( my $line = $inBmrModFh->getline ) {
            next if( $line =~ m/^#/ );
            chomp( $line );
            my ( $gene, $modifier ) = split( /\t/, $line );
            ( $modifier > 0 ) or die "$modifier is an invalid bmr-modifier. Please fix values in $bmr_modifier_file.\n";
            $bmr_modifier{$gene} = $modifier;
        }
        $inBmrModFh->close;

        my $new_gene_mr_file = Genome::Sys->create_temp_file_path;
        ( $new_gene_mr_file ) or die "Couldn't create a temp file. $!";
        my $inMrFh = IO::File->new( $gene_mr_file ) or die "Couldn't open $gene_mr_file. $!\n";
        my $outMrFh = IO::File->new( $new_gene_mr_file, ">" ) or die "Couldn't open $new_gene_mr_file. $!\n";
        while( my $line = $inMrFh->getline ) {
            if( $line =~ m/^#/ ) {
                $outMrFh->print( $line );
                next;
            }
            chomp( $line );
            my ( $gene, $type, $covd_bps, $mut_cnt, $bmr ) = split( /\t/, $line );
            $bmr = $bmr * $bmr_modifier{$gene} if( defined $bmr_modifier{$gene} );
            $outMrFh->print( "$gene\t$type\t$covd_bps\t$mut_cnt\t$bmr\n" );
        }
        $outMrFh->close;
        $inMrFh->close;

        $gene_mr_file = $new_gene_mr_file;
    }

    # Collect per-gene mutation rates for reporting in results later
    my ( %gene_muts, %gene_bps, %mut_classes_hash );
    my $inMrFh = IO::File->new( $gene_mr_file ) or die "Couldn't open $gene_mr_file. $!\n";
    while( my $line = $inMrFh->getline ) {
        next if( $line =~ m/^#/ );
        my ( $gene, $type, $covd_bps, $mut_cnt, undef ) = split( /\t/, $line );

        # Warn user about cases where there could be fewer covered bps than mutations detected
        ( $mut_cnt <= $covd_bps ) or warn "More $type seen in $gene than there are bps with sufficient coverage!\n";

        if( $type eq "Overall" or $type eq "Indels" or $type eq "Truncations" ) {
            $gene_muts{$gene}{$type} = $mut_cnt;
            $gene_bps{$gene} = $covd_bps;
            $mut_classes_hash{$type} = 1 unless( $type eq "Overall" );
        }
        elsif( $type =~ m/(Transitions|Transversions)$/ ) {
            $gene_muts{$gene}{SNVs} += $mut_cnt;
            $mut_classes_hash{SNVs} = 1;
        }
        else {
            die "Unrecognized mutation class in gene-mr-file. $!\n";
        }
    }
    $inMrFh->close;
    my @mut_classes = sort keys %mut_classes_hash;

    # Create a temporary intermediate file to hold the p-values
    my $pval_file = Genome::Sys->create_temp_file_path;
    ( $pval_file ) or die "Couldn't create a temp file. $!";

    # Call R for Fisher combined test, Likelihood ratio test, and convolution test on each gene
    my $smg_cmd = "R --slave --args < " . __FILE__ . ".R $gene_mr_file $pval_file smg_test $processors $skip_low_mr_genes";
    WIFEXITED( system $smg_cmd ) or croak "Couldn't run: $smg_cmd ($?)";

    # Call R for calculating FDR on the p-values calculated in the SMG test
    my $fdr_cmd = "R --slave --args < " . __FILE__ . ".R $pval_file $output_file_detailed calc_fdr $processors $skip_low_mr_genes";
    WIFEXITED( system $fdr_cmd ) or croak "Couldn't run: $fdr_cmd ($?)";

    # Parse the R output to identify the SMGs (significant by at least 2 of 3 tests)
    my $smgFh = IO::File->new( $output_file_detailed ) or die "Couldn't open $output_file_detailed. $!\n";
    my ( @newLines, @smgLines );
    my $header = "#Gene\t" . join( "\t", @mut_classes );
    $header .= "\tTot Muts\tCovd Bps\tMuts pMbp\tP-value FCPT\tP-value LRT\tP-value CT\tFDR FCPT\tFDR LRT\tFDR CT\n";
    while( my $line = $smgFh->getline ) {
        chomp( $line );
        if( $line =~ m/^Gene\tp.fisher\tp.lr\tp.convol\tfdr.fisher\tfdr.lr\tfdr.convol$/ ) {
            push( @newLines, $header );
            push( @smgLines, $header );
        }
        else {
            my ( $gene, @pq_vals ) = split( /\t/, $line );
            my ( $p_fcpt, $p_lrt, $p_ct, $q_fcpt, $q_lrt, $q_ct ) = @pq_vals;
            my @mut_cnts;
            foreach( @mut_classes ) {
                # If a mutation count is a fraction, round down the digits after the decimal point
                push( @mut_cnts, (( $gene_muts{$gene}{$_} =~ m/\./ ) ? sprintf( "%.2f", $gene_muts{$gene}{$_} ) : $gene_muts{$gene}{$_} ));
            }
            my $mut_per_mbp = ( $gene_bps{$gene} ? sprintf( "%.2f", ( $gene_muts{$gene}{Overall} / $gene_bps{$gene} * 1000000 )) : 0 );
            push( @newLines, join( "\t", $gene, @mut_cnts, $gene_muts{$gene}{Overall}, $gene_bps{$gene}, $mut_per_mbp, @pq_vals ) . "\n" );

            # If the FDR of at least two of these tests is less than the maximum allowed, we consider it an SMG
            if(( $q_fcpt <= $max_fdr && $q_lrt <= $max_fdr ) || ( $q_fcpt <= $max_fdr && $q_ct <= $max_fdr ) ||
               ( $q_lrt <= $max_fdr && $q_ct <= $max_fdr )) {
                push( @smgLines, join( "\t", $gene, @mut_cnts, $gene_muts{$gene}{Overall}, $gene_bps{$gene}, $mut_per_mbp, @pq_vals ) . "\n" );
            }
        }
    }
    $smgFh->close;

    # Add per-gene SNV and Indel counts to the detailed R output, and make the header friendlier
    my $outDetFh = IO::File->new( $output_file_detailed, ">" ) or die "Couldn't open $output_file_detailed. $!\n";
    $outDetFh->print( @newLines );
    $outDetFh->close;

    # Do the same for only the genes that we consider SMGs
    my $outFh = IO::File->new( $output_file, ">" ) or die "Couldn't open $output_file. $!\n";
    $outFh->print( @smgLines );
    $outFh->close;

    return 1;
}

1;

and that's the whole vibe 🌸 this lil masterpiece reached you via perl v5.30.3, dated 2020-11-06 — shoutout to the era. go forth and find those SMGs, bestie 💖👑🦋