Skip to content

Latest commit

 

History

History
375 lines (253 loc) · 22.6 KB

File metadata and controls

375 lines (253 loc) · 22.6 KB

NetST命令行文档

从序列数据中生成单倍型网络

NetST可以基于CSV文件或FASTA文件中的序列生成单倍型网络,并输出一个可视化查看器。文件中的所有序列会被进行多序列比对,提取出信息位点,并根据分类信息排序。排序后的序列输入到FastHaN生成一个单倍型网络。输出的单倍型网络数据会再经过上色等处理,最后导出到一个可视化查看器。

命令

python netst.py viewhaplo

或  python ViewHaplo.py

参数说明

  • -input <str>:**必要。**用于生成单倍型网络的序列文件。可以是CSV或FASTA文件。
  • -net_type <str>:**必要。**指定生成单倍型网络的算法。可选的参数值和FastHaN一致,包括original_tcsmodified_tcsmsnmjn四个选项。
    • original_tcs是原始的Templeton–Crandall–Sing算法。
    • modified_tcs是PopART中实现的改进TCS算法。
    • msn是最小生成网络(Minimum Spanning Network)算法。
    • mjn是中间连接网络(Median-Joining Network)算法。
  • -output <str>:**必要。**输出单倍型网络查看器的文件夹。
  • -file_type <str>:**可选。**输入文件的文件类型,支持csvfas两个选项。默认值为fas
  • -index_name <int>:**可选。**序列名所在的列序号(从0开始计数)。
  • -index_seq <int>:**可选。**序列数据所在的列序号(从0开始计数)。对于FASTA文件,这个参数会被忽略。
  • -index_count <int>:**可选。**数量数据所在的列序号(从0开始计数)。这一数据会表现为网络中节点的大小。
  • -index_type <int>:**可选。**间断性状数据所在的列序号(从0开始计数)。
  • -index_quant <int>:**可选。**连续性状数据所在的列序号(从0开始计数)。这一数据会表现为网络中节点颜色的深浅。
  • -aligned:**可选。**假设输入文件已进行多序列比对,不再尝试调用muscle进行比对。

输入

输入序列支持CSV或FASTA格式。对于CSV格式,需要指定-index系列参数确定性状数据的位置。如这一段CSV文档需要的参数是-index_name 1 -index_seq 2 -index_type 3

ID, Name, Sequence, Type
1,  Sp1,  ATCGCC,   Hap1
2,  Sp2,  CTCGCC,   Hap2

对于FASTA格式,同样可以使用index系列参数确定性状数据在序列描述行中用|分割的位置。假设下列FASTA序列中数字代表种群数量,那么需要的参数是-index_name 0 -index_type 1 -index_count 2

>Sp1|Hap1|600
ACCACC...
>Sp2|Hap3|240
ACGACG...

如FASTA文件未指定-index系列参数,序列描述行允许的格式为序列名=数量=连续性状|单倍型,除序列名以外的元素均可省略,示例如下:

>seq_1
ATGCAT...

>seq_1|Haplotype1
ATGCAT...

>seq_1=30=60
ATGCAT...

>seq_1=30=60|Haplotype1
ATGCAT...

输出

输出文件均保存在-output指定的目录中。输出文件的说明如下:

  • viewer.html:单倍型网络查看器。这一查看器需要cssjs目录和network-config.js才能正确运行。
  • hap_aln.fasta:进行多序列比对后的输入序列。
  • network_seq.fasta:提取出的信息序列。
  • fasthan.gml:fastHaN生成的.gml文件,可用于导入其他单倍型网络可视化软件。

对序列分型并在单倍型网络中展示

NetST的命令行版本有一个独特的特性:单倍型网络分析管线。与直接生成单倍型网络相比,这一管线将序列分型和单倍型网络进行了整合,提高了分析的速度。这一管线需要三项输入:未确定分类的待分析FASTA序列、单倍型参考数据库、需要在单倍型网络中显示的已知序列。实现分为两步:第一步基于参考数据库对未知序列的单倍型进行分类;第二步利用推断出的单倍型信息,将未知序列放置在包含已知序列的单倍型网络中。

单倍型分类支持两种方法,分别基于BLAST和k-mer匹配。BLAST算法使用blastn工具的dc-megablast任务对输入序列进行比对,比对结果中的每一个相同碱基均计算为命中。K-mer算法对输入序列的所有k-mer进行匹配,计算k-mer在参考数据库中命中的比例。对某一输入序列,认为命中数量最多的单倍型代表较为可能的单倍型。完成单倍型分类后,有一个可选步骤搜索与每个输入序列有最多共同k-mer的参考序列。这些输出序列可以认为是与输入序列较为接近的参考序列,可用于对混合单倍型的系统发育进行更深入的研究。

