1.下机数据质量控制
主要是针对低质量的reads和含有adapter的reads，我事先并不知道adapter序列，所以就只过滤了低质量的reads。采用的是华大开发的SOAPnuke(v1.5.6)。
SOAPnuke1.5.6 filter -d -S -n 0.1 --qualRate 0.3 --lowQual 12 -Q 2 -G -1 R1.fq.gz -2 R2.fq.gz -C out.R1.fq.gz -D out.R2.fq.gz
2.选择K-mer和基因组大小的评估
2.1选一个合适的K
kmergenie weizhi.list -k 71 -l 9 -s 2 -t 12
2.2低深度区的reads纠错
wget https://nchc.dl.sourceforge.net/project/soapdenovo2/ErrorCorrection/SOAPec_v2.01.tar.gz
tar zxvf SOAPec_v2.01.tar.gz
rm -f SOAPec_v2.01.tar.gz
KmerFreq_AR -k 17 -t 10 -p weizhi weizhi.lst > kmerfreq.log 2> kmerfreq.err
$ ll
total 89M
-rw-r--r-- 1 huangsiyuan grp3 22 Nov 2 23:07 kmerfreq.err
-rw-r--r-- 1 huangsiyuan grp3 1.5K Nov 2 23:09 kmerfreq.log
-rw-r--r-- 1 huangsiyuan grp3 89M Nov 2 23:09 weizhi.freq.cz
-rw-r--r-- 1 huangsiyuan grp3 13K Nov 2 23:09 weizhi.freq.cz.len
-rw-r--r-- 1 huangsiyuan grp3 15K Nov 2 23:09 weizhi.freq.stat
-rw-r--r-- 1 huangsiyuan grp3 867 Nov 2 23:09 weizhi.genome_estimate
-rw-r--r-- 1 huangsiyuan grp3 112 Nov 2 23:00 weizhi.lst
$ less weizhi.genome_estimate
Lowfreq cutoff is: 7
Number of lowfreq Kmers species caused by sequencing errors: 20026657
Ratio in all kmer species: 0.781208
Number of lowfreq Kmers individuals caused by sequencing errors: 21720077
Ratio in all kmer individuals: 0.0809621#这四句能告诉我们什么信息
...
...
**********Summary Table**********#k=17时，估计的G=5993946bp
kmer kmer_num kmer_depth genome_size base_num read_num base_depth(X) average_read_len unique_kmer_number
17 268274496 41.1589 5993946 319374400 3193744 53.2828 100 25635515
$ less weizhi.freq.stat
#KmerSize: 17
#SpaceSeedSize: 0
#Kmer_theory_max_num: 17179869184
#Kmer_individual_num: 268274496
#Kmer_species_num: 25635515
#Filter_kmer_num: 0
#可以用来做kmer频数分布图
#Kmer_Frequence Kmer_Species_Number Kmer_Species_Ratio Kmer_Species_Accumulate_Ratio Kmer_Individual_Number Kmer_Individual_Ratio Kmer_Individual_Accumulate_Ratio
1 18598993 0.725517 0.725517 18598993 0.0693282 0.0693282
2 1252004 0.0488387 0.774355 2504008 0.00933375 0.078662
3 125178 0.00488299 0.779238 375534 0.00139981 0.0800618
...
...
~/learn_assemble/soft/SOAPec_v2.01/bin/Corrector_AR -k 17 -l 3 -Q 33 -t 10 weizhi.freq.cz weizhi.freq.cz.len weizhi.lst >corr.cout 2>corr.cerr
$ ll -t #多了3个文件
total 89M
-rw-r--r-- 1 huangsiyuan grp3 48 Nov 2 23:39 corr.cerr
-rw-r--r-- 1 huangsiyuan grp3 1.3K Nov 2 23:39 corr.cout
-rw-r--r-- 1 huangsiyuan grp3 489 Nov 2 23:39 weizhi.lst.QC.xls
$ less weizhi.lst.QC.xls
==========*.cor.stat:==========
FileName RawReads RawBases ResultReads ResultBases TrimmedReads TrimmedBases DeletedReads BasesFilterRatio EstimatedErrorRatio
out.R1.fq.gz 1596872 159687200 1594244 157816402 138460 1607998 2628 0.0117154 0.0040411
out.R2.fq.gz 1596872 159687200 1594162 157802937 138633 1613263 2710 0.0117997 0.00403655
#然后就能在存放数据的目录下面看到reads校正后重新生成的fq文件
-rw-r--r-- 1 huangsiyuan grp3 136M Nov 2 23:39 out.R1.fq.gz.cor.pair_1.fq.gz#最好重新命个名
-rw-r--r-- 1 huangsiyuan grp3 470K Nov 2 23:39 out.R1.fq.gz.cor.single.fq.gz
-rw-r--r-- 1 huangsiyuan grp3 141M Nov 2 23:39 out.R2.fq.gz.cor.pair_2.fq.gz
-rw-r--r-- 1 huangsiyuan grp3 134M Nov 2 17:26 out.R1.fq.gz
-rw-r--r-- 1 huangsiyuan grp3 140M Nov 2 17:26 out.R2.fq.gz
-rw-r--r-- 1 huangsiyuan grp3 135M Nov 2 16:40 R2.fq.gz
-rw-r--r-- 1 huangsiyuan grp3 129M Nov 2 16:40 R1.fq.gz
现在重新画一下kmer频数分布图
~/learn_assemble/soft/SOAPec_v2.01/bin/KmerFreq_AR -k 17 -t 10 -p weizhi weizhi.lst >kmerfreq.log 2>kmerfreq.err
和前面那个同名文件相比，Ratio降低得非常明显
$ less weizhi.freq.stat
#Kmer_Frequence Kmer_Species_Number Kmer_Species_Ratio Kmer_Species_Accumulate_Ratio Kmer_Individual_Number Kmer_Individual_Ratio Kmer_Individual_Accumulate_Ratio
1 10730 0.00190359 0.00190359 10730 4.06186e-05 4.06186e-05
2 2304 0.000408749 0.00231234 4608 1.74436e-05 5.80622e-05
3 1537 0.000272676 0.00258502 4611 1.7455e-05 7.55172e-05
$ less weizhi.genome_estimate
#这次估计是G=5968689
在低深度区的reads纠错这一步中，我用的都是k=17，主要是因为SOAPec_v2.01中只提供了几种k的情况，此外我觉得k的值不影响对测序错误reads的校正。
我用校正后的fq1，fq2重新用kmergenie预测了一下，这一次k值更大了，变成了65，于是我就很疑惑，kmergenie到底有没有用？k值又应该怎么选取呢？
3.SOAPdenovo2组装
其实到这一步，我还是不确定k应该选多少，kmergenie推荐的值有些大，前面的讨论中k=17画出来的图很好且G更接近真实情况（我事先已经知道G的大概值），文献上有选35的，所以很迷。
最终决定选k=17,27,37,47,57分别组装一下再比较，好在基因组并不大。
~/learn_assemble/soft/SOAPdenovo2/SOAPdenovo-63mer pregraph -s ./my.config -p 4 -K 17 -o ./weizhi > pregraph.log
~/learn_assemble/soft/SOAPdenovo2/SOAPdenovo-63mer contig -g ./weizhi -D 1 -M 3 -p 4 > contig.log
~/learn_assemble/soft/SOAPdenovo2/SOAPdenovo-63mer map -s ./my.config -p 4 -g ./weizhi > map.log
~/learn_assemble/soft/SOAPdenovo2/SOAPdenovo-63mer scaff -p 4 -g ./weizhi -F >scaff.log
可以看出，k=47时两个N50比较大，组装的效果相比较而言更好。当然到这里组装工作并没有做完，仅仅凭借N50也不能给组装的结果评个好坏。注意到k从47到57，scaffoldN50变大了，但Scaffold_Num，Contig_Num，ContigN50均变差了，这一点是比较反常的。于是又选k=49,51,53,55组装了一下。
到这儿我猜测，k从小到大会经过一个合适值，过了这个值，N50的表现就会变差。但这个值并不一定就是最终要选取的值，拿我的这个例子来说，就不能说k=47最好，也不能说以k=47组装出来的基因组最好。还是得看生物学层面的意义，没准儿这个生物本身的基因组就是6.1k而不是5.9k。
所以才需要采取多种评价体系：
- 比如看看组装出来的基因组是否含有公认的核心基因集，涵盖多少；
- 将fq数据比对到组装出来的基因组上，看比对率是多少；
- 看看GC-depth图上的点是否聚集得比较好（理论上应该聚在一起，如果明显分成几类，说明可能存在污染）。
