生信数据预处理的Linux三大神器
> 我是hoptop，一个话很多，但努力保证写作简洁的生物信息初学者。目前在巩固Linux基础，所以会写一些Linux命令在生物信息的具体应用的文章。
阅读此文前，强烈建议先去了解一下正则表达式。
正则表达式(Regular Expression,RE,regex)是一个强大字符串操作引擎，可用于检索、替换那些符合某个模式(规则)的文本的。
grep:最快的文本搜索工具
如果你想要成功，就要做到你所在领域的最好。——中二的hoptop
grep就是在文本提取和匹配上最快的工具，因为它只有一个目标，在每一行找匹配的内容，并且在这个任务上没有其他程序比他是做的更好。
继续以拟南芥基因组和注释文件作为练习对象。
wget -c -4 -q http://www.arabidopsis.org/download_files/Genes/TAIR10_genome_release/TAIR10_chromosome_files/TAIR10_chr_all.fas &
wget -c -4 -q http://www.arabidopsis.org/download_files/Genes/TAIR10_genome_release/TAIR10_gff3/TAIR10_GFF3_genes.gff &在注释文件中查找某一个基因如AT5G25475
$ grep "AT5G25475" TAIR10_GFF3_genes.gff | head -5
Chr5 TAIR10 gene 8867797 8869821 . - . ID=AT5G25475;Note=protein_coding_gene;Name=AT5G25475
Chr5 TAIR10 mRNA 8867797 8869652 . - . ID=AT5G25475.1;Parent=AT5G25475;Name=AT5G25475.1;Index=1
Chr5 TAIR10 protein 8867920 8869492 . - . ID=AT5G25475.1-Protein;Name=AT5G25475.1;Derives_from=AT5
G25475.1
Chr5 TAIR10 five_prime_UTR 8869493 8869652 . - . Parent=AT5G25475.1
Chr5 TAIR10 CDS 8869391 8869492 . - 0 Parent=AT5G25475.1,AT5G25475.1-Protein;希望搜索AT5G25475但是不包含CDS(-v:反向选择)
$ grep "AT5G25475" TAIR10_GFF3_genes.gff | grep -v 'CDS' | head -5
Chr5 TAIR10 gene 8867797 8869821 . - . ID=AT5G25475;Note=protein_coding_gene;Name=AT5G25475
Chr5 TAIR10 mRNA 8867797 8869652 . - . ID=AT5G25475.1;Parent=AT5G25475;Name=AT5G25475.1;Index=1
Chr5 TAIR10 protein 8867920 8869492 . - . ID=AT5G25475.1-Protein;Name=AT5G25475.1;Derives_from=AT5
G25475.1
Chr5 TAIR10 five_prime_UTR 8869493 8869652 . - . Parent=AT5G25475.1
Chr5 TAIR10 exon 8869391 8869652 . - . Parent=AT5G25475.1反向选择与head,wc连用的时候，可以先排除一些注释信息(‘#’),再操作。
基因组中查找某一段特定序列并查看上下文(-A n：显示后n行-B，n：显示前n行，-C= -A,-B)
$ grep -A 2 'TTATTGTTGTTAAGAAAAAAGG' TAIR10_chr_all.fa
GTCGCCGCATTTTGTAATGCATACTTGTCTCTGTTATTGTTGTTAAGAAAAAAGGAGCACAAGTTGAGCAATGAAATAA
AATTGAATGGGCTAATGCTACAATCCCTTTTAATCAGCACAAATTGAATTAAGTTGAGGTGATTAAAAGGGATCTATCT
AGGTTTGTGGCAACAATAATAAAATGGAATCACAAACAAACTCCATAAAGGTAACCCTAAAAAAGGAGGGAAATCGCAA在基因组中统计某一段序列出现的个数 （-c : count)
$ grep 'TTATTGTTGTTAAGA' TAIR10_chr_all.fa -c
2
# 其实还可以接wc -l
$ grep 'TTATTGTTGTTAAGA' TAIR10_chr_all.fa | wc -l
2只返回匹配到的内容
$ grep 'TTATTGTTGTTAAGA' TAIR10_chr_all.fa -o
TTATTGTTGTTAAGA
TTATTGTTGTTAAGA如果你想查找一个以AT5G254开头以1结尾的基因，你就要用到强大的正则表达式。
$ grep 'AT5G254.*5$' TAIR10_GFF3_genes.gff
Chr5 TAIR10 gene 8834206 8837248 . + . ID=AT5G25415;Note=protein_coding_gene;Name=AT5G25415
Chr5 TAIR10 gene 8848549 8848986 . + . ID=AT5G25425;Note=protein_coding_gene;Name=AT5G25425
Chr5 TAIR10 gene 8867797 8869821 . - . ID=AT5G25475;Note=protein_coding_gene;Name=AT5G25475awk: 强大的文本操作工具
在上篇文章中，我使用了许多linux基本命令去完成一些琐碎的生信数据处理任务。而在awk眼中，这些任务它一个人就能搞定。
awk的工作原理
awk擅长处理表格形式的数据。它逐行从文本中读取数据，将整行数据（record)定义为$0,然后根据指定的分隔符，将各列数据（record)分别定义为$1,$2,$3
然后使用如下结构处理数据
pattern1 {action1};pattern2 {action2};....形如if函数，当满足pattern时执行后接的action.
注意：
如果没有定义pattern,则直接执行action；
如果没有提供action,则直接输出满足pattern的内容
基本用法
首先让我们用awk实现cat,cut这些命令。
cat
$ awk '{print $0}' TAIR10_GFF3_genes.gff | head -2
Chr1 TAIR10 chromosome 1 30427671 . . . ID=Chr1;Name=Chr1
Chr1 TAIR10 gene 3631 5899 . + . ID=AT1G01010;Note=protein_coding_gene;Name=AT1G01010cut
$ awk '{print $1,$4,$5}' TAIR10_GFF3_genes.gff | head -2
Chr1 1 30427671
Chr1 3631 5899awk除了能输出指定列，还能重排，甚至换一种分隔符，比cut更加强大
$ awk '{print $4","$5","$1}' TAIR10_GFF3_genes.gff | head -1
1,30427671,Chr11进阶功能
除了这些简单的功能以外，awk还支持
算术运算(+,-,*,/,%,)，
逻辑运算(==,!=,<,>,>=,<=),
或与非(&&,||,!),
还可以进行模式匹配(a~b,a!~b)。
当然知道有这些内容还不够，我们还要在实际中使用。
实例：找到长度大于10kb且在一号染色体的注释内容
$ awk '$5 - $4 > 10000 && $1 ~ /Chr1/' TAIR10_GFF3_genes.gff | head -5
Chr1 TAIR10 chromosome 1 30427671 . . . ID=Chr1;Name=Chr1
Chr1 TAIR10 gene 373335 386847 . + . ID=AT1G02080;Note=protein_coding_gene;Name=AT1G02080
Chr1 TAIR10 mRNA 373335 386847 . + . ID=AT1G02080.1;Parent=AT1G02080;Name=AT1G02080.1;Index=1
Chr1 TAIR10 protein 373335 386682 . + . ID=AT1G02080.1-Protein;Name=AT1G02080.1;Derives_from=AT1
G02080.1
Chr1 TAIR10 mRNA 373501 386846 . + . ID=AT1G02080.2;Parent=AT1G02080;Name=AT1G02080.2;Index=1练习题: 请找到长度小于1kb,在2号或者三号染色体注释
awk还有两个特殊模式BEGIN,END，顾名思义就是在操作开始或/和结束后才执行的操作。
实例： 计算长度1号染色体cds的平均长度.
$ awk 'BEGIN {s = 0;line = 0 } ;$5 - $4 > 10000 && $1 ~ /Chr1/ { s += ( $5 - $4 );line += 1}; END {print "mean=" s/line
}' TAIR10_GFF3_genes.gff | head -5
mean=227598练习题：请计算2号染色体和3号染色体的CDS的平均长度。
awk内部有许多特殊变量,如NR，表示当前所在的行数
实例： 显示第3-5行数据类似于head -n5 | tail -n2
$ awk ' NR>=3 && NR <=5 {print $0}' TAIR10_GFF3_genes.gff
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=AT1
G01010.1
Chr1 TAIR10 exon 3631 3913 . + . Parent=AT1G01010.1awk还有一个类似于Python的字典结构，叫做关联数组（associative array).使用这个我们就可以实现cat xx.file cut -f '2'| sort | uniq -c | sort类似的操作
$ awk '$1 ~ /Chr1/ { feature[$3] += 1};\
END { for (k in feature) print k "\t" feature[k]}' TAIR10_GFF3_genes.gff
pseudogene 241
ncRNA 144
exon 58048
chromosome 1
snRNA 2
protein 9263
...除此之外，他还可以使用if,for,while，还有许多常用的函数（见下表），所以说awk真的可以当做一门编程语言学习了。不过如果真的需要特别多的处理步骤，其实你可以直接使用Pyhton/perl等更强大的编程语言了。
| 函数 | 功能 |
|---|---|
| length(s) | s的长度 |
| tolower(s) | 转成小写 |
| toupper(s) | 转成大写 |
| substr(s,i,j) | 返回s的i-j部分 |
| split(s,x,d) | 根据分隔符d分隔数据，赋值给x |
| sub(f,r,s) | 根据正则f从r中提取数据到s |
番外:bioawk 生物版awk
虽然awk非常强大，但是如果它能够自动识别我们生物数据格式（如GFF/GTF,FASTQ/FASTA/,BED），自动分配record该多好。
请大家膜拜一下Heng Li，然后去下载
具体用法就不细说了。可以将bioawk代替awk重复前面的练习和实例。
sed: 流处理工具
sed的功能同样也很强大，但是这里只负责介绍两个常用的模块：文本替换和显示特定行。
假如有一个文件a.txt如下内容：
chrom1 123 456
chrom2 123 456
chrmo3 123 456我们可以很方便的使用sed将chrom替换成chr:
# sed 's/pattern/replacement/' file
sed 's/chrom/chr/' a.txt
chr1 123 456
chr2 123 456
chr3 123 456假设我们只想打印第2行，处理可以用head加tail或者awk完成，还可以用sed
#sed -n 'n1,n2p' file.txt
# -n 只输出编辑部分;p打印;n1,n2从以n1到n2
sed -n '2p' a.txt由于许多sed的功能都和前面的awk和grep冲突，所以这里只介绍了这两个功能。
如果想了解更多 awk和sed的功能，可以到进入深入的学习。
下周预告：
如何避免无谓的重复工作消耗自己对生信的热情，如何让自己的错误有迹可循，让我们下周一不见不散，学习使用脚本让生活更加美好。
客官留步：
写文章最担心就是“自嗨”，为了写出更有用的东西，让你们觉得看完这篇文章没有浪费宝贵的事情，所以麻烦大家稍稍点下的投票呗。
推荐阅读：生信分析中基本Linux命令的使用