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

1. Summary

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

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

2. 正文

目前，似乎大家都是自己写程序提取基因的最长转录本，写起来还是有些费劲的，这里是我的方案

• 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.1 提取序列

如果是NCBI的数据，可直接从FTP地址里找到所有蛋白质的序列文件，一般标注为*_protein.faa，我们这里记为all_proteins_seq变量

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

Bingo，这部分就到这里啦～