基因组覆盖度滑动窗口计算：基于bwa mem与samtools depth结果

Hey there! Let's work through how to calculate sliding window average coverage (window size 3, step 2) grouped by contig using roll-based methods—pandas' rolling functionality is exactly what you need here. Here's a step-by-step, actionable solution:

First, confirm your samtools depth output is a tab-separated file with 3 columns:

• Contig name (grouping key)
• 1-based base position
• Coverage depth (number of aligned reads)

This approach uses pandas' groupby + rolling to handle contig-specific sliding window calculations, and fixes common edge cases like non-consecutive positions.

Full Code Example #

import pandas as pd

def process_contig_coverage(group):
    # Ensure positions are sorted (critical for sliding window logic)
    group_sorted = group.sort_values("pos").reset_index(drop=True)
    
    # Fill in missing positions (samtools depth skips bases with 0 coverage by default)
    # This ensures we don't skip gaps in the genome
    full_pos_range = pd.Series(
        range(group_sorted["pos"].min(), group_sorted["pos"].max() + 1),
        name="pos"
    )
    group_filled = pd.merge(
        full_pos_range, group_sorted, on="pos", how="left"
    ).fillna({"depth": 0, "contig": group_sorted["contig"].iloc[0]})
    
    # Calculate sliding window average: window=3, step=2, keep partial windows for short contigs
    group_filled["window_mean_depth"] = group_filled["depth"].rolling(
        window=3,
        step=2,
        min_periods=1  # Compute mean even if window is smaller than 3 (e.g., end of contig)
    ).mean()
    
    # Add window position metadata (adjust based on your preferred window labeling)
    # Here we use the middle position of each window as the representative position
    group_filled["window_mid_pos"] = group_filled["pos"].shift(-1)  # For window [n, n+1, n+2], mid is n+1
    # Filter to keep only valid window rows (matches step=2 interval)
    window_results = group_filled.dropna(subset=["window_mean_depth", "window_mid_pos"]).iloc[::2].copy()
    
    # Clean up and keep only relevant columns
    return window_results[["contig", "window_mid_pos", "window_mean_depth"]]

# Load your samtools depth output
depth_df = pd.read_csv(
    "your_depth_file.tsv",
    sep="\t",
    names=["contig", "pos", "depth"]
)

# Apply processing to each contig group
final_window_df = depth_df.groupby("contig").apply(process_contig_coverage).reset_index(drop=True)

# Save the result to a new file
final_window_df.to_csv("windowed_coverage_mean.tsv", sep="\t", index=False)

• Missing Positions: The code fills in bases with 0 coverage that samtools depth omits—this ensures your sliding window doesn't skip gaps in the genome. If you don't need this, you can remove the full_pos_range and group_filled steps.
• Window Labeling: The example uses the middle position of each window as the reference point. If you prefer start/end positions instead, modify the window_mid_pos line to calculate those (e.g., window_start = group_filled["pos"] for window start positions).
• Partial Windows: The min_periods=1 parameter ensures we calculate an average even if the contig ends before a full window of 3 bases (e.g., a contig with only 2 bases will return the mean of those 2).

If you don't want to use Python, you can combine bedtools makewindows and bedtools map to achieve the same result:

1. Convert your depth file to BED format (add a 4th column for depth)
2. Generate windows with bedtools makewindows -g your_genome.fai -w 3 -s 2
3. Map coverage to windows with bedtools map -a windows.bed -b depth.bed -c 4 -o mean

But since you asked for roll-based methods, the pandas approach is the most direct fit.

内容的提问来源于stack exchange，提问作者rororo