给你几条蛋白序列，用R语言怎么比对？
安装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的位点
当然，如果比对的序列很多，还可以增加logo
# 示例代码
msaPrettyPrint(myFirstAlignment, output="asis", y=c(164, 213),subset=c(1:6), showNames="none", showLogo="top",logoColors="rasmol", shadingMode="similar",showLegend=FALSE, askForOverwrite=FALSE)
接下来才是重点
上面函数使用很简单，关键是我们得到了比对的结果，怎么提取我想要的位点呢
其实之前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)腾讯云开发者
