awk编程实战「建议收藏」

2022-09-14 12:51:27 浏览数 (1)

大家好，又见面了，我是你们的朋友全栈君。

文章目录
  • 介绍
    • 模式pattern
    • 操作action
  • awk编程
    • 常用的内置变量
    • 变量赋值
    • BEGIN模块
    • END模块
    • 重定向和管道
    • 输出print与printf
    • 条件语句
    • 循环语句
    • 数组
    • 内建函数
  • 实战演练：awk分析拟南芥gff文件
    • 下载拟南芥gff文件
    • 查看一下gff格式是什么样子的
    • 查看第1-3列的数据
    • 每个特征序列长度是多少？
    • 用模式匹配来提取CDS特征
    • 计算所有基因的累积长度
    • 计算所有CDS的累积长度
    • 计算拟南芥（Col-0）基因组的大小
    • 根据特征（features）把文件分开
    • 提取启动子区域
  • 实战演练2
    • 针对特定列的计算，比如wig文件的标准化
    • 计算某列内容出现的次数
    • 数据矩阵的格式化输出
    • 判断FASTQ文件中，输出质量值的长度是与序列长度不一致的序列ID
    • 筛选差异基因
    • ID map，常用于转换序列的ID、提取信息、合并信息等
    • 字符串匹配
    • 字符串分割

介绍

awk是linux及unix操作系统中非常优秀的数据及文本处理工具，它是一种编程语言 awk命令格式为：

代码语言：javascript复制
awk pattern { 
   action} filename

相比于sed常常作用于一整行的处理，awk则比较倾向于将一行分成数个字段来处理。awk将输入数据视为一个文本数据库，像数据库一样，它也有记录和字段的概念。默认情况下，记录的分隔符是回车，字段的分隔符是空白符(空格,t)，所以输入数据的每一行表示一个记录，而每一行中的内容被空白分隔成多个字段。利用字段和记录，awk可以非常灵活地处理文件

测试文件

代码语言：javascript复制
[sunchengquan 13:40:26 ~/test]
$ last -5 |head -5 >tmp
[sunchengquan 13:41:01 ~/test]
$ cat tmp
suncheng pts/3        106.121.73.203   Thu Oct 25 12:34 - 12:49  (00:15)    
suncheng pts/2        106.121.73.203   Thu Oct 25 12:33   still logged in   
suncheng pts/2        106.121.73.203   Thu Oct 25 12:33 - 12:33  (00:00)    
suncheng pts/1        210.13.54.234    Thu Oct 25 11:02   still logged in   
suncheng pts/0        210.13.54.234    Thu Oct 25 10:27   still logged in

别名设置了默认的分隔符是t

awk脚本是由模式和操作组成的 两者是可选的，如果没有pattern，则action应用到全部记录，如果没有action，则输出匹配全部记录

awk {action} filename

代码语言：javascript复制
[sunchengquan 14:01:58 ~/test]
$ awk -F' ' '{print $3}' tmp
106.121.73.203
106.121.73.203
106.121.73.203
210.13.54.234
210.13.54.234

awk pattern {action} filename

代码语言：javascript复制
[sunchengquan 14:08:54 ~/test]
$ awk -F' ' '$3 ~ /^1/ {print $3}' tmp
106.121.73.203
106.121.73.203
106.121.73.203
[sunchengquan 14:10:49 ~/test]
$ awk -F' ' '$3=="210.13.54.234" {print $3}' tmp
210.13.54.234
210.13.54.234
[sunchengquan 14:22:50 ~/test]
$ awk -F' ' '/in/ {print $0}' tmp
suncheng pts/2        106.121.73.203   Thu Oct 25 12:33   still logged in   
suncheng pts/1        210.13.54.234    Thu Oct 25 11:02   still logged in   
suncheng pts/0        210.13.54.234    Thu Oct 25 10:27   still logged in

模式pattern

