售后

用户工具

站点工具

您的足迹: • mmseqs2
mmseqs2

目录

MMseqs2：超快和敏感的序列搜索和聚类套件

1.软件介绍

MMseqs2（Many-against-Many序列搜索）是一个软件套件，用于搜索和集群巨大的蛋白质和核苷酸序列集。MMseqs2是开源GPL许可软件，在Linux、MacOS和（作为测试版，通过cygwin）Windows的C++中实现。该软件旨在在多个内核和服务器上运行，并具有非常好的可扩展性。MMseqs2的运行速度比BLAST快10000倍。以其速度的100倍，它的灵敏度几乎相同。它可以以与PSI-BLAST相同的灵敏度以400多倍的速度执行配置文件搜索。

2.聚类

mmseqs easy-cluster examples/DB.fasta clusterRes tmp --min-seq-id 0.5 -c 0.8 --cov-mode 1

运行时随着输入大小而线性缩放。此模式推荐用于大型数据集。 

mmseqs easy-linclust examples/DB.fasta clusterRes tmp

3.比对

mmseqs easy-search examples/QUERY.fasta examples/DB.fasta alnRes.m8 tmp

也可以预先计算目标数据库的索引。这减少了重复搜索同一数据库时的开销。 

mmseqs createdb examples/DB.fasta targetDB
mmseqs createindex targetDB tmp
mmseqs easy-search examples/QUERY.fasta targetDB alnRes.m8 tmp

databases工作流程为许多公共参考数据库提供了下载和设置程序，如Uniref、NR、NT、PFAM等（请参阅下载数据库 https://github.com/soedinglab/mmseqs2/wiki#downloading-databases ）。例如，要下载并搜索包含Swiss-Prot参考蛋白的数据库，请运行： 

mmseqs databases UniProtKB/Swiss-Prot swissprot tmp
mmseqs easy-search examples/QUERY.fasta swissprot alnRes.m8 tmp

4.物种注释

easy-taxonomy工作流程可用于分配序列分类标签。它对具有分类信息（seqTaxDb）的序列数据库进行搜索，根据不同的策略（根据–lca-mode）选择最具代表性的对齐目标序列集，并计算其中最低的共同祖先。 

mmseqs createdb examples/DB.fasta targetDB
mmseqs createtaxdb targetDB tmp
mmseqs createindex targetDB tmp
mmseqs easy-taxonomy examples/QUERY.fasta targetDB alnRes tmp

实战：如何下载搭建RefSeq数据库中的病毒Viral相关数据库内容，并进行注释

一、构建数据库

1.激活环境 

source /TJPROJ1/META_ASS/soft/mmseqs2/../anaconda3/bin/activate /TJPROJ1/META_ASS/soft/mmseqs2

2.下载需要进行构建数据库的序列内容，可以为核酸序列，可以为蛋白序列 

wget -c https://ftp.ncbi.nlm.nih.gov/refseq/release/viral/viral.1.1.genomic.fna.gz

3.构建序列索引数据库 

mmseqs createdb viral.1.1.genomic.fna refseq_viral.1.1
mmseqs createindex refseq_viral.1.1 refseq_viral.1.1_tmp/ --search-type 3

4.下载注释信息文件,Refseq数据库中可以从nt库中进行查询，因此需要nt、nr库中相关注释names.dmp或nodes.dmp文件与accession2taxid文件 

wget https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/accession2taxid/nucl_gb.accession2taxid.gz

5.构建mapping文件，该步骤的目的是为了获得Refseq数据库中的序列和taxid的对应关系 

python3 getname.py
python3 mapping.py

其中getname.py的脚本内容是： 

import sys
import io

# 确保 stdout 使用 UTF-8 编码
sys.stdout = io.TextIOWrapper(sys.stdout.buffer, encoding='utf-8')

def extract_fna_names(fna_file, output_file):
    with open(fna_file, 'r') as file:
        names = []
        for line in file:
            if line.startswith('>'):
                names.append(line[1:].strip().split(' ')[0])
    
    with open(output_file, 'w', encoding='utf-8') as outfile:
        for name in names:
            outfile.write(f"{name}\n")

# 使用示例
fna_file = 'viral.1.1.genomic.fna'
output_file = 'name.list'
extract_fna_names(fna_file, output_file)

mapping.py的脚本内容是： 

import sys
import io

# 确保 stdout 使用 UTF-8 编码
sys.stdout = io.TextIOWrapper(sys.stdout.buffer, encoding='utf-8')

def map_accession_to_taxid(name_file, accession_file, output_file):
    # 读取 name.list 文件
    with open(name_file, 'r') as nf:
        names = set(line.strip() for line in nf)

    # 创建映射字典
    mapping = {}
    with open(accession_file, 'r') as af:
        for line in af:
            columns = line.strip().split('\t')
            if len(columns) >= 4:
                taxid = columns[1]      # 第二列
                additional_info = columns[2]  # 第三列
                # 仅在名称集合中查找
                if taxid in names:
                    mapping[taxid] = additional_info

    # 将结果写入 mapping.txt
    with open(output_file, 'w', encoding='utf-8') as of:
        for taxid, additional_info in mapping.items():
            of.write(f"{taxid}\t{additional_info}\n")

