专栏首页生信宝典Linux学习 - 常用和不太常用的实用awk命令

Linux学习 - 常用和不太常用的实用awk命令

Linux学习系列文章是生信宝典最开始主推的一块,力图从一个新额视角帮助初学者快速入门Linux系统,熟悉Linux下的文件和目录文件操作文件内容操作。而且教程摒弃了完美操作,列举出常见错误和解决方式管道、标准输入输出解惑Linux下多种信息输出方式。

文件排序和FASTA文件操作中简述了awksed的使用,作为一个引子。本篇则详细列举关于awk常用的操作和一些偏门的操作。

awk基本参数解释

awk擅长于对文件按行操作,每次读取一行,然后进行相应的操作。

awk读取单个文件时的基本语法格式是awk 'BEGIN{OFS=FS="\t"}{print $0, $1;}' filename

读取多个文件时的语法是awk 'BEGIN{OFS=FS="\t"}ARGIND==1{print $0, $1;}ARGIND==2{}' file1 file2

awk后面的命令部分是用引号括起来的,可以单引号,可以双引号,但注意不能与内部命令中用到的引号相同,否则会导致最相邻的引号视为一组,引发解释错误。

OFS: 文件输出时的列分隔符 (output field separtor)

FS: 文件输入时的列分隔符 (field separtor)

BEGIN: 设置初始参数,初始化变量

END: 读完文件后做最终的处理

其它{}:循环读取文件的每一行

$0表示一行内容;$1, $2, … $NF表示第一列,第二列到最后一列。

NF (number of fields)文件多少列;NR (number of rows) 文件读了多少行: FNR 当前文件读了多少行,常用于多文件操作时。

a[$1]=1: 索引操作,类似于python中的字典,在ID map统计中有很多应用。

常见操作

  • 针对特定列的计算,比如wig文件的标准化
ct@ehbio:~/sxbd$ cat ehbio.wig
variableStep chrom=chr2
300701    12.5
300702    12.5
300703    12.5
300704    12.5
300705    12.5
ct@ehbio:~/sxbd$ awk 'BEGIN{OFS=FS="\t"}{$2=$2*10^6/(2.5*10^6); print $0}' wig 
tep chrom=chr2    0
300701    4.4
300702    4.8
300703    4
300704    4.8
300705    4.8
  • 计算某列内容出现的次数。
ct@ehbio:~/sxbd$ cat count 
ID    Type
Pou5f1    Pluripotency
Nanog    Pluripotency
Sox2    Neuron
Tet1    Epigenetic
Tet3    Epigenetic
Myc    Oncogene
ct@ehbio:~/sxbd$ awk 'BEGIN{OFS=FS="\t"}{if(FNR>1) a[$2]+=1;}END{print "Type\tCount"; for(i in a) print i,a[i];}' count
Type    Count
Neuron    1
Epigenetic    2
Oncogene    1
Pluripotency    2

# 这个也可以用下面方式代替,但不直接
ct@ehbio:~/sxbd$ tail -n +2 count | cut -f 2 | sort | uniq -c | sed -e 's/^  *//' -e 's/  */\t/' 
2    Epigenetic
1    Neuron
1    Oncogene
2    Pluripotency
  • 之前也提到过的列操作,从GTF文件中提取启动子区域
sed 's/"/\t/g' GRCh38.gtf | awk 'BEGIN{OFS=FS="\t"}{if($3=="gene") {ensn=$10; symbol=$16; if($7=="+") {start=$4-1; up=start-1000; if(up<0) up=0; dw=start+500; print $1,up, dw, ensn, symbol, $7;} else if($7=="-") {start=$5-1; up=start+1000; dw=start-500; if(dw<0) dw=0; print $1,dw,up,ensn,symbol,$7}}}' | sort -k1,1 -k2,2n >GRCh38.promoter.bed
  • 数据矩阵的格式化输出
ct@ehbio:~/sxbd$ cat numertic.matrix 
ID    A    B    C
a    1.002    1.234    1.999
b    2.333    4.232    0.889
ct@ehbio:~/sxbd$ 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
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; }}'
  • 筛选差异基因
ct@ehbio:~/sxbd$ 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

ct@ehbio:~/sxbd$ 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

ct@ehbio:~/sxbd$ 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
  • 筛选差异基因存储到不同的文件
ct@ehbio:~/sxbd$ 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 
ct@ehbio:~/sxbd$ 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、提取信息、合并信息等
ct@ehbio:~/sxbd$ 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
ct@ehbio:~/sxbd$ vim ensm
ct@ehbio:~/sxbd$ cat ensm 
ENSG00000281244
ENSG00000281165
ENSG00000282873
ct@ehbio:~/sxbd$ awk 'BEGIN{OFS=FS="\t"}ARGIND==1{if(FNR>1) ensm2entrez[$1]=$3;}ARGIND==2{print ensm2entrez[$1];}' id_map ensm
11093
11182
605

# 替代解决方案,注意 -w的使用,避免部分匹配。最稳妥的方式还是使用awk。

