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