Filtering and Processing Data

Data is retrieved from snapshots from the RCSB PDB repository.

Each file is imported into a single table raw_pdb in our PostgreSQL database using a Perl script. The script uses transactions of 10,000 inserts per chunk.

The structure of raw_pdb is this:

codecharacter varying(20)not null
line_numintegernot null
line_contcharacter varying(80)not null

The import script looks like this:

#!/usr/bin/perl -w
use FindBin '$Bin';
use DBI;
$dbName = 'bioinf'; 
$dbLogin = 'ezop'; 
$dbPass = 'XYZ';
$conn = DBI->connect("DBI:Pg:database=$dbName;host=localhost", "$dbLogin", "$dbPass", {'RaiseError' => 1, 'AutoCommit' => 0});
die "./ <table> <file>" if not defined($ARGV[0]);
$recordCount = 0;
$table = $ARGV[0];
$fName = $ARGV[1];
open F, "zcat $fName|";
while (<F>) {
    $linija = $_;
    $recordCount += 1;
    $sql = "insert into $table (code, line_num, line_cont) values (?, ?, ?)";
    $conn->do($sql, undef, $fName, $recordCount, $linija);
    $conn->commit() if ($recordCount%10000 == 0);
close F;

After lines are imported, they are parsed using functions we will define below.

From raw_pdb data, we generate the tables tspsproteinssourcessources_organela, and ss_bond by parsing the corresponding records.

The ps table has three important columns: chainlength, and sequence. Protein sequence is generated using C-alpha atoms for each chain and for each residue ordered by residue sequence, taking only the first insertion and the first alternate location. chainis taken from the TS.chain column, and length is simply the precalculated length of the sequence string. Since this article is meant to analyze only single chains and intrachain connections, multiple-chain proteins are skipped in our analysis here.

Within SSBOND records, all disulfide bonds are stored in the pdb_ssbond table, which inherits from the pdb_ssbond_extendedtable. pdb_ssbond looks like this:

idintegernot nullnextval(‘pdb_ssbond_id_seq’::regclass) 
codecharacter(4)  four-letter code
ser_numinteger  serial number of ssbond
residue1character(3)  first residue in bond
chain_id1character(1)  first chain in bond
res_seq1integer  sequential number of first residue
i_code1character(1)  insertion code of first residue
residue2character(3)  second residue in bond
chain_id2character(1)  second chain in bond
res_seq2integer  sequential number of second residue
i_code2character(1)  insertion code of second residue
sym1character(6)  first symmetry operator
sym2character(6)  second symmetry operator
distnumeric(5,2)  distance between atoms
is_reactivebooleannot nullfalsemark for reactivity (to be set)
is_consecutivebooleannot nulltruemark for consecutive bonds (to be set)
rep7booleannot nullfalsemark for set-7 structures (to be set)
rep40booleannot nullfalsemark for set-40 structures (to be set)
rep80booleannot nullfalsemark for set-80 structures (to be set)
is_from_pdbboolean trueis taken from PDB as source (to be set)

I also added these indexes:

"pdb_ssbond_pkey" PRIMARY KEY, btree (id)
"ndxcode1" btree (code, chain_id1, res_seq1)
"ndxcode2" btree (code, chain_id2, res_seq2)

It is assumed that distribution of disulfide bonds prior to the cutoff is Gaussian (without testing with, e.g., KS-test), so standard deviations were calculated on each distance between cysteines in same protein to discover the range of permitted bond lengths and compare them to the cutoff. The cutoff was same as the calculated mean plus-minus three standard deviations. We have extended the range to introduce more possible disulfide bonds which weren’t enlisted in PDB file in SSBOND rows but which we have inserted ourselves by calculating distances between ATOM records. The chosen range for ssbonds are between 1.6175344456264 and 2.48801947642267 Å, which is the mean (2.05) plus-minus four standard deviations:

select                     count(1) as     cnt
    ,      stddev(dist)             as std_dev  
    ,                     avg(dist) as avg_val
    ,     -stddev(dist) + avg(dist) as   left1
    ,      stddev(dist) + avg(dist) as  right1
    , -2 * stddev(dist) + avg(dist) as   left2
    ,  2 * stddev(dist) + avg(dist) as  right2
    , -3 * stddev(dist) + avg(dist) as   left3
    ,  3 * stddev(dist) + avg(dist) as  right3
    , -4 * stddev(dist) + avg(dist) as   left4
    ,  4 * stddev(dist) + avg(dist) as  right4
    ,         min(dist)
    ,         max(dist)
from pdb_ssbond_tmp
where dist > 0
    and dist < 20;

The TS table contains the coordinates of all atoms, but only cysteines will be used, with their sulfur named " SG ". So another staging table with " SG " sulfur atoms only is created for speeding up the process by reducing the number of records to search. When selecting sulfurs only, the number of combinations is much less than in the case of all atoms—194,574 of the former compared with 122,761,100 of the latter. Within this table joined to itself, distances are calculated using the Euclidean distance, and results are imported into the pdb_ssbond table but only where the distance is between the defined lengths calculated earlier. The reason for doing this speedup is to lessen the amount of time of running the whole process again—for checking purposes—keeping it within one day.

To clean disulfide bonds, we use the following algorithm:

·         Delete when both sides of connection point to same amino acid

·         Delete bonds whose length is not between 1.6175344456264 and 2.48801947642267

·         Remove insertions

·         Remove bonds caused by alternate atom locations, but leaving first location

The code for this (taking the pdb_ssbond table as the first argument) is:

 CREATE OR REPLACE FUNCTION pdb_ssbond_clean2( 
    clean_icodes boolean,
    clean_altloc boolean,
    mark_reactive boolean,
    mark_consecutive boolean)
RETURNS void AS $$
declare _res integer;
    delete from pdb_ssbond b 
    where exists (
        from pdb_ssbond a
        where a.code = b.code
        and >
        and (
            (   a.chain_id1 = b.chain_id1 and a.res_seq1 = b.res_seq1
            and a.chain_id2 = b.chain_id2 and a.res_seq2 = b.res_seq2)
            (   a.chain_id1 = b.chain_id2 and a.res_seq1 = b.res_seq2
            and a.chain_id2 = b.chain_id1 and a.res_seq2 = b.res_seq1)
    ) ; 
    from pdb_ssbond
    where chain_id1 = chain_id2
    and res_seq1 = res_seq2
    and i_code1 = i_code2;
    update pdb_ssbond 
    set is_consecutive = true
    , is_reactive = false;
    delete from pdb_ssbond 
    where not dist between 1.6175344456264 and 2.48801947642267;
    if clean_icodes then  
        delete from pdb_ssbond 
        where  i_code1 not in (''' ', 'A') 
        or     i_code2 not in (''' ', 'A') ;
    end if;
    if clean_altloc then  
        delete from pdb_ssbond a
        where exists (
            select 1
            from pdb_ssbond  b
            where   b.code = a.code
                and b.chain_id1 = a.chain_id1
                and b.res_seq1 = a.res_seq1
                and b.i_code1 = a.i_code1
                and b.ser_num < a.ser_num
        ) ;
        delete from pdb_ssbond a
        where exists (
            select 1
            from pdb_ssbond  b
            where   b.code = a.code
                and b.chain_id2 = a.chain_id2
                and b.res_seq2 = a.res_seq2
                and b.i_code2 = a.i_code2
                and b.ser_num < a.ser_num
        ) ;
    end if;
    if mark_reactive then  
            update pdb_ssbond 
            set is_reactive = false ; 
            update pdb_ssbond 
            set is_reactive = true
            where exists (
                select 1
                from pdb_ssbond  b
                where   b.code = pdb_ssbond.code
                    and b.chain_id1 = pdb_ssbond.chain_id1
                    and b.res_seq1 = pdb_ssbond.res_seq1
                    and b.i_code1 = pdb_ssbond.i_code1
                    and b.ser_num < pdb_ssbond.ser_num
            ) ;  
            update pdb_ssbond 
            set is_reactive = true         
            where exists (
                select 1
                from pdb_ssbond  b
                where   b.code = pdb_ssbond.code
                    and b.chain_id2 = pdb_ssbond.chain_id2
                    and b.res_seq2 = pdb_ssbond.res_seq2
                    and b.i_code2 = pdb_ssbond.i_code2
                    and b.ser_num < pdb_ssbond.ser_num
            ) ; 
            update pdb_ssbond 
            set is_reactive = true         
            where (code, chain_id1, res_seq1, i_code1)
            in (
                select code, chain_id1 as c, res_seq1 as r, i_code1 as i
                from pdb_ssbond
                select code, chain_id2 as c, res_seq2 as r, i_code2 as i
                from pdb_ssbond
            ) ; 
            update pdb_ssbond 
            set is_reactive = true         
            where (code, chain_id2, res_seq2, i_code2)
            in (
                select code, chain_id1 as c, res_seq1 as r, i_code1 as i
                from pdb_ssbond
                select code, chain_id2 as c, res_seq2 as r, i_code2 as i
                from pdb_ssbond
    end if;
    if mark_consecutive then
           update pdb_ssbond
           set is_consecutive = false;
            update pdb_ssbond
            set is_consecutive = true
            where not exists (
                select 1
                from pdb_ssbond a
                    a.code = pdb_ssbond.code
                    and (
                        (a.chain_id1 = pdb_ssbond.chain_id1 and a.chain_id2 = pdb_ssbond.chain_id2)
                        (a.chain_id1 = pdb_ssbond.chain_id2 and a.chain_id2 = pdb_ssbond.chain_id1)
                    and (
                           (a.res_seq1 between pdb_ssbond.res_seq1 and pdb_ssbond.res_seq2)
                        or (a.res_seq2 between pdb_ssbond.res_seq1 and pdb_ssbond.res_seq2)
                        or (pdb_ssbond.res_seq1 between a.res_seq1 and a.res_seq2)
                        or (pdb_ssbond.res_seq2 between a.res_seq1 and a.res_seq2)
                    and not (
                            a.code = pdb_ssbond.code
                        and a.chain_id1 = pdb_ssbond.chain_id1
                        and a.chain_id2 = pdb_ssbond.chain_id2 
                        and a.res_seq1 = pdb_ssbond.res_seq1
                        and a.res_seq2 = pdb_ssbond.res_seq2
    end if; 
$$ LANGUAGE plpgsql;

With this, the non-redundant set of proteins is imported to the nrtable which is joined to the ps and proteins tables, and sets are marked (set7set40, and set80). At the end, according to protein quantity only one set will be analyzed. Mismatched chains between PDB and NR are removed from analysis.

Proteins without disulfide bonds are excluded from research, together with proteins that don’t belong to any set. Data is processed with DSSP, and these files which had problems with resolution or too many atoms to be processed are also excluded. Only proteins with single chains are used as result for analysis because interchain connections were not kept, although they are easily calculated from the ssbond table by counting the number of connections that have different chains.

For the remaining proteins, an update is done for the total number of bonds and the number of overlapping bonds, and this is done for each of the sets.

The source organism is extracted from SOURCE records. In cases where it is unknown, synthetic, designed, artificial, or hybrid, it is discarded from research. Low-resolution proteins are also excluded only when their side chain is not visible.

SOURCE records are stored in the sources table, which contains taxonomy rows. In some cases, the taxonomy is missing or incorrect. In these cases, the manual correction of experts is needed.

From the source and taxonomy downloaded from NCBI, the class is assigned to each primary structure. In case a class is unable to be assigned, the protein is removed from the analysis list. Knowing that biological databases are being used, an extra check of all source records and NCBI taxonomy classes is recommended to be performed by a biologist, otherwise there might be problems with classifications between different domains. To match source cellular locations with taxonomy IDs, data from the source table is extracted into the table sources_organela where all data is written as code, tag, and value. Its format is shown below:

select * from sources_organela where code = '1rav';

The taxonomy archive we want to use is a ZIP file containing seven dump files. Among these files, two of the most important are names.dmp and merged.dmp. Both files are CSV tab-pipe delimited files as detailed in the documentation:

·         The file merged.dmp contains a list of previous taxonomy IDs, and the current taxonomy IDs into which each one was merged.

·         names.dmp contains these important columns in this order:

o    tax_id: The taxonomy ID

o    name_txt: The name of the species, and if applicable, the unique name (where species can be found with multiple names, this helps disambiguate)

·         division.dmp contains the names of top-level domains which we will use as our classes.

·         nodes.dmp is the table which contains a hierarchical structure of organisms using taxonomy IDs.

o    It contains a parent taxonomy ID, a foreign key which can be found in names.dmp.

o    It also contains a division ID which is important for us to correctly store the relevant top domain data.

With all this data and manual corrections (setting correct domains of life) we were able to create the structure of the taxonomy_pathtable. A sampling of data looks like this:

select * from taxonomy_path order by length(pathlimit 20 offset 2000;
142182cellular organisms;Bacteria;Gemmatimonadetesfffft
136087cellular organisms;Eukaryota;Malawimonadidaeftfff
649454Viruses;unclassified phages;Cyanophage G1168tffff
321302Viruses;unclassified viruses;Tellina virus 1tffff
649453Viruses;unclassified phages;Cyanophage G1158tffff
536461Viruses;unclassified phages;Cyanophage S-SM1tffff
536462Viruses;unclassified phages;Cyanophage S-SM2tffff
77041Viruses;unclassified viruses;Stealth virus 4tffff
77042Viruses;unclassified viruses;Stealth virus 5tffff
641835Viruses;unclassified phages;Vibrio phage 512tffff
1074427Viruses;unclassified viruses;Mouse Rosavirustffff
1074428Viruses;unclassified viruses;Mouse Mosavirustffff
480920other sequences;plasmids;IncP-1 plasmid 6-S1ffftf
2441other sequences;plasmids;Plasmid ColBM-Cl139ffftf
168317other sequences;plasmids;IncQ plasmid pIE723ffftf
536472Viruses;unclassified phages;Cyanophage Syn10tffff
536474Viruses;unclassified phages;Cyanophage Syn30tffff
2407other sequences;transposons;Transposon Tn501ffftf
227307Viruses;ssDNA viruses;Circoviridae;Gyrovirustffff
687247Viruses;unclassified phages;Cyanophage ZQS-7tffff

Before any analysis, to avoid biases, sequences have to be checked for their level of identity. Although the NR set contains sequences which are already compared between each other, an extra check is always recommended.

For each disulfide bond’s prior statistical analysis, data is marked if it is reactive or overlapping. After marking overlaps, that information automatically reveals how many consecutive and non-consecutive bonds are inside each protein, and that data is stored in the proteins table from which all protein complexes are excluded in final result.

Each disulfide bond is marked also for its association to sets, by checking both bond sides to see if they belong to the same NR set. Where that is not the case, the disulfide bond is skipped for that set analysis.

To analyze the quantity of bonds by their variance, each length has to be put in a specific class. In this case, only five classes are chosen as written in the function below.

CREATE OR REPLACE FUNCTION len2class(len integer)
    RETURNS integer AS
        when len <= 100                then 1
        when len100 and len <= 200 then 2
        when len200 and len <= 300 then 3
        when len300 and len <= 400 then 4
                                       else 5
    COST 100;

Most of the protein sizes are less than 400 amino acids, so length classification is done by splitting lengths into five ranges, each taking 100 amino acids, except the last one which takes the rest. Below you can see how to use this function to extract data for analysis:

 SELECT p.code,
    ( SELECT
           FROM src_classes s
          WHERE s.src_class::text = p.src_class::text) AS src_class_id,
    ( SELECT s.val
           FROM sources_organela s
          WHERE s.code = p.code::bpchar AND s.tag::text = 'EXPRESSION_SYSTEM_CELLULAR_LOCATION'::text) AS expression_system_cellular_location,
    ( SELECT s.val
           FROM sources_organela s
          WHERE s.code = p.code::bpchar AND s.tag::text = 'CELLULAR_LOCATION'::text) AS cellular_location,
   FROM proteins p
     JOIN nr n ON n.code::text = p.code::text AND n.rep7 = true
     JOIN ps ps ON ps.code::text = n.code::text AND ps.chain_id = n.chain_id AND ps.het = false
  WHERE p.is_excluded = false AND p.chain_cnt = 1 AND p.is_set7 = true AND p.reactive_cnt7 = 0
  ORDER BY p.code;

PostgreSQL as a Processing Intermediary

In this work, we showed how to process data, from fetching to analyzing. When working with scientific data, sometimes normalization is needed, and sometimes not. When working with small quantities of data which will not be reused for another analysis, then it is enough to leave it denormalized where processing is fast enough.

The reason why this was done all in one bioinformatics database is that PostgreSQL is able to integrate many languages. This includes R, which can do statistical analysis directly in-database—the subject for a future article on bioinformatics tools.

Related Posts

© 2023 Software Engineering - Theme by WPEnjoy · Powered by WordPress