生信流程软件-nextflow¶

//结果输出目录
params.output = "/public/home/software/test/callsnp/output"

//参考基因组
params.genome="/public/home/software/test/callsnp/data/Gbarbadense_genome.fasta"

//reads路径
params.read_file="/public/home/software/test/callsnp/data/*_{1,2}.fastq.gz"

//将参考基因组转换成file对象，以便使用getParent、getBaseName方法获取参考基因组所在路径和去掉.fa .fasta后的basename
ref_file=file(params.genome)
ref_dir=ref_file.getParent()
ref_base=ref_file.getBaseName()

//fromFilePairs默认读入一对reads，可以通过size参数设置
Channel.fromFilePairs(params.read_file)
    .ifEmpty{error "Cannot find  reads $params.read_file"}
    .into {reads}
//判断是否已经存在bwa索引文件，存在则将所有索引文件放入index_exi 这个channel；
//否则运行Run_bwa_index这process来建索引，然后输出到名字为index_exi的channel
Run_bwa_index
if(file(params.genome+".amb").exists()){
    index_exi = Channel.fromPath("${params.genome}.*")
}else{
process Run_bwa_index {

    //申请一个CPU核心
    cpus 1

    //使用PBS作业调度系统，使用local可在本地跑，也可以支持sge、slurm、lsf等调度系统
    executor = 'pbs'    

    //将output文件复制到ref_dir这个目录，mode还可以选move等
    publishDir "ref_dir",mode:'copy'

    //输入文件
    input:
    file genome from ref_file

    //输出
    output:
    file "${genome}.*" into index_exi

    //该process运行的条件，bwa index不存在时运行，其实这个地方有点重复
    when:
    !file(params.genome+".amb").exists()

    //两个"""包裹的部分是用户自己的各种脚本，一般用户需要跑的程序都写在这里，默认是shell脚本，还可以是perl、python、R等
    """
    module purge
    module load BWA
    bwa index $genome
    """
}
}

//给后面GATK运行创建genome index
process Run_dict_fai {

    cpus 1
    executor = 'pbs' 

    //pbs申请的memory，可以直接定义数字，也可以如下动态定义
    //如果这个process运行失败，可以重新开始，每次重新开始申请的内存增加2GB
    memory { 2.GB * task.attempt }

    //process运行出错时的行为，默认直接退出。因为java作业容易因为内存的原因而运行出错，
    //因此这里设置为retry，失败重新运行
    //还可以有ignore、finish、terminate
    errorStrategy 'retry'

    //process运行失败后，重新尝试的次数
    maxRetries 3

    //这个process的output的内容copy到ref_dir目录
    publishDir "ref_dir",mode:'copy'


    input:
    file genome from ref_file

    //将定义的这些文件作为一个整体输出到GATK_index这个channel
    output:
    set file("${genome}"),file("${genome}.fai"),file("${ref_base}.dict") into GATK_index

    //nextflow脚本，这里根据上面memory的值定义了一个内存值给下面的picard使用
    script:
    picardmem = "${task.memory}".replaceAll(/ GB/,"g")

    """
    module purge
    module load picard
    java -Xmx${picardmem} -jar \${EBROOTPICARD}/picard.jar CreateSequenceDictionary R=${genome} O=${ref_base}.dict
    module purge
    module load SAMtools/1.8
    samtools faidx ${genome} > ${genome}.fai
    """
}

//运行BWA比对
process Run_bwa {

    cpus 6
    executor = 'pbs'    

    input:
    //从reads中获取read id和两条read序列
    set val(id),file(readfile) from reads

    //输入bwa index文件，collect方法将这个channel中的所有item集合为一个list
    //否则就会一个个文件循环emit，不符合我们的要求
    file index from index_exi.collect()

    //将read id和生成的Sam文件放入channel bwa_sam中，给后面的process使用
    output:
    set val("${id}"),file("${id}.sam") into bwa_sam

    //process运行条件，true为运行，false为不运行，在此可以不需要这个when条件，仅用于测试
    when:
    true

    //通过输入的bwa_index文件提取bwa build生成的index文件的prefix，注意这里是Gbarbadense_genome.fasta，而不是Gbarbadense_genome
    //当然可以通过在前面的bwa build 加-p 指定生成的index文件的prefix为前面提取的ref_base,那么久可以不用这一步，下面bwa mem的index_base改为ref_base即可
    script:
    index_base = index[0].toString() - ~/.amb/

    """
    module purge
    module load BWA
    bwa mem -M -t ${task.cpus} -R '@RG\\tID:$id\\tLB:HISEQ10X\\tSM:$id' ${index_base} $readfile > ${id}.sam
    """
}

