Mercurial > repos > galaxyp > mqppep_anova
diff PhosphoPeptide_Upstream_Kinase_Mapping.pl @ 0:dbff53e6f75f draft
planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/mqppep commit 3a7b3609d6e514c9e8f980ecb684960c6b2252fe
| author | galaxyp |
|---|---|
| date | Mon, 11 Jul 2022 19:22:25 +0000 |
| parents | |
| children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/PhosphoPeptide_Upstream_Kinase_Mapping.pl Mon Jul 11 19:22:25 2022 +0000 @@ -0,0 +1,2192 @@ +#!/usr/local/bin/perl +############################################################################################################################### +# perl Kinase_enrichment_analysis_complete_v0.pl +# +# Nick Graham, USC +# 2016-02-27 +# +# Built from scripts written by NG at UCLA in Tom Graeber's lab: +# CombinePhosphoSites.pl +# Retrieve_p_motifs.pl +# NetworKIN_Motif_Finder_v7.pl +# +# Given a list of phospho-peptides, find protein information and upstream kinases. +# Output file can be used for KS enrichment score calculations using Enrichment_Score4Directory.pl +# +# Updated 2022-01-13, Art Eschenlauer, UMN on behalf of Justin Drake's lab: +# Added warnings and used strict; +# fixed some code paths resulting in more NetworKIN matches; +# applied Aho-Corasick algorithm (via external Python script because Perl implementation was still too slow) +# to speed up "Match the non_p_peptides to the @sequences array"; +# added support for SQLite-formatted UniProtKB/Swiss-Prot data as an alternative to FASTA-formatted data; +# added support for SQLite output in addition to tabular files. +# +# +############################################################################################################################### + +use strict; +use warnings 'FATAL' => 'all'; + +use Getopt::Std; +use DBD::SQLite::Constants qw/:file_open/; +use DBI qw(:sql_types); +use File::Copy; +use File::Basename; +use POSIX qw(strftime); +use Time::HiRes qw(gettimeofday); +#use Data::Dump qw(dump); + +my $USE_SEARCH_PPEP_PY = 1; +#my $FAILED_MATCH_SEQ = "Failed match"; +my $FAILED_MATCH_SEQ = 'No Sequence'; +my $FAILED_MATCH_GENE_NAME = 'No_Gene_Name'; + +my $dirname = dirname(__FILE__); +my %opts; +my ($file_in, $average_or_sum, $db_out, $file_out, $file_melt, $phospho_type); +my $dbtype; +my ($fasta_in, $networkin_in, $motifs_in, $PSP_Kinase_Substrate_in, $PSP_Regulatory_Sites_in); +my (@samples, %sample_id_lut, %ppep_id_lut, %data, @tmp_data, %n); +my $line = 0; +my @failed_match = ($FAILED_MATCH_SEQ); +my @failed_matches; +my (%all_data); +my (@p_peptides, @non_p_peptides); +my @parsed_fasta; +my (@accessions, @names, @sequences, @databases, $database); +my ($dbfile, $dbh, $stmth); +my @col_names; +my (%matched_sequences, %accessions, %names, %sites, ); +my (@tmp_matches, @tmp_accessions, @tmp_names, @tmp_sites); +my (%p_residues, @tmp_p_residues, @p_sites, $left, $right, %p_motifs, @tmp_motifs_array, $tmp_motif, $tmp_site, %residues); +my (@kinases_observed, $kinases); +my (@kinases_observed_lbl, @phosphosites_observed_lbl); +my ($p_sequence_kinase, $p_sequence, $kinase); +my (@motif_sequence, @motif_description, @motif_type_key_ary, %motif_type, %motif_count); +my (@kinases_PhosphoSite, $kinases_PhosphoSite); +my ($p_sequence_kinase_PhosphoSite, $p_sequence_PhosphoSite, $kinase_PhosphoSite); +my (%regulatory_sites_PhosphoSite_hash); +my (%domain, %ON_FUNCTION, %ON_PROCESS, %ON_PROT_INTERACT, %ON_OTHER_INTERACT, %notes, %organism); +my (%unique_motifs); +my ($kinase_substrate_NetworKIN_matches, $kinase_substrate_PhosphoSite_matches); +my %psp_regsite_protein_2; +my (%domain_2, %ON_FUNCTION_2, %ON_PROCESS_2, %ON_PROT_INTERACT_2, %N_PROT_INTERACT, %ON_OTHER_INTERACT_2, %notes_2, %organism_2); +my @timeData; +my $PhosphoSitePlusCitation; +my (%site_description, %site_id); + +my %kinase_substrate_NetworKIN_matches; +my %kinase_motif_matches; +my $regulatory_sites_PhosphoSite; +my ($seq_plus5aa, $seq_plus7aa, %seq_plus7aa_2); +my %kinase_substrate_PhosphoSite_matches; +my @formatted_sequence; +my $pSTY_sequence; +my $i; +my @a; +my $use_sqlite; +my $verbose; + +########## +## opts ## +########## + ## input files + # i : path to input file, e.g., 'outputfile_STEP2.txt' + # f : path to UniProtKB/SwissProt FASTA + # s : optional species argument + # n : path to NetworKIN_201612_cutoffscore2.0.txt + # m : path to pSTY_Motifs.txt + # p : path to 2017-03_PSP_Kinase_Substrate_Dataset.txt + # r : path to 2017-03_PSP_Regulatory_sites.txt + ## options + # P : phospho_type + # F : function + # v : verbose output + ## output files + # o : path to output file + # O : path to "melted" output file + # D : path to output SQLite file + +sub usage() + { + print STDERR <<"EOH"; + This program given a list of phospho-peptides, finds protein information and upstream kinases. + usage: $0 [-hvd] -f FASTA_file + -h : this (help) message + -v : slightly verbose + -a : use SQLite less + ## input files + -i : path to input file, e.g., 'outputfile_STEP2.txt' + -f : path to UniProtDB/SwissProt FASTA + -s : optional species filter argument for PSP records; defaults to 'human' + -n : path to NetworKIN_201612_cutoffscore2.0.txt + -m : path to pSTY_Motifs.txt + -p : path to 2017-03_PSP_Kinase_Substrate_Dataset.txt + -r : path to 2017-03_PSP_Regulatory_sites.txt + ## options + -P : phospho_type + -F : function + ## output files + -o : path to output file + -O : path to "melted" output file + -D : path to output SQLite file + example: $0 +EOH + exit; + } + +sub format_localtime_iso8601 { + # ref: https://perldoc.perl.org/Time::HiRes + my ($seconds, $microseconds) = gettimeofday; + # ref: https://pubs.opengroup.org/onlinepubs/9699919799/functions/strftime.html + return strftime("%Y-%m-%dT%H:%M:%S",localtime(time)) . sprintf(".%03d", $microseconds/1000); +} + +sub replace_pSpTpY { + my ($formatted_sequence, $phospho_type) = @_; + if ($phospho_type eq 'y') { + $formatted_sequence =~ s/pS/S/g; + $formatted_sequence =~ s/pT/T/g; + $formatted_sequence =~ s/pY/y/g; + } + elsif ($phospho_type eq "sty") { + $formatted_sequence =~ s/pS/s/g; + $formatted_sequence =~ s/pT/t/g; + $formatted_sequence =~ s/pY/y/g; + } + $formatted_sequence; +} + +sub pseudo_sed +{ + # pseudo_sed produces "UniProt_ID\tDescription\tOS\tOX\tGN\tPE\tSV" + # Comments give the sed equivalent + my ($t) = @_; + my $s = $t; + # / GN=/!{ s:\(OX=[^ \t]*\):\1 GN=N/A:; }; + unless ($s =~ m / GN=/s) + { + $s =~ s :(OX=[^ \t]*):${1} GN=N/A:s; + } + # / PE=/!{ s:\(GN=[^ \t]*\):\1 PE=N/A:; }; + unless ($s =~ m / PE=/s) + { + $s =~ s :(GN=[^ \t]*):${1} PE=N/A:s; + } + # / SV=/!{ s:\(PE=[^ \t]*\):\1 SV=N/A:; }; + unless ($s =~ m / SV=/s) + { + $s =~ s :(PE=[^ \t]*):${1} SV=N/A:s; + } + # s/^sp.//; + $s =~ s :^...::s; + # s/[|]/\t/g; + $s =~ s :[|]:\t:sg; + if ( !($s =~ m/ OX=/s) + && !($s =~ m/ GN=/s) + && !($s =~ m/ PE=/s) + && !($s =~ m/ SV=/s) + ) { + # OS= is used elsewhere, but it's not helpful without OX and GN + $s =~ s/OS=/Species /g; + # supply sensible default values + $s .= "\tN/A\t-1\tN/A\tN/A\tN/A"; + } else { + # s/ OS=/\t/; + if ($s =~ m/ OS=/s) { $s =~ s: OS=:\t:s; } else { $s =~ s:(.*)\t:$1\tN/A\t:x; }; + # s/ OX=/\t/; + if ($s =~ m/ OX=/s) { $s =~ s: OX=:\t:s; } else { $s =~ s:(.*)\t:$1\t-1\t:x; }; + # s/ GN=/\t/; + if ($s =~ m/ GN=/s) { $s =~ s: GN=:\t:s; } else { $s =~ s:(.*)\t:$1\tN/A\t:x; }; + # s/ PE=/\t/; + if ($s =~ m/ PE=/s) { $s =~ s: PE=:\t:s; } else { $s =~ s:(.*)\t:$1\tN/A\t:x; }; + # s/ SV=/\t/; + if ($s =~ m/ SV=/s) { $s =~ s: SV=:\t:s; } else { $s =~ s:(.*)\t:$1\tN/A\t:x; }; + } + return $s; +} # sub pseudo_sed + +getopts('i:f:s:n:m:p:r:P:F:o:O:D:hva', \%opts) ; + + +if (exists($opts{'h'})) { + usage(); +} +if (exists($opts{'a'})) { + $USE_SEARCH_PPEP_PY = 0; +} +if (exists($opts{'v'})) { + $verbose = 1; +} else { + $verbose = 0; +} +if (!exists($opts{'i'}) || !-e $opts{'i'}) { + die('Input File not found'); +} else { + $file_in = $opts{'i'}; +} +if (!exists($opts{'f'}) || !-e $opts{'f'}) { + die('FASTA not found'); +} else { + $fasta_in = $opts{'f'}; + $use_sqlite = 0; +} +my $species; +if ((!exists($opts{'s'})) || ($opts{'s'} eq '')) { + $species = 'human'; +} else { + $species = $opts{'s'}; + print "'-s' option is '$species'\n"; +} +print "species filter is '$species'\n"; + +if (!exists($opts{'n'}) || !-e $opts{'n'}) { + die('Input NetworKIN File not found'); +} else { + $networkin_in = $opts{'n'}; +} +if (!exists($opts{'m'}) || !-e $opts{'m'}) { + die('Input pSTY_Motifs File not found'); +} else { + $motifs_in = $opts{'m'}; +} +if (!exists($opts{'p'}) || !-e $opts{'p'}) { + die('Input PSP_Kinase_Substrate_Dataset File not found'); +} else { + $PSP_Kinase_Substrate_in = $opts{'p'}; +} +if (!exists($opts{'r'}) || !-e $opts{'r'}) { + die('Input PSP_Regulatory_sites File not found'); +} else { + $PSP_Regulatory_Sites_in = $opts{'r'}; +} +if (exists($opts{'P'})) { + $phospho_type = $opts{'P'}; +} +else { + $phospho_type = "sty"; +} +if (exists($opts{'F'})) { + $average_or_sum = $opts{'F'}; +} +else { + $average_or_sum = "sum"; +} +if (exists($opts{'D'})) { + $db_out = $opts{'D'}; +} +else { + $db_out = "db_out.sqlite"; +} +if (exists($opts{'O'})) { + $file_melt = $opts{'O'}; +} +else { + $file_melt = "output_melt.tsv"; +} +if (exists($opts{'o'})) { + $file_out = $opts{'o'}; +} +else { + $file_out = "output.tsv"; +} + + +############################################################################################################################### +# Print the relevant file names to the screen +############################################################################################################################### +# print "\nData file: $data_in\nFASTA file: $fasta_in\nSpecies: $species\nOutput file: $motifs_out\n\n"; +print "\n--- parameters:\n"; +print "Data file: $file_in\nAverage or sum identical p-sites? $average_or_sum\nOutput file: $file_out\nMelted map: $file_melt\n"; +if ($use_sqlite == 0) { + print "Motifs file: $motifs_in\nNetworKIN file: networkin_in\nPhosphosite kinase substrate data: $PSP_Kinase_Substrate_in\nPhosphosite regulatory site data: $PSP_Regulatory_Sites_in\nUniProtKB/SwissProt FASTA file: $fasta_in\nOutput SQLite file: $db_out\n"; +} else { + print "Motifs file: $motifs_in\nNetworKIN file: networkin_in\nPhosphosite kinase substrate data: $PSP_Kinase_Substrate_in\nPhosphosite regulatory site data: $PSP_Regulatory_Sites_in\nUniProtKB/SwissProt SQLIte file: $dbfile\nOutput SQLite file: $db_out\n"; +} +print "...\n\n"; + +print "Phospho-residues(s) = $phospho_type\n\n"; +if ($phospho_type ne 'y') { + if ($phospho_type ne 'sty') { + die "\nUsage error:\nYou must choose a phospho-type, either y or sty\n\n"; + } +} + +############################################################################################################################### +# read the input data file +# average or sum identical phospho-sites, depending on the value of $average_or_sum +############################################################################################################################### + +open (IN, "$file_in") or die "I couldn't find the input file: $file_in\n"; + +die "\n\nScript died: You must choose either average or sum for \$average_or_sum\n\n" if (($average_or_sum ne "sum") && ($average_or_sum ne "average")) ; + + +$line = 0; + +while (<IN>) { + chomp; + my @x = split(/\t/); + for my $n (0 .. $#x) {$x[$n] =~ s/\r//g; $x[$n] =~ s/\n//g; $x[$n] =~ s/\"//g;} + + # Read in the samples + if ($line == 0) { + for my $n (1 .. $#x) { + push (@samples, $x[$n]); + $sample_id_lut{$x[$n]} = $n; + } + $line++; + } else { + # check whether we have already seen a phospho-peptide + if (exists($data{$x[0]})) { + if ($average_or_sum eq "sum") { # add the data + # unload the data + @tmp_data = (); foreach (@{$data{$x[0]}}) { push(@tmp_data, $_); } + # add the new data and repack + for my $k (0 .. $#tmp_data) { $tmp_data[$k] = $tmp_data[$k] + $x[$k+1]; } + $all_data{$x[0]} = (); for my $k (0 .. $#tmp_data) { push(@{$all_data{$x[0]}}, $tmp_data[$k]); } + + } elsif ($average_or_sum eq "average") { # average the data + # unload the data + @tmp_data = (); foreach (@{$all_data{$x[0]}}) { push(@tmp_data, $_); } + # average with the new data and repack + for my $k (0 .. $#tmp_data) { $tmp_data[$k] = ( $tmp_data[$k]*$n{$x[0]} + $x[0] ) / ($n{$x[0]} + 1); } + $n{$x[0]}++; + $data{$x[0]} = (); for my $k (0 .. $#tmp_data) { push(@{$data{$x[0]}}, $tmp_data[$k]); } + } + } + # if the phospho-sequence has not been seen, save the data + else { + for my $k (1 .. $#x) { push(@{$data{$x[0]}}, $x[$k]); } + $n{$x[0]} = 1; + } + } +} +close(IN); + + +############################################################################################################################### +# Search the FASTA database for phospho-sites and motifs +# +# based on Retrieve_p_peptide_motifs_v2.pl +############################################################################################################################### + + +############################################################################################################################### +# +# Read in the Data file: +# 1) make @p_peptides array as in the original file +# 2) make @non_p_peptides array w/o residue modifications (p, #, other) +# +############################################################################################################################### + +foreach my $peptide (keys %data) { + $peptide =~ s/s/pS/g; $peptide =~ s/t/pT/g; $peptide =~ s/y/pY/g; + push (@p_peptides, $peptide); + $peptide =~ s/p//g; + push(@non_p_peptides, $peptide); +} + +if ($use_sqlite == 0) { + ############################################################################################################################### + # + # Read in the UniProtKB/Swiss-Prot data from FASTA; save to @sequences array and SQLite output database + # + ############################################################################################################################### + + # e.g. + # >sp|Q9Y3B9|RRP15_HUMAN RRP15-like protein OS=Homo sapiens OX=9606 GN=RRP15 PE=1 SV=2 + # MAAAAPDSRVSEEENLKKTPKKKMKMVTGAVASVLEDEATDTSDSEGSCGSEKDHFYSDD + # DAIEADSEGDAEPCDKENENDGESSVGTNMGWADAMAKVLNKKTPESKPTILVKNKKLEK + # EKEKLKQERLEKIKQRDKRLEWEMMCRVKPDVVQDKETERNLQRIATRGVVQLFNAVQKH + # QKNVDEKVKEAGSSMRKRAKLISTVSKKDFISVLRGMDGSTNETASSRKKPKAKQTEVKS + # EEGPGWTILRDDFMMGASMKDWDKESDGPDDSRPESASDSDT + # accession: Q9Y3B9 + # name: RRP15_HUMAN RRP15-like protein OS=Homo sapiens OX=9606 GN=RRP15 PE=1 SV=2 + # sequence: MAAAAPDSRVSEEENLKKTPKKKMKMVTGAVASVLEDEATDTSDSEGSCGSEKDHFYSDD DAIEADSEGDAEPCDKENENDGESSVGTNMGWADAMAKVLNKKTPESKPTILVKNKKLEK EKEKLKQERLEKIKQRDKRLEWEMMCRVKPDVVQDKETERNLQRIATRGVVQLFNAVQKH QKNVDEKVKEAGSSMRKRAKLISTVSKKDFISVLRGMDGSTNETASSRKKPKAKQTEVKS EEGPGWTILRDDFMMGASMKDWDKESDGPDDSRPESASDSDT + # + # e.g. + # >gi|114939|sp|P00722.2|BGAL_ECOLI Beta-galactosidase (Lactase) cRAP + # >gi|52001466|sp|P00366.2|DHE3_BOVIN Glutamate dehydrogenase 1, mitochondrial precursor (GDH) cRAP + # + # e.g. + # >zs|P00009.24.AR-V2_1.zs|zs_peptide_0024_AR-V2_1 + + + open (IN1, "$fasta_in") or die "I couldn't find $fasta_in\n"; + print "Reading FASTA file $fasta_in\n"; + # ref: https://perldoc.perl.org/perlsyn#Compound-Statements + # "If the condition expression of a while statement is based on any of + # a group of iterative expression types then it gets some magic treatment. + # The affected iterative expression types are readline, the <FILEHANDLE> + # input operator, readdir, glob, the <PATTERN> globbing operator, and + # `each`. If the condition expression is one of these expression types, + # then the value yielded by the iterative operator will be implicitly + # assigned to `$_`." + while (<IN1>) { + chomp; + # ref: https://perldoc.perl.org/functions/split#split-/PATTERN/,EXPR + # "If only PATTERN is given, EXPR defaults to $_." + my (@x) = split(/\|/); + # begin FIX >gi|114939|sp|P00722.2|BGAL_ECOLI Beta-galactosidase (Lactase) cRAP + if (@x > 3) { + @x = (">".$x[$#x - 2], $x[$#x - 1], $x[$#x]); + } + # end FIX >gi|114939|sp|P00722.2|BGAL_ECOLI Beta-galactosidase (Lactase) cRAP + for my $i (0 .. $#x) { + $x[$i] =~ s/\r//g; $x[$i] =~ s/\n//g; $x[$i] =~ s/\"//g; } + # Use of uninitialized value $x[0] in pattern match (m//) at /home/rstudio/src/mqppep/tools/mqppep/PhosphoPeptide_Upstream_Kinase_Mapping.pl line 411, <IN1> line 3. + if (exists($x[0])) { + if ($x[0] =~ /^>/) { + # parsing header line + $x[0] =~ s/\>//g; + push (@databases, $x[0]); + push (@accessions, $x[1]); + push (@names, $x[2]); + # format tags of standard UniProtKB headers as tab-separated values + # pseudo_sed produces "UniProt_ID\tDescription\tOS\tOX\tGN\tPE\tSV" + $_ = pseudo_sed(join "\t", (">".$x[0], $x[1], $x[2])); + # append tab as separator between header and sequence + s/$/\t/; + # parsed_fasta gets "UniProt_ID\tDescription\tOS\tOX\tGN\tPE\tSV\t" + print "push (\@parsed_fasta, $_)\n" if (0 && $x[0] ne "zs"); + push (@parsed_fasta, $_); + } elsif ($x[0] =~ /^\w/) { + # line is a portion of the sequence + if (defined $sequences[$#accessions]) { + $sequences[$#accessions] = $sequences[$#accessions].$x[0]; + } else { + $sequences[$#accessions] = $x[0]; + } + $parsed_fasta[$#accessions] = $parsed_fasta[$#accessions].$x[0]; + } + } + } + close IN1; + print "Done Reading FASTA file $fasta_in\n"; + $dbfile = $db_out; + print "Begin writing $dbfile at " . format_localtime_iso8601() . "\n"; + $dbh = DBI->connect("dbi:SQLite:$dbfile", undef, undef); + my $auto_commit = $dbh->{AutoCommit}; + print "auto_commit was $auto_commit and is now 0\n" if ($verbose); + $dbh->{AutoCommit} = 0; + + # begin DDL-to-SQLite + # --- + $stmth = $dbh->prepare(" + DROP TABLE IF EXISTS UniProtKB; + "); + $stmth->execute(); + + $stmth = $dbh->prepare(" + CREATE TABLE UniProtKB ( + Uniprot_ID TEXT PRIMARY KEY ON CONFLICT IGNORE, + Description TEXT, + Organism_Name TEXT, + Organism_ID INTEGER, + Gene_Name TEXT, + PE TEXT, + SV TEXT, + Sequence TEXT, + Database TEXT + ) + "); + $stmth->execute(); + $stmth = $dbh->prepare(" + CREATE UNIQUE INDEX idx_uniq_UniProtKB_0 on UniProtKB(Uniprot_ID); + "); + $stmth->execute(); + $stmth = $dbh->prepare(" + CREATE INDEX idx_UniProtKB_0 on UniProtKB(Gene_Name); + "); + $stmth->execute(); + # ... + # end DDL-to-SQLite + + # insert all rows + # begin store-to-SQLite "UniProtKB" table + # --- + $stmth = $dbh->prepare(" + INSERT INTO UniProtKB ( + Uniprot_ID, + Description, + Organism_Name, + Organism_ID, + Gene_Name, + PE, + SV, + Sequence, + Database + ) VALUES (?,?,?,?,?,?,?,?,?) + "); + my $row_count = 1; + my $row_string; + my (@row, @rows); + my $wrd; + while ( scalar @parsed_fasta > 0 ) { + $database = $databases[$#parsed_fasta]; + # row_string gets "UniProt_ID\tDescription\tOS\tOX\tGN\tPE\tSV\t" + # 1 2 3 4 5 6 7 sequence database + $row_string = pop(@parsed_fasta); + @row = (split /\t/, $row_string); + if ((not exists($row[4])) || ($row[4] eq "")) { + die("invalid fasta line\n$row_string\n"); + }; + if ($row[4] eq "N/A") { + print "Organism_ID is 'N/A' for row $row_count:\n'$row_string'\n"; + $row[4] = -1; + }; + for $i (1..3,5..8) { + #BIND print "bind_param $i, $row[$i]\n"; + $stmth->bind_param($i, $row[$i]); + } + #BIND print "bind_param 4, $row[4]\n"; + $stmth->bind_param(9, $database); + #BIND print "bind_param 4, $row[4]\n"; + $stmth->bind_param(4, $row[4], { TYPE => SQL_INTEGER }); + if (not $stmth->execute()) { + print "Error in row $row_count: " . $dbh->errstr . "\n"; + print "Row $row_count: $row_string\n"; + print "Row $row_count: " . ($row_string =~ s/\t/@/g) . "\n"; + } + if (0 && $database ne "zs") { + print "row_count: $row_count\n"; + #### print "row_string: $row_string\n"; + print "Row $row_count: $row_string\n"; + for $i (1..3,5..8) { + print "bind_param $i, $row[$i]\n" if (exists($row[$i])); + } + print "bind_param 4, $row[4]\n" if (exists($row[4])); + print "bind_param 9, $database\n"; + }; + $row_count += 1; + } + # ... + # end store-to-SQLite "UniProtKB" table + + print "begin commit at " . format_localtime_iso8601() . "\n"; + $dbh->{AutoCommit} = $auto_commit; + print "auto_commit is now $auto_commit\n" if ($verbose); + $dbh->disconnect if ( defined $dbh ); + print "Finished writing $dbfile at " . format_localtime_iso8601() . "\n\n"; + $dbtype = "FASTA"; +} + +if ($use_sqlite == 1) { + ############################################################################################################################### + # + # Read in the UniProtKB/Swiss-Prot data from SQLite; save to @sequences array + # + ############################################################################################################################### + + copy($dbfile, $db_out) or die "Copy $dbfile to $db_out failed: $!"; + + # https://metacpan.org/pod/DBD::SQLite#Read-Only-Database + $dbh = DBI->connect("dbi:SQLite:$dbfile", undef, undef, { + sqlite_open_flags => SQLITE_OPEN_READONLY, + }); + print "DB connection $dbh is to $dbfile\n"; + + # Uniprot_ID, Description, Organism_Name, Organism_ID, Gene_Name, PE, SV, Sequence + $stmth = $dbh->prepare(" + SELECT Uniprot_ID + , Description + || CASE WHEN Organism_Name = 'N/A' THEN '' ELSE ' OS=' || Organism_Name END + || CASE WHEN Organism_ID = -1 THEN '' ELSE ' OX=' || Organism_ID END + || CASE WHEN Gene_Name = 'N/A' THEN '' ELSE ' GN=' || Gene_Name END + || CASE WHEN PE = 'N/A' THEN '' ELSE ' PE=' || PE END + || CASE WHEN SV = 'N/A' THEN '' ELSE ' SV=' || SV END + AS Description + , Sequence + , Database + FROM + UniProtKB + "); + $stmth->execute(); + @col_names = @{$stmth->{NAME}}; + print "\nColumn names selected from UniProtKB SQLite table: " . join(", ", @col_names) . "\n\n" if ($verbose); + while (my @row = $stmth->fetchrow_array) { + push (@names, $row[1]); # redacted Description + push (@accessions, $row[0]); # Uniprot_ID + $sequences[$#accessions] = $row[2]; # Sequence + push (@databases, $row[3]); # Database (should be 'sp') + } + + $dbh->disconnect if ( defined $dbh ); + + print "Done Reading UniProtKB/Swiss-Prot file $dbfile\n\n"; + $dbtype = "SQLite"; +} + +print "$#accessions accessions were read from the UniProtKB/Swiss-Prot $dbtype file\n"; + +###################### + $dbh = DBI->connect("dbi:SQLite:$dbfile", undef, undef); + $stmth = $dbh->prepare(" + INSERT INTO UniProtKB ( + Uniprot_ID, + Description, + Organism_Name, + Organism_ID, + Gene_Name, + PE, + SV, + Sequence, + Database + ) VALUES ( + 'No Uniprot_ID', + 'NO_GENE_SYMBOL No Description', + 'No Organism_Name', + 0, + '$FAILED_MATCH_GENE_NAME', + '0', + '0', + '$FAILED_MATCH_SEQ', + 'No Database' + ) + "); + if (not $stmth->execute()) { + print "Error inserting dummy row into UniProtKB: $stmth->errstr\n"; + } + $dbh->disconnect if ( defined $dbh ); +###################### + +@timeData = localtime(time); +print "\n--- Start search at " . format_localtime_iso8601() ."\n"; + +print " --> Calling 'search_ppep' script\n\n"; +if ($verbose) { + $i = system("python $dirname/search_ppep.py -u $db_out -p $file_in --verbose"); +} else { + $i = system("python $dirname/search_ppep.py -u $db_out -p $file_in"); +} +if ($i) { + print "python $dirname/search_ppep.py -u $db_out -p $file_in\n exited with exit code $i\n"; + die "Search failed for phosphopeptides in SwissProt/SQLite file."; +} +print " <-- Returned from 'search_ppep' script\n"; + +@timeData = localtime(time); +print "... Finished search at " . format_localtime_iso8601() ."\n\n"; + + +############################################################################################################################### +# +# Match the non_p_peptides to the @sequences array: +# 1) Format the motifs +/- 10 residues around the phospho-site +# 2) Print the original data plus the phospho-motif to the output file +# +############################################################################################################################### + + +print "--- Match the non_p_peptides to the \@sequences array:\n"; + +if ($USE_SEARCH_PPEP_PY) { + print "Find the matching protein sequence(s) for the peptide using SQLite\n"; +} else { + print "Find the matching protein sequence(s) for the peptide using slow search\n"; +} + +# https://metacpan.org/pod/DBD::SQLite#Read-Only-Database +$dbh = DBI->connect("dbi:SQLite:$db_out", undef, undef, { + sqlite_open_flags => SQLITE_OPEN_READONLY, +}); +print "DB connection $dbh is to $db_out\n"; + +# CREATE VIEW uniprotid_pep_ppep AS +# SELECT deppep_UniProtKB.UniprotKB_ID AS accession +# , deppep.seq AS peptide +# , ppep.seq AS phosphopeptide +# , UniProtKB.Sequence AS sequence +# , UniProtKB.Description AS description +# FROM ppep, deppep, deppep_UniProtKB, UniProtKB +# WHERE deppep.id = ppep.deppep_id +# AND deppep.id = deppep_UniProtKB.deppep_id +# AND deppep_UniProtKB.UniprotKB_ID = UniProtKB.Uniprot_ID +# ORDER BY UniprotKB_ID, deppep.seq, ppep.seq; + +my %ppep_to_count_lut; +print "start select peptide counts " . format_localtime_iso8601() . "\n"; +my $uniprotkb_pep_ppep_view_stmth = $dbh->prepare(" + SELECT DISTINCT + phosphopeptide + , count(*) as i + FROM + uniprotkb_pep_ppep_view + GROUP BY + phosphopeptide + ORDER BY + phosphopeptide +"); +if (not $uniprotkb_pep_ppep_view_stmth->execute()) { + die "Error fetching peptide counts: $uniprotkb_pep_ppep_view_stmth->errstr\n"; +} +while (my @row = $uniprotkb_pep_ppep_view_stmth->fetchrow_array) { + $ppep_to_count_lut{$row[0]} = $row[1]; + #print "\$ppep_to_count_lut{$row[0]} = $ppep_to_count_lut{$row[0]}\n"; +} + +# accession, peptide, sequence, description, phosphopeptide, long_description, pos_start, pos_end, scrubbed, ppep_id +# 0 1 2 3 4 5 6 7 8 9 +my $COL_ACCESSION = 0; +my $COL_PEPTIDE = 1; +my $COL_SEQUENCE = 2; +my $COL_DESCRIPTION = 3; +my $COL_PHOSPHOPEPTIDE = 4; +my $COL_LONG_DESCRIPTION = 5; +my $COL_POS_START = 6; +my $COL_POS_END = 7; +my $COL_SCRUBBED = 8; +my $COL_PPEP_ID = 9; + +my %ppep_to_row_lut; +print "start select all records without qualification " . format_localtime_iso8601() . "\n"; +$uniprotkb_pep_ppep_view_stmth = $dbh->prepare(" + SELECT DISTINCT + accession + , peptide + , sequence + , description + , phosphopeptide + , long_description + , pos_start + , pos_end + , scrubbed + , ppep_id + FROM + uniprotkb_pep_ppep_view + ORDER BY + phosphopeptide +"); +if (not $uniprotkb_pep_ppep_view_stmth->execute()) { + die "Error fetching all records without qualification: $uniprotkb_pep_ppep_view_stmth->errstr\n"; +} +my $current_ppep; +my $counter = 0; +my $former_ppep = ""; +@tmp_matches = (); +@tmp_accessions = (); +@tmp_names = (); +@tmp_sites = (); +while (my @row = $uniprotkb_pep_ppep_view_stmth->fetchrow_array) { + # Identify phosphopeptide for current row; + # it is an error for it to change when the counter is not zero. + $current_ppep = $row[$COL_PHOSPHOPEPTIDE]; + + # when counter is zero, prepare for a new phosphopeptide + if (not $current_ppep eq $former_ppep) { + die "counter is $counter instead of zero" if ($counter != 0); + $ppep_id_lut{$current_ppep} = $row[$COL_PPEP_ID]; + print "next phosphpepetide: $current_ppep; id: $ppep_id_lut{$current_ppep}\n" if ($verbose); + $counter = $ppep_to_count_lut{$current_ppep}; + @tmp_matches = (); + @tmp_accessions = (); + @tmp_names = (); + @tmp_sites = (); + } + + if ($USE_SEARCH_PPEP_PY) { + push(@tmp_matches, $row[ $COL_SEQUENCE ]); + push(@tmp_accessions, $row[ $COL_ACCESSION ]); + push(@tmp_names, $row[ $COL_LONG_DESCRIPTION ]); + push(@tmp_sites, $row[ $COL_POS_START ]); + } + + # Prepare counter and phosphopeptide tracker for next row + $former_ppep = $current_ppep; + $counter -= 1; + + # Set trackers for later use after last instance of current phosphopeptide + if ($counter == 0) { + if ($USE_SEARCH_PPEP_PY) { + $matched_sequences{$current_ppep} = [ @tmp_matches ]; + $accessions{ $current_ppep} = [ @tmp_accessions ]; + $names{ $current_ppep} = [ @tmp_names ]; + $sites{ $current_ppep} = [ @tmp_sites ]; + } + } +} + + +print "end select all records without qualification " . format_localtime_iso8601() . "\n"; + +for my $j (0 .. $#p_peptides) { + + #Find the matching protein sequence(s) for the peptide using SQLite + my ($site, $sequence); + my (@row, @rows); + my $match = 0; + my $p_peptide = $p_peptides[$j]; + @tmp_matches = (); + @tmp_accessions = (); + @tmp_names = (); + @tmp_sites = (); + + #Find the matching protein sequence(s) for the peptide using slow search + $site = -1; + unless ($USE_SEARCH_PPEP_PY) { + for my $k (0 .. $#sequences) { + $site = index($sequences[$k], $non_p_peptides[$j]); + if ($site != -1) { + push(@tmp_matches, $sequences[$k]); + push(@tmp_accessions, $accessions[$k]); + push(@tmp_names, $names[$k]); + push(@tmp_sites, $site); + } + # print "Non-phosphpeptide $non_p_peptides[$j] matched accession $accessions[$k] ($names[$k]) at site $site\n"; + $site = -1; $match++; + # print "tmp_accessions @tmp_accessions \n"; + } + if ($match == 0) { # Check to see if no match was found. Skip to next if no match found. + print "Warning: Failed match for $p_peptides[$j]\n"; + $matched_sequences{$p_peptides[$j]} = \@failed_match; + push(@failed_matches,$p_peptides[$j]); + next; + } else { + $matched_sequences{$p_peptides[$j]} = [ @tmp_matches ]; + $accessions{$p_peptides[$j]} = [ @tmp_accessions ]; + $names{$p_peptides[$j]} = [ @tmp_names ]; + $sites{$p_peptides[$j]} = [ @tmp_sites ]; + } + } + +} # end for my $j (0 .. $#p_peptides) + +print "... Finished match the non_p_peptides at " . format_localtime_iso8601() ."\n\n"; + +print "--- Match the p_peptides to the \@sequences array:\n"; + +for my $peptide_to_match ( keys %matched_sequences ) { + if (grep($peptide_to_match, @failed_matches)) { + print "Failed to match peptide $peptide_to_match\n"; + } + next if (grep($peptide_to_match, @failed_matches)); + my @matches = @{$matched_sequences{$peptide_to_match}}; + @tmp_motifs_array = (); + for my $i (0 .. $#matches) { + + # Find the location of the phospo-site in the sequence(s) + $tmp_site = 0; my $offset = 0; + my $tmp_p_peptide = $peptide_to_match; + $tmp_p_peptide =~ s/#//g; $tmp_p_peptide =~ s/\d//g; $tmp_p_peptide =~ s/\_//g; $tmp_p_peptide =~ s/\.//g; + + # Find all phosphorylated residues in the p_peptide + @p_sites = (); + while ($tmp_site != -1) { + $tmp_site = index($tmp_p_peptide, 'p', $offset); + if ($tmp_site != -1) {push (@p_sites, $tmp_site);} + $offset = $tmp_site + 1; + $tmp_p_peptide =~ s/p//; + } + @tmp_p_residues = (); + for my $l (0 .. $#p_sites) { + next if not defined $sites{$peptide_to_match}[$i]; + + push (@tmp_p_residues, $p_sites[$l] + $sites{$peptide_to_match}[$i]); + + # Match the sequences around the phospho residues to find the motifs + my ($desired_residues_L, $desired_residues_R); + if ($tmp_p_residues[0] - 10 < 0) { #check to see if there are fewer than 10 residues left of the first p-site + # eg, XXXpYXX want $desired_residues_L = 3, $p_residues[0] = 3 + $desired_residues_L = $tmp_p_residues[0]; + } + else { + $desired_residues_L = 10; + } + my $seq_length = length($matched_sequences{$peptide_to_match}[$i]); + if ($tmp_p_residues[$#tmp_p_residues] + 10 > $seq_length) { #check to see if there are fewer than 10 residues right of the last p-site + $desired_residues_R = $seq_length - ($tmp_p_residues[$#tmp_p_residues] + 1); + # eg, XXXpYXX want $desired_residues_R = 2, $seq_length = 6, $p_residues[$#p_residues] = 3 + # print "Line 170: seq_length = $seq_length\tp_residue = $p_residues[$#p_residues]\n"; + } + else { + $desired_residues_R = 10; + } + + my $total_length = $desired_residues_L + $tmp_p_residues[$#tmp_p_residues] - $tmp_p_residues[0] + $desired_residues_R + 1; + my $arg2 = $tmp_p_residues[0] - $desired_residues_L; + my $arg1 = $matched_sequences{$peptide_to_match}[$i]; + + if (($total_length > 0) && (length($arg1) > $arg2 + $total_length - 1)) { + $tmp_motif = substr($arg1, $arg2, $total_length); + + # Put the "p" back in front of the appropriate phospho-residue(s). + my (@tmp_residues, $tmp_position); + for my $m (0 .. $#p_sites) { + # print "Line 183: $p_sites[$m]\n"; + if ($m == 0) { + $tmp_position = $desired_residues_L; + } else { + $tmp_position = $desired_residues_L + $p_sites[$m] - $p_sites[0]; + } + if ($tmp_position < length($tmp_motif) + 1) { + push (@tmp_residues, substr($tmp_motif, $tmp_position, 1)); + if ($tmp_residues[$m] eq "S") {substr($tmp_motif, $tmp_position, 1, "s");} + if ($tmp_residues[$m] eq "T") {substr($tmp_motif, $tmp_position, 1, "t");} + if ($tmp_residues[$m] eq "Y") {substr($tmp_motif, $tmp_position, 1, "y");} + } + } + + $tmp_motif =~ s/s/pS/g; $tmp_motif =~ s/t/pT/g; $tmp_motif =~ s/y/pY/g; + + # Comment out on 8.10.13 to remove the numbers from motifs + my $left_residue = $tmp_p_residues[0] - $desired_residues_L+1; + my $right_residue = $tmp_p_residues[$#tmp_p_residues] + $desired_residues_R+1; + $tmp_motif = $left_residue."-[ ".$tmp_motif." ]-".$right_residue; + push(@tmp_motifs_array, $tmp_motif); + $residues{$peptide_to_match}{$i} = [ @tmp_residues ]; + $p_residues{$peptide_to_match}{$i} = [ @tmp_p_residues ]; + } + } + $p_motifs{$peptide_to_match} = [ @tmp_motifs_array ]; + } # end for my $i (0 .. $#matches) ### this bracket could be in the wrong place +} + +print "... Finished match the p_peptides to the \@sequences array at " . format_localtime_iso8601() ."\n\n"; + +############################################################################################################################### +# +# Annotate the peptides with the NetworKIN predictions and HPRD / Phosida kinase motifs +# +############################################################################################################################### + + +print "--- Reading various site data:\n"; + +############################################################################################################################### +# +# Read the NetworKIN_predictions file: +# 1) make a "kinases_observed" array +# 2) annotate the phospho-substrates with the appropriate kinase +# +############################################################################################################################### +my $SITE_KINASE_SUBSTRATE = 1; +$site_description{$SITE_KINASE_SUBSTRATE} = "NetworKIN"; + +open (IN1, "$networkin_in") or die "I couldn't find $networkin_in\n"; +print "Reading the NetworKIN data: $networkin_in\n"; +while (<IN1>) { + chomp; + my (@x) = split(/\t/); + for my $i (0 .. $#x) { + $x[$i] =~ s/\r//g; $x[$i] =~ s/\n//g; $x[$i] =~ s/\"//g; + } + next if ($x[0] eq "#substrate"); + if (exists ($kinases -> {$x[2]})) { + #do nothing + } + else { + $kinases -> {$x[2]} = $x[2]; + push (@kinases_observed, $x[2]); + } + my $tmp = $x[10]."_".$x[2]; #eg, REEILsEMKKV_PKCalpha + if (exists($p_sequence_kinase -> {$tmp})) { + #do nothing + } + else { + $p_sequence_kinase -> {$tmp} = $tmp; + } +} +close IN1; + +############################################################################################################################### +# +# Read the Kinase motifs file: +# 1) make a "motif_sequence" array +# +############################################################################################################################### + +# file format (tab separated): +# x[0] = quasi-primary key (character), e.g., '17' or '23a' +# x[1] = pattern (egrep pattern), e.g., '(M|I|L|V|F|Y).R..(pS|pT)' +# x[2] = description, e.g., 'PKA_Phosida' or '14-3-3 domain binding motif (HPRD)' or 'Akt kinase substrate motif (HPRD & Phosida)' +# "counter" "pcre" "symbol" "description" "pubmed_id" "classification" "source" +# "1" "R.R..(pS|pT)(F|L)" "PKB_group" "Akt kinase" "https://pubmed.ncbi.nlm.nih.gov/?term=8985174" "kinase substrate" "HPRD" +# x[3] = old description, i.e., description in Amanchy (HPRD) and Phosida tables +# x[4] = pubmed id +# x[5] = classification +# x[6] = source (Phosida or HPRD) +my $SITE_HPRD = 2; +$site_description{$SITE_HPRD} = "HPRD"; +$site_id{$site_description{$SITE_HPRD}} = $SITE_HPRD; +my $SITE_PHOSIDA = 4; +$site_description{$SITE_PHOSIDA} = "Phosida"; +$site_id{$site_description{$SITE_PHOSIDA}} = $SITE_PHOSIDA; + +open (IN2, "$motifs_in") or die "I couldn't find $motifs_in\n"; +print "Reading the Motifs file: $motifs_in\n"; + +while (<IN2>) { + chomp; + my (@x) = split(/\t/); + my $tmp_motif_description; + if ($#x == 6) { # weirdly, a @list of length seven has $#list == 6 + # remove double-quotes which are helpful or necessary for Excel + $x[6] =~ s/\"//g; + $tmp_motif_description = $x[6]; + } else { + $tmp_motif_description = "motif"; + } + for my $i (0 .. 2) { + # remove any embedded CR or LF (none should exist) + $x[$i] =~ s/\r//g; + $x[$i] =~ s/\n//g; + # remove double-quotes which are helpful or necessary for Excel + $x[$i] =~ s/\"//g; + } + if (exists ($motif_type{$x[2]})) { + #ACE-2022.06.20 $motif_type{$x[1]} = $motif_type{$x[1]}." & ".$x[2]; + $motif_type{$x[2]} = $motif_type{$x[2]}."|".$x[2]; + } else { + $motif_type{$x[2]} = $x[2]; + $motif_count{$x[1]} = 0; + push (@motif_sequence, $x[1]); + push (@motif_description, $tmp_motif_description); + push (@motif_type_key_ary, $x[2]) + } +} +close (IN2); + + +############################################################################################################################### +# 6.28.2011 +# Read PSP_Kinase_Substrate data: +# 1) make a "kinases_PhosphoSite" array +# 2) annotate the phospho-substrates with the appropriate kinase +# +# Columns: +# (0) GENE +# (1) KINASE +# (2) KIN_ACC_ID +# (3) KIN_ORGANISM +# (4) SUBSTRATE +# (5) SUB_GENE_ID +# (6) SUB_ACC_ID +# (7) SUB_GENE +# (8) SUB_ORGANISM +# (9) SUB_MOD_RSD +# (10) SITE_GRP_ID +# (11) SITE_+/-7_AA +# (12) DOMAIN +# (13) IN_VIVO_RXN +# (14) IN_VITRO_RXN +# (15) CST_CAT# +############################################################################################################################### + +my $SITE_PHOSPHOSITE = 3; +$site_description{$SITE_PHOSPHOSITE} = "PhosphoSite"; + + +$line = 0; + +open (IN3, "$PSP_Kinase_Substrate_in") or die "I couldn't find $PSP_Kinase_Substrate_in\n"; +print "Reading the PhosphoSite Kinase-Substrate data: $PSP_Kinase_Substrate_in\n"; + +while (<IN3>) { + chomp; + my (@x) = split(/\t/); + for my $i (0 .. $#x) { + $x[$i] =~ s/\r//g; $x[$i] =~ s/\n//g; $x[$i] =~ s/\"//g; + } + if ($line != 0) { + if (($species eq $x[3]) && ($species eq $x[8])) { + if (exists ($kinases_PhosphoSite -> {$x[0]})) { + #do nothing + } + else { + $kinases_PhosphoSite -> {$x[0]} = $x[0]; + push (@kinases_PhosphoSite, $x[0]); + } + my $offset = 0; + # Replace the superfluous lower case s, t and y + my @lowercase = ('s','t','y'); + my @uppercase = ('S','T','Y'); + for my $k (0 .. 2) { + my $site = 0; + while ($site != -1) { + $site = index($x[11],$lowercase[$k], $offset); + if (($site != 7) && ($site != -1)) {substr($x[11], $site, 1, $uppercase[$k]);} + $offset = $site + 1; + } + } + my $tmp = $x[11]."_".$x[0]; #eg, RTPGRPLsSYGMDSR_PAK2 + if (exists($p_sequence_kinase_PhosphoSite -> {$tmp})) { + #do nothing + } + else { + $p_sequence_kinase_PhosphoSite -> {$tmp} = $tmp; + } + } + else { + # do nothing + #print "PSP_kinase_substrate line rejected because KIN_ORGANISM is '$x[3]' and SUB_ORGANISM is '$x[8]': $line\n"; + } + } + $line++; +} +close IN3; + + +############################################################################################################################### +# Read PhosphoSite regulatory site data: +# 1) make a "regulatory_sites_PhosphoSite" hash +# +# Columns: +# (0) GENE +# (2) PROT_TYPE +# (3) ACC_ID +# (4) GENE_ID +# (5) HU_CHR_LOC +# (6) ORGANISM --> %organism +# (7) MOD_RSD +# (8) SITE_GRP_ID +# (9) SITE_+/-7_AA --> %regulatory_sites_PhosphoSite_hash +# (10) DOMAIN --> %domain +# (11) ON_FUNCTION --> %ON_FUNCTION +# (12) ON_PROCESS --> %ON_PROCESS +# (13) ON_PROT_INTERACT --> %ON_PROT_INTERACT +# (14) ON_OTHER_INTERACT --> %ON_OTHER_INTERACT +# (15) PMIDs +# (16) LT_LIT +# (17) MS_LIT +# (18) MS_CST +# (19) NOTES --> %notes +############################################################################################################################### + + +$dbh = DBI->connect("dbi:SQLite:$db_out", undef, undef); +my $auto_commit = $dbh->{AutoCommit}; +$dbh->{AutoCommit} = 0; +print "DB connection $dbh is to $db_out, opened for modification\n"; + +# add partial PSP_Regulatory_site table (if not exists) regardless of whether SwissProt input was FASTA or SQLite +$stmth = $dbh->prepare(" +CREATE TABLE IF NOT EXISTS PSP_Regulatory_site ( + SITE_PLUSMINUS_7AA TEXT PRIMARY KEY ON CONFLICT IGNORE, + DOMAIN TEXT, + ON_FUNCTION TEXT, + ON_PROCESS TEXT, + ON_PROT_INTERACT TEXT, + ON_OTHER_INTERACT TEXT, + NOTES TEXT, + ORGANISM TEXT, + PROTEIN TEXT +) +"); +$stmth->execute(); + +# add partial PSP_Regulatory_site LUT (if not exists) regardless of whether SwissProt input was FASTA or SQLite +$stmth = $dbh->prepare(" +CREATE TABLE IF NOT EXISTS ppep_regsite_LUT +( ppep_id INTEGER REFERENCES ppep(id) +, site_plusminus_7AA TEXT REFERENCES PSP_Regulatory_site(site_plusminus_7AA) +, PRIMARY KEY (ppep_id, site_plusminus_7AA) ON CONFLICT IGNORE +); +"); +$stmth->execute(); + +# $stmth = $dbh->prepare(" +# CREATE UNIQUE INDEX idx_PSP_Regulatory_site_0 +# ON PSP_Regulatory_site(site_plusminus_7AA); +# "); +# $stmth->execute(); + + +# add Citation table (if not exists) regardless of whether SwissProt input was FASTA or SQLite +my $citation_sql; +$citation_sql = " +CREATE TABLE IF NOT EXISTS Citation ( + ObjectName TEXT REFERENCES sqlite_schema(name) ON DELETE CASCADE, + CitationData TEXT, + PRIMARY KEY (ObjectName, CitationData) ON CONFLICT IGNORE +) +"; +$stmth = $dbh->prepare($citation_sql); +$stmth->execute(); + + +open (IN4, "$PSP_Regulatory_Sites_in") or die "I couldn't find $PSP_Regulatory_Sites_in\n"; +print "Reading the PhosphoSite regulatory site data: $PSP_Regulatory_Sites_in\n"; + + +$line = -1; +while (<IN4>) { + $line++; + chomp; + if ($_ =~ m/PhosphoSitePlus/) { + #$PhosphoSitePlusCitation = ($_ =~ s/PhosphoSitePlus/FooBar/g); + $PhosphoSitePlusCitation = $_; + $PhosphoSitePlusCitation =~ s/\t//g; + $PhosphoSitePlusCitation =~ s/\r//g; + $PhosphoSitePlusCitation =~ s/\n//g; + $PhosphoSitePlusCitation =~ s/""/"/g; + $PhosphoSitePlusCitation =~ s/^"//g; + $PhosphoSitePlusCitation =~ s/"$//g; + print "$PhosphoSitePlusCitation\n"; + next; + } + my (@x) = split(/\t/); + for my $i (0 .. $#x) { + $x[$i] =~ s/\r//g; $x[$i] =~ s/\n//g; $x[$i] =~ s/\"//g; + } + my $found_GENE=0; + if ( (not exists($x[0])) ) { + next; + } + elsif ( ($x[0] eq "GENE") ) { + $found_GENE=1; + next; + } + if ( (not exists($x[9])) || ($x[9] eq "") ) { + if (exists($x[8]) && (not $x[8] eq "")) { + die "$PSP_Regulatory_Sites_in line $line has no SITE_+/-7_AA: $_\n"; + } else { + if ( (not exists($x[1])) || (not $x[1] eq "") ) { + print "$PSP_Regulatory_Sites_in line $line (".length($_)." characters) has no SITE_+/-7_AA: $_\n" + if $found_GENE==1; + } + next; + } + } + elsif ($line != 0) { + if ($species ne $x[6]) { + # Do nothing - this record was filtered out by the species filter + } + elsif (!exists($regulatory_sites_PhosphoSite_hash{$x[9]})) { + if (!defined $domain{$x[9]} || $domain{$x[9]} eq "") { + $regulatory_sites_PhosphoSite_hash{$x[9]} = $x[9]; + $domain{$x[9]} = $x[10]; + $ON_FUNCTION{$x[9]} = $x[11]; + $ON_PROCESS{$x[9]} = $x[12]; + $ON_PROT_INTERACT{$x[9]} = $x[13]; + $ON_OTHER_INTERACT{$x[9]} = $x[14]; + $notes{$x[9]} = $x[19]; + $organism{$x[9]} = $x[6]; + } + } + else { + # $domain + if (!defined $domain{$x[9]} || $domain{$x[9]} eq "") { + if ($x[10] ne "") { + $domain{$x[9]} = $domain{$x[10]}; + } + else { + # do nothing + } + } + else { + if ($domain{$x[9]} =~ /$x[10]/) { + # do nothing + } + else { + $domain{$x[9]} = $domain{$x[9]}." / ".$x[10]; + #print "INFO line $line - compound domain for 7aa: GENE $x[0] PROTEIN $x[1] PROT_TYPE $x[2] ACC_ID $x[3] GENE_ID $x[4] HU_CHR_LOC $x[5] ORGANISM $x[6] MOD_RSD $x[7] SITE_GRP_ID $x[8] SITE_+/-7_AA $x[9] DOMAIN $domain{$x[9]}\n"; + } + } + + # $ON_FUNCTION + if (!defined $ON_FUNCTION{$x[9]} || $ON_FUNCTION{$x[9]} eq "") { + $ON_FUNCTION{$x[9]} = $ON_FUNCTION{$x[10]}; + } elsif ($x[10] eq "") { + # do nothing + } + else { + $ON_FUNCTION{$x[9]} = $ON_FUNCTION{$x[9]}." / ".$x[10]; + } + + # $ON_PROCESS + if (!defined $ON_PROCESS{$x[9]} || $ON_PROCESS{$x[9]} eq "") { + $ON_PROCESS{$x[9]} = $ON_PROCESS{$x[10]}; + } elsif ($x[10] eq "") { + # do nothing + } + else { + $ON_PROCESS{$x[9]} = $ON_PROCESS{$x[9]}." / ".$x[10]; + } + + # $ON_PROT_INTERACT + if (!defined $ON_PROT_INTERACT{$x[9]} || $ON_PROT_INTERACT{$x[9]} eq "") { + $ON_PROT_INTERACT{$x[9]} = $ON_PROT_INTERACT{$x[10]}; + } elsif ($x[10] eq "") { + # do nothing + } + else { + $ON_PROT_INTERACT{$x[9]} = $ON_PROT_INTERACT{$x[9]}." / ".$x[10]; + } + + # $ON_OTHER_INTERACT + if (!defined $ON_OTHER_INTERACT{$x[9]} || $ON_OTHER_INTERACT{$x[9]} eq "") { + $ON_OTHER_INTERACT{$x[9]} = $ON_OTHER_INTERACT{$x[10]}; + } elsif ($x[10] eq "") { + # do nothing + } + else { + $ON_OTHER_INTERACT{$x[9]} = $ON_OTHER_INTERACT{$x[9]}." / ".$x[10]; + } + + # $notes + if (!defined $notes{$x[9]} || $notes{$x[9]} eq "") { + $notes{$x[9]} = $notes{$x[10]}; + } elsif ($x[10] eq "") { + # do nothing + } + else { + $notes{$x[9]} = $notes{$x[9]}." / ".$x[10]; + } + + # $organism + if (!defined $organism{$x[9]} || $organism{$x[9]} eq "") { + $organism{$x[9]} = $organism{$x[10]}; + } elsif ($x[10] eq "") { + # do nothing + } + else { + $organism{$x[9]} = $organism{$x[9]}." / ".$x[10]; + } + } + } +} +close IN4; + +print "... Finished reading various site data at " . format_localtime_iso8601() ."\n\n"; + +$stmth = $dbh->prepare(" +INSERT INTO Citation ( + ObjectName, + CitationData +) VALUES (?,?) +"); + +sub add_citation { + my ($cit_table, $cit_text, $cit_label) = @_; + $stmth->bind_param(1, $cit_table); + $stmth->bind_param(2, $cit_text); + if (not $stmth->execute()) { + print "Error writing $cit_label cit for table $cit_table: $stmth->errstr\n"; + } +} +my ($citation_text, $citation_table); + +# PSP regulatory or kinase/substrate site +$citation_text = 'PhosphoSitePlus(R) (PSP) was created by Cell Signaling Technology Inc. It is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License. When using PSP data or analyses in printed publications or in online resources, the following acknowledgements must be included: (a) the words "PhosphoSitePlus(R), www.phosphosite.org" must be included at appropriate places in the text or webpage, and (b) the following citation must be included in the bibliography: "Hornbeck PV, Zhang B, Murray B, Kornhauser JM, Latham V, Skrzypek E PhosphoSitePlus, 2014: mutations, PTMs and recalibrations. Nucleic Acids Res. 2015 43:D512-20. PMID: 25514926."'; +$citation_table = "PSP_Regulatory_site"; +add_citation($citation_table, $citation_text, "PSP_Kinase_Substrate"); +$citation_table = "psp_gene_site"; +add_citation($citation_table, $citation_text, "PSP_Kinase_Substrate"); +$citation_table = "psp_gene_site_view"; +add_citation($citation_table, $citation_text, "PSP_Regulatory_site"); +$citation_text = 'Hornbeck, 2014, "PhosphoSitePlus, 2014: mutations, PTMs and recalibrations.", https://pubmed.ncbi.nlm.nih.gov/22135298, https://doi.org/10.1093/nar/gkr1122'; +$citation_table = "PSP_Regulatory_site"; +add_citation($citation_table, $citation_text, "PSP_Regulatory_site"); +$citation_table = "psp_gene_site"; +add_citation($citation_table, $citation_text, "PSP_Kinase_Substrate"); +$citation_table = "psp_gene_site_view"; +add_citation($citation_table, $citation_text, "PSP_Kinase_Substrate"); + +# NetworKIN site +$citation_text = 'Linding, 2007, "Systematic discovery of in vivo phosphorylation networks.", https://pubmed.ncbi.nlm.nih.gov/17570479, https://doi.org/10.1016/j.cell.2007.05.052'; +$citation_table = "psp_gene_site"; +add_citation($citation_table, $citation_text, "NetworkKIN"); +$citation_table = "psp_gene_site_view"; +add_citation($citation_table, $citation_text, "NetworkKIN"); +$citation_text = 'Horn, 2014, "KinomeXplorer: an integrated platform for kinome biology studies.", https://pubmed.ncbi.nlm.nih.gov/24874572, https://doi.org/10.1038/nmeth.296'; +$citation_table = "psp_gene_site"; +add_citation($citation_table, $citation_text, "NetworkKIN"); +$citation_table = "psp_gene_site_view"; +add_citation($citation_table, $citation_text, "NetworkKIN"); +$citation_text = 'Aken, 2016, "The Ensembl gene annotation system.", https://pubmed.ncbi.nlm.nih.gov/33137190, https://doi.org/10.1093/database/baw093'; +$citation_table = "psp_gene_site"; +add_citation($citation_table, $citation_text, "NetworkKIN"); +$citation_table = "psp_gene_site_view"; +add_citation($citation_table, $citation_text, "NetworkKIN"); + +# pSTY motifs +$citation_text = 'Amanchy, 2007, "A curated compendium of phosphorylation motifs.", https://pubmed.ncbi.nlm.nih.gov/17344875, https://doi.org/10.1038/nbt0307-285'; +$citation_table = "psp_gene_site"; +add_citation($citation_table, $citation_text, "Amanchy_pSTY_motifs"); +$citation_table = "psp_gene_site_view"; +add_citation($citation_table, $citation_text, "Amanchy_pSTY_motifs"); +$citation_text = 'Gnad, 2011, "PHOSIDA 2011: the posttranslational modification database.", https://pubmed.ncbi.nlm.nih.gov/21081558, https://doi.org/10.1093/nar/gkq1159'; +$citation_table = "psp_gene_site"; +add_citation($citation_table, $citation_text, "Phosida_pSTY_motifs"); +$citation_table = "psp_gene_site_view"; +add_citation($citation_table, $citation_text, "Phosida_pSTY_motifs"); + + +############################################################################################################################### +# +# Read the data file: +# 1) find sequences that match the NetworKIN predictions +# 2) find motifs that match the observed sequences +# +############################################################################################################################### + +print "--- Find sequences that match the NetworKIN predictions and find motifs that match observed sequences\n"; + +my $ppep_regsite_LUT_stmth; +$ppep_regsite_LUT_stmth = $dbh->prepare(" + INSERT INTO ppep_regsite_LUT ( + ppep_id, + site_plusminus_7AA + ) VALUES (?,?) +"); + +my ($start_seconds, $start_microseconds) = gettimeofday; + +foreach my $peptide (keys %data) { + # find the unique phospho-motifs for this $peptide + my @all_motifs = (); + my $have_all_motifs = 0; + for my $i (0 .. $#{ $matched_sequences{$peptide} } ) { + my $tmp_motif = $p_motifs{$peptide}[$i]; + push(@all_motifs, $tmp_motif); + $have_all_motifs = 1; + } + if ($have_all_motifs == 1) { + for my $j (0 .. $#all_motifs) { + if (defined $all_motifs[$j]) { + $all_motifs[$j] =~ s/\d+-\[\s//; + $all_motifs[$j] =~ s/\s\]\-\d+//; + } + } + } + my %seen = (); + if ($have_all_motifs == 1) { + foreach my $a (@all_motifs) { + if (defined $a) { + if (exists($seen{$a})) { + next; + } else { + push(@{$unique_motifs{$peptide}}, $a); + $seen{$a} = 1; + } + } + print "push(\@{\$unique_motifs{$peptide}}, $a);\n" if ($verbose); + } + } + + # count the number of phospo-sites in the motif + my $number_pY = 0; + my $number_pSTY = 0; + if ($phospho_type eq 'y') { + if (defined(${$unique_motifs{$peptide}}[0])) { + while (${$unique_motifs{$peptide}}[0] =~ /pY/g) { + $number_pY++; + } + } + } + if ($phospho_type eq 'sty') { + print "looking for unique_motifs for $peptide\n" if ($verbose); + if (defined(${$unique_motifs{$peptide}}[0])) { + while (${$unique_motifs{$peptide}}[0] =~ /(pS|pT|pY)/g) { + $number_pSTY++; + print "We have found $number_pSTY unique_motifs for $peptide\n" if ($verbose); + } + } + } + + + # search each of the unique motifs for matches + print "searching $#{$unique_motifs{$peptide}} motifs for peptide $peptide\n" if ($verbose); + for my $i (0 .. $#{$unique_motifs{$peptide}}) { + print "\$i = $i; peptide = $peptide; unique_motif = ${$unique_motifs{$peptide}}[$i]\n" if ($verbose); + my $tmp_motif = ${$unique_motifs{$peptide}}[$i]; + print " --- matching unique motif $tmp_motif for peptide $peptide at " . format_localtime_iso8601() ."\n" if ($verbose); + my $formatted_sequence; + if (($number_pY == 1) || ($number_pSTY == 1)) { + my $seq_plus5aa = ""; + my $seq_plus7aa = ""; + $formatted_sequence = &replace_pSpTpY($tmp_motif, $phospho_type); + print " a #pY $number_pY; #pSTY $number_pSTY; matching formatted motif $formatted_sequence for peptide $peptide at " . format_localtime_iso8601() ."\n" if ($verbose); + if ($phospho_type eq 'y') { + $seq_plus5aa = (split(/(\w{0,5}y\w{0,5})/, $formatted_sequence))[1]; + $seq_plus7aa = (split(/(\w{0,7}y\w{0,7})/, $formatted_sequence))[1]; + } + elsif ($phospho_type eq "sty") { + $seq_plus5aa = (split(/(\w{0,5}(s|t|y)\w{0,5})/, $formatted_sequence))[1]; + $seq_plus7aa = (split(/(\w{0,7}(s|t|y)\w{0,7})/, $formatted_sequence))[1]; + } + + if (defined $seq_plus7aa) { + # commit the 7aa LUT records + $ppep_regsite_LUT_stmth->bind_param( 1, $ppep_id_lut{$peptide} ); + $ppep_regsite_LUT_stmth->bind_param( 2, $seq_plus7aa ); + if (not $ppep_regsite_LUT_stmth->execute()) { + print "Error writing tuple ($ppep_id_lut{$peptide},$seq_plus7aa) for peptide $peptide to ppep_regsite_LUT: $ppep_regsite_LUT_stmth->errstr\n"; + } + } + for my $i (0 .. $#kinases_observed) { + if (defined $seq_plus5aa) { + my $tmp = $seq_plus5aa."_".$kinases_observed[$i]; #eg, should be PGRPLsSYGMD_PKCalpha + if (exists($p_sequence_kinase -> {$tmp})) { + $kinase_substrate_NetworKIN_matches{$peptide}{$kinases_observed[$i]} = "X"; #ACE + } + } + } + for my $i (0 .. $#motif_sequence) { + print "matching $motif_sequence[$i]" if ($verbose); + if ($peptide =~ /$motif_sequence[$i]/) { + $kinase_motif_matches{$peptide}{$motif_type{$motif_type_key_ary[$i]}} = "X"; + } + } + for my $i (0 .. $#kinases_PhosphoSite) { + if (defined $seq_plus7aa) { + my $tmp = $seq_plus7aa."_".$kinases_PhosphoSite[$i]; #eg, should be RTPGRPLsSYGMDSR_PAK2 + if (exists($p_sequence_kinase_PhosphoSite -> {$tmp})) { + $kinase_substrate_PhosphoSite_matches{$peptide}{$kinases_PhosphoSite[$i]} = "X"; + } + } + } + if (exists($regulatory_sites_PhosphoSite_hash{$seq_plus7aa})) { + $seq_plus7aa_2{$peptide} = $seq_plus7aa; + $domain_2{$peptide} = $domain{$seq_plus7aa}; + $ON_FUNCTION_2{$peptide} = $ON_FUNCTION{$seq_plus7aa}; + $ON_PROCESS_2{$peptide} = $ON_PROCESS{$seq_plus7aa}; + $ON_PROT_INTERACT_2{$peptide} = $ON_PROT_INTERACT{$seq_plus7aa}; + $ON_OTHER_INTERACT_2{$peptide} = $ON_OTHER_INTERACT{$seq_plus7aa}; + $notes_2{$peptide} = $notes{$seq_plus7aa}; + $organism_2{$peptide} = $organism{$seq_plus7aa}; + } else { + } + } + elsif (($number_pY > 1) || ($number_pSTY > 1)) { #eg, if $x[4] is 1308-[ VIYFQAIEEVpYpYDHLRSAAKKR ]-1329 and $number_pY == 2 + $formatted_sequence = $tmp_motif; + $seq_plus5aa = ""; + $seq_plus7aa = ""; + #Create the sequences with only one phosphorylation site + #eg, 1308-[ VIYFQAIEEVpYpYDHLRSAAKKR ]-1329, which becomes 1308-[ VIYFQAIEEVpYYDHLRSAAKKR ]-1329 and 1308-[ VIYFQAIEEVYpYDHLRSAAKKR ]-1329 + + my (@sites, $offset, $next_p_site); + $sites[0] = index($tmp_motif, "p"); + $offset = $sites[0] + 1; + $next_p_site = 0; + while ($next_p_site != -1) { + $next_p_site = index($tmp_motif, "p", $offset); + if ($next_p_site != -1) { + push (@sites, $next_p_site); + } + $offset = $next_p_site+1; + } + + my @pSTY_sequences; + for my $n (0 .. $#sites) { + $pSTY_sequences[$n] = $tmp_motif; + for (my $m = $#sites; $m >= 0; $m--) { + if ($m != $n) {substr($pSTY_sequences[$n], $sites[$m], 1) = "";} + } + } + + my @formatted_sequences; + for my $k (0 .. $#sites) { + $formatted_sequences[$k] = &replace_pSpTpY($pSTY_sequences[$k], $phospho_type); + } + + for my $k (0 .. $#formatted_sequences) { + print " b #pY $number_pY; #pSTY $number_pSTY; matching formatted motif $formatted_sequences[$k] for peptide $peptide at " . format_localtime_iso8601() ."\n" if ($verbose); + if ($phospho_type eq 'y') { + $seq_plus5aa = (split(/(\w{0,5}y\w{0,5})/, $formatted_sequences[$k]))[1]; + $seq_plus7aa = (split(/(\w{0,7}y\w{0,7})/, $formatted_sequences[$k]))[1]; + } + elsif ($phospho_type eq "sty") { + $seq_plus5aa = (split(/(\w{0,5}(s|t|y)\w{0,5})/, $formatted_sequences[$k]))[1]; + $seq_plus7aa = (split(/(\w{0,7}(s|t|y)\w{0,7})/, $formatted_sequences[$k]))[1]; + } + for my $i (0 .. $#kinases_observed) { + my $tmp = $seq_plus5aa."_".$kinases_observed[$i]; #eg, should look like REEILsEMKKV_PKCalpha + if (exists($p_sequence_kinase -> {$tmp})) { + $kinase_substrate_NetworKIN_matches{$peptide}{$kinases_observed[$i]} = "X"; + } + } + $pSTY_sequence = $formatted_sequences[$k]; + for my $i (0 .. $#motif_sequence) { + if ($pSTY_sequence =~ /$motif_sequence[$i]/) { + $kinase_motif_matches{$peptide}{$motif_type{$motif_type_key_ary[$i]}} = "X"; + } + } + for my $i (0 .. $#kinases_PhosphoSite) { + my $tmp = $seq_plus7aa."_".$kinases_PhosphoSite[$i]; #eg, should be RTPGRPLsSYGMDSR_PAK2 + #print "seq_plus7aa._.kinases_PhosphoSite[i] is $tmp"; + if (exists($p_sequence_kinase_PhosphoSite -> {$tmp})) { + $kinase_substrate_PhosphoSite_matches{$peptide}{$kinases_PhosphoSite[$i]} = "X"; + } + } + if (exists($regulatory_sites_PhosphoSite -> {$seq_plus7aa})) { + $seq_plus7aa_2{$peptide} = $seq_plus7aa; + + # $domain + if ($domain_2{$peptide} eq "") { + $domain_2{$peptide} = $domain{$seq_plus7aa}; + } + elsif ($domain{$seq_plus7aa} eq "") { + # do nothing + } + else { + $domain_2{$peptide} = $domain_2{$peptide}." / ".$domain{$seq_plus7aa}; + } + + + # $ON_FUNCTION_2 + if ($ON_FUNCTION_2{$peptide} eq "") { + $ON_FUNCTION_2{$peptide} = $ON_FUNCTION{$seq_plus7aa}; + } + elsif ($ON_FUNCTION{$seq_plus7aa} eq "") { + # do nothing + } + else { + $ON_FUNCTION_2{$peptide} = $ON_FUNCTION_2{$peptide}." / ".$ON_FUNCTION{$seq_plus7aa}; + } + + # $ON_PROCESS_2 + if ($ON_PROCESS_2{$peptide} eq "") { + $ON_PROCESS_2{$peptide} = $ON_PROCESS{$seq_plus7aa}; + } + elsif ($ON_PROCESS{$seq_plus7aa} eq "") { + # do nothing + } + else { + $ON_PROCESS_2{$peptide} = $ON_PROCESS_2{$peptide}." / ".$ON_PROCESS{$seq_plus7aa}; + } + + # $ON_PROT_INTERACT_2 + if ($ON_PROT_INTERACT_2{$peptide} eq "") { + $ON_PROT_INTERACT_2{$peptide} = $ON_PROT_INTERACT{$seq_plus7aa}; + } + elsif ($ON_PROT_INTERACT{$seq_plus7aa} eq "") { + # do nothing + } + else { + $ON_PROT_INTERACT_2{$peptide} = $ON_PROT_INTERACT_2{$peptide}." / ".$ON_PROT_INTERACT{$seq_plus7aa}; + } + + # $ON_OTHER_INTERACT_2 + if ($ON_OTHER_INTERACT_2{$peptide} eq "") { + $ON_OTHER_INTERACT_2{$peptide} = $ON_OTHER_INTERACT{$seq_plus7aa}; + } + elsif ($ON_OTHER_INTERACT{$seq_plus7aa} eq "") { + # do nothing + } + else { + $ON_OTHER_INTERACT_2{$peptide} = $ON_OTHER_INTERACT_2{$peptide}." / ".$ON_OTHER_INTERACT{$seq_plus7aa}; + } + + # $notes_2 + if ($notes_2{$peptide} eq "") { + $notes_2{$peptide} = $notes{$seq_plus7aa}; + } + elsif ($notes{$seq_plus7aa} eq "") { + # do nothing + } + else { + $notes_2{$peptide} = $notes_2{$peptide}." / ".$notes{$seq_plus7aa}; + } + $notes_2{$peptide} = $notes{$seq_plus7aa}; + + # $organism_2 + if ($organism_2{$peptide} eq "") { + $organism_2{$peptide} = $organism{$seq_plus7aa}; + } + elsif ($organism{$seq_plus7aa} eq "") { + # do nothing + } + else { + $organism_2{$peptide} = $organism_2{$peptide}." / ".$organism{$seq_plus7aa}; + } + $organism_2{$peptide} = $organism{$seq_plus7aa}; + } else { + } # if (exists($regulatory_sites_PhosphoSite -> {$seq_plus7aa})) + } # for my $k (0 .. $#formatted_sequences) + } # if/else number of phosphosites + } # for each motif i # for my $i (0 .. $#{$unique_motifs{$peptide}}) +} # for each $peptide + +my ($end_seconds, $end_microseconds) = gettimeofday; + +my $delta_seconds = $end_seconds - $start_seconds; +my $delta_microseconds = $end_microseconds - $start_microseconds; +$delta_microseconds += 1000000 * $delta_seconds; +my $key_count = keys(%data); +print sprintf("Average search time is %d microseconds per phopshopeptide\n", ($delta_microseconds / $key_count)); + +($start_seconds, $start_microseconds) = gettimeofday; + +print "Writing PSP_Regulatory_site records\n"; + +my $psp_regulatory_site_stmth = $dbh->prepare(" + INSERT INTO PSP_Regulatory_site ( + DOMAIN, + ON_FUNCTION, + ON_PROCESS, + ON_PROT_INTERACT, + ON_OTHER_INTERACT, + NOTES, + SITE_PLUSMINUS_7AA, + ORGANISM + ) VALUES (?,?,?,?,?,?,?,?) + "); + +foreach my $peptide (keys %data) { + if (exists($domain_2{$peptide}) and (defined $domain_2{$peptide}) and (not $domain_2{$peptide} eq "") ) { + $psp_regulatory_site_stmth->bind_param(1, $domain_2{$peptide}); + $psp_regulatory_site_stmth->bind_param(2, $ON_FUNCTION_2{$peptide}); + $psp_regulatory_site_stmth->bind_param(3, $ON_PROCESS_2{$peptide}); + $psp_regulatory_site_stmth->bind_param(4, $ON_PROT_INTERACT_2{$peptide}); + $psp_regulatory_site_stmth->bind_param(5, $ON_OTHER_INTERACT_2{$peptide}); + $psp_regulatory_site_stmth->bind_param(6, $notes_2{$peptide}); + $psp_regulatory_site_stmth->bind_param(7, $seq_plus7aa_2{$peptide}); + $psp_regulatory_site_stmth->bind_param(8, $organism_2{$peptide}); + if (not $psp_regulatory_site_stmth->execute()) { + print "Error writing PSP_Regulatory_site for one regulatory site with peptide '$domain_2{$peptide}': $psp_regulatory_site_stmth->errstr\n"; + } else { + } + } elsif (exists($domain_2{$peptide}) and (not defined $domain_2{$peptide})) { + print "\$domain_2{$peptide} is undefined\n"; #ACE + } +} + +$dbh->{AutoCommit} = $auto_commit; +# auto_commit implicitly finishes psp_regulatory_site_stmth, apparently # $psp_regulatory_site_stmth->finish; +$dbh->disconnect if ( defined $dbh ); + + +($end_seconds, $end_microseconds) = gettimeofday; + +$delta_seconds = $end_seconds - $start_seconds; +$delta_microseconds = $end_microseconds - $start_microseconds; +$delta_microseconds += 1000000 * $delta_seconds; +$key_count = keys(%data); +print sprintf("Write time is %d microseconds\n", ($delta_microseconds)); + +print "... Finished find sequences that match the NetworKIN predictions and find motifs that match observed sequences at " . format_localtime_iso8601() ."\n\n"; + +############################################################################################################################### +# +# Print to the output file +# +############################################################################################################################### + + +open (OUT, ">$file_out") || die "could not open the fileout: $file_out"; +open (MELT, ">$file_melt") || die "could not open the fileout: $file_melt"; + +# print the header info +print MELT "phospho_peptide\tgene_names\tsite_type\tkinase_map\n"; +print OUT "p-peptide\tProtein description\tGene name(s)\tFASTA name\tPhospho-sites\tUnique phospho-motifs, no residue numbers\tAccessions\tPhospho-motifs for all members of protein group with residue numbers\t"; + +# print the PhosphoSite regulatory data +print OUT "Domain\tON_FUNCTION\tON_PROCESS\tON_PROT_INTERACT\tON_OTHER_INTERACT\tPhosphoSite notes\t"; + +# print the sample names +for my $i (0 .. $#samples) { print OUT "$samples[$i]\t"; } + +# print the kinases and groups +for my $i (0 .. $#kinases_observed) { + my $temp = $kinases_observed[$i]."_NetworKIN"; + print OUT "$temp\t"; + push(@kinases_observed_lbl, $temp); +} +my @motif_type_keys = keys %motif_type; +for my $i (1 .. $#motif_type_keys) { + print OUT "$motif_type{$motif_type_keys[$i]}\t"; +} +for my $i (0 .. $#kinases_PhosphoSite) { + my $temp = $kinases_PhosphoSite[$i]; # ."_PhosphoSite"; + if ($i < $#kinases_PhosphoSite) { print OUT "$temp\t"; } + if ($i == $#kinases_PhosphoSite) { print OUT "$temp\n"; } + push(@phosphosites_observed_lbl, $temp); +} + +# begin DDL-to-SQLite +# --- +$dbh = DBI->connect("dbi:SQLite:$db_out", undef, undef); +$auto_commit = $dbh->{AutoCommit}; +$dbh->{AutoCommit} = 0; +print "DB connection $dbh is to $db_out, opened for modification\n"; + +my $sample_stmth; +$sample_stmth = $dbh->prepare(" + INSERT INTO sample ( + id, + name + ) VALUES (?,?) +"); + +my $ppep_intensity_stmth; +$ppep_intensity_stmth = $dbh->prepare(" + INSERT INTO ppep_intensity ( + ppep_id, + sample_id, + intensity + ) VALUES (?,?,?) +"); + +my $site_type_stmth; +$site_type_stmth = $dbh->prepare(" + insert into site_type ( + id, + type_name + ) values (?,?) +"); + +my $ppep_gene_site_stmth; +$ppep_gene_site_stmth = $dbh->prepare(" + insert into ppep_gene_site ( + ppep_id, + gene_names, + kinase_map, + site_type_id + ) values (?,?,?,?) +"); + +my $ppep_metadata_stmth; +$ppep_metadata_stmth = $dbh->prepare(" + INSERT INTO ppep_metadata + ( ppep_id + , protein_description + , gene_name + , FASTA_name + , phospho_sites + , motifs_unique + , accessions + , motifs_all_members + , domain + , ON_FUNCTION + , ON_PROCESS + , ON_PROT_INTERACT + , ON_OTHER_INTERACT + , notes + ) VALUES ( + ?,?,?,?,?,?,? + , ?,?,?,?,?,?,? + ) +"); +# end DDL-to-SQLite +# ... + +# begin store-to-SQLite "sample" table +# --- +# %sample_id_lut maps name -> ID +for my $sample_name (keys %sample_id_lut) { + $sample_stmth->bind_param( 2, $sample_name ); + $sample_stmth->bind_param( 1, $sample_id_lut{$sample_name} ); + if (not $sample_stmth->execute()) { + print "Error writing tuple ($sample_name,$sample_id_lut{$sample_name}): $sample_stmth->errstr\n"; + } +} +# end store-to-SQLite "sample" table +# ... + +# begin store-to-SQLite "site_type" table +# --- +sub add_site_type { + my ($site_type_id, $site_type_type_name) = @_; + $site_type_stmth->bind_param( 2, $site_type_type_name ); + $site_type_stmth->bind_param( 1, $site_type_id ); + if (not $site_type_stmth->execute()) { + die "Error writing tuple ($site_type_id,$site_type_type_name): $site_type_stmth->errstr\n"; + } +} +add_site_type($SITE_KINASE_SUBSTRATE, $site_description{$SITE_KINASE_SUBSTRATE}); +add_site_type($SITE_HPRD , $site_description{$SITE_HPRD }); +add_site_type($SITE_PHOSIDA , $site_description{$SITE_PHOSIDA }); +add_site_type($SITE_PHOSPHOSITE , $site_description{$SITE_PHOSPHOSITE }); +# end store-to-SQLite "site_type" table +# ... + +foreach my $peptide (sort(keys %data)) { + next if (grep($peptide, @failed_matches)); + my $ppep_id = $ppep_id_lut{$peptide}; + my @ppep_metadata = (); + my @ppep_intensity = (); + my @gene = (); + my $gene_names; + my $j; + # Print the peptide itself + # column 1: p-peptide + print OUT "$peptide\t"; + push (@ppep_metadata, $ppep_id); + push (@ppep_intensity, $peptide); + + my $verbose_cond = 0; # $peptide eq 'AAAAAAAGDpSDpSWDADAFSVEDPVR' || $peptide eq 'KKGGpSpSDEGPEPEAEEpSDLDSGSVHSASGRPDGPVR'; + # skip over failed matches + print "\nfirst match for '$peptide' is '$matched_sequences{$peptide}[0]' and FAILED_MATCH_SEQ is '$FAILED_MATCH_SEQ'\n" if $verbose_cond; + if ($matched_sequences{$peptide}[0] eq $FAILED_MATCH_SEQ) { + # column 2: Protein description + # column 3: Gene name(s) + # column 4: FASTA name + # column 5: phospho-residues + # Column 6: UNIQUE phospho-motifs + # Column 7: accessions + # Column 8: ALL motifs with residue numbers + # 2 3 4 5 6 7 8 + print OUT "Sequence not found in FASTA database\tNA\tNA\tNA\tNA\tNA\tNA\t"; + print "No match found for '$peptide' in sequence database\n"; + $gene_names = '$FAILED_MATCH_GENE_NAME'; + } else { + my @description = (); + my %seen = (); + # Print just the protein description + for $i (0 .. $#{$names{$peptide}}) { + my $long_name = $names{$peptide}[$i]; + my @naming_parts = split(/\sOS/, $long_name); + my @front_half = split(/\s/, $naming_parts[0]); + push(@description, join(" ", @front_half[1..($#front_half)])); + } + # column 2: Protein description + print OUT join(" /// ", @description), "\t"; + push (@ppep_metadata, join(" /// ", @description)); + + # Print just the gene name + for $i (0 .. $#{$names{$peptide}}) { + my $tmp_gene = $names{$peptide}[$i]; + $tmp_gene =~ s/^.*GN=//; + $tmp_gene =~ s/\s.*//; + if (!exists($seen{$tmp_gene})) { + push(@gene, $tmp_gene); + $seen{$tmp_gene} = $tmp_gene; + } + } + # column 3: Gene name(s) + $gene_names = join(" /// ", @gene); + print OUT $gene_names, "\t"; + push (@ppep_metadata, join(" /// ", @gene)); + + # column 4: FASTA name + print OUT join(" /// ", @{$names{$peptide}}), "\t"; + push (@ppep_metadata, join(" /// ", @{$names{$peptide}})); + + # column 5: phospho-residues + my $tmp_for_insert = ""; + my $foobar; + for my $i (0 .. $#{ $matched_sequences{$peptide} } ) { + print "match $i for '$peptide' is '$matched_sequences{$peptide}[$i]'\n" if $verbose_cond; + if ($i < $#{ $matched_sequences{$peptide} }) { + if (defined $p_residues{$peptide}{$i}) { + @tmp_p_residues = @{$p_residues{$peptide}{$i}}; + for $j (0 .. $#tmp_p_residues) { + if ($j < $#tmp_p_residues) { + my $tmp_site_for_printing = $p_residues{$peptide}{$i}[$j] + 1; # added 12.05.2012 for Justin's data + print OUT "p$residues{$peptide}{$i}[$j]$tmp_site_for_printing, "; + $tmp_for_insert .= "p$residues{$peptide}{$i}[$j]$tmp_site_for_printing, "; + } + elsif ($j == $#tmp_p_residues) { + my $tmp_site_for_printing = $p_residues{$peptide}{$i}[$j] + 1; # added 12.05.2012 for Justin's data + print OUT "p$residues{$peptide}{$i}[$j]$tmp_site_for_printing /// "; + $tmp_for_insert .= "p$residues{$peptide}{$i}[$j]$tmp_site_for_printing /// "; + } + } + } + } + elsif ($i == $#{ $matched_sequences{$peptide} }) { + my $there_were_sites = 0; + if (defined $p_residues{$peptide}{$i}) { + @tmp_p_residues = @{$p_residues{$peptide}{$i}}; + if ($#tmp_p_residues > 0) { + for my $j (0 .. $#tmp_p_residues) { + if ($j < $#tmp_p_residues) { + if (defined $p_residues{$peptide}{$i}[$j]) { + my $tmp_site_for_printing = $p_residues{$peptide}{$i}[$j] + 1; # added 12.05.2012 for Justin's data + $foobar = $residues{$peptide}{$i}[$j]; + if (defined $foobar) { + print OUT "$foobar"; + print OUT "$tmp_site_for_printing, "; + $tmp_for_insert .= "p$residues{$peptide}{$i}[$j]$tmp_site_for_printing, "; + $there_were_sites = 1; + } + } + } + elsif ($j == $#tmp_p_residues) { + if (defined $p_residues{$peptide}{$i}[$j]) { + $foobar = $residues{$peptide}{$i}[$j]; + if (defined $foobar) { + my $tmp_site_for_printing = $p_residues{$peptide}{$i}[$j] + 1; # added 12.05.2012 for Justin's data + print OUT "$foobar"; + print OUT "$tmp_site_for_printing\t"; + $tmp_for_insert .= "p$residues{$peptide}{$i}[$j]$tmp_site_for_printing"; + $there_were_sites = 1; + } + } + } + } + } + } + if (0 == $there_were_sites) { + print OUT "\t"; + } + } + } + print "tmp_for_insert '$tmp_for_insert' for '$peptide'\n" if $verbose_cond; + push (@ppep_metadata, $tmp_for_insert); + + # Column 6: UNIQUE phospho-motifs + print OUT join(" /// ", @{$unique_motifs{$peptide}}), "\t"; + push (@ppep_metadata, join(" /// ", @{$unique_motifs{$peptide}})); + + # Column 7: accessions + if (defined $accessions{$peptide}) { + print OUT join(" /// ", @{$accessions{$peptide}}), "\t"; + push (@ppep_metadata, join(" /// ", @{$accessions{$peptide}})); + } else { + print OUT "\t"; + push (@ppep_metadata, ""); + } + + # Column 8: ALL motifs with residue numbers + if (defined $p_motifs{$peptide}) { + print OUT join(" /// ", @{$p_motifs{$peptide}}), "\t"; + push (@ppep_metadata, join(" /// ", @{$p_motifs{$peptide}})); + } else { + print OUT "\t"; + push (@ppep_metadata, ""); + } + + } + + # Print the PhosphoSite regulatory data + + if (defined $domain_2{$peptide}) { print OUT "$domain_2{$peptide}\t"; } else { print OUT "\t"; } + if (defined $ON_FUNCTION_2{$peptide}) { print OUT "$ON_FUNCTION_2{$peptide}\t"; } else { print OUT "\t"; } + if (defined $ON_PROCESS_2{$peptide}) { print OUT "$ON_PROCESS_2{$peptide}\t"; } else { print OUT "\t"; } + if (defined $ON_PROT_INTERACT_2{$peptide}) { print OUT "$ON_PROT_INTERACT_2{$peptide}\t"; } else { print OUT "\t"; } + if (defined $ON_OTHER_INTERACT_2{$peptide}) { print OUT "$ON_OTHER_INTERACT_2{$peptide}\t"; } else { print OUT "\t"; } + if (defined $notes_2{$peptide}) { print OUT "$notes_2{$peptide}\t"; } else { print OUT "\t"; } + + if (defined $domain_2{$peptide}) { push (@ppep_metadata, $domain_2{$peptide}); } else { push(@ppep_metadata, ""); } + if (defined $ON_FUNCTION_2{$peptide}) { push (@ppep_metadata, $ON_FUNCTION_2{$peptide}); } else { push(@ppep_metadata, ""); } + if (defined $ON_PROCESS_2{$peptide}) { push (@ppep_metadata, $ON_PROCESS_2{$peptide}); } else { push(@ppep_metadata, ""); } + if (defined $ON_PROT_INTERACT_2{$peptide}) { push (@ppep_metadata, $ON_PROT_INTERACT_2{$peptide}); } else { push(@ppep_metadata, ""); } + if (defined $ON_OTHER_INTERACT_2{$peptide}) { push (@ppep_metadata, $ON_OTHER_INTERACT_2{$peptide}); } else { push(@ppep_metadata, ""); } + if (defined $notes_2{$peptide}) { push (@ppep_metadata, $notes_2{$peptide}); } else { push(@ppep_metadata, ""); } + + # begin store-to-SQLite "ppep_metadata" table + # --- + for $i (1..14) { + $ppep_metadata_stmth->bind_param($i, $ppep_metadata[$i-1]); + } + if (not $ppep_metadata_stmth->execute()) { + print "Error writing ppep_metadata row for phosphopeptide $ppep_metadata[$i]: $ppep_metadata_stmth->errstr\n"; + } + # ... + # end store-to-SQLite "ppep_metadata" table + + # Print the data + @tmp_data = (); + foreach (@{$data{$peptide}}) { + push(@tmp_data, $_); + } + print OUT join("\t", @tmp_data), "\t"; + + # begin store-to-SQLite "ppep_intensity" table + # --- + # commit the sample intensities + $i = 0; + foreach (@{$data{$peptide}}) { + my $intense = $_; + $ppep_intensity_stmth->bind_param( 1, $ppep_id ); + $ppep_intensity_stmth->bind_param( 2, $sample_id_lut{$samples[$i]} ); + $ppep_intensity_stmth->bind_param( 3, $intense ); + if (not $ppep_intensity_stmth->execute()) { + print "Error writing tuple ($peptide,$samples[$i],$intense): $ppep_intensity_stmth->errstr\n"; + } + $i += 1; + } + # ... + # end store-to-SQLite "ppep_intensity" table + + # print the kinase-substrate data + for my $i (0 .. $#kinases_observed) { + if (exists($kinase_substrate_NetworKIN_matches{$peptide}{$kinases_observed[$i]})) { + print OUT "X\t"; + my $NetworKIN_label = $kinases_observed[$i]; #."_NetworKIN"; + print MELT "$peptide\t$gene_names\t$site_description{$SITE_KINASE_SUBSTRATE}\t$NetworKIN_label\n"; + # begin store-to-SQLite "ppep_gene_site" table + # --- + $ppep_gene_site_stmth->bind_param(1, $ppep_id); # ppep_gene_site.ppep_id + $ppep_gene_site_stmth->bind_param(2, $gene_names); # ppep_gene_site.gene_names + $ppep_gene_site_stmth->bind_param(3, $NetworKIN_label); # ppep_gene_site.kinase_map + $ppep_gene_site_stmth->bind_param(4, $SITE_KINASE_SUBSTRATE); # ppep_gene_site.site_type_id + if (not $ppep_gene_site_stmth->execute()) { + print "Error writing tuple ($peptide,$gene_names,$kinases_observed[$i]): $ppep_gene_site_stmth->errstr\n"; + } + # ... + # end store-to-SQLite "ppep_gene_site" table + } + else { print OUT "\t";} + } + my %wrote_motif; + my $motif_parts_0; + my @motif_split; + my $one_motif; + + for my $i (0 .. $#motif_type_keys) { + if (exists($kinase_motif_matches{$peptide}{$motif_type_keys[$i]})) { + print OUT "X\t"; + #ACE-2022.06.20 $motif_parts_0 = $motif_type{$motif_sequence[$i]}." ".$motif_sequence[$i]; + $motif_parts_0 = $motif_type{$motif_type_keys[$i]}; + @motif_split = split("[|]", $motif_parts_0); + #ACE-2022.06.20 my $key = "$peptide\t$gene_names\t$motif_parts_0"; + for my $j (0 .. $#motif_split) { + $one_motif = $motif_split[$j]; + #ACE-2022.06.20 my $key = "$peptide\t$gene_names\t$motif_parts_0"; + my $key = "$peptide\t$gene_names\t$one_motif"; + if (!exists($wrote_motif{$key})) { + $wrote_motif{$key} = $key; + print MELT "$peptide\t$gene_names\t$motif_description[$i]\t$one_motif\n"; + # print "Line 657: i is $i\t$kinase_motif_matches{$peptide}{$motif_sequence[$i]}\n"; #debug + # begin store-to-SQLite "ppep_gene_site" table + # --- + $ppep_gene_site_stmth->bind_param(1, $ppep_id); # ppep_gene_site.ppep_id + $ppep_gene_site_stmth->bind_param(2, $gene_names); # ppep_gene_site.gene_names + $ppep_gene_site_stmth->bind_param(3, $one_motif); # ppep_gene_site.kinase_map + $ppep_gene_site_stmth->bind_param(4, $site_id{$motif_description[$i]}); # ppep_gene_site.site_type_id + if (not $ppep_gene_site_stmth->execute()) { + print "Error writing tuple ($peptide,$gene_names,$one_motif): $ppep_gene_site_stmth->errstr\n"; + } + # ... + # end store-to-SQLite "ppep_gene_site" table + } + } + } + else { print OUT "\t";} + } + for my $i (0 .. $#kinases_PhosphoSite) { + if (exists($kinase_substrate_PhosphoSite_matches{$peptide}{$kinases_PhosphoSite[$i]})) { + print MELT "$peptide\t$gene_names\t$site_description{$SITE_PHOSPHOSITE}\t$phosphosites_observed_lbl[$i]\n"; + if ($i < $#kinases_PhosphoSite) { + print OUT "X\t"; + } + else { + print OUT "X\n"; + } + # begin store-to-SQLite "ppep_gene_site" table + # --- + $ppep_gene_site_stmth->bind_param(1, $ppep_id); # ppep_gene_site.ppep_id + $ppep_gene_site_stmth->bind_param(2, $gene_names); # ppep_gene_site.gene_names + $ppep_gene_site_stmth->bind_param(3, $phosphosites_observed_lbl[$i]); # ppep_gene_site.kinase_map + $ppep_gene_site_stmth->bind_param(4, $SITE_PHOSPHOSITE); # ppep_gene_site.site_type_id + if (not $ppep_gene_site_stmth->execute()) { + print "Error writing tuple ($peptide,$gene_names,$phosphosites_observed_lbl[$i]): $ppep_gene_site_stmth->errstr\n"; + } + # ... + # end store-to-SQLite "ppep_gene_site" table + } + else { + if ($i < $#kinases_PhosphoSite) { + print OUT "\t"; + } + elsif ($i == $#kinases_PhosphoSite) { + print OUT "\n"; + } + } + } +} + +close OUT; +close MELT; +$ppep_gene_site_stmth->finish; +print "begin DB commit at " . format_localtime_iso8601() . "\n"; +$dbh->{AutoCommit} = $auto_commit; +$dbh->disconnect if ( defined $dbh ); + +print "\nFinished writing output at " . format_localtime_iso8601() ."\n\n"; + +###############################################################################################################################