参考

github: 李恒：bioawk
简书：刘小泽：一个神奇的小软件bioawk

会一点awk很容易上手

1. 介绍

https://github.com/lh3/bioawk

2. 安装

  1. conda install bioawk

3. 参数

3.1 -t

  1. bioawk -t # 相当于
  2. bioawk -F '\t' -v OFS="\t" # 即分隔符、输出分隔符设置为制表符

3.2 -c

3.2.1 help

  1. $ bioawk -c help
  2. bed:
  3.         1:chrom 2:start 3:end 4:name 5:score 6:strand 7:thickstart 8:thickend 9:rgb 10:blockcount 11:blocksizes 12:blockstarts
  4. sam:
  5.         1:qname 2:flag 3:rname 4:pos 5:mapq 6:cigar 7:rnext 8:pnext 9:tlen 10:seq 11:qual
  6. vcf:
  7.         1:chrom 2:pos 3:id 4:ref 5:alt 6:qual 7:filter 8:info
  8. gff:
  9.         1:seqname 2:source 3:feature 4:start 5:end 6:score 7:filter 8:strand 9:group 10:attribute
  10. fastx:
  11.         1:name 2:seq 3:qual 4:comment

在这些有规范格式的文件中，对各成分定义了变量，并且添加了一些子命令，非常方便！

3.2.2 fastx

对最熟悉的fasta/fastq文件进行操作

  1. $ cat a.fa
  2. >gene1
  3. actcagatcgatcca
  4. cagatcaggcagtta
  5. cag
  6. >gene2
  7. cagatcagatcagcg
  8. cagatcgcgacagat
  9. cagtctt
  1. $ bioawk -t -c fastx 'BEGIN{print "#geneid\tlength\tGC_content\tsequence"}{print $name,length($seq),gc($seq),$seq}' a.fa
  2. #geneid length  GC_content      sequence
  3. gene1   33      0.484848        actcagatcgatccacagatcaggcagttacag
  4. gene2   37      0.513514        cagatcagatcagcgcagatcgcgacagatcagtctt

gc(): 返回序列gc含量
length()：返回长度
revcomp(): 返回互补链
…

3.2.3 hdr or header

甚至可以按自己的第一行进行变量定义

  1. $ cat b.txt
  2. geneid  T1_1    T1_2    T1_3    T2_1    T2_2    T2_3
  3. gene1   1       3       9       2       8       5
  4. gene2   3       4       5       8       0       8
  5. $ bioawk -t -c header '{print $geneid,$T1_1,$T2_1}' b.txt
  6. geneid  T1_1    T2_1
  7. gene1   1       2
  8. gene2   3       8