Skip to content

Add optional parallel genome alignment support #86

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
SergeWielhouwer wants to merge 3 commits into
base: development
Choose a base branch
from
Open

Add optional parallel genome alignment support #86

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
ab594b2
Initial plan
Copilot May 6, 2026
3669fa9
Add --parallel_genomes option with Parallel::ForkManager support
Copilot May 6, 2026
44fd25c
Fix typo: amibiguous -> ambiguous in error message
Copilot May 6, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
212 changes: 147 additions & 65 deletions fastq_screen
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -72,33 +72,35 @@ my $pass;
my $inverse = 0;
my $get_genomes;
my $add_genome;
my $parallel_genomes = 1;

my $config_result = GetOptions(
"subset=i" => \$subset_count,
"outdir=s" => \$outdir,
"illumina1_3" => \$illumina,
"quiet" => \$quiet,
"help" => \$help,
"version" => \$version,
"conf=s" => \$conf,
"bowtie=s" => \$bowtie_opts,
"bowtie2=s" => \$bowtie2_opts,
"bismark=s" => \$bismark_opts,
"bwa=s" => \$bwa_opts,
"minimap2=s" => \$minimap2_opts,
"threads=i" => \$threads,
"nohits" => \$nohits,
"aligner=s" => \$aligner,
"force" => \$force,
"paired" => \$paired,
"bisulfite" => \$bisulfite,
"tag" => \$tag,
"filter=s" => \$filter,
"top=s" => \$top,
"pass=i" => \$pass,
"inverse" => \$inverse,
"get_genomes" => \$get_genomes,
"add_genome=s" => \$add_genome,
"subset=i" => \$subset_count,
"outdir=s" => \$outdir,
"illumina1_3" => \$illumina,
"quiet" => \$quiet,
"help" => \$help,
"version" => \$version,
"conf=s" => \$conf,
"bowtie=s" => \$bowtie_opts,
"bowtie2=s" => \$bowtie2_opts,
"bismark=s" => \$bismark_opts,
"bwa=s" => \$bwa_opts,
"minimap2=s" => \$minimap2_opts,
"threads=i" => \$threads,
"nohits" => \$nohits,
"aligner=s" => \$aligner,
"force" => \$force,
"paired" => \$paired,
"bisulfite" => \$bisulfite,
"tag" => \$tag,
"filter=s" => \$filter,
"top=s" => \$top,
"pass=i" => \$pass,
"inverse" => \$inverse,
"get_genomes" => \$get_genomes,
"add_genome=s" => \$add_genome,
"parallel_genomes=i" => \$parallel_genomes,
);

die "Could not parse options, please adjust configuration.\n" unless ($config_result);
Expand Down Expand Up @@ -239,6 +241,8 @@ if( (defined $bismark_opts) and (!defined $bisulfite) ){
die "Option --bismark may not be specified without --bisulfite\n";
}

die "Option --parallel_genomes must be a positive integer.\n" if ($parallel_genomes < 1);

$bowtie_opts = '' unless (defined $bowtie_opts); # Get undef warnings otherwise
$bowtie2_opts = '' unless (defined $bowtie2_opts);
$bwa_opts = '' unless (defined $bwa_opts);
Expand Down Expand Up @@ -339,6 +343,13 @@ if ( ( ( 1 << 32 ) != 4294967296 ) and ( scalar @libraries > 15 ) ) { #32-bit
die "Maximum number of reference genomes exceeded, please adjust configuration and specify at most 32 libraries.\n";
}

if ($parallel_genomes > 1) {
my $total_cpus = $parallel_genomes * $number_of_threads;
warn "Running $parallel_genomes genome alignments in parallel (total CPU usage approximately $total_cpus threads)\n" unless ($quiet);
eval { require Parallel::ForkManager }
or die "Parallel::ForkManager is required when --parallel_genomes > 1.\n"
. "Install it with: cpanm Parallel::ForkManager\n";
}

die "No files to process\n" unless (@files);

