仅作个人笔记 使用4DNvestigator计算VNE

% ==== 固定参数（按需改动路径） ====
repo      = 'F:\Matlab.workspace\4DNvestigator';
juicerJar = fullfile(repo, 'juicer_tools.jar');        % 建议 Juicer Tools 2.20.00
hicFile   = 'E:\workspace\FenshuChrVar\48.heart.HiC\Erothschildi.heart.hic';
outDir    = 'E:\workspace\FenshuChrVar\48.heart.HiC\VNE_allchr_fixed';
bpFrag    = 'BP';
normName  = 'KR';
binSize   = 1e5;                                       % 固定分辨率：100 kb

% ==== 环境 ====
cd(repo); addpath(genpath(repo));
if ~exist(outDir,'dir'), mkdir(outDir); end

% ==== 主流程：Chr1..Chr31，整条染色体一次性计算 ====
chroms = 1:29;
VNE    = nan(numel(chroms),1);

for idx = 1:numel(chroms)
    c = chroms(idx);
    chrArg = sprintf('Chr%d', c);                      % 你的 .hic 命名为 ChrN

    dumpTxt = fullfile(outDir, sprintf('dump_%s_%dkb.txt', chrArg, binSize/1000));
    cmd = sprintf('java -jar "%s" dump oe %s "%s" %s %s %s %d "%s"', ...
        juicerJar, upper(normName), hicFile, chrArg, chrArg, upper(bpFrag), binSize, dumpTxt);
    st = system(cmd);
    if st~=0 || ~exist(dumpTxt,'file') || dir(dumpTxt).bytes==0
        % 若极少数样本并非 ChrN 命名，可试备用命名（chrN）
        chrArg2 = sprintf('chr%d', c);
        dumpTxt = fullfile(outDir, sprintf('dump_%s_%dkb.txt', chrArg2, binSize/1000));
        cmd = sprintf('java -jar "%s" dump oe %s "%s" %s %s %s %d "%s"', ...
            juicerJar, upper(normName), hicFile, chrArg2, chrArg2, upper(bpFrag), binSize, dumpTxt);
        st = system(cmd);
        if st~=0 || ~exist(dumpTxt,'file') || dir(dumpTxt).bytes==0
            warning('Chr%d 无法导出（命名或分辨率不存在），VNE=NaN', c);
            VNE(idx) = NaN; continue
        end
    end

    M = readmatrix(dumpTxt);
    if isempty(M) || size(M,2)<3
        warning('Chr%d dump 文件异常，VNE=NaN', c);
        VNE(idx) = NaN; continue
    end

    i0 = double(M(:,1)); j0 = double(M(:,2)); v = double(M(:,3));
    if max([i0; j0]) > 1e7
        % dump 输出为坐标，转为 bin 索引
        i = floor(i0./binSize) + 1;
        j = floor(j0./binSize) + 1;
    else
        % dump 输出为 0-based bin 索引，转 1-based
        i = i0 + 1;
        j = j0 + 1;
    end

    nBins = max([i; j]);
    if nBins < 2
        warning('Chr%d nBins<2，VNE=NaN', c);
        VNE(idx) = NaN; continue
    end

    S = sparse(i, j, v, nBins, nBins);
    S = S + S.' - spdiags(diag(S), 0, nBins, nBins);   % 对称化

    % 4DNvestigator 默认：先做相关矩阵，再算 VNE
    H = full(S);                                       % 注意：这里会占用 O(nBins^2) 内存
    VNE(idx) = hicEntropy(H, [], [], 'corr');

    fprintf('Chr%d 完成，%dkb，nBins=%d，VNE=%.4f\n', c, binSize/1000, nBins, VNE(idx));
end

T = table(chroms(:), VNE, 'VariableNames', {'chrom','VNE'});
outTsv = fullfile(outDir, sprintf('VNE_allchr_%dkb.tsv', binSize/1000));
writetable(T, outTsv, 'FileType','text','Delimiter','\t');

disp('完成，结果表：'); disp(outTsv);