Analiza ekspresije gena - Radna verzija

In [125]:
import GEOparse
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

from scipy.stats import ttest_ind
import stringdb
In [4]:
gse = GEOparse.get_GEO(geo="GSE33267", destdir="./")
04-Mar-2021 09:28:43 DEBUG utils - Directory ./ already exists. Skipping.
04-Mar-2021 09:28:43 INFO GEOparse - File already exist: using local version.
04-Mar-2021 09:28:43 INFO GEOparse - Parsing ./GSE33267_family.soft.gz: 
04-Mar-2021 09:28:43 DEBUG GEOparse - DATABASE: GeoMiame
04-Mar-2021 09:28:43 DEBUG GEOparse - SERIES: GSE33267
04-Mar-2021 09:28:43 DEBUG GEOparse - PLATFORM: GPL4133
04-Mar-2021 09:28:44 WARNING GEOTypes - Detected duplicated columns in d GPL4133. Correcting.

04-Mar-2021 09:28:44 DEBUG GEOparse - SAMPLE: GSM823177
04-Mar-2021 09:28:44 DEBUG GEOparse - SAMPLE: GSM823178
04-Mar-2021 09:28:44 DEBUG GEOparse - SAMPLE: GSM823179
04-Mar-2021 09:28:44 DEBUG GEOparse - SAMPLE: GSM823180
04-Mar-2021 09:28:45 DEBUG GEOparse - SAMPLE: GSM823181
04-Mar-2021 09:28:45 DEBUG GEOparse - SAMPLE: GSM823182
04-Mar-2021 09:28:45 DEBUG GEOparse - SAMPLE: GSM823183
04-Mar-2021 09:28:45 DEBUG GEOparse - SAMPLE: GSM823184
04-Mar-2021 09:28:45 DEBUG GEOparse - SAMPLE: GSM823185
04-Mar-2021 09:28:45 DEBUG GEOparse - SAMPLE: GSM823186
04-Mar-2021 09:28:45 DEBUG GEOparse - SAMPLE: GSM823187
04-Mar-2021 09:28:45 DEBUG GEOparse - SAMPLE: GSM823188
04-Mar-2021 09:28:45 DEBUG GEOparse - SAMPLE: GSM823189
04-Mar-2021 09:28:45 DEBUG GEOparse - SAMPLE: GSM823190
04-Mar-2021 09:28:46 DEBUG GEOparse - SAMPLE: GSM823191
04-Mar-2021 09:28:46 DEBUG GEOparse - SAMPLE: GSM823192
04-Mar-2021 09:28:46 DEBUG GEOparse - SAMPLE: GSM823193
04-Mar-2021 09:28:46 DEBUG GEOparse - SAMPLE: GSM823194
04-Mar-2021 09:28:46 DEBUG GEOparse - SAMPLE: GSM823195
04-Mar-2021 09:28:46 DEBUG GEOparse - SAMPLE: GSM823196
04-Mar-2021 09:28:46 DEBUG GEOparse - SAMPLE: GSM823197
04-Mar-2021 09:28:46 DEBUG GEOparse - SAMPLE: GSM823198
04-Mar-2021 09:28:46 DEBUG GEOparse - SAMPLE: GSM823199
04-Mar-2021 09:28:47 DEBUG GEOparse - SAMPLE: GSM823200
04-Mar-2021 09:28:47 DEBUG GEOparse - SAMPLE: GSM823201
04-Mar-2021 09:28:47 DEBUG GEOparse - SAMPLE: GSM823202
04-Mar-2021 09:28:47 DEBUG GEOparse - SAMPLE: GSM823203
04-Mar-2021 09:28:47 DEBUG GEOparse - SAMPLE: GSM823204
04-Mar-2021 09:28:47 DEBUG GEOparse - SAMPLE: GSM823205
04-Mar-2021 09:28:47 DEBUG GEOparse - SAMPLE: GSM823206
04-Mar-2021 09:28:47 DEBUG GEOparse - SAMPLE: GSM823207
04-Mar-2021 09:28:47 DEBUG GEOparse - SAMPLE: GSM823208
04-Mar-2021 09:28:47 DEBUG GEOparse - SAMPLE: GSM823209
04-Mar-2021 09:28:48 DEBUG GEOparse - SAMPLE: GSM823210
04-Mar-2021 09:28:48 DEBUG GEOparse - SAMPLE: GSM823211
04-Mar-2021 09:28:48 DEBUG GEOparse - SAMPLE: GSM823212
04-Mar-2021 09:28:48 DEBUG GEOparse - SAMPLE: GSM823213
04-Mar-2021 09:28:48 DEBUG GEOparse - SAMPLE: GSM823214
04-Mar-2021 09:28:48 DEBUG GEOparse - SAMPLE: GSM823215
04-Mar-2021 09:28:48 DEBUG GEOparse - SAMPLE: GSM823216
04-Mar-2021 09:28:48 DEBUG GEOparse - SAMPLE: GSM823217
04-Mar-2021 09:28:48 DEBUG GEOparse - SAMPLE: GSM823218
04-Mar-2021 09:28:49 DEBUG GEOparse - SAMPLE: GSM823219
04-Mar-2021 09:28:49 DEBUG GEOparse - SAMPLE: GSM823220
04-Mar-2021 09:28:49 DEBUG GEOparse - SAMPLE: GSM823221
04-Mar-2021 09:28:49 DEBUG GEOparse - SAMPLE: GSM823222
04-Mar-2021 09:28:49 DEBUG GEOparse - SAMPLE: GSM823223
04-Mar-2021 09:28:49 DEBUG GEOparse - SAMPLE: GSM823224
04-Mar-2021 09:28:49 DEBUG GEOparse - SAMPLE: GSM823225
04-Mar-2021 09:28:49 DEBUG GEOparse - SAMPLE: GSM823226
04-Mar-2021 09:28:49 DEBUG GEOparse - SAMPLE: GSM823227
04-Mar-2021 09:28:50 DEBUG GEOparse - SAMPLE: GSM823228
04-Mar-2021 09:28:50 DEBUG GEOparse - SAMPLE: GSM823229
04-Mar-2021 09:28:50 DEBUG GEOparse - SAMPLE: GSM823230
04-Mar-2021 09:28:50 DEBUG GEOparse - SAMPLE: GSM823231
04-Mar-2021 09:28:50 DEBUG GEOparse - SAMPLE: GSM823232
04-Mar-2021 09:28:50 DEBUG GEOparse - SAMPLE: GSM823233
04-Mar-2021 09:28:50 DEBUG GEOparse - SAMPLE: GSM823234
04-Mar-2021 09:28:50 DEBUG GEOparse - SAMPLE: GSM823235
04-Mar-2021 09:28:51 DEBUG GEOparse - SAMPLE: GSM823236
04-Mar-2021 09:28:51 DEBUG GEOparse - SAMPLE: GSM823237
04-Mar-2021 09:28:51 DEBUG GEOparse - SAMPLE: GSM823238
04-Mar-2021 09:28:51 DEBUG GEOparse - SAMPLE: GSM823239
04-Mar-2021 09:28:51 DEBUG GEOparse - SAMPLE: GSM823240
04-Mar-2021 09:28:51 DEBUG GEOparse - SAMPLE: GSM823241
04-Mar-2021 09:28:51 DEBUG GEOparse - SAMPLE: GSM823242
04-Mar-2021 09:28:51 DEBUG GEOparse - SAMPLE: GSM823243
04-Mar-2021 09:28:51 DEBUG GEOparse - SAMPLE: GSM823244
04-Mar-2021 09:28:52 DEBUG GEOparse - SAMPLE: GSM823245
04-Mar-2021 09:28:52 DEBUG GEOparse - SAMPLE: GSM823246
04-Mar-2021 09:28:52 DEBUG GEOparse - SAMPLE: GSM823247
04-Mar-2021 09:28:52 DEBUG GEOparse - SAMPLE: GSM823248
04-Mar-2021 09:28:52 DEBUG GEOparse - SAMPLE: GSM823249
04-Mar-2021 09:28:52 DEBUG GEOparse - SAMPLE: GSM823250
04-Mar-2021 09:28:52 DEBUG GEOparse - SAMPLE: GSM823251
04-Mar-2021 09:28:52 DEBUG GEOparse - SAMPLE: GSM823252
04-Mar-2021 09:28:52 DEBUG GEOparse - SAMPLE: GSM823253
04-Mar-2021 09:28:53 DEBUG GEOparse - SAMPLE: GSM823254
04-Mar-2021 09:28:53 DEBUG GEOparse - SAMPLE: GSM823255
04-Mar-2021 09:28:53 DEBUG GEOparse - SAMPLE: GSM823256
04-Mar-2021 09:28:53 DEBUG GEOparse - SAMPLE: GSM823257
04-Mar-2021 09:28:53 DEBUG GEOparse - SAMPLE: GSM823258
04-Mar-2021 09:28:53 DEBUG GEOparse - SAMPLE: GSM823259
04-Mar-2021 09:28:53 DEBUG GEOparse - SAMPLE: GSM823260
04-Mar-2021 09:28:53 DEBUG GEOparse - SAMPLE: GSM823261
04-Mar-2021 09:28:53 DEBUG GEOparse - SAMPLE: GSM823262
04-Mar-2021 09:28:54 DEBUG GEOparse - SAMPLE: GSM823263
04-Mar-2021 09:28:54 DEBUG GEOparse - SAMPLE: GSM823264
04-Mar-2021 09:28:54 DEBUG GEOparse - SAMPLE: GSM823265
04-Mar-2021 09:28:54 DEBUG GEOparse - SAMPLE: GSM823266
04-Mar-2021 09:28:54 DEBUG GEOparse - SAMPLE: GSM823267
04-Mar-2021 09:28:54 DEBUG GEOparse - SAMPLE: GSM823268
04-Mar-2021 09:28:54 DEBUG GEOparse - SAMPLE: GSM823269
04-Mar-2021 09:28:54 DEBUG GEOparse - SAMPLE: GSM823270
04-Mar-2021 09:28:54 DEBUG GEOparse - SAMPLE: GSM823271
04-Mar-2021 09:28:55 DEBUG GEOparse - SAMPLE: GSM823272
04-Mar-2021 09:28:55 DEBUG GEOparse - SAMPLE: GSM823273
04-Mar-2021 09:28:55 DEBUG GEOparse - SAMPLE: GSM823274
04-Mar-2021 09:28:55 DEBUG GEOparse - SAMPLE: GSM823275
In [7]:
dir(gse)
Out[7]:
['__class__',
 '__delattr__',
 '__dict__',
 '__dir__',
 '__doc__',
 '__eq__',
 '__format__',
 '__ge__',
 '__getattribute__',
 '__gt__',
 '__hash__',
 '__init__',
 '__init_subclass__',
 '__le__',
 '__lt__',
 '__metaclass__',
 '__module__',
 '__ne__',
 '__new__',
 '__reduce__',
 '__reduce_ex__',
 '__repr__',
 '__setattr__',
 '__sizeof__',
 '__str__',
 '__subclasshook__',
 '__weakref__',
 '_get_metadata_as_string',
 '_get_object_as_soft',
 '_phenotype_data',
 'database',
 'download_SRA',
 'download_supplementary_files',
 'geotype',
 'get_accession',
 'get_metadata_attribute',
 'get_type',
 'gpls',
 'gsms',
 'merge_and_average',
 'metadata',
 'name',
 'phenotype_data',
 'pivot_and_annotate',
 'pivot_samples',
 'relations',
 'show_metadata',
 'to_soft']
In [8]:
gse.show_metadata()
!Series_title = SCL005: icSARS CoV Urbani or icSARS deltaORF6 infections of the 2B4 clonal derivative of Calu-3 cells - Time course
!Series_geo_accession = GSE33267
!Series_status = Public on Nov 01 2011
!Series_submission_date = Oct 26 2011
!Series_last_update_date = Feb 22 2018
!Series_pubmed_id = 23365422
!Series_summary = Purpose of experiment was to compare transcriptomics of 2B4 cells (clonal derivative of Calu-3 cells) infected with either icSARS CoV or the icSARS deltaORF6 mutant.
!Series_overall_design = Calu-3 cells were infected with either icSARS CoV or the icSARS deltaORF6 mutant at MOI of 5.0.  Cells samples were collected at 0, 3, 7, 12, 24, 30, 36, 48, 54, 60 or 72h post infection. Each infected sample was done in triplicate. (Triplicates are defined as 3 different wells, plated at the same time using the same cell stock for all replicates.)There are triplicate time-matched mock for each time point from the same cell stock as rest of samples. Culture medium (the same as what the virus stock is in) was used for the mock infections.
!Series_overall_design = The NIAID Systems Virology Center
!Series_type = Expression profiling by array
!Series_contributor = Armand,,Bankhead
!Series_contributor = Jean,,Chang
!Series_contributor = Amy,,Sims
!Series_contributor = Michael,,Katze
!Series_contributor = Ralph,,Baric
!Series_contributor = Shannon,,McWeeney
!Series_sample_id = GSM823177
!Series_sample_id = GSM823178
!Series_sample_id = GSM823179
!Series_sample_id = GSM823180
!Series_sample_id = GSM823181
!Series_sample_id = GSM823182
!Series_sample_id = GSM823183
!Series_sample_id = GSM823184
!Series_sample_id = GSM823185
!Series_sample_id = GSM823186
!Series_sample_id = GSM823187
!Series_sample_id = GSM823188
!Series_sample_id = GSM823189
!Series_sample_id = GSM823190
!Series_sample_id = GSM823191
!Series_sample_id = GSM823192
!Series_sample_id = GSM823193
!Series_sample_id = GSM823194
!Series_sample_id = GSM823195
!Series_sample_id = GSM823196
!Series_sample_id = GSM823197
!Series_sample_id = GSM823198
!Series_sample_id = GSM823199
!Series_sample_id = GSM823200
!Series_sample_id = GSM823201
!Series_sample_id = GSM823202
!Series_sample_id = GSM823203
!Series_sample_id = GSM823204
!Series_sample_id = GSM823205
!Series_sample_id = GSM823206
!Series_sample_id = GSM823207
!Series_sample_id = GSM823208
!Series_sample_id = GSM823209
!Series_sample_id = GSM823210
!Series_sample_id = GSM823211
!Series_sample_id = GSM823212
!Series_sample_id = GSM823213
!Series_sample_id = GSM823214
!Series_sample_id = GSM823215
!Series_sample_id = GSM823216
!Series_sample_id = GSM823217
!Series_sample_id = GSM823218
!Series_sample_id = GSM823219
!Series_sample_id = GSM823220
!Series_sample_id = GSM823221
!Series_sample_id = GSM823222
!Series_sample_id = GSM823223
!Series_sample_id = GSM823224
!Series_sample_id = GSM823225
!Series_sample_id = GSM823226
!Series_sample_id = GSM823227
!Series_sample_id = GSM823228
!Series_sample_id = GSM823229
!Series_sample_id = GSM823230
!Series_sample_id = GSM823231
!Series_sample_id = GSM823232
!Series_sample_id = GSM823233
!Series_sample_id = GSM823234
!Series_sample_id = GSM823235
!Series_sample_id = GSM823236
!Series_sample_id = GSM823237
!Series_sample_id = GSM823238
!Series_sample_id = GSM823239
!Series_sample_id = GSM823240
!Series_sample_id = GSM823241
!Series_sample_id = GSM823242
!Series_sample_id = GSM823243
!Series_sample_id = GSM823244
!Series_sample_id = GSM823245
!Series_sample_id = GSM823246
!Series_sample_id = GSM823247
!Series_sample_id = GSM823248
!Series_sample_id = GSM823249
!Series_sample_id = GSM823250
!Series_sample_id = GSM823251
!Series_sample_id = GSM823252
!Series_sample_id = GSM823253
!Series_sample_id = GSM823254
!Series_sample_id = GSM823255
!Series_sample_id = GSM823256
!Series_sample_id = GSM823257
!Series_sample_id = GSM823258
!Series_sample_id = GSM823259
!Series_sample_id = GSM823260
!Series_sample_id = GSM823261
!Series_sample_id = GSM823262
!Series_sample_id = GSM823263
!Series_sample_id = GSM823264
!Series_sample_id = GSM823265
!Series_sample_id = GSM823266
!Series_sample_id = GSM823267
!Series_sample_id = GSM823268
!Series_sample_id = GSM823269
!Series_sample_id = GSM823270
!Series_sample_id = GSM823271
!Series_sample_id = GSM823272
!Series_sample_id = GSM823273
!Series_sample_id = GSM823274
!Series_sample_id = GSM823275
!Series_contact_name = Michael,,Katze
!Series_contact_email = data@viromics.washington.edu
!Series_contact_laboratory = Michael G. Katze, Ph.D
!Series_contact_department = Microbiology
!Series_contact_institute = University of Washington
!Series_contact_address = Rosen Building 960 Republican St.
!Series_contact_city = Seattle
!Series_contact_state = WA
!Series_contact_zip/postal_code = 98109-4325
!Series_contact_country = USA
!Series_supplementary_file = ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE33nnn/GSE33267/suppl/GSE33267_RAW.tar
!Series_supplementary_file = ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE33nnn/GSE33267/suppl/GSE33267_marray_matrix_readme.txt
!Series_supplementary_file = ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE33nnn/GSE33267/suppl/GSE33267_scl005_EXPRS.txt.gz
!Series_platform_id = GPL4133
!Series_platform_taxid = 9606
!Series_sample_taxid = 9606
!Series_relation = BioProject: https://www.ncbi.nlm.nih.gov/bioproject/PRJNA149059
In [12]:
dir(gse.gpls['GPL4133'])
Out[12]:
['__class__',
 '__delattr__',
 '__dict__',
 '__dir__',
 '__doc__',
 '__eq__',
 '__format__',
 '__ge__',
 '__getattribute__',
 '__gt__',
 '__hash__',
 '__init__',
 '__init_subclass__',
 '__le__',
 '__lt__',
 '__metaclass__',
 '__module__',
 '__ne__',
 '__new__',
 '__reduce__',
 '__reduce_ex__',
 '__repr__',
 '__setattr__',
 '__sizeof__',
 '__str__',
 '__subclasshook__',
 '__weakref__',
 '_get_columns_as_string',
 '_get_metadata_as_string',
 '_get_object_as_soft',
 '_get_table_as_string',
 'columns',
 'database',
 'geotype',
 'get_accession',
 'get_metadata_attribute',
 'get_type',
 'gses',
 'gsms',
 'head',
 'metadata',
 'name',
 'relations',
 'show_columns',
 'show_metadata',
 'show_table',
 'table',
 'to_soft']
