作图seq 比对

bed文件转成等bin的bedgraph格式

2020-04-07  本文已影响0人  caokai001

目的:

当你想可视化你所研究的数据在染色体上分布,可能会用到RIdeogram来画类似的图。如何整理成要求的输入文件格式呢?

image.png

标准输入文件: 一定长度为bins 的count 数统计文件(类似bedgraph格式)

假设你手里有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

结果

image.png

补充代码:

image.png

思考:

1.bedtools makewindows 也可以分bin.
2.输入文件是普通的bedgraph 可能也可以,为了保险起见,开始用bedtools 来按照每个bin分别统计count 数。

上一篇 下一篇

猜你喜欢

热点阅读