一、ALLPATH简介
ALLPATHS-LG是一个基因组组装软件,适合于组装short reads数据,由Computational Research and Development group at the Broad Institute开发。ALLPATHS-LG是现在行业内公认进行基因组De novo组装效果最好的软件。
二. 基础注意事项
1. 不能只使用一个library数据进行组装; 2. 必须有一个"overlapping"的片段文库的paired-reads数据。比如,reads长度~ 100bp,插入片段库长度~180bp3. 必须有jumping library数据; 4. 基因组组装需要100x或以上基因组覆盖度的碱基,这个覆盖度是指raw reads数据(在 error correction和filtering之前)的覆盖度; 5. 可以使用PacBio数据; 6. 不能使用454数据和Torrent数据。主要是这两者测序太贵,如果什么时候价格降低,有 需求的话,会写出相应的代码来满足要求; 7. 官方提供了测试用数据; 8. 不支持在整个计算机集群上进行运算; 9. 需要消耗的内存峰值大约是1.7bytes每个碱基,即输入10G的碱基数据量,大约需要17 G内存; 10. 对于试探性的参数,比如K,原则上可以调整。但是我们不会自行调整,并也不推荐。AL LPATHS-LG不像其它De novo一样,Kmer大小的参数K和read大小之间没有直接的联系, ALLPATHS-LG会在运行过程中运用一系列的K值。
三. ALLPATHS-LG使用方法
1. 基础的使用方法和命令
使用RunAllPathsLG这个命令来运行。虽然有很多参数,但是在没有指导的情况下不要随意使用,使用默认设置即可。其使用方法为:
$ RunAllPathsLG arg1=value1 arg2=value2 ...
参数主要是设置程序辨别的一些目录,在程序的运行过程,会输入相应目录中的数据,将结果输入到指定的目录。一个简单的命令使用例子:
#!/bin/sh # ALLPATHS-LG needs 100 MB of stack space. In 'csh' run 'limit stacksize 100000'. ulimit -s 100000 # ALLPATHS-LG命令的写法与一般的linux参数写法不是很一样。采用 ‘参数=值’ 的方法,并使之成每行一个参数,使用'\'来连接各个参数,这样看起来直观易懂。初始接触的人可能会不适应。 RunAllPathsLG \ PRE=$PWD\ REFERENCE_NAME=species.genome\ DATA_SUBDIR=data\ RUN=run\ SUBDIR=test\ EVALUATION=STANDARD\ TARGETS=standard\ OVERWRITE=True\ MAXPAR=8 | tee -a assemble.out
2. 详细的参数说明
必须的参数 PRE (String) 程序运行的根目录,所有的其它目录全在该目录下REFERENCE_NAME (String) 参考基因组目录名称,位于PRE目录下。如果有一个参考基因组,可将参考基因组放到该 目录中;若没有,则创建该文件夹用于基因组组装DATA_SUBDIR (String) DATA子目录名称,位于REFERENCE_NAME目录下。程序从该目录中读取数据。 RUN (String) 运行目录名称,位于DATA_SUBDIR下。程序将生成的中间文件和结果文件存储于该目录 。比如组装结果是一个名为ASSEMBLES的目录,位于该目录下。 部分可选参数: SUBDIR (String) default: test 子目录名,在REF/DATA/RUN/ASSEMBLIES目录下创建的存放基因组组装结果的目录 名。 K (int) default: 96 核心Kmer大小,只有K=96能很好地运行。 EVALUATION (String: {NONE,BASIC,STANDARD,FULL,CHEAT})default:BASIC 给定一个参考基因组,pipeline能在基因组组装的不同阶段对组装过程和结果进行评估。 BASIC:基础评估,不需要参考基因组; STANDARD:使用参考基因组来运行评估模块; FULL:在某些组装模块下打开in-place评估,不会影响组装结果; CHEAT:稍微使用参考基因组指导组装,产生更详细的分析,能对组装结果产生小的(好方 向的)改变。REFERENCE_FASTA (String) default: REF/genome.fasta 评估中使用的参考基因组。 MAXPAR (int) default: 1 有些模块的运行是独立的,不相互依赖,能同时运行。该参数设定能同时运行的模块的最 大数目。由于pipeline中的绝大部分模块都能多线程运行,因此将该值设定大于1,效果不明 显。 THREADS (String) default: max 有些模块能多线程程运行,默认使用最大线程数运行。 OVERWRITE (Bool) default: False 是否覆盖存在的文件。可以设置该选项为True,在每次运行程序的时候设定RUN参数为 一个新的目录名,则比较好。 TARGETS (vec) default: standard pipeline会生成一系列的文件,不同的文件的生成需要call不同的模块。如果某文件 已经存在了并且是最新的,则跳过相应的模块的运行。本参数指定生成哪些拟定的目标文件(p seudo targets)。若目标文件没有相应的模块能生成,则会得到报错。 none:没有拟定的目标文件,仅仅生成指定的目标文件; standard:生成组装文件和选定的评估文件; full_eval:生成组装文件和额外的评估文件。TARGETS_REF (String) 在ref_dir目录中生成的目标文件。 多个目标文件的书写方法为: TARGETS_REF="{target1,target2,target3}" 。 TARGETS_DATA (String) 在data目录中生成的目标文件。 TARGETS_RUN (String) 在run目录中生成的目标文件。 TARGETS_SUBDIR (String) 在subdir中生成的目标文件。FORCE_TARGETS (Bool) default: False 生成目标文件,即使文件已经存在并且看起来是很新的。
3. 输入文件与目录的准备
两个文库:插入片段长度为180bp和3000bp,illumina测序文件结果为fastq格式。以此为例来准备ALLPATHS-LG运行所需的文件和目录。
(1) 准备 in_groups.csv 和 in_libs.csv 文件。
这两个文件内容由逗号隔开,in_groups.csv文件内容如下:
group_name, library_name, file_name firest, Illumina_180bp, seq/species_500bp_read?.fastq second, Illumina_3000bp, seq/species_3000bp_read?.fastq
in_groups.csv文件的解释:
group_name:数据独特的代号,每一份数据有一个代号; library_name:数据所属文库的名字,体现出该; filename:数据文件所存放位置。可以为相对位置,文件名可以包含'*'和'?'(但是扩展名 中不能有该符号,因为要根据扩展名识别文件类型),从而代表paired数据。支持的文件类型有 '.bam','fasta','fa','fastq','fq','fastq.gz'和'fq.gz'。
in_libs.csv文件内容如下:
library_name, project_name, organism_name, type, paired, frag_size, frag_stddev, insert_size, insert_stddev, read_orientation, genomic_start, genomic_end Illumina_180bp, species, species.genome, fragment, 1, 180, 10, , , inward, 0, 0 Illumina_3000bp, species, species.genome, jumping, 1, , , 3000, 500, outward, 0, 0
in_libs.csv文件的解释:
library_name:和in_groups.csv中的相匹配; project_name:project的名字; organism_name:测序物种的名字; type:仅仅只是一个信息; paired:0:Unpaired reads1:paired readsfrag_size:小片段文库插入片段长度的均值; frag_stddev:小片段文库的插入片段长度估算的标准偏差; insert_size:大片段文库插入片段长度的均值; insert_stddev:大片段文库插入片段长度估算的标准偏差; read_orientation:reads的方向,小片段文库为inward,大片段文库为outward; genomic_start:reads从该位置开始,读入数据,如果不为0,之前的碱基都被剪掉; genomic_end:reads从该位置开始,停止读入数据,如果不为0,之后的碱基都被剪掉。
(2) 使用PrepareAllPathsInputs.pl来对数据进行转换
ALLPATHS-LG接受的输入数据要求如下:
1. ALLPATHS-LG的输入数据支持小片段文库(fragment library)、大片段文库(jum ping library)和超大片段文库(long jumping library)。并且前两种文库至少各有 一个才能进行基因组组装。超大片段文库是只插入片段>20kb的文库,其测序方向和小片段文 库一致,为inward。 2. ALLPATHS-LG的输入数据放置在//文件夹下,包含3种文件:碱基文件,质量文件和配 对信息文件 frag_reads_orig.fastb frag_reads_orig.qualb frag_reads_orig.pairs jump_reads_orig.fastb jump_reads_orig.qualb jump_reads_orig.pairs 以下是可选的超大插入片段文库对应的数据文件(非必须): long_jump_reads_orig.fastb long_jump_reads_orig.qualb long_jump_reads_orig.pairs
使用PrepareAllPathsInputs.pl来将fastq等格式的测序结果转换成ALLPATHS-LG可接受的文件。以下是该程序的参数:
DATA_DIR 将转换后的数据文件放到此文件夹下。 PICARD_TOOLS_DIR 若输入数据为bam格式,则需要用到Picard软件,该参数Picard的路径 IN_GROUPS_CSV 输入的in_groups.csv文件名 IN_LIBS_CSV 输入的in_libs.csv文件名INCLUDE_NON_PF_READS default: 1 1:包含non-PF reads;0:仅仅只包含PF reads. PHRED_64 default: 0 0:碱基质量是ASCII的33到126,一般情况下Illumina数据的最低碱基质量是'B'1:碱基质量的ASCII码是从64到126,一般情况下Illumina数据的最低碱基质量是'#'。 PLOIDY 生成ploidy文件。该文件就包含一个数字 1 或者 2 。1表示基因组为单倍体型,2表 示双倍体型。 HOSTS 列出平行forking的host主机(这些主机必须要能无密码直接ssh连上)。比如“2,3. host2,4.host3"表示使用本地机器的2个CPU线程,host2机器的3个CPU线程和host3机 器的4个CPU线程。 以下是不常用的参数,主要用来选择转换的数据量的大小。当测序数据量太多,而只想使用其 中一部分数据的时候,可以用到 FRAG_FRAC 使用小片段库reads的比例。比如 30% 或 0.3 。如果设定了此值,则不能同时设定 FRAG_COVERAGE。 JUMP_FRAC 使用大片段库reads的比例。比如 20% 或 0.2 。如果设定了此值,则不能同时设定 JUMP_COVERAGE。 LONG_JUMP_FRAC 使用超大片段库reads的比例。 比如 90% 或 0.9 。如果设定了此值,则不能同时 设定LONG_JUMP_COVERAGE。 GENOME_SIZE 估计的基因组大小,用来计算对应覆盖度所对应的reads数 FRAG_COVERAGE 所期望的小片度库的覆盖度,比如 45. 要求GENOME_SIZE有设定 JUMP_COVERAGE 所期望的大片度库的覆盖度,比如 45. 要求GENOME_SIZE有设定 LONG_JUMP_COVERAGE 所期望的超大片度库的覆盖度,比如 1. 要求GENOME_SIZE有设定
(一)安装bowtieBowtie可以在个人计算机上使用,也可以在CSC服务器上使用终端连接。请参阅以下文档的第一部分,了解如何在笔记本电脑上安装Bowtie。特别是对他们的计算机没有管理员权限的那些应该确保软件的正确安装和功能。Bowtie也可以在服务器计算机上远程使用。我们将提供临时帐户访问CSC,但你将需要一个安全Shell终端程序进行通信。默认情况下,Mac和Linux上都有这样的程序,但需要安装Windows。普遍的实现是PuTTY。即使终端程序不用于读取映射,也将需要其他练习,并且应该可用。Bowtie的安装:从下载页面下载相应的版本(Linux,Mac或Win,小编使用的是在Linux下进行)。将zip文件解压缩到新的目录中,并转到该目录。下载的bowtie包装包含大肠杆菌基因组的预先建立的指数,以及从该基因组模拟的一组1000个35bp的读数。要使用Bowtie对齐这些读取,请键入以下命令。bowtiee_colireads/e_coli_1000.fqmap_result.txt
如果你收到错误消息"commandnotfound",请尝试在"bowtie"(./bowtie)之前添加"./"。
(二)使用Bowtie
(1)Mapping
要使用Bowtie对齐示例读取,请发出以下命令。bowtiee_colireads/e_coli_1000.fqmap_result.txt
如果你收到错误消息"commandnotfound",请尝试在"bowtie"(./bowtie)之前添加"./"。"e_coli"与"indexes/e_coli"相同。你可以在文本编辑器中打开map_result.txt。每行都是一个读取对齐。对齐读取的名称显示在第一列中。对于Mac和Linux,使用"少"会更好。
lessmap_result.txt#extrareading
ReadthemanualinthefolderorwebsitetogetadeeperunderstandinghowBowtieworksandfurtheroptionsinBowtie.
我们来看看Bowtie在1中使用的一些不同的选项,报告所有有效的对齐方式与一些不匹配。
./bowtie-a-v2e_coli--suppress1,5,6,7-cATGCATCATGCGCCAT-a/--all报告每个读取或对的所有有效对齐(默认值:off)
-v
最多不相匹配的报告对齐
-c
查询序列在命令行
--suppress
上以默认输出模式抑制输出列
2限制对齐
$./bowtie-k3-v2e_coli--suppress1,5,6,7-cATGCATCATGCGCCAT-k
每次读取或配对时报告有效对齐(默认值:1)。
3不匹配排名
$./bowtie-a--best-v2e_coli--suppress1,5,6,7-cATGCATCATGCGCCAT
所有相同的对齐方式按最佳到最坏的顺序进行报告
4只有最不匹配
$./bowtie-a--best--strata-v2--suppress1,5,6,7e_coli-cATGCATCATG
(2)配对对齐
当使用-1和-2选项指定正确配对的读取文件时,Bowtie可以对齐配对端读取(对于原始,FASTA或FASTQ读取文件)
./bowtiee_coli-1reads/e_coli_1000_1.fq-2reads/e_coli_1000_2.fqmap_paired.txt
SAMtools(http://samtools.sf.net)是一套用于存储, *** 纵和分析对齐方式的工具,例如Bowtie输出的对齐方式。bowtie-Se_colireads/e_coli_1000.fqec.sam
我们可以再次检查sam文件以查看与txt文件的区别(也是在r4,r5中未映射的读取)。接下来,我们将SAM文件转换为BAM以准备排序。
samtoolsview-bS-oec.bamec.sam
接下来,我们对BAM文件进行排序,
samtoolssortec.bamec.sorted
这样我们就简单的对bam文件中的基因组进行配对对齐。
欢迎分享,转载请注明来源:内存溢出
评论列表(0条)