view release on metacpan or  search on metacpan

lib/Bio/Tradis/RunTradis.pm  view on Meta::CPAN

      glob("$temporary_directory/$outfile.*.insert_site_plot.gz");

    foreach my $plotname (@current_plots) {
        print STDERR "Reversing $plotname\n" if($self->verbose);

        #my $plotname = $self->_plotname($sn);
        system("gunzip -c $plotname > $temporary_directory/tmp.plot");
        system(
"awk '{ t = \$1; \$1 = \$2; \$2 = t; print; }' $temporary_directory/tmp.plot > rv_plot"
        );
        system("gzip -c rv_plot > $plotname");
    }
    unlink("$temporary_directory/tmp.plot");
    unlink("rv_plot");
}

sub _stats {
    my ($self)                = @_;
    my $outfile               = $self->outfile;
    my $temporary_directory   = $self->_temp_directory;
    my $fq                    = $self->_unzipped_fastq;
    my $seq_info              = $self->_sequence_info;

    #write header to stats file
    $self->_write_stats_header;

    # Add file name and number of reads in it
    my @fql         = split( "/", $fq );
    my $stats       = "$fql[-1],";
    my $total_reads = `wc $fq | awk '{print \$1/4}'`;
    chomp($total_reads);
    $stats .= "$total_reads,";

    # Matching reads
    my $matching = $total_reads;
    if (defined($self->tag)) {
        $matching = `wc $temporary_directory/filter.fastq | awk '{print \$1/4}' `;
        chomp($matching);
    }
    $stats .= "$matching,";
    $stats .= ( $matching / $total_reads ) * 100 . ",";

    # Mapped reads
    my $mapped = $self->_number_of_mapped_reads;
    $stats .= "$mapped,";
    $stats .= ( $mapped / $matching ) * 100 . ",";

    # Unique insertion sites
    my ( $total_uis, $total_seq_len );
    foreach my $si ( keys %{$seq_info} ) {
        my $plotname = $self->_plotname($si);
        system(
"gunzip -c $temporary_directory/$plotname > $temporary_directory/tmp.plot"
        );
        my $uis = `grep -c -v "0 0" $temporary_directory/tmp.plot`;
        chomp($uis);
        $total_uis += $uis;
        $stats .= "$uis,";
        my $seqlen = ${$seq_info}{$si};
        $total_seq_len += $seqlen;
        my $uis_per_seqlen = "NaN";
        $uis_per_seqlen = $seqlen / $uis if ( $uis > 0 );
        chomp($uis_per_seqlen);
        $stats .= "$uis_per_seqlen,";
    }
    $stats .= "$total_uis,";
    my $t_uis_p_l = "NaN";
    $t_uis_p_l = $total_seq_len / $total_uis if ( $total_uis > 0 );
    $stats .= "$t_uis_p_l\n";
    print { $self->_stats_handle } $stats;
}

sub _write_stats_header {
    my ($self)   = @_;
    my @seqnames = keys %{ $self->_sequence_info };
    my @fields   = (
        "File",
        "Total Reads",
        "Reads Matched",
        "\% Matched",
        "Reads Mapped",
        "\% Mapped"
    );
    print { $self->_stats_handle } join( ",", @fields ) . ",";
    foreach my $sn (@seqnames) {
        print { $self->_stats_handle } "Unique Insertion Sites : $sn,";
        print { $self->_stats_handle } "Seq Len/UIS : $sn,";
    }
    print { $self->_stats_handle } "Total Unique Insertion Sites,";
    print { $self->_stats_handle } "Total Seq Len/Total UIS\n";
}

sub _plotname {
    my ( $self, $seq_name ) = @_;
    my $outfile = $self->outfile;

    $seq_name =~ s/[^\w\d\.]/_/g;
    my $plotfile_name = "$outfile.$seq_name.insert_site_plot.gz";
    return $plotfile_name;
}

sub _number_of_mapped_reads {
    my ($self) = @_;
    my $temporary_directory = $self->_temp_directory;

    my $pars =
      Bio::Tradis::Parser::Bam->new(
        file => "$temporary_directory/mapped.bam" );
    my $c = 0;
    while ( $pars->next_read ) {
        if ( $pars->is_mapped ) {
            $c++;
        }
    }
    return $c;
}

__PACKAGE__->meta->make_immutable;
no Moose;
1;

__END__

=pod

=encoding UTF-8