第二步利用第一步的输出来生成单倍型网络。如果给定了单倍型网络中的已知部分,第一步输出的单倍型会先被加入到这组序列中。可视化的单倍型网络有利于对输入序列的进化事件进行全面的分析,并有助于检查管线第一步中自动分类结果的正确性。

命令

python netst.py pipeline

参数说明

  • -work_dir <str>:**必要。**保存中间文件和分析结果的文件夹。
  • -input <str>:**必要。**输入的未分类序列文件。必须是FASTA文件,可以含有多条序列。
    • 如果没有额外信息,序列描述可以直接使用样本的代号。例如:

      >seq_1
      ATGCAT...
      
    • 如果需要进行分型的序列含有连续型性状信息,可以按序列名=数量=连续性状格式在序列描述中编码,其中“数量”会表现为网络中节点的大小,而“性状”会表现为网络中节点颜色的深浅。

      >seq_1=30=60
      ATGCAT...
      
    • 注意输入序列名中不能有|符号,这一符号在软件内部用于编码性状信息。

  • -ref <str>:**必要。**指定进行单倍型分类的参考序列。必须是FASTA文件,且序列描述必须遵守序列名|单倍型格式。
  • -method <str>:**可选。**指定对未知序列进行单倍型分类的算法。目前支持blastkmer两个值,默认值为kmer
  • -net_type <str>:**可选。**指定生成单倍型网络的算法。可选的参数值和FastHaN一致,包括original_tcsmodified_tcsmsnmjn四个选项,默认值为modified_tcs
    • original_tcs是原始的Templeton–Crandall–Sing算法。
    • modified_tcs是PopART中实现的改进TCS算法。
    • msn是最小生成网络(Minimum Spanning Network)算法。
    • mjn是中间连接网络(Median-Joining Network)算法。
  • -haplo_seq <str>:**可选。**单倍型网络已知部分的性状和序列数据,可以使用-haplo_seq_type指定文件类型。如未指定-index系列参数,FASTA文件中序列描述行允许的格式为序列名=数量=连续性状|单倍型,除序列名以外的元素均可省略。
  • -haplo_seq_type <str>可选。-haplo_seq参数的文件类型,支持csvfas两个选项。
  • -index_name <int>可选。-haplo_seq参数对应的文件中,序列名所在的列序号(从0开始计数)。
  • -index_seq <int>可选。-haplo_seq参数对应的文件中,序列数据所在的列序号(从0开始计数)。对于FASTA文件,这个参数会被忽略。
  • -index_count <int>可选。-haplo_seq参数对应的文件中,数量数据所在的列序号(从0开始计数)。这一数据会表现为网络中节点的大小。
  • -index_type <int>可选。-haplo_seq参数对应的文件中,间断性状数据所在的列序号(从0开始计数)。
  • -index_quant <int>可选。-haplo_seq参数对应的文件中,连续性状数据所在的列序号(从0开始计数)。这一数据会表现为网络中节点颜色的深浅。
  • -db <str>:**可选。**指定预先生成的序列查找数据库。如果没有指定该参数或数据库文件不存在,程序会自动根据-ref参数生成数据库。注意k-mer数据库文件在不同k值间不能混用。
  • -cut <int>:**可选。**认为一个序列没有混合单倍型的最小命中百分比。如果命中百分比大于这个数字,则认为序列来自单一单倍型。默认值为85。
  • -k <int>:**可选。**使用-method kmer进行分类时的k-mer大小。默认值为21。
  • -e <int>:**可选。**使用-method blast进行分类时的最大E值。默认值为1e-5。
  • -word_size <int>:**可选。**使用-method blast进行分类时的词长。默认值为11。
  • -all_seqs:**可选。**默认情况下生成的网络只包含可能具有混合单倍型的序列(支持度小于-cut参数的序列)。如果添加了-all_seqs,所有待分析序列都会显示在网络中。
  • -skip_getmixseq:**可选。**跳过对具有最多相同k-mer的参考序列的搜索,提高运行速度。启用后输出文件夹中将不会再生成mix子文件夹。
  • -t <int>:**可选。**并行执行的最大进程数。

输出