模式可以是以下任意一个：

  • /正则表达式/：使用通配符的扩展集。 awk -F' ' '/in/ {print $0}' tmp
  • 关系表达式：可以用下面运算符表中的关系运算符进行操作，可以是字符串或数字的比较，如$2>%1选择第二个字段比第一个字段长的行。 awk -F' ' '$3=="210.13.54.234" {print $3}' tmp
  • 模式匹配表达式：用运算符(匹配)和!(不匹配)。 awk -F' ' '$3 ~ /^1/ {print $3}' tmp
  • 模式，模式：指定一个行的范围。该语法不能包括BEGIN和END模式。
  • BEGIN：让用户指定在第一条输入记录被处理之前所发生的动作，通常可在这里设置全局变量。
  • END：让用户在最后一条输入记录被读取之后发生的动作。

操作action

操作由一个或多个命令、函数、表达式组成，之间由换行符或分号隔开，并位于大括号内，主要有四个部分：

  • 变量或数组赋值
  • 输出命令
  • 内置函数
  • 控制流命令

awk编程

典型的awk语法如下：

代码语言：javascript复制
awk '{ BEGIN{stat1} BEGIN{stat2} pattern1{action1} pattern2{action2} ... patternn{actionn} {默认动作,无条件,始终执行} END{stat1} END{stat2} }'

常用的内置变量

内置变量

解释

$0

当前所有字段

$1 – $n

当前第n个字段

FS

输入字段分隔符 默认是空格或t，相当于-F

RS

输入的记录分隔符， 默认为换行符 。awk ‘BEGIN{RS=”:”}{print $0}’ /etc/passwd

NF

当前处理行的列个数，字段个数。awk -F: ‘{print NF}’ /etc/passwd

NR 行号

awk ‘{print NR}’ /etc/passwd

FNR

当前记录数，与NR不同的是，这个值会是各个文件自己的行号

OFS

输出字段分隔符， 默认也是空格

ORS

输出的记录分隔符，默认为换行符

FILENAME

当前输入文件的名字

统计/etc/passwd:文件名，每行的行号，每行的列数，对应的完整行内容:

代码语言：javascript复制
[sunchengquan 14:23:27 ~/test]
$ awk  -F ':'  '{print "filename:" FILENAME ",linenumber:" NR ",columns:" NF ",linecontent:"$0}' /etc/passwd |head -3
filename:/etc/passwd,linenumber:1,columns:7,linecontent:root:x:0:0:root:/root:/bin/bash
filename:/etc/passwd,linenumber:2,columns:7,linecontent:bin:x:1:1:bin:/bin:/sbin/nologin
filename:/etc/passwd,linenumber:3,columns:7,linecontent:daemon:x:2:2:daemon:/sbin:/sbin/nologin

变量赋值

Variable = expression 变量不需要定义就可以直接使用，变量类型可以是数字或字符串

代码语言：javascript复制
[sunchengquan 14:51:25 ~/test]
$ awk '$1 ~/suncheng/{count = "SCQ"; print count}' tmp 
SCQ
SCQ
SCQ
SCQ
SCQ

域变量也可被赋值和修改

代码语言：javascript复制
[sunchengquan 14:54:44 ~/test]
$ awk -F' ' '$1 ~/suncheng/{ 
   $2=100 $6 ; print $0}' tmp 
suncheng    125    106.121.73.203    Thu    Oct    25    12:34    -    12:49    (00:15)

自定义及外部变量

代码语言：javascript复制
awk -F: '$7 ~ /^/bin/{print $0}' /etc/passwd
awk -F: -v reg='^/bin.*' '$7 ~ reg {print $0}' /etc/passwd

awk -v n=$HOSTNAME '{print n}' /etc/passwd
echo $HOSTNAME   输出hostname
awk -v n=$HOSTNAME 'BEGIN{print n}'

BEGIN模块

BEGIN模块后紧跟着动作块，这个动作块在awk处理任何输入文件之前执行。所以它可以在没有任何输入的情况下进行测试。它通常用来改变内建变量的值，如OFS,RS和FS等，以及打印标题。

