kkonganti@1: #!/usr/bin/env perl kkonganti@1: kkonganti@1: # Kranti Konganti kkonganti@1: # 10/12/2022 kkonganti@1: kkonganti@1: use strict; kkonganti@1: use warnings; kkonganti@1: use Getopt::Long; kkonganti@1: use Data::Dumper; kkonganti@1: use Pod::Usage; kkonganti@1: use File::Basename; kkonganti@1: use File::Spec::Functions; kkonganti@1: kkonganti@1: my $tbl = {}; kkonganti@1: my $snp_2_serovar = {}; kkonganti@1: my $acc_2_serovar = {}; kkonganti@1: my $acc_2_target = {}; kkonganti@1: my $snp_count = {}; kkonganti@1: my $snp_2_acc = {}; kkonganti@1: my $acc_2_snp = {}; kkonganti@1: my $multi_cluster_acc = {}; kkonganti@1: my ( kkonganti@1: $serovar_limit, $serovar_or_type_col, $min_asm_size, kkonganti@1: $complete_serotype_name, $PDG_file, $table_file, kkonganti@1: $not_null_pdg_serovar, $snp_cluster, $help, kkonganti@1: $out_prefix kkonganti@1: ); kkonganti@1: my @custom_serovars; kkonganti@1: kkonganti@1: GetOptions( kkonganti@1: 'help' => \$help, kkonganti@1: 'pdg=s' => \$PDG_file, kkonganti@1: 'tbl=s' => \$table_file, kkonganti@1: 'snp=s' => \$snp_cluster, kkonganti@1: 'min_contig_size=i' => \$min_asm_size, kkonganti@1: 'complete_serotype_name' => \$complete_serotype_name, kkonganti@1: 'serocol:i' => \$serovar_or_type_col, kkonganti@1: 'not_null_pdg_serovar' => \$not_null_pdg_serovar, kkonganti@1: 'num_serotypes_per_serotype:i' => \$serovar_limit, kkonganti@1: 'include_serovar=s' => \@custom_serovars, kkonganti@1: 'op=s' => \$out_prefix kkonganti@1: ) or pod2usage( -verbose => 2 ); kkonganti@1: kkonganti@1: if ( defined $help ) { kkonganti@1: pod2usage( -verbose => 2 ); kkonganti@1: } kkonganti@1: kkonganti@1: if ( !defined $serovar_limit ) { kkonganti@1: $serovar_limit = 1; kkonganti@1: } kkonganti@1: kkonganti@1: if ( !defined $serovar_or_type_col ) { kkonganti@1: $serovar_or_type_col = 49; kkonganti@1: } kkonganti@1: kkonganti@1: if ( !defined $min_asm_size ) { kkonganti@1: $min_asm_size = 0; kkonganti@1: } kkonganti@1: kkonganti@1: if ( defined $out_prefix ) { kkonganti@1: $out_prefix .= '_'; kkonganti@1: } kkonganti@1: else { kkonganti@1: $out_prefix = ''; kkonganti@1: } kkonganti@1: kkonganti@1: pod2usage( -verbose => 2 ) if ( !$PDG_file || !$table_file || !$snp_cluster ); kkonganti@1: kkonganti@1: open( my $pdg_file, '<', $PDG_file ) kkonganti@1: || die "\nCannot open PDG file $PDG_file: $!\n\n"; kkonganti@1: open( my $tbl_file, '<', $table_file ) kkonganti@1: || die "\nCannot open tbl file $table_file: $!\n\n"; kkonganti@1: open( my $snp_cluster_file, '<', $snp_cluster ) kkonganti@1: || die "\nCannot open $snp_cluster: $!\n\n"; kkonganti@1: open( my $acc_fh, '>', 'acc2serovar.txt' ) kkonganti@1: || die "\nCannot open acc2serovar.txt: $!\n\n"; kkonganti@1: open( my $Stdout, '>&', STDOUT ) || die "\nCannot pipe to STDOUT: $!\n\n"; kkonganti@1: open( my $Stderr, '>&', STDERR ) || die "\nCannot pipe to STDERR: $!\n\n"; kkonganti@1: open( my $accs_snp_fh, '>', $out_prefix . 'accs_snp.txt' ) kkonganti@1: || die "\nCannnot open " . $out_prefix . "accs_snp.txt for writing: $!\n\n"; kkonganti@1: open( my $genome_headers_fh, '>', $out_prefix . 'mash_snp_genome_list.txt' ) kkonganti@1: || die "\nCannnot open " kkonganti@1: . $out_prefix kkonganti@1: . "mash_snp_genome_list.txt for writing: $!\n\n"; kkonganti@1: kkonganti@1: my $pdg_release = basename( $PDG_file, ".metadata.tsv" ); kkonganti@1: kkonganti@1: while ( my $line = <$pdg_file> ) { kkonganti@1: chomp $line; kkonganti@1: next if ( $line =~ m/^\#/ ); kkonganti@1: kkonganti@1: # Relevent columns (Perl index): kkonganti@1: # 9: asm_acc kkonganti@1: # 33: serovar kkonganti@1: # 48: computed serotype kkonganti@1: kkonganti@1: my @cols = split( /\t/, $line ); kkonganti@1: my $serovar_or_type = $cols[ $serovar_or_type_col - 1 ]; kkonganti@1: my $acc = $cols[9]; kkonganti@1: my $serovar = $cols[33]; kkonganti@1: my $target_acc = $cols[41]; kkonganti@1: kkonganti@1: $serovar_or_type =~ s/\"//g; kkonganti@1: kkonganti@1: my $skip = 1; kkonganti@1: foreach my $ser (@custom_serovars) { kkonganti@1: $skip = 0, next if ( $serovar_or_type =~ qr/\Q$ser\E/ ); kkonganti@1: } kkonganti@1: kkonganti@1: if ( defined $complete_serotype_name ) { kkonganti@1: next kkonganti@1: if ( $skip kkonganti@1: && ( $serovar_or_type =~ m/serotype=.*?\-.*?\,antigen_formula.+/ ) kkonganti@1: ); kkonganti@1: } kkonganti@1: kkonganti@1: next kkonganti@1: if ( kkonganti@1: $skip kkonganti@1: && ( $serovar_or_type =~ m/serotype=\-\s+\-\:\-\:\-/ kkonganti@1: || $serovar_or_type =~ m/antigen_formula=\-\:\-\:\-/ ) kkonganti@1: ); kkonganti@1: kkonganti@1: # next kkonganti@1: # if ( kkonganti@1: # ( kkonganti@1: # $serovar_or_type =~ m/serotype=\-\s+\-\:\-\:\-/ kkonganti@1: # || $serovar_or_type =~ m/antigen_formula=\-\:\-\:\-/ kkonganti@1: # ) kkonganti@1: # ); kkonganti@1: kkonganti@1: if ( defined $not_null_pdg_serovar ) { kkonganti@1: $acc_2_serovar->{$acc} = $serovar_or_type, kkonganti@1: $acc_2_target->{$acc} = $target_acc, kkonganti@1: print $acc_fh "$acc\t$serovar_or_type\n" kkonganti@1: if ( $acc !~ m/NULL/ kkonganti@1: && $serovar !~ m/NULL/ kkonganti@1: && $serovar_or_type !~ m/NULL/ ); kkonganti@1: } kkonganti@1: else { kkonganti@1: $acc_2_serovar->{$acc} = $serovar_or_type, kkonganti@1: $acc_2_target->{$acc} = $target_acc, kkonganti@1: print $acc_fh "$acc\t$serovar_or_type\n" kkonganti@1: if ( $acc !~ m/NULL/ && $serovar_or_type !~ m/NULL/ ); kkonganti@1: } kkonganti@1: kkonganti@1: # $snp_count->{$serovar_or_type} = 0; kkonganti@1: } kkonganti@1: kkonganti@1: # kkonganti@1: # SNP to ACC kkonganti@1: # kkonganti@1: kkonganti@1: while ( my $line = <$snp_cluster_file> ) { kkonganti@1: chomp $line; kkonganti@1: my @cols = split( /\t/, $line ); kkonganti@1: kkonganti@1: # Relevant columns kkonganti@1: # 0: SNP Cluster ID kkonganti@1: # 3: Genome Accession belonging to the cluster (RefSeq or GenBank) kkonganti@1: my $snp_clus_id = $cols[0]; kkonganti@1: my $acc = $cols[3]; kkonganti@1: kkonganti@1: next if ( $acc =~ m/^NULL/ || $snp_clus_id =~ m/^PDS_acc/ ); kkonganti@1: next if ( !exists $acc_2_serovar->{$acc} ); kkonganti@1: push @{ $snp_2_acc->{$snp_clus_id} }, $acc; kkonganti@1: if ( exists $acc_2_snp->{$acc} ) { kkonganti@1: print $Stderr kkonganti@1: "\nGot a duplicate assembly accession. Cannot proceed!\n\n$line\n\n"; kkonganti@1: exit 1; kkonganti@1: } kkonganti@1: $acc_2_snp->{$acc} = $snp_clus_id; kkonganti@1: $snp_count->{$snp_clus_id} = 0; kkonganti@1: } kkonganti@1: kkonganti@1: while ( my $line = <$tbl_file> ) { kkonganti@1: chomp $line; kkonganti@1: kkonganti@1: my @cols = split( /\t/, $line ); kkonganti@1: kkonganti@1: # .tbl file columns (Perl index): kkonganti@1: # kkonganti@1: # 0: Accession kkonganti@1: # 1: AssemblyLevel kkonganti@1: # 2: ScaffoldN50 kkonganti@1: # 3: ContigN50 kkonganti@1: kkonganti@1: my $acc = $cols[0]; kkonganti@1: my $asm_lvl = $cols[1]; kkonganti@1: my $scaffold_n50 = $cols[2]; kkonganti@1: my $contig_n50 = $cols[3]; kkonganti@1: kkonganti@1: # my $idx0 = $acc_2_serovar->{$cols[0]}; kkonganti@1: my $idx0 = $acc_2_snp->{$acc} if ( exists $acc_2_snp->{ $cols[0] } ); kkonganti@1: kkonganti@1: if ( not_empty($acc) && defined $idx0 ) { kkonganti@1: my $fna_rel_loc = kkonganti@1: "$pdg_release/ncbi_dataset/data/$acc/" kkonganti@1: . $acc kkonganti@1: . '_scaffolded_genomic.fna.gz'; kkonganti@1: kkonganti@1: if ( not_empty($scaffold_n50) ) { kkonganti@1: next if ( $scaffold_n50 <= $min_asm_size ); kkonganti@1: push @{ $snp_2_serovar->{$idx0}->{ sort_asm_level($asm_lvl) } kkonganti@1: ->{$scaffold_n50} }, "$acc_2_serovar->{$acc}|$fna_rel_loc"; kkonganti@1: } kkonganti@1: elsif ( not_empty($contig_n50) ) { kkonganti@1: next if ( $contig_n50 <= $min_asm_size ); kkonganti@1: push @{ $snp_2_serovar->{$idx0}->{ sort_asm_level($asm_lvl) } kkonganti@1: ->{$contig_n50} }, "$acc_2_serovar->{$acc}|$fna_rel_loc"; kkonganti@1: } kkonganti@1: } kkonganti@1: } kkonganti@1: kkonganti@1: foreach my $snp_cluster_id ( keys %$snp_2_acc ) { kkonganti@1: my $count = $snp_count->{$snp_cluster_id}; kkonganti@1: foreach my $asm_lvl ( kkonganti@1: sort { $a cmp $b } kkonganti@1: keys %{ $snp_2_serovar->{$snp_cluster_id} } kkonganti@1: ) kkonganti@1: { kkonganti@1: if ( $asm_lvl =~ m/Complete\s+Genome/i ) { kkonganti@1: $count = kkonganti@1: print_dl_metadata( $asm_lvl, kkonganti@1: \$snp_2_serovar->{$snp_cluster_id}->{$asm_lvl}, kkonganti@1: $count, $snp_cluster_id ); kkonganti@1: } kkonganti@1: if ( $asm_lvl =~ m/Chromosome/i ) { kkonganti@1: $count = kkonganti@1: print_dl_metadata( $asm_lvl, kkonganti@1: \$snp_2_serovar->{$snp_cluster_id}->{$asm_lvl}, kkonganti@1: $count, $snp_cluster_id ); kkonganti@1: } kkonganti@1: if ( $asm_lvl =~ m/Scaffold/i ) { kkonganti@1: $count = kkonganti@1: print_dl_metadata( $asm_lvl, kkonganti@1: \$snp_2_serovar->{$snp_cluster_id}->{$asm_lvl}, kkonganti@1: $count, $snp_cluster_id ); kkonganti@1: } kkonganti@1: if ( $asm_lvl =~ m/Contig/i ) { kkonganti@1: $count = kkonganti@1: print_dl_metadata( $asm_lvl, kkonganti@1: \$snp_2_serovar->{$snp_cluster_id}->{$asm_lvl}, kkonganti@1: $count, $snp_cluster_id ); kkonganti@1: } kkonganti@1: printf $Stderr "%-17s | %s\n", $snp_cluster_id, $count kkonganti@1: if ( $count > 0 ); kkonganti@1: last if ( $count == $serovar_limit ); kkonganti@1: } kkonganti@1: } kkonganti@1: kkonganti@1: close $pdg_file; kkonganti@1: close $tbl_file; kkonganti@1: close $snp_cluster_file; kkonganti@1: close $acc_fh; kkonganti@1: close $accs_snp_fh; kkonganti@1: kkonganti@1: #------------------------------------------- kkonganti@1: # Main ends kkonganti@1: #------------------------------------------- kkonganti@1: # Routines begin kkonganti@1: #------------------------------------------- kkonganti@1: kkonganti@1: sub print_dl_metadata { kkonganti@1: my $asm_lvl = shift; kkonganti@1: my $acc_sizes = shift; kkonganti@1: my $curr_count = shift; kkonganti@1: my $snp_cluster_id = shift; kkonganti@1: kkonganti@1: $asm_lvl =~ s/.+?\_(.+)/$1/; kkonganti@1: kkonganti@1: foreach my $acc_size ( sort { $b <=> $a } keys %{$$acc_sizes} ) { kkonganti@1: foreach my $serovar_url ( @{ $$acc_sizes->{$acc_size} } ) { kkonganti@1: my ( $serovar, $url ) = split( /\|/, $serovar_url ); kkonganti@1: return $curr_count if ( exists $multi_cluster_acc->{$url} ); kkonganti@1: $multi_cluster_acc->{$url} = 1; kkonganti@1: $curr_count++; kkonganti@1: my ( $final_acc, $genome_header ) = kkonganti@1: ( split( /\//, $url ) )[ 3 .. 4 ]; kkonganti@1: print $accs_snp_fh "$final_acc\n"; kkonganti@1: print $genome_headers_fh catfile( 'scaffold_genomes', kkonganti@1: $genome_header ) kkonganti@1: . "\n"; kkonganti@1: print $Stdout "$serovar|$asm_lvl|$acc_size|$url|$snp_cluster_id\n" kkonganti@1: if ( $curr_count > 0 ); kkonganti@1: last if ( $curr_count == $serovar_limit ); kkonganti@1: } kkonganti@1: last if ( $curr_count == $serovar_limit ); kkonganti@1: } kkonganti@1: return $curr_count; kkonganti@1: } kkonganti@1: kkonganti@1: sub sort_asm_level { kkonganti@1: my $level = shift; kkonganti@1: kkonganti@1: $level =~ s/(Complete\s+Genome)/a\_$1/ kkonganti@1: if ( $level =~ m/Complete\s+Genome/i ); kkonganti@1: $level =~ s/(Chromosome)/b\_$1/ if ( $level =~ m/Chromosome/i ); kkonganti@1: $level =~ s/(Scaffold)/c\_$1/ if ( $level =~ m/Scaffold/i ); kkonganti@1: $level =~ s/(Contig)/d\_$1/ if ( $level =~ m/Contig/i ); kkonganti@1: kkonganti@1: return $level; kkonganti@1: } kkonganti@1: kkonganti@1: sub not_empty { kkonganti@1: my $col = shift; kkonganti@1: kkonganti@1: if ( $col !~ m/^$/ ) { kkonganti@1: return 1; kkonganti@1: } kkonganti@1: else { kkonganti@1: return 0; kkonganti@1: } kkonganti@1: } kkonganti@1: kkonganti@1: __END__ kkonganti@1: kkonganti@1: =head1 SYNOPSIS kkonganti@1: kkonganti@1: This script will take in a PDG metadata file, a C<.tbl> file and generate kkonganti@1: the final list by B> priority. kkonganti@1: kkonganti@1: See complete description: kkonganti@1: kkonganti@1: perldoc waterfall_per_snp_cluster.pl kkonganti@1: kkonganti@1: or kkonganti@1: kkonganti@1: waterfall_per_snp_cluster.pl --help kkonganti@1: kkonganti@1: Examples: kkonganti@1: kkonganti@1: waterfall_per_snp_cluster.pl kkonganti@1: kkonganti@1: =head1 DESCRIPTION kkonganti@1: kkonganti@1: We will retain up to N number of genome accessions per SNP cluster. kkonganti@1: It prioritizes SNP Cluster participation over serotype coverage. kkonganti@1: Which N genomes are selected depends on (in order): kkonganti@1: kkonganti@1: 1. Genome assembly level, whose priority is kkonganti@1: kkonganti@1: a: Complete Genome kkonganti@1: b: Chromosome kkonganti@1: c: Scaffold kkonganti@1: d: Contig kkonganti@1: kkonganti@1: 2. If the genomes are of same assembly level, then kkonganti@1: scaffold N50 followed by contig N50 is chosen. kkonganti@1: kkonganti@1: 3. If the scaffold or contig N50 is same, then all kkonganti@1: of them are included kkonganti@1: kkonganti@1: =head1 OPTIONS kkonganti@1: kkonganti@1: =over 3 kkonganti@1: kkonganti@1: =item -p PDGXXXXX.XXXX.metadata.tsv kkonganti@1: kkonganti@1: Absolute UNIX path pointing to the PDG metadata file. kkonganti@1: Example: PDG000000002.2505.metadata.tsv kkonganti@1: kkonganti@1: =item -t asm.tbl kkonganti@1: kkonganti@1: Absolute UNIX path pointing to the file from the result kkonganti@1: of the C script, which is the C kkonganti@1: file. kkonganti@1: kkonganti@1: =item -snp PDGXXXXXXX.XXXX.reference_target.cluster_list.tsv kkonganti@1: kkonganti@1: Absolute UNIX path pointing to the SNP Cluster metadata file. kkonganti@1: Examples: PDG000000002.2505.reference_target.cluster_list.tsv kkonganti@1: kkonganti@1: kkonganti@1: =item --serocol (Optional) kkonganti@1: kkonganti@1: Column number (non 0-based index) of the PDG metadata file kkonganti@1: by which the serotypes are collected. Default: 49 kkonganti@1: kkonganti@1: =item --complete_serotype_name (Optional) kkonganti@1: kkonganti@1: Skip indexing serotypes when the serotype name in the column kkonganti@1: number 49 (non 0-based) of PDG metadata file consists a "-". For example, if kkonganti@1: an accession has a I> string as such in column kkonganti@1: number 49 (non 0-based): C<"serotype=- 13:z4,z23:-","antigen_formula=13:z4,z23:-"> kkonganti@1: then, the indexing of that accession is skipped. kkonganti@1: Default: False kkonganti@1: kkonganti@1: =item --not_null_pdg_serovar (Optional) kkonganti@1: kkonganti@1: Only index the B> column i.e. column number 49 (non 0-based) kkonganti@1: if the B> column is not C. kkonganti@1: kkonganti@1: =item -i (Optional) kkonganti@1: kkonganti@1: Make sure the following serotype is included. Mention C<-i> multiple kkonganti@1: times to include multiple serotypes. kkonganti@1: kkonganti@1: =item -num (Optional) kkonganti@1: kkonganti@1: Number of genome accessions per SNP Cluster. Default: 1 kkonganti@1: kkonganti@1: =back kkonganti@1: kkonganti@1: =head1 AUTHOR kkonganti@1: kkonganti@1: Kranti Konganti kkonganti@1: kkonganti@1: =cut