ct@ehbio:~/sxbd$ grep -w -f ensm id_map | cut -f 3
11093
11182
605
  • 转换大小写, toupper, tolower
ct@ehbio:~/sxbd$ cat symbol 
Tgm4
Dgkd
Abo

ct@ehbio:~/sxbd$ 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
  • awk数值操作
# log2对数
awk 'BEGIN{OFS="\t";FS="\t"}{print log($0)/log(2)}' file

# 取整,四舍五入
awk 'BEGIN{OFS="\t";FS="\t"}{print int($1+0.5);}' file
  • awk定义函数
awk 'function abs(x){return ((x < 0.0) ? -x : x)}BEGIN{OFS="\t";FS="\t"}{pos[1]=$1;pos[2]=$2;pos[3]=$3;pos[4]=$4; len=asort(pos);for(i=len;i>1;i--) print abs(pos[i]-pos[i-1]);}' file
  • 字符串匹配
ct@ehbio:~/sxbd$ cat ens.bed
1    100    105
2    100    105
3    100    105
Mt    100    105
X    100    105
ct@ehbio:~/sxbd$ 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
  • 字符串分割
ct@ehbio:~/sxbd$ 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
ct@ehbio:~/sxbd$ 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
  • awk脚本
cat <<END >grade.awk
f ( avg >= 90 ) grade="A";
else if ( avg >= 80) grade ="B";
else if (avg >= 70) grade ="C";
else grade="D";
END

awk -f grade.awk grade
  • awk给每行增加行号,使其变为唯一
awk 'BEGIN{OFS="\t";FS="\t"}NR!=1{$4=$4"_"NR;print $0}' file

糅合操作

  • awk中执行系统命令 (注意引号的使用)
# 系统命令组成字符串,交给system函数运行
awk 'BEGIN{OFS=FS="\t"}{system("mv "$1".fq "$2".fq");}' input_mat
  • awk 引用系统变量
ct@ehbio:~/sxbd$ echo 1 | awk -v ehbio="shengxinbaodian" -v ehbio2="SXBD" '{print ehbio, ehbio2;}'
shengxinbaodian SXBD

学会了基本命令,生信分析中还有一大块是使用已经安装好的工具,针对软件安装中遇到的问题,推出了系列文章:环境变量和可执行属性彻底释义环境变量的概念;列举出Linux下软件安装的各种方法,并针对DockerConda分别发文介绍。

本文分享自微信公众号 - 生信宝典(Bio_data),作者:陈同

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2017-09-25

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 父爱无言!科学家首次发现来自爸爸的基因,可以通过胎儿来控制妈妈对宝宝的爱

    近期,卡迪夫大学研究者发表在PLOS Biology上的一项研究发现,宝宝体内一个名为Phlda2的基因竟然能够影响妈妈是否关爱宝宝[1]!当这个来自鼠爸的印记...

    生信宝典
  • Python学习极简教程 (一)

    Python 教程 欢迎来到Python的世界,本教程将带你遨游Python,领悟Python的魅力。本教程专注于帮助初学者,尤其是生物信息分析人员快速学会P...

    生信宝典
  • 增强火山图,要不要试一下?

    包作者:Kevin Blighe撰文:协和医学院 苑晓梅编辑:生信宝典 时间:2019-06-03

    生信宝典
  • NodeJs —— 在Visual Studio中开发C++插件之环境配置

         最新版地址:  https://www.nasm.us/pub/nasm/releasebuilds/2.13.03/

    strawberry
  • 最全常见算法工程师面试题目整理(二)

    接着上回写的《最全常见算法工程师面试题目整理(一)》,继续填接下来的坑。 11boost算法的思路是什么样的?讲一下你对adaboost 和 gbdt的了解? ...

    用户1332428
  • 总结:常见算法工程师面试题目整理(二)

    答: boost的核心思想不同于bagging,它在基于样本预测结果对照与真实值得差距,进行修正,再预测再修正,逐步靠近正确值。

    sladesal
  • 小程序富文本解析的「伪需求」,从wxParse到towxml的坑

    其实有很多场景会用到富文本框「通常后台维护一长串html文本,前台进行渲染展示」。但由于小程序的一些特殊性,无法直接渲染html,因此类似wxParse的开源组...

    Bug生活2048
  • 博客网页导致电脑CPU飙升的问题解决记录

    已经有好几个访客朋友匿名反馈只要打开我的博客电脑的 CPU 就狂转: ? 因为忙一直也没当一回事,一是我自己的 MacbookAir 打开并没有异常,二是因为我...

    张戈
  • Docker搭建自己的Gitlab CI Runner

    目录: Gitlab CI介绍 环境、软件准备 安装、注册并启动Gitlab Runner 定义项目构建流程 FAQ 1、Gitlab CI介绍 CI:持续...

    哎_小羊
  • 基础篇章:关于 React Native 的props,state,style的讲解

    (友情提示:RN学习,从最基础的开始,大家不要嫌弃太基础,会的同学请自行略过,希望不要耽误已经会的同学的宝贵时间) React Native看起来很像React...

    非著名程序员

扫码关注云+社区

领取腾讯云代金券