搜索
查看: 2499|回复: 4

[Other] Biostar:课程15、16-samtools简介&利用readseq对GenBank文件格式转换

[复制链接]

33

主题

46

帖子

230

积分

管理员

Rank: 9Rank: 9Rank: 9

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

翻译:生物女博士
注:本系列课程翻译已获得原作者授权。

#为作者的注释
#*为译者的注释


Lecture 15 - BAM文件和samtools

# 安装一个短read模拟器。Heng Li写的wgsim。
cd ~/src
cd wgsim
gcc -g -O2 -Wall -o wgsim wgsim.c -lz -lm
ln -s ~/src/wgsim/wgsim ~/bin/wgsim

# 检查一下是否可以用。
wgsim

# 现在,下载和安装samtools。
# 把samtools链接到~/bin目录。
cd ~/src
tar jxvf samtools-1.1.tar.bz2
cd samtools-1.1
make
ln -s ~/src/samtools-1.1/samtools ~/bin/

# 检查一下是否可以用。
samtools

# 查看samtools的手册。
man ~/src/samtools-1.1/samtools.1

# 我们将会利用一个突变的参考基因组来生成一些短reads,并map(回贴?)回去。
mkdir ~/edu/lec15
cd ~/edu/lec15

# 生成一些数据,并把输出保存到一个文件。
wgsim -N 10000 ~/refs/852/ebola-1999.fa read1.fq read2.fq > mutations.txt

# 运行比对。
bwa mem ~/refs/852/ebola-1999.fa read1.fq read2.fq > results.sam

# 开始转换成BAM文件。
#*运行的时候请注意,这里用的samtools版本比较旧,新版本有的命令已经改变了。
samtools view -Sb results.sam > temp.bam

# 对比对结果进行排序。
samtools sort -f temp.bam results.bam

# 索引比对结果。
samtools index results.bam
#*做完了这一步,如果想方便查看一下比对结果,可以下载一个IGV软件,并将参考序列、bam文件和bam文件的索引后生成的一个文件导入IGV即可。

# 通常,我们会把所有的命令放到一个脚本中并运行。
# 查看文末的脚本。
#*原文并未将这个脚本放上来。但是后边的课程中会有更加完善的脚本。
#*本行命令暂时不运行:bash align.sh read1.fq read2.fq results.bam

# 使用samtools过滤。
# -f 匹配标识(flag)(保留匹配的标识对应的reads)
# -F 过滤标识(flag)(去除匹配标识对应的reads,保留剩余的)
#*在上两个lecture中有说过sam格式,sam每一列有不同含义,其中一列(第二列)叫标识(flag),不同数字有不同含义
#*下边内容引用自参考资料1

FLAG, 概括出一个合适的标记,各个数字分别代表
1     序列是一对序列中的一个
2     比对结果是一个pair-end比对的末端
4     没有找到位点
8     这个序列是pair中的一个但是没有找到位点
16   在这个比对上的位点,序列与参考序列反向互补
32   这个序列在pair-end中的的mate序列与参考序列反响互补
64   序列是 mate 1
128 序列是 mate 2
假如说标记为以上列举出的数目,就可以直接推断出匹配的情况。假如说标记不是以上列举出的数字,比如说83=(64+16+2+1),就是这几种情况值和。

# 比对到反向互补链的reads。
samtools view -f 16 results.bam

# 比对到正向链的reads。
samtools view -F 16 results.bam

# -c 可以用来计数
samtools view -F 16 results.bam | wc -l
samtools view -c -F 16 results.bam

# 用质量来过滤。BWA 的mapping质量为0时,表示reads是map到多个位置。而 q>=1表示该reads是单一mapping。
samtools view -c -q 1 results.bam

# 统计高质量的比对数目。
samtools view -c -q 40 results.bam

# 参考基因组的名字为 gi|10313991|ref|NC_002549.1|
# 总是手动输入会有些乏味,我们可以创建一个简称。
CHR='gi|10313991|ref|NC_002549.1|'

# 截取部分数据。
samtools view results.bam $CHR:1-100

# 查看文件的特定某列
samtools view results.bam $CHR:1-100 | cut -f 4 |  more


Lecture 16 - Converting data with ReadSeq

# 安装文件格式转换器。
# 为之创建一个目录。
mkdir -p x
cd ~/src/readseq

