我有一个 fasta 文件,如下所示。我想将三个字母代码转换为一个字母代码。我怎样才能用 python 或 R 做到这一点?
>2ppo
ARGHISLEULEULYS
>3oot
METHISARGARGMET
期望的输出
>2ppo
RHLLK
>3oot
MHRRM
您的建议将不胜感激!
我有一个 fasta 文件,如下所示。我想将三个字母代码转换为一个字母代码。我怎样才能用 python 或 R 做到这一点?
>2ppo
ARGHISLEULEULYS
>3oot
METHISARGARGMET
期望的输出
>2ppo
RHLLK
>3oot
MHRRM
您的建议将不胜感激!
BioPython 已经内置了字典来帮助进行此类翻译。以下命令将向您显示可用字典的完整列表:
import Bio
help(Bio.SeqUtils.IUPACData)
您要查找的预定义字典:
Bio.SeqUtils.IUPACData.protein_letters_3to1['Ala']
使用字典查找单字母代码:
d = {'CYS': 'C', 'ASP': 'D', 'SER': 'S', 'GLN': 'Q', 'LYS': 'K',
'ILE': 'I', 'PRO': 'P', 'THR': 'T', 'PHE': 'F', 'ASN': 'N',
'GLY': 'G', 'HIS': 'H', 'LEU': 'L', 'ARG': 'R', 'TRP': 'W',
'ALA': 'A', 'VAL':'V', 'GLU': 'E', 'TYR': 'Y', 'MET': 'M'}
还有一个简单的函数可以将三个字母代码与整个字符串的一个字母代码匹配:
def shorten(x):
if len(x) % 3 != 0:
raise ValueError('Input length should be a multiple of three')
y = ''
for i in range(len(x) // 3):
y += d[x[3 * i : 3 * i + 3]]
return y
测试你的例子:
>>> shorten('ARGHISLEULEULYS')
'RHLLK'
这是在 R 中执行此操作的一种方法:
# Variables:
foo <- c("ARGHISLEULEULYS","METHISARGARGMET")
# Code maps:
code3 <- c("Ala", "Arg", "Asn", "Asp", "Cys", "Glu", "Gln", "Gly", "His",
"Ile", "Leu", "Lys", "Met", "Phe", "Pro", "Ser", "Thr", "Trp",
"Tyr", "Val")
code1 <- c("A", "R", "N", "D", "C", "E", "Q", "G", "H", "I", "L", "K",
"M", "F", "P", "S", "T", "W", "Y", "V")
# For each code replace 3letter code by 1letter code:
for (i in 1:length(code3))
{
foo <- gsub(code3[i],code1[i],foo,ignore.case=TRUE)
}
结果是 :
> foo
[1] "RHLLK" "MHRRM"
请注意,我更改了变量名称,因为变量名称不允许以 R 中的数字开头。
>>> src = "ARGHISLEULEULYS"
>>> trans = {'ARG':'R', 'HIS':'H', 'LEU':'L', 'LYS':'K'}
>>> "".join(trans[src[x:x+3]] for x in range(0, len(src), 3))
'RHLLK'
您只需要将其余条目添加到trans
字典中。
编辑:
要完成其余部分trans
,您可以执行此操作。文件table
:
Ala A
Arg R
Asn N
Asp D
Cys C
Glu E
Gln Q
Gly G
His H
Ile I
Leu L
Lys K
Met M
Phe F
Pro P
Ser S
Thr T
Trp W
Tyr Y
Val V
阅读:
trans = dict((l.upper(), s) for l, s in
[row.strip().split() for row in open("table").readlines()])
您可以尝试查看并安装Biopython,因为您正在解析 .fasta 文件,然后转换为单字母代码。不幸的是,Biopython 只有函数 seq3(在包 Bio::SeqUtils 中),它与你想要的相反。IDLE 中的示例输出:
>>>seq3("MAIVMGRWKGAR*")
>>>'MetAlaIleValMetGlyArgTrpLysGlyAlaArgTer'
不幸的是,没有“seq1”功能(还...),但我认为这可能对您将来有所帮助。至于你的问题,Junuxx 是正确的。创建一个字典并使用 for 循环以三个为一组读取字符串并进行翻译。这是一个与他提供的功能类似的功能,它包罗万象,也可以处理小写字母。
def AAcode_3_to_1(seq):
'''Turn a three letter protein into a one letter protein.
The 3 letter code can be upper, lower, or any mix of cases
The seq input length should be a factor of 3 or else results
in an error
>>>AAcode_3_to_1('METHISARGARGMET')
>>>'MHRRM'
'''
d = {'CYS': 'C', 'ASP': 'D', 'SER': 'S', 'GLN': 'Q', 'LYS': 'K',
'ILE': 'I', 'PRO': 'P', 'THR': 'T', 'PHE': 'F', 'ASN': 'N',
'GLY': 'G', 'HIS': 'H', 'LEU': 'L', 'ARG': 'R', 'TRP': 'W', 'TER':'*',
'ALA': 'A', 'VAL':'V', 'GLU': 'E', 'TYR': 'Y', 'MET': 'M','XAA':'X'}
if len(seq) %3 == 0:
upper_seq= seq.upper()
single_seq=''
for i in range(len(upper_seq)/3):
single_seq += d[upper_seq[3*i:3*i+3]]
return single_seq
else:
print("ERROR: Sequence was not a factor of 3 in length!")
Biopython 有一个很好的解决方案
>>> from Bio.PDB.Polypeptide import *
>>> three_to_one('ALA')
'A'
对于您的示例,我将通过这一班轮解决它
>>> from Bio.PDB.Polypeptide import *
>>> str3aa = 'ARGHISLEULEULYS'
>>> "".join([three_to_one(aa3) for aa3 in [ "".join(g) for g in zip(*(iter(str3aa),) * 3)]])
>>> 'RHLLK'
他们可能会批评我这种类型的一个班轮:),但在我内心深处,我仍然爱着 PERL。
使用 R:
convert <- function(l) {
map <- c("A", "R", "N", "D", "C", "E", "Q", "G", "H", "I",
"L", "K", "M", "F", "P", "S", "T", "W", "Y", "V")
names(map) <- c("ALA", "ARG", "ASN", "ASP", "CYS", "GLU", "GLN",
"GLY", "HIS", "ILE", "LEU", "LYS", "MET", "PHE",
"PRO", "SER", "THR", "TRP", "TYR", "VAL")
sapply(strsplit(l, "(?<=[A-Z]{3})", perl = TRUE),
function(x) paste(map[x], collapse = ""))
}
convert(c("ARGHISLEULEULYS", "METHISARGARGMET"))
# [1] "RHLLK" "MHRRM"
# install.packages("seqinr")
# source("https://bioconductor.org/biocLite.R")
# biocLite("iPAC")
library(seqinr)
library(iPAC)
#read in file
fasta = read.fasta(file = "test_fasta.fasta", seqtype = "AA", as.string = T, set.attributes = F)
#split string
n = 3
fasta1 = lapply(fasta, substring(x,seq(1,nchar(x),n),seq(n,nchar(x),n)))
#convert the three letter code for each element in the list
fasta2 = lapply(fasta1, function(x) paste(sapply(x, get.SingleLetterCode), collapse = ""))
# > fasta2
# $`2ppo`
# [1] "RHLLK"
#
# $`3oot`
# [1] "MHRRM"
my %aa_hash=(
Ala=>'A',
Arg=>'R',
Asn=>'N',
Asp=>'D',
Cys=>'C',
Glu=>'E',
Gln=>'Q',
Gly=>'G',
His=>'H',
Ile=>'I',
Leu=>'L',
Lys=>'K',
Met=>'M',
Phe=>'F',
Pro=>'P',
Ser=>'S',
Thr=>'T',
Trp=>'W',
Tyr=>'Y',
Val=>'V',
Sec=>'U', #http://www.uniprot.org/manual/non_std;Selenocysteine (Sec) and pyrrolysine (Pyl)
Pyl=>'O',
);
while(<>){
chomp;
my $aa=$_;
warn "ERROR!! $aa invalid or not found in hash\n" if !$aa_hash{$aa};
print "$aa\t$aa_hash{$aa}\n";
}
使用此 perl 脚本将三元组 aa 代码转换为单字母代码。
Python 3 解决方案。
在我的工作中,令人烦恼的部分是氨基酸代码可以引用经常出现在 PDB/mmCIF 文件中的修改过的代码,例如
'Tih'-->'A'。
所以映射可以超过22对。Python 中的第 3 方工具,例如
Bio.SeqUtils.IUPACData.protein_letters_3to1
无法处理。我最简单的解决方案是使用http://www.ebi.ac.uk/pdbe-srv/pdbechem查找映射并将异常映射添加到我自己的函数中的字典,每当我遇到它们时。
def three_to_one(three_letter_code):
mapping = {'Aba':'A','Ace':'X','Acr':'X','Ala':'A','Aly':'K','Arg':'R','Asn':'N','Asp':'D','Cas':'C',
'Ccs':'C','Cme':'C','Csd':'C','Cso':'C','Csx':'C','Cys':'C','Dal':'A','Dbb':'T','Dbu':'T',
'Dha':'S','Gln':'Q','Glu':'E','Gly':'G','Glz':'G','His':'H','Hse':'S','Ile':'I','Leu':'L',
'Llp':'K','Lys':'K','Men':'N','Met':'M','Mly':'K','Mse':'M','Nh2':'X','Nle':'L','Ocs':'C',
'Pca':'E','Phe':'F','Pro':'P','Ptr':'Y','Sep':'S','Ser':'S','Thr':'T','Tih':'A','Tpo':'T',
'Trp':'W','Tyr':'Y','Unk':'X','Val':'V','Ycm':'C','Sec':'U','Pyl':'O'} # you can add more
return mapping[three_letter_code[0].upper() + three_letter_code[1:].lower()]
另一种解决方案是在线检索映射(但 url 和 html 模式可能会随着时间而改变):
import re
import urllib.request
def three_to_one_online(three_letter_code):
url = "http://www.ebi.ac.uk/pdbe-srv/pdbechem/chemicalCompound/show/" + three_letter_code
with urllib.request.urlopen(url) as response:
single_letter_code = re.search('\s*<td\s*>\s*<h3>One-letter code.*</h3>\s*</td>\s*<td>\s*([A-Z])\s*</td>', response.read().decode('utf-8')).group(1)
return single_letter_code
为了简单起见,这里我直接使用 re 而不是 html 解析器。
希望这些能有所帮助。
对于那些在 2017 年及以后登陆这里的人:
这是一个单行 Linux bash 命令,用于将蛋白质氨基酸三字母代码转换为文本文件中的单字母代码。我知道这不是很优雅,但我希望这可以帮助搜索相同内容并希望使用单行命令的人。
sed 's/ALA/A/g;s/CYS/C/g;s/ASP/D/g;s/GLU/E/g;s/PHE/F/g;s/GLY/G/g;s/HIS/H/g;s/HID/H/g;s/HIE/H/g;s/ILE/I/g;s/LYS/K/g;s/LEU/L/g;s/MET/M/g;s/ASN/N/g;s/PRO/P/g;s/GLN/Q/g;s/ARG/R/g;s/SER/S/g;s/THR/T/g;s/VAL/V/g;s/TRP/W/g;s/TYR/Y/g;s/MSE/X/g' < input_file_three_letter_code.txt > output_file_single_letter_code.txt
上面原始问题的解决方案,作为单个命令行:
sed 's/.\{3\}/& /g' | sed 's/ALA/A/g;s/CYS/C/g;s/ASP/D/g;s/GLU/E/g;s/PHE/F/g;s/GLY/G/g;s/HIS/H/g;s/HID/H/g;s/HIE/H/g;s/ILE/I/g;s/LYS/K/g;s/LEU/L/g;s/MET/M/g;s/ASN/N/g;s/PRO/P/g;s/GLN/Q/g;s/ARG/R/g;s/SER/S/g;s/THR/T/g;s/VAL/V/g;s/TRP/W/g;s/TYR/Y/g;s/MSE/X/g' | sed 's/ //g' < input_file_three_letter_code.txt > output_file_single_letter_code.txt
解释:
[1]sed 's/.\{3\}/& /g'
将拆分序列。它会在每 3 个字母后添加一个空格。
[2] 管道中的第二个 'sed'
命令将获取上面的输出并转换为单字母代码。s/XYZ/X/g;
为该命令添加任何非标准残基。
[3] 第三个 ' sed
' 命令,sed 's/ //g'
将删除空格。