# 使用示例
name_file = 'name.list'  # 名称列表文件
accession_file = 'nucl_gb.accession2taxid'  # accession 到 taxid 的映射文件
output_file = 'mapping.txt'  # 输出文件名
map_accession_to_taxid(name_file, accession_file, output_file)

6.构建注释索引库 

mkdir ncbi-taxdump && cd ncbi-taxdump
wget https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxdump.tar.gz
tar xzvf taxdump.tar.gz
cd -

7.构建mmseqs专用的数据库格式，并输入指定名称 

mmseqs createtaxdb refseq_viral.1.1 refseq_viral.1.1_tmp/ --ncbi-tax-dump ncbi-taxdump --tax-mapping-file taxidmapping
#检查有多少标识符具有以下代码的指定分类群
awk 'FNR==NR{f[$1]=$2; next} $1 in f{ print $2" has taxid "f[$1];} !($1 in f){print $2" has no taxid";} ' refseq_viral.1.1_mapping refseq_viral.1.1.lookup

二、物种注释

1.激活环境 

source /TJPROJ1/META_ASS/soft/mmseqs2/../anaconda3/bin/activate /TJPROJ1/META_ASS/soft/mmseqs2

2.将需要被比对的、注释的序列转换为数据库格式，如果不转化，则下面的search和taxonomy提示数据格式不可用，也可以使用easy-模式 

mmseqs createdb CK.1.1.fasta queryDB

3.序列比对 

mmseqs search queryDB /TJPROJ1/META_ASS/soft/mmseqs2/database/test/refseq_viral.1.1 resultDB tmp

4.物种注释 

mmseqs taxonomy queryDB /TJPROJ1/META_ASS/soft/mmseqs2/database/test/refseq_viral.1.1 taxonomyResult tmp

5.导出注释结果 

mmseqs createtsv queryDB taxonomyResult taxonomyResult.tsv
mmseqs taxonomyreport /TJPROJ1/META_ASS/soft/mmseqs2/database/test/refseq_viral.1.1 taxonomyResult taxonomyResult_report
mmseqs taxonomyreport /TJPROJ1/META_ASS/soft/mmseqs2/database/test/refseq_viral.1.1 taxonomyResult report.html --report-mode 1

拓展：如何将Silva数据库转换为NCBI node/names 的格式

拓展1.普通列表项目为SILVA创建一个seqTaxDB

如果想根据SILVA数据库对核糖体RNA（16S、18S、SSU）序列进行分类，您可以通过databases工作流程进行下载： 

mmseqs databases SILVA silvadb tmp

要了解SILVA MMseqs2分类序列数据库是如何构建的，请查看下面的示例脚本。我们首先需要从SILVA分类法中创建一个像分类法一样的NCBI。 

# build name.dmp, node.dmp from SILVA taxonomy
mkdir taxonomy/ && cd "$_"
wget ftp://ftp.arb-silva.de/current/Exports/taxonomy/tax_slv_ssu_*.txt.gz
buildNCBITax=$(cat << 'EOF'
BEGIN{
  ids["root"]=1; 
  print "1\t|\t1\t|\tno rank\t|\t-\t|" > "nodes.dmp"
  print "1\t|\troot\t|\t-\t|\tscientific name\t|" > "names.dmp";
} 
{ n=split($1, a, ";"); 
  gsub("domain", "superkingdom", $3);
  ids[$1]=$2;
  gsub(/[^,;]*;$/,"",$1); 
  pname=$1; 
  if(n==2){ 
    pname="root"
  }
  pid=ids[pname];  
  printf("%s\t|\t%s\t|\t%s\t|\t-\t|\n", $2, pid, $3) > "nodes.dmp";
  printf("%s\t|\t%s\t|\t-\t|\tscientific name\t|\n",$2,a[n-1]) >"names.dmp"; 
}
EOF
)
awk -F'\t' "$buildNCBITax" <(gunzip -c tax_slv_ssu_*.txt.gz)
touch merged.dmp 
touch delnodes.dmp
cd .. 

# create the database SILVA database from Nr99 fasta
wget ftp://ftp.arb-silva.de/current/Exports/SILVA_*_SSURef_NR99_tax_silva_full_align_trunc.fasta.gz
mmseqs createdb SILVA_*_SSURef_NR99_tax_silva_full_align_trunc.fasta.gz SILVA_DB

# add taxonomy to SILVA_DB
wget ftp://ftp.arb-silva.de/current/Exports/taxonomy/tax_slv_ssu_*.acc_taxid
mmseqs createtaxdb SILVA_DB tmp --ncbi-tax-dump taxonomy/ --tax-mapping-file tax_slv_ssu_*.acc_taxid

拓展2.为GTDB创建一个seqTaxDB

基因组分类数据库（GTDB）是一个系统发育一致的数据库，它重新定义了分类树。MMseqs2可以搜索GTDB，但它需要一些预处理步骤。 

