从蛋白质数据库到Cosmic或Uniprot的蛋白质序列比对

Question

我想将蛋白质数据库中的PDB文件与Cosmic或Uniprot中显示的蛋白质的典型AA序列相匹配。具体来说，我需要做的是从pdb文件中拉出主干中的碳原子和它们的xyz位置。我还需要在蛋白质序列中提取他们的实际顺序。对于结构3GFT（克拉斯 - Uniprot登录号P01116），这很容易，我可以采用ResSeq编号。但是，对于其他一些蛋白质，我无法弄清楚这是可能的。例如，对于结构（2ZHQ）（蛋白质F2- Uniprot登录号P00734），Seqres具有针对数字“1”和“14”重复的ResSeq数字，并且仅在Icode条目中有所不同。进一步icode条目不是在词法顺序，因此很难说出要提取的顺序。

如果考虑结构3V5Q（Uniprot登录号Q16288），情况会更糟糕。对于大多数蛋白质，ResSeq编号与来自COSMIC或UNIPROT等来源的实际氨基酸相匹配。在位置711之后，它跳到位置730.当查看REMARK 465（缺失的原子）时，它显示对于链A，缺失了726-729。然而，匹配到蛋白质后，那些AA实际上是712-715。

我附上了一些简单的3GFT示例代码，但如果有人是pdb文件的专家，并且可以帮助我解决其余的问题，我会非常感激。

library(gdata) 

#answer<- get.positions("http://www.pdb.org/pdb/files/2ZHQ.pdb","L") 
answer<- get.positions("http://www.pdb.org/pdb/files/3GFT.pdb","A") 


#This function reads a pdb file and returns the appropriate data structure 
get.positions <- function(sourcefile, chain_required = "A"){ 
N <- 10^5 
AACount <- 0 
positions = data.frame(Residue=rep(NA, N),AtomCount=rep(0, N),SideChain=rep(NA, N),XCoord=rep(0, N),YCoord=rep(0, N),ZCoord=rep(0, N),stringsAsFactors=FALSE)  

visited = list() 
filedata <- readLines(sourcefile, n= -1) 
for(i in 1: length(filedata)){ 
input = filedata[i] 
id = substr(input,1,4) 
if(id == "ATOM"){ 
    type = substr(input,14,15) 
    if(type == "CA"){ 
    resSerial = strtoi(substr(input, 7,11)) 
    residue = substr(input,18,20) 
    type_of_chain = substr(input,22,22) 
    resSeq = strtoi(substr(input, 23,26)) 
    altLoc = substr(input,17,17) 

    if(resSeq >=1){ #does not include negative residues 
     if(type_of_chain == chain_required && !(resSerial %in% visited) && (altLoc == " " || altLoc == "A")){ 
     visited <- c(visited, resSerial) 
     AACount <- AACount + 1 
     position_string =list() 
     position_string[[1]]= as.numeric(substr(input,31,38)) 
     position_string[[2]] =as.numeric(substr(input,39,46)) 
     position_string[[3]] =as.numeric(substr(input,47,54)) 
     #print(input) 
     positions[AACount,]<- c(residue, resSeq, type_of_chain, position_string[[1]], position_string[[2]], position_string[[3]]) 
     } 
    } 
    } 
} 

    } 
    positions<-positions[1:AACount,] 
    positions[,2]<- as.numeric(positions[,2]) 
    positions[,4]<- as.numeric(positions[,4]) 
    positions[,5]<- as.numeric(positions[,5]) 
    positions[,6]<- as.numeric(positions[,6]) 
    return (positions) 
} 

Answer 1

写作时你可能想这个问题转移到www.biostars.org和写入help@uniprot.org（你知道，这些序列在数据库级别就足够了联系？）在任何情况下， help@uniprot.org请求Jules Jacobsen，因为他是UniProt常驻专家，负责将PDB结构链接到uniprot.org规范序列。

Answer 2

这是一种方法。它需要bio3d R软件包（http://thegrantlab.org/bio3d/），并安装肌肉对齐可执行文件。我按照这里的'其他实用程序'的说明http://thegrantlab.org/bio3d/tutorials/installing-bio3d

下面的示例代码似乎适用于您列出的三种情况。

library(bio3d) 

## Download PDB file with given 'id' (can also read from online) 
id <- "3GFT" #"3V5Q" 
file.name <- get.pdb(id) 
pdb <- read.pdb(file.name) 
pdb.seq <- pdbseq(pdb, atom.select(pdb, chain="A", elety="CA")) 

## Get UniProt identifier and sequence 
pdb.ano <- pdb.annotate(id, "db_id") 
uni.seq <- get.seq(pdb.ano) 

## Align sequences to define corespondences 
aln <- seqaln(seqbind(pdb.seq, uni.seq), id=c(file.name, unlist(pdb.ano))) 

## Read aligned coordinate data (all the info you want is in here) 
pdbs <- read.fasta.pdb(aln) 

answer2 <- cbind(1:ncol(pdbs$ali), t(pdbs$ali), 
      pdbs$resno[1,], matrix(pdbs$xyz[1,], ncol=3, byrow=T)) 

head(answer2) 

[1,] "1" "M"  "M" "1" "62.935" "97.579" "30.223" 
[2,] "2" "T"  "T" "2" "63.155" "95.525" "27.079" 
[3,] "3" "E"  "E" "3" "65.289" "96.895" "24.308" 
[4,] "4" "Y"  "Y" "4" "64.899" "96.22" "20.615" 
[5,] "5" "K"  "K" "5" "67.593" "96.715" "18.023" 
[6,] "6" "L"  "L" "6" "65.898" "97.863" "14.816" 

，如果你希望你的氨基酸3字母代码上市有一个在bio3d一个aa321（）函数。