inferCNV相关思考（配对样本、多样本问题 ）

由于层次聚类的问题（数据量太大的时候，层次聚类的速度很慢，同时会报堆错误，出现内存溢出），inferCNV的不能处理大批量的数据。另外，对于来源于不同病人的样本，inferCNV的issue的讨论也是分开跑。
首先是对于配对的病人和多样本的的处理问题
官网的解释： https://github.com/broadinstitute/infercnv/issues/429
对于同时具癌旁和癌paired的病人，可以使用正常的epithelial cel作为reference，tumor的epi cell作为observation。但是对于同时提供多个Normal作为ref，inferCNV并不能把他们和tumor的epi cell做到一一对应。
另外，单独对一名患者的数据运行 infercnv 确实会更快，并且如果您的不同样本具有相似的类型，则结果应该在很大程度上相似，尽管它们可能不会完全相同，因为低表达/不表达的初始过滤基因会受到细胞总数以及不同细胞类型之间混合比例的影响。由于你对亚克隆结构感兴趣，官网强烈建议使用HMM=T,analysis_mode="subclusters"否则预测可能不准确。
如果您决定对多个患者的数据一起运行 infercnv，除了cluster_by_groups您提到的选项之外，plot_per_group()函数还可以将每个组绘制在不同的图形上，并且dynamic_resize绘图方法的选项可以使热图的可读性更好，得到的heatmap对于大型数据集友好。

对于单样本分析完成样本，如果想做可比的heatmap图，有一些问题：

1. plot_cnv默认的参数x轴的scale，是直接根据测到gene数目进行scale，需要使用plot_chr_scale 参数来保证染色体的宽度是可比的

chr_df = read.chromInfo(species = 'hg38')$df
  chr_df = chr_df[chr_df$chr %in% paste0("chr", 1:22), ]
  chrLen <- chr_df$end
  names(chrLen) <- chr_df$chr
  
  cnvobj <- readRDS(file = paste0('./result/run.final.infercnv_obj'))
  fileDir <- paste0('./result_ChromeScale/')
  dir_mk(fileDir)
  plot_cnv(infercnv_obj = cnvobj,
           out_dir = fileDir,
           title = sampleid,
           plot_chr_scale = T,
           chr_lengths = chrLen,
           #old version.
           x.range = c(0.9,1.1),
           x.center=1,
           png_res = 300,
           output_format = 'png',
           useRaster = T)

2. 每个样本测到的基因数目是不一样的，为了不同样本的CNV可比，需要使用complexHeatmap的提供的帮助函数，从而获得chromesome level的CNV结果。具体的操作如下：

library(circlize)
library(GenomicRanges)
chr_df = read.chromInfo()$df
chr_df = chr_df[chr_df$chr %in% paste0("chr", 1:22), ]
chr_gr = GRanges(seqnames = chr_df[, 1], ranges = IRanges(chr_df[, 2] + 1, chr_df[, 3]))

library(EnrichedHeatmap)
chr_window = makeWindows(chr_gr, w = 1e6)
chr_window

average_in_window = function(window, gr, v, method = "weighted", empty_v = NA) {

    if(missing(v)) v = rep(1, length(gr))
    if(is.null(v)) v = rep(1, length(gr))
    if(is.atomic(v) && is.vector(v)) v = cbind(v)

    v = as.matrix(v)
    if(is.character(v) && ncol(v) > 1) {
        stop("`v` can only be a character vector.")
    }

    if(length(empty_v) == 1) {
        empty_v = rep(empty_v, ncol(v))
    }

    u = matrix(rep(empty_v, each = length(window)), nrow = length(window), ncol = ncol(v))

    mtch = as.matrix(findOverlaps(window, gr))
    intersect = pintersect(window[mtch[,1]], gr[mtch[,2]])
    w = width(intersect)
    v = v[mtch[,2], , drop = FALSE]
    n = nrow(v)

    ind_list = split(seq_len(n), mtch[, 1])
    window_index = as.numeric(names(ind_list))
    window_w = width(window)

    if(is.character(v)) {
        for(i in seq_along(ind_list)) {
            ind = ind_list[[i]]
            if(is.function(method)) {
                u[window_index[i], ] = method(v[ind], w[ind], window_w[i])
            } else {
                tb = tapply(w[ind], v[ind], sum)
                u[window_index[i], ] = names(tb[which.max(tb)])
            }
        }
    } else {
        if(method == "w0") {
            gr2 = reduce(gr, min.gapwidth = 0)
            mtch2 = as.matrix(findOverlaps(window, gr2))
            intersect2 = pintersect(window[mtch2[, 1]], gr2[mtch2[, 2]])

            width_intersect = tapply(width(intersect2), mtch2[, 1], sum)
            ind = unique(mtch2[, 1])
            width_setdiff = width(window[ind]) - width_intersect

            w2 = width(window[ind])

            for(i in seq_along(ind_list)) {
                ind = ind_list[[i]]
                x = colSums(v[ind, , drop = FALSE]*w[ind])/sum(w[ind])
                u[window_index[i], ] = (x*width_intersect[i] + empty_v*width_setdiff[i])/w2[i]
            }

        } else if(method == "absolute") {
            for(i in seq_along(ind_list)) {
                u[window_index[i], ] = colMeans(v[ind_list[[i]], , drop = FALSE])
            }
            
        } else if(method == "weighted") {
            for(i in seq_along(ind_list)) {
                ind = ind_list[[i]]
                u[window_index[i], ] = colSums(v[ind, , drop = FALSE]*w[ind])/sum(w[ind])
            }
        } else {
            if(is.function(method)) {
                for(i in seq_along(ind_list)) {
                    ind = ind_list[[i]]
                    u[window_index[i], ] = method(v[ind], w[ind], window_w[i])
                }
            } else {
                stop("wrong method.")
            }
        }
    }

    return(u)
}

bed1 = generateRandomBed(nr = 1000, nc = 10) # generateRandomBed() is from circlize package
# convert to a GRanes object
gr1 = GRanges(seqnames = bed1[, 1], ranges = IRanges(bed1[, 2], bed1[, 3]))
#获取统一的结果
num_mat = average_in_window(chr_window, gr1, bed1[, -(1:3)])
#行为genomic ranges， 列为cell
dim(num_mat)

获得num_mat后，自己根据cellid作图即可
参考： https://jokergoo.github.io/ComplexHeatmap-reference/book/genome-level-heatmap.html
我的数据得到的图：