美文网首页编程学习
Linux学习 - awk使用

Linux学习 - awk使用

作者: 生信宝典 | 来源:发表于2018-04-26 20:18 被阅读37次

    欢迎关注天下博客：http://blog.genesino.com/2017/09/linux-awk/
    Jump to...

    1. awk基本参数解释
    2. 常见操作
    3. 糅合操作
    4. 生信宝典

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

    在文件排序和FASTA文件操作中简述了awk和sed的使用，作为一个引子。本篇则详细列举关于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下软件安装的各种方法，并针对Docker和Conda分别发文介绍。

    生信宝典

    <footer class="entry-meta" style="box-sizing: border-box; display: block; font-size: 0.75rem; text-transform: uppercase; color: rgba(187, 187, 187, 0.8); margin: 50px 30px 30px; text-align: center; font-family: Lato, Calibri, Arial, sans-serif; font-style: normal; font-variant-ligatures: normal; font-variant-caps: normal; font-weight: 400; letter-spacing: normal; orphans: 2; text-indent: 0px; white-space: normal; widows: 2; word-spacing: 0px; -webkit-text-stroke-width: 0px; text-decoration-style: initial; text-decoration-color: initial;">LINUXBIOINFOCHENTONG
    版权声明：本文为博主原创文章，转载请注明出处。

    image.png image.png

    </footer>

    相关文章

      网友评论

      本文标题：Linux学习 - awk使用

      本文链接：https://www.haomeiwen.com/subject/pfnmlftx.html

      延伸阅读

      深度阅读

      • 您也可以注册成为美文阅读网的作者，发表您的原创作品、分享您的心情！

      栏目导航

      热点阅读

      编程学习

      关于我们|服务条款|联系我们|Linux学习 - awk使用|投稿指南|网站地图|RSS订阅|排版工具|手机版

      提供经典美文摘抄,优美散文欣赏,现代诗歌精选,短篇小说,心情随笔,表白情书范文,故事会在线阅读欣赏

      Copyright © 2014-2023 Haomeiwen.com All Rights Reserved. 好美文阅读网 版权所有

      备案信息：桂公网安备 45052102000051号 · 桂ICP备13007215号-3

      本站所收录作品、热点评论等信息部分来源互联网，目的只是为了系统归纳学习和传递资讯

      所有作品版权归原创作者所有，与本站立场无关，如不慎侵犯了你的权益，请联系我们告知，我们将做删除处理！