首先,对不起,如果这是有点长,但我想充分描述我有什么问题,我已经尝试过了。
我正在尝试在多个条件下连接(合并)两个数据框架对象。如果要满足的条件都是'等于'操作符,我知道如何做到这一点,但是,我需要使用LESS THAN和MORE THAN。
数据框代表遗传信息:一个是基因组中的突变列表(称为snp),另一个提供了基因在人类基因组上的位置信息。对这些语句执行df.head()将返回以下结果:
SNP DataFrame (snp_df):
chromosome SNP BP
0 1 rs3094315 752566
1 1 rs3131972 752721
2 1 rs2073814 753474
3 1 rs3115859 754503
4 1 rs3131956 758144
显示SNP引用ID及其位置。'BP'代表'Base-Pair'的位置。
基因数据帧(gene_df):
chromosome chr_start chr_stop feature_id
0 1 10954 11507 GeneID:100506145
1 1 12190 13639 GeneID:100652771
2 1 14362 29370 GeneID:653635
3 1 30366 30503 GeneID:100302278
4 1 34611 36081 GeneID:645520
该数据框显示了所有感兴趣的基因的位置。
我想找出的是基因组中基因区域内的所有snp,并丢弃这些区域外的snp。
如果我想基于多个(等于)条件合并两个数据帧,我将做如下操作:
merged_df = pd.merge(snp_df, gene_df, on=['chromosome', 'other_columns'])
然而,在这种情况下,我需要找到染色体值与基因数据框中值匹配的snp, BP值落在'chr_start'和'chr_stop'之间。这一挑战在于这些数据帧非常大。在当前数据集中,snp_df有6795021行,gene_df有34362行。
我试图通过分别观察染色体或基因来解决这个问题。由于没有使用性染色体,有22个不同的染色体值(整数1-22)。这两种方法都需要很长时间。一种方法是使用pandasql
模块,而另一种方法是通过单独的基因进行循环。
import pandas as pd
import pandasql as psql
pysqldf = lambda q: psql.sqldf(q, globals())
q = """
SELECT s.SNP, g.feature_id
FROM this_snp s INNER JOIN this_genes g
WHERE s.BP >= g.chr_start
AND s.BP <= g.chr_stop;
"""
all_dfs = []
for chromosome in snp_df['chromosome'].unique():
this_snp = snp_df.loc[snp_df['chromosome'] == chromosome]
this_genes = gene_df.loc[gene_df['chromosome'] == chromosome]
genic_snps = pysqldf(q)
all_dfs.append(genic_snps)
all_genic_snps = pd.concat(all_dfs)
基因迭代法
all_dfs = []
for line in gene_df.iterrows():
info = line[1] # Getting the Series object
this_snp = snp_df.loc[(snp_df['chromosome'] == info['chromosome']) &
(snp_df['BP'] >= info['chr_start']) & (snp_df['BP'] <= info['chr_stop'])]
if this_snp.shape[0] != 0:
this_snp = this_snp[['SNP']]
this_snp.insert(len(this_snp.columns), 'feature_id', info['feature_id'])
all_dfs.append(this_snp)
all_genic_snps = pd.concat(all_dfs)
有谁能给出更有效的方法吗?
我刚刚想到了一个解决这个问题的方法——将我的两个方法结合起来:
首先,关注单个染色体,然后遍历这些较小数据框中的基因。这也不需要使用任何SQL查询。我还包括了一个部分,以立即识别任何多余的基因,没有任何snp在其范围内。这利用了我通常尽量避免的双for循环,但在这种情况下,它工作得很好。
all_dfs = []
for chromosome in snp_df['chromosome'].unique():
this_chr_snp = snp_df.loc[snp_df['chromosome'] == chromosome]
this_genes = gene_df.loc[gene_df['chromosome'] == chromosome]
# Getting rid of redundant genes
min_bp = this_chr_snp['BP'].min()
max_bp = this_chr_snp['BP'].max()
this_genes = this_genes.loc[~(this_genes['chr_start'] >= max_bp) &
~(this_genes['chr_stop'] <= min_bp)]
for line in this_genes.iterrows():
info = line[1]
this_snp = this_chr_snp.loc[(this_chr_snp['BP'] >= info['chr_start']) &
(this_chr_snp['BP'] <= info['chr_stop'])]
if this_snp.shape[0] != 0:
this_snp = this_snp[['SNP']]
this_snp.insert(1, 'feature_id', info['feature_id'])
all_dfs.append(this_snp)
all_genic_snps = pd.concat(all_dfs)
虽然它运行得不是特别快,但它确实运行了,所以我可以得到一些答案。我还是想知道有没有人有什么建议可以让它更有效地运行。
您可以使用以下命令来完成您要查找的内容:
merged_df=snp_df.merge(gene_df,on=['chromosome'],how='inner')
merged_df=merged_df[(merged_df.BP>=merged_df.chr_start) & (merged_df.BP<=merged_df.chr_stop)][['SNP','feature_id']]
注意:您的示例数据帧不符合您的连接条件。下面是一个使用修改过的数据框架的例子:
snp_df
Out[193]:
chromosome SNP BP
0 1 rs3094315 752566
1 1 rs3131972 30400
2 1 rs2073814 753474
3 1 rs3115859 754503
4 1 rs3131956 758144
gene_df
Out[194]:
chromosome chr_start chr_stop feature_id
0 1 10954 11507 GeneID:100506145
1 1 12190 13639 GeneID:100652771
2 1 14362 29370 GeneID:653635
3 1 30366 30503 GeneID:100302278
4 1 34611 36081 GeneID:645520
merged_df
Out[195]:
SNP feature_id
8 rs3131972 GeneID:100302278