用awk进行简单的编程

  • A+
所属分类:Script

# 用awk进行简单的编程

# 我们之前已经做过这个了,但是为了让大家我们了解我们正在做什么,我们来再做一遍。

# 提取埃博拉基因组的编码特征、基因以及编码序列。

efetch -db nucleotide -id NC_002549.1 -format gb > NC.gb

readseq -format=GFF -o NC.gff NC.gb

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

#*curl -O ftp://ftp.arabidopsis.org/home/tair/Genes/TAIR10_genome_release/TAIR10_gff3/TAIR10_GFF3_genes.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文件中分离提取基因名字。

$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 } '

avatar

发表评论

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen: