234-给你几条蛋白序列,用R语言怎么比对?

刘小泽写于2021.2.2 今天我们以两条蛋白序列为例,来看看用R怎么比对,并且重点是要找到目标位点的比对结果。当然,今天这个R包不仅可以比对蛋白序列,还支持DNA, RNA

安装R包

rm(list=ls())
options(stringsAsFactors = F)

# BiocManager::install("msa") 
suppressMessages(library(msa))
library(stringr)

下载蛋白序列

首页(https://www.uniprot.org/)输入:基因名_物种,获取ID号

然后就可以下载到小鼠该基因的蛋白序列

注意网页链接的规律

https://www.uniprot.org/uniprot/P48754.fasta

同理,我们也获取到human的该基因蛋白序列: https://www.uniprot.org/uniprot/P38398.fasta

R中继续操作

官方文档在:https://bioconductor.org/packages/release/bioc/vignettes/msa/inst/doc/msa.pdf

# 读取序列
hg <- readAAStringSet('P38398.fasta')
mm <- readAAStringSet('P48754.fasta')
# 结果保存为AAStringSet对象。当然,还支持读入 DNAStringSet, and RNAStringSet对象

# 进行比对
mySequences=c(mm,hg)
myFirstAlignment <- msa(mySequences)

可以查看全部比对结果
print(myFirstAlignment, show="complete")

还可以输出为pdf
msaPrettyPrint(myFirstAlignment, output="pdf", showNames="none",showLogo="none", askForOverwrite=FALSE, verbose=FALSE)

其中感叹号!就是全部比对上的,*就是存在gap的位点

# 示例代码
msaPrettyPrint(myFirstAlignment, output="asis", y=c(164, 213),subset=c(1:6), showNames="none", showLogo="top",logoColors="rasmol", shadingMode="similar",showLegend=FALSE, askForOverwrite=FALSE)

image-20210202143631003

接下来才是重点

上面函数使用很简单,关键是我们得到了比对的结果,怎么提取我想要的位点呢

其实之前biostar上也有人提出了这个问题:https://www.biostars.org/p/467870/

也就是说,他之前感兴趣的位点是82位的AA,想要看看比对后这个AA对应着哪个AA。但是他发现,比对后原来的AA位置变成了95。所以不知道应该怎么去统计

这里给出我自己的理解

**首先明确,比对前后的差异在什么地方?**=》原来可以比对的当然能对应,但是不能比对的分了两种情况:一个是硬比对(比如D比对到了N);一个是软比对,即增加了gap(比如D 比对到. )。

麻烦就出在gap身上,原来序列没有,比对过程加上去。

因此,我们在统计的时候,就需要考虑进来gap的数量

这里以小鼠为query,人为target。也就是说,小鼠的位点可以按照之前下载的序列,但是人的位点,就需要加上gap。

以下面👇这个位点Q为例

  • 这里显示的一条序列长度是120,正因为混入了两个gap(蓝色的星号表示),所以原本238的位置,现在处于240【所以,比对结果右侧数字记录的就是:原始的位点所处的位置

  • 我们可以轻松得到小鼠原来第238位的位点:

    > pos=238
    > str_sub(as.character(mm),pos,pos)
    [1] "Q"
    
  • 那human呢?用同样的方法能行吗?

    > str_sub(as.character(hg),pos,pos)
    [1] "H"
    

    结果显然不对,它应该是Q,但这里显示H 。因此,我们应该把gap考虑进来

  • 在小鼠的第238位点之前,一共有几个gap?

    # 这里需要从比对结果(myFirstAlignment@unmasked)中提取,即myFirstAlignment@unmasked[1]
    > ngap=str_sub(as.character(myFirstAlignment@unmasked[1]),1,pos) %>% str_count('-')
    > ngap
    [1] 2
    
  • 这2个gap也符合我们的肉眼验证,因此human此时对应的位置应该是:

    > hg_site=str_sub(as.character(myFirstAlignment@unmasked[2]),pos+ngap,pos+ngap)
    > hg_site
    [1] "Q"
    
原以为这样就结束了,其实并没有

我测试了好多个,结果都对,但又随机挑选了一个502位点,发现了错误:

正确的应该是:S<=>T,但我得到的是:S<=>R

> pos=502
# 小鼠的
> str_sub(as.character(mm),pos,pos)
[1] "S"
# 人的
>         ngap=str_sub(as.character(myFirstAlignment@unmasked[1]),1,pos+ngap)%>% str_count('-')
> str_sub(as.character(myFirstAlignment@unmasked[2]),pos+ngap,pos+ngap)
[1] "R"

错误原因就是ngap计算错误:

# 此时计算的是
> ngap
# 但其实数一数,在502位点前,其实还有2个,因此是7个
但为什么gap会统计错误呢?

因为我们这里给出的pos=502,在比对结果中,是落在了真实502位置的前面,而且恰巧也落在了那2个新的gap的前面,所以没有统计上。这里我想了一种解决方案,就是增加一步while循环,来探索在比对结果的502位点之后,有没有新的gap出现

ngap=str_sub(as.character(myFirstAlignment@unmasked[1]),1,pos) %>% str_count('-')
new_ngap=str_sub(as.character(myFirstAlignment@unmasked[1]),1,pos+ngap) %>% str_count('-')
while (new_ngap != ngap) {
       ngap=new_ngap
       next
}

mm_site=str_sub(as.character(mm),pos,pos)
hg_site=str_sub(as.character(myFirstAlignment@unmasked[2]),pos+ngap,pos+ngap)
Yunze Liu
Yunze Liu
Bioinformatics Sharer

Co-founder of Bioinfoplanet(生信星球)

Next
Previous

Related