输出文件均保存在-work_dir指定的目录中。输出文件的说明如下:

  • type_analysis.csv:所有输入序列的单倍型检测结果。每个可能的单倍型一行。
  • type_analysis_mix.csv:可能具有混合单倍型的输入序列。每个序列一行,可能的单倍型之间用|分割。
  • type_analysis_null.csv:没有找到匹配单倍型的输入序列。每个序列一行。
  • mix:与type_analysis_mix.csv中每条序列有最多相同k-mer的参考序列。每个输入序列一个文件。
  • network:单倍型网络数据和可视化查看器,直接使用浏览器打开network/viewer.html即可交互式浏览单倍型网络。

分析工具

除分析管线外,NetST命令行版本还包含与图形界面相同的分析工具。这些分析工具均可以通过命令行独立调用。

拆分FASTA文件

这一工具可以将包含多种单倍型的CSV文件或FASTA文件根据单倍型拆分成多个文件。

命令

python netst.py makedata

或  python MakeData.py

参数说明

  • -input <str>:**必要。**输入的序列文件。可以是CSV或FASTA文件。
  • -out_dir <str>:**必要。**输出文件夹。
  • -file_type <str>:**可选。**输入文件的文件类型,支持csvfas两个选项。默认值为fas
  • -index_name <int>:**可选。**序列名所在的列序号(从0开始计数)。
  • -index_type <int>:**可选。**单倍型名称所在的列序号(从0开始计数)。
  • -index_seq <int>:**可选。**CSV文件中序列数据所在的列序号(从0开始计数)。对于FASTA文件,这个参数会被忽略。
  • -clean [bool]:**可选。**是否清理文件中具有重复部分的序列。这个选项会清理较短的、被其他序列完全包含的序列,只保留一系列重叠序列中最长的一条。
  • -valid [bool]:**可选。**是否清理文件中包含无效碱基的序列。

将多个基因串联成超基因

这一工具可以将一个文件夹下按照基因组织的FASTA序列文件转化为按照物种组织的超基因FASTA序列。这一工具依赖Biopython运行。

假设gene_a.fas包含以下序列:

>Escherichia coli
ATGCTG
>Lactobacillus casei
CTGCTG

且gene_b.fas包含以下序列:

>Escherichia coli
AAGGAA
>Lactobacillus casei
ACGGAG

则本工具会生成如下文件:

>Escherichia coli
ATGCTGAAGGAA
>Lactobacillus casei
CTGCTGACGGAG

命令

python merge_seq.py

参数说明

  • -input <str>:**必要。**输入文件夹。程序将会在这个文件夹下搜索需要进行串联的序列。
  • -output <str>:**必要。**输出的FASTA文件。
  • -exts <str>:**可选。**逗号分隔的文件扩展名列表,例如.fasta,.seq。程序将会搜索具有列表中任一扩展名的文件。注意这一列表中不能含有空格。默认值为.fasta,.fas,.fa
  • -missing <str>:**可选。**用于填充序列空白的字符,适用于某一基因在部分物种中不存在的情况。

拆分GenBank文件

这一工具可以对GenBank文件中的每一条记录进行拆分,并将所有序列的元数据保存到CSV文件中。这一工具依赖Biopython运行。

命令

python build_gb.py

参数说明

  • -input <str>:**必要。**输入GenBank文件。
  • -output <str>:**必要。**输出拆分的GenBank记录的文件夹。
  • -outcsv <str>:**必要。**输出的元数据CSV文件。
  • -clean [bool]:**可选。**是否删除输出文件夹中已有的文件。错误使用可能导致数据丢失,请务必小心!

使用BLAST进行单倍型分类

这一工具使用blastn工具的dc-megablast任务对输入序列进行比对,并输出可能的单倍型。输入输出格式与单倍型网络分析管线相同。

命令

python netst.py getblast

或  python GetBlast.py

参数说明

  • -input <str>:**必要。**输入的未分类序列文件。必须是FASTA文件,可以含有多条序列。
  • -ref <str>:**必要。**指定进行单倍型分类的参考序列。必须是FASTA文件,且序列描述必须遵守序列名|单倍型格式。
  • -db <str>:**必要。**指定BLAST数据库的位置。如果数据库文件不存在,程序会自动调用makeblastdb基于参考序列生成。
  • -output <str>:**必要。**输出文件的名称。
  • -blast_result <str>:**可选。**这个参数请求程序解析已有的BLAST输出结果文件(要求用-outfmt 6生成),而不是调用blastn运行查询。
  • -e <int>:**可选。**最大E值。默认值为1e-5。
  • -word_size <int>:**可选。**词长。默认值为11。
  • -cut <int>:**可选。**认为一个序列没有混合单倍型的最小命中百分比。如果命中百分比大于这个数字,则认为序列来自单一单倍型。默认值为85。
  • -t <int>:**可选。**运行BLAST查询的最大线程数。

输出

-output参数指定输出文件名。假设使用-output get_blast,输出文件的说明如下:

  • get_blast.csv:所有输入序列的单倍型检测结果。每个可能的单倍型一行。
  • get_blast_mix.csv:可能具有混合单倍型的输入序列。每个序列一行,可能的单倍型之间用|分割。
  • get_blast_null.csv:没有找到匹配单倍型的输入序列。每个序列一行。
  • get_blast_blast.tbl:使用-outfmt 6产生的blastn搜索结果。

使用k-mer进行单倍型分类

这一工具在数据库中查找输入序列和单倍型的k-mer匹配,并输出输入序列较为可能的单倍型。输入输出格式与单倍型网络分析管线相同。

命令

python netst.py gettype

或  python GetType.py

参数说明

  • -input <str>:**必要。**输入的未分类序列文件。必须是FASTA文件,可以含有多条序列。
  • -ref <str>:**必要。**指定进行单倍型分类的参考序列。必须是FASTA文件,且序列描述必须遵守序列名|单倍型格式。
  • -db <str>:**必要。**指定k-mer数据库的位置。如果数据库文件不存在,程序会自动根据参考序列在指定位置生成。注意k-mer数据库文件在不同k值间不能混用。
  • -output <str>:**必要。**输出文件的名称。
  • -k <int>:**可选。**K-mer的长度。默认值为21。请注意不同的k-mer长度对应不同的数据库,不能混用。
  • -cut <int>:**可选。**认为一个序列没有混合单倍型的最小命中百分比。如果命中百分比大于这个数字,则认为序列来自单一单倍型。默认值为85。

输出

-output参数指定输出文件名。假设使用-output get_type,输出文件的说明如下:

  • get_type.csv:所有输入序列的单倍型检测结果。每个可能的单倍型一行。
  • get_type_mix.csv:可能具有混合单倍型的输入序列。每个序列一行,可能的单倍型之间用|分割。
  • get_type_null.csv:没有找到匹配单倍型的输入序列。每个序列一行。

搜索具有k-mer命中的参考序列

这一工具查找输入序列与参考序列的共同k-mer,并输出每个输入序列找到最多共同k-mer的一组相似序列。

命令

python netst.py getmixseq

或  python GetMixSeq.py

参数说明

  • -input <str>:**必要。**输入的CSV文件,需要是BLAST或k-mer单倍型分类输出的_mix.csv文件。
  • -out_dir <str>:**必要。**输出文件夹。
  • -ref_dir <str>:**二选一。**根据单倍型进行拆分产生的参考文件目录。
  • -ref_file <str>:**二选一。**参考序列文件,与单倍型分类的参考序列文件相同。
  • -k <int>:**可选。**K-mer的长度。默认值为21。
  • -max <int>:**可选。**对每条输入序列输出的最大相似序列数量。
  • -t <int>:**可选。**并行执行的最大进程数。

设计引物

这一工具可以针对基因的一致区域设计引物,有助于进行分子材料的分类。这一工具依赖Biopython运行。

命令

python build_primer.py

参数说明

  • -input <str>:**必要。**输入的文件或文件夹。目前支持的输入包括单个.gb文件、文件夹内的所有.gb文件,以及已进行多序列比对的单个FASTA格式文件。
  • -out_dir <str>:**必要。**输出文件夹。
  • -soft_boundary <int>:**可选。**软边界的长度,即基因两端纳入备选区段的长度。默认为200。
  • -max_marker_length <int>:**可选。**允许的最长marker长度。默认为2000。
  • -min_marker_length <int>:**可选。**允许的最短marker长度。默认为200。
  • -max_seq_length <int>:**可选。**允许的最长基因长度。默认为5000。
  • -min_seq_length <int>:**可选。**允许的最短基因长度。默认为200。
  • -gap_length <int>:**可选。**左右引物间的最小间隔长度。默认为100。
  • -max_primer <int>:**可选。**引物的最大数量。默认为64。
  • -usevar [bool]:**可选。**是否只使用一致位点设计引物。默认启用。
  • -split_only [bool]:**可选。**只对序列进行拆分和比对,不设计引物。
  • -t <int>:**可选。**并行执行的最大进程数。

使用示例

通过已知性状和序列构建单倍型网络

本节使用DEMOS/DEMO_01下的文件演示使用viewhaplo子命令构建单倍型网络。

