学一点儿基因组组装

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

推荐值k=59
k=59时的kmer频数分布

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
...
...

k=17

~/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

k=17, 好看多了

在低深度区的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图上的点是否聚集得比较好（理论上应该聚在一起，如果明显分成几类，说明可能存在污染）。