使用Crossbow作SNP分析

Crossbow:http://bowtie-bio.sourceforge.net/crossbow/index.shtml

在这里讲述本地化使用single computer运行Crossbow的命令详解和流程。

1. Crossbow的安装

在上面提到的Crossbow主页上点击下载链接去sourceforge上下载Crossbow的压缩包,解压即可。

Crossbow的使用需要依赖于bowtie和soapsnp。当使用sra数据的时候,还需要依赖于SRA toolkit中的fastq-dump命令。

$ wget http://kaz.dl.sourceforge.net/project/bowtie-bio/crossbow/1.2.1/crossbow-1.2.1.zip
$ unzip crossbow-1.2.1.zip
$ cd crossbow-1.2.1/soapsnp
$ make
$ sudo cp soapsnp /usr/local/bin/
$ cd ..
$ ./cb_local --test

cb_local是本地化的Crossbow运行命令,–test用于检测 bowtie, soapsnp 和 fastq-dump 这三个程序是否能正常使用。

2. cb_local 命令主要参数

Crossbow可以运行在Amazon Elastic MapReduce (EMR)上,应该是按计算量收费;也可以运行在计算机集群上;或运行在单个的计算机上。其cb_local特有的参数有:

--reference <path>
    参考序列的文件夹。该文件夹由一个参考序列jar文件使用unzip解压缩出来的,里面主
要包括3个子文件夹和一个列表文件:1."index":里面包含bowtie-bulid产生的index文
件;2."sequences":里面包含参考序列;3."snps":已有的参考snp结果;4."cmap.
txt":参考序名和对应的sequences文件夹中的fasta的header名各占一列。

--input <path>
输入的文件或文件夹路径。如果不指定--preprocess 参数,则该参数后接的是一个
文件夹,该文件夹包含了reads的预处理结果,即 --preprocess-output 的输出结果;
如果指定了 --preprocess 或 --just-preprocess 参数,则是一个文本文件,内容
为输入的序列文件的一个清单(Labeled manifest files)。该文件描述了输入序列文件
和样本的信息。序列文件是FASTQ或sra格式文件,FASTQ格式可以是gzip或bzip2压缩。
该文件每一行代表一个测序的结果,是以下是 unpaired reads 和 paired reads 的
书写方法:
<URL>(tab)<Optional MD5>
<URL-1>(tab)<Optional MD5>(tab)<URL-2>(tab)<Optional MD5>

其中URL是测序结果文件的路径,可以是ftp站点上的文件;MD5是可选的,以利于下载文件后
用于文件的完整度检验,若不想进行检验,则设置为0.

--output <path>
输出文件的路径。如果 --just-preprocess 参数给定,则只是对reads文件进行预
处理,给出结果。默认情况下,则是输出最终结果: 由SOAPsnp对每个参考序列进行SNP calls
的计算结果。

--intermediate <path> default: /tmp/myrna/intermediate/<PID>
中间结果文件存放的临时路径。如果指定了 --keep-intermediates 或 --keep-
all参数,则该文件则会一直存在。默认情况下会使用完毕后删除。

--preprocess-output <path>
对reads进行预处理的输出路径。由于对reads的预处理会消耗很多时间,因此保留数据
的预处理结果,对再次运行程序来说,可以节约很多时间。

--keep-intermediates <path>
保留程序运行中生成的结果文件夹和文件,默认情况下这些文件会被尽快的自动删除。

--keep-all
保留所有的临时文件,默认情况下这些文件会被尽快的自动删除。

--cpus <int> default: 1
设定程序运行使用的CPU线程数。

--bowtie <path>
    设定bowtie的路径

--fastq-dump <path>
    sra toolkit中fastq-dump的路径

--soapsnp <path>
    SOAPsnp的路径

以上输入参数有个注意事项: 当 –input 参数的输入为一个manifest文件时,文件中指向某一个fastq或sra文件的 URL 必须是绝对路径;如果是相对路径则会报错。

在EMR、Hadoop和local上使用Crossbow共有的参数如下:

--quality { phred33 | phred64 | solexa64 }  default: phred33
    设置输入reads的碱基质量格式。phred33:输入的碱基质量等于ASCII码值加上33;最
低碱基质量是“#”. phred64:输入的碱基质量等于ASCII码值加上64;最低碱基质量是“B”.

--preprocess
    当 --input 的参数是一个文本文件(含有输入序列信息的清单)时,则必须加入该参数
;该参数指定出了文本文件中的序列信息来,Crossbow则进行这些reads的预处理,将生成的预处
理reads结果输出到 --preprocess-output 指定的文件夹;如果不指定--preprocess
-output参数,则输出到 --intermediate 指定的文件夹。

--just-preprocess
    和 --preprocess 参数类似,但不同的是进行reads预处理后,即停止pipeline的
运行,并将结果输出到 --output指定的文件夹。

--just-align
    不将Myrna的pipeline运行到底,只是运行到比对结束,并将结果存放到 --output 
参数指定的URL.

--resume-align
    使用此命令来运行上一个命令(比对完毕)过后的阶段。 此时 --input 的参数则是上
一个命令中 --output的URL

--bowtie-args "<args>"  default: -m 1
    设置bowtie的参数。点击查看Bowtie manual

--discard-reads <fraction>  default: 0.0
    丢弃一定比例输的reads来进行程序的运行。比如 0.5 则代表丢弃50%的reads。该参
数应用于所有输入的reads。对于程序的调试比较有用。

--discard-ref-bins <fraction> default:0.0
    在进行SNP calling之前,丢弃一定比例的参考序列,再进行SNP calling。有利于
debugging。

--discard-all <fraction>
    等同于将上两个参数同时设定到该值。

--soapsnp-args "<args>"  default: -2 -u -n -q
    设置SOAPsnp的参数。点击查看SOAPsnp manual

