售后

编辑本页后请点击“保存”。请参阅 syntax 了解维基语法。只有在您能改进该页面的前提下才编辑它。如果您想尝试一些东西，请先到 playground 热身。

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

===== 1.软件介绍 =====
MMseqs2（Many-against-Many序列搜索）是一个软件套件，用于搜索和集群巨大的蛋白质和核苷酸序列集。MMseqs2是开源GPL许可软件，在Linux、MacOS和（作为测试版，通过cygwin）Windows的C++中实现。该软件旨在在多个内核和服务器上运行，并具有非常好的可扩展性。MMseqs2的运行速度比BLAST快10000倍。以其速度的100倍，它的灵敏度几乎相同。它可以以与PSI-BLAST相同的灵敏度以400多倍的速度执行配置文件搜索。
===== 2.聚类 =====
<code>
mmseqs easy-cluster examples/DB.fasta clusterRes tmp --min-seq-id 0.5 -c 0.8 --cov-mode 1
</code>
运行时随着输入大小而线性缩放。此模式推荐用于大型数据集。
<code>
mmseqs easy-linclust examples/DB.fasta clusterRes tmp
</code>
===== 3.比对 =====
<code>
mmseqs easy-search examples/QUERY.fasta examples/DB.fasta alnRes.m8 tmp
</code>
也可以预先计算目标数据库的索引。这减少了重复搜索同一数据库时的开销。
<code>
mmseqs createdb examples/DB.fasta targetDB
mmseqs createindex targetDB tmp
mmseqs easy-search examples/QUERY.fasta targetDB alnRes.m8 tmp
</code>
databases工作流程为许多公共参考数据库提供了下载和设置程序，如Uniref、NR、NT、PFAM等（请参阅下载数据库 https://github.com/soedinglab/mmseqs2/wiki#downloading-databases ）。例如，要下载并搜索包含Swiss-Prot参考蛋白的数据库，请运行：
<code>
mmseqs databases UniProtKB/Swiss-Prot swissprot tmp
mmseqs easy-search examples/QUERY.fasta swissprot alnRes.m8 tmp
</code>
===== 4.物种注释 =====
easy-taxonomy工作流程可用于分配序列分类标签。它对具有分类信息（seqTaxDb）的序列数据库进行搜索，根据不同的策略（根据--lca-mode）选择最具代表性的对齐目标序列集，并计算其中最低的共同祖先。
<code>
mmseqs createdb examples/DB.fasta targetDB
mmseqs createtaxdb targetDB tmp
mmseqs createindex targetDB tmp
mmseqs easy-taxonomy examples/QUERY.fasta targetDB alnRes tmp
</code>

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

----
==== 一、构建数据库 ====
1.激活环境
<code>
source /TJPROJ1/META_ASS/soft/mmseqs2/../anaconda3/bin/activate /TJPROJ1/META_ASS/soft/mmseqs2
</code>
2.下载需要进行构建数据库的序列内容，可以为核酸序列，可以为蛋白序列

{{:yuxi:viral.png?500|}}
<code>
wget -c https://ftp.ncbi.nlm.nih.gov/refseq/release/viral/viral.1.1.genomic.fna.gz
</code>
3.构建序列索引数据库
<code>
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
</code>
4.下载注释信息文件,Refseq数据库中可以从nt库中进行查询，因此需要nt、nr库中相关注释names.dmp或nodes.dmp文件与accession2taxid文件
{{:yuxi:accession2taxid.png?400|}}
<code>
wget https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/accession2taxid/nucl_gb.accession2taxid.gz
</code>
5.构建mapping文件，该步骤的目的是为了获得Refseq数据库中的序列和taxid的对应关系
<code>
python3 getname.py
python3 mapping.py
</code>
其中getname.py的脚本内容是：
<code>
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)
</code>
mapping.py的脚本内容是：
<code>
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)
</code>
6.构建注释索引库

{{:yuxi:taxonomy.png?400|}}
<code>
mkdir ncbi-taxdump && cd ncbi-taxdump
wget https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxdump.tar.gz
tar xzvf taxdump.tar.gz
cd -
</code>
7.构建mmseqs专用的数据库格式，并输入指定名称
<code>
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
</code>

