chatGPT辅助学员造轮子去解析fasta文件来搜索指定的多个基因的CDS长度

起

前一段时间，同事提了一个需求，大致内容是：给出若干个基因名，去基因数据库匹配，求出这些基因的CDS序列的长度。他说：

你写个for循环就行了。

当时月初才开始学R语言，显然没法完成。我试图用excel的vlookup解决，但是该CDS文件太长（用R读取后发现有138w+行），excel都无法完整读取，所以彼时这个问题我没有解决。现在月末学完R基础之后，这个问题可以解决了。

承

先去ensembl官网下载目的文件——大鼠的基因CDS文件，文件是fa格式，先用R读取；目的基因list是存储为一列的txt文档，一并读取。

fa <- read.csv("Rattus_norvegicus.mRatBN7.2.cds.all.fa", header = F)
mono <- read.csv("mono-716.txt", header = F)

先观察fa文件的结构，先有一行“>”起始的注释，后面有基因类型、等注释信息，有些有gene_symbol，有些则没有。起始行往下另起一行是序列，有若干行，直到下一个“>”为止。（因为是完全零基础，所以才开始学习生物信息学常见的数据文件结构，包括fasta和fastq ）

我的思路

1. 找出所有的序列起始行的行号
2. 筛选出有gene_symbol的行号
3. 确定序列所在的行，统计字符数
4. 读取目的基因字段，去上述行查找

其中第二步不是必须的。由于有些序列行不带gene_symbol，我这里直接简单粗暴删掉；如果进行注释会更好，不会遗漏。

1&2.找出目的行的行号

用for循环来执行如果fa中某行带有“>”，那就是某序列起始行，输出行号；如还带有“gene_symbol”，那就再输出一次行号并输出该行内容；循环直至最后一行。使用到的函数是str_starts()、str_detect()、str_split()和length()。

因为要执行两次if判断，如果分开写会明显增加计算量，于是，我选择了写嵌套的if，没想到引发了后面一系列问题。确实如果第一次把这一段代码写得冗余一些，后面某些问题将不存在...

遇到的第一个问题

没错，刚起头就出问题。准确说不算问题，代码能跑通但输出结果不太满意。起初代码是这样的：

fa <- read.csv("Rattus_norvegicus.mRatBN7.2.cds.all.fa", header = F)
mono <- read.csv("mono-716.txt", header = F)
i = 1
x = 1
library(stringr)
ENS_row = data.frame()
for(i in (1:nrow(fa))){
  if (str_starts(fa[i,1],">ENSRN")){  #找出所有起始行
    if (str_detect(fa[i,1],"gene_symbol")){    #找出所有带gene_symbol的起始行
      ENS_row[x,2] = i   #gene_symbol的起始行行号输出到第2列
      ENS_row[x,3:(2+length(str_split(fa[i,1]," ", simplify = T)))] = (str_split(fa[i,1]," ", simplify = T))    #输出其他信息，genen_symbol为第9列
    }
    ENS_row[x,1] = i    #所有起始行行号输出到第1列
  }  
  i = i + 1
  x = x + 1
}

查看运行结果，产生了无数空行。一个解决办法是，先生成，再删除。但是这样凭空增加了很多无效计算量。那么是哪里出现了问题？排查并尝试之后，发现是循环变量自增，也就是 i = i + 1的位置。在我这个场景下，如果放在for循环里就会产生空行。这里的执行逻辑是，每读一行就进行一次判断，放在for里将导致每次都自增一行；而我的需求是遇到序列起始行才输出。解决办法是把i = i + 1放在外层if下，于是就如期呈现了。

for(i in (1:nrow(fa))){
  if (str_starts(fa[i,1],">ENSRN")){  #找出所有起始行
    if (str_detect(fa[i,1],"gene_symbol")){    #找出所有带gene_symbol的起始行
      ENS_row[x,2] = i   #gene_symbol的起始行行号输出到第2列
      ENS_row[x,3:(2+length(str_split(fa[i,1]," ", simplify = T)))] = (str_split(fa[i,1]," ", simplify = T))    #输出其他信息，genen_symbol为第9列
    }
    ENS_row[x,1] = i    #所有起始行行号输出到第1列
    i = i + 1
    x = x + 1
  }   
}

3.确定序列所在的行，统计字符数