In [91]:
gse.gpls['GPL4133'].table
gpls = gse.gpls['GPL4133'].table
In [92]:
gpls
Out[92]:
ID COL ROW NAME SPOT_ID CONTROL_TYPE REFSEQ GB_ACC GENE GENE_SYMBOL ... ENSEMBL_ID TIGR_ID ACCESSION_STRING CHROMOSOMAL_LOCATION CYTOBAND DESCRIPTION GO_ID SEQUENCE SPOT_ID.1 ORDER
0 1 266 170 GE_BrightCorner GE_BrightCorner pos NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN 1
1 2 266 168 DarkCorner DarkCorner pos NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN 2
2 3 266 166 DarkCorner DarkCorner pos NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN 3
3 4 266 164 DarkCorner DarkCorner pos NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN 4
4 5 266 162 DarkCorner DarkCorner pos NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN 5
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
45215 45216 1 9 DarkCorner DarkCorner pos NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN 45216
45216 45217 1 7 DarkCorner DarkCorner pos NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN 45217
45217 45218 1 5 DarkCorner DarkCorner pos NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN 45218
45218 45219 1 3 GE_BrightCorner GE_BrightCorner pos NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN 45219
45219 45220 1 1 GE_BrightCorner GE_BrightCorner pos NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN 45220

