思路：
1.拟配对个体ROH片段取交集
2.提取出交集的染色体编号和位置坐标
3.从总vcf文件中提取出拟配对个体的vcf文件
4.根据第2步中的染色体编号和位置坐标提取出每个ROH片段对应的vcf文件
5.制定规则，预测拟配对个体配对后产生的新的vcf文件
6.新的vcf文件能否继续被评估为ROH片段，预测ROH的保留和消除
7.计算每一个配对ROH保留和消除的概率
代码实现：
1.配对个体ROH片段取交集
import pandas as pd
# 定义读取数据的函数
def read_data(file_path):
return pd.read_csv(file_path, sep='\t', usecols=[3, 6, 7], header=0, names=['CHR', 'POS1', 'POS2'])
# 定义寻找相交ROH的函数
def find_intersections(df1, df2):
intersections = []
for chr_num in pd.unique(df1['CHR']):
df1_chr = df1[df1['CHR'] == chr_num]
df2_chr = df2[df2['CHR'] == chr_num]
for _, row1 in df1_chr.iterrows():
for _, row2 in df2_chr.iterrows():
start_max = max(row1['POS1'], row2['POS1'])
end_min = min(row1['POS2'], row2['POS2'])
if start_max < end_min:
kb = end_min - start_max
intersections.append((chr_num, start_max, end_min, kb))
return intersections
# 加载数据
data_1 = read_data('467.txt') # 请替换为第一个文件的实际路径
data_2 = read_data('522.txt') # 请替换为第二个文件的实际路径
# 查找相交区域
intersections = find_intersections(data_1, data_2)
# 将结果转换为DataFrame
intersections_df = pd.DataFrame(intersections, columns=['CHR', 'POS1', 'POS2', 'KB'])
# 显示或保存结果
print(intersections_df)
intersections_df.to_csv('intersections.txt', sep='\t', index=False)
2.提取出交集的染色体编号和位置坐标
awk '{print $1 "\t" $2"-"$3}' intersections.txt > 1.txt
3.从总vcf文件中提取出拟配对个体的vcf文件
## 1.提取出配对个体的snp位点组成新的vcf文件
vcftools --vcf chr1.vcf --recode --recode-INFO-all --stdout --indv Pta_467 --indv Pta_522 > 467_522.vcf
## 2.vcf文件转换为vcf.gz
bcftools view 467_522.vcf -Oz -o 467_522.vcf.gz
## 3.加标签
bcftools index 467_522.vcf.gz
## 4.使用bcftools提取某一范围内的snp位点
bcftools filter 467_522.vcf.gz --regions 1:75520805-81296155 > position1.vcf
4. 根据第2步中的染色体编号和位置坐标提取出每个ROH片段对应的vcf文件 (利用脚本实现批量化处理）
这个脚本会占用服务器全部线程，这可能会导致其他任务崩溃，所以我在下面更新了代码
#!/bin/bash
# 假定原始 VCF 文件名为 source.vcf.gz
input_vcf="467_522.vcf.gz"
# 读取1.txt文件中的每一行
while IFS=$'\t' read -r chr range; do
if [ "$chr" != "CHR" ]; then # 跳过标题行
# 生成输出文件名
output_file="CHR${chr}.${range}.vcf"
# 格式化区域字符串
region="${chr}:${range}"
# 使用 bcftools 提取指定区域的数据
bcftools filter "$input_vcf" --regions "$region" > "$output_file"
fi
done < "1.txt"
4.更新后提取出每个ROH片段对应的vcf文件（利用snakemake）
自动设置使用线程，不会导致服务器崩溃
import pandas as pd
# 载入位置信息
positions = pd.read_csv("chr_position.txt", sep="\t", header=0)
# 生成明确的输出文件列表
output_files = ["CHR{0}.{1}.vcf".format(row['CHR'], row['POS1-POS2']) for index, row in positions.iterrows()]
rule all:
input:
output_files
rule filter_vcf:
input:
vcf="467_334.vcf.gz"
output:
vcf="CHR{chr}.{range}.vcf"
params:
region=lambda wildcards: f"{wildcards.chr}:{wildcards.range}" # Simplified region specification
shell:
"""
bcftools filter {input.vcf} --regions {params.region} > {output.vcf}
"""
5. 制定规则，预测拟配对个体配对后产生的新的基因型的vcf文件
这一步需要注意，我将建立一个新的文件夹，将生成的结果导入到新的文件夹中，同时生成chr.txt文件，这个文件是上一步生成的vcf文件的路径和文件名
ls CHR* | tr "\t" "\n" > CHR.txt
awk '{print "../" $1}' CHR.txt > chr.txt
生成新基因型vcf文件的代码
def recombine_genotypes(gt1, gt2):
# 根据两个基因型来重组新的基因型
if gt1 == "0/1" or gt2 == "0/1":
return "0/1"
elif (gt1 == "0/0" and gt2 == "1/1") or (gt1 == "1/1" and gt2 == "0/0"):
return "0/1"
elif gt1 == "1/1" and gt2 == "1/1":
return "1/1"
else:
return "0/0"
def process_vcf(input_filename, output_filename):
with open(input_filename, 'r') as file_in, open(output_filename, 'w') as file_out:
for line in file_in:
if line.startswith('#'):
if line.startswith('#CHROM'):
# 修改列标题
headers = line.strip().split('\t')
headers = headers[:9] + ['forecast']
file_out.write('\t'.join(headers) + '\n')
else:
# 复制注释行
file_out.write(line)
else:
# 处理基因型数据
parts = line.strip().split('\t')
gt1 = parts[9]
gt2 = parts[10]
forecast_gt = recombine_genotypes(gt1, gt2)
new_line = parts[:9] + [forecast_gt]
file_out.write('\t'.join(new_line) + '\n')
def batch_process(file_list):
with open(file_list, 'r') as files:
for file_path in files:
file_path = file_path.strip()
output_filename = file_path.split('/')[-1] # 提取文件名
process_vcf(file_path, output_filename)
# 输入文件名列表文件
file_list = 'chr.txt'
# 调用批量处理函数
batch_process(file_list)
6.新的vcf文件能否继续被评估为ROH片段，预测ROH的保留和消除
### 新建文件夹ped
ls *vcf | tr "\t" "\n" | awk -F "." '{print $1"."$2}' | while read id ; do plink --vcf $id.vcf --make-bed --out ped/$id ; done
ls *bed | tr "\t" "\n" | awk -F "." '{print $1"."$2}' | while read id ; do plink --bfile $id --recode --out $id ; done
### 新建文件夹ROH
ls *bed | tr "\t" "\n" | awk -F "." '{print $1"."$2}' | while read id ; do plink --bfile $id --homozyg --homozyg-density 20 --homozyg-gap 100 --homozyg-kb 500 --homozyg-snp 50 --homozyg-window-het 3 --homozyg-window-snp 50 --homozyg-window-threshold 0.05 --homozyg-window-missing 5 --out ROH/ROH$id; done
rm ROH*summary ROH*log ROH*nosex ROH*indiv
7.计算每一个配对ROH保留和消除的概率
find . -type f -size +135c | wc -l