Mercurial > repos > jvolkening > b2b_fq_interleave

comparison frag_lens @ 0:6530b8a73225 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/jvolkening/galaxy-tools/tree/master/tools/b2b_utils commit 9260aa02a5f703bce63d2db5b69003df9be371ac
author jvolkening
date Fri, 08 Mar 2024 00:48:17 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:6530b8a73225
1 #!/usr/bin/env perl
2
3 use strict;
4 use warnings;
5 use 5.012;
6
7 use Cwd qw/abs_path/;
8 use File::Temp qw/tempdir tempfile/;
9 use IPC::Cmd qw/can_run/;
10 use List::Util qw/sum/;
11 use Getopt::Long;
12 use Pod::Usage;
13
14 my @good_codes = ( 0x0002, 0x0040 );
15 my @bad_codes = ( 0x0004, 0x0100, 0x0800 );
16
17 #-inputs---------------------------------------------------------------------#
18 my $fasta;
19 my $forward;
20 my $reverse;
21 my $sam;
22 #-knobs----------------------------------------------------------------------#
23 my $threads = 1;
24 my $max_align = 10000;
25
26 my $PROGRAM = 'frag_lens';
27 my $VERSION = 0.001;
28
29 GetOptions(
30 #-inputs-----------------------------------------------------------------#
31 'sam=s' => \$sam,
32 'forward=s' => \$forward,
33 'reverse=s' => \$reverse,
34 'ref=s' => \$fasta,
35 #-knobs------------------------------------------------------------------#
36 'threads=i' => \$threads,
37 'max_aln=i' => \$max_align,
38 'help' => sub{ pod2usage(-verbose => 2); },
39 'version' => sub{ say "This is $PROGRAM v$VERSION";exit; },
40
41 ) or pod2usage( -verbose => 1);
42
43 my $fh_sam;
44 my $tmp_fasta;
45
46 if (defined $sam) {
47 open $fh_sam, '<', $sam or die "failed to open SAM\n";
48 }
49
50 else {
51
52 my $BWA = can_run('bwa')
53 // die "BWA is required but not found\n";
54
55 my ($tmp_dir) = tempdir( CLEANUP => 1);
56
57 die "specify forward and reverse read files and reference\n"
58 if (! defined $forward || ! defined $reverse || ! defined $fasta);
59
60 $fasta = abs_path($fasta);
61
62 my $res = system(
63 'ln',
64 '-s',
65 $fasta,
66 "$tmp_dir/tmp.fasta"
67 );
68 die "link failed" if ($res);
69 $res = system(
70 $BWA,
71 'index',
72 "$tmp_dir/tmp.fasta"
73 );
74 die "index failed" if ($res);
75 open $fh_sam, '-|', $BWA,
76 'mem',
77 '-t' => $threads,
78 '-v' => 1,
79 "$tmp_dir/tmp.fasta",
80 $forward,
81 $reverse
82 ;
83 }
84
85 my $c = 0;
86 while (my $line = <$fh_sam>) {
87 next if ($line =~ /^\@/);
88 chomp $line;
89 my @parts = split "\t", $line;
90 my $flags = $parts[1];
91 my $sum1 = sum map {$_ & $flags ? 1 : 0} @good_codes;
92 my $sum2 = sum map {$_ & $flags ? 1 : 0} @bad_codes;
93 if ($sum1 == scalar @good_codes && $sum2 == 0) {
94 say abs($parts[8]);
95 last if (++$c >= $max_align);
96 }
97 }
98 close $fh_sam;
99
100 __END__
101
102 =head1 NAME
103
104 frag_lens - Calculate paired end fragment lengths from read alignment
105
106 =head1 SYNOPSIS
107
108 frag_lens [--sam <in.sam>] OR [--ref <cons.fa> --forward <R1.fq> --reverse <R2.fq>] [options] > frag_lens.txt
109
110 =head1 DESCRIPTION
111
112 Calculates library fragment lengths based on paired-end read alignment.
113 Takes as input either a preprepared SAM alignment or a reference and read
114 files from which it produces an alignment. Outputs calculated fragment
115 lengths, one per line.
116
117 =head1 PREREQUISITES
118
119 Requires the following binaries:
120
121 =over 1
122
123 =item * bwa
124
125 =back
126
127 =head1 OPTIONS
128
129 =head2 Input option one
130
131 =over 4
132
133 =item B<--sam> I<filename>
134
135 Path to input SAM alignment.
136
137 =back
138
139 =head2 Input option two
140
141 =over 4
142
143 =item B<--ref> I<filename>
144
145 Path to reference sequence (e.g. assembly)
146
147 =item B<--forward> I<filename>
148
149 Forward reads in FASTQ format
150
151 =item B<--reverse> I<filename>
152
153 Reverse reads in FASTQ format
154
155 =back
156
157 =head2 Configuration
158
159 =over 4
160
161 =item B<--max_align>
162
163 Maximum number of alignment records to read as input. Used to limit run times.
164
165 =item B<--threads>
166
167 Number of threads to use for alignment (ignored if --sam is given)
168
169 =back
170
171 =head1 CAVEATS AND BUGS
172
173 Please submit bug reports to the issue tracker in the distribution repository.
174
175 =head1 AUTHOR
176
177 Jeremy Volkening (jdv@base2bio.com)
178
179 =head1 LICENSE AND COPYRIGHT
180
181 Copyright 2014-19 Jeremy Volkening
182
183 This program is free software: you can redistribute it and/or modify
184 it under the terms of the GNU General Public License as published by
185 the Free Software Foundation, either version 3 of the License, or
186 (at your option) any later version.
187
188 This program is distributed in the hope that it will be useful,
189 but WITHOUT ANY WARRANTY; without even the implied warranty of
190 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
191 GNU General Public License for more details.
192
193 You should have received a copy of the GNU General Public License
194 along with this program. If not, see <http://www.gnu.org/licenses/>.
195
196 =cut
197