根据华南虎的谱系，我得到了每个华南虎个体所属的世代，我将老虎按照各自的世代进行分群，分析随着世代的增加，各项遗传参数的变化
1.计算当前有效种群规模
- 根据Hill and Weir (1988) 的公式，Ne可以通过r²值计算：
image.png
第一步：计算连锁不平衡
plink --bfile /home/DL_2/CUTech/comput_2/Tech_No1/group4/yjc/workspace/SCT_data/var/01.combine_gvcf/04.deep_allele_maf/01.plink/SCT1_SNP_ind_maf0.05_autosome --r2 --ld-window-r2 0 --ld-window 99999 --ld-window-kb 100 --out ld_output
#--ld-window-r2 0 设置无最小r^2阈值。
#--ld-window 99999 表示窗口内最大考虑的SNP对数量。
#--ld-window-kb 100 设置窗口大小为100千碱基。
当然，在进行第一步之前，我需要将群体按照代数进行拆分，将不同的代数的个体单独生成文件
plink --bfile data --keep keep_individuals.txt --make-bed --out extracted_708
#keep_individuals.txt文件格式
HNH 708
### 很明显在生成keep_individuals.txt这个文件，我需要截取fam文件中前两列，并将每一代的个体从fam文件的前两列中提出出来
awk 'NR==FNR{ids[$1]; order[FNR]=$1; next} $1 in ids {data[$1]=$2} END {for (i in order) print order[i], data[order[i]]}' B.txt A.txt > result.txt
image.png
我发现第一步计算的结果文件，id文件中的R2值存在很多极度接近1的值，我想删除文件中R2大于0.8的值，但文件很大，计算很慢，所以我选择了一种更快的方式，将调用更多的线程来完成这项工作
#第一步，拆分文件
split -l 5000000 eight_output.ld part_
#第二步，创建一个 filter.awk 文件，内容如下：
$7 <= 0.8 {print}
#使用 parallel 和 awk 处理所有分割后的文件，并将结果合并到一个新文件中
parallel "awk -f filter.awk {} > {.}_filtered.ld" ::: part_*
#步骤 4: 合并过滤后的文件
cat part_*_filtered.ld > filtered_output.ld
步骤 5: 清理中间文件
rm part_* part_*_filtered.ld
在完成最后一步之前，我要将上一步生成的文件filtered_output.ld的第一行增加列名，脚本如下：
import pandas as pd
# 定义列名
column_names = ['CHR_A', 'BP_A', 'SNP_A', 'CHR_B', 'BP_B', 'SNP_B', 'R2']
# 读取文件，假设字段是由空格分隔的
data = pd.read_csv('te_output.ld', delim_whitespace=True, header=None)
# 设置列名
data.columns = column_names
# 将带有列名的数据写回到一个新的文件
data.to_csv('nin_filtered_output.ld', sep='\t', index=False)
最后一步，利用Python脚本计算Ne，其中公式中c这个值是重组率，引用自一篇文章，这篇文章计算了猫的重组率
Dumont, B.L. & Payseur, B.A. Evolution of the genomic rate of recombination in mammals. Evolution: International Journal of Organic Evolution 62, 276-294 (2008).
import pandas as pd
# 读取PLINK输出的LD结果文件
ld_data = pd.read_csv('nin_filtered_output.ld', delim_whitespace=True)
# 计算平均r^2值
average_r2 = ld_data['R2'].mean()
# 设置标记间的平均物理距离（这里假设为50Kb，需根据实际情况调整）
#1Mb = 1.1 cM
c = 0.011 # 转换成摩尔根
# 计算短期有效种群大小
Ne = (1 / (3 * c * average_r2))
print(f'Estimated short-term effective population size (Ne): {Ne}')
