0.脚本由大模型生成
- 鉴定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就是后续要用到的文件
- 统计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)