45220 rows × 22 columns

In [16]:
gse.gsms['GSM823177'].table
Out[16]:
ID_REF VALUE
0 1 104000.00
1 2 3.61
2 3 3.62
3 4 3.63
4 5 3.63
... ... ...
45010 45216 3.29
45011 45217 3.29
45012 45218 3.55
45013 45219 111000.00
45014 45220 108000.00

45015 rows × 2 columns

In [32]:
def process_gsms(gsms):
    df = None
    initialized = False
    for k in gsms:
        curr_df = gsms[k].table
        if not initialized:
            df = pd.DataFrame(data=curr_df.values, columns=['ID_REF', k])
            initialized = True
        else:
            new_df = pd.DataFrame(data=curr_df.loc[:,['VALUE']].values, columns=[k])
            df = df.join(new_df)
    return df
In [48]:
expression_data = process_gsms(gse.gsms)
expression_data.iloc[:, 1:] = np.log2(expression_data.iloc[:, 1:])
In [49]:
expression_data
Out[49]:
ID_REF GSM823177 GSM823178 GSM823179 GSM823180 GSM823181 GSM823182 GSM823183 GSM823184 GSM823185 ... GSM823266 GSM823267 GSM823268 GSM823269 GSM823270 GSM823271 GSM823272 GSM823273 GSM823274 GSM823275
0 1.0 16.666224 16.919981 16.376872 16.133704 16.459239 16.508042 16.595141 15.476746 16.356383 ... 16.236313 15.649481 13.813781 16.129685 16.269565 16.487777 16.333854 16.285908 15.467223 16.298492
1 2.0 1.851999 1.575312 1.516015 1.163499 2.367371 1.042644 0.903038 0.565597 0.704872 ... 0.941106 0.895303 0.765535 1.438293 0.847997 0.925999 1.169925 1.250962 1.014355 1.063503
2 3.0 1.855990 1.575312 0.526069 1.157044 1.042644 1.042644 0.903038 0.565597 0.704872 ... 0.941106 0.895303 1.327687 1.443607 0.847997 0.925999 1.169925 1.250962 0.722466 1.063503
3 4.0 1.859970 1.575312 0.526069 1.157044 1.042644 1.042644 0.910733 0.565597 0.704872 ... 0.941106 0.903038 1.090853 1.443607 0.847997 0.925999 1.169925 1.250962 0.722466 1.063503
4 5.0 1.859970 1.570463 0.526069 1.150560 1.042644 1.042644 0.910733 0.565597 0.704872 ... 0.941106 1.744161 0.748461 1.443607 0.847997 0.933573 1.169925 1.250962 0.722466 1.063503
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
45010 45216.0 1.718088 1.536053 3.058316 1.124328 1.049631 1.130931 0.839960 0.782409 0.545968 ... 0.910733 0.799087 0.650765 1.367371 0.790772 0.925999 1.189034 1.077243 0.739848 1.049631
45011 45217.0 1.718088 1.536053 3.237258 1.124328 5.133399 1.130931 0.831877 1.333424 0.545968 ... 0.918386 0.799087 0.650765 1.367371 0.790772 0.925999 1.182692 1.077243 1.220330 1.049631
45012 45218.0 1.827819 1.541019 2.304511 1.130931 1.042644 1.130931 0.823749 1.416840 1.207893 ... 0.918386 0.799087 0.659925 1.367371 0.790772 0.918386 1.182692 1.077243 0.731183 1.049631
45013 45219.0 16.760200 16.919981 16.267738 16.376872 16.289515 16.186888 16.500282 16.078484 16.451211 ... 16.153631 15.556746 16.232571 16.161526 16.163492 16.457637 16.219395 16.131696 15.615400 16.192678
45014 45220.0 16.720672 16.943064 16.309192 16.253035 16.215609 16.265908 16.544723 15.894155 16.471993 ... 16.184953 15.603858 16.208006 16.204189 16.194603 16.476746 16.159556 16.562719 15.660664 16.251187

