写在前面：根据“白烟囱假说”，一切的一切都从深海热液开始。从最初的单细胞生物到如今繁杂的世界，为了适应不断变化的外界环境，生物体需要不断且缓慢地改变自己。46亿年的演化，基因的消亡、出现不断发生，给生物体带来了不一样的机遇······
1. Summary
基因家族扩张收缩分析可细分为六部分：
- 获取每个基因对应的最长编码区转录本
- OrthoFinder聚类基因家族
- 系统发育树推断
- 物种分歧时间推断
- 基因家族扩张收缩分析
- 功能富集
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，这部分就到这里啦～
