illumina CASAVA 序列比对使用说明

我们假设任务是序列比对,而不包括Bcl文件的转换,也不包括后期的数据处理,只单单是把序列比对的程序跑起来。

(一)设置参考序列文件

在安装好CASAVA后,我们需要做的是准备好reference文件。它有一定的目录格式和文件格式要求,具体需要参见其说明手册。但是主要的一些参考序列文件在illumina的网站上可以下载到。可以下载到的文件其中包括:

  1. Arabidopsis_thaliana/
  2. Bacillus_cereus_ATCC_10987/
  3. Bos_taurus/
  4. Caenorhabditis_elegans/
  5. Canis_familiaris
  6. Drosophila_melanogaster/
  7. Equus_caballus/
  8. Escherichia_coli_K_12_DH10B/
  9. Escherichia_coli_K_12_MG1655/
  10. Gallus_gallus/
  11. Homo_sapiens/
  12. Mus_musculus/
  13. Mycobacterium_tuberculosis_H37RV/
  14. Oryza_sativa_japonica/2
  15. Pan_troglodytes/
  16. PhiX/
  17. Pseudomonas_aeruginosa_PAO1/
  18. Rattus_norvegicus/
  19. Rhodobacter_sphaeroides_2.4.1/
  20. Saccharomyces_cerevisiae/
  21. Sorghum_bicolor/
  22. Staphylococcus_aureus_NCTC_8325/
  23. Sus_scrofa/
  24. Zea_mays/

使用ftp下载,地址:ftp.illumina.com

用户名: igenome

密码:G3nom3s4u

以人类基因组为例,我们使用命令:

wget –ftp-user=igenome –ftp-password=G3nom3s4u ftp://ftp.illumina.com/Homo_sapiens/UCSC/hg18/Homo_sapiens_UCSC_hg18.tar.gz

即可。文件4.28GB下载需要一段时间。

(二)设置配置文件

接下来的工作就是设置好配置文件。之所以使用配置文件,一般都是因为参数太多,如果都写在命令行的话,命令行显得冗长。

关于输入文件。对于CASAVA 1.8来说,输入文件有了较大的变化:

  1. 使用FASTQ文件作为输入源
  2. 可以通过读取样品描述文件SampleSheet.csv来识为输入源。

输入源期待有以下内容:

  1. Unaligned目录,其目录下以Project名为子目录,子目录当中存放需要比对的fastq.gz文件。注意,这些文件是gzip压缩的。
  2. 配置文件config.txt。
  3. 多路复用总结文件summary.xml
  4. 多路复用配制文件config.xml
  5. FASTA格式的参考序列文件
  6. 如果是RNA比对,还需要一些相关文件, refFlat.txt.gz 或者 seq_gene.md.gz。这些文件可以在UCSC或者NCBI上下载到。还有非基因组序列(contaminat sequences),比如线粒体基因组等。

对于config.txt的内容,可以参考下面的文件,但使用前最好把注释的所有内容都去除。

######################################################################
#####核心参数
## Unaligned所在路径
EXPT_DIR $HOME/data/expt_name/Unaligned
 
## 短序列处理方式
## 几个例子:
## USE_BASES nYYY	#忽略第一个碱基,使用每2至4个碱基
## USE_BASES Y30	#比对前30个碱基
## USE_BASES nY30	#忽略第一个碱基,比对基后的30个碱基
## USE_BASES nY30n	#忽略第一个碱基,比对基后的30个碱基,并且忽略最后一个碱基
## USE_BASES nY*n	#忽略第一个碱基,和最后一个碱基,碱基数设定为短序数-2
## USE_BASES Y*n	#忽略最后一个碱基,碱基数设定为短序数-1
## USE_BASES Y*n,Y*n	#成对比对,并且均忽略最后一个碱基,且成对的两短序不一定等长
## USE_BASES nY*,nY*	#成对比对,并且均忽略第一个碱基,且成对的两短序不一定等长
## USE_BASES nY*	#忽略第一个碱基,碱基数设定为短序数-1
## USE_BASES n*,Y*n	#成对读入,但是第一组全部忽略,第二组忽略最后一个碱基
## USE_BASES Y*n,n*	#成对读入,但是第一组忽略最后一个碱基,第二组全部忽略
USE_BASES Y*	#使用全部

## 参考序列文件所在目录,如果想使用ELANDv2e的话,就指向中一个单一的FASTA参考序列文件
ELAND_GENOME $HOME/iGenome/Homo_sapiens/UCSC/hg19/Sequence/Chromosomes/
SAMTOOLS_GENOME	# Direct CASAVA to a multi-sequence FASTA reference file.