45015 rows × 100 columns

In [100]:
id_ref = expression_data.iloc[:,0]
id_ref
Out[100]:
0            1.0
1            2.0
2            3.0
3            4.0
4            5.0
          ...   
45010    45216.0
45011    45217.0
45012    45218.0
45013    45219.0
45014    45220.0
Name: ID_REF, Length: 45015, dtype: float64
In [50]:
control = expression_data.iloc[:, 1:34]
infected = expression_data.iloc[:, 34:]
In [45]:
plt.hist(control.iloc[:, 1]) # !!!
plt.hist(control.iloc[:, 2], color='#FF0000AA')
Out[45]:
(array([4.4159e+04, 5.1800e+02, 1.7200e+02, 1.0400e+02, 4.6000e+01,
        7.0000e+00, 3.0000e+00, 2.0000e+00, 3.0000e+00, 1.0000e+00]),
 array([1.08000000e+00, 5.61009720e+04, 1.12200864e+05, 1.68300756e+05,
        2.24400648e+05, 2.80500540e+05, 3.36600432e+05, 3.92700324e+05,
        4.48800216e+05, 5.04900108e+05, 5.61000000e+05]),
 <BarContainer object of 10 artists>)
In [61]:
plt.hist(control.iloc[17, :])
plt.hist(infected.iloc[17, :], color='#FF0000AA')
Out[61]:
(array([ 6.,  3., 12.,  3.,  9.,  5.,  5.,  9.,  8.,  6.]),
 array([5.93781517, 6.1382851 , 6.33875504, 6.53922497, 6.7396949 ,
        6.94016484, 7.14063477, 7.3411047 , 7.54157464, 7.74204457,
        7.94251451]),
 <BarContainer object of 10 artists>)
In [81]:
a = [-1,0,0,1,2,2,3,3,3,3,4,4,5,6,7,7,8]
b = np.array([1,2,2,3,3,3,3,4,4,5]) +1
plt.hist(a)
plt.hist(b, color='#FF0000AA')
Out[81]:
(array([1., 0., 2., 0., 0., 4., 0., 2., 0., 1.]),
 array([2. , 2.4, 2.8, 3.2, 3.6, 4. , 4.4, 4.8, 5.2, 5.6, 6. ]),
 <BarContainer object of 10 artists>)
In [82]:
ttest_ind(a, b)
Out[82]:
Ttest_indResult(statistic=-0.7353361753098311, pvalue=0.4689776626892319)
In [86]:
np.array([[1,2],[3,4]]).ravel()
Out[86]:
array([1, 2, 3, 4])
In [102]:
i = 0
id_ref.iloc[i]
Out[102]:
1.0
In [103]:
sig_data = []

for i in range(expression_data.shape[0]):
    control_expression_i = control.iloc[i,:].values.ravel()
    infected_expression_i = infected.iloc[i,:].values.ravel()
    ref = id_ref.iloc[i]

    
    control_mean_i = control_expression_i.mean()
    infected_mean_i = infected_expression_i.mean()
    
    log2FC = infected_mean_i - control_mean_i
    
    _, p = ttest_ind(control_expression_i, infected_expression_i, equal_var=False)
    
    sig = False
    
    if p < 0.05 and abs(log2FC) > 1.5:
        sig = True
        
    sig_data.append([ref, control_mean_i, infected_mean_i, log2FC, p, sig])  

sig_df = pd.DataFrame(data=sig_data, columns=['ID_REF', 'CONTROL_MEAN', 'INFECTED_MEAN', 'LOG2FC', 'P', 'SIG'])
In [104]:
sig_df
Out[104]:
ID_REF CONTROL_MEAN INFECTED_MEAN LOG2FC P SIG
0 1.0 16.473359 16.253439 -0.219920 0.009117 False
1 2.0 1.293028 1.382012 0.088984 0.581833 False
2 3.0 1.003003 0.947770 -0.055233 0.448868 False
3 4.0 0.981373 0.963976 -0.017397 0.816628 False
4 5.0 0.983190 0.955862 -0.027329 0.718231 False
... ... ... ... ... ... ...
45010 45216.0 1.127839 0.969169 -0.158670 0.108337 False
45011 45217.0 1.370958 1.104557 -0.266401 0.214096 False
45012 45218.0 1.200348 1.021117 -0.179231 0.058063 False
45013 45219.0 16.517345 16.276750 -0.240596 0.001317 False
45014 45220.0 16.498123 16.285789 -0.212335 0.012087 False

45015 rows × 6 columns

In [105]:
gene_data = gpls.loc[:,['ID', 'GENE_SYMBOL']]
gene_data
Out[105]:
ID GENE_SYMBOL
0 1 NaN
1 2 NaN
2 3 NaN
3 4 NaN
4 5 NaN
... ... ...
45215 45216 NaN
45216 45217 NaN
45217 45218 NaN
45218 45219 NaN
45219 45220 NaN

45220 rows × 2 columns

