/usr/share/perl5/Bio/Tradis/Map.pm is in bio-tradis 1.3.3+dfsg-3.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 | package Bio::Tradis::Map;
# ABSTRACT: Perform mapping
=head1 SYNOPSIS
Takes a reference genome and indexes it.
Maps given fastq files to ref.
use Bio::Tradis::Map;
my $pipeline = Bio::Tradis::Map->new(fastqfile => 'abc', reference => 'abc');
$pipeline->index_ref();
$pipeline->do_mapping();
=head1 PARAMETERS
=head2 Required
=over
=item * C<fastqfile> - path to/name of file containing reads to map to the reference
=item * C<reference> - path to/name of reference genome in fasta format (.fa)
=back
=head2 Optional
=over
=item * C<refname> - name to assign to the reference index files. Default = ref.index
=item * C<outfile> - name to assign to the mapped SAM file. Default = mapped.sam
=back
=head1 METHODS
=over
=item * C<index_ref> - create index files of the reference genome. These are required
for the mapping step. Only skip this step if index files already
exist. -k and -s options for referencing are calculated based
on the length of the reads being mapped as per table:
=begin html
<table>
<tr><th>Read length</th><th>k</th><th>s</th></tr>
<tr><td><70</td><td>13</td><td>4<td></tr>
<tr><td>>70 and <100</td><td>13</td><td>6<td></tr>
<tr><td>>100</td><td>20</td><td>6<td></tr>
</table>
=end html
=item * C<do_mapping> - map C<fastqfile> to C<reference>. Options used for mapping are: C<-r -1 -x -y 0.96>
=back
For more information on the mapping and indexing options discussed here, see the L<SMALT manual|ftp://ftp.sanger.ac.uk/pub4/resources/software/smalt/smalt-manual-0.7.4.pdf>
=cut
use Moose;
use Bio::Tradis::Parser::Fastq;
has 'fastqfile' => ( is => 'rw', isa => 'Str', required => 1 );
has 'reference' => ( is => 'rw', isa => 'Str', required => 1 );
has 'refname' =>
( is => 'rw', isa => 'Str', required => 0, default => 'ref.index' );
has 'outfile' =>
( is => 'rw', isa => 'Str', required => 0, default => 'mapped.sam' );
has 'smalt_k' => ( is => 'rw', isa => 'Maybe[Int]', required => 0 );
has 'smalt_s' => ( is => 'rw', isa => 'Maybe[Int]', required => 0 );
has 'smalt_y' => ( is => 'rw', isa => 'Maybe[Num]', required => 0, default => 0.96 );
has 'smalt_r' => ( is => 'rw', isa => 'Maybe[Int]', required => 0, default => -1 );
has 'smalt_n' => ( is => 'rw', isa => 'Maybe[Int]', required => 0, default => 1 );
sub index_ref {
my ($self) = @_;
my $ref = $self->reference;
my $refname = $self->refname;
# Calculate index parameters
my $pars = Bio::Tradis::Parser::Fastq->new( file => $self->fastqfile );
$pars->next_read;
my @read = $pars->read_info;
my $read_len = length($read[1]);
my ( $k, $s ) = $self->_calculate_index_parameters($read_len);
my $cmd = "smalt index -k $k -s $s $refname $ref > /dev/null 2>&1";
system($cmd);
return $cmd;
}
sub _calculate_index_parameters {
my ($self, $read_len) = @_;
my ( $k, $s );
if( defined $self->smalt_k ){ $k = $self->smalt_k; }
else{ $k = $self->_smalt_k_default($read_len); }
if( defined $self->smalt_s ){ $s = $self->smalt_s; }
else{ $s = $self->_smalt_s_default($read_len); }
return ( $k, $s );
}
sub _smalt_k_default {
my ($self, $read_len) = @_;
if($read_len < 100){ return 13; }
else{ return 20; }
}
sub _smalt_s_default {
my ( $self, $read_len ) = @_;
if( $read_len < 70 ){ return 4; }
elsif( $read_len > 100 ){ return 13; }
else{ return 6; }
}
sub do_mapping {
my ($self) = @_;
my $fqfile = $self->fastqfile;
my $refname = $self->refname;
my $outfile = $self->outfile;
my $y = $self->smalt_y;
my $r = $self->smalt_r;
my $n = $self->smalt_n;
my $smalt = "smalt map -n $n -x -r $r -y $y $refname $fqfile 1> $outfile 2> smalt.stderr";
system($smalt);
unlink('smalt.stderr');
return $smalt;
}
__PACKAGE__->meta->make_immutable;
no Moose;
1;
|