0.脚本由大模型生成

1. 鉴定ROH

plink --bfile judaei -allow-extra-chr --chr-set 30 --homozyg-window-snp 1000 --homozyg-window-het 10 --homozyg-window-threshold 0.001  --homozyg-kb 100 --out judaei.ROH
#输出文件judaei.ROH.hom

awk '{print "Chr"$4"\t"$9"\t"$2}' judaei.phased.hom > judaei.tmp
#judaei.tmp就是后续要用到的文件

2. 统计ROH

python stat.ROH.py -i1 judaei.tmp -i2 Chr.length -o1 judaei.tmp.FROH.boxplot.txt -o2 judaei.tmp.FROH.heatmap.txt
#Chr.length这个文件即参考基因组，每条染色体的长度

#第一个输出文件格式要是这样的（正好用于ggplot2画箱线图，比较不同染色体的FROH）：
#第一列是染色体号，第二列是FROH，第三列是样本名

#第二个输出文件是这样的，第一列是染色体号，第二列、第三列...是样本名，
#然后整个表格中的值就是不同样本在不同染色体上的染色体号。
#这个文件正好可以用于R里面的pheatmap直接画图。

cat stat.ROH.py

import argparse
import pandas as pd

def calculate_froh(input_file_a, input_file_b, output_file1, output_file2):
    # 读取染色体长度文件
    chr_length_df = pd.read_csv(input_file_b, sep="\t", header=None, names=["Chr", "Length"])
    chr_length_dict = dict(zip(chr_length_df["Chr"], chr_length_df["Length"]))
    
    # 读取ROH文件
    roh_df = pd.read_csv(input_file_a, sep="\t", header=None, names=["Chr", "Length", "Sample"])
    
    # 计算FROH
    froh_data = roh_df.groupby(["Chr", "Sample"]).sum().reset_index()
    froh_data["FROH"] = froh_data.apply(lambda row: row["Length"] * 1000 / chr_length_dict.get(row["Chr"], 1), axis=1)
    froh_data = froh_data[["Chr", "FROH", "Sample"]]
    
    # 按照染色体顺序排序
    froh_data.sort_values(by="Chr", key=lambda x: x.str.extract(r'(\d+)')[0].astype(int), inplace=True)
    froh_data.to_csv(output_file1, sep="\t", index=False, header=False)
    
    # 生成heatmap格式
    froh_pivot = froh_data.pivot(index="Chr", columns="Sample", values="FROH").fillna(0)
    froh_pivot = froh_pivot.reindex(sorted(chr_length_dict.keys(), key=lambda x: int(x[3:])))
    froh_pivot.to_csv(output_file2, sep="\t")
    
if __name__ == "__main__":
    parser = argparse.ArgumentParser()
    parser.add_argument("-i1", required=True, help="Input file A (ROH file)")
    parser.add_argument("-i2", required=True, help="Input file B (Chromosome length file)")
    parser.add_argument("-o1", required=True, help="Output file 1 (for boxplot)")
    parser.add_argument("-o2", required=True, help="Output file 2 (for heatmap)")
    
    args = parser.parse_args()
    calculate_froh(args.i1, args.i2, args.o1, args.o2)