Pandas，让Python像R一样处理数据，但快

%matplotlib inline

#import plotly
#plotly.offline.init_notebook_mode()

import matplotlib
matplotlib.style.use('ggplot')
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import os
from glob import glob

dir_1 = "pandas_data/"
glob(dir_1+'*')

['pandas_data/ENCFF289HGQ.tsv',
 'pandas_data/gencode.v24.ENS2SYN',
 'pandas_data/ENCFF262OBL.tsv',
 'pandas_data/Gene_metadata_primary_wt_whole_cell.tsv',
 'pandas_data/ENCFF673KYR.tsv',
 'pandas_data/ENCFF060LPA.tsv']

!head -n 4 pandas_data/gencode.v24.ENS2SYN

gene_id    gene_symbol

ENSG00000223972.5    DDX11L1

ENSG00000227232.5    WASH7P

ENSG00000278267.1    MIR6859-1

!head -n 4 pandas_data/ENCFF556YSD.tsv

transcript_id    gene_id    length    effective_length    expected_count    TPM    FPKM

ENST00000373020.4    ENSG00000000003.10    2206    1925.57    997.23    1.71    7.21

ENST00000494424.1    ENSG00000000003.10    820    539.58    24.77    0.15    0.64

ENST00000496771.1    ENSG00000000003.10    1025    744.57    0.00    0.00    0.00

ens2syn_file = "pandas_data/gencode.v24.ENS2SYN"

# pandas中的计数都是从0开始的
# header=0: 指定第一行包含列的名字
# index_col=0: 指定第一列为行的名字
ens2syn = pd.read_table(ens2syn_file, header=0, index_col=0)

ens2syn[:3]

ens2syn[ens2syn.index=="ENSG00000227232.5"]

ens2syn[ens2syn['gene_symbol'].isin(['DDX11L1','MIR6859-1'])]

# head: 只展示部分数据
ens2syn[ens2syn.index.str.contains(r'ENSG0000022')].head()

tsvL = glob(dir_1+'ENC*.tsv')
tsvL

['pandas_data/ENCFF289HGQ.tsv',
 'pandas_data/ENCFF262OBL.tsv',
 'pandas_data/ENCFF673KYR.tsv',
 'pandas_data/ENCFF060LPA.tsv']

index = 0
tsvFile = tsvL[index]
expr = pd.read_table(tsvFile, header=0, index_col=0)
expr.head(3)

expr[['FPKM','TPM']].head(3)

# 因为要把多个文件的同一类型表达值合并到一个文件，我们使用文件名作为列的名字
name = os.path.split(tsvFile)[-1][:-4]
print name
expr_tpm = expr.loc[:,['TPM']] # 取出所有的行和名字为TPM的列
expr_tpm.columns=[name]  
expr_tpm[:3]

ENCFF289HGQ

# 为了读取多个文件，定义一个函数简化操作
def readExpr_1(tsvFileL, typeL=['TPM','FPKM']):
    '''
    tsvFileL: lists of files waiting for reading
    resultD: a dictionary to save data matrix
            {'TPM':[mat1, mat2,...]
             'FPKM':[mat1, mat2, ...]}
    typeL; list of names for columns to be extracted
    '''
    resultD = {}
    for _type in typeL: resultD[_type] = []

    for tsvFile in tsvFileL:
        expr = pd.read_table(tsvFile, header=0, index_col=0)
        name = os.path.split(tsvFile)[-1][:-4]  #this options is very arbitary
        for _type in typeL: # add _ to type to avoid override Python inner function `type` 
            expr_type = expr.loc[:,[_type]]
            expr_type.columns = [name]
            resultD[_type].append(expr_type)
    return resultD
#-----------------------------------------------------

exprD = readExpr_1(tsvL)
TPM_mat = exprD['TPM']
FPKM_mat = exprD['FPKM']

# 选取第一个矩阵
_idL = ['ENSG00000000003.14', 'ENSG00000000005.5','ENSG00000000419.12','ENSG00000000457.13']
mat1 = TPM_mat[0]
mat1 = mat1[mat1.index.isin(_idL)]
mat1

# 选取第二个矩阵
_idL = ['ENSG00000001561.6','ENSG00000000003.14', 'ENSG00000000419.12','ENSG00000001036.13']
mat2 = TPM_mat[1]
mat2 = mat2[mat2.index.isin(_idL)]
mat2

* outer: 合并所有的索引，缺失值填充NA
* inner：保留共有的索引
* left：使用第一个矩阵的索引
* right：使用第二个矩阵的索引