--soapsnp-hap-args "<args>"  default: -r 0.0001
    对单核体测序数据设置的参数

--soapsnp-dip-args "<args>"  default: -r 0.00005 -e 0.0001
    对双核体测序数据设置的参数

--haploids <chromosome-list>
    后接逗号分开的参考序列名。这些序列在运用SOAPsnp进行SNP calling的时候按
hapolid进行处理,其余的则按diploid进行处理。默认下所有的参考序列都是按diploid
处理

--all-haploids
    设定该参数,表明在使用SOAPsnp的时候所有的参考序列都按hapolid进行处理。

--partition-len <int> default: 1,000,000
    The bin size to use when binning alignments into partitions
prior to SNP calling. If load imbalance occurrs in the SNP 
calling step (some tasks taking far longer than others), try 
decreasing this.

--dry-run
    Just generate a script containing the commands needed to 
launch the job, but don't run it. The script's location will be 
printed so that you may run it later.

--test
    搜索Crossbow需要的工具(bowtie,SOAPsnp和fastq-dump),来检测Crossbow
能否正常运行

--tempdir `<path>` default: /tmp/Crossbow/invoke.scripts
    Local directory where temporary files (e.g. dynamically 
generated scripts) should be deposited.

3. 输入的reference jar文件配置

该jar文件使用unzip解压后,包括3个文件夹和1个文件,具体如下:

1. sequences 文件夹: 该文件夹中包含很多fasta文件。fasta文件的取名一定要按 chrX.fa 来作为文件名,其中 X 是从0开始的连续的数字;每个fasta文件中只含有一条序列,fasta文件的header命名必须为 X.
2. cmap.txt 文件:该文件记录着参考序列的真实名称和其对应着的数字。第一列为参考序列名,第二列为数字;中间使用tab分隔。
3. index 文件夹:使用bowtie-build将sequences中的所有的fasta文件构建一个索引,并且索引文件的前缀必须为 index;并且索引的顺序必须按数字的id进行排序。
4. snps 文件夹:如果有已知的snp结果,放置在此文件夹下;SNP 结果的文件名必须和 sequences 文件夹下的 FASTA 文件的前缀一致,并加上后缀 .snps. 比如 chrX.snps。SNP结果文件的格式为文本文件,包含由tab分隔的10列内容。

1. 参考序列名
2. 1-based 位置
3. SNP是否有频率信息(1 = yes, 0 = no)
4. SNP是否有实验验证(1 = yes, 0 = no)
5. SNP是否实际是indel(1 = yes, 0 = no)
6. A allele的频率
7. C allele的频率
8. T allele的频率
9. G allele的频率
10. SNP id

将上述文件和问夹就准备好后,即可打包成jar文件。

$ jar cf ref-XXX.jar sequences snps index cmap.txt

4. Crossbow运行过程和结果

Crossbow运行有4步,分别是:

1. 对reads进行预处理。通过$CrossbowHome/Copy.pl进行reads的预处理。
2. 进行序列的比对。通过$CrossbowHome/Align.pl调用bowtie来进行序列比对。
3. 进行SNPs Calling。通过$CrossbowHome/Soapsnp.pl调用soapsnp来Call SNPs。
4. 后处理。通过$CrossbowHome/CBFinish.pl来处理结果,使结果和原始的参考序列名一致。

Crossbow得到的结果在 –output 参数设置的文件夹的 crossbow_results 子文件夹中。每个参考序列得到一个使用gzip压缩后的结果。比如:
/crossbow_results/chr1.gz
/crossbow_results/chr2.gz
/crossbow_results/chr3.gz

每一个文件包含着SOAPsnp格式的SNPs,每行记录一个SNP。该SOAPsnp格式的包括18列,其意义如下:

1. 参考序列名
2. 1-based位置
3. 参考序列的基因型
4. 样品的基因型
5. 样品基因型的质量分数
6. 最优的碱基
7. 最优碱基的平均碱基质量
8. 唯一比对到该位点的对应着最优碱基的reads数
9. 所有比对到该位点的对应着最优碱基的reads数
10. 次优碱基
11. 次优碱基的平均碱基质量
12. 唯一比对到该位点的对应着次优碱基的reads数
13. 所有比对到该位点的对应着次优碱基的reads数
14. 该位点的测序深度
15. 该位点对应着paired比对的测序深度
16. Rank sum test P-value
17. Average copy number of nearby region
18. 该位点是否是一个已知的SNP(通过和 -s 参数中指定文件的对比)

注意: 第15列是Crossbow对SOAPsnp修饰后才添加的。

5. 点评:Crossbow和samtools进行SNPs calling

个人通过使用以上两种软件对实际数据的运算结果,发现:samtools比Crossbow得到的结果要好。表现为如下几点:

1. 这两个软件都是使用bowtie进行的比对;而Crossbow是使用SOAPsnp对比对结果进行SNP calling。而samtools结合其一起发布的bcftools进行SNPs calling。
2. 通过samtool tview对bowtie的比对结果进行浏览,来人工观察并比较这两个软件的给出的结果。发现samtools的结果比较正确。很多在samtools tview中明显显示为SNP的确在Crossbow中没有给出;但是,一些map质量很差,或比对上数目很少的SNPs却只在Crossbow中给出。
3. Crossbow只能进行SNPs calling,不能进行Indels calling。而samtools则能进行short Indels calling。
4. Crossbow使用bowtie作为比对软件;而samtools则可以运用bowtie2的比对结果。

综上,可以优先考虑使用samtools来进行SNP或Indel的calling。

发表评论

您的电子邮箱地址不会被公开。 必填项已用*标注

此站点使用Akismet来减少垃圾评论。了解我们如何处理您的评论数据