代码语言：javascript复制
[sunchengquan 14:55:51 ~/test]
$ awk 'BEGIN{FS=" ";OFS="t";ORS="nn" ;print "start it "}{print $1,$3}' tmp
start it 

suncheng    106.121.73.203

suncheng    106.121.73.203

suncheng    106.121.73.203

suncheng    210.13.54.234

suncheng    210.13.54.234

END模块

END不匹配任何的输入文件，但是执行动作块中的所有动作，它在整个输入文件处理完成后被执行。

代码语言：javascript复制
[sunchengquan 15:12:17 ~/test]
$ awk 'END{print "The number of records is " NR}' tmp
The number of records is 5
上式将打印所有被处理的记录数。

重定向和管道

awk可使用shell的重定向符进行重定向输出

代码语言：javascript复制
[sunchengquan 15:12:32 ~/test]
$ awk '$1=100 {print $1 > "output_file"}' tmp
[sunchengquan 15:15:27 ~/test]
$ cat output_file |head -2
100
100

输出重定向需用到getline函数 getline从标准输入、管道或者当前正在处理的文件之外的其他输入文件获得输入。

代码语言：javascript复制
[sunchengquan 15:15:46 ~/test]
$ awk 'BEGIN{ "date" | getline d; print d}' 
2018年 10月 25日 星期四 15:17:26 CST
[sunchengquan 15:17:26 ~/test]
$ awk 'BEGIN{ "date" | getline d; print d}' tmp
2018年 10月 25日 星期四 15:17:42 CST
[sunchengquan 15:18:10 ~/test]
$ awk 'BEGIN{ "date" | getline d; print d} {print $0}' tmp
2018年 10月 25日 星期四 15:18:30 CST
suncheng pts/3        106.121.73.203   Thu Oct 25 12:34 - 12:49  (00:15)    
suncheng pts/2        106.121.73.203   Thu Oct 25 12:33   still logged in   
suncheng pts/2        106.121.73.203   Thu Oct 25 12:33 - 12:33  (00:00)    
suncheng pts/1        210.13.54.234    Thu Oct 25 11:02   still logged in   
suncheng pts/0        210.13.54.234    Thu Oct 25 10:27   still logged in

执行shell的date命令，并通过管道输出给getline，然后getline从管道中读取并将输入赋值给d，split函数把变量d转化成数组mon，然后打印数组mon的第二个元素

代码语言：javascript复制
[sunchengquan 15:21:09 ~/test]
$ awk 'BEGIN{"date" | getline d; split(d,mon," "); print mon[2]}'
10月

输出print与printf

代码语言：javascript复制
[sunchengquan 21:07:38 ~]
$ awk -F: '{print $1 ":" $2}' /etc/passwd|tail -5
mysql:x
tomcat:x
apache:x
mongod:x
avahi:x
代码语言：javascript复制
[sunchengquan 21:07:47 ~]
$ awk -F: '{print "hello "$1}' /etc/passwd | tail -5
hello mysql
hello tomcat
hello apache
hello mongod
hello avahi

printf格式化输出，默认不换行

代码语言：javascript复制
[sunchengquan 21:08:23 ~]
$ awk -F: '{printf("hello %sn", $1)}' /etc/passwd |tail -5
hello mysql
hello tomcat
hello apache
hello mongod
hello avahi

条件语句

awk中的条件语句是从C语言中借鉴过来的，可控制程序的流程。

if 语句

代码语言：javascript复制
格式：
        { 
   if (expression){ 
   
                   statement; statement; ...
                     }
        }

awk -F: '{if(NF=="/bin/bash" &&

if/else语句，用于双重判断

代码语言：javascript复制
格式：
        { 
   if (expression){ 
   
                   statement; statement; ...
                       }
        else { 
   
                   statement; statement; ...
                       }
        }
代码语言：javascript复制
seq 10 | awk '{if($0%2==0){print "OK"}else{print "no"}}’ 对，标准写法，最好不要省略{} seq 10 | awk '{ 
   if($0%2==0) print "OK";else print "no”}’ 对 seq 10 |awk '{if ($0%2==0){ count  ; print "Y"} else {count--; print "N"}}’ 对 seq 10 |awk '{if ($0%2==0) count  ; print "Y"; else count--; print "N”}’    错
