生物技术

3D基因组之Hi-C数据分析(三)

这期主要讲Hi-C数据分析(protocol)

1. 数据过滤(质量筛选,NGS必做^_^)

2.分析HiC数据,得到自己设定分辨率(例如100kb)的热图(由于这部分需要下载基因组,网速较慢,下期介绍)

3.Hi-C数据可视化以及Find Topologically Associating Domains

(1)需要数据:Hi-C交互矩阵:N * N 的矩阵

由第二步得到的Hi-C数据是这样的格式(中间以tab分割)

#tab chrT_001 chrT_002 chrT_003 chrT_004chrT_001 629     164     88      105chrT_002 164      612     175     110

chrT_003 88        175     437     105

chrT_004 105      110     105     278

(2)软件包:tadbit

https://github.com/zilhua/tadbit

官方网站:http://3dgenomes.github.io/TADbit/

(3)计算代码:

#!/user/bin/python
# -*- coding:UTF-8 -*-
'''
Created on ${date}
@author: www.zilhua.com
'''
from pytadbit import Chromosome
# initiate a chromosome object that will store all Hi-C data and analysis
#设置环境,初始化
my_chrom = Chromosome(name='chr19', centromere_search=True,
					  species='Homo sapiens', assembly='NCBI36')
###my_chrom
'''
Chromosome chr19:
   0  experiment loaded:
   0  alignment loaded:
   species		 : Homo sapiens
   assembly version: NCBI36
'''
#加载数据以BingRen的数据为例
# load Hi-C data
my_chrom.add_experiment('k562', cell_type='wild type', exp_type='Hi-C',
	identifier='k562',project='TADbit tutorial',
		hic_data="scripts/sample_data/HIC_k562_chr19_chr19_100000_obs.txt",
	resolution=100000)
my_chrom.add_experiment('gm06690',cell_type='cancer', exp_type='Hi-C',
	identifier='gm06690',project='TADbit tutorial',
	hic_data="scripts/sample_data/HIC_gm06690_chr19_chr19_100000_obs.txt",
	resolution=100000)
#参数解释:k562 样本名;cell_type:细胞类型
#exp_type:实验类型;hic_data:HiC数据;
#resolution:分辨率
#对两个样本分别可视化
my_chrom.experiments["k562"].view()
my_chrom.experiments["gm06690"].view()

figure_1-1024x532gm06690-1024x532

#如果这两个样本是对照实验,则可以这样exp = my_chrom.experiments["k562"] + my_chrom.experiments["gm06690"]
my_chrom.add_experiment(exp)
my_chrom.visualize([('k562', 'gm06690'), 'k562+gm06690'])

k562_gm06690-1024x388

#Find Topologically Associating Domainsmy_chrom.find_tad('k562', n_cpus=4)
my_chrom.find_tad('gm06690', n_cpus=4)
#参数解释:函数find_tad第一个参数k652样本名;
#n_cpus线程数,个人电脑一般是1--4,服务器看配置
#获取计算结果
my_chrom.experiments["k562"].tads
my_chrom.experiments["gm06690"].tads
#将数组转换成txt文本
my_chrom.experiments["k562"].write_tad_borders()
#结果中,第一列就是TAD的编号,
#第二列第三列分别是TAD的起始和终止位置
#第四列打分
#将TAD的边界在热图上显示
my_chrom.visualize('k562', paint_tads=True)

k562_bd-1024x532

#如果需要对比两个样本的TAD情况my_chrom.visualize([('k562', 'gm06690')],
paint_tads=True, focus=(490,620), normalized=True)
#参数解释:paint_tads是否描述TAD边界;
#focus 显示染色体位置;normalized是否标准化

tutorial_1_general_45_11

#TAD density
my_chrom.tad_density_plot('k562')

tutorial_1_general_50_0

#保存和重载数据my_chrom.save_chromosome("some_path.tdb", force=True)
from pytadbit import load_chromosome
my_chrom = load_chromosome("some_path.tdb")
print my_chrom.experiments

发表评论