三维基因组: SELFISH 差异分析

引言

本系列主要讲解 3D-Genome (Hi-C) 系列的分析，主要涉及三维基因组分析中的数据处理，重复性评估，Compartment/TAD/Loop 检测，差异分析等，欢迎关注！

SELFISH

SELFISH 是一款用于检测 Hi-C 矩阵间差异相互作用的软件，支持 MATLAB 和 Python。

原理

基于矩阵间局部自相似性原理。如果矩阵 A 和 B 之间，基因组 bin_i 和 j 的相互作用接触频率存在显著差异，那么这种差异在该相互作用周围的像素（即“影响半径”）内仍会显著。为了降低远离测试相互作用的相互作用的影响，会以 i, j 为中心，对影响半径施加一个逐渐增大的高斯滤波器。同时，由于基因组位点的线性邻近性和布朗运动，大多数相互作用发生在矩阵对角线附近，因此会根据与对角线的接近程度对相互作用进行归一化处理。最终，每个相互作用用一组向量表示，这些向量包含其邻近相互作用在不同影响半径下的频率，而这组向量的一阶导数将用于检测矩阵 A 和 B 之间的差异。

Example

本文讨论的 SELFISH 的 Python 版本是一个命令行工具，可接受 .hic、. cool 和 HiC-Pro 矩阵作为输入。它输出的文件格式有两种：一种是 numpy 格式的差异相互作用概率 p 值矩阵；另一种是更易读的制表符分隔文件，其中包含两个差异相互作用 bin 的起始坐标、p 值以及差异相互作用的对数变化倍数（logFC）。

SELFISH 的运行命令如下：

selfish -f1 HiC_Rep1.hic -f2 HiC_Rep1.hic -ch chr2 -r 5kb -o SELFISH_chr2_5kb_D00vsD15.tsv -t 0.05

在这里，-f1 和 -f2 是需要分析的矩阵，-ch 和 -r 是从 .hic 文件中提取的染色体和分辨率，-t 参数用于指定输出结果为制表符分隔文件，同时设置记录结果的 p 值阈值（此处设定得较为宽松，<0.05，后续可以对结果进行进一步筛选），而 -o 参数则用于指定结果文件的路径和文件名。尽管可以通过参数 -p 绘制每个染色体的差异相互作用图，但对于分辨率更高的分析（例如，bin 尺寸小于 50 kb），强烈不推荐使用此选项。

在将结果转换为 R 中的 .bedpe 格式后，SELFISH 的结果可以在 JuiceBox 中轻松可视化：

# R version 3.6.2

options(scipen=999)

bin.size <- 5000
pval_threshold <- 0.01
logFC_threshold <- 1
chrs <- c("chr1","chr2","chr3","chr4","chr5","chr6","chr7","chr8","chr9","chr10","chr11","chr12","chr13","chr14","chr15","chr16","chr17","chr18","chr19","chr20","chr21","chr22","chrX")
all_interactions <- NULL

for (chr in chrs) {
    print(chr)
    diffInts <- read.table(paste("SELFISH", chr, "D00vsD15_5kb.tsv",sep="_"), sep= "\t", header= T, stringsAsFactors= F)
    diffInts_filtered <- diffInts[diffInts$P_VAL < pval_threshold & (diffInts$LOG_FOLD_CHANGE > logFC_threshold | diffInts$LOG_FOLD_CHANGE < -(logFC_threshold)),]
    Nrow <- nrow(diffInts_filtered)
    diffInts_filtered_bedpe <- data.frame(
    "chr1"= rep(chr, Nrow),
    "chr1_bin_start"= diffInts_filtered$LOC1, "chr1_bin_end"= diffInts_fil-
    tered$LOC1 + bin.size,
    "chr2"= rep(chr, Nrow),
    "chr2_bin_start"= diffInts_filtered$LOC2,
    "chr2_bin_end"= diffInts_filtered$LOC2 + bin.size,
    "name"= rep(".", Nrow),
    "score"= rep(".", Nrow),
    "strand1"= rep(".", Nrow),
    "strand2"= rep(".", Nrow),
    "color"= rep("0,255,255", Nrow),
    "p_value"= diffInts_filtered$P_VAL,
    "logFC"= diffInts_filtered$LOG_FOLD_CHANGE, stringsAsFactors=F)
    all_interactions <- rbind(all_interactions, diffInts_filtered_bedpe, stringsAsFactors=F)
}

all_interactions_noEmpty <- all_interactions[!is.na(all_interactions$p_value),]

# RGB color code for dark blue
all_interactions_noEmpty[all_interactions_noEmpty$logFC>0,11] <- "0,0,139" 

# RGB color code for orange
all_interactions_noEmpty[all_interactions_noEmpty$logFC<0,11] <- "255,140,0" 

write.table(all_interactions_noEmpty, "SELFISH_D00vsD15_5kb.bedpe", sep= "\t", col.names= F, row.names= F, quote=F)

在这个过程中，通过 bin. size 设置分辨率，并利用 pval_threshold 和 logFC_threshold 来确定 p 值和 logFC 的阈值，以便筛选出符合条件的相互作用。使用 read. table 函数读取 SELFISH 输出的 .tsv 文件，然后根据设定的 p 值阈值（<0.01）和对数变化倍数阈值（logFC >1 或 < -1）对相互作用列表进行筛选。为了使结果符合 .bedpe 格式，还添加了名称、得分、strand1、strand2 和相互作用颜色列。最后，根据对数变化倍数的正负，调整相互作用的颜色，并通过 write. table 函数将筛选后的差异相互作用以 .bedpe 格式保存到一个制表符分隔文件中。根据这些阈值，SELFISH 在染色体 2 上共发现了 169,703 个差异相互作用。

本文由mdnice多平台发布