//运行picar，给sam文件排序、去重
process Run_picard {

    executor = 'pbs'    
    cpus 1

    //失败重试，动态资源申请
    memory { 2.GB * task.attempt }
    errorStrategy 'retry'
    maxRetries 3

    //这2个结果后面的process不会用到，但结果可能有用，可以move到自己指定的目录
    publishDir "$params.output/metrics",pattern:"*_metrics.txt",mode:'move'

    //删除不需要的文件
    afterScript     "rm ${samfile};rm ${id}_srt.bam"

    input:

    //Run_bwa生成的output channel作为这一步的input channel
    set val(id),file(samfile) from bwa_sam

    //输入参考基因组
    file genome from ref_file

    output:

    //输出排序后的bam
    file "${id}_srt.bam" into srtbam

    //输出metrics文件
    file "${id}_metrics.txt" into metrics

    //将id和相应的排序去重后的bam文件输入到2个channel中，因为一个channel只能使用一次，后面有2个process需要用到这个output channel的内容
    //所以需要输出2个内容相同、名称不同的channel
    set val("${id}"), file("${id}_srt_redup.bam") into samtools_input,GATK_input

    //所用同上，调试之用
    when:
    true

    //作用同上，用于给java程序设定内存参数
    script:
    picardmem = "${task.memory}".replaceAll(/ GB/,"g")

    """
    module purge
    module load picard
    java -Xmx${picardmem} -jar \${EBROOTPICARD}/picard.jar SortSam I=${samfile} O=${id}_srt.bam SORT_ORDER=coordinate
    java -Xmx${picardmem} -jar \${EBROOTPICARD}/picard.jar MarkDuplicates INPUT=${id}_srt.bam OUTPUT=${id}_srt_redup.bam CREATE_INDEX=true VALIDATION_STRINGENCY=SILENT METRICS_FILE=${id}_metrics.txt
    """
}

//samtools callsnp
process Run_samtools {

    cpus 1
    executor = 'pbs'    

    //output channel中的文件move到指定目录
    publishDir "$params.output/samtools",mode:'move'

    input:
    set val(id),file(bam) from samtools_input
    file genome from ref_file

    output:
    file "${id}.samtools.vcf" into samtoolsvcf

    when:
    true

    """
    module purge
    module load SAMtools/1.8
    module load BCFtools/1.6
    samtools mpileup -t DP -t SP -ugf ${genome} ${bam} |bcftools call --threads ${task.cpus} -mo ${id}.samtools.vcf
    """
}

process Run_GATK {

    cpus 4
    executor = 'pbs'    
    memory { 20.GB * task.attempt }
    errorStrategy 'retry'
    maxRetries 3

    publishDir "$params.output/GATK",mode:'move'

    input:

    //输入bam
    set val(id),file(bam) from GATK_input

    //输入index
    set file(genome),file(fai),file(dict) from GATK_index

    //将所有生成文件输出到GATKvcf channel中
    output:
    file "*" into GATKvcf

    when:
    true

    script:
    gatkmem = "${task.memory}".replaceAll(/ GB/,"g")

    """
    module purge
    module load GATK/3.8-0-Java-1.8.0_92
    java -Xmx${gatkmem} -jar \$EBROOTGATK/GenomeAnalysisTK.jar -R ${genome} -T HaplotypeCaller --emitRefConfidence GVCF -I ${bam} -o ${id}.gatk.vcf --variant_index_type LINEAR --variant_index_parameter 128000 -nct ${task.cpus}
    """
}