查看原文
其他

Biostar:课程13、14

2017-07-20 István Albert 生信媛

翻译:生物女博士

注:本系列课程翻译已获得原作者授权。


#为作者的注释

#*为译者的注释


Lecture 13 - 安装和使用bwa


#从SourceForge上下载bwa

# 从SourceForge获取的URL是比较复杂的,因为SourceForge希望你去访问网址并从网址上下载,这样你可以看到网页投放的广告。

# 你可以访问网址,并手动下载。

#*最新版本是bwa-0.7.15

#*响应时间可能会比较久

cd ~/src

curl -OL http://downloads.sourceforge.net/project/bio-bwa/bwa-0.7.15.tar.bz2


# 这个安装包是用bzip2压缩的,所以需要加-j 来解压。

tar jxvf bwa.0.7.15.tar.bz2


# 我们需要编译一下程序

cd bwa-0.7.15/


# 用make命令来编译。

make


# 把bwa程序链接到bin目录下。

ln -s ~/src/bwa-0.7.15/bwa ~/bin


# 运行一下程序,确保其可用。

mkdir -p ~/edu/lec13; cd ~/edu/lec13


# 运行这个工具,看看它是怎样用的。

bwa


# 手册不错呦!

man ~/src/bwa-0.7.15/bwa.1

#*点击"q"可以退出。


# 创建一个总目录来存储参考序列。

mkdir -p ~/refs/852


# 获取1999年大爆发的埃博拉病毒的参考基因组序列。

#*貌似最近用efetch工具下载NCBI上的数据有问题。大家可以手动下载,然后放到目的文件夹。

efetch -db nucleotide -id NC_002549 -format fasta > ~/refs/852/ebola-1999.fa


# 对基因组创建索引,这样bwa就可以进行基因组搜索。

bwa index ~/refs/852/ebola-1999.fa


# 来看一下运行结果。

ls ~/refs/852/


# 使之也能通过blast进行搜索。

makeblastdb -in ~/refs/852/ebola-1999.fa -dbtype nucl -parse_seqids


# 现在,来看看发生了什么。

ls ~/refs/852/


# 从一个文件中衍生了16个文件,我们最好还是把它们放在一个独立文件夹。

ls -l ~/refs/852/ | wc -l


# 获取埃博拉基因组序列的第一行。

head -2 ~/refs/852/ebola-1999.fa > query.fa


# 用bwa-mem命令将其map回它的基因组。

bwa mem ~/refs/852/ebola-1999.fa query.fa > results.sam


# 与该序列的blast过程比较一下

blastn -db ~/refs/852/ebola-1999.fa -query query.fa > results.blastn



Lecture 14 - 了解SAM格式


# 从lecture13结束的地方继续讲。

head -2 ~/refs/852/ebola-1999.fa  > query.fa


# 做序列比对。

bwa mem ~/refs/852/ebola-1999.fa query.fa > results.sam


# 注意,bwa还是会在屏幕上输出“流”(stream)。

#* stream:这里是指软件运行过程中屏幕上不断出现的一些文字。

# 有两种输出流。它们被称为:标准输出(standard out)和标准错误输出(standard error)。

# 下边的代码是告诉你,你可以如何重定向这两个输出到文件中。

bwa mem ~/refs/852/ebola-1999.fa query.fa 1> results.sam 2> errors.txt


# 有很多技巧可以将一个流重定向到文件,尽管在绝大多数情况下,这是不必要的。但有时确实需要存储错误消息。

# http://www.tldp.org/LDP/abs/html/io-redirection.html


# 根据你有的序列创建“查询序列”。

# 添加或删除一些碱基来创造“插入”和“缺失”。

# 做序列比对。

bwa mem ~/refs/852/ebola-1999.fa query.fa > results.sam


# 查看某些列。

# 了解每列的含义很重要。

#*这里作者讲的过于简略,推荐大家看参考阅读1和2或者官网。

# (这三列是)query id, 比对标识和目标id

#*query id其实就是reads序列的名字。

cat results.sam | cut -f 1,2,3


# 比对起始位置,mapping质量, 简要比对信息表达式

cat results.sam | cut -f 4,5,6


# 双端reads信息:名字,位置和模板长度

cat results.sam | cut -f 7,8,9


# Query序列, query质量, 其他属性

#*Query序列(也就是reads序列)

cat results.sam | cut -f 10,11,12,13,14


# 如果你列两个文件给bwa-mem,它会以配对模式来比对它们。

# 从课程网站上获取示例数据集。

curl -O http://www.personal.psu.edu/iua1/courses/files/2014/data-14.tar.gz


# 解压和比对这个数据集,根据结果回答作业上的问题

bwa mem ~/refs/852/ebola-1999.fa read1.fq read2.fq > results.sam



参考阅读1:http://www.cnblogs.com/emanlee/p/5366610.html

参考阅读2:http://blog.sina.com.cn/s/blog_670445240101l30k.html



欢迎关注我们




您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存