## 指定比对方式,有四种,分别是
## ANALYSIS eland_extended	逐一比对
## ANALYSIS eland_pair		成对比对
## ANALYSIS eland_rna		RNA序列比对
## ANALYSIS none		未设定
ANALYSIS eland_extended
 
## 并行处理的文件数。为使CASAVA占用内存小于2GB,所以需要设定同时处理的数据大小。
## 其计算公式为:
## (ELAND_FASTQ_FILES_PER_PROCESS value) × (--fastq-cluster-count value) ≤ ## 16 million
## 要注意,该值的设定影响着内存的控制,如果太大就会占用更多的内存,
## 从而导致内存不足而影响效率。
## 以下是参考数值:
## --fastq-cluster-count	|	ELAND_FASTQ_FILES_PER_PROCESS
## 12,000,000			|	1
##  6,000,000			|	2
##  4,000,000			|	3
##  3,000,000			|	4
##  2,000,000			|	6
##  1,000,000			|	12
## 关于如何获取--fastq-cluster-cout,这其实是Bcl文件转换过程当中的一个参数,
## 可以简单的想象成fastq文件行数除以4
ELAND_FASTQ_FILES_PER_PROCESS 3
 
######################################################################
#####可选参数,使用前去除注释答符号
## --singleseed模式,只针对ANALYSIS eland_extended和ANALYSIS eland_pair有效
#SINGLESEED --singleseed

## --ungapped模式,不允许空格
#UNGAPPED --ungapped

## --sensitive增强灵敏度,以速度为代价
#INCREASED_SENSITIVITY --sensitive

## 输出文件夹,默认为当前目录
#OUT_DIR $HOME/output/

## 后续命令行。
#DATASET_POST_RUN_COMMAND /yourPath/yourCommand yourArgs

## 程序运行完全结束之后,发邮件通知
#EMAIL_LIST user@example.com user2@example.com
#EMAIL_SERVER mailserver
#EMAIL_DOMAIN example.com

## 在超链接中运行目录前使用web前缀,
## NUM_LEADING_DIRS_TO_STRIP,在使用前缀前运行目录需要截取的字符数
#WEB_DIR_ROOT file://server.example.com/share
#NUM_LEADING_DIRS_TO_STRIP 0

## 非基因组序列文件所在目录
#ELAND_RNA_GENOME_CONTAM $HOME/iGenome/Homo_sapiens/UCSC/hg19/Sequence/Contaminants/

## 基因注释文件
#ELAND_RNA_GENOME_ANNOTATION $HOME/iGenome/Homo_sapiens/UCSC/hg19/Annotation/Genes/refFlat.txt.gz

##不解释
#ELAND_RNA_GENE_MD_GROUP_LABEL GRCh37.p2-Primary Assembly

##不解释
#KAGU_PARAMS

######################################################################
#####paired模式参数,不解释,省略。

如果有多个不同的project的话,可以使用SampleSheet.csv来细化分析。对于SampleSheet.csv的内容,可以参考

FCID Lane SampleID SampleRef Index Description Control Recipe Operator Project
12345AAXX 1 sample1 human ATCACG desc1 N R1 name Proj1
12345AAXX 1 sample2 human CGATGT desc2 N R1 name Proj1
12345AAXX 2 sample3 rat TTAGGC desc3 N R1 name Proj2
12345AAXX 2 sample4 mouse TGACCA desc4 N R1 name Proj3

在使用SampleSheet.csv之后,config.txt也需要做出相应的变化。如果针对某个具体的Project,reference,sample,barcode,lane的话,需要使用前缀来特别指定。在读取config.txt时,CASAVA会依照PROJECT,REFERENCE,SAMPLE,BARCODE,Lane,Global settings的优先顺序来进行参数设置。如果全局设置的ANALYSIS为eland_pair的话,你如果设置PROJECT proj1 ANALYSIS eland_rna的话,对于proj1就会使用eland_rna来运算。具体示例如下:

PROJECT projx REFERENCE SampleRef_x SAMPLE SampleID_x 2: ANALYSIS eland_pair

(三)运行

在config.txt所在的目录试运行命令(不加–make)

/path-to-CASAVA/bin/configureAlignment.pl config.txt –EXPT_DIR path_to_Unaligned_folder

在config.txt所在的目录运行命令

/path-to-CASAVA/bin/configureAlignment.pl config.txt –EXPT_DIR path_to_Unaligned_folder –make

make all -j n

其中-j n代表使用多个CPU,有几个CPU,n就是几。

Leave a Reply

  

  

  

%d 博主赞过: