bed文件转成等bin的bedgraph格式
2020-04-07 本文已影响0人
caokai001
目的:
当你想可视化你所研究的数据在染色体上分布,可能会用到RIdeogram来画类似的图。如何整理成要求的输入文件格式呢?
![](https://img.haomeiwen.com/i9589088/bd8ab3f2e7ab3609.png)
标准输入文件: 一定长度为bins 的count 数统计文件(类似bedgraph格式)
![](https://img.haomeiwen.com/i9589088/dc45ebfd67b5ba70.png)
假设你手里有ChIP-seq测序结果的bed 文件,如何得到等bin区间的bedgraph 结果呢。
你可以选择bed 转成bam再转成bedgraph:
工具列表 bedtools bedToBam ;deeptools bamCoverage .
但是你是否想过这样得到的结果,也就是bedgraph 文件存在一个问题,相同的value的区间会自动合并,你可以通过下面操作让bin 不合并)
实践:
输入文件
csi.chromosome.fa.fai : 基因组samtools faidx 索引文件
50bpC5_1.id0.9co50bp.nochrUn : 比如ChIP-seq uniq.bamtobed 文件
代码
### 将基因组以10k 为bin进行分割
awk '{n=int($2/10000);for(i=0;i<=n+1;i++){print $1"\t"i*10000"\t"(1+i)*10000}}' csi.chromosome.fa.fai > csi.chromosome.10k.bedgraph
### 将bed 文件进行排序,注意strand 正负向问题
$ less -S 50bpC5_1.id0.9co50bp.nochrUn |cut -f 2,9,10 |awk 'BEGIN{FS=OFS="\t"}{if($3>$2)print $1,$2,$3;else print $1,$3,$2}' > 50bpC5_1.id0.9co50bp.nochrUn.sort.bed
### 利用bedtools coverage 得到bedgraph 文件(相同value的bin不会合并)
$ bedtools coverage -a csi.chromosome.10k.bedgraph -b 50bpC5_1.id0.9co50bp.nochrUn.sort.bed |cut -f 1-4 >50bpC5_1.id0.9co50bp.nochrUn.sort.bedgraph
结果
![](https://img.haomeiwen.com/i9589088/a2d4f103c1d11b8a.png)
补充代码:
![](https://img.haomeiwen.com/i9589088/622e037dc53b9029.png)
思考:
1.bedtools makewindows 也可以分bin.
2.输入文件是普通的bedgraph 可能也可以,为了保险起见,开始用bedtools 来按照每个bin分别统计count 数。