假设

新版本的实验设计:

1.杂合易位系由于两个原因，形成了跨染色体的ld

2. 裸子植物中的杂合易位系很有可能造成了非断点层面的重组的抑制。

从日本文献和浙大的文献来看，裸子植物中，杂合易位系形成的四价体是一个松散的四价体，很快会打碎形成不同的二价体，在这种情况下，重组会抑制，减弱，LD会增加。但是另外一方面，我们仍然应当能在跨染色体尺度上看到较强信号的LD。两者是不冲突的。

实操

跨染色体的ld计算

抽样

全局的SNP总数太大了，不可能能计算完，所以必须要用抽样的办法。

#!/bin/bash
set -euo pipefail

vcf="11_27_merged_vcf.gz"
out="merged.10000snps.vcf"
n_snp=10000
threads=24   # 并行线程数，可根据CPU调整

# 获取所有染色体列表
chroms=$(bcftools query -f '%CHROM\n' "$vcf" | sort -u)

# 并行抽样每条染色体
parallel -j $threads '
  chr={} 
  bcftools view -r "$chr" -v snps '"$vcf"' \
  | bcftools view -H \
  | shuf -n '"$n_snp"' \
  > tmp_chr_"$chr".vcf
' ::: $chroms

# 生成最终合并文件
bcftools view -h "$vcf" > "$out"   # 先输出头部
for chr in $chroms; do
    cat tmp_chr_"$chr".vcf >> "$out"
done

# 清理临时文件
rm tmp_chr_*.vcf

echo "抽样完成，结果保存在 $out"ww

VCF → bfile

plink \
  --vcf merged.10000snps.vcf \
  --make-bed \
  --double-id \
  --out merged.10000snps

计算

plink --bfile merged_1000snps       --r2 inter-chr       --ld-window-r2 0       --out ld_all_inter_full

注：ld_all_inter_full这个文件有359GB大小

统计

awk '
BEGIN {
    # 设置输出格式为制表符分隔
    OFS = "\t";
    print "ChrA", "ChrB", "Total_Pairs", "Avg_R2"
}
# 跳过标题行 (NR > 1)
NR > 1 {
    # 提取关键字段 (CHR_A在第1列，CHR_B在第4列，R2在第7列)
    chrA = $1
    chrB = $4
    r2 = $7
    
    # 创建标准化染色体对标识
    key = (chrA <= chrB) ? chrA ":" chrB : chrB ":" chrA
    
    # 累加统计值
    sum_r2[key] += r2
    count[key]++
}
END {
    # 输出各染色体对结果
    for (pair in sum_r2) {
        split(pair, chrs, ":")
        avg = sum_r2[pair] / count[pair]
        print chrs[1], chrs[2], count[pair], avg
    }
    
    # 添加全局平均值行
    total_pairs = 0
    total_r2 = 0
    for (pair in sum_r2) {
        total_pairs += count[pair]
        total_r2 += sum_r2[pair]
    }
    if (total_pairs > 0) {
        global_avg = total_r2 / total_pairs
        print "Global", "Average", total_pairs, global_avg
    }
}' ld_all_inter_full.ld > chromosome_pair_avg_r2.tsv

像这样的输出结果应该可以绘制一个三角热图。就是有点慢了，总的数据量还是很大的。

相当显著

染色体内的LD计算

 plink --vcf 11_27_merged_vcf.gz --make-bed --out 11_27_merged_vcf --double-id

plink --bfile 11_27_merged_vcf --r2 --ld-window 99999 --ld-window-kb 1000 --ld-window-r2 0

LD衰减曲线的计算

这个时候应该要计算不同染色体的LD衰减曲线，工具考虑使用这个https://github.com/BGI-shenzhen/PopLDdecay文章原文doi.org/10.1093/bioinformatics/bty875。

文章原文与文件储存在老猫个人计算机的裸子植物染色体臂易位固定机制\何伟明ld计算工具路径下。