awk: cmd. line:1: { 
   if ($0%2==0) count  ; print "Y"; else count--; print "N"}
awk: cmd. line:1:                                   ^ syntax error

if/else else if语句，用于多重判断

代码语言：javascript复制
格式：
        { 
   if (expression){ 
   
                    statement; statement; ...
                   }
        else if (expression){ 
   
                    statement; statement; ...
                   }
        else if (expression){ 
   
                    statement; statement; ...
                   }
        else { 
   
                   statement; statement; ...
             }
        }

循环语句

常用的两种循环：while，for

代码语言：javascript复制
[sunchengquan 15:48:13 ~/test]
$ awk 'BEGIN{FS=" "}{i=1; while ( i <= NF ) { print NF,$i; i  }}' tmp
[sunchengquan 15:50:08 ~/test]
$ awk 'BEGIN{FS=" "}{for(i=1; i <= NF;i   ) { print NF,$i}}' tmp

break：用于在满足条件的情况下跳出循环 continue：用于在满足条件的情况下忽略后面的语句 next exit

数组

定义或添加数组元素 array[1]="hello" array["name"]=jake 解释：1和name为索引，hello和jake为键值 awk的数组中的索引可以是数字或字符串

代码语言：javascript复制
awk BEGIN'{a[5]="Jack";a["name"]="lilei";print a[5],a["name"]}’ awk 'BEGIN{ 
   a[5]="Jack";a["name"]="lilei";print a[5],a["name"]}' awk BEGIN'{ 
   a[5]="Jack";a["name"]="lilei";for(i in a){ 
   print a[i]}}' awk BEGIN'{ 
   a[5]="Jack";a["name"]="lilei";for(i in a){ 
   print i}}' awk BEGIN'{ 
   a[5]="Jack";a["name"]="lilei";for(i in a){ 
   print i":"a[i]}}’

