R编程反面教材：1646992行的gff文件处理了6小时

最新推荐文章于 2024-11-12 23:43:13 发布

学生信的大叔

最新推荐文章于 2024-11-12 23:43:13 发布

阅读量139

点赞数

文章标签： r语言开发语言

本文链接：https://blog.csdn.net/qq_33408922/article/details/134110302

版权

耐心看下背景

前段时间帮客户做转录组上游分析定量，要求定量到基因水平，遇到一个非模式真核生物。它的gff注释文件示例如下：

Chr4	.	gene	215480503	215499887	.	-	.	ID=evm.TU.ctg388.13;Name=EVM%20prediction%20ctg388.13
Chr4	.	mRNA	215480503	215499887	.	-	.	ID=evm.model.ctg388.13;Parent=evm.TU.ctg388.13;Name=EVM%20prediction%20ctg388.13
Chr4	.	five_prime_UTR	215499673	215499887	.	-	.	ID=evm.model.ctg388.13.utr5p1;Parent=evm.model.ctg388.13
Chr4	.	exon	215499406	215499887	.	-	.	ID=evm.model.ctg388.13.exon1;Parent=evm.model.ctg388.13
Chr4	.	CDS	215499406	215499672	.	-	0	ID=cds.evm.model.ctg388.13;Parent=evm.model.ctg388.13
Chr4	.	exon	215480503	215481095	.	-	.	ID=evm.model.ctg388.13.exon2;Parent=evm.model.ctg388.13
Chr4	.	CDS	215480751	215481095	.	-	0	ID=cds.evm.model.ctg388.13;Parent=evm.model.ctg388.13
Chr4	.	three_prime_UTR	215480503	215480750	.	-	.	ID=evm.model.ctg388.13.utr3p1;Parent=evm.model.ctg388.13

从上面这个基因的注释来看，feature为exon 的行中attributes注释部分，并没有给出基因IDevm.TU.ctg388.13，而是只给到了转录本IDevm.model.ctg388.13 。由于整个注释文件中存在特殊情况，并不能直接将转录本ID中的model部分替换为TU。因此，我们需要做一些工作使我们能定量到基因水平。

目前的大方向有两个：

提取基因ID与转录本ID对应关系。一般非模式真核生物的注释，很少有基因能注释出可变剪切。也就是一个基因只注释出一个转录本，我们可以从gtf文件中提取基因ID和转录本ID对应关系。提取的代码可以参考我之前的代码

探索GRCh38的转录组本、基因和SYMBOL等的对应关系

给到对应的five_prime_UTR、exon、CDS和three_prime_UTR注释部分添加基因ID。

我当时用的就是第二种方法。一开始是用R写的代码，写代码10分钟，运行6小时。后来嫌弃太慢，花40分钟写了个perl版本代码，然后4秒钟解决了。。。

今天带大家看下我写的反面教材R代码（改进思路见小结部分，有兴趣的可以找一个模式生物验证下自己的思路和代码是否正确。）。

R反面教材代码

library(vroom)
library(stringr)

gff <- vroom_lines(file = "./test.gff.gz")

new_gff <- vector()
for (i in 1:length(gff)) {
  message(i)
  gff_line <- gff[i]
  gff_tmp <- str_split(gff_line, pattern = "\t")[[1]]
  feature <- c("mRNA","exon","CDS")
 if(length(gff_tmp) == 9){
  if(gff_tmp[3] == "gene"){gene <- gff_tmp[9] %>% str_extract(pattern = "ID=.*?;") %>% str_remove("ID=") %>% str_remove(";")}
  if(gff_tmp[3] %in% feature) gff_line <- paste0(gff_line,";gene_id=",gene)
 }
  new_gff <- c(new_gff, gff_line)
}

vroom_write_lines(new_gff, file = "new_test.gff.gz", eol = "\n")
save.image(file = "new_gff.Rdata")

原注释文件有1646992行，这里是尝试用R按行处理的数据，还使用了for循环。我是知道for循环慢，但是也没有想到for循环这么慢啊！单纯的处理gff就花费6小时。