一种有效的划分大范围fi子集的方法

Question

我以前曾发布过一个问题，用R:subset recursively a data.frame来解决，但是这个文件太大了，我需要大量的时间和RAM内存来读取它。我想知道我是否可以在python中使用pandas来做同样的事情，因为我是python的新手，pandas看起来更像R，至少在sintax中是这样。以下是我前一篇文章的摘要：

上一篇文章： 我有一个以制表符分隔的文件，它的行数接近1500万，大小为27GB。我需要一种基于两个标准的有效的数据子集方法。我可以这样做是一个for循环，但我想知道是否有一种更优雅的方式来完成这项工作，而且显然效率更高。这个数据帧看起来像这样：

SNP         CHR     BP          P
rs1000000   chr1    126890980   0.000007
rs10000010  chr4    21618674    0.262098    
rs10000012  chr4    1357325     0.344192
rs10000013  chr4    37225069    0.726325    
rs10000017  chr4    84778125    0.204275    
rs10000023  chr4    95733906    0.701778
rs10000029  chr4    138685624   0.260899
rs1000002   chr3    183635768   0.779574
rs10000030  chr4    103374154   0.964166    
rs10000033  chr2    139599898   0.111846    
rs10000036  chr4    139219262   0.564791
rs10000037  chr4    38924330    0.392908    
rs10000038  chr4    189176035   0.971481    
rs1000003   chr3    98342907    0.000004
rs10000041  chr3    165621955   0.573376
rs10000042  chr3    5237152     0.834206    
rs10000056  chr4    189321617   0.268479
rs1000005   chr1    34433051    0.764046
rs10000062  chr4    5254744     0.238011    
rs10000064  chr4    127809621   0.000044
rs10000068  chr2    36924287    0.000003
rs10000075  chr4    179488911   0.100225    
rs10000076  chr4    183288360   0.962476
rs1000007   chr2    237752054   0.594928
rs10000081  chr1    17348363    0.517486    
rs10000082  chr1    167310192   0.261577    
rs10000088  chr1    182605350   0.649975
rs10000092  chr4    21895517    0.000005
rs10000100  chr4    19510493    0.296693    

首先我要做的是选择那些p值低于阈值的SNP，然后按CHR和BP对这个子集进行排序。一旦我有了这个子集，我需要从重要的SNP上下获取500000个窗口中的所有SNP，这个步骤将定义一个区域。我需要为所有重要的SNP做这件事，并将每个区域存储到一个列表或类似的东西中进行进一步的分析。例如，在所显示的数据帧中，CHR==chr1的最高有效SNP（即低于阈值0.001）为rs1000000，CHR==chr4的最高有效SNP为rs1000002。因此，这两个SNP将定义两个区域，我需要从每个最重要的SNP的POS上下提取属于500000个区域的SNP。在

@eddi和@rafaelpereira提供的R代码解决方案如下：

^{pr2}$

Answer 1

首先，我强烈建议从CSV文件切换到PyTables（HDF存储），并存储按['SNP','BP']排序的DF，因为它速度更快，允许条件选择（参见where参数），并且通常占用较少的空间-请参见this comparison。在

HDF store recipes

下面是一个工作示例脚本，它执行以下操作：

1. 生成样本DF（20M行，8列：'SNP', 'CHR', 'BP', 'P', 'SNP2', 'CHR2', 'BP2', 'P2'）。我故意将列数翻倍，因为我认为您的CSV有更多的列，因为我生成的具有20M行和8列的CSV文件的大小只有1.7GB。在
2. 将生成的DF保存到CSV文件（文件大小：1.7 GB）
3. 以1M行分块读取CSV文件（前4列只读）
4. 按['CHR','BP']排序数据并将结果另存为PyTable（.h5）
5. 从HDF读取只存储那些P < threshold
6. 从HDF读取存储min(SNP) - 500K和max(SNP) + 500K之间的所有行-您可能需要改进这一部分

代码：

import numpy as np
import pandas as pd

##############################
# generate sample DF
##############################
rows = 2*10**7
chr_lst = ['chr{}'.format(i) for i in range(1,10)]

df = pd.DataFrame({'SNP': np.random.randint(10**6, 10**7, rows).astype(str)})
df['CHR'] = np.random.choice(chr_lst, rows)
df['BP'] = np.random.randint(10**6, 10**9, rows)
df['P'] = np.random.rand(rows)
df.SNP = 'rs' + df.SNP

"""
# NOTE: sometimes it gives me MemoryError !!!
# because of that i did it "column-by-column" before
df = pd.DataFrame({
    'SNP': np.random.randint(10**6, 10**7, rows).astype(str),
    'CHR': np.random.choice(chr_lst, rows),
    'BP':  np.random.randint(10**6, 10**9, rows),
    'P':   np.random.rand(rows)
}, columns=['SNP','CHR','BP','P'])

df.SNP = 'rs' + df.SNP
"""

# make 8 columns out of 4 ...
df = pd.concat([df]*2, axis=1)
df.columns = ['SNP', 'CHR', 'BP', 'P', 'SNP2', 'CHR2', 'BP2', 'P2']

##############################
# store DF as CSV file
##############################
csv_path = r'c:/tmp/file_8_cols.csv'
df.to_csv(csv_path, index=False)

##############################
# read CSV file (only needed cols) in chunks
##############################
csv_path = r'c:/tmp/file_8_cols.csv'
cols = ['SNP', 'CHR', 'BP', 'P']
chunksize = 10**6

df = pd.concat([x for x in pd.read_csv(csv_path, usecols=cols,
                                       chunksize=chunksize)],
               ignore_index=True )

##############################
# sort DF and save it as .h5 file
##############################
store_path = r'c:/tmp/file_sorted.h5'
store_key = 'test'
(df.sort_values(['CHR','BP'])
   .to_hdf(store_path, store_key, format='t', mode='w', data_columns=True)
)


##############################
# read HDF5 file in chunks
##############################
store_path = r'c:/tmp/file_sorted.h5'
store_key = 'test'
chunksize = 10**6

store = pd.HDFStore(store_path)

threshold = 0.001
store_condition = 'P < %s' % threshold

i = store.select(key=store_key, where=store_condition)

# select all rows between `min(SNP) - 500K` and `max(SNP) + 500K`
window_size = 5*10**5
start = max(0, i.index.min() - window_size)
stop = min(store.get_storer(store_key).nrows, i.index.max() + window_size)
df = pd.concat([
        x for x in store.select(store_key, chunksize=chunksize,
                                start=start, stop=stop, )
               ])

# close the store before exiting...
store.close()

样本数据：

^{pr2}$