生信数据预处理的Linux三大神器

原创 hoptop 生信媛

> 我是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=AT5G25475

awk: 强大的文本操作工具

在上篇文章中，我使用了许多linux基本命令去完成一些琐碎的生信数据处理任务。而在awk眼中，这些任务它一个人就能搞定。

awk的工作原理

awk擅长处理表格形式的数据。它逐行从文本中读取数据，将整行数据（record)定义为$0,然后根据指定的分隔符，将各列数据（record)分别定义为$1,$2,$3

然后使用如下结构处理数据

pattern1 {action1};pattern2 {action2};....

形如if函数，当满足pattern时执行后接的action.
注意：

如果没有定义pattern,则直接执行action；
如果没有提供action,则直接输出满足pattern的内容

基本用法

首先让我们用awk实现cat,cut这些命令。

$ 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=AT1G01010

$ awk '{print $1,$4,$5}' TAIR10_GFF3_genes.gff | head -2
Chr1 1 30427671
Chr1 3631 5899

awk除了能输出指定列，还能重排，甚至换一种分隔符，比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.1

awk还有一个类似于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
具体用法就不细说了。可以将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的功能，可以到http://dongweiming.github.io/sed_and_awk进入深入的学习。

下周预告：

如何避免无谓的重复工作消耗自己对生信的热情，如何让自己的错误有迹可循，让我们下周一不见不散，学习使用脚本让生活更加美好。

客官留步：

写文章最担心就是“自嗨”，为了写出更有用的东西，让你们觉得看完这篇文章没有浪费宝贵的事情，所以麻烦大家稍稍点下的投票呗。

推荐阅读：生信分析中基本Linux命令的使用

继续滑动看下一个