1. 前期准备
背景信息:
- GC含量 和 GC分布
- 基因组重复程度
- 基因组大小估计
- 杂合情况
最好的情况是对方能提供已经发表的近源物种。根据近源物种分析以上信息,尤其是GC含量以及对应的GC分布,重复程度。
2. 测序策略
根据基因组大小和具体情况选择个大概的k值,根据“测序X数推导说明.pdf”制定用于构建contig所需的数据量以及所需的构建的文库数量。对于植物基因组一般考虑的是大kmer(>31),动物的话一般在27左右,具体根据基因组情况调整。需要在短片段数据量达到20X左右的时候进行kmer分析。Kmer分析正常后,继续加测数据以达到最后期望的数据量。
3. 组装流程
原始数据-数据过滤-纠错-kmer分析-denovo组装
3.1 数据过滤:
/nas/GAG_01A/assembly/Database/Assembly/Package/Filter_data/目录下有程序、源代码、使用文档、test实例
/nas/GAG_01A/assembly/yanglinfeng/Filter_gz/目录下程序用法和上面一样,读gz压缩文件,省去解压缩
3.2 数据纠错:
/ifs1/GAG/assemble/fanw/Assembly/source_codes/correct_error/correct_error_v1.0/有先使用多线程版本correct_error_pread
说明文档以及算法详见“Error_correction_algorithm.doc”
3.3 kmer分析:
/nas/GAG_01A/assembly/Database/Assembly/Package/kmerfreq/目录下有程序、源代码、使用文档、
/ifs1/GAG/assemble/lizhenyu/kmer/kmerfreq2buff/kmerfreq 多线程版本,原理与上面程序一致,程序帮助信息包含使用实例
原理说明可参见“kmer分析.docx”
Kmer 分析中估计基因组大小采用纠错后的数据。
3.4 SOAP denovo(grape)组装:
/nas/RD_01A/xieyinlong/bin/目录下grape63mer大kmer版本,可设k>31,grape1123设置K<=31,其他在使用没有差别
/nas/RD_01A/zhuhm/bin/
使用说明以及一些参数选择上具体参见“SOAP denovo使用说明”。
grape组装并非一次就能得到理想的结果,会根据以有的组装结果做分析,调整参数,处理数据,加测少量数据等策略来得到比较理想的结果,在下面的说明中进行详细解释。
4. 基因组的grape组装输出以及soap比对分析
4.1 Contig构建质量
Contig构建的质量情况是组装的基础,较为理想的contig结果才能得到比较好的scaffold结果,对于正常的项目得到的congtig N50在1k左右,而N90在200bp左右(100PE测序)。如果是正常的基因组在kmer分析中没有出现异常的情况,那么需要查看是否由于过多的测序错误造成的,即kmer频数<3所占的比例过高造成的,那么对于数据纠错是否有效执行,对于一些情况下,经过纠错后的数据仍然可以设置-d参数,另外对于其他参数的选择是否合理,k值大小(与数据量、基因组特性相关),M、R参数的使用。
4.2 文库插入片段大小的校正
利用每个文库的一个lane与组装得到的基因组进行soap比对,利用程序画出比对的插入长度分布,调整配置文件中的avg_ins
另外对于这个分析还能得到大片段文库的可用数据的比率,有助于估计实际组装可用的大片段数据量,详见“大片段文库质量标准.doc”。
4.3 基因组GC_depth分布图
根据GC_depth分布图可以看出测序是否有明显的GC偏向,是否有存在样品被细菌污染等情况。做GC_depth图需要用到所有短片段reads的soap结果。
污染:高GC或者GC与正常整体有较大差异,覆盖度低。取出异常区域
与nt库做blast比对,确定可疑的污染细菌的全基因组,用测序得到的reads和可疑的污染细菌的全基因组做soap比对,取出能比对上的reads后重新组装。
/ifs1/GAG/assemble/wangzhiwen/raid6/bin/gc_vs_depth/coverage_gc_dis.sh
运行脚本会得到使用说明
4.4 基因组的单碱基深度分布
较为理想的分布是接近泊松分布,这样说明测序的随机性较好。计算基因组单碱基深度用到所有reads的soap的结果,程序如下,直接运行的到帮助信息。
/nas/GAG_01A/assembly/Database/Assembly/chenyan/01.Data_analysis/soap.coverage
4.5 确定污染,去除污染reads
将在疑似污染区域中的scaffold,scaffold 80%落在疑似区域,将scaffold断成contig,取片段>1k或者更长的部分(根据实际情况选择)与nt库做blast比对,取比对结果最好的以确定是否为污染,污染物为什么物种,使用注释组流程:
/nas/GAG_02/huangqf/GACP-8.0/02.gene-annotation/Dbxref/bin/blast_database.pl
下载确定的污染物的基因组序列作为reference,将所有的reads与之比对,能有比上的就去掉对应的一对reads,过滤后的reads再进行重新组装。
5. 确定组装版本后的处理以及检验
5.1 原始contig定位
主要是针对原始contig的准确度远高于补洞部分的序列,这样在后期注释、进化可以区别对待。
将原始contig定位到补洞后scaffold上,选取长度大于cutoff的contig,利用的是nucmer的比对,根据比对结果进行定位,比上位置为大写的ATCG,其余位置为小写atcg。输出文件为*.OUT
相应脚本和程序目录:
/nas/GAG_01A/assembly/Database/Assembly/yanglinfeng/confirm_contig2scaffold/
使用说明请参考目录下README.doc
5.2 针对关键区域的PE关系数的确定
针对注释或者进化对某些scaffold区域的分析得到较为重要的结果,需要进一步确定是否存在组装错误问题。
从之前的所有reads的soap结果中选取对应的scaffold比对结果,利用程序
/ifs1/GAG/assemble/wangzhiwen/raid6/bin/pair_ends_map/version3.0/pair_end_in_scaffold.pl
得到可是的PE关系数的图:
因为得到的整条scaffold的图片比较大,window自带的照片库浏览器看不了,推荐一个看图片的软件,“ACDSee 9 Photo Manager”
6. 基因组的组装版本的EST和BAC评价
6.1 EST评价
参考/nas/GAG_01A/assembly/yanglinfeng/EST_and_BAC/EST_work.sh
6.2 BAC评价
参考/nas/GAG_01A/assembly/yanglinfeng/EST_and_BAC/BAC_pipeline.pl
使用文档请查看:
/nas/GAG_01A/assembly/yanglinfeng/EST_and_BAC/BAC_pipeline.README.txt
7. 其他使用到的程序以及说明
7.1 Connect Mated Read(CMR)
用于将测通的PE reads连接成一条长reads。
/nas/GAG_01A/assembly/Database/Assembly/Package/cmr/cmr
运行得到帮助信息,说明文档在目录:
/nas/GAG_01A/assembly/Database/Assembly/Package/cmr/
7.2 simulate_solexa_reads
用于模拟solexa reads,运行得到帮助信息。
/nas/GAG_01A/assembly/Database/Assembly/Package/simulate_solexa_reads/simulate_solexa_reads
说明文档在目录:
/nas/GAG_01A/assembly/Database/Assembly/Package/simulate_solexa_reads/
7.3 项目目录结构
对于项目的目录结构,每个项目下分6个子目录。如:
/ifs1/GAG/assemble/chenyan/cucumber/目录里面包含有
00.data 01.analysis 02.assembly 03.soap 04.evaule 05.super-scaffold 06.backup
00.data主要放置原始数据链接,处理后数据,Sanger测序数据
01.analysis主要是针对初步产出的20X的分析,包括Kmer分析和细菌污染比对分析
02.assembly 主要是保存组装的版本
03.soap主要是需要分析文库的插入片段,大片段数据去除小峰,reads数据去除污染,还有GC—Depth的分布等等
04.evaluate 主要是对组装结果的评价,包括EST,BAC评价,全基因组线性化的比对(有reference序列) 和两个组装版本间的比较
05.super-scaffold 主要是针对有fosmid-end和bac-end数据的项目
06.backup 主要是备份中期传给客户的报告,序列,还有上述需要备份的脚本和配置文件。
文章来源:基因组所培训