搜索
查看: 3162|回复: 3

[Other] Biostar:课程17、18-awk进行简单的编程&序列比对软件bwa和bowtie2

[复制链接]

33

主题

46

帖子

230

积分

管理员

Rank: 9Rank: 9Rank: 9

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

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

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

Lecture 17 - 用awk进行简单的编程

# 用awk进行简单的编程
# 我们之前已经做过这个了,但是为了让大家我们了解我们正在做什么,我们来再做一遍。
# 提取埃博拉基因组的编码特征、基因以及编码序列。
efetch -db nucleotide -id NC_002549.1 -format gb > NC.gb
readseq -format=GFF -o NC.gff NC.gb

#*可以跳过上边两步,使用拟南芥gff文件。

#*我们可以先来查看一下gff格式是什么样子的
less -S TAIR10_GFF3_genes.gff |head

#*gff文件是tab分隔的文件。第1列是染色体信息;第2列是gff注释数据来源;第3列为特征(feature)即属于gene还是mRNA还是CDS等等;第4和5列分别是这个特征序列的起始和终止位置,第6列是得分,可以是序列相似性比对时的E-values值或者基因预测是的P-values值。”.”表示为空(来自参考资料1);第7列是表示序列的方向:正义链为+,反义链为-;第8列仅为对CDS的注释,表示起始编码的位置,有效值为0、1、2(来自参考资料1);第9列为注释信息。

# 找到每个特征的长度,注意“魔法”。
#*如果使用的是拟南芥的gff文件为操作对象,请自行将之后所有命令的NC.gff改为TAIR10_GFF3_genes.gff

#*$1, $2, $3指的就是第1-3列的数据
cat NC.gff | awk ' { print $1, $2, $3 } ' | head -5

# 这(基本上)等同于截取列。
cat NC.gff | cut -f 1,2,3 | head -5

# 我们可以进行计算。每个特征序列长度是多少?
#*第5列的数值减去第4列的数值后+1,即得到特征序列的长度
cat NC.gff | awk ' { print $3, $5-$4 + 1 } ' | head -5

# 我们可以利用模式匹配来提取CDS特征。
cat NC.gff | awk '$3 =="gene" { print $3, $5-$4 + 1, $9 } '

# 计算所有基因的累积长度。
cat NC.gff | awk '$3 =="gene" { len=$5-$4 + 1; size += len; print "Size:", size } '


# 计算所有CDS的累积长度。
cat NC.gff | awk '$3 =="CDS" { len=$5-$4 + 1; size += len; print "Size:", size } '

# 基因组有多大?有很多方法可以知道,比如从SAM文件的头文件(header)。
samtools view -h ../lec16/results.bam | head -2

# 基因组上有多少是基因。

#*如果是拟南芥,把perc=size/18959对应改成perc=size/119667750,119667750是拟南芥(Col-0)基因组的大小。
#*可以用下边的代码自行计算。
#*cat TAIR10_GFF3_genes.gff |awk '$3 == "chromosome"{len=$5-$4 + 1; size += len; print "Size:", size } '

cat NC.gff | awk '$3 =="CDS" { len=$5-$4 + 1; size += len; perc=size/18959; print "Size:", size, perc } '

# 我们需要使awk只用tab分隔符进行分隔。
# 用空格符分隔虽然有用,但是常常导致烦人的bugs
# 告诉awk,以tab作为F(字段分隔符)和OFS(输出字段分隔符)
# 把下边这行命令放到 .profile 或 .bashrc 文件中
alias awk="awk -F '\t' -v OFS='\t'"

# 使设置生效。
source ~/.profile

# 生成新的含有基因和CDS的gff文件。
# The first line specifies what kind of file this is.
head -1 NC.gff  > NC-genes.gff
head -1 NC.gff  > NC-cds.gff

# 根据特征(features)把文件分开。
cat NC.gff | awk ' $3=="gene" { print $0 }' >> NC-genes.gff
cat NC.gff | awk ' $3=="CDS" { print $0 }'  >> NC-cds.gff

# Sam文件是tab分隔符的,可以用awk来处理。
# 上周的数据中,有多少碱基被覆盖了超过50次?
samtools depth ../lec16/results.bam | awk '$3 > 50 { print $0 } ' | wc -l

# 有多少的模板长度超过50 bp。
samtools view ../lec16/results.bam | awk ' $9 > 50 { print $0 } '  | wc -l

从GFF文件中分离提取基因名字。


[Shell] 纯文本查看 复制代码
$3 == "gene" {
# 通过 ; 分离提取第9列

split($9, x, ";")
# 基因名字是第一个元素。

# 通过空格符切分。基因名字是第二个元素

split(x[1], y, " ")

# 去除基因名字旁边的双引号

name = y[2]

#*对于拟南芥的gff,通过“=”切分。基因名字是第二个元素:
#*split(x[1], y, "=")
#*name = y[2]# 用空格全局替换 " 符号

# 由于 " 是一个特殊字符,我们必须写成 \"

#*反斜杠\表示转义符。

gsub("\"", "", name)

# 打印特征类型、基因名字以及大小。

print $3, name, $5 - $4 + 1}

#*最后,我们可以写成下边这条命令

#*cat TAIR10_GFF3_genes.gff |awk '$3 == "gene"{split($9,x,";");split(x[1], y, " ");name = y[2];gsub("\"", "", name); print $3, name, $5 - $4 + 1 } '


