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

Question

我想将蛋白质数据库中的PDB文件与Cosmic或Uniprot中显示的蛋白质的标准AA序列进行匹配。具体地说，我需要做的是从pdb文件中提取主干中的碳α原子及其xyz位置。我还需要提取它们在蛋白质序列中的实际顺序。对于structure 3GFT (Kras - Uniprot注册号P01116)，这很简单，我只需获取ResSeq号即可。然而，对于其他一些蛋白质，我不知道这是如何实现的。

例如，对于结构(2ZHQ) (蛋白质F2 - Uniprot登录号P00734)，序列具有重复的数字"1“和"14”的ResSeq编号，并且仅在Icode条目上不同。此外，icode条目不是按词法顺序排列的，因此很难判断提取的顺序。

如果考虑结构3V5Q (Uniprot注册号Q16288)，情况会变得更糟。对于大多数蛋白质，ResSeq数与来自COSMIC或UNIPROT等来源的实际氨基酸相匹配。然而，在位置711之后，它跳到位置730。当查看备注465 (缺少的原子)时，它显示对于链A，缺少726-729。然而，在将其与蛋白质进行匹配后，这些氨基酸实际上是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，因为他是将PDB结构链接到uniprot.org规范序列的常驻UniProt专家。

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()函数。