In [107]:
gene_diff = gene_data.set_index('ID').join(sig_df.set_index('ID_REF'))
In [109]:
gene_diff.iloc[11]
Out[109]:
GENE_SYMBOL        APOBEC3B
CONTROL_MEAN        7.79091
INFECTED_MEAN       8.41112
LOG2FC             0.620215
P                0.00314838
SIG                   False
Name: 12, dtype: object
In [120]:
sig_exp = gene_diff[gene_diff['SIG'] == True].dropna()
In [121]:
sorted_sig_exp = sig_exp.sort_values(by=['LOG2FC'])
sorted_sig_exp
Out[121]:
GENE_SYMBOL CONTROL_MEAN INFECTED_MEAN LOG2FC P SIG
ID
4701 CAPNS2 10.361622 6.964108 -3.397514 3.567899e-18 True
6494 ANXA13 8.921412 6.294918 -2.626493 5.017886e-15 True
6748 CLRN3 7.397385 4.859611 -2.537774 1.952626e-14 True
991 TTR 8.851363 6.335043 -2.516320 4.249409e-17 True
27545 ANXA13 9.821740 7.325120 -2.496620 7.539615e-15 True
... ... ... ... ... ... ...
30072 PAEP 3.327891 7.303118 3.975227 1.057617e-19 True
1093 CSF3 1.538147 6.317711 4.779565 7.796504e-16 True
17062 FADS2 5.073183 9.929027 4.855844 1.382949e-22 True
39118 PRNT 1.659269 7.161729 5.502460 1.985034e-23 True
28337 EDAR 0.779498 6.825019 6.045520 1.172127e-24 True

433 rows × 6 columns

In [166]:
sorted_sig_exp[sorted_sig_exp['GENE_SYMBOL'] == 'DEFA4']
Out[166]:
GENE_SYMBOL CONTROL_MEAN INFECTED_MEAN LOG2FC P SIG
ID
In [122]:
up_regulated = sorted_sig_exp[-10:]
down_regulated = sorted_sig_exp[:10]
In [123]:
up_regulated
Out[123]:
GENE_SYMBOL CONTROL_MEAN INFECTED_MEAN LOG2FC P SIG
ID
44003 HHIP 3.137443 7.000847 3.863404 7.245861e-20 True
11787 RGR 5.132149 9.006631 3.874481 3.813493e-19 True
29465 CCR8 0.787747 4.728019 3.940272 6.288462e-20 True
35433 CREB5 6.609011 10.580027 3.971016 2.748210e-18 True
28117 IL28B 2.135472 6.107292 3.971820 1.915863e-17 True
30072 PAEP 3.327891 7.303118 3.975227 1.057617e-19 True
1093 CSF3 1.538147 6.317711 4.779565 7.796504e-16 True
17062 FADS2 5.073183 9.929027 4.855844 1.382949e-22 True
39118 PRNT 1.659269 7.161729 5.502460 1.985034e-23 True
28337 EDAR 0.779498 6.825019 6.045520 1.172127e-24 True
In [124]:
down_regulated
Out[124]:
GENE_SYMBOL CONTROL_MEAN INFECTED_MEAN LOG2FC P SIG
ID
4701 CAPNS2 10.361622 6.964108 -3.397514 3.567899e-18 True
6494 ANXA13 8.921412 6.294918 -2.626493 5.017886e-15 True
6748 CLRN3 7.397385 4.859611 -2.537774 1.952626e-14 True
991 TTR 8.851363 6.335043 -2.516320 4.249409e-17 True
27545 ANXA13 9.821740 7.325120 -2.496620 7.539615e-15 True
5020 LOC100132167 11.964410 9.672812 -2.291598 4.835899e-18 True
26799 MS4A8B 6.113347 3.858874 -2.254473 6.549172e-11 True
2169 F5 5.390926 3.260689 -2.130237 4.100002e-16 True
5853 LEFTY1 6.228430 4.161496 -2.066934 7.948616e-15 True
21599 CA9 8.475047 6.424245 -2.050802 3.440158e-12 True
In [162]:
genes = ['EDAR']
string_ids = stringdb.get_string_ids(genes)
enrichment_df = stringdb.get_enrichment(string_ids.queryItem)
In [163]:
string_ids
Out[163]:
queryItem queryIndex stringId ncbiTaxonId taxonName preferredName annotation
0 EDAR 0 9606.ENSP00000258443 9606 Homo sapiens EDAR Tumor necrosis factor receptor superfamily mem...
In [164]:
enrichment_df
Out[164]:
category term number_of_genes number_of_genes_in_background ncbiTaxonId inputGenes preferredNames p_value fdr description
0 Process GO:0042476 5 113 9606 WNT10A,EDAR,SHH,PAX9,EDA WNT10A,EDAR,SHH,PAX9,EDA 3.280000e-09 3.530000e-06 odontogenesis
1 Process GO:0060662 3 5 9606 EDAR,SHH,EDA EDAR,SHH,EDA 7.380000e-09 3.970000e-06 salivary gland cavitation
2 Process GO:0001942 4 72 9606 WNT10A,EDAR,SHH,EDA WNT10A,EDAR,SHH,EDA 6.780000e-08 1.820000e-05 hair follicle development
3 Process GO:0042481 3 26 9606 WNT10A,SHH,PAX9 WNT10A,SHH,PAX9 4.780000e-07 5.150000e-05 regulation of odontogenesis
4 Process GO:0009887 6 865 9606 WNT10A,EDAR,SHH,WNT10B,PAX9,EDA WNT10A,EDAR,SHH,WNT10B,PAX9,EDA 2.910000e-06 2.200000e-04 animal organ morphogenesis
... ... ... ... ... ... ... ... ... ... ...
198 RCTM HSA-373080 3 93 9606 WNT10A,SHH,WNT10B WNT10A,SHH,WNT10B 1.830000e-05 4.200000e-04 Class B/2 (Secretin family receptors)
199 RCTM HSA-3238698 2 26 9606 WNT10A,WNT10B WNT10A,WNT10B 1.100000e-04 1.200000e-03 WNT ligand biogenesis and trafficking
200 NetworkNeighborAL CL:11051 3 8 9606 EDAR,EDARADD,EDA EDAR,EDARADD,EDA 2.170000e-08 8.040000e-07 mixed, incl. salivary gland cavitation, and Ta...
201 NetworkNeighborAL CL:19550 2 10 9606 SLC45A2,SLC24A5 SLC45A2,SLC24A5 1.890000e-05 1.700000e-04 Melanin biosynthesis, and Sodium/potassium/cal...
202 NetworkNeighborAL CL:5745 2 12 9606 WNT10A,WNT10B WNT10A,WNT10B 2.600000e-05 1.900000e-04 WNT ligand biogenesis and trafficking

