Significant TMRCA discrepancy between tsdate and Relate on 1KGP YRI (Chr 22)

[Jesson-mark]

Hi,

I am currently using tsinfer + tsdate to reconstruct ARGs and estimate coalescence times. I recently ran this pipeline on Chromosome 22 of the Yoruba (YRI) samples from the 1000 Genomes Project and compared the results with those inferred by Relate.

I noticed a substantial discrepancy in the TMRCA estimates for local trees between the two methods. Specifically, tsdate produces significantly older estimates compared to Relate, especially in the deep past.

Observations

Here is the distribution of TMRCAs across Chromosome 22:

• tsdate:

  • Median TMRCA: ~2.57 million years (assuming 25 years/gen)
  • Max TMRCA: ~40 million years

• Relate:

  • Median TMRCA: ~1.46 million years (assuming 25 years/gen)
  • Max TMRCA: ~15 million years

As shown, the maximum TMRCA inferred by tsdate reaches ~40 MYA, whereas Relate caps around ~15 MYA. The median age is also nearly double in tsdate.

Parameters & Reproducibility

I used the following parameters:

tsinfer:

• mismatch_ratio: 1
• recombination_rate: 1.25e-8

tsdate:

• mutation_rate: 1.25e-8
• Ne (Effective Population Size): I didn't specify Ne in this pipeline.

Questions:

Is such a large discrepancy expected for African populations which have deep coalescence times?

Could this be related to the prior $N_e$ used in tsdate? If I haven't specified a demographic history, would the default prior cause this overestimation in deep time?

Any insights or suggestions on how to align these estimates would be greatly appreciated.

Thanks!

[hyanwong]

Out of interest, what if you don't specify the recombination rate in tsinfer (or set the mismatch ratio to 0, which is the same thing?)

What is the correlation between log(tsdate) and log(tsinfer) time in both cases?

[Jesson-mark]

Hi Yan,

Thanks for the reply!

Following your suggestion, I set both the mismatch_ratio and recombination_rate to None in tsinfer and reran the pipeline. However, the results remain largely unchanged from the initial run (where the mismatch ratio was 1). The local trees are still showing very deep coalescence times, with a median TMRCA of approximately 2.60 million years (assuming 25 years/gen):

I also compared the tree heights between tsinfer and tsdate. I noticed that the raw tree heights produced by tsinfer are almost exclusively clustered around 1, as shown below (the height_infer column):

Here is the scatter plot comparing the two:

Is this distribution of tsinfer tree heights (clustering at 1) expected behavior? I wanted to confirm if this looks correct to you.

Best,

Jie

[hyanwong]

Hmm, how are you getting the tree heights? The height of the tree might not be the height of the oldest MRCA: the oldest node could be older than that. E.g. in the tree below, the oldest node is a little greater than 1, but the oldest MRCA is at freq 0.67

If you want the MRCA times then you need to remove nodes above the local roots, e.g. by calling simplify() first.

The tsinfer times for nodes or mutations should simply be the frequency of each of the alleles used for inference. What happens if you compare the histogram of tsinfer node times with e.g. the histogram of allele frequencies?

from matplotlib import pyplot as plt
# Just the nodes
plt.hist(inferred_ts.nodes_time)

# The mutation freqs
tables = inferred_ts.dump_tables()
tables.compute_mutation_times()
plt.hist(tables.tree_sequence().mutations_time, label="mutation times (equivalent to freqs)")

# Could also use the AFS
plt.plot(np.arange(inferred_ts.num_samples + 1) / inferred_ts.num_samples, inferred_ts.allele_frequency_spectrum(), label="AFS")

# assuming you are using the variant_data format:
freqs = [np.mean(v.genotypes != a) for v, a in zip(vdata.variants(), vdata.sites_ancestral_allele)]
plt.hist(freqs, label="Derived freqs")

I assume that you have called tsdate.preprocess_ts(...) on the output from tsinfer?