# build name.dmp, node.dmp from GTDB taxonomy
mkdir taxonomy/ && cd "$_"
wget https://data.ace.uq.edu.au/public/gtdb/data/releases/latest/ssu.fna
buildNCBITax=$(cat << 'EOF'
BEGIN{
  ids["root"]=1; 
  rank["c"]="class"l
  rank["d"]="superkingdom";
  rank["f"]="family";
  rank["g"]="genus";
  rank["o"]="order";
  rank["p"]="phylum";
  rank["s"]="species";
  taxCnt=1;
  print "1\t|\t1\t|\tno rank\t|\t-\t|" > "nodes.dmp"
  print "1\t|\troot\t|\t-\t|\tscientific name\t|" > "names.dmp";
} 
/^>/{
  str=$2
  for(i=3; i<=NF; i++){ str=str" "$i} 
  n=split(str, a, ";"); 
  prevTaxon=1;
  for(i = 1; i<=n; i++){ 
    if(a[i] in ids){
      prevTaxon=ids[a[i]];
    }else{
      taxCnt++;
      split(a[i],b,"_");
      printf("%s\t|\t%s\t|\t%s\t|\t-\t|\n", taxCnt, prevTaxon, rank[b[1]]) > "nodes.dmp";
      printf("%s\t|\t%s\t|\t-\t|\tscientific name\t|\n", taxCnt, b[3]) >"names.dmp"; 
      ids[a[i]]=taxCnt;
      prevTaxon=ids[a[i]];
    }
  }
  gsub(">", "", $1);
  printf("%s\t%s\n", $1, ids[a[n]]) > "mapping";
}
EOF
)
awk -F'\\[loc' '{ print $1}' ssu.fna | awk "$buildNCBITax" 
touch merged.dmp 
touch delnodes.dmp
cd .. 

mmseqs createdb ssu.fna ssu

# add taxonomy to GTDB
mmseqs createtaxdb ssu tmp --ncbi-tax-dump taxonomy/ --tax-mapping-file taxonomy/mapping

拓展3.通过序列数据库的手动注释创建seqTaxDB

这是一个如何用分类信息手动注释序列数据库的示例。该示例使用Uniprot标识符。 作为第一步，将FAST[A/Q]文件转换为mmseqs序列数据库，使用createdb 

# Turn the sequences into a MMseqs2 database (this also creates sequenceDB.lookup)
# Skip this step if you already created a database
mmseqs createdb sequence.fasta sequenceDB

createdb生成一个制表符分隔的sequenceDB.lookup文件，其中包含数字-db-id，Accession（例如Uniprot Accession Q6GZX4）和文件。ID是从输入数据库的标头中解析的（请参阅从标头解析的id）。 

0 Q6GZX4  0
1 Q6GZX3  0
2 Q197F8  0
3 P0A031  0
4 Q197F7  0

作为下一步，我们创建一个制表符分隔的映射，每个目标数据库标识符都映射到NCBI分类群标识符。映射文件应采用Accession numeric-ncbi-tax-id的格式。 

Q6GZX4    654924
Q6GZX3    654924
Q197F8    345201
Q197F7    345201

以下是如何将Uniprot映射文件转换为制表符分隔的映射文件的示例。其他数据库可以参考我上面写的Refseq数据库构建 

# The taxidmapping file should be in the format 
# Accession numeric-ncbi-tax-id
# Q6GZX4    654924
# Q6GZX3    654924
# Q197F8    345201
# Q197F7    345201

# e.g. download the uniprot mapping file and convert it to the taxidmapping mapping format
URL="ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/idmapping/idmapping.dat.gz"
wget -nv -O - "$URL" | zcat | awk '$2 == "NCBI_TaxID" {print $1"\t"$3 }' > taxidmapping

我们需要NCBI分类法taxdump.tar.gz。它在NCBI FTP服务器上可用： 

mkdir ncbi-taxdump && cd ncbi-taxdump
wget ftp://ftp.ncbi.nih.gov/pub/taxonomy/taxdump.tar.gz
tar xzvf taxdump.tar.gz
cd -

作为最后一步，我们现在可以使用createtaxdb来注释我们的序列数据库。 

# now we can use createtaxdb with our own mapping.
mmseqs createtaxdb sequenceDB tmp --ncbi-tax-dump ncbi-taxdump --tax-mapping-file taxidmapping

可以检查有多少标识符具有以下代码的指定分类群 

awk 'FNR==NR{f[$1]=$2; next} $1 in f{ print $2" has taxid "f[$1];} !($1 in f){print $2" has no taxid";} ' sequenceDB_mapping sequenceDB.lookup

MMseqs2超快和敏感的蛋白质搜索，Linclust在线性时间内聚集巨大的蛋白质序列集和Plass蛋白水平组装，以增加复杂元基因组的蛋白质序列回收。

mmseqs2.txt · 最后更改: 2024/10/22 03:03 由 yuxi

页面工具

除额外注明的地方外，本维基上的内容按下列许可协议发布： CC Attribution-Share Alike 4.0 International
CC Attribution-Share Alike 4.0 International Donate Powered by PHP Valid HTML5 Valid CSS Driven by DokuWiki