发布时间 : 2024/5/15 20:38:25 星期三 文章组装组参考文档version1更新完毕开始阅读6320b48f9b89680202d82527

附录二 纠错流程

1 纠错策略

建立高频kmer表.。

遍历reads, 在一个read上找到一个连续高频kmer最多的一个区域，做为下一步动态规划的起始点S。

从S向左向右遍历，如下图所示：

read S Kmer ,X=A kmer??,X??=A, X??=A … kmer??,X??=T, X??=A … kmer??,X=T kmer??,X??=A, X??=T … kmer??,X??=T, X??=T

以向左遍历为例：在S左端取n-1长度的序列Ni，其中n为kmer的长度。当Xi依次

为A,C,G,T四种碱基时，考察kmeri= Xi+Ni是否为高频kmer,， 若kmeri是高频则存入table，初始化其上游指针为NULL，若当前Xi的值与对应位点read上的碱基一致，则其scorei=0，否则scorei=1。在kmeri左端取n-1长度的序列Ni?1，考察kmeri?1= Xi?1+Ni?1

是否为高频 ,， 若kmeri?1是高频则存入table，初始化其上游指针指向kmeri，若当前Xi?1的值与对应位点read上的碱基一致，则其scorei?1= scorei，否则scorei?1= scorei+1。继续上述迭代，迭代结束后，遍历table找到一条分数最低的路径作为最优解。

2 纠错流程

程序路径：

/share/data2/shizhb/workspace/MySoftware/assembly/correction/ 流程：kmerfreq -> unicorn -> merge_pair_lst.pl 1. kmerfreq

使用kmerfreq程序读入所有需纠错的文件（尽量挑取高质量无偏向性数据），统计kmer频数，生成频数表。 Usage:

kmerfreq [options]

-i 需要纠错的reads文件列表 -o 输出的文件前缀

-q quality cutoff (默认为5,可以不变) 碱基质量值低于该值的kmer 去除

-s seed length (默认为17,需16G内存，可以不变)

-n output kmer index?: 0: no, 1: yes. (默认为0,输出一个16G大小的kmer频数表，unicorn要用，如果是小数据集可以设为1，注意unicorn也有这个参数，必须一致。)

-f file format: 1: fq, 2: fa.

-l is there have file list?: 0: no, 1: yes. (默认为1,把需要纠错的文件全路径保存在这个file list文件中，如果只有一个文件需要纠错，可以设为0)

-h -? Help

脚本 /share/data2/shizhb/workspace/MySoftware/assembly/correction/kmerfreq -i kmer.lst -o ant -f 1 -l 1 -n 0 >kmer.log 输出文件 ant.stat 和ant.freq 樊老师新纠错程序：

/ifs1/GAG/assemble/fanw/Assembly/source_codes/correct_error/correct_error_v1.0/

2. unicorn

使用unicorn程序读入需纠错的文件和频数表进行纠错。 Usage:

unicorn [options]

-i 需要纠错的reads 文件list

-r 指定kmerfreq 输出的kmer频数表*freq

-n input kmer index?: 0: no, 1: yes. (和kmerfreq保持一致)

-k start of kmer frequence cutoff (默认为5，可以不变) 纠错前，去除kmer频数低于5的kmer

-e end of kmer frequence cutoff (默认为5，可以不变) 纠错后去除kmer频数低于5的kmer

-d set delta's value (默认为2，可以不变) reads中有超出2个以上碱基错误的不进行纠错

-s seed length (默认为17,16G内存，和kmerfreq保持一致，可以不变)

-t thread number （默认为4个） -f file format: 1: fq, 2: fa.

-l is there have file list?: 0: no, 1: yes. (这个file list可以和前面一样，也可以将前面的list拆分成几份，并行的纠错。) -h -? Help

脚本：/share/data2/shizhb/workspace/MySoftware/assembly/correction/unicorn -i corr.lst -r ant.freq -n 0 -f 1 -l 1 >corr.log 输出文件：corr.lst中的文件纠错后的*.corr文件。

3. 合并

使用merge_pair.pl或merge_pair_lst.pl程序将pair reads合并成一个文件。用merge_pair_lst.pl进行合并时 merge_pair_lst.pl和merge_pair.pl须放在同一目录下。 脚本：perl /share/data2/shizhb/workspace/MySoftware/assembly/correction/merge_pair_lst.pl corr.lst

corr.lst文件中，纠错后的reads1和reads2文件放在一起，read1在前，read2在后。

3 内存及时间损耗

kmerfreq程序kmer等于17mer的时候占用内存16G。 unicorn程序kmer等于17mer的时候占用内存16G。另外每个线程在处理一个文件的时候需要将该文件的所有reads读入内存。现在reads一般长度75bp， reads name 75个字节，每个文件至少25M个reads，那么一个文件要占4G左右内存。每个线程还有单独的动态规划表占1G内存。一个线程5G，unicorn程序默认开设4个线程要占36G。 merge_pair_lst.pl 程序所耗内存很少

kmerfreq程序统计kmer频数，输出频数表的耗时跟文件的多少和io状况有关。处理一个文件约需100s，african总共606个文件耗时15h。

unicorn程序可将所有需纠错文件拆分处理，处理一个文件约需1000s，100个文件4个线程耗时1000s*100/4 = 25000s = 7h。african 606个文件拆分成6份耗时7h。 African纠错总共耗时22h。

参考文献：http://soap.genomics.org.cn/about.html

附录三 Kmer分析总论 1 引言 Kmerfreq程序说明 程序功能： 统计测序read中kmer频数。 参数介绍： -k 设定kmer长度，建议设定为奇数，默认17，占用内存16G。 -a 将read从头截取长度。 -d 将read从尾截取长度。这两个参数根据实际测序质量进行调整。 -r 设定固定read长度，每个read仅仅截取该固定长度，用于提取kmer。建议-r和-d不要同时使用。 -n 设定输出最小kmer深度。只输出深度大于该值的深度频数列表。 -t 设定总碱基数上限，读取碱基数达到该值时不再读取文件。 基因组大小估计： 获得深度分布曲线和深度乘积曲线各自峰值深度，综合考虑估计准确峰值位置，然后根据公式：G=kmer_num/kmer_depth，估计基因组大小。 通过深度1处频数估计错误率，特异kmer数参考估计基因组大小。 输出图表： xx 17-mer depth distribution543210161116212631364146 1num kmer_num pkdepth genome_size used_base X node_num Xxx 306358 999579000 8 124947375 1189975000 9.5 101166945 其中1num为深度为1处的kmer频数，其余均由程序自身生成（*.log文件）。