二代测序技术的迅速发展大幅度的降低了高通量测序价格且提高了测序质量，同时高通量测序相关的数据处理和分析的软件也是越来越好用。SPAdes在二代测序小基因组的组装上表现优异，使用方便，真是爱不释手。Unicycler在二代测序小基因组装上则是依赖于SPAdes，在此基础上进行了优化；此外在二代三代测序数据混合组装和纯三代测序数据组装上也表现较好；该项目在github上被标星210次。今天使用SPAdes和Unicycler对一个350bp小片段建库测序1G数据量的细菌菌株进行组装，熟悉一下流程并且看看两个软件的差别。

1. 数据质控

公司给的建库测序的数据分为raw_data和clean_data, 我先直接使用fastqc对clean_data的数据进行质控。双端测序数据是N46-1_add_BDMS190627350-1a_1.clean.fq.gz和N46-1_add_BDMS190627350-1a_2.clean.fq.gz。

(denovo_genome) czh@ubuntu:~/Desktop/01_denovo_genome/N46_DENOVO$ fastqc '/home/czh/Desktop/01_denovo_genome/N46_DENOVO/N46-1_add_BDMS190627350-1a_1.clean.fq.gz' '/home/czh/Desktop/01_denovo_genome/N46_DENOVO/N46-1_add_BDMS190627350-1a_2.clean.fq.gz'  -o '/home/czh/Desktop/01_denovo_genome/N46_DENOVO/fastqc_out'

正向序列

反向序列

fastqc结果看，测序长度是150bp，序列数量约为8M，数据量约为1.2G，序列GC含量和duplication水平不是特别完美，总体来看双端序列质量很好，数据的质量和数量满足下游分析。

2. 使用spades对Illumina双端序列进行组装。

(denovo_genome) czh@ubuntu:~/Desktop/01_denovo_genome/N46_DENOVO$ spades.py -1 N46-1_add_BDMS190627350-1a_1.clean.fq.gz -2 N46-1_add_BDMS190627350-1a_2.clean.fq.gz --careful -o spades_out -t 12 -m 24
Command line: /home/czh/miniconda3/envs/denovo_genome/bin/spades.py -1  /home/czh/Desktop/01_denovo_genome/N46_DENOVO/N46-1_add_BDMS190627350-1a_1.clean.fq.gz  -2  /home/czh/Desktop/01_denovo_genome/N46_DENOVO/N46-1_add_BDMS190627350-1a_2.clean.fq.gz  --careful   -o  /home/czh/Desktop/01_denovo_genome/N46_DENOVO/spades_out    -t  12  -m  24  

System information:
  SPAdes version: 3.13.1
  Python version: 3.7.3
  OS: Linux-5.0.0-23-generic-x86_64-with-debian-buster-sid

Output dir: /home/czh/Desktop/01_denovo_genome/N46_DENOVO/spades_out
Mode: read error correction and assembling
Debug mode is turned OFF

Dataset parameters:
  Multi-cell mode (you should set '--sc' flag if input data was obtained with MDA (single-cell) technology or --meta flag if processing metagenomic dataset)
  Reads:
    Library number: 1, library type: paired-end
      orientation: fr
      left reads: ['/home/czh/Desktop/01_denovo_genome/N46_DENOVO/N46-1_add_BDMS190627350-1a_1.clean.fq.gz']
      right reads: ['/home/czh/Desktop/01_denovo_genome/N46_DENOVO/N46-1_add_BDMS190627350-1a_2.clean.fq.gz']
      interlaced reads: not specified
      single reads: not specified
      merged reads: not specified
Read error correction parameters:
  Iterations: 1
  PHRED offset will be auto-detected
  Corrected reads will be compressed
Assembly parameters:
  k: automatic selection based on read length
  Repeat resolution is enabled
  Mismatch careful mode is turned ON
  MismatchCorrector will be used
  Coverage cutoff is turned OFF
Other parameters:
  Dir for temp files: /home/czh/Desktop/01_denovo_genome/N46_DENOVO/spades_out/tmp
  Threads: 12
  Memory limit (in Gb): 24


