文本处理是完成生信项目的基本功,能实现的方法绝对不止一种,最先想到的肯定是你比较熟悉的,这个时候谁能抵过单行命令的诱惑呢?
awk是由Alfred Aho、Peter Weinberger 和Brain Kernighan这3个人创造的,awk是这三个人的姓氏的首字母:Aho、Weinberger、Kernighan缩写而成。
awk最早是在unix上实现的,因此现在linux上用的是gawk,及GNU awk,简称gawk。
说到awk我们就不得不提linux三剑客中的另外两者:
grep: 查找和匹配文本(类似R语言中的grep
、grepl
、str_detect
和excel的查找
功能)
sed:编辑匹配到的文本(类似R中的str_replace
、sub
、gsub
和excel的查找和替换
)
awk:格式化文本,对文本进行复杂处理
由此可见awk绝对算是三剑客中的大拇指!
awk基础之action
awk基本语法:
awk [options] 'program' file1, file2
其中program部分又可分为pattern和action,即awk基本语法:
awk [options] 'Pattern{action}' file
字面意思理解,action就是动作,awk擅长文本格式化,并将格式化后的文本输出,所以最常用的动作是print
和printf
。
先不用[options]和pattern,直接用常用的action:print
,尝试一下:
$ echo abc > tested
$ awk '{print}' tested
abc
上面只用awk执行了print动作,将tested打印了出来。
类似的情景:
上述内容仍然没有用[options]和pattern
df | awk '{ print $5} '
$5:表示当前行按照分割符分割后的第5列,不指定分割符时,默认使用空格为分割符。
其实上面信息的有多个地方有连续空格,awk是怎么处理的呢?
这就是awk牛×的地方所在,自动将连续的空格理解为一个分割符。
awk是逐行处理文本的,即:当awk处理文本时,会一行一行处理,处理完当前行,再处理下一行,awk默认以“换行符”为标记,识别每一行。到这里是不是感觉有点意思?这TMD就是我们人类识别文本的流程嘛。当然,awk还会按照用户指定的分割符去分割当前行,默认空格是分割符。
$0
:显示整行$NF
:当前行分割后的最后一列
注意:
$0
和$NF
均为内置变量,但$NF
和NF
要表达的意思不一样,NF表示当前行被分割符切开后,一共有几个字段。
假如一行文本被空格分为7段,那么NF = 7,$NF
的值就是$7
,而$7
表示当前行的第7个字段,即最后一列,那么倒数第二列可以写成$(NF-1)
。
返回多列:awk '{print $1, $2}' tested
, 中间用
$ cat tested
adc 123 iuy ddd
8ua 456 auv ppp 7y7
$ awk '{print $1,$2}' tested
adc 123
8ua 456
$ awk '{print $1,$2,$5}' tested
adc 123
8ua 456 7y7
可以看出输出多列时,用逗号隔开要输出的列即可。此外,还能添加自定义的字段,将给定的字段与文件列结合起来:
#在文本后方插入
$ awk '{print $1,$2, "string"}' tested
adc 123 string
8ua 456 string
#文本前方插入
$ awk '{print "xiaoming:"$1,"xiaoqiang:"$2}' tested
xiaoming:adc xiaoqiang:123
xiaoming:8ua xiaoqiang:456
#任何位置
$ awk '{print "xiaoming:"$1,666,"xiaoqiang:"$2}' tested
xiaoming:adc 666 xiaoqiang:123
xiaoming:8ua 666 xiaoqiang:456
注意:
$1
外侧不能加人双引号,否则$1
会被当做文本处理。例如:
$ awk '{print $1}' tested
adc
8ua
$ awk '{print "$1"}' tested
$1
awk基础之pattern
先来学习下AWK的两种特殊模式:BEGIN
和END
。
BEGIN
:处理文本之前需执行的操作
END
:处理完术后行之后所需要执行的操作
例子:
$ awk 'BEGIN{print "aaa","bbb"} {print $1,$2}' tested
aaa bbb
adc 123
8ua 456
awk '{print $1,$2} END{print "ccc", "ddd"}' tested
adc 123
8ua 456
ccc ddd
awk 'BEGIN{print "aaa","bbb"}{print $1,$2} END{print "ccc", "ddd"}' tested
aaa bbb
adc 123
8ua 456
ccc ddd
可见返回的结果是不是有点像一张“报表”,有“表头”、“表内容”、“表尾”,足可见awk对文本的格式化能力有多强。
实际操作
下载测试数据
wget https://raw.github.com/nachocab/nachocab.github.io/master/assets/transcriptome.gtf
刘小泽金句:任何时候,一旦有放弃的念头,就要想想,世界上这么多做相似的事情,为什么他们没有放弃?是不是有更好的解决办法?试着解决而不是放弃,这样你会发现越来越有趣。
less -SN filename
强大的预览功能,竟然还有行号
$ less -SN transcriptome.gtf
1 ##description: evidence-based annotation of the human genome (GRCh37), version 18 (Ensembl 73)
2 ##provider: GENCODE
3 ##contact: gencode@sanger.ac.uk
4 ##format: gtf
5 ##date: 2013-09-02
6 chr1 HAVANA exon 173753 173862 . - . gene_id "ENSG00000241860.2"; transcript_id "ENST00000466557.2"; gene_type "processed_transcript"; gene_stat>
7 chr1 HAVANA transcript 1246986 1250550 . - . gene_id "ENSG00000127054.14"; transcript_id "ENST00000478641.1"; gene_type "protein_coding"; gene_s>
8 chr1 HAVANA CDS 1461841 1461911 . + 0 gene_id "ENSG00000197785.9"; transcript_id "ENST00000378755.5"; gene_type "protein_coding"; gene_status "KN>
9 chr1 HAVANA exon 1693391 1693474 . - . gene_id "ENSG00000008130.11"; transcript_id "ENST00000341991.3"; gene_type "protein_coding"; gene_status "K>
10 chr1 HAVANA CDS 1688280 1688321 . - 0 gene_id "ENSG00000008130.11"; transcript_id "ENST00000497186.1"; gene_type "protein_coding"; gene_status "K>
11 chr1 ENSEMBL exon 2125078 2125349 . - . gene_id "ENSG00000162585.12"; transcript_id "ENST00000400919.3"; gene_type "protein_coding"; gene_status "K>
12 chr1 HAVANA exon 2334474 2334503 . + . gene_id "ENSG00000157916.14"; transcript_id "ENST00000306256.9"; gene_type "protein_coding"; gene_status "K>
13 chr1 HAVANA stop_codon 2453049 2453051 . - 0 gene_id "ENSG00000157881.9"; transcript_id "ENST00000502770.1"; gene_type "protein_coding"; gene_st>
【↑ 和 ↓
:上下翻页;shift + g
跳到底部;gg
返回顶部】
这个GTF文件共有9列,分别是:
chromosome, annotation source, feature type, start, end, score, strand, and phase
最后一列phase非常长,但很重要,它是以键值对的形式给出的,同一键值对的键与值用空格分隔,不同键值对用分号隔开。【这些信息都可能为以后准确提取大量文本提供帮助】
Q:提取只有一个外显子(exon)的蛋白编码基因?
先对第3列进行筛选
$ awk '$3=="gene"' transcriptome.gtf | head |less -SN
1 chr1 HAVANA gene 32826871 32829879 . - . gene_id "ENSG00000225828.1"; transcript_id "ENSG00000225828.1"; gene_type "protein_coding";>
2 chr1 HAVANA gene 43241484 43242461 . - . gene_id "ENSG00000228452.1"; transcript_id "ENSG00000228452.1"; gene_type "antisense"; gene>
3 chr1 HAVANA gene 109512836 109584850 . - . gene_id "ENSG00000085433.11"; transcript_id "ENSG00000085433.11"; gene_type "protein_coding>
4 chr1 HAVANA gene 155121020 155121295 . + . gene_id "ENSG00000223452.2"; transcript_id "ENSG00000223452.2"; gene_type "pseudogene"; gen>
5 chr1 HAVANA gene 179559507 179560440 . + . gene_id "ENSG00000261060.1"; transcript_id "ENSG00000261060.1"; gene_type "sense_overlappin>
6 chr1 HAVANA gene 237205505 237997288 . + . gene_id "ENSG00000198626.10"; transcript_id "ENSG00000198626.10"; gene_type "protein_coding>
7 chr2 HAVANA gene 24123563 24124019 . + . gene_id "ENSG00000238111.1"; transcript_id "ENSG00000238111.1"; gene_type "pseudogene"; gen>
8 chr2 ENSEMBL gene 29726773 29726859 . - . gene_id "ENSG00000266310.1"; transcript_id "ENSG00000266310.1"; gene_type "miRNA"; gene_sta>
9 chr2 ENSEMBL gene 122360649 122360740 . - . gene_id "ENSG00000222467.1"; transcript_id "ENSG00000222467.1"; gene_type "misc_RNA"; gene_>
10 chr2 ENSEMBL gene 203156055 203156144 . + . gene_id "ENSG00000271852.1"; transcript_id "ENSG00000271852.1"; gene_type "snoRNA"; gene_st>
上面提到过第9列phase中键和键值之间是用空格分隔的,虽然列之间默认只用空格分隔的,但是列内没有默认,因此需要制定分隔符【options: -F
】Field。
$ awk -F "\t" '$3=="gene" {print $9}' transcriptome.gtf |head |less -S
gene_id "ENSG00000225828.1"; transcript_id "ENSG00000225828.1"; gene_type "protein_coding"; gene_status "KNOWN"; gene_name "FAM229A"; transcript_type "protein_coding"; transcript_>
gene_id "ENSG00000228452.1"; transcript_id "ENSG00000228452.1"; gene_type "antisense"; gene_status "NOVEL"; gene_name "RP5-994D16.9"; transcript_type "antisense"; transcript_statu>
gene_id "ENSG00000085433.11"; transcript_id "ENSG00000085433.11"; gene_type "protein_coding"; gene_status "KNOWN"; gene_name "WDR47"; transcript_type "protein_coding"; transcript_>
gene_id "ENSG00000223452.2"; transcript_id "ENSG00000223452.2"; gene_type "pseudogene"; gene_status "KNOWN"; gene_name "HMGN2P18"; transcript_type "pseudogene"; transcript_status >
gene_id "ENSG00000261060.1"; transcript_id "ENSG00000261060.1"; gene_type "sense_overlapping"; gene_status "NOVEL"; gene_name "RP11-545A16.4"; transcript_type "sense_overlapping";>
gene_id "ENSG00000198626.10"; transcript_id "ENSG00000198626.10"; gene_type "protein_coding"; gene_status "KNOWN"; gene_name "RYR2"; transcript_type "protein_coding"; transcript_s>
gene_id "ENSG00000238111.1"; transcript_id "ENSG00000238111.1"; gene_type "pseudogene"; gene_status "KNOWN"; gene_name "AC066692.3"; transcript_type "pseudogene"; transcript_statu>
gene_id "ENSG00000266310.1"; transcript_id "ENSG00000266310.1"; gene_type "miRNA"; gene_status "NOVEL"; gene_name "AC106899.1"; transcript_type "miRNA"; transcript_status "NOVEL";>
gene_id "ENSG00000222467.1"; transcript_id "ENSG00000222467.1"; gene_type "misc_RNA"; gene_status "NOVEL"; gene_name "Y_RNA"; transcript_type "misc_RNA"; transcript_status "NOVEL">
gene_id "ENSG00000271852.1"; transcript_id "ENSG00000271852.1"; gene_type "snoRNA"; gene_status "KNOWN"; gene_name "SNORD11B"; transcript_type "snoRNA"; transcript_status "KNOWN";>
番外
less常用命令选项
选项 | 功能 |
---|---|
-N | 显示每行的行号 |
-S | 行过长时将超出部分舍弃 |
-e | 当文件显示结束后,自动离开 |
-i | 忽略搜索时的大小写 |
-m | 显示类似 more 命令的百分比 |
常用交互指令
交互指令 | 功能 |
---|---|
b | 向上移动一页 |
d | 向下移动半页 |
q 或 Q | 退出 less 命令 |
空格键 | 向下移动一页 |
回车键 | 向下移动一行 |
Ctrl+g | 到页面底部 |
gg | 回到页面顶端 |
推荐一本学习awk的书籍:
由AWK语言作者亲自主笔编写,非常贴近实际应用
说明:本文为下面参考链接的学习实践笔记,如有侵权,请联系删除。
参考链接:
生信单行命令的美,慢慢体会(一)
awk从放弃到入门(1):awk基础 (通俗易懂,快进来看)
【笔记】less命令 - 查看文件内容