====== 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.下载需要进行构建数据库的序列内容，可以为核酸序列，可以为蛋白序列 {{:yuxi:viral.png?500|}} 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文件 {{:yuxi:accession2taxid.png?400|}} 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.构建注释索引库 {{:yuxi:taxonomy.png?400|}} 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 ---- * 官方网站：[[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蛋白水平组装，以增加复杂元基因组的蛋白质序列回收。