我们知道.hic 文件是高度压缩的二进制文件,便于存储和分析。那么如果我们想要从.hic提取某一区域的交互信息的话,该如何操作呢?这就涉及到了juicer dump。
https://github.com/aidenlab/juicer/wiki/Data-Extraction
Juicer dump 有以下参数:
Usage:
juicebox dump <observed/oe> <NONE/VC/VC_SQRT/KR> <hicFile(s)> <chr1>[:x1:x2] <chr2>[:y1:y2] <BP/FRAG> <binsize> [outfile]
dump <norm/expected> <NONE/VC/VC_SQRT/KR> <hicFile(s)> <chr> <BP/FRAG> <binsize> [outfile]
dump <loops/domains> <hicFile URL> [outfile]
示例:
juicer_tools dump observed NONE sam1.chr20.hic 20:32679500:32680500 20 BP 10000 extract_matrix.txt
提取的矩阵主要有三列:(start,end,contacts)
提取矩阵示例:
120000 32680000 1.0
350000 32680000 2.0
370000 32680000 1.0
560000 32680000 1.0
850000 32680000 1.0
980000 32680000 1.0
1190000 32680000 2.0
1270000 32680000 1.0
1300000 32680000 1.0
1800000 32680000 1.0
那么如果我们想要进行可视化的话,可以参照以下代码转换成HiTC格式的矩阵:
def reform_matrix(self):
#-----------HiTC matrix---------------------
chr=self.chr;start=self.start;end=self.end;bin=self.bin;genome=self.genome
self.hitc_matrix="{}/{}_{}_{}_{}_hitc.mat".format(self.outdir,self.prefix,chr,start,end)
mat=pd.read_table(self.juicer_dump_mat,names=['frag1','frag2','contacts'])
min=math.ceil(int(start)/bin)*bin
max=int(int(end)/bin)*bin
N=int(end/bin)-math.ceil(start/bin)+1
#---------------------- add header --------------------------
inddf=np.arange(N)
headers_ref=[genome for x in inddf]
bin_num_df=pd.Series(inddf).apply(lambda x : str(x))
headers_ref=pd.Series(headers_ref)
chromdf=pd.Series([chr for x in list(range(N))])
startdf=pd.Series(inddf*bin+min)
enddf=pd.Series((inddf+1)*bin+min)
headers_suf=chromdf.str.cat(startdf.apply(lambda x :str(x)),sep=':')
headers_suf=headers_suf.str.cat(enddf.apply(lambda x:str(x)),sep="-")
headers=bin_num_df.str.cat([headers_ref,headers_suf],sep="|")
headers=list(headers)
mat['b1']=mat['frag1'].apply(lambda x: (x-min)/bin)
mat['b2']=mat['frag2'].apply(lambda x: (x-min)/bin)
counts=sparse.coo_matrix((mat['contacts'],(mat['b1'],mat['b2'])),shape=(N, N),dtype=float).toarray()
diag_matrix=np.diag(np.diag(counts))
counts=counts.T + counts
counts=counts-diag_matrix
df=pd.DataFrame(counts)
df.columns=headers
df.index=headers
#print('DataFrame.....')
#print(df.head())
df.to_csv(self.hitc_matrix,sep="\t")
return df
来查看一下结果。
image.png
网友评论