======= SPAdes pipeline started. Log can be found here: /home/czh/Desktop/01_denovo_genome/N46_DENOVO/spades_out/spades.log


===== Read error correction started. 


== Running read error correction tool: /home/czh/miniconda3/envs/denovo_genome/bin/spades-hammer /home/czh/Desktop/01_denovo_genome/N46_DENOVO/spades_out/corrected/configs/config.info

  0:00:00.000     4M / 4M    INFO    General                 (main.cpp                  :  75)   Starting BayesHammer, built from N/A, git revision N/A
  0:00:00.000     4M / 4M    INFO    General                 (main.cpp                  :  76)   Loading config from /home/czh/Desktop/01_denovo_genome/N46_DENOVO/spades_out/corrected/configs/config.info
  0:00:00.000     4M / 4M    INFO    General                 (main.cpp                  :  78)   Maximum # of threads to use (adjusted due to OMP capabilities): 12
  0:00:00.000     4M / 4M    INFO    General                 (memory_limit.cpp          :  49)   Memory limit set to 24 Gb
  0:00:00.000     4M / 4M    INFO    General                 (main.cpp                  :  86)   Trying to determine PHRED offset
  0:00:00.000     4M / 4M    INFO    General                 (main.cpp                  :  92)   Determined value is 33

耗时大约70分钟。

3. 使用unicycler对Illumina双端序列进行组装。

(denovo_genome) czh@ubuntu:~/Desktop/01_denovo_genome/N46_DENOVO$ unicycler -1 N46-1_add_BDMS190627350-1a_1.clean.fq.gz -2 N46-1_add_BDMS190627350-1a_2.clean.fq.gz -o unicycler_out -t 12
Starting Unicycler (2019-11-03 22:46:05)
    Welcome to Unicycler, an assembly pipeline for bacterial genomes. Since you
provided only short reads, Unicycler will essentially function as a SPAdes-
optimiser. It will try many k-mer sizes, choose the best based on contig length
and graph connectivity, and scaffold the graph using SPAdes repeat resolution.
    For more information, please see https://github.com/rrwick/Unicycler

Command: /home/czh/miniconda3/envs/denovo_genome/bin/unicycler -1 /home/czh/Desktop/01_denovo_genome/N46_DENOVO/N46-1_add_BDMS190627350-1a_1.clean.fq.gz -2 /home/czh/Desktop/01_denovo_genome/N46_DENOVO/N46-1_add_BDMS190627350-1a_2.clean.fq.gz -o /home/czh/Desktop/01_denovo_genome/N46_DENOVO/add_out -t 12

Unicycler version: v0.4.8
Using 12 threads

The output directory already exists:
  /home/czh/Desktop/01_denovo_genome/N46_DENOVO/add_out

Dependencies:
  Program         Version           Status  
  spades.py       3.13.1            good    
  racon                             not used
  makeblastdb     2.9.0+            good    
  tblastn         2.9.0+            good    
  bowtie2-build   2.3.5             good    
  bowtie2         2.3.5             good    
  samtools        1.9               good    
  java            11.0.1-internal   good    
  pilon           1.23              good    
  bcftools                          not used


SPAdes read error correction (2019-11-03 22:47:16)
    Unicycler uses the SPAdes read error correction module to reduce the number
of errors in the short read before SPAdes assembly. This can make the assembly
faster and simplify the assembly graph structure.

Command: /home/czh/miniconda3/envs/denovo_genome/bin/spades.py -1 /home/czh/Desktop/01_denovo_genome/N46_DENOVO/N46-1_add_BDMS190627350-1a_1.clean.fq.gz -2 /home/czh/Desktop/01_denovo_genome/N46_DENOVO/N46-1_add_BDMS190627350-1a_2.clean.fq.gz -o /home/czh/Desktop/01_denovo_genome/N46_DENOVO/add_out/spades_assembly/read_correction --threads 12 --only-error-correction

Corrected reads:
  /home/czh/Desktop/01_denovo_genome/N46_DENOVO/add_out/spades_assembly/corrected_1.fastq.gz
  /home/czh/Desktop/01_denovo_genome/N46_DENOVO/add_out/spades_assembly/corrected_2.fastq.gz


