「轮子」实现人类基因id和蛋白序列的对应

我也不知道为啥又在造轮子，就是一个需求，我用肉眼可见知道用ticyverse一定可以实现，我就做了。大概搜索还是我的短板吧。 有时候想想或许哪里有现成的只是我没找到，包括函数也是，我基于已有的知识来写的代码，经常发现有现成的函数只是我不知道。

需求

从ensembl下载的pep序列文件形式如下，奇怪的是他的name是tianscript，而不是基因id。所以需要吧geneid从Annot里面提出来放到第一列。或许这个工作在linux里面也可以实现，可我沉迷于R无法自拔啦啦啦。

思路

（1）提取到整个Annot
（2）使用空格作为分隔符，然后留下gene（也就是分隔完成后的第四个元素）。
（3）把gene和冒号去掉，顺便把id后面的小数点一二也去掉
（4）数据框专项fasta

代码

#读取人类pep的fasta格式
library(seqinr)
fastafile <- read.fasta(file = "Homo_sapiens.GRCh38.pep.all.fa", 
                        as.string = TRUE,
                        forceDNAtolower = FALSE)
#修改fasta的第一行为基因名（讨厌linux，喜欢用R）
gene <- vector("character")
sequence <- vector("character")
for (i in 1:length(fastafile)){
  gene[[i]] =attr(fastafile[[i]],'Annot')
  sequence[i] =fastafile[[i]]
}
df <- data.frame(gene = gene,sequence = sequence) 
(test2 <- apply(df,
                1,
                function(x){
                        str_split(x[1],' ',simplify=T)
                }))
#现在生成的第一列形如“gene:ENSG00000211923.1”，我们不需要小数点，也不要前面的“gene”，都丢掉
id <- vector("character")
for(i in 1:length(test2)){
  id[[i]] <- test2[[i]][4]
}
df$gene <- id
df2 <- df%>%
  separate(gene,into = c("drop","gene","drop2")) %>%
  select(-c(1,3))
#看看是否有重复
dump <- function(x){
  count <- vector("integer")
  for(i in 1:ncol(x)){
    count[i]=nrow(x[!duplicated(x[,i]),])
  }
  print(c(nrow(x),count))
}
dump(df2)
df2=df2[!duplicated(df2[,1]),]
all_recs=paste(
  apply(
    df2,
    1,
    function(x){
      paste0('>',x[1],'\n',x[2])
    }
  ),collapse = '\n') #生成fasta格式的字符串
save(all_recs,file = 'df_gene.Rdata')