使用 Biostrings 编写 R 程序以将核苷酸序列转换为蛋白质序列答案

【问题标题】：Writing R program using Biostrings to convert nucleotide sequences to protein sequence使用 Biostrings 编写 R 程序以将核苷酸序列转换为蛋白质序列
【发布时间】：2021-08-10 04:03:11
【问题描述】：

我有一个包含许多不同序列的 .fasta 文件。我的目标是使用 Biostrings 包将每个单独的序列转换为它的氨基酸序列。 .fasta 文件如下所示：

>Sequence 1
AAATTTGGGCCC
>Sequence 2
TTTGGGCCCAAA

任何帮助将不胜感激。谢谢。

【问题讨论】：

标签： r bioinformatics

【解决方案1】：

translate function 会做你想做的事：

library(Biostrings)

`Sequence 1` <- DNAString("AAATTTGGGCCC")
`Sequence 2` <- DNAString("TTTGGGCCCAAA")

seq_1 <- translate(`Sequence 1`, no.init.codon=TRUE)
seq_1
#> 4-letter AAString object
#> seq: KFGP

seq_2 <- translate(`Sequence 2`, no.init.codon=TRUE)
seq_2
#> 4-letter AAString object
#> seq: FGPK

读入整个fasta文件：

seqs <- Biostrings::readDNAStringSet("file.fasta", format = "fasta", use.names = TRUE)
seqs_translated <- translate(seqs, no.init.codon = TRUE)
seqs_translated
#> AAStringSet object of length 2:
#>     width seq                               names               
#> [1]     4 KFGP                              Sequence 1
#> [2]     4 FGPK                              Sequence 2

编辑

您翻译 fasta 文件的问题是序列使用 'full' alphabet，而不仅仅是 ATCG - 您有“无呼叫”（“N”）、间隙（“-”）和矛盾/未解决的呼叫，例如“K”（鸟嘌呤或胸腺嘧啶）。我使用 sed 找到了这些：

grep -v ">" SEQUENCE_orf1ab.fasta | sed 's/[ATCG]//g' | sed '/^$/d'
# explanation: remove lines beginning with ">"
# then remove all A/T/C/G's and blank lines
# what you have left is causing the "not a base" error

如果您使用例如删除这些“非碱基”碱基

sed '/^>/! s/[-NYRKW]//g' SEQUENCE_orf1ab.fasta > test.fasta 
#explanation: in lines not beginning with ">", substitute all of the characters "-NYRKW" with nothing (i.e. delete them)

然后文件被翻译没有问题：

seqs <- Biostrings::readDNAStringSet("test.fasta", format = "fasta", use.names = TRUE)
seqs_translated <- translate(seqs, no.init.codon = TRUE)
seqs_translated

#>  AAStringSet object of length 91:
#>       width seq                              names               
#>  [1]  7095 MESLVPGFNEKTHVQ...KTTELLFLVMFLLT MZ505877.1 |Sever...
#>  [2]  7095 MESLVPGFNEKTHVQ...KTTELLFLVMFLLT MZ020653.1 |Sever...
#>  [3]  7092 MESLVPGFNEKTHVQ...KTTELLFLVMFLLT MW988268.1 |Sever...
#>  [4]  7095 MESLVPGFNEKTHVQ...KTTELLFLVMFLLT MW928277.1 |Sever...
#>  [5]  7095 MESLVPGFNEKTHVQ...KTTELLFLVMFLLT MW885875.1 |Sever...
#>  ...   ... ...
#> [87]  7095 MESLVPGFNEKTHVQ...KTTELLFLVMFLLT MN996529.1 |Sever...
#> [88]  7095 MESLVPGFNEKTHVQ...KTTELLFLVMFLLT MN996530.1 |Sever...
#> [89]  7095 MESLVPGFNEKTHVQ...KTTELLFLVMFLLT MN996531.1 |Sever...
#> [90]  7094 MESLVPGFNEKTHVQ...KTTELLFLVMFLLT MN988713.1 |Sever...
#> [91]  7095 MESLVPGFNEKTHVQ...KTTELLFLVMFLLT MN975262.1 |Sever...

【讨论】：

感谢您的回复。我遇到的问题是序列比这长得多。我希望从 fasta 中加载它们，而不是手动输入每个。
我不断收到以下错误：Error in .Call2("DNAStringSet_translate", x, skip_code, dna_codes[codon_alphabet], : in 'x[[3]]': not a base at pos 11023 我不知道如何解决这个问题，因为我检查了序列并且位置 11023 有碱基 C。谢谢您的帮助。
每当我运行sed '/^>/! s/[-NYRKW]//g' SEQUENCE_orf1ab.fasta > test.fasta 命令时，我都会收到以下错误：Error: unexpected string constant in "sed '/^>/! s/[-NYRKW]//g'" 我不确定如何解决此问题，我们将不胜感激。
很奇怪。所以命令grep -v ">" SEQUENCE_orf1ab.fasta | sed 's/[ATCG]//g' | sed '/^$/d' 起作用了，但是当您尝试运行sed '/^>/! s/[-NYRKW]//g' SEQUENCE_orf1ab.fasta > test.fasta 时会出现意外的字符串常量错误？
第一个命令给了我以下错误：Error: unexpected string constant in "grep -v ">""