pd.merge(mat1, mat2, left_index=True, right_index=True, how="outer")

pd.merge(mat1, mat2, left_index=True, right_index=True, how="inner")

pd.merge(mat1, mat2, left_index=True, right_index=True, how="left")

pd.concat([mat1, mat2], axis=1)

pd.concat([mat1, mat2], axis=1, join="inner")

mat3 = mat1.join(mat2, how="outer")
mat3

mat3 = mat3.fillna(0)
mat3

#Both works well here
#mat3[(mat3>0).any(axis=1)]
mat3.loc[(mat3>0).any(axis=1)]

%timeit test_merge = reduce(lambda left,right: pd.merge(left,right,left_index=True,right_index=True,how='outer'), TPM_mat)

100 loops, best of 3: 3.36 ms per loop

%timeit test_merge = pd.concat(TPM_mat, axis=1)

1000 loops, best of 3: 1.21 ms per loop

%timeit TPM_mat[0].join(TPM_mat[1:], how="outer")

1000 loops, best of 3: 1.25 ms per loop

%load_ext memory_profiler

%memit test_merge = reduce(lambda left,right: pd.merge(left,right,left_index=True,right_index=True,how='outer'), TPM_mat)

peak memory: 107.32 MiB, increment: 0.01 MiB

%memit test_merge = pd.concat(TPM_mat, axis=1)

peak memory: 107.30 MiB, increment: 0.10 MiB

%memit TPM_mat[0].join(TPM_mat[1:], how="outer")

peak memory: 107.32 MiB, increment: 0.00 MiB

# 读取多个文件，并且合并矩阵，定义一个函数简化操作
def concatExpr(tsvFileL, typeL=['TPM','FPKM']):
    '''
    tsvFileL: lists of files waiting for reading
    resultD: a dictionary to save data matrix
            {'TPM':[mat1, mat2,...]
             'FPKM':[mat1, mat2, ...]}
    typeL; list of names for columns to be extracted
    '''
    resultD = {}
    for _type in typeL: resultD[_type] = []

    for tsvFile in tsvFileL:
        expr = pd.read_table(tsvFile, header=0, index_col=0)
        name = os.path.split(tsvFile)[-1][:-4]  #this options is very arbitary
        for _type in typeL: # add _ to type to avoid override Python inner function `type` 
            expr_type = expr.loc[:,[_type]]
            expr_type.columns = [name]
            resultD[_type].append(expr_type)
    #-------------------------------------------
    mergeD = {}
    for _type in typeL:
        mergeM = pd.concat(resultD[_type], axis=1)
        mergeM = mergeM.fillna(0) # Substitute all NA with 0
        mergeM = mergeM.loc[(mergeM>0).any(axis=1)] # Delete aoo zero rows.
        mergeD[_type] = mergeM
    return mergeD
#-----------------------------------------------------

typeL = ['TPM','FPKM']
exprD = concatExpr(tsvL, typeL)
TPM_mat = exprD['TPM']
FPKM_mat = exprD['FPKM']

TPM_mat.head()

# 回顾下数据格式
ens2syn.head(3)

ens2syn.shape

(60725, 1)

ens2syn = ens2syn[ens2syn.index.isin(TPM_mat.index)]

ens2syn.shape

(48, 1)

ens2syn.head(3)

meta = "pandas_data/meta.tsv"
metaM = pd.read_table(meta, header=0, index_col=0)
# 重名了列的名字
oriColnames = metaM.columns.values
nameD = dict([(i,i.replace(' ','_')) for i in oriColnames])
metaM.rename(columns=nameD, inplace=True)
metaM.head(3)

sampleL = TPM_mat.columns.values
metaM = metaM[metaM.index.isin(sampleL)]
# 同时索引行和列
metaM.ix[:4,:5]

# 假如只提取`Biosample`开头的列
#meta_colL = ['Biosample term id', 'Biosample term name']

# Extract columns matching specific patterns
# Both works well, filter is more simple
#metaM.loc[:,metaM.columns.str.contains(r'^Biosample')]
metaM = metaM.filter(regex=("^Biosample"))
metaM

metaM.fillna('')

metaM.to_csv("pandas_data/meta2.tsv", sep="\t")

ens2syn.to_csv("pandas_data/gencode.v24.ENS2SYN", sep="\t")

TPM_mat.to_csv("pandas_data/TPM", sep=b'\t', float_format="%.2f")

# 写入模式打开一个HDF5文件,使用压缩格式以节省空间
store = pd.HDFStore("pandas_data/ENCODE.hdf5", "w", complib=str("zlib"), complevel=9)

