探索猪基因组的秘密——从GTF文件开始_技术教程

作为一个非人类、小鼠等经典模式生物的生物信息工作者，常常感叹猪的基因组注释信息不够详细，不够规范。

尤其是在进行ID转换的时候，经常有一大批的基因或者蛋白找不到对应的基因名、GeneID或者其他重要的信息。

这就导致我们在进行功能分析的时候，那些研究的越充分的基因，得到的信息越多。反而是那些可能具有重要功能的新基因，得到的信息聊聊无几。

为了寻找更完整的信息，我决定从GTF文件中提取相关的ID信息，自己创建一个猪的基因ID注释表，并定期进行更新。

今天先从Ensembl数据库中猪参考基因组的GTF文件开始。该文件的下载地址是：

或者按照下图顺序登录Ensembl数据库依次点击。

下载后保存在本地，解压。

参考基因组GTF文件是该物种参考基因组的注释文件，包含该物种的已知的基因组相关的信息。它在格式上是用制表符分隔为9列的TSV（Tab-Separated Values，制表符分隔的值）文件。

其中：

第一列是seqid, 通常是染色体编号；
第二列是source，表示信息的来源；
第三列是feature,表示类型，如gene，transcript，exon，CDS等；
第四和第五列分别是区间的起始位置与终止位置；
第六列是score, 软件提供的统计值；
第七列是strand, 代表正负链的信息, +表示正链，-表示负链，?表示不清楚正负链的信息；
第八列是phase，当第三列是CDS时，需要指定翻译时开始的位置，取值范围有0,1,2；
第九列是attributes, 表示属性，必须有gene_id和transcript_id这两个属性，多个属性用分号分隔

这个文件可以用Rstudio导入并提取其中的信息。

library(stringr)
library(readr)
#导入gtf文件
gtf_in <- read_delim("./data/ensembl/Sus_scrofa.Sscrofa11.1.113.chr.gtf",
                     "\t", 
                     escape_double = FALSE, 
                     col_names = FALSE, 
                     comment = "#", 
                     trim_ws = TRUE) 
#导入时间会比较长
#看一下
View(gtf_in)

这个GTF文件共有1339823行，我想提取全部基因的ID信息，它主要储存在第9列，也就是X9这一列。

但是首先要过滤掉其他feature，只保留gene类型。

gtf_in <- gtf_ingtf_inx3='="gene",]'> dim(gtf_in)
[1] 34307     9
#说明包含了34307个基因信息

继续提取第9列中的信息

gene_info <- data.frame(str_split_fixed(unique(gtf_in$X9),";",5))
# 将第9列根据分号分隔成多个列，str_split_fix(string, pattern,n)输入向量，返回矩阵
View(gene_info)#看看得到了什么

得到了一个新的数据框，第一列是ensembl ID，第三列是gene name，第5列是gene biotype，继续提取这三列的信息

gene_info2 <- unique(data.frame("gene_id"=str_split_fixed(gene_info$X1,'"',3)[,2],
                         "gene_name"=str_split_fixed(gene_info$X3,'"',3)[,2],
                         "gene_biotype"=str_split_fixed(gene_info$X5,'"',3)[,2]))
#再用str_split_fixed获取引号内的内容，这样就获得了gene_id、gene_name、gene_biotype构成的基因注释表格
View(gene_info)#看看得到了什么

这个表就是我想要的样子了。检查一下ensembl ID和gene name有没有重复

length(unique(gene_info2$gene_id)) 
[1] 34307
# 34307行，与geneinfo2的行数一致，说明没有重复的ensembl ID
length(unique(gene_info2$gene_name)) #16171行，其余无symbol
[1] 16171
# 少于geneinfo2的行数

值得注意的事，gene_info2中唯一的基因名只有16171个，说明有将近2万个基因还没有基因名，这个数量比已知的基因数量还多。猪的基因组研究果然还是任重道远啊。

再看看Ensembl数据库中都收集了哪些类型的基因

length(unique(gene_info2$gene_biotype)) 
[1] "protein_coding"         ""                      
[3] "snRNA"                  "miRNA"                 
[5] "misc_RNA"               "snoRNA"                
[7] "ribozyme"               "pseudogene"            
[9] "Y_RNA"                  "vault_RNA"             
[11] "rRNA"                   "TR_V_gene"             
[13] "unprocessed_pseudogene"

合计有12种类型，果然Ensembl数据库整理的数据很全面嘛。

保存一下数据。

完成。下次在看看NCBI数据库的猪的GTF文件。

更多优质内容，关注微信公众号：育种与生信

探索猪基因组的秘密——从GTF文件开始

相关文章