需要做一个感兴趣的基因的测序数据在基因组上的分布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还有更多参数,参考官网即可。
今天的文章
基因分析器怎么完成_DNA基因图谱分享到此就结束了,感谢您的阅读。
版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 举报,一经查实,本站将立刻删除。
如需转载请保留出处:https://bianchenghao.cn/61507.html