[hyanwong]

BTW, I think it's standard to use 27 or 29 years per generation in humans nowadays?

[hyanwong]

(also, could you show some log-time plots of tsdate vs tsinfer ages for the mutations?)

[Jesson-mark]

Hi Yan,

here is the codes that I used to extract tree height:

def create_tree_height_df(tree_seq):
    data = []
    
    for tree in tree_seq.trees():
        if tree.num_roots > 0:
            height = max(tree.time(root) for root in tree.roots)
        else:
            height = 0
            
        data.append({
            'tree_index': tree.index, 
            'tree_height': height
        })

    tree_height_df = pd.DataFrame(data)
    
    return tree_height_df

I didn't call ts.simplify() in previous version. However, when I add ts.simplify() and then extracted the tree height, the plot is consistent with the previous one.

I also plotted the frequency using the codes modified from the codes you provided:

def plot_diagnostics(trees_path, freqs_path, output_path):
    print(f"Loading tree sequence: {trees_path}")
    ts = tskit.load(trees_path)
    
    print(f"Loading pre-calculated frequencies: {freqs_path}")
    vdata_freqs = np.load(freqs_path)
    
    # 创建画布
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    fig.suptitle(f"tsinfer Diagnostics\n(File: {trees_path})", fontsize=16)
    
    # --- 1. Node Times (节点时间分布) ---
    ax = axes[0, 0]
    # 过滤掉时间为0的节点（通常是样本），只看内部节点
    internal_nodes = ts.nodes_time[ts.nodes_time > 0]
    ax.hist(internal_nodes, bins=50, color='skyblue', edgecolor='black', alpha=0.7)
    ax.set_title("1. Node Times (inferred_ts.nodes_time)")
    ax.set_xlabel("Time (Proxy / Frequency)")
    ax.set_ylabel("Count")
    
    # --- 2. Mutation Times (突变时间分布) ---
    # 作者代码: tables.compute_mutation_times() -> mutations_time
    ax = axes[0, 1]
    tables = ts.dump_tables()
    # compute_mutation_times 会根据突变所在的节点时间来估算突变时间
    tables.compute_mutation_times() 
    mut_times = tables.mutations.time
    # 过滤掉 NaN (如果有)
    mut_times = mut_times[~np.isnan(mut_times)]
    
    ax.hist(mut_times, bins=50, color='lightgreen', edgecolor='black', alpha=0.7, label="Mutation Times")
    ax.set_title("2. Mutation Times (derived from node times)")
    ax.set_xlabel("Time (Proxy)")
    ax.legend()

    # --- 3. Allele Frequency Spectrum (AFS) ---
    # 作者代码: plt.plot(..., ts.allele_frequency_spectrum())
    ax = axes[1, 0]
    afs = ts.allele_frequency_spectrum()
    # 归一化 x 轴：从 0 到 1
    x_axis = np.arange(len(afs)) / (len(afs) - 1)
    
    # 通常 AFS 的首尾 (0 和 1) 计数极大，为了看清中间趋势，我们切掉首尾
    # 或者使用 log scale
    ax.plot(x_axis[1:-1], afs[1:-1], color='purple', marker='.', linestyle='-', linewidth=1, label="AFS")
    ax.set_title("3. Allele Frequency Spectrum (AFS)")
    ax.set_xlabel("Frequency")
    ax.set_ylabel("Count (Log Scale)")
    ax.set_yscale('log') # AFS 通常跨度很大，用对数坐标更清晰
    ax.legend()

    # --- 4. Derived Allele Frequencies (原始数据频率) ---
    # 作者代码: plt.hist(freqs, label="Derived freqs")
    # 这里使用你从服务器下载的 .npy 文件
    ax = axes[1, 1]
    # 清洗数据：去除 NaN
    clean_freqs = vdata_freqs[~np.isnan(vdata_freqs)]
    
    ax.hist(clean_freqs, bins=50, color='salmon', edgecolor='black', alpha=0.7, label="Derived Freqs (vdata)")
    ax.set_title("4. Raw Derived Allele Frequencies (vdata)")
    ax.set_xlabel("Frequency")
    ax.set_ylabel("Count")
    ax.legend()

    plt.tight_layout(rect=[0, 0.03, 1, 0.95]) # 调整布局以容纳总标题
    plt.savefig(output_path, dpi=300)
    print(f"Diagnostic plot saved to: {output_path}")
    plt.show()

