三维基因组学基因组

【HiC】数据如何选择合适分辨率?

2020-12-11  本文已影响0人  caokai001

参考:

比如原位HiC,DLO-HiC 实验流程不同,reads利用率存在差别。但是你的数据量适合多大分辨率,是否满足TAD(40Kb),Loop(10Kb/5Kb)呢?


image.png

上面是一个方法,计算该分辨率下,大于1000交互的bin所占的比例,是否大于80%。
粗略的计算一下,10Kb分辨率需要多大的数据量(valid interaction)?

image.png

相反,我们数据量是否达到10Kb分辨率呢,需要计算统计交互大于1000的bin所占比例了。


操作(HiC-Pro)

Step0: 将bedpe 格式修改为allValidirs格式

# awk 语法和echo "" 冒号冲突
for sample in *.bedpe; do 
name=${sample/bedpe/bedpe2hicpro.txt}
echo $name; 
( cat $sample | awk 'BEGIN{FS=OFS="\t"}{print $5,$1,$2,$7,$3,$4,$8}' > ${name} & )
done

Step1 :对hic-pro结果allValidirs格式转换成raw matrix

需要计算的分辨率:1000000 250000 100000 50000 25000 20000 15000 10000

### 转换成不同分辨率raw matrix
(base) [13:29:10] kcao@comput4:~/Work
$ 
sample=clean.bedpe2hicpro.txt

echo $sample ;
for res in 1000000  250000  100000  50000  25000 20000 15000 10000; do \
( nohup sh ./kcao_raw_maps.sh -i $sample -o ${sample/bedpe2hicpro.txt/${res}.raw} \
-r $res -s /ChrSize.bed & );
done 

Step2:计算满足大于1000 contact分辨率bin所占比例

### 分组求和,计算合适的分辨率
(py2_env) [14:40:58] kcao@comput4:~/Work
$ for f in *raw.matrix; do 
# echo $f
 high_bin=$(cat $f | awk  '{x[$1] += $3}END{for (i in x){print i, x[i]}}' | sort -k 1,1n  | awk  '($2>1000){sum=sum+1}END{print sum}');
total_bin=$(cat ${f/raw.matrix/raw_abs.bed} | wc -l);
awk -v high_bin=$high_bin -v total_bin=$total_bin -v name=$f 'BEGIN{print name,high_bin/total_bin}' >> ./hic_res_suitable.txt ;
done

结果:

1000000.raw.matrix 0.933826
250000.raw.matrix 0.920303
100000.raw.matrix 0.90917
50000.raw.matrix 0.899155
25000.raw.matrix 0.878731
20000.raw.matrix 0.860007
15000.raw.matrix 0.81036
10000.raw.matrix 0.668495

从上面结果表明,该数据量满足15Kb 分辨率。

总结

上一篇 下一篇

猜你喜欢

热点阅读