脚本由GPT生成
运行完nucmer后，想要知道每条染色体的比对情况
nucmer --mum -t 6 -b 500 ref.fa query.fa -p ref.query
delta-filter -1 -l 500 ref.query.delta > ref.query.filtered.delta
show-coords -H -c -l -o -r -T ref.query.filtered.delta > ref.query.filtered.coords
#ref.query.filtered.coords 就是输入文件
cat merge.gap.analyze_alignment.py
#统计nucmer的比对结果
#nucmer --mum -t 6 -b 500 ref.fa query.fa -p ref.query
#delta-filter -1 -l 500 ref.query.delta > ref.query.filtered.delta
#show-coords -H -c -l -o -r -T ref.query.filtered.delta > ref.query.filtered.coords
#ref.query.filtered.coords 就是本脚本的输入文件
#Usage: python merge.gap.analyze_alignment.py $i.filtered.coords -rag=200000 > $i.200kgap.stat.txt
#-rag=的含义是如果相邻的比对片段在两个基因组上的距离都小于 rag，则合并
#比如-rag=30000的时候
#704513 705497 72524266 72525272 --> 初始比对片段
#721485 723852 72549371 72551732 --> 15988 和 25105 都 < 30000，所以合并
#
#合并后
#704513 723852 72524266 72551732
#Ref_Chr ref基因组 染色体
#Qry_Chr query基因组 染色体
#Total_Aligned_Length 合并后的总比对长度
#Ref_Coverage(%) ref基因组 染色体被 query基因组 覆盖的比例
#Qry_Coverage(%) query基因组 染色体被 ref基因组 覆盖的比例
#
import sys
from collections import defaultdict
def parse_coords(file_path, rag_threshold):
""" 解析coords文件，并根据rag阈值合并片段 """
chr_stats = defaultdict(lambda: defaultdict(list)) # 按 (ref_chr, qry_chr) 存储片段
with open(file_path, "r") as f:
for line in f:
cols = line.strip().split("\t")
if len(cols) < 13:
continue # 确保数据完整
# 提取必要信息
ref_chr, qry_chr = cols[11], cols[12] # 染色体名称
ref_start, ref_end = int(cols[0]), int(cols[1]) # 参考基因组位置
qry_start, qry_end = int(cols[2]), int(cols[3]) # 查询基因组位置
ref_length, qry_length = int(cols[7]), int(cols[8]) # 染色体总长
chr_stats[ref_chr][qry_chr].append((ref_start, ref_end, qry_start, qry_end, ref_length, qry_length))
# 合并片段
merged_stats = defaultdict(lambda: defaultdict(list))
for ref_chr, qry_dict in chr_stats.items():
for qry_chr, segments in qry_dict.items():
segments.sort() # 按 ref_start 排序
merged = []
current_start, current_end, qry_start, qry_end, ref_len, qry_len = segments[0]
for i in range(1, len(segments)):
ref_s, ref_e, qry_s, qry_e, _, _ = segments[i]
ref_gap = ref_s - current_end
qry_gap = qry_s - qry_end
if ref_gap <= rag_threshold and qry_gap <= rag_threshold:
# 合并
current_end = max(current_end, ref_e)
qry_end = max(qry_end, qry_e)
else:
merged.append((current_start, current_end, qry_start, qry_end, ref_len, qry_len))
current_start, current_end, qry_start, qry_end = ref_s, ref_e, qry_s, qry_e
# 添加最后一个片段
merged.append((current_start, current_end, qry_start, qry_end, ref_len, qry_len))
merged_stats[ref_chr][qry_chr] = merged
return merged_stats
def compute_coverage(merged_stats):
""" 计算合并后的覆盖度和比对长度 """
results = []
for ref_chr, qry_dict in merged_stats.items():
for qry_chr, segments in qry_dict.items():
total_aligned_length = sum(seg[1] - seg[0] + 1 for seg in segments)
ref_length = segments[0][4] if segments else 1 # 参考染色体长度
qry_length = segments[0][5] if segments else 1 # 查询染色体长度
coverage_ref = (total_aligned_length / ref_length) * 100 if ref_length > 0 else 0
coverage_qry = (total_aligned_length / qry_length) * 100 if qry_length > 0 else 0
results.append((ref_chr, qry_chr, total_aligned_length, coverage_ref, coverage_qry))
return results
def main():
# 检查命令行参数
if len(sys.argv) != 3:
print("Usage: python merge_and_analyze.py <filtered.coords> -rag <threshold>")
sys.exit(1)
coords_file = sys.argv[1]
rag_threshold = int(sys.argv[2].replace("-rag=", ""))
# 解析并合并
merged_stats = parse_coords(coords_file, rag_threshold)
results = compute_coverage(merged_stats)
# 输出结果（支持 `>` 重定向）
print("Ref_Chr\tQry_Chr\tTotal_Aligned_Length\tRef_Coverage(%)\tQry_Coverage(%)")
for ref_chr, qry_chr, aligned_length, coverage_ref, coverage_qry in results:
print(f"{ref_chr}\t{qry_chr}\t{aligned_length}\t{coverage_ref:.2f}\t{coverage_qry:.2f}")
if __name__ == "__main__":
main()
cat analyze_alignment.py
#统计nucmer的比对结果
#nucmer --mum -t 6 -b 500 ref.fa query.fa -p ref.query
#delta-filter -1 -l 500 ref.query.delta > ref.query.filtered.delta
#show-coords -H -c -l -o -r -T ref.query.filtered.delta > ref.query.filtered.coords
#ref.query.filtered.coords 就是本脚本的输入文件
#Usage: python analyze_alignment.py $i.filtered.coords > $i.200kgap.stat.txt
#输出解释
#Ref_Chr ref基因组染色体名称
#Qry_Chr query基因组染色体名称
#Aligned_Regions 比对区段数量
#Total_Aligned_Length 参考染色体和目标染色体的比对总长度（bp）
#Ref_Coverage(%) ref基因组 该染色体比对到 query基因组 的比例 (aligned_length / ref_length * 100)
#Qry_Coverage(%) query基因组 该染色体被 ref基因组 覆盖的比例 (aligned_length / qry_length * 100)
#Average_Identity 平均序列相似度
import sys
from collections import defaultdict
# 检查是否传入了文件路径
if len(sys.argv) != 2:
print("Usage: python analyze_alignment.py <filtered.coords>")
sys.exit(1)
coords_file = sys.argv[1] # 从命令行参数获取输入文件
# 统计数据结构
chr_stats = defaultdict(lambda: defaultdict(lambda: {"count": 0, "aligned_length": 0, "total_identity": 0, "ref_length": 0, "qry_length": 0}))
# 读取 filtered.coords 文件并解析数据
with open(coords_file, "r") as f:
for line in f:
cols = line.strip().split("\t")
if len(cols) < 13:
continue # 确保数据完整
ref_chr = cols[11] # Chicken.T2T 染色体
qry_chr = cols[12] # Falco_naumanni 染色体
length = int(cols[5]) # 比对长度
identity = float(cols[6]) # 相似度
ref_length = int(cols[7]) # Chicken.T2T 染色体总长度
qry_length = int(cols[8]) # Falco_naumanni 染色体总长度
# 更新参考基因组（Chicken.T2T）比对情况
chr_stats[ref_chr][qry_chr]["count"] += 1
chr_stats[ref_chr][qry_chr]["aligned_length"] += length
chr_stats[ref_chr][qry_chr]["total_identity"] += identity * length # 加权相似度
chr_stats[ref_chr][qry_chr]["ref_length"] = ref_length # 参考染色体长度
chr_stats[ref_chr][qry_chr]["qry_length"] = qry_length # 目标染色体长度
# 打印结果（支持重定向到文件）
print("Ref_Chr\tQry_Chr\tAligned_Regions\tTotal_Aligned_Length\tRef_Coverage(%)\tQry_Coverage(%)\tAverage_Identity")
for ref_chr, qry_dict in chr_stats.items():
for qry_chr, stats in qry_dict.items():
avg_identity = stats["total_identity"] / stats["aligned_length"] if stats["aligned_length"] > 0 else 0
coverage_ref = (stats["aligned_length"] / stats["ref_length"]) * 100 if stats["ref_length"] > 0 else 0
coverage_qry = (stats["aligned_length"] / stats["qry_length"]) * 100 if stats["qry_length"] > 0 else 0
print(f"{ref_chr}\t{qry_chr}\t{stats['count']}\t{stats['aligned_length']}\t{coverage_ref:.2f}\t{coverage_qry:.2f}\t{avg_identity:.2f}")