the plot:

Besides, I also plotted the log-time of tsdate vs tsinfer ages for the mutations, with codes below:

def plot_mutation_comparison(tsinfer_path, tsdate_path, output_file="mutation_time_comparison.png"):
    # 1. 加载两个 TreeSequence
    print(f"Loading tsinfer: {tsinfer_path}")
    ts_infer = tskit.load(tsinfer_path)
    
    print(f"Loading tsdate: {tsdate_path}")
    ts_date = tskit.load(tsdate_path)
    
    # 检查突变数量是否一致（通常必须一致）
    if ts_infer.num_mutations != ts_date.num_mutations:
        print(f"Warning: Mutation counts differ! Infer: {ts_infer.num_mutations}, Date: {ts_date.num_mutations}")
        # 如果不一致，可能需要取交集，但通常流程下是一致的
        return

    # 2. 计算突变时间 (Compute Mutation Times)
    # tsinfer 和 tsdate 的原始 node time 并不直接等于 mutation time
    # tskit 提供了一个方法，根据突变所在的枝长位置估算其精确时间
    
    print("Computing mutation times for tsinfer...")
    tables_infer = ts_infer.dump_tables()
    tables_infer.compute_mutation_times()
    # 这里的 time 其实就是频率 (Frequency) 或 代理时间
    infer_mut_times = tables_infer.mutations.time
    
    print("Computing mutation times for tsdate...")
    tables_date = ts_date.dump_tables()
    tables_date.compute_mutation_times()
    # 这里的 time 是世代数 (Generations)
    date_mut_times = tables_date.mutations.time
    
    # 3. 数据清洗
    # 因为要画 Log 图，必须去掉 <= 0 的值
    # tsinfer 的时间通常是频率 (0~1)，tsdate 是世代 (0~inf)
    mask = (infer_mut_times > 0) & (date_mut_times > 0)
    
    x = infer_mut_times[mask]
    y = date_mut_times[mask]
    
    print(f"Plotting {len(x)} mutations (filtered {len(infer_mut_times) - len(x)} <= 0)...")

    # 4. 绘图
    plt.figure(figsize=(10, 8))
    
    # 使用 hexbin (六边形分箱图) 而不是 scatter
    # 因为突变通常有几十万个，scatter 会糊成一团，hexbin 能看到密度分布
    hb = plt.hexbin(x, y, gridsize=50, cmap='inferno', xscale='log', yscale='log', mincnt=1)
    
    # 添加颜色条
    cb = plt.colorbar(hb, label='Count of Mutations')
    
    # 添加对角线参考（虽然单位不同，但斜率应该大致为正）
    # 这里不画 y=x，因为单位完全不同
    
    plt.title("Mutation Ages Comparison: tsinfer vs tsdate", fontsize=15)
    plt.xlabel("tsinfer Mutation Time (Frequency / Proxy)\n(Log Scale)", fontsize=12)
    plt.ylabel("tsdate Mutation Time (Generations)\n(Log Scale)", fontsize=12)
    
    plt.grid(True, which="both", ls="--", alpha=0.3)
    
    plt.tight_layout()
    plt.savefig(output_file, dpi=300)
    print(f"Plot saved to {output_file}")
    plt.show()

the plot:

It seems the times between tsinfer and tsdate are well correlated.

