基因家族收缩扩张分析1/6：利用gff注释文件获得最长转录本id(R代码

本篇推文是由 山 投稿，小明负责排版整理

写在前面：“白烟囱假说”认为一切的一切都从深海热液开始。从最初的电闪雷鸣到如今纷繁的世界，为了适应不断变化的外界环境，生物体需要持续且缓慢地改变自己。46亿年的演化，基因的消亡、出现不断发生，给生物体带来了不一样的机遇······

1. Summary

基因家族扩张收缩分析可细分为六部分：

• 1、获取每个基因对应的最长编码区转录本
• 2、OrthoFinder聚类基因家族
• 3、系统发育树推断
• 4、物种分歧时间推断
• 5、基因家族扩张收缩分析
• 6、功能富集

2. 正文

本篇推文提供一段R语言代码用来根据gff格式注释文件提取基因的最长转录本

• 1、从注释文件拿到最长转录本信息
• 2、从基因组中提取序列

2.1 提取最长转录本ID

• 最后的信息输出在注释文件所在目录的longest_protein_list.txt中
• 如果是NCBI里的数据，identifier_type选CDSNAME更便于后续操作

     # 1. packages and external scripts ---------------------------------------- TODO:
    suppressWarnings(suppressMessages(library(GenomicFeatures)))

    # 2. functions ------------------------------------------------------------ TODO:


    # 3. input ---------------------------------------------------------------- TODO:
    Args <- commandArgs()
    # anno_filename：注释文件绝对路径，这里用的是gff文件
    anno_filename <- Args[6] 
    # identifier_type：最后要提取的id水平，比如GENEID，TXNAME，CDSNAME
    identifier_type <- Args[7] 

    # 4. variable setting of test module--------------------------------------- TODO:

    # 5. process -------------------------------------------------------------- TODO:
    anno <- makeTxDbFromGFF(anno_filename, format = "auto")
    tx_lens <- transcriptLengths(anno, with.cds_len = TRUE)
    tx_lens <- tx_lens[tx_lens$cds_len != 0, ] # cds长度为0的情况：假基因、编码RNA的基因
    data_split_by_gene <- split(tx_lens, f = tx_lens$gene_id, drop = TRUE)
    longest_isoform <- lapply(data_split_by_gene, function(x) {
    return(x[which.max(x$cds_len), c("tx_id", "tx_name", "gene_id")])
    })
    longest_isoform <- do.call(rbind, longest_isoform)
    mapdata <- select(anno, keys = as.character(longest_isoform$tx_id), columns = identifier_type, keytype = "TXID")
    protein_name <- mapdata[, 2]
    protein_name <- protein_name[!is.na(protein_name)]
    write.table(protein_name, paste0(dirname(anno_filename), "longest_protein_list.txt"), quote = F, row.names = F, col.names = F)

2.2 提取序列

如果是NCBI的数据，可直接从FTP地址里找到所有蛋白质的序列文件，一般标注为*_protein.faa，我们这里记为all_proteins_seq变量。如果是自己的基因功能注释文件，相应替换一下就好。

seqkit grep -f longest_protein_list.txt $all_proteins_seq >longest_protein.fa

Bingo，这部分就到这里啦～

小明补充

利用gff注释文件获得最长转录本 实现这个功能的方法有很多，本篇推文提供了一个备选项是使用R语言的代码。

需要安装GenomicFeatures这个R包，我在服务器上用conda安装

conda install bioconductor-GenomicFeatures

因为我是在linux系统下运行，上面代码中的最后一行要改成/

我从NCBI下载拟南芥的基因组注释gff文件

对应的蛋白文件

第一步是利用脚本获取最长转录本对应的蛋白id

Rscript getLongestTranscriptID.R GCA_000001735.2_TAIR10.1_genomic.gff CDSNAME

最后一个参数这个设置了3个备选，分别是 GENEID，TXNAME，CDSNAME，这里大部分时候要选择CDSNAME，不确定用哪个就把三个分别试一下，然后看那个输出和自己的蛋白id相匹配

Rscript getLongestTranscriptID.R GCA_000001735.2_TAIR10.1_genomic.gff GENEID
Rscript getLongestTranscriptID.R GCA_000001735.2_TAIR10.1_genomic.gff TXNAME

最后是用seqkit这个工具利用id列表从总的蛋白文件中提取自己想要的蛋白

seqkit grep -f longest_protein_list.txt GCA_000001735.2_TAIR10.1_protein.faa -o longest_protein.fa

欢迎大家关注我的公众号

小明的数据分析笔记本

小明的数据分析笔记本

分享R语言和python在生物信息领域做数据分析和数据可视化的简单小例子；偶尔会分享一些组学数据处理相关的内容

703篇原创内容

公众号

小明的数据分析笔记本 公众号 主要分享：1、R语言和python做数据分析和数据可视化的简单小例子；2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记；3、生物信息学入门学习资料及自己的学习笔记！