参考:
比如原位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
分辨率。
总结
- 数据量达不到要求,得到的Loop结果数量差别就很大。
- 如果有其他相关工具,可以直接计算最佳分辨率,欢迎评论留言~