| 1 |
#!/usr/bin/perl |
| 2 |
use strict; |
| 3 |
use Getopt::Std; |
| 4 |
#use LWP::Simple; |
| 5 |
use FindBin;use lib $FindBin::Bin; |
| 6 |
use dbSession; |
| 7 |
my $usage = qq/Usage: |
| 8 |
uniprot2bcp.pl [-t <taxons.tab>] [-F] uniref100_clusters.pfa [<uniprot_data_file_or_stream>] |
| 9 |
|
| 10 |
Parses the UNIPROT data into bcp files for loading into geanno db. |
| 11 |
|
| 12 |
Also creates a uniprot fasta file (uniprot.fa) with all the |
| 13 |
uniprot entries found in the input data\/stream |
| 14 |
|
| 15 |
If no database connection is possible, use -t option to provide a file |
| 16 |
with taxon information; <taxons.tab> can be built with |
| 17 |
pullsql -b common\@CBCBMYSQL -o taxons.tab |
| 18 |
.. with the query: |
| 19 |
select tax_id, name, c_name from taxon |
| 20 |
|
| 21 |
If -F option is given, ONLY the fasta file will be created. |
| 22 |
/; |
| 23 |
umask 0002; |
| 24 |
getopts('t:F') || die($usage."\n"); |
| 25 |
my $fastaOnly=$Getopt::Std::opt_F; |
| 26 |
my $ftaxon=$Getopt::Std::opt_t; |
| 27 |
my @tax; |
| 28 |
if ($ftaxon) { |
| 29 |
open(TAX, $ftaxon)||die("Error opening $ftaxon\n"); |
| 30 |
while (<TAX>) { |
| 31 |
chomp; |
| 32 |
my @t=split(/\t/); |
| 33 |
$tax[int($t[0])]=[$t[1],$t[2]]; |
| 34 |
} |
| 35 |
close(TAX); |
| 36 |
} |
| 37 |
|
| 38 |
#create the bcp files for loading |
| 39 |
|
| 40 |
my $funiref=shift(@ARGV) || die("$usage\n"); |
| 41 |
|
| 42 |
my $ds; |
| 43 |
unless ($ftaxon) { |
| 44 |
$ds=dbSession->new('common@CBCBMYSQL'); |
| 45 |
} |
| 46 |
|
| 47 |
my %uniref; # uniprot_id => uniref_cluster_id |
| 48 |
open(UNIREF, $funiref) || die("Error opening $funiref file!\n"); |
| 49 |
while (<UNIREF>) { |
| 50 |
chomp; |
| 51 |
s/^>//; |
| 52 |
my @ids=split; |
| 53 |
unless (@ids>1) { |
| 54 |
print STDERR "WARNING: cluster: '$_' has no members!\n"; |
| 55 |
next; |
| 56 |
} |
| 57 |
my $clname=shift(@ids); |
| 58 |
foreach my $upid (@ids) { |
| 59 |
$uniref{$upid}=$clname; |
| 60 |
} |
| 61 |
} |
| 62 |
close(UNIREF); |
| 63 |
print STDERR "UniRef clusters loaded, now parsing the main uniprot KB data..\n"; |
| 64 |
|
| 65 |
open(UPFA, '>uniprot.fa') || die ("Error creating file uniprot.fa!\n"); |
| 66 |
unless ($fastaOnly) { |
| 67 |
open(UNIPROT, '>uniprot.bcp') || die ("Error creating file uniprot.bcp!\n"); |
| 68 |
open(UNAMES, '>uniprot_names.bcp') || die ("Error creating file uniprot_names.bcp!\n"); |
| 69 |
open(UXREF, '>uniprot_xref.bcp') || die ("Error creating file uniprot_xref.bcp!\n"); |
| 70 |
} |
| 71 |
|
| 72 |
#open(FGO, $file) || die ("Error opening $file\n"); |
| 73 |
my $sth; |
| 74 |
unless ($ftaxon) { |
| 75 |
$sth=$ds->prep('select name, c_name from taxon where tax_id=?'); |
| 76 |
} |
| 77 |
my $c_id; #current uniprot ID |
| 78 |
my $c_taxon; |
| 79 |
my $c_rev=0; |
| 80 |
my $c_len; |
| 81 |
my $c_priacc; |
| 82 |
my $c_descr; |
| 83 |
my $c_gene; |
| 84 |
my $c_uniref; |
| 85 |
my $c_htaxon; |
| 86 |
my $seq; |
| 87 |
my $inSeq=0; |
| 88 |
my $xcode=0; |
| 89 |
my ($c_orgname, $c_comname); |
| 90 |
while (<>) { |
| 91 |
if (m{^//}) { |
| 92 |
writeUP() if $c_id; |
| 93 |
next; |
| 94 |
} |
| 95 |
if ($inSeq) { |
| 96 |
#within SQ record |
| 97 |
chomp; |
| 98 |
tr/\t //d; |
| 99 |
$seq.=$_; |
| 100 |
next; |
| 101 |
} |
| 102 |
if (m/^SQ\s+SEQUENCE\s+(\d+)/) { |
| 103 |
my $seqlen=$1; |
| 104 |
die("Error: invalid length for $c_id ($c_len vs $seqlen)!\n") if ($seqlen != $c_len); |
| 105 |
$inSeq=1; |
| 106 |
next; |
| 107 |
} |
| 108 |
if (m/^ID\s+(\S+)\s+(\S+)\s+(\d+)/) { |
| 109 |
my ($id, $rev, $len)=($1,$2,$3); |
| 110 |
writeUP() if $c_id; |
| 111 |
($c_id, $c_len)=($id, $len); |
| 112 |
$c_rev=($rev=~/^Reviewed/)?'1':'0'; |
| 113 |
next; |
| 114 |
}# ID |
| 115 |
if (m/^AC\s+(\S.+)/) { |
| 116 |
my $a=$1; |
| 117 |
chomp; |
| 118 |
my @accs=split(/;\s*/, $a); |
| 119 |
$c_priacc=$accs[0] unless $c_priacc; |
| 120 |
unless ($fastaOnly) { |
| 121 |
foreach my $acc (@accs) { |
| 122 |
next unless $acc; |
| 123 |
print UXREF join("\t",$c_id, 'UP_ACC', '', $acc, '')."\n"; |
| 124 |
} |
| 125 |
} |
| 126 |
next; |
| 127 |
}# AC |
| 128 |
if (m/^DE\s+(\S.+)/) { |
| 129 |
my $d=$1; |
| 130 |
$c_descr.=' '.$d; |
| 131 |
next; |
| 132 |
} |
| 133 |
if (m/^PE\s+(\d+)\:/) { |
| 134 |
$xcode=$1; |
| 135 |
next; |
| 136 |
} |
| 137 |
if (m/^GN\s+(\S.+)/) { |
| 138 |
my $d=$1; |
| 139 |
my ($gname)=($d=~m/\bName\s*=\s*([\w\-\.]+)/); |
| 140 |
if ($gname) { |
| 141 |
$c_gene=$gname unless $c_gene; |
| 142 |
unless ($fastaOnly) { |
| 143 |
print UXREF join("\t", $c_id, 'UP_GENE', '', $gname, '')."\n"; |
| 144 |
my ($gsyns)=($d=~m/\bSynonyms\s*=\s*(\S[\w\-\,\.\s]+)\;/); |
| 145 |
if ($gsyns) { |
| 146 |
my @gsyn=split(/\,\s*/,$gsyns); |
| 147 |
foreach my $gs (@gsyn) { |
| 148 |
print UNAMES join("\t",$gname, 'GSYN', $gs)."\n"; |
| 149 |
} |
| 150 |
} #synonyms found |
| 151 |
} |
| 152 |
}#a gname found |
| 153 |
next; |
| 154 |
}#GN |
| 155 |
if (s/^OS\s+//) { |
| 156 |
($c_orgname)=(m/^([\w\. \-]+)/); |
| 157 |
$c_orgname=~s/\s+$//; |
| 158 |
$c_orgname=~s/\.$//; |
| 159 |
($c_comname)=(m/\(([\w\. \-]+)/); |
| 160 |
next; |
| 161 |
} |
| 162 |
if (m/^OX\s+NCBI_TaxID\s*=\s*(\d+)/) { |
| 163 |
my $t=$1; |
| 164 |
$c_taxon=$t unless $c_taxon; |
| 165 |
print UXREF join("\t", $c_id, 'UP_Taxon', $t, '', '')."\n" |
| 166 |
unless $fastaOnly; #taxon |
| 167 |
next; |
| 168 |
}#OX (Taxon) |
| 169 |
|
| 170 |
if (m/^OH\s+NCBI_TaxID\s*=\s*(\d+)/) { |
| 171 |
my $t=$1; |
| 172 |
$c_htaxon=$t unless $c_htaxon; |
| 173 |
print UXREF join("\t", $c_id, 'UP_Host', $t, '', '')."\n" |
| 174 |
unless $fastaOnly; #host taxon |
| 175 |
next; |
| 176 |
} #OH (Host taxon) |
| 177 |
|
| 178 |
if (m/^DR\s+(\S.+)/) { # db xref |
| 179 |
my $d=$1; |
| 180 |
$d=~s/\-?\.\s*$//; |
| 181 |
$d=~s/\;\s*$//; |
| 182 |
my @xd=split(/\;\s*/,$d); |
| 183 |
@xd=grep { !/^\-/ } @xd; |
| 184 |
my $xdb=shift(@xd); |
| 185 |
if ($xdb eq 'GO') { |
| 186 |
my $goid=shift(@xd); |
| 187 |
my ($gonum)=($goid=~/GO:(\d+)/); |
| 188 |
$gonum=int($gonum); |
| 189 |
die("Error parsing GO ID from DR line: $_\n") unless $gonum>0; |
| 190 |
my $gobranch=shift(@xd); |
| 191 |
my $go_ann_type=shift(@xd); |
| 192 |
print UXREF join("\t", $c_id, $xdb, $gonum, $goid, $go_ann_type.' '.$gobranch)."\n" |
| 193 |
unless $fastaOnly; |
| 194 |
} |
| 195 |
else { |
| 196 |
my $acc=shift(@xd); |
| 197 |
my $data=join(' ', @xd); |
| 198 |
print UXREF join("\t", $c_id, $xdb, '', $acc, $data)."\n" unless $fastaOnly; |
| 199 |
} |
| 200 |
next; |
| 201 |
} #DR |
| 202 |
} #while |
| 203 |
|
| 204 |
writeUP() if $c_id; |
| 205 |
unless ($fastaOnly) { |
| 206 |
close(UNIPROT); |
| 207 |
close(UNAMES); |
| 208 |
close(UXREF); |
| 209 |
} |
| 210 |
close(UPFA); |
| 211 |
unless ($ftaxon) { |
| 212 |
$ds->logout(); |
| 213 |
} |
| 214 |
#-- now run bcpin |
| 215 |
|
| 216 |
#my $bcpcmd='bcpin -TI -b geanno@NEOSYBASE'; |
| 217 |
#system($bcpcmd. ' uniprot.bcp uniprot_names.bcp uniprot_xref') |
| 218 |
# && die("Error at bcpin!\n"); |
| 219 |
|
| 220 |
sub writeUP { |
| 221 |
# c_descr now can contain strange RecName: Full=...; and AltName: Full=...; |
| 222 |
# keep only the RecName Full, remove anything else |
| 223 |
if ($c_descr=~m/Name:\s+Full=\s*([^\;]+)/) { # RecName or SubName |
| 224 |
$c_descr=$1; |
| 225 |
} |
| 226 |
# elsif ($c_descr=~m/SubName:\s+Full=\s*([^\;]+)/) { |
| 227 |
# $c_descr=$1; |
| 228 |
# } |
| 229 |
|
| 230 |
prepFields($c_descr); |
| 231 |
# -- c_descr can probably be insanely long.. truncate it here to the max field length |
| 232 |
$c_uniref=$uniref{$c_id} || $uniref{$c_priacc}; |
| 233 |
unless ($fastaOnly) { |
| 234 |
print UNIPROT join("\t", |
| 235 |
$c_id, $c_taxon, $c_rev, $c_len, $c_priacc, $c_descr, $c_gene, $xcode, $c_uniref, $c_htaxon)."\n"; |
| 236 |
} |
| 237 |
my $upre=($c_rev==1)?'UPr|':'UP|'; |
| 238 |
#$upre.='e' if $xcode<3; |
| 239 |
my $defline=''; |
| 240 |
$defline.='|gid|'.$c_gene if $c_gene; |
| 241 |
$defline.=' '.$c_uniref if $c_uniref; |
| 242 |
$defline.=' '.$c_descr; |
| 243 |
my ($oname, $coname); |
| 244 |
if ($ftaxon) { |
| 245 |
my $td=$tax[$c_taxon]; |
| 246 |
if ($td) { |
| 247 |
($oname, $coname)=@$td; |
| 248 |
} |
| 249 |
#else { |
| 250 |
# print STDERR "Warning: no data for taxon $c_taxon". |
| 251 |
# " was loaded ($c_id|acc|$c_priacc $defline)!\n"; |
| 252 |
# } |
| 253 |
} |
| 254 |
else { |
| 255 |
$ds->exec($sth, $c_taxon); |
| 256 |
while (my $rr=$ds->fetch($sth)) { ($oname, $coname)=@$rr; }; |
| 257 |
} |
| 258 |
unless ($oname) { |
| 259 |
print STDERR "Warning: no data for taxon $c_taxon". |
| 260 |
" ($c_id|acc|$c_priacc $defline) - using '$c_orgname' instead.\n"; |
| 261 |
$oname=$c_orgname; |
| 262 |
} |
| 263 |
$defline.=' {'.$oname.'}'; |
| 264 |
$defline=~tr/;"\t/.' /; #" |
| 265 |
print UPFA '>'.$upre.$c_id.'|acc|'.$c_priacc.$defline."\n"; |
| 266 |
print UPFA fastafmt(\$seq); |
| 267 |
$inSeq=0; |
| 268 |
($c_id, $c_taxon, $c_rev, $c_len, $c_priacc, $c_descr, $c_gene, $xcode, $c_uniref, $c_htaxon)= |
| 269 |
(undef, undef, 0, undef, undef, undef, undef, 0, undef, undef ); |
| 270 |
($c_orgname, $c_comname)=(undef,undef); |
| 271 |
$seq=''; |
| 272 |
} |
| 273 |
|
| 274 |
sub fastafmt { |
| 275 |
my $s=shift; |
| 276 |
my $slen=length($$s); |
| 277 |
my @lines=unpack("A60" x (int(($slen-1)/60)+1),$$s); |
| 278 |
return join("\n",@lines)."\n"; |
| 279 |
} |
| 280 |
|
| 281 |
sub prepFields { |
| 282 |
foreach (@_) { |
| 283 |
s/^\s+//; |
| 284 |
s/\s+$//; |
| 285 |
s/\.$//; |
| 286 |
tr/\t \n/ /s; |
| 287 |
} |
| 288 |
} |
