Mercurial > repos > dereeper > pangenome_explorer
view COG/bac-genomics-scripts/rod_finder/blast_rod_finder.pl @ 15:dbde253606c5 draft default tip
Uploaded
author | dereeper |
---|---|
date | Wed, 11 Dec 2024 08:25:06 +0000 |
parents | e42d30da7a74 |
children |
line wrap: on
line source
#!/usr/bin/perl ####### # POD # ####### =pod =head1 NAME blast_rod_finder.pl 07-11-2011 =head1 SYNOPSIS C<perl blast_rod_finder.pl -q query.embl -r blastn.out -m 2000> =head1 DESCRIPTION This script is intended to identify region of differences (RODs) between a nucleotide query and a nucleotide subject/reference sequence. In order to do so, a BLASTN (L<http://blast.ncbi.nlm.nih.gov/Blast.cgi>) needs to be performed beforehand with the query and the subject sequences. I<blast_rod_finder.pl> is mainly designed to work with bacterial genomes, while a query genome can be blasted against several subject sequences to detect RODs over a number of references. Although the results are optimized towards a complete query genome, both the reference(s) as well as the query can be used in draft form. To create artificial genomes use I<cat_seq.pl> or the EMBOSS application union (L<http://emboss.sourceforge.net/>). The BLASTN report file, the query sequence file (preferably in embl or genbank format with annotation, see below) and a minimum size for ROD detection have to be provided for I<blast_rod_finder.pl> to run. Subsequently, RODs are summarized in a tab-separated result file, a gff3 (usable e.g. in Artemis/DNAPlotter, L<http://www.sanger.ac.uk/resources/software/artemis/>) and a BRIG (BLAST Ring Image Generator, L<http://brig.sourceforge.net/>) output file. Nucleotide sequences of each ROD are written to a multi-fasta file. The query sequence can be provided in embl or genbank format, but has to correspond to the fasta file used in querying the BLAST database (the accession numbers have to correspond to the fasta headers). Use I<seq_format-converter.pl> to create a corresponding fasta file from embl|genbank files for BLASTN if needed. With annotated query files additional info is given in the result summary and the amino acid sequences of all non-pseudo CDSs, which are contained or overlap a ROD, are written to a result file. Furthermore, all detected RODs are saved in individual sequence files in the corresponding query sequence format. Run BLASTN and the script I<blast_rod_finder.pl> manually or use the bash shell scripts I<blast_rod_finder*.sh> (see examples below) to perform the pipeline consecutively in one folder. The Perl script runs on BioPerl (L<http://www.bioperl.org>). =head1 OPTIONS =head2 Mandatory options =over 21 =item B<-m>=I<int>, B<-min>=I<int> Minimum size of RODs that are reported =item B<-q>=I<str>, B<-query>=I<str> Query sequence file [fasta, embl, or genbank format] =item B<-r>=I<str>, B<-report>=I<str> BLASTN report/output file =back =head2 Optional options =over =item B<-h>, B<-help> Help (perldoc POD) =back =head1 OUTPUT =over 17 =item F<./results> All output files are stored in this result folder =item F<rod_summary.txt> Summary of detected ROD regions (for embl/genbank queries includes annotation), tab-separated =item F<rod.gff> GFF3 file with ROD coordinates to use in Artemis/DNAPlotter etc. =item F<rod_BRIG.txt> ROD coordinates to use in BRIG (BLAST Ring Image Generator), tab-separated =item F<rod_seq.fasta> Nucleotide sequences of ROD regions (>ROD# size start..stop), multi-fasta =item (F<rod_aa_fasta.txt>) Only present if query is given with annotation, i.e. embl|genbank format Amino acid sequences of all CDSs that are contained in or overlap a ROD region in multi-fasta format (>locus_tag gene product). RODs are seperated in the file via '~~' (~~ROD# size start..stop). =item (F<ROD#.embl|gbk>) Only present if query is given with annotation, i.e. embl|genbank format. Each identified ROD is written to an individual sequence file (in the same format as the query). =back =head1 EXAMPLES =head2 Legacy BLASTN =over =item C<formatdb -p F -i subject.fasta -n blast_db> =item C<blastall -p blastn -d blast_db -i query.fasta -o blastn.out -e 2e-11 -F F> =back =head2 BLASTN+ =over =item C<makeblastdb -in subject.fasta -input_type fasta -dbtype nucl -out blast_db> =item C<blastn -db blast_db -query query.fasta -out blastn.out -evalue 2e-11 -soft_masking false> =back =head2 blast_rod_finder.pl =over =item C<perl blast_rod_finder.pl -q query.embl|gbk|fasta -r blastn.out -m 5000> =back =head2 All-in-one with unix bash-shell scripts =over =item C<./blast_rod_finder_legacy.sh subject.fasta query.fasta query.embl|gbk|fasta 5000> =item C<./blast_rod_finder_plus.sh subject.fasta query.fasta query.embl|gbk|fasta 5000> =back =head1 VERSION 0.4 update: 13-02-2013 =head1 AUTHORS Andreas Leimbach aleimba[at]gmx[dot]de David Studholme D[dot]J[dot]Studholme[at]exeter[dot]ac[dot]uk =head1 LICENSE This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3 (GPLv3) of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see L<http://www.gnu.org/licenses/>. =cut ######## # MAIN # ######## use strict; use warnings; use Getopt::Long; # module to get options from the command line use Bio::SearchIO; # bioperl module to handle blast reports use Bio::SeqIO; # bioperl module to handle sequence input/output use Bio::Seq; # bioperl module to handle sequences with features use Bio::SeqFeatureI; # bioperl module to handle features in a sequence use Bio::SeqUtils; # bioperl module with additional methods (including features) for Bio::Seq objects (e.g. revcom, truncate, concatenate) ### Get the options with Getopt::Long, works also abbreviated and with two "--": -r, --r, -report ... my $usage = "perldoc $0"; my $blast_report = ''; # blastn report/output file my $q_file = ''; # query sequence file to get the RODs and possible including annotations my $minimum_size = ''; # minimum size of uncovered regions (RODs) that are reported my $help = ''; # run perldoc on the POD GetOptions ('report=s' => \$blast_report, 'query=s' => \$q_file, 'min=s' => \$minimum_size, 'help' => \$help); ### Run perldoc on POD if -h|--h|--help is given as argument or arguments are not given if (!$blast_report || !$q_file || !$minimum_size) { die system($usage); } elsif ($help) { die system($usage); } ### Parse the blast report/output file my %q_covered; # stores for each seq position (key) of a query accession if it falls within a hsp hit (value set to 1) my @blast_q_acc; # store query accessions to test if identical in blast-parse and query-seq (seqio below) my $blast_db; # add to GFF3 print out my $parser = new Bio::SearchIO(-file => "<$blast_report", -format => 'blast'); # Bio::SearchIO object $|++; # turn on autoflush, forces STDOUT to flush right away print "\nParsing blast report file."; # status display (with autoflush) while(my $result = $parser->next_result) { # several query sequences possible (result = entire analysis for a single query seq) my $query_acc = $result->query_accession; $query_acc =~ s/\.$//; # rm a '.' if present at the end of the string (for non-NCBI fasta headers) $query_acc =~ s/\.\d$//; # rm version nr. from NCBI query acc (to fit it to the acc.-nr. from Bio::SeqIO below) push(@blast_q_acc, $query_acc); $blast_db = $result->database_name; chop $blast_db; # always an empty character after the name? while (my $hit = $result->next_hit) { # several subject sequences in the database might have hits to a single query print '.'; # status display while(my $hsp = $hit->next_hsp) { # each hit might have one or more hsps (the alignments shown in a blast report) my ($query_start, $query_end) = $hsp->range('query'); # query hsp start/stop coords foreach my $i ($query_start .. $query_end) { $q_covered{$query_acc}{$i}++; # changes for each hsp hit the value from 'undef' to '1' } } } } ### Read the query sequence file into RAM, for Bio::Seq::RichSeq files (e.g. embl/genbank) also read features my %q_seqobj; # combine all seqobj (values; for multi-fasta/embl/gbk queries) to acc.-nr.s (keys ) my %features; # stores all features of a query RichSeq file my $multi = 0; # test if multi-seq query (multi-fasta/embl/genbank) >= 2; and counter my $seqio_obj = Bio::SeqIO->new(-file => "<$q_file"); # no '-format' to leave it to bioperl guessing print "\nParsing query sequence file and potential annotation."; # status display (with autoflush) while (my $seq_obj = $seqio_obj->next_seq) { # Bio::Seq object, query might be multi-fasta/embl/genbank print '.'; # status display my $query_acc; if (ref($seqio_obj) =~ m/\:\:fasta$/i) { # $seqio_obj blessed, thus ref returns package name (here= Bio::SeqIO::fasta) from bioperl guessing $query_acc = $seq_obj->display_id; # method '->accession_number' doesn't work with fasta files $query_acc =~ s/gi\|\d+\|(emb|gb|dbj|ref)\|(.+)\|/$2/; # get acc. nr. from NCBI fasta ID lines $query_acc =~ s/\.\d$//; # rm version nr., to fit it to blast acc. nr. above } else { # Bio::Seq::RichSeq file (embl, genbank) inherited from Bio::SeqIO::embl/genbank $query_acc = $seq_obj->accession_number; @{$features{$query_acc}} = $seq_obj->get_all_SeqFeatures; # store all Bio::SeqFeatureI objects in anonymous array } if ($query_acc ne $blast_q_acc[$multi]) { # check if acc.s fit; SAME order of acc.s in blast-report and seq-obj die "\n\n###Fatal error:\nBlast query accession \'$blast_q_acc[$multi]\' doesn't fit to query sequence accession \'$query_acc\'!\n\n"; } $q_seqobj{$query_acc} = $seq_obj; $multi++; } ### Get the contiguous unaligned regions of the query seq(s) my %uncovered_regions; # stores coords of uncovered regions (RODs) print "\nLooking for RODs."; # status display (with autoflush) foreach my $acc (keys %q_seqobj) { print "."; # status display my $start; my $seq = $q_seqobj{$acc}->seq; # Bio::SeqIO object method 'seq' my $previous_state = 1; # begin at base #1 of query seq foreach my $i (1 .. length($seq)) { # loop through the whole query seq and look for undef positions my $current_state = $q_covered{$acc}{$i}; if (!defined $current_state and defined $previous_state) { # warn "We have just entered an unaligned region $acc: $i\n"; $start = $i; } elsif (defined $current_state and !defined $previous_state) { # warn "We have just exited an unaligned region $acc: ", $i-1, "\n"; my $end = $i - 1; # position before current_state was undef my $length = $end - $start + 1; push @{$uncovered_regions{$length}{$acc}}, $start; # anonymous array in hash-in-hash data structure undef $start; } elsif ($i == length($seq) and !defined $previous_state) { # warn "We have just exited an unaligned region at end of contig $acc: $i\n"; my $end = $i; my $length = $end - $start + 1; push @{$uncovered_regions{$length}{$acc}}, $start; undef $start; } $previous_state = $current_state; } } ### Create results directory, where output files are written to my $out_dir = './results/'; if (-e $out_dir) { print "\n\n###Directory \'$out_dir\' already exists! Replace the directory and all its contents? [y|n] "; my $user_ask = <STDIN>; if ($user_ask =~ /y/i) { unlink glob $out_dir."*"; # remove all files in results directory rmdir $out_dir; # remove the empty directory } else { die "Script abborted!\n"; } } mkdir $out_dir or die "Can't create directory \"$out_dir\": $!\n"; ### List the longest uncovered regions and print results in output files my $rod_result = 'rod_summary.txt'; # overview of parsed ROD results and possible ROD annotations open(ROD, ">$out_dir"."$rod_result") or die "Failed to create $rod_result file: $!\n"; my $seq_out = 'rod_seq.fasta'; # multi-fasta nucleotide file of ROD seqs open(SEQ, ">$out_dir"."$seq_out") or die "Failed to create $seq_out file: $!\n"; my $gff_out = 'rod.gff'; # GFF3 file of ROD regions for use in artemis/DNAPlotter etc. open(GFF, ">$out_dir"."$gff_out") or die "Failed to create $gff_out file: $!\n"; my $brig_out = 'rod_BRIG.txt'; # tab-separated ROD file to load into BRIG open (BRIG, ">$out_dir"."$brig_out") or die "Failed to create $brig_out file: $!\n"; print ROD "$q_file regions with no blast hits (RODs) in $blast_report\n"; print ROD "Rank\t"; if ($multi >=2) { # for multi-seq query files include acc. nr. column print ROD "Accession\t"; } print ROD "Length\tROD position"; my $cds_out; # only for RichSeq query files (embl/genbank) if (ref($seqio_obj) !~ m/\:\:fasta$/i) { # query NOT a fasta file print ROD "\tFeature primary_tag\tFeature locus_tag\tgene\tFeature product\tFeature position"; $cds_out = 'rod_aa_fasta.txt'; # CDS amino acid sequence output file open(CDS, ">$out_dir"."$cds_out") or die "Failed to create $cds_out file: $!\n"; } print ROD "\n"; print GFF "##gff-version 3\n"; # header, shows it's a gff version 3 file print GFF "#$q_file regions with no blast hits (RODs) in $blast_report\n"; # comment line for description print BRIG "#Start\tStop\tLabel\n"; # column headers (commented) print "\nPreparing results."; # status display (with autoflush) my $rank = 0; # number of ROD regions foreach my $length (sort {$b<=>$a} keys %uncovered_regions) { # sort RODs from large length to small length if ($length >= $minimum_size) { print "."; # status display foreach my $acc (sort keys %{$uncovered_regions{$length}}) { # de-reference hash-in-hash, sort by acc. nr. foreach my $start (sort @{$uncovered_regions{$length}{$acc}}) { # de-reference anonymous array my $end = $start + $length - 1; my $pos = "$start..$end"; $rank++; print GFF "$acc\tBLASTN\tsequence_difference\t$start\t$end\t.\t+\t.\tName=ROD$rank;Target=$blast_db;color=2\n"; print BRIG "$start\t$end\tROD$rank\n"; print ROD "$rank\t"; if ($multi >=2) { # for multi-seq query files include acc. nr. print ROD "$acc\t"; } print ROD "$length\t$pos"; print SEQ ">ROD$rank $length $pos\n"; print SEQ $q_seqobj{$acc}->subseq($start,$end), "\n\n"; # subseq method of Bio::Seq object ### fill ROD feature columns in 'rod_result_summary.txt', print CDS aa seqs to 'rod_aa_fasta.txt' for RichSeq files, ### and print each ROD as a separate sequence file (with Bio::Sequtils->'trunc_with_features') if (ref($seqio_obj) !~ m/\:\:fasta$/i) { # RichSeq query file if ($rank >1) { # blank line in front of the next ROD aa seq block, if it's not the first block print CDS "\n"; } my $truncseq = Bio::SeqUtils->trunc_with_features($q_seqobj{$acc}, $start, $end); # truncate current Bio::Seq object to coordinates my $seqout; # variable to hold Bio::SeqIO object if (ref($seqio_obj) =~ /genbank/) { # print the truncated ROD sequences in the same '-format' as $q_file; Bio::SeqIO::genbank object $seqout = Bio::SeqIO->new(-file => ">$out_dir"."ROD$rank".".gbk", -format => 'genbank'); } elsif (ref($seqio_obj) =~ /embl/) { # Bio::SeqIO::embl object $seqout = Bio::SeqIO->new(-file => ">$out_dir"."ROD$rank".".embl", -format => 'embl'); } $seqout->write_seq($truncseq); # write truncated sequence object my $loop = 0; # indicates that more than one feature in the ROD range print CDS "~~ROD$rank $length $pos\n"; foreach my $feature (@{$features{$acc}}) { # de-reference anonymous array in hash if ($feature->location->start >= $start && $feature->location->end <= $end) { # features that are fully within ROD region $loop = print_CDSs($feature, $loop); # subroutine 'print_CDSs' } elsif (($feature->location->start <= $start && ($feature->location->end > $start && $feature->location->end <= $end)) || (($feature->location->start >= $start && $feature->location->start < $end) && $feature->location->end > $end)) { # features that overlap the ROD region $loop = print_CDSs($feature, $loop, 1); # 1 defines $overlap in subroutine } } if ($loop == 0) { # print "\n" if no features exist in the ROD print ROD "\n"; } } else { # fasta files don't contain feature information print ROD "\n"; } } } } } close ROD; close SEQ; close GFF; close BRIG; if (ref($seqio_obj) !~ m/\:\:fasta$/i) { close CDS; } $| = 0; # turn off autoflush ### Print which files have been created! print "\n\nThe following files were created in the \"./results\" directory:\n"; print "- $rod_result: Tab-separated summary of the found regions of difference (RODs)\n"; print "- $gff_out: GFF3 file with ROD coordinates to use in Artemis/DNAPlotter ...\n"; print "- $brig_out: Tab-separated file for BRIG (BLAST Ring Image Generator)\n"; print "- $seq_out: Nucleotide sequences of ROD regions in multi-fasta format\n"; if (!(ref($seqio_obj) =~ m/\:\:fasta$/i)) { print "- $cds_out: Amino acid sequences of CDSs in ROD regions\n"; print "- ROD#.format: Each ROD as a separate sequence file in the query file format\n"; } print "\n"; exit; ############### # Subroutines # ############### ### Subroutine to print the position of features according to leading and lagging strand sub position { my $feature = shift; if ($feature->strand == 1) { # leading strand print ROD "\t", $feature->location->start, "..", $feature->location->end; } elsif ($feature->strand == -1) { # lagging strand print ROD "\t", $feature->location->end, "..", $feature->location->start; } return 1; } ### Subroutine that prints feature information from RichSeq files sub print_CDSs { my $feature = shift; my $loop = shift; my $overlap = shift; # defined if feature overlaps ROD edge my $primary_tag = $feature->primary_tag; if (($primary_tag eq 'CDS') || ($primary_tag eq 'misc_RNA') || ($primary_tag eq 'ncRNA') || ($primary_tag eq 'rRNA') || ($primary_tag eq 'tRNA') || ($primary_tag eq 'tmRNA') || ($primary_tag eq 'misc_feature') || ($primary_tag eq 'repeat_region')) { skip_columns($loop); # sub to skip to correct column in print out print ROD "\t$primary_tag"; if ($feature->has_tag('locus_tag')) { # TRUE if tag 'locus_tag' exists for this feature print ROD "\t", $feature->get_tag_values('locus_tag'); print CDS ">", $feature->get_tag_values('locus_tag') if ($primary_tag eq 'CDS' && !$feature->has_tag('pseudo')); # include only CDSs with "/translation" in file 'rod_aa_fasta.txt' (exclude pseudo-genes, RNAs ...) } else { print ROD "\t"; print CDS ">" if ($primary_tag eq 'CDS' && !$feature->has_tag('pseudo')); } if ($feature->has_tag('gene')) { print ROD "\t", $feature->get_tag_values('gene'); print CDS " ", $feature->get_tag_values('gene') if ($primary_tag eq 'CDS' && !$feature->has_tag('pseudo')); } else { print ROD "\t"; print CDS " " if ($primary_tag eq 'CDS' && !$feature->has_tag('pseudo')); } if ($feature->has_tag('product')) { my ($product) = $feature->get_tag_values('product'); print ROD "\t", $product; print CDS " ", replace_white($product) if ($primary_tag eq 'CDS' && !$feature->has_tag('pseudo')); # subroutine 'replace_white' } elsif ($feature->has_tag('rpt_family')) { # don't include feature 'rpt_family' in the aa output print ROD "\t", $feature->get_tag_values('rpt_family'); } elsif ($feature->has_tag('note')) { # only include feature 'note' if 'product' doesn't exist my ($note) = $feature->get_tag_values('note'); print ROD "\t", $note; print CDS " ", replace_white($note) if ($primary_tag eq 'CDS' && !$feature->has_tag('pseudo')); } position($feature); # subroutine 'position' to print coordinates of features if (defined $overlap) { # include '(overlap)' if feature overlaps edge of ROD print ROD " (overlap)\n"; print CDS " (overlap)\n" if ($primary_tag eq 'CDS' && !$feature->has_tag('pseudo')); } else { print ROD "\n"; print CDS "\n" if ($primary_tag eq 'CDS' && !$feature->has_tag('pseudo')); } if ($primary_tag eq 'CDS' && $feature->has_tag('translation')) { # print translations of CDSs print CDS $feature->get_tag_values('translation'), "\n"; } $loop++; # Indicates that more than one feature in the ROD range } elsif ($primary_tag eq 'mobile_element') { # include 'mobile_element' features only in ROD skip_columns($loop); print ROD "\t$primary_tag"; if ($feature->has_tag('mobile_element_type')) { print ROD "\t\t\t", $feature->get_tag_values('mobile_element_type'); } else { print ROD "\t\t\t"; } position($feature); # sub to print positions of features if (defined $overlap) { print ROD " (overlap)\n"; } else { print ROD "\n"; } $loop++; } return $loop; } ### Subroutine to replace whitespaces in a string sub replace_white { my $string = shift; $string =~ s/\s/_/g; return $string; } ### Subroutine to skip to the correct column for features in tab-separated result file 'rod_result_summary.txt' sub skip_columns { my $loop = shift; # $multi is in the same scope ... (actually also $loop) if ($loop >= 1) { # for the features after the first ROD line skip to the 'feature primary_tag' column in the tab-separated summary output print ROD "\t\t"; if ($multi >=2) { # multi-seq query files include acc. nr. column print ROD "\t"; } } return 1; }