----
==== 二、物种注释 ====
1.激活环境
<code>
source /TJPROJ1/META_ASS/soft/mmseqs2/../anaconda3/bin/activate /TJPROJ1/META_ASS/soft/mmseqs2
</code>
2.将需要被比对的、注释的序列转换为数据库格式，如果不转化，则下面的search和taxonomy提示数据格式不可用，也可以使用easy-模式
<code>
mmseqs createdb CK.1.1.fasta queryDB
</code>
3.序列比对
<code>
mmseqs search queryDB /TJPROJ1/META_ASS/soft/mmseqs2/database/test/refseq_viral.1.1 resultDB tmp
</code>
4.物种注释
<code>
mmseqs taxonomy queryDB /TJPROJ1/META_ASS/soft/mmseqs2/database/test/refseq_viral.1.1 taxonomyResult tmp
</code>
5.导出注释结果
<code>
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
</code>

----

===== 拓展：如何将Silva数据库转换为NCBI node/names 的格式 =====
==== 拓展1.普通列表项目为SILVA创建一个seqTaxDB ====
如果想根据SILVA数据库对核糖体RNA（16S、18S、SSU）序列进行分类，您可以通过databases工作流程进行下载：
<code>
mmseqs databases SILVA silvadb tmp
</code>
要了解SILVA MMseqs2分类序列数据库是如何构建的，请查看下面的示例脚本。我们首先需要从SILVA分类法中创建一个像分类法一样的NCBI。
<code>
# 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
</code>
==== 拓展2.为GTDB创建一个seqTaxDB ====
基因组分类数据库（GTDB）是一个系统发育一致的数据库，它重新定义了分类树。MMseqs2可以搜索GTDB，但它需要一些预处理步骤。
<code>
# 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
</code>
==== 拓展3.通过序列数据库的手动注释创建seqTaxDB ====
这是一个如何用分类信息手动注释序列数据库的示例。该示例使用Uniprot标识符。
作为第一步，将FAST[A/Q]文件转换为mmseqs序列数据库，使用createdb
<code>
# 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
</code>
createdb生成一个制表符分隔的sequenceDB.lookup文件，其中包含数字-db-id，Accession（例如Uniprot Accession Q6GZX4）和文件。ID是从输入数据库的标头中解析的（请参阅从标头解析的id）。
<code>
0 Q6GZX4  0
1 Q6GZX3  0
2 Q197F8  0
3 P0A031  0
4 Q197F7  0
</code>
作为下一步，我们创建一个制表符分隔的映射，每个目标数据库标识符都映射到NCBI分类群标识符。映射文件应采用Accession numeric-ncbi-tax-id的格式。
<code>
Q6GZX4    654924
Q6GZX3    654924
Q197F8    345201
Q197F7    345201
</code>
以下是如何将Uniprot映射文件转换为制表符分隔的映射文件的示例。其他数据库可以参考我上面写的Refseq数据库构建
<code>
# 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
</code>
我们需要NCBI分类法taxdump.tar.gz。它在NCBI FTP服务器上可用：
<code>
mkdir ncbi-taxdump && cd ncbi-taxdump
wget ftp://ftp.ncbi.nih.gov/pub/taxonomy/taxdump.tar.gz
tar xzvf taxdump.tar.gz
cd -
</code>
作为最后一步，我们现在可以使用createtaxdb来注释我们的序列数据库。
<code>
# now we can use createtaxdb with our own mapping.
mmseqs createtaxdb sequenceDB tmp --ncbi-tax-dump ncbi-taxdump --tax-mapping-file taxidmapping
</code>
可以检查有多少标识符具有以下代码的指定分类群
<code>
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
</code>
----

* 官方网站：[[https://github.com/soedinglab/MMseqs2]]
  * 用户使用说明：[[https://github.com/soedinglab/mmseqs2/wiki]]
  * 实战教程：[[https://github.com/soedinglab/MMseqs2/wiki/Tutorials]]
  * 软件路径：/TJPROJ1/META_ASS/soft/mmseqs2/
  * 数据库构建路径：/TJPROJ1/META_ASS/soft/mmseqs2/database/refseq_viral.1.1
  * 注释测试路径：/TJPROJ1/META_ASS/soft/mmseqs2/test

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

编辑摘要

当您选择开始编辑本页，即寓示你同意将你贡献的内容按下列许可协议发布： CC Attribution-Share Alike 4.0 International

售后

用户工具

站点工具

页面工具