Choosing k-mer range for assembly (2019-11-03 23:47:00)
    Unicycler chooses a k-mer range for SPAdes based on the length of the input
reads. It uses a wide range of many k-mer sizes to maximise the chance of
finding an ideal assembly.

SPAdes maximum k-mer: 127
Median read length: 150
K-mer range: 27, 47, 63, 77, 89, 99, 107, 115, 121, 127


SPAdes assemblies (2019-11-03 23:48:07)
    Unicycler now uses SPAdes to assemble the short reads. It scores the
assembly graph for each k-mer using the number of contigs (fewer is better) and
the number of dead ends (fewer is better). The score function is 1/(c*(d+2)),
where c is the contig count and d is the dead end count.

K-mer   Contigs   Dead ends   Score   
   27                           failed
   47                           failed
   63                           failed
   77                           failed
   89                           failed
   99                           failed
  107                           failed
  115                           failed
  121                           failed
  127       198           0   2.53e-03 <-best

Read depth filter: removed 216 contigs totalling 72842 bp
Deleting /home/czh/Desktop/01_denovo_genome/N46_DENOVO/add_out/spades_assembly/


Determining graph multiplicity (2019-11-04 00:07:15)
    Multiplicity is the number of times a sequence occurs in the underlying
sequence. Single-copy contigs (those with a multiplicity of one, occurring only
once in the underlying sequence) are particularly useful.

Saving /home/czh/Desktop/01_denovo_genome/N46_DENOVO/add_out/001_best_spades_graph.gfa


Cleaning graph (2019-11-04 00:07:15)
    Unicycler now performs various cleaning procedures on the graph to remove
overlaps and simplify the graph structure. The end result is a graph ready for
bridging.

Graph overlaps removed

Removed zero-length segments:
    111, 112, 116, 117, 122, 123, 125, 126, 134, 135, 140, 143, 144, 145, 146,
149, 151, 153, 158, 161, 169, 170, 171, 172, 179, 181, 183

Removed zero-length segments:
    110, 165

Removed zero-length segments:
    195

Merged small segments:
    184, 185, 186, 187, 192, 196, 197

Saving /home/czh/Desktop/01_denovo_genome/N46_DENOVO/add_out/002_overlaps_removed.gfa

    Unicycler now selects a set of anchor contigs from the single-copy contigs.
These are the contigs which will be connected via bridges to form the final
assembly.

43 anchor segments (5,368,430 bp) out of 162 total segments (5,423,141 bp)


Creating SPAdes contig bridges (2019-11-04 00:07:15)
    SPAdes uses paired-end information to perform repeat resolution (RR) and
produce contigs from the assembly graph. SPAdes saves the graph paths
corresponding to these contigs in the contigs.paths file. When one of these
paths contains two or more anchor contigs, Unicycler can create a bridge from
the path.

                                                                                                                                         Bridge
Start                                                            Path                                                           End     quality
  -35                                                            -58                                                            44         18.2
  -19       109 -> 114 -> 95 -> -136 -> 106 -> 112 -> 123 -> 80 -> 130 -> 65 -> 136 -> 110 -> -145 -> -103 -> -109 -> 131       43         12.0
    2                                                            -42                                                            37          9.3
    3                                                            102                                                            10         62.2
    5                                                     -96 -> 152 -> -96                                                     -12        37.8
    8                                                            -102                                                           -15        61.8
   10                                                      -98 -> 46 -> -98                                                     1           7.3
   16   -150 -> 113 -> 133 -> -103 -> -109 -> 131 -> 45 -> 132 -> 109 -> 114 -> 124 -> 129 -> 83 -> 121 -> -117 -> -88 -> 156   -28         7.5
   22                                                            -57                                                            -40        11.8
   37                                                            -42                                                            4           8.7


Creating loop unrolling bridges (2019-11-04 00:07:15)
    When a SPAdes contig path connects an anchor contig with the middle contig
of a simple loop, Unicycler concludes that the sequences are contiguous (i.e.
the loop is not a separate piece of DNA). It then uses the read depth of the
middle and repeat contigs to guess the number of times to traverse the loop and
makes a bridge.

                                  Loop count   Loop count    Loop    Bridge
