我们假设任务是序列比对,而不包括Bcl文件的转换,也不包括后期的数据处理,只单单是把序列比对的程序跑起来。
(一)设置参考序列文件
在安装好CASAVA后,我们需要做的是准备好reference文件。它有一定的目录格式和文件格式要求,具体需要参见其说明手册。但是主要的一些参考序列文件在illumina的网站上可以下载到。可以下载到的文件其中包括:
- Arabidopsis_thaliana/
- Bacillus_cereus_ATCC_10987/
- Bos_taurus/
- Caenorhabditis_elegans/
- Canis_familiaris
- Drosophila_melanogaster/
- Equus_caballus/
- Escherichia_coli_K_12_DH10B/
- Escherichia_coli_K_12_MG1655/
- Gallus_gallus/
- Homo_sapiens/
- Mus_musculus/
- Mycobacterium_tuberculosis_H37RV/
- Oryza_sativa_japonica/2
- Pan_troglodytes/
- PhiX/
- Pseudomonas_aeruginosa_PAO1/
- Rattus_norvegicus/
- Rhodobacter_sphaeroides_2.4.1/
- Saccharomyces_cerevisiae/
- Sorghum_bicolor/
- Staphylococcus_aureus_NCTC_8325/
- Sus_scrofa/
- 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来说,输入文件有了较大的变化:
- 使用FASTQ文件作为输入源
- 可以通过读取样品描述文件SampleSheet.csv来识为输入源。
输入源期待有以下内容:
- Unaligned目录,其目录下以Project名为子目录,子目录当中存放需要比对的fastq.gz文件。注意,这些文件是gzip压缩的。
- 配置文件config.txt。
- 多路复用总结文件summary.xml
- 多路复用配制文件config.xml
- FASTA格式的参考序列文件
- 如果是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就是几。