在前面的直播基因组系列，我们讲解过那些比对不少我们人类的参考基因组序列的数据，其实可以细致的进行探究。
直播】我的基因组（十五）:提取未比对的测序数据
这里主要参考这篇文章的图4：http://www./ng/journal/v42/n11/figtab/ng.691 F4.html
组装的contig注释到物种 这是2010年发表于nature genetics杂志的Whole-genome sequencing and comprehensive variant analysis of a Japanese individual using massively parallel sequencing 虽然文章选择的是SOAPdenovo,ABySS,Velvet这3款软件
选择Minia工具来组装 Minia软件也是基于de Bruijn图原理的短序列组装工具，优于以前的ABySS和SOAPdenovo，关键是速度非常快，十几分钟就OK了，不消耗计算机资源，所以这里就选择它啦。
下载安装Minia 安装官网的指导说明书下载二进制版本即可，代码如下：
## Download and install Minia
# http://minia./
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
使用 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
Input Information
Input file(s): unmapped.fq Input format(s): FASTQ # Sequences: 7,870,271 Total bases: 1,180,540,650
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./
cd ~/ biosoft
mkdir KmerGenie && cd KmerGenie
wget http : //kmergenie.bx./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软件也简单的给出了一个污染物种分布情况表，但是这个原理不一样。 以后再讲
