生信编程11.根据gtf格式的基因注释文件提取人所有基因的染色体坐标
有一些五六年前的学生们都成长为了各个生物信息学相关公司的小领导，而且他们都有了自己的公众号，知乎号，也算是一番人物。最近他们跟我反馈面试找不到或者说很难直接考核筛选到认真干活的生信工程师，挺有意思的。让我想起来了早在生信技能树论坛创立之初我为了引流，而规划的200个生信工程师面试题。值得继续分享：
首先我们要知道gtf文件中所包含的信息
GTF文件介绍
GTF全程为gene transfer format，主要用对基因进行注释。GTF文件共由9列组成，这两种数据的前8列都是相同的，存在一些细微的差别。
原文链接：https://blog.csdn.net/sinat_38163598/article/details/72851239
| 列 | 列名 | 包含的信息 |
|---|---|---|
| 1 | seq_id | 序列的编号，一般为染色体号或者scaffold的号 |
| 2 | source | 注释的来源，一般为数据库或者注释的机构，如果是未知来源，用'.'代替 |
| 3 | type | 注释信息的类型，一般为Gene，cDNA，mRNA，CDS等 |
| 4 | start | 该基因或者转录本在参考序列上的起始位置 |
| 5 | end | 该基因或者转录本在参考序列上的终止位置 |
| 6 | score | 对注释信息的准确性进行打分，可以是序列比对时的相似性的份E-value，也可以是基因预测是的p-value |
| 7 | strand | 仅该基因或者转录本位于参考序列的正链或者负链上 |
| 8 | phase | 该注释仅对CDs有效，仅对注释类型为“CDS”有效，表示起始编码的位置，有效值为0、1、2(对于编码蛋白质的CDS来说，本列指定下一个密码子开始的位置。每3个核苷酸翻译一个氨基酸，从0开始，CDS的起始位置，除以3，余数就是这个值，，表示到达下一个密码子需要跳过的碱基个数。该编码区第一个密码子的位置，取值0,1,2。0表示该编码框的第一个密码子第一个碱基位于其5'末端；1表示该编码框的第一个密码子的第一个碱基位于该编码区外；2表示该编码框的第一个密码子的第一、二个碱基位于该编码区外；如果Feature为CDS时，必须指明具体值。 |
| 9 | attributes | 不同属性的列表 |
这里以GRCh38为例
##description: evidence-based annotation of the human genome (GRCh38), version 32 (Ensembl 98)
##provider: GENCODE
##contact: gencode-help@ebi.ac.uk
##format: gtf
##date: 2019-09-05
chr1 HAVANA gene 11869 14409 . + . gene_id "ENSG00000223972.5"; gene_type "transcribed_unprocessed_pseudogene"; gene_name "DDX11L1"; level 2; hgnc_id "HGNC:37102"; havana_gene "OTTHUMG00000000961.2";
chr1 HAVANA transcript 11869 14409 . + . gene_id "ENSG00000223972.5"; transcript_id "ENST00000456328.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_name "DDX11L1"; transcript_type "lncRNA"; transcript_name "DDX11L1-202"; level 2; transcript_support_level "1"; hgnc_id "HGNC:37102"; tag "basic"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000362751.1";
我们的要求是提取基因的染色体坐标
首先我们要明确这个注释文件到底有几个来源
awk '{if(NR>5) print $2}' gencode.v32.annotation.gtf | sort | uniq
#245472 ENSEMBL #2654302 HAVANA
该文件中的注释信息有两个来源，我们挑选数据多的HAVANA
挑选gene信息
awk '{if((NR>5) && ($2=="HAVANA") && ($3=="gene")) print $0}' gencode.v32.annotation.gtf | perl -alne '{next unless $F[2] eq "gene" ;/gene_name \"(.*?)\";/; print "$F[0]\t$F[3]\t$F[4]\t$1" }' | head
结果文件如下：
chr1 11869 14409 DDX11L1
chr1 14404 29570 WASH7P
chr1 29554 31109 MIR1302-2HG
chr1 34554 36081 FAM138A
chr1 52473 53312 OR4G4P
chr1 57598 64116 OR4G11P
chr1 65419 71585 OR4F5
chr1 89295 133723 AL627309.1
chr1 89551 91105 AL627309.3
chr1 131025 134836 CICP27