逻辑：如要统计某基因CDS长度，则需要获取该CDS所在行的区间，然后统计区间内的字符数。举个栗子：求上图第二行基因的CDS长度。目前已知，该基因起始行是8，下一个基因的起始行是16；则该基因的CDS行的区间是9-15，那么统计第9-15行字符数即是CDS长度。尝试增加这部分代码：

ENS_row[x,3] = ENS_row[x+1,1] - 1    #获取某基因的下一个基因行号再减1，即某基因CDS末行

遇到的第二个问题

看上去很合理，加到for循环里，跑一下：

for(i in (1:nrow(fa))){
  if (str_starts(fa[i,1],">ENSRN")){  #找出所有起始行
    ENS_row[x,1] = i    #所有起始行行号输出到第1列
      if (str_detect(fa[i,1],"gene_symbol")){    #找出所有带gene_symbol的起始行
      ENS_row[x,3] = i+1  #gene_symbol的起始行的下一行输出到第3列
      ENS_row[x,4:(2+length(str_split(fa[i,1]," ", simplify = T)))] = (str_split(fa[i,1]," ", simplify = T))    #输出其他信息，genen_symbol为第9列
    }   
    ENS_row[x,2] = ENS_row[x+1,1] - 1
    i = i + 1
    x = x + 1
  }   
}

第二列输出为空。从代码上来看，似乎能说得通？“ENS_row[x,3] = ENS_row[x+1,1] - 1”，class查看属性显示是“numeric”，可以排除字符型不能运算的原因。并且这段代码单独执行是正常的，唯独放到for循环里就会出错。仔细思考排查之后找到了原因——这段代码引入了一个未来函数。未来函数在金融上用得多，我觉得我这个符合它的定义（具体以专业人士）：

······通达信官方对未来函数的解释：某一量依赖另一量，如量A和量B，B变化使A改变，那么A是B的函数，如果是稍后的量，A是稍早的量，A跟着B变，A是B的未来函数。

未来函数，简单的说，就是带有未来信息的函数。

未来函数就像开了上帝之眼，拥有了预知的超能力。在我这个实例中，因为是逐行输出，当输出x行的时候，x+1行还没有生成，因此"ENS_row[x,3] = ENS_row[x+1,1] - 1“将返回空值。而这一切的源头，就是我使用了if的嵌套以及只写了一个for······如果多写几个for的话，那么可以在第一列和第三列生成后，再生成第二列，可以避免上述问题。那么在现有场景下如何改？

我的办法是，待下一行生成之后在生成上一行第二列的值。具体操作是引入y，使y的值落后x一个循环。即y不自增，并使y=x，整合之后代码如下：

