使用HicPro處理fastq data达舒，并進行了數(shù)據(jù)格式轉(zhuǎn)換以后哭懈，就可以進行一些常規(guī)的HiC下游計算了，最常規(guī)的包括AB compartment挫鸽，TAD相關(guān)計算缚柳，loop相關(guān)計算埃脏。本文整理了我自己使用的AB區(qū)室計算的相關(guān)方法和代碼。
一秋忙、HiC ValidPair修飾
hicpro 軟件輸出的validPair文件彩掐，經(jīng)過簡單修飾后，再使用homer進行compartment計算PC1值灰追，PC1>0的區(qū)域為A區(qū)室堵幽，PC1<0為B區(qū)室。
homer軟件只需要validPair的前7列：
cat CT.allValidPairs |awk 'BEGIN{FS=OFS="\t"}{print $1,$2,$3,$4,$5,$6,$7}' > CT.allValidPairs.homer
二监嗜、call PC1 value —— homer
## step 1 : makeTagDir
tagdir=~/data/hic/homer/tagDir
makeTagDirectory ${tagdir}/CT -format HiCsummary CT.allValidPairs.homer
## step 2 : calculate PC1 value
## 如果有ATAC 測序數(shù)據(jù)谐檀，則ATAC narrowPeak文件可以作為 seed 輸入抡谐，輔助計算PC1
mm10=~/data/database/Mus_musculus.GRCm38.dna.toplevel.fa
seed=~/data/hic/homer/seed_active/
output=~/data/hic/homer/compartment/
runHiCpca.pl ${output}-500k ${tagdir}/CT -cpu 20 -res 500000 -genome $mm10 -pc 1 -active ${seed}/CT.narrowPeak
## output data,兩個output
head -n 5 CT-500k.PC1.bedGraph
track name="..." yLineMark="0.0" alwaysZero=on maxHeightPixels=100:75:11 visibility=full viewLimits=-1:1 autoScale=on type=bedGraph
chr11 3080000 3120000 -0.220452952221998
chr11 3120000 3160000 -0.200718202287763
chr11 3160000 3200000 0.0416765017206624
chr11 3200000 3240000 0.343580286630125
[medhdl@login2 compartment]$ head -n 5 CT-500k.PC1.txt
#peakID chr start end strand PC1
chr11-3080000 chr11 3080000 3120000 + -0.510452952221998
chr11-3120000 chr11 3120000 3160000 + -0.314718202287763
chr11-3160000 chr11 3160000 3200000 + 0.0806765017206624
chr11-3200000 chr11 3200000 3240000 + 0.894380286630125
以上運行所得的矩陣可以導入R中進行相關(guān)下游分析裁奇，非常方便。
三麦撵、使用PC1 value 進行后續(xù)分析
變成文本了以后刽肠，用R進行統(tǒng)計就可以。
1免胃、標記 AB compartment
通過PC1的正負符號標記 AB compartment音五。
#為了方便理解，這里盡量了使用最簡單的Rcode
pc1.ct <- read.table("CT-500k.PC1.txt") #讀入計算所得的PC1矩陣
colnames(pc1.ct) <- c("peakID","chr","start","end","strand","PC1")
pc1.A <- read.table("A-500k.PC1.txt") #讀入計算所得的PC1矩陣
colnames(pc1.A) <- c("peakID","chr","start","end","strand","PC1")
pc1.dat <- merge(pc1.ct,pc1.A,by=c("peakID","chr","start","end","strand"))
colnames(pc1.dat)[6:7] <- c("CT","A") #合并矩陣
pc1.dat$ABtrans <- ifelse(pc1.dat$CT < 0 & pc1.dat$A < 0,"BB",
ifelse(pc1.dat$CT > 0 & pc1.dat$A > 0,"AA",
ifelse(pc1.dat$CT < 0 & pc1.dat$A > 0,"BA","AB"))) #標記區(qū)室信息
2羔沙、計算區(qū)室轉(zhuǎn)換的比例
統(tǒng)計AB區(qū)室轉(zhuǎn)換的數(shù)量和比例躺涝。
pc1.trans <- pc1.dat %>% group_by(ABtrans) %>%
summarise(n=n()) %>% as.data.frame() #計算四類區(qū)室轉(zhuǎn)換的染色質(zhì)區(qū)域個數(shù)
pc1.trans$percentage <- round(x=(pc1.trans$n)/sum(pc1.trans$n),digits = 2) #計算四類區(qū)室轉(zhuǎn)換的占比
pc1.trans$label <- paste0(pc1.trans$ABtrans,"(",pc1.trans$percentage,"%)") #設(shè)置畫圖標記
3、畫圖
使用ggplot2包簡單畫個餅圖扼雏，存成svg格式坚嗜，再在PPT中修改格式細節(jié)。
p <- pc1.transW7vsM25.pctrans %>%
ggplot(aes(x="",y=percentage,fill=ABtrans))+
geom_bar(stat="identity")+
coord_polar(theta = "y")+
labs(title = "ABtrans_portion")+
scale_fill_manual(values = brewer.pal(9,"Blues")[c(2,4,5,3)],labels=pc1.trans$label)+
theme_bw(base_line_size = NA)+
theme(axis.text = element_text(size = 12,color="black"),legend.position = "top")
svg(file="ABtrans.svg",height=8,width=8) #保存svg格式诗充，方便在window上進行修改
p
dev.off()
svg格式的圖片靈活性還是蠻高的苍蔬，可以根據(jù)需求隨便改【只要不改數(shù)據(jù)】，一個簡單的效果圖：
也可以畫成類似文獻中的bar圖蝴蜓，這個就看大家的喜好啦碟绑，比如：
圖片來源發(fā)表文獻
區(qū)室分析內(nèi)容差不多就這些俺猿，一些其他的聯(lián)合分析，還需要視具體情況而定格仲。例如可以在IGV或WashU上進行Track觀察押袍，聯(lián)合ATAC-seq和RNA-seq數(shù)據(jù)觀察區(qū)室轉(zhuǎn)換與基因表達調(diào)控的聯(lián)系等。
下篇繼續(xù)總結(jié)TAD相關(guān)的分析凯肋。
公眾號原文: HI-C 數(shù)據(jù)分析pipeline（三：AB compartment 區(qū)室計算）
