基于Pandas/Python按Position阈值划分可变窗口计算斜率

我需要在基因数据集中按可变窗口计算斜率，窗口大小由Position列的数值距离阈值决定。数据集df包含Scaffold、centiMorgan、Position三列，其中centiMorgan和Position为连续变量。

规则为：当相邻行Position的差值超过设定阈值（如20）时，当前窗口结束，计算该窗口首尾点的斜率（公式：m=(Y2-Y1)/(X2-X1)，Y为centiMorgan，X为Position），并开启新窗口；若差值未超过阈值，则继续扩展窗口。

示例数据集：

Scaffold    centiMorgan    Position
Scaffold01     0.0           10004
Scaffold01     0.1           10006 
Scaffold01     0.5           10008
Scaffold02     1.5           10450
Scaffold02     2.9           11100
Scaffold02     3.0           11102
Scaffold03     3.8           12600
Scaffold04     4.6           12610

我尝试用Pandas的.diff方法计算Position差值，再通过嵌套循环划分窗口，代码如下：

import pandas as pd

df = pd.dataframe()

threshold = 20


# Subset position column:
posCol = ['position']

# create a new column 'Dist' in original df with results of .diff
# use .fillna(0) to fill in 0 where there is no previous row

df['Dist'] = posCol.diff().fillna(0)

# use a nested for/while loop to read through every line in df[Dist], the result column from .diff

for dist in df[Dist]:

     while dist < threshold 
        # write entire row to a new dataframe called window1
        # extract the first line of window1 
        # extract the last line of window1 
        # calculate the slope of window1 using the first and last points
        if dist > 20:
             break 
   
# Move on to the next window and repeat. 

尚未运行代码，想确认该方案是否合理或有更优实现，恳请提供建议。

原方案的问题 #

1. 语法错误：pd.dataframe()应为pd.DataFrame()；posCol = ['position']是列表，无法直接调用.diff()，需使用df['Position']；循环缺少冒号且无更新逻辑，会陷入死循环。
2. 效率低下：嵌套循环在处理大型基因数据集时性能极差，Pandas的向量化操作远优于循环。

优化实现代码 #

import pandas as pd

# 示例数据集（实际使用时替换为你的数据集）
data = {
    'Scaffold': ['Scaffold01', 'Scaffold01', 'Scaffold01', 'Scaffold02', 'Scaffold02', 'Scaffold02', 'Scaffold03', 'Scaffold04'],
    'centiMorgan': [0.0, 0.1, 0.5, 1.5, 2.9, 3.0, 3.8, 4.6],
    'Position': [10004, 10006, 10008, 10450, 11100, 11102, 12600, 12610]
}
df = pd.DataFrame(data)

threshold = 20

# 1. 按Scaffold和Position排序，确保窗口划分的顺序正确
df = df.sort_values(['Scaffold', 'Position']).reset_index(drop=True)

# 2. 计算同一Scaffold内相邻Position的差值，不同Scaffold的差值设为超过阈值（强制分割窗口）
df['diff_pos'] = df.groupby('Scaffold')['Position'].diff().fillna(threshold + 1)

# 3. 生成窗口ID：当差值超过阈值时标记为分割点，累加得到唯一窗口ID
df['window_id'] = (df['diff_pos'] > threshold).cumsum()

# 4. 定义斜率计算函数，处理单一行窗口的边界情况
def calculate_slope(group):
    if len(group) == 1:
        # 窗口仅一行时无法计算斜率，返回NaN，可根据需求调整
        return pd.Series({
            'slope': float('nan'),
            'start_pos': group['Position'].iloc[0],
            'end_pos': group['Position'].iloc[0],
            'start_cM': group['centiMorgan'].iloc[0],
            'end_cM': group['centiMorgan'].iloc[0]
        })
    x1, y1 = group['Position'].iloc[0], group['centiMorgan'].iloc[0]
    x2, y2 = group['Position'].iloc[-1], group['centiMorgan'].iloc[-1]
    # 避免除以0的情况
    slope = (y2 - y1) / (x2 - x1) if (x2 - x1) != 0 else float('nan')
    return pd.Series({
        'slope': slope,
        'start_pos': x1,
        'end_pos': x2,
        'start_cM': y1,
        'end_cM': y2
    })

# 按Scaffold和窗口ID分组计算斜率
window_slopes = df.groupby(['Scaffold', 'window_id']).apply(calculate_slope).reset_index()

print(window_slopes)

代码优势 #

• 高效性：使用Pandas向量化和分组操作，避免循环，处理百万级数据时性能提升明显。
• 准确性：自动处理不同Scaffold的窗口分割（基因数据中不同Scaffold属于独立序列），同时处理单一行窗口的边界情况。
• 可读性：结果清晰展示每个窗口的所属Scaffold、窗口ID、斜率及首尾点的Position和centiMorgan值，便于后续分析。

内容的提问来源于stack exchange，提问作者Daniela Zárate