for(i in (1:nrow(fa))){
  if (str_starts(fa[i,1],">ENSRN")){      #列出所有起始行行号
    ENS_row[x,1] = i
    if (str_detect(fa[i,1],"gene_symbol")){   #列出带有symbol的起始行行号    
      ENS_row[x,3:(2+length(str_split(fa[i,1]," ", simplify = T)))] = str_split(fa[i,1]," ", simplify = T)
      ENS_row[x,3] = i+1   #找到起始行的下一行
      ENS_row[,4] = ENS_row[,9]
       y = x      
    }   
   ENS_row[y,2] = ENS_row[y+1,1] - 1
   x = x + 1
   i = i + 1
  }

4.读取目的基因字段，去上述行查找

现在获取了CDS所在行行号，再按行号读取fa中的内容。用nchrc()统计字符数，再用sum()求和，输出到第5列。

ENS_row[y,5] = sum(nchar(fa[(ENS_row[y,3]:ENS_row[y,2]),1]))

遇到的第三个问题

这句代码在单独运行时没有问题，但是放到循环里就报错，提示参数的长度为零。查看生成的ENS_row：

从图得知，在输出到第一行第三列时出错。上述代码里”(ENS_row[y,3]:ENS_row[y,2]“在第一行为空，因此无法继续。查了下可以用is.na()处理空值，加上if判断，重开一个for循环：

a = 1
for(a in (1:nrow(mono))){
   if (!is.na(ENS_row[a,2])){
      ENS_row[a,5] = sum(nchar(fa[(ENS_row[a,3]:ENS_row[a,2]),1]))
   }
   a = a + 1
}

终于正确获得了想要的结果！

造轮子解析了fasta文件

前面我们下载的 Rattus_norvegicus.mRatBN7.2.cds.all.fa 文件是 100多万行，几万个基因的序列信息。通过我们自己的冗余的代码，终于给每个基因一个cds长度信息啦。

第四列是基因名，第五列是CDS长度，接下来就是从mono.txt中提取元素挨个比对，直到匹配成功。此时脑子里晃过好几个可行的函数，比如x%in%y啥的；不过脑瓜子一转想起一个好东西（柯南灵光一闪，头上有一条白线闪过），可以使用merge()连接！先将gene_symbol全部改为大写，以免merge的时候因大小写不一致而遗漏：

ENS_row$V4 = str_remove_all(ENS_row$V4,"gene_symbol:") 
colnames(mono) <- "name"
colnames(ENS_row) <- c("NO.","CDS_t","CDS_h","name","CDS")  #赋列名以便merge
mono[,1] = c(str_to_upper(mono$name))  
ENS_row[,4] = c(str_to_upper(ENS_row$name))
mononew = merge(mono, ENS_row, by.x = "name")

最终结果如图。这个解题到这基本结束，为了美观可以单独保存第一列和第五列即可。唯一不足的是，没有使用注释，有可能造成部分基因没有匹配到。不过目测占比不是太大，下次可以尝试加上注释。

转

扩展一下

尝试用chatGPT解题。对话是这样的：

尝试了下这玩意是跑不通的。但是它提到用grep()解题（不过查看帮助文档，R语言里grep没有 | 这个用法），确实好用，比使用for和if快多了。

ENSRNROW = data.frame(
  row = c(grep(">ENSRN",fa$V1))
  )

只要上面这段代码就能之间输出目的行号...比我自己写的不知高到哪里去了......

其实如果合理的询问chatGPT，会有更好的结果：

这段代码首先加载了seqinr库，然后使用read.fasta函数加载了fasta文件。接下来，我们创建了一个存储每个基因序列长度的向量gene_lengths。最后，我们使用循环来遍历fasta文件中的每个基因，并使用nchar函数获取每个基因的序列长度，然后将长度存储在gene_lengths向量中。最后，我们使用gene_lengths来查看每个基因的长度。

合

算是一次战斗力为五的实操（文章最后要扣题），中间连续卡了几次bug，还好都解决了。总结一下：1.新手使用if嵌套要小心，以免掉坑; 2.chatGPT可以尝试一下，会有收获，但不能全信；3.bug一定可以解决。

附我写的全部代码：

library(dplyr)
library(stringr)
fa <- read.csv("Rattus_norvegicus.mRatBN7.2.cds.all.fa", header = F)
mono <- read.csv("mono-716.txt", header = F)

i = 1  #fa
x = 1  #ens_row
y = 2 #很重要，没有y就是空值

ENS_row = data.frame()
#sym = data.frame()
#class(sym)
for(i in (1:nrow(fa))){
  if (str_starts(fa[i,1],">ENSRN")){      #列出所有起始行行号
    ENS_row[x,1] = i
    if (str_detect(fa[i,1],"gene_symbol")){   #列出带有symbol的起始行行号
      #print(i) 
      ENS_row[y,2] = ENS_row[y+1,1] - 1
      ENS_row[x,3:(2+length(str_split(fa[i,1]," ", simplify = T)))] = str_split(fa[i,1]," ", simplify = T)
      ENS_row[x,3] = i+1   #找到起始行的下一行
      ENS_row[,4] = ENS_row[,9]
      y = x
    }
   x = x + 1
   i = i + 1 
  }
}
a = 1
   for(a in (1:nrow(fa))){
   as.numeric(ENS_row[a,2])
   if (!is.na(ENS_row[a,2])){
      ENS_row[a,5] = sum(nchar(fa[(ENS_row[a,3]:ENS_row[a,2]),1]))
   }
   a = a + 1
}
ENS_row$V4 = str_remove_all(ENS_row$V4,"gene_symbol:") 
colnames(mono) <- "name"
colnames(ENS_row) <- c("NO.","CDS_t","CDS_h","name","CDS")  #赋列名以便merge
mono[,1] = c(str_to_upper(mono$name))  
ENS_row[,4] = c(str_to_upper(ENS_row$name))
mononew = merge(mono, ENS_row, by.x = "name")

总结一下

其实上面的代码仍然是很烂，毕竟是才学一个月的R，而且没有生物信息学基础不知道fasta格式的技巧，也不熟悉大量的成熟的包和工具，只能说自己的凭思维慢慢摸索。

首先需要找到对应物种的每个基因的cds序列的fasta文件地址并且下载，然后需要解析fasta文件计算每个基因的cds序列长度，最后筛选自己需要的基因即可。