利用deeptools来做基因在genome上的分布图

需要做一个感兴趣的基因的测序数据在基因组上的分布profile，deeptools可以解决这个问题，记录一下过程。
1、安装deeptools
pip install deeptools

2、准备bw和bed文件
（1）bw
bw即样本测序数据的bigwig文件，之前做过，也是用deeptools的bamCoverage就可以实现。

bamCoverage -b sample.bam -o sample.bw

（2）bed
bed即感兴趣的基因在genome上的位置，我先用10个基因测试，10个基因的ensembl_id放入10.txt中。
①利用bedtools来将人类基因组注释文件中的‘gene’区域提取出来，并按照bed格式要求转换为txt。
bed格式要求：
第一列：染色体
第二列：起始位置
第三列：终止位置
第四列：基因名
第五列：基因长度
第六列：正负链
这六列在gtf中分别对应1、4、5、9、5-4、7列。其中长度通过5-4计算。

awk -F ";” '{print $1}' Homo_sapiens.GRCh38.91.gtf| awk -F "\t" '{if($3~/gene/) print  $1"\t"$4"\t"$5"\t"$9"\t"$5"-"$4"\t"$7" }'|sed 's/gene_id//g'|sed 's/"//g'> hg38.txt

②利用R，将hg38.txt和10.txt进行关联合并，从而获得10个基因在染色体上的位置信息。

hg38.bed <- read.delim("hg38.txt",sep ="",header = FALSE)
rownames(hg38.bed)<-hg38.bed$V4
10.txt <- read.delim("10.txt", header = FALSE)
rownames(10.txt)<- 10.txt$V1

10.bed <-hg38.bed[match(rownames(10.txt),rownames(hg38.bed)),c(1:6)]
10.bed$V4 <- rownames(I.bed)
write.table(10.bed,file ="10.bed",sep ="\t", quote = FALSE, col.names = FALSE,row.names=FALSE)

到此为止，获得了sample.bw以及10.bed，即可以开始作图。

3、利用deeptools作图
（1）先计算computeMatrix
输出matrix.mat.gz进入下一步

computeMatrix scale-regions -S sample.bw  -R 10.bed --beforeRegionStartLength 1000 --regionBodyLength 5000 --afterRegionStartLength 1000 --skipZeros -o matrix.mat.gz

（2）利用plotProfile绘图

plotProfile -m matrix.mat.gz  -out ExampleProfile1.png --plotTitle "Test data profile"

deeptools还有更多参数，参考官网即可。