问题
我有多条mRNA序列，进行体外蛋白表达。由于成熟的蛋白是切除了信号肽的，我在体外表达蛋白时就需要去除信号肽，方法就是设计引物时就把信号肽那部分去掉。
一下子拿了九条序列，一个个手动去数很容易错。我就比较懒，想有没有办法批量去除信号肽呢？
只要想总会有的吧！
思路
1.预测信号肽的软件是signalIP，非常简单好用的在线工具。首先挨个预测下信号肽的氨基酸个数。这个似乎没什么捷径
2.首先要从genebank下载序列，查了一下这个可以批量！查到的结果是一个网址，ncbi的Batch Entrez。
https://www.ncbi.nlm.nih.gov/sites/batchentrez 可以批量查找并导出fasta序列。
3.然后要将fasta读取到R中，我本来已经打算当做字符串来做了，谁知刚好生信技能树发了一篇推文：《把fasta序列读入到R里面去》，喜出望外！
(被曾老师嫌弃说可以直接下载到R里,琢磨了一会没搞懂，可能涉及到了什么转录本。放弃治疗，反正我跑出来了)
4.然后就是按照字符串那样处理，只要知道了信号肽氨基酸个数，X3就是核苷酸个数了，设为n，然后用substring来提取第n+1位到最后一位就可以啦！
解决方案
（1）下载序列
将想要下载的mRNA序列号放进一个text文档ids.txt，每行一个号。
如果有错误的id号就会被列在Reject里
（2）导出fasta文件
点击上一页的蓝字会跳转到这一页
挨个选中（吐槽这里竟然没有全选按钮），拉到最后有个send to
导出
（3）读取到R
导出了一个名为sequence.fasta的多序列文件。
试了《把fasta序列读入到R里面去》列出的几种方法，有的因为NA报错了，最终发现对多序列比较友好的是seqinr。
BiocManager::install('seqinr')
library(seqinr)
fastafile <- read.fasta(file = "sequence.fasta",
as.string = TRUE,
forceDNAtolower = FALSE)
此时fasta被读取为一个列表了
想要得到一个数据框，一列是id，一列是原始的mRNA序列
seq_before <- vector("character")
id <- vector("character")
for (i in 1:length(fastafile)){
seq_before[i] =fastafile[[i]]
id[[i]] =names(fastafile[i])
}
df <- data.frame(id = id,seq_before = seq_before)
然后将信号肽数目加到最后一列，我还顺便加上了基因名。这个信息储存在一个xls文档中，三列，分别是基因名，id和信号肽数。数据不便公开我就不写了。
signal <- readxl::read_xlsx('id&signal peptide.xlsx')
seq <- left_join(signal,df,by = "id")
最重要的一步就是提取了，顺便导出文件
after <- vector()
for (i in 1:nrow(seq)) {
after[[i]]=substring(seq$seq_before[[i]],
seq$number_of_signal_peptide[[i]]+1,
str_length(seq_before)[[i]])
}
output <- mutate(seq,after =after)
library(xlsx)
write.xlsx(output,"no_signal_peptide.xlsx")
看着这些马赛克好想笑
网友评论