# 跟调用trimmomatic很相似
# 将调用trimmomatic脚本修改一下。
#*在lecture8时,我们创建了一个文件叫trimmomatic来简化trimmomatic的命令执行。详情参见Biostar:课程7、8
cp ~/bin/trimmomatic ~/bin/readseq

# 用下边程序来替换:
# java -jar ~/src/readseq/readseq.jar $@
# 你也可以加上:
# java -cp ~/src/readseq/readseq.jar app $@
# 用图形用户界面来运行。

# 获取完整的GenBank记录的1999年的埃博拉基因组。

# 获取完整genbank文件。
mkdir -p ~/edu/lec16
cd -p ~/edu/lec16

# 获取完整数据。
#*efetch貌似有问题,可以用网页直接下载。详情参见Biostar:课程3、4
efetch -db nucleotide -id NC_002549.1 -format gb > NC.gb

# 提取数据未GFF(Generic Feature Format)格式。
# 自动获取输入文件的格式。
readseq -format=GFF NC.gb

# 你也可以设置输出文件的名字。
readseq -format=GFF -o NC-all.gff NC.gb

# 提取为fasta文件。
readseq -format=FASTA -o NC.fa NC.gb

# 从GFF文件中提取“gene”的行
#* \t表示的是tab分隔符
cat NC-all.gff | egrep '\tgene\t'

# 同时,前两行也是我们要的。
echo '##gff-version 2' > NC-genes.gff
cat NC-all.gff | egrep '\tgene\t' >> NC-genes.gff

# 染色体坐标不匹配!我们需要去把它修复一下。
# 我们需要重做比对,或者交换现有比对或特征的坐标。
# 重新索引和比对数据。
cp NC.fa ~/refs/852/
bwa index ~/refs/852/NC.fa
cp ../lec15/*.fq .

# 编辑align.sh使其能用于新文件,并重新运行。
bash align.sh read1.fq read2.fq results.bam

# 把数据导入IGV并查看。
# 查看100-500bp区域的深度。
# 输出为序列号、基因组索引和覆盖度。
samtools depth -r NC_002549:100-500 results.bam


参考资料1:http://www.cnblogs.com

欢迎到微信公众号订阅我们
生信媛
bio_sxy







上一篇:samtools安装报错
下一篇:Biostar:课程17、18-awk进行简单的编程&序列比对软件bwa和bowtie2
回复

使用道具 举报

0

主题

7

帖子

51

积分

注册会员

Rank: 2

积分
51
发表于 2017-10-12 17:18:42 | 显示全部楼层
ubuntu 17.04的samtools 1.6在安装时会报错,需要安装
apt-get install libncurses5-dev
apt-get install libbz2-dev
apt-get install liblzma-dev
回复 支持 反对

使用道具 举报

0

主题

7

帖子

51

积分

注册会员

Rank: 2

积分
51
发表于 2017-10-12 19:48:22 | 显示全部楼层
align.sh文件是?
回复 支持 反对

使用道具 举报

33

主题

46

帖子

230

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
230
 楼主| 发表于 2017-10-13 09:16:08 | 显示全部楼层

原文没提供这个脚本。不过按照文意是,把如下几步放入脚本,命名未align.sh
# 运行比对。
bwa mem ~/refs/852/ebola-1999.fa read1.fq read2.fq > results.sam

# 开始转换成BAM文件。
#*运行的时候请注意,这里用的samtools版本比较旧,新版本有的命令已经改变了。
samtools view -Sb results.sam > temp.bam

# 对比对结果进行排序。
samtools sort -f temp.bam results.bam

# 索引比对结果。
samtools index results.bam

并对应修改 read1.fq read2.fq results.bam 为变量(类似Lec22那个脚本那样修改)
回复 支持 反对

使用道具 举报

0

主题

7

帖子

51

积分

注册会员

Rank: 2

积分
51
发表于 2017-10-13 10:54:17 | 显示全部楼层
生信媛 发表于 2017-10-13 09:16
原文没提供这个脚本。不过按照文意是,把如下几步放入脚本,命名未align.sh
# 运行比对。
bwa mem ~/refs ...

谢谢
回复 支持 反对

使用道具 举报

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

本版积分规则

QQ|手机版|小黑屋|生信技能树    

GMT+8, 2019-5-22 17:54 , Processed in 0.052821 second(s), 27 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.