搜索
查看: 3937|回复: 0

[Other] Biostar:课程25、26:bedtools的安装和使用

[复制链接]

33

主题

46

帖子

230

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
230
发表于 2017-9-24 14:22:09 | 显示全部楼层 |阅读模式
注:本文是生信媛微信公众号原创文章
作者:István Albert
原文链接:https://mp.weixin.qq.com/s/Gx7edC6re8gb2ffuZhP4GA


Lecture 25 - Bedtools简介

# 安装bedtools
#* 后台有人发消息来问,你怎么知道你要找的这个软件的链接在哪里?
#* 其实,很多官网上,都是有的。跟你下个office、PS等这些软件是一样的道理。
#* 比如今天要学习的bedtools
#* 首先,百度或者google关键词:bedtools
#* 搜索得到的词条,前两个都是官网的(有时候不一定排在很前面)。
#* 进入官网以后,就找到了。可以看到是在github上存放的。

#* 先大致看一下网页,对于github,可以点击“clone or download”点击下载或者获得下载的链接。


#* 具体下载方式见下图:

#* 在本软件的下载中,在页面下方,还有另外一个安装包下载链接:

#* 留给有兴趣的读者自行探索。
#* 其他也是相类似的搜索。如果是R包,则可以到CRAN(https://cran.r-project.org/)上的Packages(https://cran.r-project.org/web/packages/index.html)或者是bioconductor(http://www.bioconductor.org/)网站去找。


#* 好了,进入今天的正题。

# 安装bedtools

cd ~/srccurl -OL https://github.com/arq5x/bedtools2/releases/download/v2.22.0/bedtools-2.22.0.tar.gz
tar zxvf bedtools-2.22.0.tar.gz
cd bedtools2
make
ln -sf ~/src/bedtools2/bin/bedtools ~/bin/bedtools


# 演示版的bed文件 (demo.bed)

[AppleScript] 纯文本查看 复制代码
KM034562    100    200    one    0    +KM034562    400    500    two    0    -


# 我们的基因组文件(genome.txt)
[AppleScript] 纯文本查看 复制代码
KM034562    18959


# 没有链特异性的运算

[AppleScript] 纯文本查看 复制代码
bedtools slop -i demo.bed -g genome.txt -l 10 -r 0


# 有链特异性的运算

[AppleScript] 纯文本查看 复制代码
bedtools slop -i demo.bed -g genome.txt -l 10 -r 0 -s


# 把BED转成对应的GFF
# 这并非是真的正确地把BED转成GFF。

[AppleScript] 纯文本查看 复制代码
cat demo.bed | bioawk -c bed '{print $chrom, ".", ".", $start+1, $end, $score, $strand, ".", "." }' > demo.gff


# 看!它与其他格式可以很好地协同工作!那怎么可能呢?好神奇!
[AppleScript] 纯文本查看 复制代码
bedtools slop -i demo.gff -g genome.txt -l 10 -r 0 -s
# 两侧的运算。把结果重定向到一个文件中。
[AppleScript] 纯文本查看 复制代码
bedtools flank -i demo.gff -g genome.txt -l 10 -r 0 -s > flank.gff
# 填充运算。
[AppleScript] 纯文本查看 复制代码
bedtools complement -i demo.gff -g genome.txt > complement.gff
# 获取埃博拉基因组。
#* KM034562.gb文件在之前的课程其实已经有给过大家了。
#* 或者也可在这里下载,链接:https://share.weiyun.com/e880cca38636c8d90780fb487050bf3d (密码:HhCUk9)

[AppleScript] 纯文本查看 复制代码
efetch -id KM034562 -format gb -db nucleotide > KM034562.gbreadseq -format=FASTA -o ~/refs/852/KM034562.fa KM034562.gbreadseq -format=GFF -o KM034562.gff KM034562.gb
# 看一下提示,这个工具可以用于哪些格式。
[AppleScript] 纯文本查看 复制代码
bioawk -c help
# 从整个文件中提取基因。
[AppleScript] 纯文本查看 复制代码
cat KM034562.gff | bioawk -c gff ' $feature=="gene" { print $0 } ' > genes.gff
# 序列提取。 获取与间隔相对应的序列。
[AppleScript] 纯文本查看 复制代码
bedtools flank -i genes.gff -g genome.txt -l 10 -r 0 -s > flank.gff
#* 在这一步,使用最新版本的bedtools会有奇怪的报错。而用作者提供的这个旧版本则不会。

# 提取与genes.gff的间隔相对应的序列。
[AppleScript] 纯文本查看 复制代码
bedtools getfasta -fi ~/refs/852/KM034562.fa -bed genes.gff -fo genes.fa
# 来看一下结果文件。
[AppleScript] 纯文本查看 复制代码
head genes.fa

Lecture 26 - Bedtools指南
# 来到ecture 23,我们进行了埃博拉的短reads比对
cd ~/edu/lec23

# 创建一个间隔文件,叫region.bed
# KM034562        1000        2000
# 我们可以用这个间隔文件去分割数据。
bedtools intersect -a bam/SRR1553593.bam -b region.bed > region.bamsamtools index region.bam

# 我们也可以把数据作为bed文件输出。
bedtools intersect -a bam/SRR1553593.bam -b region.bed -bed > region.bed

# 又或者,我们可以把bam文件转成bed文件。
bedtools bamtobed -i bam/SRR1553595.bam > SRR1553595.bed

# 我们将跟着bedtools的指南来做
http://quinlanlab.org/tutorials/cshl2014/bedtools.htmlmkdir ~/edu/lec26cd ~/edu/lec26

# 获取所有的数据,并解包。

# 这个数据集含有更多的数据。
# 胎儿的组织样品,包括脑、心脏、肠道、肾脏、肺、肌肉、皮肤以及胃
tar -zxvf maurano.dnaseI.tgz

# 我们将在课堂上跟着指南做一点实际操作
# 并给作业留一些建议
# 从注释文件中,选取启动子。
cat hesc.chromHmm.bed | grep  Promoter > promoters.bed

# 找到跟每个exon最近的启动子。
# -d 表示打印列
#* 多的一列数值是-a 和 -b 两者最近的距离。(见下图)
bedtools closest -a exons.bed -b promoters.bed  -d | head -1


# 找到覆盖了最多外显子的CPG岛
bedtools intersect -a cpg.bed -b exons.bed -c | sort -k5,5nr | head -1

# 以5Kb一个窗口把人类基因组以覆盖。
bedtools makewindows -g genome.txt -w 50000 > windows.bed
欢迎到微信公众号订阅我们
生信媛
bio_sxy






本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x



上一篇:推荐本书《统计学习方法》-李航
下一篇:Biostar:课程27、28:RNAseq分析
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-9-21 14:15 , Processed in 0.042292 second(s), 29 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.