# 写入表达矩阵
store["TPM"] = TPM_mat
store["FPKM"] = FPKM_mat

# 写入注释文件
store['ens2syn'] = ens2syn
store['meta'] = metaM

# 关闭HDF5句柄
store.close()

/MPATHB/soft/anacond/lib/python2.7/site-packages/IPython/core/interactiveshell.py:3035: PerformanceWarning: 
your performance may suffer as PyTables will pickle object types that it cannot
map directly to c-types [inferred_type->mixed,key->block0_values] [items->['Biosample term id', 'Biosample term name', 'Biosample type', 'Biosample life stage', 'Biosample sex', 'Biosample organism', 'Biosample Age']]

  exec(code_obj, self.user_global_ns, self.user_ns)

# 写入模式打开一个HDF5文件,使用压缩格式已节省空间
store = pd.HDFStore("pandas_data/ENCODE.hdf5", "w", complib=str("zlib"), complevel=9)

# 写入表达矩阵
store["TPM"] = TPM_mat
store["FPKM"] = FPKM_mat

# 写入注释文件
store['ens2syn'] = ens2syn
store['meta'] = metaM.fillna('')

# 关闭HDF5句柄
store.close()

store = pd.HDFStore("pandas_data/ENCODE.hdf5")

# 列出HDF5文件的索引名字
store.keys()

['/FPKM', '/TPM', '/ens2syn', '/meta']

TPM_mat = store['TPM']
TPM_mat.head(3)

ens2syn = store['ens2syn']
meta = store['meta']

store.close()

targetL = ['KRIT1','AK2']

ensID = ens2syn[ens2syn["gene_symbol"].isin(targetL)]
ensID

targetExpr = TPM_mat[TPM_mat.index.isin(ensID.index)]
targetExpr

ensID_dict = ensID.to_dict()
ensID_dict

{'gene_symbol': {'ENSG00000001631.14': 'KRIT1', 'ENSG00000004455.16': 'AK2'}}

targetExpr = targetExpr.rename(index=ensID_dict['gene_symbol'])
targetExpr

targetExpr_t = targetExpr.T
targetExpr_t

meta_type = ["Biosample term name","Biosample type", "Biosample life stage", "Biosample sex"]

meta = meta[meta_type]
meta

meta.loc['ENCFF673KYR',"Biosample life stage"] = "fetal"
# Much faster
meta = meta.set_value('ENCFF673KYR','Biosample sex','female')
meta = meta.set_value('ENCFF289HGQ','Biosample sex','female')
meta

target_expr_meta = targetExpr_t.join(meta, how="left")
target_expr_meta

oriColnames = target_expr_meta.columns.values
nameD = dict([(i,i.replace(' ','_')) for i in oriColnames])
target_expr_meta.rename(columns=nameD, inplace=True)
target_expr_meta

target_expr_meta.plot.scatter(x='KRIT1', y='AK2')

<matplotlib.axes._subplots.AxesSubplot at 0x7fbcaefc0c10>

/MPATHB/soft/anacond/lib/python2.7/site-packages/matplotlib/collections.py:590: FutureWarning:

elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison

a = target_expr_meta.boxplot(["KRIT1", "AK2"])

/MPATHB/soft/anacond/lib/python2.7/site-packages/IPython/kernel/__main__.py:1: FutureWarning:

The default value for 'return_type' will change to 'axes' in a future release.
 To use the future behavior now, set return_type='axes'.
 To keep the previous behavior and silence this warning, set return_type='dict'.

fig = {
    'data': [
        {
            'x': target_expr_meta.KRIT1,
            'y': target_expr_meta.AK2,
            'text': target_expr_meta.Biosample_term_name,
            'mode': 'markers+texts',
            'name': 'Legend name',
            'marker': {
                'symbol':"circle",
                "opacity": "0.9"
            }
        },
        {
            'x': [0,40],
            'y': [0,40],
            'text': "Diagonal",
            'mode': 'lines',
            'name': 'Diagonal line',
            'showlegend': False,
            'line': {
                'color': ('rgb(192,192,192)')
            }
        }
    ],
    'layout': {
        'xaxis': {'title':'KRIT1 expression value','range':[0,40]},
        'yaxis': {'title':'AK2 expression value','range':[0,40]},
        'height':500,
        'width':600,
        'showlegend': True,
         "legend": {
            "x": 0.6,
            "y": 1
        }
    }
}
plotly.offline.iplot(fig)