Start   Repeat   Middle     End    by repeat    by middle   count   quality
    2      -42       37       4         1.07         0.97       1      40.5
    5      -96      152     -12         0.67         0.60       1      36.4
    9       97      147      12         0.82         0.71       1      40.1
   10      -98       46       1         0.12         0.88       1      37.5


Applying bridges (2019-11-04 00:07:15)
    Unicycler now applies to the graph in decreasing order of quality. This
ensures that when multiple, contradictory bridges exist, the most supported
option is used.

Bridge type   Start -> end   Path                                       Quality
SPAdes           3 -> 10     102                                         62.210
SPAdes           8 -> -15    -102                                        61.764
SPAdes           5 -> -12    -96, 152, -96                               37.778
SPAdes         -35 -> 44     -58                                         18.167
SPAdes         -19 -> 43     109, 114, 95, -136, 106, 112, 123, 80,      11.994
                             130, 65, 136, 110, -145, -103, -109, 131          
SPAdes          22 -> -40    -57                                         11.829
loop             2 -> 4      -42, 37, -42                                40.460
loop             9 -> 12     97, 147, 97                                 40.120
loop            10 -> 1      -98, 46, -98                                37.498

Saving /home/czh/Desktop/01_denovo_genome/N46_DENOVO/add_out/003_bridges_applied.gfa


Bridged assembly graph (2019-11-04 00:07:16)
    The assembly is now mostly finished and no more structural changes will be
made. Ideally the assembly graph should now have one contig per replicon and no
erroneous contigs (i.e a complete assembly). If there are more contigs, then
the assembly is not complete.

Saving /home/czh/Desktop/01_denovo_genome/N46_DENOVO/add_out/004_final_clean.gfa

Component   Segments   Links   Length      N50         Longest segment   Status    
        1        131     184   5,430,795   1,298,586         1,645,591   incomplete


Polishing assembly with Pilon (2019-11-04 00:07:16)
    Unicycler now conducts multiple rounds of Pilon in an attempt to repair any
remaining small-scale errors with the assembly.

Aligning reads to find appropriate insert size range...
Insert size 1st percentile:  190
Insert size 99th percentile: 513

Pilon polish round 1
Total number of changes: 29

Pilon polish round 2
Total number of changes: 4

Pilon polish round 3
No Pilon changes

Saving /home/czh/Desktop/01_denovo_genome/N46_DENOVO/add_out/005_polished.gfa


Assembly complete (2019-11-04 02:13:11)
Saving /home/czh/Desktop/01_denovo_genome/N46_DENOVO/add_out/assembly.gfa
Saving /home/czh/Desktop/01_denovo_genome/N46_DENOVO/add_out/assembly.fasta

总共耗时3小时左右，比较慢。

4. 使用QUAST评估两个软件组装结果。

(qc_genome) czh@ubuntu:~/Desktop/01_denovo_genome/N46_DENOVO$ quast -o quast_out(qc_genome) czh@ubuntu:~/Desktop/01_denovo_genome/N46_DENOVO$ quast -o quast_out/ -t 8 spades_assembly.fasta unicycler_assembly.fasta 
/home/czh/miniconda3/envs/qc_genome/bin/quast -o quast_out/ -t 8 spades_assembly.fasta unicycler_assembly.fasta

Version: 5.0.2

System information:
  OS: Linux-5.0.0-23-generic-x86_64-with-debian-buster-sid (linux_64)
  Python version: 3.6.7
  CPUs number: 16

Started: 2019-11-04 06:14:57

Logging to /home/czh/Desktop/01_denovo_genome/N46_DENOVO/quast_out/quast.log

CWD: /home/czh/Desktop/01_denovo_genome/N46_DENOVO
Main parameters: 
  MODE: default, threads: 8, minimum contig length: 500, minimum alignment length: 65, \
  ambiguity: one, threshold for extensive misassembly size: 1000

Contigs:
  Pre-processing...
  1  spades_assembly.fasta ==> spades_assembly
  2  unicycler_assembly.fasta ==> unicycler_assembly

