/
call.nf
149 lines (104 loc) · 3.53 KB
/
call.nf
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
#!/usr/bin/env nextflow
log.info ""
log.info '--------------------------------------------------'
log.info 'NEXTFLOW Variant Calling PAIRED END'
log.info '--------------------------------------------------'
if(params.help){
log.info ''
log.info 'Usage: '
log.info 'nextflow run --fastq_folder FASTQ/ --genome genome.fa --species species --genus genus'
log.info ''
log.info 'Mandatory arguments:'
log.info ' --fastq_folder FOLDER Folder containing FASTQ paired files'
log.info ' --genome FILE File containing the reference genome'
log.info ' --species STRING The species you are investigating e.g. chinensis'
log.info ' --genus String The genus (first name) of the species e.g. Actinidia'
log.info 'Options:'
log.info ''
exit 1
}
/*
* Files
*/
/*
* The reference genome file
* Config in nextflow.config:
*/
label_variant_calling_gatk = 'variant_calling_gatk'
genome_file = file(params.genome)
/*
* Creates the `read_pairs` channel that emits for each read-pair a tuple containing
* three elements: the pair ID, the first read-pair file and the second read-pair file
* This is the start of the pipeline and needs to be configured in nextflow.config (params.read).
*
* The raw read files
* Config in nextflow.config:
*/
Channel
.fromPath("${params.output_dir}/240.add_read_group_id/*.bam")
.set{ aligned_bam }
Channel
.fromPath("${params.output_dir}/240.add_read_group_id/*.bai")
.set{ aligned_bai }
label_variant_calling_gatk = 'variant_calling_gatk'
process index_faidx_unzip{
input:
file genome from genome_file
output:
file 'kiwitest.fasta' into genome_raw
"""
gunzip -c ${genome} > kiwitest.fasta
"""
}
process gatk_haplotype_calling{
module = params.variant_calling_gatk_module
publishDir "${params.publish_dir}/260.${tag}"
input:
file genome from genome_raw
file bwa_mdup_rg_bam from aligned_bam
file bwa_mdup_rg_bai from aligned_bai
output:
file "${bwa_mdup_rg_bam}.vcf" into gatk_gvcf
script:
"""
samtools faidx ${genome}
tt=${genome}
java -jar \$PICARD CreateSequenceDictionary \
R=${genome} O=\${tt%.fasta}.dict
java -Djava.io.tmpdir=/workspace/cfphxd/tmp -jar /software/bioinformatics/gatk-1.0/GenomeAnalysisTK.jar \
-T HaplotypeCaller \
-R ${genome} \
-I ${bwa_mdup_rg_bam} \
--genotyping_mode DISCOVERY \
-ploidy ${params.ploidy} \
-o ${bwa_mdup_rg_bam}.vcf
ln -s ${params.publish_dir}/260.${label_variant_calling_gatk} ${params.ouput_dir}
"""
}
gvcf_list = file("VCFlist.list")
process gatk_joint_calling{
module = params.joint_calling_module
input:
file genome from genome_file
file 'vcf_*' from gatk_gvcf.toList()
file gvcf_list from gvcf_list
output:
file 'gvcf_joint.g.vcf' into gvcf_joint
publishDir "${params.publish_dir}/270.joint_calling"
script:
"""
samtools faidx ${genome}
tt=${genome}
java -jar \$PICARD CreateSequenceDictionary \
R=${genome} O=\${tt%.fasta}.dict
cat ${gvcf_list} > list.list
ls vcf_* >> list.list
java -Djava.io.tmpdir=/workspace/cfphxd/tmp -jar /software/bioinformatics/gatk-1.0/GenomeAnalysisTK.jar \
-T GenotypeGVCFs \
-nt ${params.joint_calling_nt} \
-R ${genome} \
-V list.list \
-o gvcf_joint.g.vcf
ln -s ${params.publish_dir}/270.joint_calling ${params.ouput_dir}
"""
}