删除数组元素 delete array['"name"]

代码语言：javascript复制
[sunchengquan 15:55:34 ~/test]
$ awk 'BEGIN{a[5]="Jack";a["name"]="lilei";delete a["name"];for(i in a){print i":"a[i]}}'   
5:Jack

内建函数

算术函数

算术函数

解释

int(x)

返回x的整数部分的值

sqrt(x)

返回x的平方根

rand()

返回伪随机数r，其中0<=r<1

srand(x)

建立rand()新的种子数。如果没有指定就用当天的时间

代码语言：javascript复制
[sunchengquan 21:18:23 ~]
$ awk BEGIN'{print rand()}'
0.237788
[sunchengquan 21:20:37 ~]
$ awk BEGIN'{print rand();srand();print rand()}'
0.237788
0.774373
[sunchengquan 21:20:44 ~]
$ awk 'BEGIN{print rand();srand();print rand()}'
0.237788
0.483444

字符串函数 sub, gsub() 替换函数 gsub全局替换 sub (regular expression, substitution string, target string)

代码语言：javascript复制
[sunchengquan 21:21:12 ~]
$ echo "hello world" |awk '{sub("world","sunchengquan");print $0} '
hello sunchengquan
[sunchengquan 21:26:25 ~]
$ echo "hello world world" |awk '{sub("world","sunchengquan");print $0} '
hello sunchengquan world
[sunchengquan 21:26:37 ~]
$ echo "hello world world" |awk '{gsub("world","sunchengquan");print $0} '
hello sunchengquan sunchengquan

index(s,t) 返回子串t在字符串s中的位置，如果没有则返回0

代码语言：javascript复制
[sunchengquan 21:27:20 ~]
$ echo "helloworld" |awk '{print index($0,"world")} '
6
[sunchengquan 21:27:24 ~]
$ echo "hello world" |awk '{print index($0,"world")} '
7

length(s) 返回字符串长度，当没有给出s时，返回$0 的长度 match(s,r) 如果正则表达式r在s中匹配到，则返回出现的起始位置，否则返回0 split(s,a,sep)使用sep将字符串s分解到数组a中。默认sep为FS。 split( string, array, field separator )

代码语言：javascript复制
[sunchengquan 21:27:37 ~]
$ echo "00-11-22-33-44"|awk '{split($0,a,sep="-");for(i in a){print i":"a[i]}}'
4:33
5:44
1:00
2:11
3:22

tolower(s) 将字符串s中的所有大写字符转换为小写

代码语言：javascript复制
[sunchengquan 21:28:28 ~]
$ echo "HELLO"|awk '{print tolower($0)}'
hello

toupper(s) 将字符串s中的所有小写字符转换为大写

代码语言：javascript复制
[sunchengquan 21:28:34 ~]
$ echo "hello"|awk '{print toupper($0)}'
HELLO

substr(s,p) 返回字符串s中从p开始的后缀部分 substr(s,p,n) 返回字符串s中从p开始长度为n的后缀部分 substr(3,10,8) —> 表示是从第3个字段里的第10个字符开始，截取8个字符结束. substr(

代码语言：javascript复制
[sunchengquan 17:42:58 ~/local/app/bowtie2-2.2.4/example/reads]
$ echo "hello" | awk '{print substr($0,1,2)}'
he
[sunchengquan 17:45:46 ~/local/app/bowtie2-2.2.4/example/reads]
$ echo "hello" | awk '{print substr($0,2)}'
ello

自定义函数

代码语言：javascript复制
[sunchengquan 21:09:24 ~]
$ awk 'function sum(n,m){total=n m;return total} BEGIN{print sum(1,2)}'
3

实战演练：awk分析拟南芥gff文件

下载拟南芥gff文件

代码语言：javascript复制
curl -O ftp://ftp.arabidopsis.org/home/tair/Genes/TAIR10_genome_release/TAIR10_gff3/TAIR10_GFF3_genes.gff

查看一下gff格式是什么样子的

代码语言：javascript复制
less -S TAIR10_GFF3_genes.gff |head
Chr1	TAIR10	chromosome	1	30427671	.	.	.	ID=Chr1;Name=Chr1
Chr1	TAIR10	gene	3631	5899	.	 	.	ID=AT1G01010;Note=protein_coding_gene;Name=AT1G01010
Chr1	TAIR10	mRNA	3631	5899	.	 	.	ID=AT1G01010.1;Parent=AT1G01010;Name=AT1G01010.1;Index=1
Chr1	TAIR10	protein	3760	5630	.	 	.	ID=AT1G01010.1-Protein;Name=AT1G01010.1;Derives_from=AT1G01010.1
Chr1	TAIR10	exon	3631	3913	.	 	.	Parent=AT1G01010.1
Chr1	TAIR10	five_prime_UTR	3631	3759	.	 	.	Parent=AT1G01010.1
Chr1	TAIR10	CDS	3760	3913	.	 	0	Parent=AT1G01010.1,AT1G01010.1-Protein;
Chr1	TAIR10	exon	3996	4276	.	 	.	Parent=AT1G01010.1
Chr1	TAIR10	CDS	3996	4276	.	 	2	Parent=AT1G01010.1,AT1G01010.1-Protein;
Chr1	TAIR10	exon	4486	4605	.	 	.	Parent=AT1G01010.1
gff文件是tab分隔的文件
第1列是染色体信息
第2列是gff注释数据来源
第3列为特征（feature）即属于gene还是mRNA还是CDS等等
第4和5列分别是这个特征序列的起始和终止位置
第6列是得分，可以是序列相似性比对时的E-values值或者基因预测是的P-values值, ”.”表示为空
第7列是表示序列的方向：正义链为 ，反义链为-
第8列仅为对CDS的注释，表示起始编码的位置，有效值为0、1、2
第9列为注释信息

查看第1-3列的数据

1, 2,

代码语言：javascript复制
$ cat TAIR10_GFF3_genes.gff | awk ' { print $1, $2, $3 } ' | head -5
Chr1	TAIR10	chromosome
Chr1	TAIR10	gene
Chr1	TAIR10	mRNA
Chr1	TAIR10	protein
Chr1	TAIR10	exon

这（基本上）等同于截取列
$ cat TAIR10_GFF3_genes.gff | cut -f 1,2,3 | head -5
Chr1	TAIR10	chromosome
Chr1	TAIR10	gene
Chr1	TAIR10	mRNA
Chr1	TAIR10	protein
Chr1	TAIR10	exon

每个特征序列长度是多少？

第5列的数值减去第4列的数值后 1，即得到特征序列的长度

代码语言：javascript复制
$ cat TAIR10_GFF3_genes.gff | awk ' { print $3, $5-$4   1 } ' | head -8
chromosome	30427671
gene	2269
mRNA	2269
protein	1871
exon	283
five_prime_UTR	129
CDS	154
exon	281

用模式匹配来提取CDS特征

代码语言：javascript复制
cat TAIR10_GFF3_genes.gff | awk '$3 =="gene" { print $3, $5-$4   1, $9 } '| head -5
cat TAIR10_GFF3_genes.gff | awk '{if($3 =="gene") { print $3, $5-$4   1, $9 }} '| head -5

gene	2269	ID=AT1G01010;Note=protein_coding_gene;Name=AT1G01010
gene	2810	ID=AT1G01020;Note=protein_coding_gene;Name=AT1G01020
gene	2066	ID=AT1G01030;Note=protein_coding_gene;Name=AT1G01030
gene	8082	ID=AT1G01040;Note=protein_coding_gene;Name=AT1G01040
gene	207	ID=AT1G01046;Note=miRNA;Name=AT1G01046

计算所有基因的累积长度

代码语言：javascript复制
cat TAIR10_GFF3_genes.gff  | awk '$3 =="gene" { len=$5-$4   1; size  = len; print "Size:", size } '

.....................
Size: 61219702
Size: 61220038
Size: 61220156
Size: 61222091
Size: 61222409
Size: 61223024

计算所有CDS的累积长度

代码语言：javascript复制
cat TAIR10_GFF3_genes.gff | awk '$3 =="CDS" { len=$5-$4   1; size  = len; print "Size:", size } '

............................
Size: 43543450
Size: 43543888
Size: 43545472
Size: 43545808
Size: 43546126
Size: 43546741

计算拟南芥（Col-0）基因组的大小

119667750是拟南芥（Col-0）基因组的大小 可以用下边的代码自行计算

代码语言：javascript复制
$ cat TAIR10_GFF3_genes.gff |awk '{if($3 == "chromosome"){len=$5-$4   1; size  = len; print "Size:", size,$9 } }'
Size:	30427671	ID=Chr1;Name=Chr1
Size:	50125960	ID=Chr2;Name=Chr2
Size:	73585790	ID=Chr3;Name=Chr3
Size:	92170846	ID=Chr4;Name=Chr4
Size:	119146348	ID=Chr5;Name=Chr5
Size:	119300826	ID=ChrC;Name=ChrC
Size:	119667750	ID=ChrM;Name=ChrM

根据特征（features）把文件分开

代码语言：javascript复制
cat NC.gff | awk ' $3=="gene" { print $0 }' >> NC-genes.gff
less -S NC-genes.gff|head -5

##gff-version 2
NC_002549    -    gene    56    3026    .         .    gene "NP" ; locus_tag "ZEBOVgp1" ; db_xref "GeneID:911830"
NC_002549    -    gene    3032    4407    .         .    gene "VP35" ; locus_tag "ZEBOVgp2" ; db_xref "GeneID:911827"
NC_002549    -    gene    4390    5894    .         .    gene "VP40" ; locus_tag "ZEBOVgp3" ; db_xref "GeneID:911825"
NC_002549    -    gene    5900    8305    .         .    gene "GP" ; locus_tag "ZEBOVgp4" ; db_xref "GeneID:911829"

cat NC.gff | awk ' $3=="CDS" { print $0 }'  >> NC-cds.gff
less -S NC-cds.gff|head -5

提取启动子区域

代码语言：javascript复制
cat TAIR10_GFF3_genes.gff| awk 'BEGIN{OFS=FS="t"}{if($3=="gene"){if($7==" "){start=$4-1;up=start-1000;if(up<0) up=0;dw=start 500;print $1,up,dw,$7;} else if($7=="-"){start=$5-1; up=start 1000; dw=start-500; if(dw<0) dw=0; print $1,dw,up,$7}}}' >TAIR10.promoter.bed

实战演练2

针对特定列的计算，比如wig文件的标准化

代码语言：javascript复制
$ cat ehbio.wig 
variableStep chrom=chr2
300701	12.5
300702	12.5
300703	12.5
300704	12.5
300705	12.0


$ awk 'BEGIN{OFS=FS="t"}{ 
   $2=$2*10^6/(2.5*10^6); print $0}' ehbio.wig
variableStep chrom=chr2	0
300701	5
300702	5
300703	5
300704	5
300705	4.8

计算某列内容出现的次数

代码语言：javascript复制
$ cat count 
ID	Type
Pou5f1	Pluripotency
Nanog	Pluripotency
Sox2	Neuron
Tet1	Epigenetic
Tet3	Epigenetic
Myc	Oncogene

$ awk 'BEGIN{OFS=FS="t"}{if(FNR>1) a[$2] =1;}END{print "CounttType"; for(i in a) print a[i],i;}' count
Count	Type
2	Pluripotency
1	Oncogene
1	Neuron
2	Epigenetic

数据矩阵的格式化输出

代码语言：javascript复制
 cat numertic.matrix 
ID	A	B	C
a	1.002	1.234	1.999
b	2.333	4.232	0.889

$ awk '{if(FNR==1) print $0;else {printf "%s%s",$1,FS;for (i=2;i<=NF;i  ) printf "%.1f%s",$i, (i==NF?RS:FS)}}' numertic.matrix 
ID	A	B	C
a 1.0 1.2 2.0
b 2.3 4.2 0.9

判断FASTQ文件中，输出质量值的长度是与序列长度不一致的序列ID

代码语言：javascript复制
zcat Test_2.fq.gz | awk '{if(FNR%4==1) ID=$0; else if(FNR%4==2) seq_len=length($0); else if(FNR%4==0) {quality_len=length($0); if(seq_len!=quality_len) print ID; }}'

筛选差异基因

代码语言：javascript复制
 cat de_gene
ID	log2fc	padj
A	1	0.001
B	-1	0.001
C	1	0.001
D	2	0.0001
E	-0.51	0.051
F	0.1	0.1
G	1	0.1

awk '$3<0.05 || NR==1' de_gene 
ID	log2fc	padj
A	1	0.001
B	-1	0.001
C	1	0.001
D	2	0.0001


awk 'BEGIN{OFS=FS="t"}{if(FNR==1) print $0;else{abs_log2fc=($2<0?$2*(-1):$2);if(abs_log2fc>=1 && $3<0.05) print $0;}}' de_gene 
ID	log2fc	padj
A	1	0.001
B	-1	0.001
C	1	0.001
D	2	0.0001


存储到不同的文件
awk 'BEGIN{OFS=FS="t"; up="up"; dw="dw";}{if(FNR==1) {print $0 >up; print $0 >dw;} else if ($3<0.05) {if($2>=1) print $0>up; else if($2<=-1) print $0>dw}}' de_gene 

$ head up dw
==> up <==
ID	log2fc	padj
A	1	0.001
C	1	0.001
D	2	0.0001

==> dw <==
ID	log2fc	padj
B	-1	0.001

ID map，常用于转换序列的ID、提取信息、合并信息等

代码语言：javascript复制
$ cat id_map 
ENSM	Symbol	Entrez
ENSG00000280516	TMEM42	693149
ENSG00000281886	TGM4	7047
ENSG00000280873	DGKD	8527
ENSG00000281244	ADAMTS13	11093
ENSG00000280701	RP11-272D20.2	
ENSG00000280674	ZDHHC3	51304
ENSG00000281623	Y_RNA	
ENSG00000280479	CACFD1	11094
ENSG00000281165	SLC2A6	11182
ENSG00000281879	ABO	28
ENSG00000282873	BCL7A	605
ENSG00000280651	AC156455.1	100506691
[sunchengquan 23:19:33 ~]
$ cat ensm 
ENSG00000281244
ENSG00000281165
ENSG00000282873

$ awk 'BEGIN{OFS=FS="t"} ARGIND==1{if(FNR>1) {map[$1]=$3};} ARGIND==2{print map[$1];}' id_map ensm 
11093
11182
605


转换大小写, toupper, tolower
$ cat symbol 
Tgm4
Dgkd
Abo
$ awk '{print toupper($1)}' symbol 
TGM4
DGKD
ABO

$ awk 'BEGIN{OFS=FS="t"}ARGIND==1{if(FNR>1) ensm2entrez[$2]=$3;}ARGIND==2{print ensm2entrez[toupper($1)];}' id_map symbol 
7047
8527
28

字符串匹配

代码语言：javascript复制
$ cat ens.bed
1	100	105
2	100	105
3	100	105
Mt	100	105
X	100	105
[sunchengquan 16:11:25 ~]
$ awk 'BEGIN{OFS=FS="t"}{if($1~/^[0-9XY]/) $1="chr"$1; else if($1~/M.*/)gsub(/M.*/, "chrM", $1); print $0}' ens.bed
chr1	100	105
chr2	100	105
chr3	100	105
chrM	100	105
chrX	100	105

字符串分割

代码语言：javascript复制
$ cat trinity_id 
Trinity_C1_g1_i1
Trinity_C1_g1_i2
Trinity_C1_g1_i3
Trinity_C2_g1_i1
Trinity_C3_g1_i1
Trinity_C3_g3_i2
[sunchengquan 16:40:11 ~]
$  awk 'BEGIN{OFS=FS="t"}{count=split($1, geneL, "_"); gene=geneL[1];print gene}' trinity_id
Trinity
Trinity
Trinity
Trinity
Trinity
Trinity
[sunchengquan 16:40:24 ~]
$  awk 'BEGIN{OFS=FS="t"}{count=split($1, geneL, "_"); gene=geneL[1];print geneL[1],geneL[2];}' trinity_id
Trinity	C1
Trinity	C1
Trinity	C1
Trinity	C2
Trinity	C3
Trinity	C3
[sunchengquan 16:41:27 ~]
$ awk 'BEGIN{OFS=FS="t"}{count=split($1, geneL, "_"); gene=geneL[1]; for(i=2;i<count;i  ) gene=gene"_"geneL[i]; print gene,$1;}' trinity_id 
Trinity_C1_g1	Trinity_C1_g1_i1
Trinity_C1_g1	Trinity_C1_g1_i2
Trinity_C1_g1	Trinity_C1_g1_i3
Trinity_C2_g1	Trinity_C2_g1_i1
Trinity_C3_g1	Trinity_C3_g1_i1
Trinity_C3_g3	Trinity_C3_g3_i2
[sunchengquan 16:47:16 ~]
$ awk 'BEGIN{OFS=FS="t"}{count=split($1, geneL, "_"); gene=geneL[1]; for(i=2;i<=count;i  ) gene=gene"_"geneL[i]; print gene,$1;}' trinity_id 
Trinity_C1_g1_i1	Trinity_C1_g1_i1
Trinity_C1_g1_i2	Trinity_C1_g1_i2
Trinity_C1_g1_i3	Trinity_C1_g1_i3
Trinity_C2_g1_i1	Trinity_C2_g1_i1
Trinity_C3_g1_i1	Trinity_C3_g1_i1
Trinity_C3_g3_i2	Trinity_C3_g3_i2

发布者：全栈程序员栈长，转载请注明出处：https://javaforall.cn/159733.html原文链接：https://javaforall.cn

linux unix 编程算法

0 人点赞