桑基图/冲击图/河流图（如下图I），可以使用变化的宽条带和堆叠条形图来表示数据在不同分类下的流量变化，通常应用于能源、材料成分、金融等数据的可视化分析。在上一篇文章中我们使用Liger v2整合了单细胞转录组和染色质开放两个组学的数据，就很适合使用桑基图进行可视化。

以下的例子中绘制了三列的桑基图，从左到右分别是scRNA-seq中使用seurat进行分群、liger整合后分群、以及signac分群的结果。因为这并不是在同一个细胞中测了RNA和可及性，所以中间一列实际上是左右两列细胞之和。

接上篇的内容，首先将Liger对象中添加scRNA和scATAC分群的结果。

数据格式整理

• 匹配原始scRNA-seq聚类结果

rownames(Liger@cellMeta)
paste0("rna_",rownames(scRNA@meta.data))
[14881] "rna_TTTGTTGCAAGGCTTT-1"         "rna_TTTGTTGCAGAGCGTA-1"
[14883] "rna_TTTGTTGCATATCTGG-1"         "rna_TTTGTTGGTAAGGTCG-1"        
[14885] "rna_TTTGTTGGTCACTAGT-1"         "rna_TTTGTTGGTCTAACTG-1"        
[14887] "rna_TTTGTTGGTGTCCAAT-1"         "rna_TTTGTTGTCATCTGTT-1"        
[14889] "rna_TTTGTTGTCATGACAC-1"

Liger@cellMeta$rna_clusters = scRNA@meta.data$seurat_clusters[match(rownames(Liger@cellMeta),paste0("rna_",rownames(scRNA@meta.data)))]
pdf("atac-rna_integrated_raw_rna_clusters.pdf")
plotClusterDimRed(Liger,useCluster = "rna_clusters")
dev.off()

• 匹配原始scATAC-seq聚类结果

在这里我的scATAC数据barcode前已经添加了样本信息，下面我去掉了样本信息粘贴上“atac_”

Liger@cellMeta$atac_clusters = signac@meta.data$seurat_clusters[match(rownames(Liger@cellMeta),gsub("sample1_", "atac_", rownames(signac@meta.data)))]
pdf("atac-rna_integrated_raw_atac_clusters.pdf")
plotClusterDimRed(Liger,useCluster = "atac_clusters")
dev.off()

如果采用常规绘制桑基图的方法，比如ggalluvial包，可以展现不同clsuter来源不同数据集的细胞占比情况，但是不能很好的表示原来数据集的分群。使用to_lodes_form函数将数据转换为lodes形式。该函数可以对数据框进行重塑，构成不同分组之间的关系。

data <- Liger@cellMeta[,c("seurat_clusters", "leiden_cluster", "signac_seurat_clusters")]
data = data.frame(data)

corLodes <- to_lodes_form(data, axes = 1:ncol(data), id = "id")

# 绘制桑葚图
pdf(file = "ggalluvial.pdf" , width = 7, height = 6)
ggplot(corLodes, aes(x = x, stratum = stratum, alluvium = Cohort, fill = stratum, label = stratum)) +
    scale_x_discrete(expand = c(0, 0)) +  
    geom_flow(width = 0.2, aes.flow = "forward") + 
    geom_stratum(alpha = 0.9, width = 0.2) +
    scale_fill_manual(values = mycol) +
    geom_text(stat = "stratum", size = 2, color = "black") +
    xlab("") + 
    ylab("") + 
    theme_bw() + 
    theme(
        axis.line = element_blank(),
        axis.ticks = element_blank(),
        axis.text.y = element_blank(),
        panel.grid = element_blank(),
        panel.border = element_blank()
    ) + 
    ggtitle("") + 
    guides(fill = FALSE)                            
dev.off()

就像之前说的，由于左右两边细胞数之和加起来才是中间一列，所以这种方法会产生NA值，也就是白色方块区域。

可以通过去掉NA值进行处理，或者代码中加入na.rm

corLodes <- corLodes[complete.cases(corLodes),]

pdf("alluvial_plot.pdf")
ggplot(corLodes, aes(x = x, stratum = stratum, alluvium = id, fill = stratum, label = stratum)) +
    scale_x_discrete(expand = c(0, 0)) +  
    geom_flow(width = 0.2, aes.flow = "forward",na.rm = TRUE) + 
    geom_stratum(alpha = 0.9, width = 0.2,na.rm = TRUE) +
    scale_fill_brewer(type = "qual", 
                    palette = "Set1")+
    geom_text(stat = "stratum", size = 2, color = "black",na.rm = TRUE) +
    xlab("") + 
    ylab("") + 
    theme_bw() + 
    theme(
        axis.line = element_blank(),
        axis.ticks = element_blank(),
        axis.text.y = element_blank(),
        panel.grid = element_blank(),
        panel.border = element_blank()
    ) + 
    ggtitle("") + 
    guides(fill = FALSE)
dev.off() 

新的问题是，去掉NA值后ggalluvial包作图效果表现不佳。ggalluvial包做出的桑基图矩形之间没有空隙，高度似乎也不好单独调整。因此换别的工具。最终效果为：

代码如下：