生信单行命令的美，慢慢体会（二）

刘小泽写于19.3.2

上一次说到了怎样找到特定的列（比如：$3 == "gene"以及返回特定的列（{print $1}）。今天继续探究上次的文件

我们看到GTF文件中存在许多的键值对，上一次我们知道了怎么将全部的第九列提取出来：

awk -F '\t' '$3 == "gene" {print $9}' transcriptome.gtf | head | less -S

那么现在，如果想提取第9列中的关于gene_type 信息（也就是第三列），需要在上面代码的基础上用awk再次指定分隔符为; ，接着提取第三列就好

awk -F '\t' '$3 == "gene" {print $9}' transcriptome.gtf |awk -F ';' '{print $3}'| head | less -S
# 得到如下
 gene_type "protein_coding"
 gene_type "antisense"
 gene_type "protein_coding"
 gene_type "pseudogene"
 gene_type "sense_overlapping"
 gene_type "protein_coding"
 gene_type "pseudogene" 
 ...

因为存在许多的gene_type，而我们就是想要蛋白编码区（protein_coding）的，于是可以指定筛选条件：【需要注意一点："protein_coding"的前后引号需要用反斜线进行转义，便于awk识别筛选条件】

awk -F "\t" '$3 == "gene" { print $9 }' transcriptome.gtf | \
    awk -F "; " '$3 == "gene_type \"protein_coding\""' | \
    head | less -S

或者：如果感觉反斜线的使用不美观或者经常忘记使用，那么也可以用~表示部分匹配，只要搜索的列中包含你输入的字符，那么就会输出结果

awk -F '\t' '$3 == "gene" {print $9}' transcriptome.gtf | \
    awk -F ';' '$3 ~ "protein_coding"' |\
    head | less -S

到这里我们已经输出了所有的protein_coding基因，那么之前的问题又来了：

如何输出只包含一个exon的protein_coding基因呢？

这里发现：开始我们选择的提取条件是第三列的gene，但是第三列除了gene还有其他信息，统计一下：

awk '{print $3}' transcriptome.gtf | sort | uniq -c | sort -r
# 有用信息如下
   1203 exon
    730 CDS
    273 UTR
    228 transcript
     94 start_codon
     75 stop_codon
     63 gene

最多的还是exon，于是就可以将$3 == "gene"改成$3 == "exon" ，挑出全部的exon基因名信息（在分号隔开的第5列）

awk -F '\t' '$3 == "exon" {print $9}'  transcriptome.gtf | \
    awk -F ';' '$3 ~ "protein_coding"' |\
    head | less -S
# 得到的结果
      1  gene_name "NADK"
      2  gene_name "C1orf86"
      3  gene_name "RER1"
      4  gene_name "SLC2A5"
      5  gene_name "LZIC"
      6  gene_name "FBXO6"
      7  gene_name "NPPA"
      8  gene_name "PLOD1"
      9  gene_name "FHAD1"
     10  gene_name "CLCNKA"
     ...
 # 还有另外一种实现方式(因为发现GTF第9列中既有分号又有空格分隔，还有引号。因此可以直接将分号和引号去掉，然后直接指定分隔符是空格)，接下来就可以直接指定第6列为protein_gene时将第10列name打印出来
 awk -F '\t' '$3 == "exon" {print $9}' transcriptome.gtf |\
    tr -d ";\"" |awk -F ' ' '$6 == "protein_coding" {print $10}'|\      head |less -S
# 结果如此，更加方便下面进行 
NADK
C1orf86
RER1
SLC2A5
LZIC
FBXO6
NPPA
PLOD1
FHAD1

这样就得到了exon与protein_gene的对应关系

接下来到了计算每个基因对应的exon出现多少次的时候了，其实就是统计上面的结果

awk -F '\t' '$3 == "exon" {print $9}' transcriptome.gtf | tr -d ";\"" |awk -F ' ' '$6 == "protein_coding" {print $10}' | sort | uniq -c | sort -r >gene_exon_num.txt
# 结果
4 TTN
4 KIF12
3 LAT2
3 ITSN2
3 CBWD6
2 UBXN11
...

然后看看有多少基因是只有一个exon的：

awk '$2 == 1' gene_exon_num.txt | wc -l 

---

欢迎关注我们的公众号～_～　　
我们是两个农转生信的小硕，打造生信星球，想让它成为一个不拽术语、通俗易懂的生信知识平台。需要帮助或提出意见请后台留言或发送邮件到Bioplanet520@outlook.com

Welcome to our bioinfoplanet!