diff --git a/src/SNP_calling.config b/src/SNP_calling.config index 4b693893951517f76dd13ab6bc51448a8575c488..37cda85457cc041339d319eb49d7979f738f0e17 100644 --- a/src/SNP_calling.config +++ b/src/SNP_calling.config @@ -9,6 +9,9 @@ profiles { withName: trimming { container = "urqt:d62c1f8" } + withName: filter_fasta { + container = "bioawk:1.0" + } withName: index_fasta { container = "bwa:0.7.17" } @@ -21,9 +24,6 @@ profiles { withName: sort_bam { container = "sambamba:0.6.7" } - withName: name_fasta { - container = "samtools:1.7" - } withName: index_bam { container = "sambamba:0.6.7" } @@ -36,6 +36,18 @@ profiles { withName: HaplotypeCaller { container = "gatk:4.0.8.1" } + withName: GetPileupSummaries { + container = "gatk:4.0.8.1" + } + withName: CalculateContamination { + container = "gatk:4.0.8.1" + } + withName: CollectSequencingArtifactMetrics { + container = "gatk:4.0.8.1" + } + withName: filter_SNP { + container = "gatk:4.0.8.1" + } } } sge { diff --git a/src/SNP_calling.nf b/src/SNP_calling.nf index e52db5972c994bf7be7e5ad201b530b1c4f6eeca..eead4083d9466a0ea10d9efb3ebfa8d1b2921c9a 100644 --- a/src/SNP_calling.nf +++ b/src/SNP_calling.nf @@ -1,7 +1,9 @@ params.fastq = "$baseDir/data/*.fastq" params.fasta = "$baseDir/data/*.fasta" +params.seq_length = 800000 log.info "fastq files : ${params.fastq}" log.info "fasta files : ${params.fasta}" +log.info "fasta length to retain : ${params.seq_length}" def normal_sample = Eval.me(params.normal) def tumor_sample = Eval.me(params.tumor) log.info "normal : ${normal_sample}" @@ -11,11 +13,7 @@ Channel .fromPath( params.fasta ) .ifEmpty { error "Cannot find any fasta files matching: ${params.fasta}" } .map { it -> [(it.baseName =~ /([^\.]*)/)[0][1], it]} - .into { fasta_file; - indel_fasta_file; - recalibration_fasta_file; - haplotypecaller_fasta_file - } + .set { fasta_file } Channel .fromFilePairs( params.fastq ) .ifEmpty { error "Cannot find any fastq files matching: ${params.fastq}" } @@ -61,13 +59,39 @@ UrQt --t 20 --m ${task.cpus} --gz \ """ } +process filter_fasta { + tag "$fasta_id" + cpus 4 + publishDir "results/fasta/", mode: 'copy' + + input: + set fasta_id, file(fasta) from fasta_file + + output: + set fasta_idf, "*_filtered.fasta" into filter_fasta_files + + script: + fasta_idf = "${fasta_id}_filtered" +""" +bioawk -c fastx '{ if(length(\$seq) > $params.seq_length) { print ">"\$name; print \$seq }}' ${fasta} > \ +${fasta_id}_filtered.fasta +""" +} + +filter_fasta_files.into{ + filtered_fasta_files; + indel_fasta_file; + recalibration_fasta_file; + haplotypecaller_fasta_file +} + process index_fasta { tag "$fasta_id" cpus 4 publishDir "results/mapping/index/", mode: 'copy' input: - set fasta_id, file(fasta) from fasta_file + set fasta_id, file(fasta) from filtered_fasta_files output: set fasta_id, "${fasta.baseName}.*" into index_files @@ -205,7 +229,7 @@ index_bam_files.into{ } haplotypecaller_fasta_file.into{ - haplo_fasta_file; + final_fasta_file; index2_fasta_file index3_fasta_file } @@ -242,20 +266,49 @@ samtools faidx ${fasta} """ } -haplotypecaller_bam_files_norm = haplo_bam_files_norm +final_bam_files_norm = haplo_bam_files_norm .filter{ "normal_sample" == it[0] } -haplotypecaller_bam_files_tumor = haplo_bam_files_tumor +final_bam_files_tumor = haplo_bam_files_tumor .filter{ "tumor_sample" == it[0] } indexed_bam_files.into { index_bam_files_norm; index_bam_files_tumor } -indexed_bam_files_norm = index_bam_files_norm +final_indexed_bam_files_norm = index_bam_files_norm .filter{ "normal_sample" == it[0] } -indexed_bam_files_tumor = index_bam_files_tumor +final_indexed_bam_files_tumor = index_bam_files_tumor .filter{ "tumor_sample" == it[0] } +final_bam_files_norm.set{ + haplotypecaller_bam_files_norm +} +final_bam_files_tumor.into{ + haplotypecaller_bam_files_tumor; + artifact_bam_files_tumor; + pileup_bam_files_tumor +} +final_indexed_bam_files_norm.set{ + haplo_index_bam_files_norm +} +final_indexed_bam_files_tumor.into{ + haplo_index_bam_files_tumor; + artifact_index_bam_files_tumor; + pileup_index_bam_files_tumor +} +final_fasta_file.into{ + haplo_fasta_file + artifact_fasta_file +} +indexed2_fasta_file.into{ + haplo_indexed2_fasta_file + artifact_indexed2_fasta_file +} +indexed3_fasta_file.into{ + haplo_indexed3_fasta_file; + artifact_indexed3_fasta_file +} + process HaplotypeCaller { tag "$file_id_norm" cpus 10 @@ -263,21 +316,22 @@ process HaplotypeCaller { input: set file_id_norm, file(bam_norm) from haplotypecaller_bam_files_norm - set file_ididx_norm, file(bamidx_norm) from indexed_bam_files_norm + set file_ididx_norm, file(bamidx_norm) from haplo_index_bam_files_norm set file_id_tumor, file(bam_tumor) from haplotypecaller_bam_files_tumor - set file_ididx_tumor, file(bamidx_tumor) from indexed_bam_files_tumor + set file_ididx_tumor, file(bamidx_tumor) from haplo_index_bam_files_tumor set genome_id, file(fasta) from haplo_fasta_file - set genome2_idx, file(fasta2idx) from indexed2_fasta_file - set genome3_idx, file(fasta3idx) from indexed3_fasta_file + set genome2_idx, file(fasta2idx) from haplo_indexed2_fasta_file + set genome3_idx, file(fasta3idx) from haplo_indexed3_fasta_file output: set file_id_norm, "*.vcf" into vcf_files + set file_id_norm, "*.vcf.idx" into index_vcf_files set file_id_norm, "*.bam" into realigned_bams_files - set "*_mutect2_report.txt" into mutect2_report + file "*_mutect2_report.txt" into mutect2_report script: """ -gatk Mutect2 --native-pair-hmm-threads ${task.cpus} -R ${fasta} \ +gatk --java-options "-Xmx32G" Mutect2 --native-pair-hmm-threads ${task.cpus} -R ${fasta} \ -I ${bam_tumor} -tumor ${file_id_tumor} \ -I ${bam_norm} -normal ${file_id_norm} \ -O ${file_id_norm}_raw_calls.g.vcf \ @@ -285,75 +339,106 @@ gatk Mutect2 --native-pair-hmm-threads ${task.cpus} -R ${fasta} \ """ } +vcf_files.into{ + pileup_vcf_files; + filter_vcf_files +} +index_vcf_files.into{ + pileup_index_vcf_files; + filter_index_vcf_files +} + /* -process filter_SNP { - tag "$file_id" - cpus 4 +process GetPileupSummaries { + tag "$file_id_norm" + cpus 1 publishDir "results/SNP/vcf/", mode: 'copy' input: + set file_id_norm, file(vcf) from pileup_vcf_files + set fileidx_id_norm, file(vcfidx) from pileup_index_vcf_files + set file_id_tumor, file(bam_tumor) from pileup_bam_files_tumor output: - set file_id, "*.vcf" into vcf_files_filtered + set file_id_tumor, "*.table" into pileup_files + file "*_pileup_report.txt" into pileup_report script: """ -gatk --java-options "-Xmx2g" Mutect2 \ --R hg38/Homo_sapiens_assembly38.fasta \ --I tumor.bam \ --I normal.bam \ --tumor HCC1143_tumor \ --normal HCC1143_normal \ --pon resources/chr17_pon.vcf.gz \ ---germline-resource resources/chr17_af-only-gnomad_grch38.vcf.gz \ ---af-of-alleles-not-in-resource 0.0000025 \ ---disable-read-filter MateOnSameContigOrNoMappedMateReadFilter \ --L chr17plus.interval_list \ --O 1_somatic_m2.vcf.gz \ --bamout 2_tumor_normal_m2.bam - -gatk Mutect2 \ --R ~/Documents/ref/hg38/Homo_sapiens_assembly38.fasta \ --I HG00190.bam \ --tumor HG00190 \ ---disable-read-filter MateOnSameContigOrNoMappedMateReadFilter \ --L chr17plus.interval_list \ --O 3_HG00190.vcf.gz - -gatk CreateSomaticPanelOfNormals \ --vcfs 3_HG00190.vcf.gz \ --vcfs 4_NA19771.vcf.gz \ --vcfs 5_HG02759.vcf.gz \ --O 6_threesamplepon.vcf.gz - -gatk GetPileupSummaries \ --I tumor.bam \ --V resources/chr17_small_exac_common_3_grch38.vcf.gz \ --O 7_tumor_getpileupsummaries.table - -gatk CalculateContamination \ --I 7_tumor_getpileupsummaries.table \ --O 8_tumor_calculatecontamination.table +gatk --java-options "-Xmx32G" GetPileupSummaries \ +-I ${bam_tumor} \ +-V ${vcf} \ +-O ${file_id_tumor}_pileup.table \ +2> ${file_id_tumor}_pileup_report.txt +""" +} -gatk FilterMutectCalls \ --V somatic_m2.vcf.gz \ ---contamination-table tumor_calculatecontamination.table \ --O 9_somatic_oncefiltered.vcf.gz +process CalculateContamination { + tag "$file_id_tumor" + cpus 1 + publishDir "results/SNP/vcf/", mode: 'copy' + + input: + set file_id_tumor, file(pileup_table) from pileup_files + + output: + set file_id_tumor, "*.table" into contamination_files + file "*_contamination_report.txt" into contamination_report + + script: +""" +gatk --java-options "-Xmx32G" CalculateContamination \ +-I ${pileup_table} \ +-O $file_id_tumor}_contamination.table \ +2> ${file_id_tumor}_contamination_report.txt +""" +} +*/ + +process CollectSequencingArtifactMetrics { + tag "$file_id_tumor" + cpus 1 + publishDir "results/SNP/vcf/", mode: 'copy' + + input: + set file_id_tumor, file(bam_tumor) from artifact_bam_files_tumor + set genome_id, file(fasta) from artifact_fasta_file + set genome2_idx, file(fasta2idx) from artifact_indexed2_fasta_file + set genome3_idx, file(fasta3idx) from artifact_indexed3_fasta_file + + output: + set file_id_tumor, "*_artifact.*" into artifact_files + file "*_artifact_report.txt" into artifact_report + script: +""" gatk CollectSequencingArtifactMetrics \ --I tumor.bam \ --O 10_tumor_artifact \ -–-FILE_EXTENSION ".txt" \ --R ~/Documents/ref/hg38/Homo_sapiens_assembly38.fasta - -gatk FilterByOrientationBias \ --A G/T \ --A C/T \ --V 9_somatic_oncefiltered.vcf.gz \ --P tumor_artifact.pre_adapter_detail_metrics.txt \ --O 11_somatic_twicefiltered.vcf.gz +-I ${bam_tumor} \ +-O ${file_id_tumor}_artifact \ +-R ${fasta} \ +2> ${file_id_tumor}_artifact_report.txt +""" +} + +process filter_SNP { + tag "$file_id_norm" + cpus 1 + publishDir "results/SNP/vcf/", mode: 'copy' + input: + set file_id_norm, file(vcf) from filter_vcf_files + set fileidx_id_norm, file(vcfid) from filter_index_vcf_files + + output: + set file_id_norm, "*.vcf" into vcf_files_filtered + file "*_filter_report.txt" into filter_report + + script: +""" +gatk FilterMutectCalls \ +-V ${vcf} \ +-O ${file_id_norm}_filtered.vcf \ +2> ${file_id_norm}_filter_report.txt """ } -*/