Besides, I used 25 years per generation for convenience, although it's more common to use 29 years per generation.

Pipeline of tsinfer + tsdate

I do called tsdate.preprocess_ts(...) on the output from tsinfer in tsdate.

Here is the pipeline that I used to run tsinfer:

def main():
    args = get_args()

    vcz_path = args.vcz_path
    recombination_rate = args.recombination_rate
    mismatch_ratio = args.mismatch_ratio
    threads = args.threads
    out_file = args.out_file # with .trees suffix or .tree.tsz

    check_files(vcz_path)

    print(f'vcf_path is {vcz_path}')
    print(f'recombination_rate is {recombination_rate}')
    print(f'mismatch_ratio is {mismatch_ratio}')
    print(f'threads is {threads}')
    print(f'out_file is {out_file}')

    print(f'Loading vcz file: {vcz_path}')
    vdata = tsinfer.VariantData(vcz_path, ancestral_state="ancestral_state", individuals_population="individuals_population")

    if recombination_rate is None:
        mismatch_ratio = None
        print("Recombination rate is None, so mismatch ratio will be set to None as well.")

    print(f'Running inference with mismatch ratio: {mismatch_ratio}, recombination rate: {recombination_rate}, threads: {threads}')
    inferred_ts = tsinfer.infer(
        vdata,
        recombination_rate=recombination_rate,
        mismatch_ratio=mismatch_ratio,
        num_threads=threads
    )

    print(f'Saving inferred tree sequence to file: {out_file}')
    if out_file.endswith('.tsz'):
        tszip.compress(inferred_ts, out_file)
    else:
        inferred_ts.dump(out_file)

    print(f'Output file is {out_file}')

    print("Program is done")

and tsdate:

def main():
    args = get_args()

    ts_file = args.ts_file
    mutation_rate = args.mutation_rate
    out_file = args.out_file

    check_files(ts_file)
    print(f'ts file: {ts_file}')
    print(f'mutation rate: {mutation_rate}')
    print(f'output file: {out_file}')

    print(f"Load tree sequence from {ts_file}")
    ts = tszip.load(ts_file)
    print(f"Original: {ts.num_nodes} nodes, {ts.num_mutations} mutations")

    print("Simplifying the tree sequence")
    ts_clean = ts.simplify(filter_populations=False)
    print(f"Simplified: {ts_clean.num_nodes} nodes, {ts_clean.num_mutations} mutations")

    ## preprocess ts
    print("Preprocessing the tree sequence")
    simplified_ts = tsdate.preprocess_ts(
        ts_clean,
        remove_telomeres=True
    )

    ## infer the time
    print("Infering the time")
    dated_ts = tsdate.date(simplified_ts,
                           mutation_rate=mutation_rate)

    print("Restoring population metadata...")
    tables = dated_ts.dump_tables()

    # 检查 ts_clean 中是否有种群信息
    if len(ts_clean.tables.populations) > 0:
        # 直接用 ts_clean 的种群表替换 dated_ts 的种群表
        # 因为 simplify(filter_populations=False) 保证了索引没乱
        tables.populations.replace_with(ts_clean.tables.populations)

        # since the metadata_schema is JSON in the tree outputed by tsinfer, we need to set it to None here which will use bytes string
        tables.populations.metadata_schema = tskit.MetadataSchema(None)

        # 重新构建最终的 TreeSequence
        final_ts = tables.tree_sequence()
        print(f"Restored metadata for {len(final_ts.tables.populations)} populations.")
    else:
        print("Warning: No population metadata found in input.")
        final_ts = dated_ts

    if out_file.endswith(".tsz"):
        print(f"Save the dated tree sequence to {out_file}")
        tszip.compress(final_ts, out_file)
    else:
        print(f"Save the dated tree sequence to {out_file}")
        final_ts.dump(out_file)

    print("Program is done")
    print(f'Dated tree sequence saved to {out_file}')