主要参考这篇文章的图4：http://www.nature.com/ng/journal/v42/n11/fig_tab/ng.691_F4.html
这是2010年发表于nature genetics杂志的Whole-genome sequencing and comprehensive variant analysis of a Japanese individual using massively parallel sequencing 虽然文章选择的是SOAPdenovo,ABySS,Velvet这3款软件来进行组装，但毕竟是2010年的文章了，现在其实有更好的选择，比如Minia
选择Minia工具来组装
Minia软件也是基于de Bruijn图原理的短序列组装工具，优于以前的ABySS和SOAPdenovo，所以这里就选择它啦。
下载安装Minia
安装官网的指导说明书下载二进制版本即可，代码如下：
## Download and install Minia
# http://minia.genouest.org/
cd ~/biosoft
mkdir Minia && cd Minia
wget https://github.com/GATB/minia/releases/download/v2.0.7/minia-v2.0.7-bin-Linux.tar.gz
tar -zxvf minia-v2.0.7-bin-Linux.tar.gz
~/biosoft/Minia/minia-v2.0.7-bin-Linux/bin/minia --help
## eg: ./minia -in reads.fa -kmer-size 31 -abundance-min 3 -out output_prefix
软件使用方法也非常简单，就一行命令，其中最佳-kmer-size需要用KmerGenie来确定。
使用
step1:提取比对失败的reads
samtools view -f4 jmzeng_recal.bam |perl -alne '{print "\@$F[0]\n$F[9]\n+\n$F[10]" }' >unmapped.fq
perl ~/biosoft/PRINSEQ/prinseq-lite-0.20.4/prinseq-lite.pl -verbose -fastq unmapped.fq -graph_data unmapped.gd -out_good null -out_bad null
perl ~/biosoft/PRINSEQ/prinseq-lite-0.20.4/prinseq-graphs.pl -i unmapped.gd -png_all -o unmapped
perl ~/biosoft/PRINSEQ/prinseq-lite-0.20.4/prinseq-graphs.pl -i unmapped.gd -html_all -o unmapped
cd ~/data/project/myGenome/gatk/jmzeng/unmapped
共31481084/4=7870271，仅仅是7.8M的reads
step2: 用KmerGenie确定kmer值
KmerGenie estimates the best k-mer length for genome de novo assembly.
KmerGenie predictions can be applied to single-k genome assemblers (e.g. Velvet, SOAPdenovo 2, ABySS, Minia).
## http://kmergenie.bx.psu.edu/
cd ~/biosoft
mkdir KmerGenie && cd KmerGenie
wget http://kmergenie.bx.psu.edu/kmergenie-1.7044.tar.gz
tar zxvf kmergenie-1.7044.tar.gz
cd kmergenie-1.7044
make
python setup.py install --user
~/.local/bin/kmergenie --help
cd ~/data/project/myGenome/gatk/jmzeng/unmapped
~/.local/bin/kmergenie unmapped.fq
step3: 运行Minia
cd ~/data/project/myGenome/gatk/jmzeng/unmapped
~/biosoft/Minia/minia-v2.0.7-bin-Linux/bin/minia -in unmapped.fq -kmer-size 31 -abundance-min 3 -out output_prefix
7.8M的reads组装之后有272007条contigs
组装之后：
Prinseq v0.20.4 was used to calculate assembly statistics, including N50 contig size, GC content
cd ~/data/project/myGenome/gatk/jmzeng/unmapped
perl ~/biosoft/PRINSEQ/prinseq-lite-0.20.4/prinseq-lite.pl -verbose -fasta output_prefix.contigs.fa -graph_data contigs.gd -out_good null -out_bad null
perl ~/biosoft/PRINSEQ/prinseq-lite-0.20.4/prinseq-graphs.pl -i contigs.gd -png_all -o contigs
perl ~/biosoft/PRINSEQ/prinseq-lite-0.20.4/prinseq-graphs.pl -i contigs.gd -html_all -o contigs
perl ~/biosoft/PRINSEQ/prinseq-lite-0.20.4/prinseq-lite.pl -verbose -fasta output_prefix.contigs.fa -stats_assembly
就是给出一些指标，如下；
stats_assembly N50 176 stats_assembly N75 113 stats_assembly N90 78 stats_assembly N95 70
Input Information
| Input file(s): | output_prefix.contigs.fa |
|---|---|
| Input format(s): | FASTA |
| # Sequences: | 272,007 |
| Total bases: | 44,868,011 |
Length Distribution
| Mean sequence length: | 164.95 ± 204.44 bp |
|---|---|
| Minimum length: | 63 bp |
| Maximum length: | 10,187 bp |
| Length range: | 10,125 bp |
| Mode length: | 150 bp with 16,461 sequences |
然后用RNA-SEQ数据来比对验证！ 以后再讲
把组装好的contigs拿去NCBI做blast看看物种分布,Distribution of top nucleotide BLAST hits by species from the NCBI nr database for 1000 random contigs in the assembly！其实上面的prinseq软件也简单的给出了一个污染物种分布情况表，但是这个原理不一样。以后再讲