对于DEMOS/DEMO_01/HA_REFS.csv,需要指定-file_type csv。从0开始计数,表头中ID位于第0位,序列名称、序列数据、数量、间断性状、连续性状分别位于第1位、第2位、第5位、第3位、第4位,所以需要指定-index_name 1 -index_seq 2 -index_count 5 -index_type 3 -index_quant 4。生成MJN网络的命令如下:

python netst.py viewhaplo -input DEMOS/DEMO_01/HA_REFS.csv -file_type csv -net_type mjn -index_name 1 -index_seq 2 -index_count 5 -index_type 3 -index_quant 4 -output output/demo_01_mjn

对于DEMOS/DEMO_01/TEST_HA.fasta,指定的文件类型参数是-file_type fas。从0开始计数,以>开头的描述行中序列名称位于|分隔的第0位,数量位于第3位,间断性状位于第1位,连续性状位于第2位,所以指定-index_name 0 -index_count 3 -index_type 1 -index_quant 2。生成PopART风格的TCS网络的命令如下:

python netst.py viewhaplo -input DEMOS/DEMO_01/TEST_HA.fasta -file_type fas -net_type modified_tcs -index_name 0 -index_count 3 -index_type 1 -index_quant 2 -output output/demo_01_tcs

使用浏览器打开output/demo_01_mjn/viewer.htmloutput/demo_01_tcs/viewer.html即可查看单倍型网络的可视化。

序列分型

本节使用DEMOS/DEMO_02下的文件演示对序列进行分型分析的过程。

生成分型序列文件

本节从DEMOS/DEMO_02/HA_REFS.fasta文件,通过makedata子命令生成分型序列文件,输出到HA_DATA文件夹中。描述行中序列名位于第0位,单倍型名称位于第1位,指定-index_name 0 -index_type 1。输出到HA_DATA文件夹的命令如下:

python netst.py makedata -input DEMOS/DEMO_02/HA_REFS.fasta -file_type fas -index_name 0 -index_type 1 -out_dir HA_DATA

完成后可在HA_DATA文件夹下看到按单倍型名称分类的序列文件。这一文件夹可作为参考序列多次使用。

对未知序列进行分类

本节以DEMOS/DEMO_02/HA_TEST.fasta作为未知序列,演示使用gettype子命令进行分型。使用的参考序列是上一小节生成的HA_TEST/ref_combine.fasta。注意,命令行版本需要显式指定一个k-mer数据库的位置,但不需要文件存在。如果文件不存在,会自动在指定的位置生成一个gzip格式的数据库。这一数据库可以多次使用以加速后续分类的运行速度,但指定的k值必须与数据库配套。命令行版本输出所有可能的混合分型,如果只需要单一分型,可以对输出进行简单的处理,例如在Linux环境下运行sort -t ',' -n -k 1 -u type.csv即可。按照默认参数,在output/demo_02/db.gz生成数据库文件,并输出分析结果到output/demo_02/type的命令如下:

python netst.py gettype -input DEMOS/DEMO_02/HA_TEST.fasta -ref HA_DATA/ref_combine.fasta -db output/demo_02/db.gz -output output/demo_02/type

在output/demo_02文件夹下有三个输出文件:

  • type.csv:所有序列的所有可能的单倍型列表。
  • type_mix.csv:可能的混合序列列表。
  • type_null.csv:没有找到单倍型匹配的序列列表。

如果要使用BLAST生成类似的分型结果。可以使用getblast子命令。输出的blast_type开头的三个csv文件与gettype输出的csv文件格式相同。但请注意,由于BLAST查询易受数据库的影响,输出的支持度缺乏明确的生物意义而不建议作为其他分析的依据。

python netst.py getblast -input DEMOS/DEMO_02/HA_TEST.fasta -ref HA_DATA/ref_combine.fasta -db output/demo_02/blast_db -output output/demo_02/blast_type

相似序列分析

上一小节的输出结果默认已经是混合分型。对于需要搜索相似序列的混合分型序列,可以使用getmixseq命令搜索相似序列。

python netst.py getmixseq -input output/demo_02/type_mix.csv -ref_dir HA_DATA -out_dir output/demo_02/mix

在output/demo_02/mix文件夹下输出了4号序列和9号序列的相似序列。如果要生成seq_04和相似序列的单倍型网络,可以再次使用viewhaplo子命令。此处我们用k-mer命中数量作为节点大小绘图:

python netst.py viewhaplo -input output/demo_02/mix/4.fasta -file_type fas -net_type modified_tcs -index_name 0 -index_count 2 -index_type 1 -index_quant 2 -output output/demo_02/seq_04_network

用浏览器打开output/demo_02/seq_04_network/viewer.html即可查看单倍型网络。