刘小泽写于19.3.1
最近一直忙毕业论文，感觉过去的一周像是过去的一个月一样，快接近尾声了，知道我的朋友都知道，我英语单词、配音、Peak连续打卡几百天，现在与其称之为打卡，倒不如对我是一种生活方式，是不做就睡不着觉的事情，我也很开心有这样的一个习惯。因此，虽然这段时间很忙，但是推送也在一直坚持写，anyway，等熬过这段时间就爽了~
今天在火车上用手机闲来无事，一边听着“潘吉Jenny告诉你”（嗯，是我喜欢的英语播客，也推荐给你哦），一边重新复习了下Shell单行命令
写在前面
文本文件想必是日常中完成一个大的生信计算项目，方法绝对不止一种，但最先想到的肯定是最简单的方法，那么有什么方法比以一行命令来的更优雅呢？
神奇的Awk
不止提取某列这么简单
Awk应该是功能最强大的命令行工具了，厉害的大牛一般不需要编程写脚本，awk就是他们的脚本，而且一般都是单行完成，我们没那么厉害，先从基础学起咯，一点点来
awk中一般数字就代表列号（除了0，因为0表示打印整行）
想输出第4列：awk '{print $4} a.gtf' ，这样会直接输出到屏幕，当然可以直接将结果保存到一个文件中，这个操作叫做重定向awk '{print $4} a.gtf' >a.txt ，又或者可以利用管道符号将结果进行下一步操作awk '{print $4} a.gtf' | COMMAND
Awk本身就是一门编程语言，因此也要有输入文件、命令、输出文件。需要注意的是，awk的全部需求都是要写在{}中，并且是对文件一行一行进行操作。
默认情况下，awk会把文件的每一行按照空格分成几部分（Filed），当然按什么分隔可以人为设定，然后将每一个部分赋给$1, $2,... $N这样的变量【$N表示当前行的最后一部分，另外前面也提到了，$0表示输出整行】，然后默认的操作是print
实际操作
其中有许多实用小操作哦
#下载测试数据
wget https://raw.github.com/nachocab/nachocab.github.io/master/assets/transcriptome.gtf
下载好一个文件后，首先你肯定要看看它包含了什么，于是可能你会用cat或者head查看，保证让你看过就想放弃生信~
第一步
但任何时候，一旦有放弃的念头，就要想想，世界上这么多人做相似的事情，为什么他们没有放弃？是不是有更好的解决办法？试着解决而不是逃避，这样会发现越来越有趣
less -SN transcriptome.gtf 再看一下，发现是不是好多了？并且还带上了行号？这就是-SN的作用
调整格式
【这样就可以上下翻页了，并且shift + g跳到底部，gg返回顶部，非常好用】
这个GTF文件有9列（虽然最后一列很长，但它相当重要），意思分别是
chromosome, annotation source, feature type, start, end, score, strand, and phase
其中值得一提的是最后一列phase：它以键值对的形式给出，同一键值对的键与值用空格分隔，不同键值对用分号隔开，有一些键（Key）是可选的，一些键是大部分GTF文件中都有的，值（Value）是用引号包围的【这些信息都可能为以后准确提取大量文本提供帮助】
问题
我么知道编码区有的存在许多外显子（Exon），也有的只有一个。那么现在要从GTF中提取只有一个外显子（exon）的蛋白编码基因，该怎么做？
首先我们要找到基因有哪些：
awk '$3 == "gene"' transcriptome.gtf | head | less -SN
# awk所有的规则制定都在单引号''中完成，其中$3指定第三列，后面==表示要满足什么要求，具体的匹配条件需要用双引号包围，例如“gene”
对第3列进行过滤后，我们只想要过滤后的第9列
awk '$3 == "gene" {print $9}' transcriptome.gtf | head | less -S
# 前面提到了，awk需要三大块：输入内容、运行程序、输入内容。如果不写{}的话，即不指定输出什么内容，awk默认输出全部内容
gene_id
gene_id
gene_id
gene_id
...
但是第9列不是很多键值对组成的吗？这列怎么只有一列gene_id呢？
这里又需要学习一下awk对列（column）的理解【有没有发现，生信中许多知识都是交叉学习的，学这个知识点往往能带出来其他知识，逻辑能力就是这样锻炼的吧】
我们知道，awk默认是按照空白进行分隔。空白包括两种：空格和tab，那么awk是按照什么分隔的呢？看个例子就知道：
echo "1 2 3" | awk '{ print $2 }'
# 2
echo "1\t2\t3" | awk '{ print $2 }'
# 结果为空
虽然awk默认用空格（space）分隔，但是我们也可以进行指定：
echo "1,2,3" | awk -F "," '{ print $2 }'
# 设置按，分隔，结果输出2
echo "1 2 3" | awk -F "\t" '{ print $2 }'
# 设置按\t分隔，结果为空
echo -e "1\t2\t3"" | awk -F "\t" '{ print $2 }'
# 设置按\t分隔，结果为2
看到这里是-F起的作用，它可以选择Field。看到第二行命令没有输出结果，因为它找不到任何的\t进行分隔，因此只能将全部内容输出，但全部内容却作为了一整列，我们想输出第二列$2，因此没有内容，但如果命令将$2改为$1，结果就是全部的1 2 3
好了，知道了-F的重要性以后，我们就知道为什么上面的awk '$3 == "gene" {print $9}' transcriptome.gtf命令只输出了gene_id：没有指定分隔符，对整个文件来讲这9列是tab分隔，而在第9列内部，键与值之间又是空格分隔。因此如果不指定分隔符的话，默认就是找空格所在位置，也就是第九列的起始gene_id处；指定了分隔符以后，就可以按照tab找到真正的第九列
awk -F "\t" '$3 == "gene" { print $9 }' transcriptome.gtf | head | less -S
gene_id "ENSG00000225828.1"; transcript_id "ENSG00000225828.1"; gene_type "protein_coding"; gene_status "KNOWN"; gene_name "FAM229A" ...
欢迎关注我们的公众号～_～
我们是两个农转生信的小硕，打造生信星球，想让它成为一个不拽术语、通俗易懂的生信知识平台。需要帮助或提出意见请后台留言或发送邮件到Bioplanet520@outlook.com
Welcome to our bioinfoplanet!
网友评论