203 rows × 10 columns

In [167]:
# up_regulated.loc[:,'GENE_SYMBOL'].values.ravel() + down_regulated.loc[:,'GENE_SYMBOL'].values.ravel()

up_regulated.loc[:,'GENE_SYMBOL'].values

network = stringdb.get_network(['PGLYRP1', 'LEF1', 'CXCR5', 'CXCR1', 'CAMP', 'DEFA4', 'ELANE', 'LTF', 'CR2', 'ARG1', 'SLPI', 'MS4A1'])
network
Out[167]:
stringId_A stringId_B preferredName_A preferredName_B ncbiTaxonId score nscore fscore pscore ascore escore dscore tscore
0 9606.ENSP00000008938 9606.ENSP00000349446 PGLYRP1 ARG1 9606 0.917 0 0 0 0.213 0.000 0.9 0.000
1 9606.ENSP00000008938 9606.ENSP00000349446 PGLYRP1 ARG1 9606 0.917 0 0 0 0.213 0.000 0.9 0.000
2 9606.ENSP00000008938 9606.ENSP00000342082 PGLYRP1 SLPI 9606 0.922 0 0 0 0.000 0.000 0.9 0.255
3 9606.ENSP00000008938 9606.ENSP00000342082 PGLYRP1 SLPI 9606 0.922 0 0 0 0.000 0.000 0.9 0.255
4 9606.ENSP00000008938 9606.ENSP00000466090 PGLYRP1 ELANE 9606 0.932 0 0 0 0.191 0.000 0.9 0.227
5 9606.ENSP00000008938 9606.ENSP00000466090 PGLYRP1 ELANE 9606 0.932 0 0 0 0.191 0.000 0.9 0.227
6 9606.ENSP00000008938 9606.ENSP00000231751 PGLYRP1 LTF 9606 0.951 0 0 0 0.149 0.000 0.9 0.473
7 9606.ENSP00000008938 9606.ENSP00000231751 PGLYRP1 LTF 9606 0.951 0 0 0 0.149 0.000 0.9 0.473
8 9606.ENSP00000008938 9606.ENSP00000297435 PGLYRP1 DEFA4 9606 0.964 0 0 0 0.542 0.000 0.9 0.297
9 9606.ENSP00000008938 9606.ENSP00000297435 PGLYRP1 DEFA4 9606 0.964 0 0 0 0.542 0.000 0.9 0.297
10 9606.ENSP00000008938 9606.ENSP00000296435 PGLYRP1 CAMP 9606 0.968 0 0 0 0.560 0.000 0.9 0.337
11 9606.ENSP00000008938 9606.ENSP00000296435 PGLYRP1 CAMP 9606 0.968 0 0 0 0.560 0.000 0.9 0.337
12 9606.ENSP00000231751 9606.ENSP00000349446 LTF ARG1 9606 0.914 0 0 0 0.096 0.000 0.9 0.125
13 9606.ENSP00000231751 9606.ENSP00000349446 LTF ARG1 9606 0.914 0 0 0 0.096 0.000 0.9 0.125
14 9606.ENSP00000231751 9606.ENSP00000342082 LTF SLPI 9606 0.941 0 0 0 0.101 0.000 0.9 0.398
15 9606.ENSP00000231751 9606.ENSP00000342082 LTF SLPI 9606 0.941 0 0 0 0.101 0.000 0.9 0.398
16 9606.ENSP00000231751 9606.ENSP00000296435 LTF CAMP 9606 0.946 0 0 0 0.231 0.000 0.9 0.355
17 9606.ENSP00000231751 9606.ENSP00000296435 LTF CAMP 9606 0.946 0 0 0 0.231 0.000 0.9 0.355
18 9606.ENSP00000231751 9606.ENSP00000466090 LTF ELANE 9606 0.952 0 0 0 0.155 0.067 0.9 0.473
19 9606.ENSP00000231751 9606.ENSP00000466090 LTF ELANE 9606 0.952 0 0 0 0.155 0.067 0.9 0.473
20 9606.ENSP00000231751 9606.ENSP00000297435 LTF DEFA4 9606 0.963 0 0 0 0.284 0.000 0.9 0.526
21 9606.ENSP00000231751 9606.ENSP00000297435 LTF DEFA4 9606 0.963 0 0 0 0.284 0.000 0.9 0.526
22 9606.ENSP00000292174 9606.ENSP00000356024 CXCR5 CR2 9606 0.680 0 0 0 0.062 0.000 0.0 0.673
23 9606.ENSP00000292174 9606.ENSP00000356024 CXCR5 CR2 9606 0.680 0 0 0 0.062 0.000 0.0 0.673
24 9606.ENSP00000292174 9606.ENSP00000295683 CXCR5 CXCR1 9606 0.910 0 0 0 0.055 0.000 0.9 0.630
25 9606.ENSP00000292174 9606.ENSP00000295683 CXCR5 CXCR1 9606 0.910 0 0 0 0.055 0.000 0.9 0.630
26 9606.ENSP00000296435 9606.ENSP00000349446 CAMP ARG1 9606 0.917 0 0 0 0.137 0.000 0.9 0.121
27 9606.ENSP00000296435 9606.ENSP00000349446 CAMP ARG1 9606 0.917 0 0 0 0.137 0.000 0.9 0.121
28 9606.ENSP00000296435 9606.ENSP00000342082 CAMP SLPI 9606 0.964 0 0 0 0.088 0.000 0.9 0.645
29 9606.ENSP00000296435 9606.ENSP00000342082 CAMP SLPI 9606 0.964 0 0 0 0.088 0.000 0.9 0.645
30 9606.ENSP00000296435 9606.ENSP00000466090 CAMP ELANE 9606 0.980 0 0 0 0.320 0.379 0.9 0.586
31 9606.ENSP00000296435 9606.ENSP00000466090 CAMP ELANE 9606 0.980 0 0 0 0.320 0.379 0.9 0.586
32 9606.ENSP00000296435 9606.ENSP00000297435 CAMP DEFA4 9606 0.980 0 0 0 0.749 0.000 0.9 0.292
33 9606.ENSP00000296435 9606.ENSP00000297435 CAMP DEFA4 9606 0.980 0 0 0 0.749 0.000 0.9 0.292
34 9606.ENSP00000297435 9606.ENSP00000349446 DEFA4 ARG1 9606 0.914 0 0 0 0.168 0.000 0.9 0.050
35 9606.ENSP00000297435 9606.ENSP00000349446 DEFA4 ARG1 9606 0.914 0 0 0 0.168 0.000 0.9 0.050
36 9606.ENSP00000297435 9606.ENSP00000342082 DEFA4 SLPI 9606 0.916 0 0 0 0.000 0.000 0.9 0.195
37 9606.ENSP00000297435 9606.ENSP00000342082 DEFA4 SLPI 9606 0.916 0 0 0 0.000 0.000 0.9 0.195
38 9606.ENSP00000297435 9606.ENSP00000466090 DEFA4 ELANE 9606 0.987 0 0 0 0.784 0.000 0.9 0.456
39 9606.ENSP00000297435 9606.ENSP00000466090 DEFA4 ELANE 9606 0.987 0 0 0 0.784 0.000 0.9 0.456
40 9606.ENSP00000342082 9606.ENSP00000349446 SLPI ARG1 9606 0.906 0 0 0 0.061 0.000 0.9 0.081
41 9606.ENSP00000342082 9606.ENSP00000349446 SLPI ARG1 9606 0.906 0 0 0 0.061 0.000 0.9 0.081
42 9606.ENSP00000342082 9606.ENSP00000466090 SLPI ELANE 9606 0.997 0 0 0 0.061 0.870 0.9 0.833
43 9606.ENSP00000342082 9606.ENSP00000466090 SLPI ELANE 9606 0.997 0 0 0 0.061 0.870 0.9 0.833
44 9606.ENSP00000349446 9606.ENSP00000466090 ARG1 ELANE 9606 0.928 0 0 0 0.108 0.000 0.9 0.260
45 9606.ENSP00000349446 9606.ENSP00000466090 ARG1 ELANE 9606 0.928 0 0 0 0.108 0.000 0.9 0.260
46 9606.ENSP00000356024 9606.ENSP00000433277 CR2 MS4A1 9606 0.541 0 0 0 0.232 0.000 0.0 0.427
47 9606.ENSP00000356024 9606.ENSP00000433277 CR2 MS4A1 9606 0.541 0 0 0 0.232 0.000 0.0 0.427
In [168]:
genes = ['PGLYRP1', 'LEF1', 'CXCR5', 'CXCR1', 'CAMP', 'DEFA4', 'ELANE', 'LTF', 'CR2', 'ARG1', 'SLPI', 'MS4A1']
string_ids = stringdb.get_string_ids(genes)
enrichment_df = stringdb.get_enrichment(string_ids.queryItem)
In [169]:
enrichment_df
Out[169]:
category term number_of_genes number_of_genes_in_background ncbiTaxonId inputGenes preferredNames p_value fdr description
0 Process GO:0045321 12 894 9606 PGLYRP1,LTF,LEF1,CXCR5,CXCR1,CAMP,DEFA4,SLPI,A... PGLYRP1,LTF,LEF1,CXCR5,CXCR1,CAMP,DEFA4,SLPI,A... 9.000000e-17 9.490000e-14 leukocyte activation
1 Process GO:0006955 12 1560 9606 PGLYRP1,LTF,LEF1,CXCR5,CXCR1,CAMP,DEFA4,SLPI,A... PGLYRP1,LTF,LEF1,CXCR5,CXCR1,CAMP,DEFA4,SLPI,A... 6.910000e-14 2.430000e-11 immune response
2 Process GO:0006959 8 252 9606 PGLYRP1,LTF,CAMP,DEFA4,SLPI,CR2,MS4A1,ELANE PGLYRP1,LTF,CAMP,DEFA4,SLPI,CR2,MS4A1,ELANE 4.110000e-13 1.080000e-10 humoral immune response
3 Process GO:0002252 10 927 9606 PGLYRP1,LTF,LEF1,CXCR1,CAMP,DEFA4,SLPI,ARG1,CR... PGLYRP1,LTF,LEF1,CXCR1,CAMP,DEFA4,SLPI,ARG1,CR... 3.640000e-12 7.680000e-10 immune effector process
4 Process GO:0035821 7 182 9606 PGLYRP1,LTF,LEF1,CAMP,DEFA4,SLPI,ELANE PGLYRP1,LTF,LEF1,CAMP,DEFA4,SLPI,ELANE 5.310000e-12 9.350000e-10 modification of morphology or physiology of ot...
... ... ... ... ... ... ... ... ... ... ...
231 RCTM HSA-380108 2 48 9606 CXCR5,CXCR1 CXCR5,CXCR1 4.200000e-04 3.400000e-03 Chemokine receptors bind chemokines
232 RCTM HSA-977606 2 47 9606 CR2,ELANE CR2,ELANE 4.000000e-04 3.400000e-03 Regulation of Complement cascade
233 NetworkNeighborAL CL:10260 6 11 9606 PGLYRP1,LTF,CAMP,DEFA4,SLPI,ELANE PGLYRP1,LTF,CAMP,DEFA4,SLPI,ELANE 1.460000e-16 5.690000e-15 specific granule lumen
234 NetworkNeighborAL CL:10261 5 6 9606 PGLYRP1,CAMP,DEFA4,SLPI,ELANE PGLYRP1,CAMP,DEFA4,SLPI,ELANE 1.520000e-14 1.490000e-13 mixed, incl. Cathelicidin, and neutrophil medi...
235 NetworkNeighborAL CL:7364 2 24 9606 CXCR5,CXCR1 CXCR5,CXCR1 1.100000e-04 4.800000e-04 CXC Chemokine domain, and dendritic cell chemo...

236 rows × 10 columns