Lecture 18 - 序列比对工具的对比

# 下载并安装比对软件bowtie2
cd ~/src

# Mac OSX上用:

# Linux上用:

# 这个已经是个二进制文件了,所以可以直接执行。
unzip bowtie2-2.2.4-macos-x86_64.zip

# 把比对软件以及相关程序链接到bin文件夹。
ln -s ~/src/bowtie2-2.2.4/bowtie2 ~/bin/
ln -s ~/src/bowtie2-2.2.4/bowtie2-align ~/bin/
ln -s ~/src/bowtie2-2.2.4/bowtie2-build ~/bin/

# 对参考基因组建立索引。
# 这是bowtie2的索引。
bowtie2-build ~/refs/852/NC.fa NC.fa

# 把突变弄成一定的格式,使得可以在浏览器上展示。
# 把它变为gff文件。
#*mutations.txt这个文件是在lecture15中生成的。并且存放在目录~/edu/lec15下。
cat mutations.txt | awk ' {print $1, "wgsim", "mutation", $2, $2, ".", "+",  ".", "." }' > mutations.gff

# 运行比对程序。
#*脚本见文末
bash compare.sh

# 修改脚本使之生成有较大的错误率的fastq文件。
# 当你修改了错误率后,计算一下比对上的reads数目。
samtools view -cF 4 bwa.bam
samtools view -cF 4 bow.bam

最终脚本compare.sh

[Shell] 纯文本查看 复制代码
# 对比两个比对软件的输出。
# 使用方法: bash compare.sh
# 在出错的地方停止运行。
set -ue

# 数据文件名。
READ1=r1.fqREAD2=r2.fq

# 参考序列文件。两个比对软件需要分别对它进行索引。
REFS=~/refs/852/NC.fa

# 从基因组中生成模拟reads。
# 编辑错误率并重新运行这个pipeline。
wgsim -N 10000 -e 0.1 $REFS $READ1 $READ2 > mutations.txt

# 将mutations文件转成GFF文件。
cat mutations.txt | awk -F '\t' ' BEGIN { OFS="\t"; print "##gff-version 2" } { print $1, "wgsim", "mutation", $2, $2, ".", "+",  ".", "name " $3 "/" $4 }' > mutations.gff

# 我们需要添加一个read组到mapping中。

GROUP='@RG\tID:123\tSM:liver\tLB:library1'

# 运行bwa并创建bwa.sam文件。
bwa mem -R $GROUP $REFS $READ1 $READ2 > bwa.sam

# 运行bowtie2并创建bow.sam文件。
bowtie2 -x $REFS -1 $READ1 -2 $READ2 > bow.sam

# 调整bowtie2# bowtie2 -D 20 -R 3 -N 1 -L 20 -i S,1,0.50 -x $REFS -1 $READ1 -2 $READ2 > bow.sam
# 对每个sam文件,都转换成bam(sam的二进制)格式。
for name in *.sam;
do
samtools view -Sb $name > tmp.bamsamtools sort -f tmp.bam $name.bam
done

# 删除中间文件。
rm -f tmp.sam tmp.bam

# 修改名字,使得他们不再含有两个后缀。
mv bwa.sam.bam bwa.bammv bow.sam.bam bow.bam

# 对数据进行索引。
for name in *.bam
do
samtools index $name
done

# 删除这个程序生成的所有文件。
# rm -f *.bam *.bai *fq *.txt *.sam *.gff



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



本帖子中包含更多资源

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

x



上一篇:Biostar:课程15、16-samtools简介&利用readseq对GenBank文件格式转换
下一篇:Biostar:课程19、20-利用samtools mpileup和bcftools进行SNP calling
回复

使用道具 举报

0

主题

7

帖子

51

积分

注册会员

Rank: 2

积分
51
发表于 2017-10-13 13:12:04 | 显示全部楼层
从SAM文件的头文件(header)中提取results.ban的NC的大小是18959,但是根据cat NC.gff | awk '$3 =="gene" { len=$5-$4 + 1; size += len; print "Size:", size } '计算所有基因的累积长度是18127,不知道怎么是以哪个为准?

  
回复 支持 反对

使用道具 举报

0

主题

7

帖子

51

积分

注册会员

Rank: 2

积分
51
发表于 2017-10-13 14:58:47 | 显示全部楼层
命令ln -s ~/src/bowtie2-2.2.4/bowtie2-align ~/bin/;
可是在bowtie目录下,存在bowtie2-align-l和bowtie-aligh-s及*-debug文件,不知道这边软链接需要链接的是哪一个align文件?
回复 支持 反对

使用道具 举报

33

主题

46

帖子

230

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
230
 楼主| 发表于 2017-10-21 10:42:43 | 显示全部楼层
二口 发表于 2017-10-13 14:58
命令ln -s ~/src/bowtie2-2.2.4/bowtie2-align ~/bin/;
可是在bowtie目录下,存在bowtie2-align-l和bowtie ...

我个人觉得作者这里可能少了个通配符*     改为ln -s ~/src/bowtie2-2.2.4/bowtie2-align* ~/bin/可能好些。我自己安装的时候是直接在根目录的.bashrc(苹果系统是.profile)里修改的,加上这一句export PATH=$PATH:~/src/bowtie2-2.2.4
回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-9-20 19:25 , Processed in 0.035885 second(s), 26 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.