2019-11-04 06:15:06
Running Basic statistics processor...
  Contig files: 
    1  spades_assembly
    2  unicycler_assembly
  Calculating N50 and L50...
    1  spades_assembly, N50 = 472283, L50 = 5, Total length = 5427777, GC % = 57.74, # N's per 100 kbp =  15.92
    2  unicycler_assembly, N50 = 1298586, L50 = 2, Total length = 5419702, GC % = 57.76, # N's per 100 kbp =  0.00
  Drawing Nx plot...
    saved to /home/czh/Desktop/01_denovo_genome/N46_DENOVO/quast_out/basic_stats/Nx_plot.pdf
  Drawing cumulative plot...
    saved to /home/czh/Desktop/01_denovo_genome/N46_DENOVO/quast_out/basic_stats/cumulative_plot.pdf
  Drawing GC content plot...
    saved to /home/czh/Desktop/01_denovo_genome/N46_DENOVO/quast_out/basic_stats/GC_content_plot.pdf
  Drawing spades_assembly GC content plot...
    saved to /home/czh/Desktop/01_denovo_genome/N46_DENOVO/quast_out/basic_stats/spades_assembly_GC_content_plot.pdf
  Drawing unicycler_assembly GC content plot...
    saved to /home/czh/Desktop/01_denovo_genome/N46_DENOVO/quast_out/basic_stats/unicycler_assembly_GC_content_plot.pdf
  Drawing Coverage histogram (bin size: 24x)...
    saved to /home/czh/Desktop/01_denovo_genome/N46_DENOVO/quast_out/basic_stats/coverage_histogram.pdf
  Drawing spades_assembly coverage histogram (bin size: 24x)...
    saved to /home/czh/Desktop/01_denovo_genome/N46_DENOVO/quast_out/basic_stats/spades_assembly_coverage_histogram.pdf
Done.

NOTICE: Genes are not predicted by default. Use --gene-finding or --glimmer option to enable it.

2019-11-04 06:15:08
Creating large visual summaries...
This may take a while: press Ctrl-C to skip this step..
  1 of 2: Creating Icarus viewers...
  2 of 2: Creating PDF with all tables and plots...
Done

2019-11-04 06:15:08
RESULTS:
  Text versions of total report are saved to /home/czh/Desktop/01_denovo_genome/N46_DENOVO/quast_out/report.txt, report.tsv, and report.tex
  Text versions of transposed total report are saved to /home/czh/Desktop/01_denovo_genome/N46_DENOVO/quast_out/transposed_report.txt, transposed_report.tsv, and transposed_report.tex
  HTML version (interactive tables and plots) is saved to /home/czh/Desktop/01_denovo_genome/N46_DENOVO/quast_out/report.html
  PDF version (tables and plots) is saved to /home/czh/Desktop/01_denovo_genome/N46_DENOVO/quast_out/report.pdf
  Icarus (contig browser) is saved to /home/czh/Desktop/01_denovo_genome/N46_DENOVO/quast_out/icarus.html
  Log is saved to /home/czh/Desktop/01_denovo_genome/N46_DENOVO/quast_out/quast.log

Finished: 2019-11-04 06:15:08
Elapsed time: 0:00:11.039108
NOTICEs: 1; WARNINGs: 0; non-fatal ERRORs: 0

Thank you for using QUAST!

quast评估结果可以看出unicycler组装结果优于spades，unicycler应该是矫正和剔除了长度太短和一些冗余的contig，减少了contig的数量，提高了N50的长度，但基因组总长度小于spades组装。

5. unicycler调用spades初步组装后对contigs的处理。

区分单拷贝的contigs和重复的contigs，选择合适的单拷贝contigs进行下一步处理。

剔除contig间的重叠区域，为了尽可能的环化contigs。

单拷贝contigs桥接图

unicycler相较于spades耗时较长，组装结果较好，这只是一个菌基因组的组装比较结果，可能其他菌株使用unicycler组装并不会有较大提升，反而会降低N50，不过N50不是评估基因组的唯一标准，过于苛求更长的N50可能会导致错误组装，可能会误导后续基因注释和功能研究。