Expand Down Expand Up @@ -488,12 +499,58 @@ sub process_file {

$file = $temp_file;

my $library_index = -1; # Make lists in the same order as @libraries
# Check bisulfite/aligner compatibility once before the loop (does not vary per library)
if ($bisulfite) {
die "The aligners BWA or minimap2 may not be used in --bisulfite mode, please adjust configuration.\n"
if (($aligner eq 'bwa') or ($aligner eq 'minimap2'));
}

# Set up Parallel::ForkManager when parallel_genomes > 1.
# When parallel_genomes == 1 (default), $pm is undef and the loop runs sequentially,
# preserving byte-for-byte identical output to the previous implementation.
my $pm;
if ($parallel_genomes > 1) {
$pm = Parallel::ForkManager->new($parallel_genomes);

# Merge per-library results returned by each child into the shared arrays.
# Each library owns a fixed, non-overlapping set of bits in @index_genomes
# determined by $library_index, so bitwise OR is safe regardless of child
# completion order.
$pm->run_on_finish(sub {
my ($pid, $exit_code, $ident, $exit_signal, $core_dump, $data) = @_;
return unless defined $data;
my ($lib_idx, $hits_ref, $bs_counts) = @$data;

# Merge per-read hit bitmask contributions via bitwise OR
while (my ($seqname, $val) = each %$hits_ref) {
$index_genomes[$seqname] = ($index_genomes[$seqname] // 0) | $val;
}

foreach my $library (@libraries) {
if ($bisulfite) {
$bisulfite_orientation[$lib_idx] = $bs_counts;
}
});
}

#Write Bowtie/Bowtie2 Standard Error to a temporary output file
#Generate a random filename and place in $outdir (if specified)
for (my $library_index = 0; $library_index < @libraries; $library_index++) {
my $library = $libraries[$library_index];

warn "Searching " . basename($file) . " against $library->[0]\n" unless ($quiet);

my $illumina_flag = '';
if ($illumina) {
$illumina_flag = '--phred64-quals';
}

if ($parallel_genomes > 1) {
$pm->start and next; # parent resumes loop; child continues below
}

# === This block runs in the child (parallel) or inline (sequential) ===

# Create per-invocation error temp file after the fork so each child
# gets its own unique file; tempfile() generates a random name which is
# safe even without forking.
my $error_fh;
my $error_filename;
if ($outdir) {
Expand All @@ -502,24 +559,34 @@ sub process_file {
( $error_fh, $error_filename ) = tempfile( 'aligner_standard_error.XXXXXXXX', SUFFIX => '.txt' );
}

warn "Searching " . basename($file) . " against $library->[0]\n" unless ($quiet);

my $illumina_flag = '';
if($illumina){
$illumina_flag = '--phred64-quals';
my ($hits_ref, $bs_counts);
if ($bisulfite) {
# Bismark writes output files prefixed with $library->[0] (the genome/database
# name), e.g. "Human.sample_bismark_bt2.bam". Because every library has a
# unique name in the configuration, parallel children for different libraries
# will never write to the same file.
($hits_ref, $bs_counts) = bisulfite_mapping( $illumina_flag, $library, $file, $error_filename, $library_index, $error_fh );
} else {
$hits_ref = conventional_mapping( $illumina_flag, $read_length, $library, $file, $error_filename, $library_index, $error_fh );
}

#Count the index of the library being used
$library_index++;

if ($bisulfite) {
die "The aligners BWA or minimap2 may not be used in --bisulfite mode, please adjust configuration.\n" if(($aligner eq 'bwa') or ($aligner eq 'minimap2'));
bisulfite_mapping( $illumina_flag, $library, $file, $error_filename, $library_index, $error_fh, \@index_genomes, \@bisulfite_orientation );
if ($parallel_genomes > 1) {
$pm->finish(0, [$library_index, $hits_ref, $bs_counts]);
} else {
conventional_mapping( $illumina_flag, $read_length, $library, $file, $error_filename, \@index_genomes, $library_index, $error_fh );
# Sequential case: merge results immediately
while (my ($seqname, $val) = each %$hits_ref) {
$index_genomes[$seqname] = ($index_genomes[$seqname] // 0) | $val;
}
if ($bisulfite) {
$bisulfite_orientation[$library_index] = $bs_counts;
}
}
}

if ($parallel_genomes > 1) {
$pm->wait_all_children;
}

# Collate the hit results from the Bowtie searches.
# Result categories are:
# 0 - read not map to library
Expand Down Expand Up @@ -1119,11 +1186,17 @@ sub datestampGenerator {



#Uses Bismark to map an input file to a specified genome
#Results stored in the @index_genomes array
#Uses Bismark to map an input file to a specified genome.
#Returns a hashref of per-read hit values { seqname => bitmask_value } and
#an arrayref of bisulfite orientation counts [OT, CTOT, CTOB, OB].
#The bitmask values use only the two bits owned by $library_index; the caller
#merges them into the shared @index_genomes array via bitwise OR.
sub bisulfite_mapping {

my ( $illumina_flag, $library, $file, $error_filename, $library_index, $error_fh, $index_genomes_ref, $bisulfite_orientation_ref ) = @_;
my ( $illumina_flag, $library, $file, $error_filename, $library_index, $error_fh ) = @_;
my %local_hits;
my @bs_counts = (0, 0, 0, 0); # [OT, CTOT, CTOB, OB]

my $bismark_command;
my $sam_output_option = '';
$sam_output_option = '--sam' unless ( defined $path_to_samtools );
Expand Down Expand Up @@ -1186,10 +1259,8 @@ sub bisulfite_mapping {
next if ( substr( $line, 0, 1 ) eq '@' );

my ($seqname) = split( /\./, $line ); #Extract the index id
unless ( defined ${$index_genomes_ref}[$seqname] ) {
${$index_genomes_ref}[$seqname] = 0; #Initialise - array may have 'gaps'
}
${$index_genomes_ref}[$seqname] = record_hit( ${$index_genomes_ref}[$seqname], $library_index + 1 );
$local_hits{$seqname} //= 0;
$local_hits{$seqname} = record_hit( $local_hits{$seqname}, $library_index + 1 );

#Now determine the directionality of the read
my $read_conversion;
Expand All @@ -1208,13 +1279,13 @@ sub bisulfite_mapping {
}

if ($read_conversion eq 'CT' and $genome_conversion eq 'CT') { ## original top strand OT
${ $bisulfite_orientation_ref }[$library_index]->[0]++;
$bs_counts[0]++;
} elsif ($read_conversion eq 'GA' and $genome_conversion eq 'CT') { ## complementary to original top strand CTOT
${ $bisulfite_orientation_ref }[$library_index]->[1]++;
$bs_counts[1]++;
} elsif ($read_conversion eq 'GA' and $genome_conversion eq 'GA') { ## complementary to original bottom strand CTOB
${ $bisulfite_orientation_ref }[$library_index]->[2]++;
$bs_counts[2]++;
} elsif ($read_conversion eq 'CT' and $genome_conversion eq 'GA') { ## original bottom strand OB
${ $bisulfite_orientation_ref }[$library_index]->[3]++;
$bs_counts[3]++;
} else {
die "Unexpected combination of read and genome conversion: '$read_conversion' / '$genome_conversion'\n";
}
Expand All @@ -1232,12 +1303,9 @@ sub bisulfite_mapping {
$seqname = substr($seqname, 1); #Ignore @

#Record twice, so reads from the ambiguous file are considered as multi-mapping
unless ( defined ${$index_genomes_ref}[$seqname] ) {
${$index_genomes_ref}[$seqname] = 0; #Initialise - array may have 'gaps'
}

${$index_genomes_ref}[$seqname] = record_hit( ${$index_genomes_ref}[$seqname], $library_index + 1 );
${$index_genomes_ref}[$seqname] = record_hit( ${$index_genomes_ref}[$seqname], $library_index + 1 );
$local_hits{$seqname} //= 0;
$local_hits{$seqname} = record_hit( $local_hits{$seqname}, $library_index + 1 );
$local_hits{$seqname} = record_hit( $local_hits{$seqname}, $library_index + 1 );

scalar <AMBIGUOUS_FILE>; #Ignore rest of FASTQ read
scalar <AMBIGUOUS_FILE>;
Expand All @@ -1259,14 +1327,20 @@ sub bisulfite_mapping {
my $bismark_report_file = $mapped_file;
$bismark_report_file =~ s/\.sam$|\.bam$/_SE_report.txt/;
unlink $bismark_report_file or die "Could not delete Bismark report file '$bismark_report_file'.\n";
unlink $ambiguous_file or die "Could not delete Bismark amibiguous reads outputfile '$ambiguous_file'.\n";
unlink $ambiguous_file or die "Could not delete Bismark ambiguous reads output file '$ambiguous_file'.\n";

return (\%local_hits, \@bs_counts);
}

#Uses Bowtie or Bowtie2 to map an input file to a specified genome
#Results stored in the @index_genomes array
#Uses Bowtie or Bowtie2 to map an input file to a specified genome.
#Returns a hashref of per-read hit values { seqname => bitmask_value } where the
#bitmask values use only the two bits owned by $library_index. The caller merges
#them into the shared @index_genomes array via bitwise OR.
sub conventional_mapping {

my ( $illumina_flag, $read_length, $library, $file, $error_filename, $index_genomes_ref, $library_index, $error_fh ) = @_;
my ( $illumina_flag, $read_length, $library, $file, $error_filename, $library_index, $error_fh ) = @_;
my %local_hits;
my $aligner_command;

#Determine whether to execute bowtie1 or bowtie2
Expand Down Expand Up @@ -1297,10 +1371,8 @@ sub conventional_mapping {
($seqname) = split(/\./, $seqname);
next if $samFlag & 0x4; #Did not align

unless ( defined ${$index_genomes_ref}[$seqname] ) {
${$index_genomes_ref}[$seqname] = 0; #Initialise - array may have 'gaps'
}
${$index_genomes_ref}[$seqname] = record_hit( ${$index_genomes_ref}[$seqname], $library_index + 1 );
$local_hits{$seqname} //= 0;
$local_hits{$seqname} = record_hit( $local_hits{$seqname}, $library_index + 1 );
}

#Check the Standard Error output file and report any errors
Expand All @@ -1315,6 +1387,8 @@ sub conventional_mapping {
}
close $error_fh;
unlink $error_filename or die "Could not delete temporary Standard Error file '$error_filename' : $!";

return \%local_hits;
}

sub remove_duplicates {
Expand Down Expand Up @@ -2676,6 +2750,14 @@ Options
If no directory is specified then output files
are saved in the current working directory.

--parallel_genomes <int>
Align against this many reference genomes in parallel.
Defaults to 1 (sequential, identical to previous
behaviour). When set above 1, total CPU usage is
approximately parallel_genomes * threads cores.
Requires the Perl module Parallel::ForkManager when
set above 1 (install with: cpanm Parallel::ForkManager).

--pass <int> Used in conjunction with --filter. By default all
genome filters must be passed for a read to pass
the --filter option. However, a minimum number
Expand Down