In [1]:
import sys
sys.version_info
Out[1]:
This case study illustrates the use of the petl package for doing some simple profiling and comparison of data from two tables.
The files used in this case study can be downloaded from the following link:
Download and unzip the files:
In [2]:
!wget http://aliman.s3.amazonaws.com/petl/petl-case-study-1-files.zip
!unzip -o petl-case-study-1-files.zip
The first file (snpdata.csv) contains a list of locations in the
genome of the malaria parasite P. falciparum, along with some basic
data about genetic variations found at those locations.
The second file (popdata.csv) is supposed to contain the same list
of genome locations, along with some additional data such as allele
frequencies in different populations.
The main point for this case study is that the first file
(snpdata.csv) contains the canonical list of genome locations, and
the second file (popdata.csv) contains some additional data about
the same genome locations and therefore should be consistent with the
first file. We want to check whether this second file is in fact
consistent with the first file.
Start by importing the petl package:
In [3]:
import petl as etl
etl.__version__
Out[3]:
To save some typing, let a be the table of data extracted from the
first file (snpdata.csv), and let b be the table of data extracted
from the second file (popdata.csv), using the fromcsv()
function:
In [4]:
a = etl.fromtsv('snpdata.csv')
b = etl.fromtsv('popdata.csv')
Examine the header from each file:
In [5]:
a.header()
Out[5]:
In [6]:
b.header()
Out[6]:
There is a common set of 9 fields that is present in both tables, and
we would like focus on comparing these common fields, however
different field names have been used in the two files. To simplify
comparison, use rename() to rename some fields in the
second file:
In [7]:
b_renamed = b.rename({'Chromosome': 'Chr',
'Coordinates': 'Pos',
'Ref. Allele': 'Ref',
'Non-Ref. Allele': 'Nref',
'Derived Allele': 'Der',
'Mutation type': 'Mut',
'Gene': 'GeneId',
'Gene Aliases': 'GeneAlias',
'Gene Description': 'GeneDescr'})
b_renamed.header()
Out[7]:
Use cut() to extract only the fields we're interested in
from both tables:
In [8]:
common_fields = ['Chr', 'Pos', 'Ref', 'Nref', 'Der', 'Mut', 'GeneId', 'GeneAlias', 'GeneDescr']
a_common = a.cut(common_fields)
b_common = b_renamed.cut(common_fields)
Inspect the data:
In [9]:
a_common
Out[9]:
In [10]:
b_common
Out[10]:
The fromucsv() function does not attempt to parse any of the
values from the underlying CSV file, so all values are reported as
strings:
In [11]:
b_common.display(vrepr=repr)
However, the 'Pos' field should be interpreted as an integer.
Also, the 'Mut' field has a different representation in the two tables, which needs to be converted before the data can be compared:
In [12]:
a_common.valuecounts('Mut')
Out[12]:
In [13]:
b_common.valuecounts('Mut')
Out[13]:
Use the convert() function to convert the type of the 'Pos'
field in both tables and the representation of the 'Mut' field in
table b:
In [14]:
a_conv = a_common.convert('Pos', int)
b_conv = (
b_common
.convert('Pos', int)
.convert('Mut', {'SYN': 'S', 'NON': 'N'})
)
In [15]:
highlight = 'background-color: yellow'
a_conv.display(caption='a', vrepr=repr, td_styles={'Pos': highlight})
In [16]:
b_conv.display(caption='b', vrepr=repr, td_styles={'Pos': highlight, 'Mut': highlight})
Now the tables are ready for comparison.
Because both tables should contain the same list of genome locations,
they should have the same number of rows. Use nrows() to
compare:
In [17]:
a_conv.nrows()
Out[17]:
In [18]:
b_conv.nrows()
Out[18]:
There is some discrepancy. First investigate by comparing just the
genomic locations, defined by the 'Chr' and 'Pos' fields, using
complement():
In [19]:
a_locs = a_conv.cut('Chr', 'Pos')
b_locs = b_conv.cut('Chr', 'Pos')
locs_only_in_a = a_locs.complement(b_locs)
locs_only_in_a.nrows()
Out[19]:
In [20]:
locs_only_in_a.displayall(caption='a only')
In [21]:
locs_only_in_b = b_locs.complement(a_locs)
locs_only_in_b.nrows()
Out[21]:
So it appears that 29 locations are missing from table b. Export
these missing locations to a CSV file using toucsv():
In [22]:
locs_only_in_a.tocsv('missing_locations.csv')
An alternative method for finding rows in one table where some key
value is not present in another table is to use the antijoin()
function:
In [23]:
locs_only_in_a = a_conv.antijoin(b_conv, key=('Chr', 'Pos'))
locs_only_in_a.nrows()
Out[23]:
We'd also like to compare the values given in the other fields, to find any discrepancies between the two tables.
The simplest way to find conflicts is to merge() both tables under a given key:
In [24]:
ab_merge = etl.merge(a_conv, b_conv, key=('Chr', 'Pos'))
ab_merge.display(caption='ab_merge',
td_styles=lambda v: highlight if isinstance(v, etl.Conflict) else '')
From a glance at the conflicts above, it appears there are discrepancies in the 'GeneAlias' and 'GeneDescr' fields. There may also be conflicts in other fields, so we need to investigate further.
Note that the table ab_merge will contain all rows, not only those containing conflicts. To find only conflicting rows, use cat() then conflicts(), e.g.:
In [25]:
ab = etl.cat(a_conv.addfield('source', 'a', index=0),
b_conv.addfield('source', 'b', index=0))
ab_conflicts = ab.conflicts(key=('Chr', 'Pos'), exclude='source')
ab_conflicts.display(10)
Finally, let's find conflicts in a specific field:
In [26]:
ab_conflicts_mut = ab.conflicts(key=('Chr', 'Pos'), include='Mut')
ab_conflicts_mut.display(10, caption='Mut conflicts', td_styles={'Mut': highlight})
In [27]:
ab_conflicts_mut.nrows()
Out[27]:
For more information about the petl package see the petl online documentation.