Skip to content

Commit 59228bb

Browse files
jbv2jbv2
committed
getting manipulate damage from eager/dev
1 parent 5f5c10a commit 59228bb

1 file changed

Lines changed: 127 additions & 121 deletions

File tree

‎subworkflows/local/manipulate_damage.nf‎

Lines changed: 127 additions & 121 deletions
Original file line numberDiff line numberDiff line change
@@ -2,171 +2,177 @@
22
// Calculate PMD scores, trim, or rescale DNA damage from mapped reads.
33
//
44

5-
include { addNewMetaFromAttributes } from '../../subworkflows/local/utils_nfcore_eager_pipeline/main'
5+
include { addNewMetaFromAttributes } from '../../subworkflows/local/utils_nfcore_eager_pipeline/main'
66

7-
include { BEDTOOLS_MASKFASTA } from '../../modules/nf-core/bedtools/maskfasta/main'
8-
include { MAPDAMAGE2 } from '../../modules/nf-core/mapdamage2/main'
9-
include { PMDTOOLS_FILTER } from '../../modules/nf-core/pmdtools/filter/main'
10-
include { BAMUTIL_TRIMBAM } from '../../modules/nf-core/bamutil/trimbam/main'
11-
include { SAMTOOLS_INDEX as SAMTOOLS_INDEX_DAMAGE_RESCALED } from '../../modules/nf-core/samtools/index/main'
12-
include { SAMTOOLS_INDEX as SAMTOOLS_INDEX_DAMAGE_FILTERED } from '../../modules/nf-core/samtools/index/main'
13-
include { SAMTOOLS_INDEX as SAMTOOLS_INDEX_DAMAGE_TRIMMED } from '../../modules/nf-core/samtools/index/main'
14-
include { SAMTOOLS_FLAGSTAT as SAMTOOLS_FLAGSTAT_DAMAGE_FILTERED } from '../../modules/nf-core/samtools/flagstat/main'
7+
include { BEDTOOLS_MASKFASTA } from '../../modules/nf-core/bedtools/maskfasta/main'
8+
include { MAPDAMAGE2 } from '../../modules/nf-core/mapdamage2/main'
9+
include { PMDTOOLS_FILTER } from '../../modules/nf-core/pmdtools/filter/main'
10+
include { BAMUTIL_TRIMBAM } from '../../modules/nf-core/bamutil/trimbam/main'
11+
include { SAMTOOLS_INDEX as SAMTOOLS_INDEX_DAMAGE_RESCALED } from '../../modules/nf-core/samtools/index/main'
12+
include { SAMTOOLS_INDEX as SAMTOOLS_INDEX_DAMAGE_FILTERED } from '../../modules/nf-core/samtools/index/main'
13+
include { SAMTOOLS_INDEX as SAMTOOLS_INDEX_DAMAGE_TRIMMED } from '../../modules/nf-core/samtools/index/main'
14+
include { SAMTOOLS_FLAGSTAT as SAMTOOLS_FLAGSTAT_DAMAGE_FILTERED } from '../../modules/nf-core/samtools/flagstat/main'
1515

1616
// TODO: Add required channels and channel manipulations for reference-dependent bed masking before pmdtools. Requires multi-ref support before implementation.
1717
workflow MANIPULATE_DAMAGE {
1818
take:
19-
ch_bam_bai // [ [ meta ], bam , bai ]
20-
ch_fasta // [ [ meta ], fasta ]
21-
ch_fasta_elongated // [ [ meta ], fasta ]
22-
ch_pmd_masking // [ [ meta ], masked_fasta, bed_for_masking ]
19+
ch_bam_bai // [ [ meta ], bam , bai ]
20+
ch_fasta // [ [ meta ], fasta ]
21+
ch_pmd_masking // [ [ meta ], masked_fasta, bed_for_masking ]
2322

2423
main:
25-
ch_versions = Channel.empty()
26-
ch_rescaled_bams = Channel.empty()
27-
ch_pmd_filtered_bams = Channel.empty()
28-
ch_trimmed_bams = Channel.empty()
29-
ch_pmd_filtered_flagstat = Channel.empty()
30-
// Only run flagstat on pmd filtered bam, since rescaling and trimming does not change the number of reads
24+
ch_versions = Channel.empty()
25+
ch_rescaled_bams = Channel.empty()
26+
ch_pmd_filtered_bams = Channel.empty()
27+
ch_trimmed_bams = Channel.empty()
28+
ch_pmd_filtered_flagstat = Channel.empty() // Only run flagstat on pmd filtered bam, since rescaling and trimming does not change the number of reads
3129

3230
// Ensure correct reference is associated with each bam_bai pair
33-
if (params.mapping_tool == 'circularmapper') {
34-
ch_refs = ch_fasta_elongated.map {
35-
// Prepend a new meta that contains the meta.id value as the new_meta.reference attribute
36-
addNewMetaFromAttributes(it, "id", "reference", false)
37-
}
38-
}
39-
else {
40-
ch_refs = ch_fasta.map {
31+
ch_refs = ch_fasta
32+
.map {
4133
// Prepend a new meta that contains the meta.id value as the new_meta.reference attribute
42-
addNewMetaFromAttributes(it, "id", "reference", false)
34+
addNewMetaFromAttributes( it, "id" , "reference" , false )
4335
}
44-
}
45-
4636

4737
ch_input_for_damage_manipulation = ch_bam_bai
4838
.map {
4939
// Prepend a new meta that contains the meta.reference value as the new_meta.reference attribute
50-
addNewMetaFromAttributes(it, "reference", "reference", false)
51-
}
52-
.combine(ch_refs, by: 0)
53-
// [ [combine_meta], [meta], bam, bai, [ref_meta], fasta ]
54-
55-
if (params.run_mapdamage_rescaling) {
56-
ch_mapdamage_prep = ch_input_for_damage_manipulation.branch { ignore_me, meta, bam, bai, ref_meta, fasta ->
57-
skip: meta.damage_treatment == 'full'
58-
no_skip: true
40+
addNewMetaFromAttributes( it, "reference" , "reference" , false )
5941
}
42+
.combine( ch_refs, by: 0 ) // [ [combine_meta], [meta], bam, bai, [ref_meta], fasta ]
43+
44+
if ( params.run_mapdamage_rescaling ) {
45+
ch_mapdamage_prep = ch_input_for_damage_manipulation
46+
.branch {
47+
ignore_me, meta, bam, bai, ref_meta, fasta ->
48+
skip: meta.damage_treatment == 'full'
49+
no_skip: true
50+
}
6051

61-
ch_skip_rescale = ch_mapdamage_prep.skip.map { ignore_me, meta, bam, bai, ref_meta, fasta ->
62-
[meta, bam, bai]
63-
}
52+
ch_skip_rescale = ch_mapdamage_prep.skip
53+
.map {
54+
ignore_me, meta, bam, bai, ref_meta, fasta ->
55+
[ meta, bam, bai ]
56+
}
6457

65-
ch_mapdamage_input = ch_mapdamage_prep.no_skip.multiMap { ignore_me, meta, bam, bai, ref_meta, fasta ->
66-
bam: [meta, bam]
67-
fasta: fasta
68-
}
58+
ch_mapdamage_input = ch_mapdamage_prep.no_skip
59+
.multiMap {
60+
ignore_me, meta, bam, bai, ref_meta, fasta ->
61+
bam: [ meta, bam ]
62+
fasta: fasta
63+
}
6964

70-
MAPDAMAGE2(ch_mapdamage_input.bam, ch_mapdamage_input.fasta)
71-
ch_versions = ch_versions.mix(MAPDAMAGE2.out.versions.first())
65+
MAPDAMAGE2( ch_mapdamage_input.bam, ch_mapdamage_input.fasta )
66+
ch_versions = ch_versions.mix( MAPDAMAGE2.out.versions.first() )
7267

73-
SAMTOOLS_INDEX_DAMAGE_RESCALED(MAPDAMAGE2.out.rescaled)
74-
ch_versions = ch_versions.mix(SAMTOOLS_INDEX_DAMAGE_RESCALED.out.versions.first())
68+
SAMTOOLS_INDEX_DAMAGE_RESCALED( MAPDAMAGE2.out.rescaled )
69+
ch_versions = ch_versions.mix( SAMTOOLS_INDEX_DAMAGE_RESCALED.out.versions.first() )
7570
ch_rescaled_index = params.fasta_largeref ? SAMTOOLS_INDEX_DAMAGE_RESCALED.out.csi : SAMTOOLS_INDEX_DAMAGE_RESCALED.out.bai
7671

7772
// TODO When adding library-level data merging pre-genotyping, make sure that rescaled bams are not merged in any way as the underlying damage model could differ between libraries
78-
ch_rescaled_bams = MAPDAMAGE2.out.rescaled
79-
.join(ch_rescaled_index)
80-
.mix(ch_skip_rescale)
73+
ch_rescaled_bams = MAPDAMAGE2.out.rescaled.join(ch_rescaled_index)
74+
.mix(ch_skip_rescale) // Should these be mixed actually, or excluded? Might not make sense to take rescaled and non-rescaled bams togetehr for anything downstream.
8175
}
8276

83-
if (params.run_pmd_filtering) {
77+
if ( params.run_pmd_filtering ) {
8478
ch_masking_prep = ch_pmd_masking
85-
.combine(ch_fasta, by: 0)
86-
.branch { meta, masked_fasta, bed, fasta ->
87-
alreadymasked: masked_fasta != ""
88-
tobemasked: masked_fasta == "" && bed != ""
89-
nomasking: masked_fasta == "" && bed == ""
90-
}
91-
92-
ch_masking_input = ch_masking_prep.tobemasked.multiMap { meta, masked_fasta, bed, fasta ->
93-
bed: [meta, bed]
94-
fasta: fasta
95-
}
96-
97-
BEDTOOLS_MASKFASTA(ch_masking_input.bed, ch_masking_input.fasta)
79+
.combine( ch_fasta, by: 0 ) // [ [meta], masked_fasta, bed, fasta ]
80+
.branch {
81+
meta, masked_fasta, bed, fasta ->
82+
alreadymasked: masked_fasta != ""
83+
tobemasked: masked_fasta == "" && bed != ""
84+
nomasking: masked_fasta == "" && bed == ""
85+
}
86+
87+
ch_masking_input = ch_masking_prep.tobemasked
88+
.multiMap{
89+
meta, masked_fasta, bed, fasta ->
90+
bed: [ meta, bed ]
91+
fasta: fasta
92+
}
93+
94+
BEDTOOLS_MASKFASTA( ch_masking_input.bed, ch_masking_input.fasta )
9895
ch_masking_output = BEDTOOLS_MASKFASTA.out.fasta
99-
ch_versions = ch_versions.mix(BEDTOOLS_MASKFASTA.out.versions.first())
100-
101-
ch_already_masked = ch_masking_prep.alreadymasked.map { meta, masked_fasta, bed, fasta ->
102-
[meta, masked_fasta]
103-
}
104-
105-
ch_no_masking = ch_masking_prep.nomasking.map { meta, masked_fasta, bed, fasta ->
106-
[meta, fasta]
107-
}
108-
109-
ch_pmd_fastas = ch_masking_output
110-
.mix(ch_already_masked, ch_no_masking)
111-
.map {
112-
// Prepend a new meta that contains the meta.id value as the new_meta.reference attribute
113-
addNewMetaFromAttributes(it, "id", "reference", false)
114-
}
96+
ch_versions = ch_versions.mix( BEDTOOLS_MASKFASTA.out.versions.first() )
97+
98+
ch_already_masked = ch_masking_prep.alreadymasked
99+
.map {
100+
meta, masked_fasta, bed, fasta ->
101+
[ meta, masked_fasta ]
102+
}
103+
104+
ch_no_masking = ch_masking_prep.nomasking
105+
.map {
106+
meta, masked_fasta, bed, fasta ->
107+
[ meta, fasta ]
108+
}
109+
110+
ch_pmd_fastas = ch_masking_output.mix( ch_already_masked, ch_no_masking )
111+
.map {
112+
// Prepend a new meta that contains the meta.id value as the new_meta.reference attribute
113+
addNewMetaFromAttributes( it, "id" , "reference" , false )
114+
}
115115

116116
ch_pmdtools_input = ch_bam_bai
117-
.map { addNewMetaFromAttributes(it, "reference", "reference", false) }
118-
.combine(ch_pmd_fastas, by: 0)
119-
.multiMap { combine_meta, meta, bam, bai, ref_meta, fasta ->
120-
bam: [meta, bam, bai]
121-
fasta: fasta
122-
}
123-
124-
PMDTOOLS_FILTER(ch_pmdtools_input.bam, params.damage_manipulation_pmdtools_threshold, ch_pmdtools_input.fasta)
125-
ch_versions = ch_versions.mix(PMDTOOLS_FILTER.out.versions.first())
126-
127-
SAMTOOLS_INDEX_DAMAGE_FILTERED(PMDTOOLS_FILTER.out.bam)
128-
ch_versions = ch_versions.mix(SAMTOOLS_INDEX_DAMAGE_FILTERED.out.versions.first())
117+
.map { addNewMetaFromAttributes( it, "reference" , "reference" , false ) }
118+
.combine( ch_pmd_fastas, by: 0 ) // [ [combine_meta], [meta], bam, bai, [ref_meta] fasta ]
119+
.multiMap {
120+
combine_meta, meta, bam, bai, ref_meta, fasta ->
121+
bam: [ meta, bam, bai ]
122+
fasta: fasta
123+
}
124+
125+
PMDTOOLS_FILTER( ch_pmdtools_input.bam, params.damage_manipulation_pmdtools_threshold, ch_pmdtools_input.fasta )
126+
ch_versions = ch_versions.mix( PMDTOOLS_FILTER.out.versions.first() )
127+
128+
SAMTOOLS_INDEX_DAMAGE_FILTERED( PMDTOOLS_FILTER.out.bam )
129+
ch_versions = ch_versions.mix( SAMTOOLS_INDEX_DAMAGE_FILTERED.out.versions.first() )
129130
ch_filtered_index = params.fasta_largeref ? SAMTOOLS_INDEX_DAMAGE_FILTERED.out.csi : SAMTOOLS_INDEX_DAMAGE_FILTERED.out.bai
130131

131-
ch_pmd_filtered_bams = PMDTOOLS_FILTER.out.bam.join(ch_filtered_index)
132+
ch_pmd_filtered_bams = PMDTOOLS_FILTER.out.bam.join( ch_filtered_index )
132133

133-
SAMTOOLS_FLAGSTAT_DAMAGE_FILTERED(ch_pmd_filtered_bams)
134+
SAMTOOLS_FLAGSTAT_DAMAGE_FILTERED( ch_pmd_filtered_bams )
134135
ch_pmd_filtered_flagstat = SAMTOOLS_FLAGSTAT_DAMAGE_FILTERED.out.flagstat
135-
ch_versions = ch_versions.mix(SAMTOOLS_FLAGSTAT_DAMAGE_FILTERED.out.versions.first())
136+
ch_versions = ch_versions.mix( SAMTOOLS_FLAGSTAT_DAMAGE_FILTERED.out.versions.first() )
136137
}
137138

138-
if (params.run_trim_bam) {
139-
if (params.run_pmd_filtering) {
140-
ch_to_trim = ch_pmd_filtered_bams.map { meta, bam, bai ->
141-
[meta, bam]
142-
}
143-
}
144-
else {
145-
ch_to_trim = ch_bam_bai.map { meta, bam, bai ->
146-
[meta, bam]
147-
}
139+
if ( params.run_trim_bam ) {
140+
if ( params.run_pmd_filtering ) {
141+
ch_to_trim = ch_pmd_filtered_bams
142+
.map{
143+
meta, bam, bai ->
144+
[ meta, bam ]
145+
}
146+
} else {
147+
ch_to_trim = ch_bam_bai
148+
.map {
149+
meta, bam, bai ->
150+
[ meta, bam ]
151+
}
148152
}
149153

150-
ch_trimbam_input = ch_to_trim.map { meta, bam ->
151-
def trim_left = meta.strandedness == 'single' ? (meta.damage_treatment == 'none' ? params.damage_manipulation_bamutils_trim_single_stranded_none_udg_left : meta.damage_treatment == 'half' ? params.damage_manipulation_bamutils_trim_single_stranded_half_udg_left : 0) : (meta.damage_treatment == 'none' ? params.damage_manipulation_bamutils_trim_double_stranded_none_udg_left : meta.damage_treatment == 'half' ? params.damage_manipulation_bamutils_trim_double_stranded_half_udg_left : 0)
152-
def trim_right = meta.strandedness == 'single' ? (meta.damage_treatment == 'none' ? params.damage_manipulation_bamutils_trim_single_stranded_none_udg_right : meta.damage_treatment == 'half' ? params.damage_manipulation_bamutils_trim_single_stranded_half_udg_right : 0) : (meta.damage_treatment == 'none' ? params.damage_manipulation_bamutils_trim_double_stranded_none_udg_right : meta.damage_treatment == 'half' ? params.damage_manipulation_bamutils_trim_double_stranded_half_udg_right : 0)
153-
[meta, bam, trim_left, trim_right]
154-
}
154+
ch_trimbam_input = ch_to_trim
155+
.map {
156+
meta, bam ->
157+
trim_left = meta.strandedness == 'single' ? ( meta.damage_treatment == 'none' ? params.damage_manipulation_bamutils_trim_single_stranded_none_udg_left : meta.damage_treatment == 'half' ? params.damage_manipulation_bamutils_trim_single_stranded_half_udg_left : 0 ) : ( meta.damage_treatment == 'none' ? params.damage_manipulation_bamutils_trim_double_stranded_none_udg_left : meta.damage_treatment == 'half' ? params.damage_manipulation_bamutils_trim_double_stranded_half_udg_left : 0 )
158+
trim_right = meta.strandedness == 'single' ? ( meta.damage_treatment == 'none' ? params.damage_manipulation_bamutils_trim_single_stranded_none_udg_right : meta.damage_treatment == 'half' ? params.damage_manipulation_bamutils_trim_single_stranded_half_udg_right : 0 ) : ( meta.damage_treatment == 'none' ? params.damage_manipulation_bamutils_trim_double_stranded_none_udg_right : meta.damage_treatment == 'half' ? params.damage_manipulation_bamutils_trim_double_stranded_half_udg_right : 0 )
159+
[ meta, bam, trim_left, trim_right ]
160+
}
155161

156-
BAMUTIL_TRIMBAM(ch_trimbam_input)
157-
ch_versions = ch_versions.mix(BAMUTIL_TRIMBAM.out.versions.first())
162+
BAMUTIL_TRIMBAM( ch_trimbam_input )
163+
ch_versions = ch_versions.mix( BAMUTIL_TRIMBAM.out.versions.first() )
158164

159-
SAMTOOLS_INDEX_DAMAGE_TRIMMED(BAMUTIL_TRIMBAM.out.bam)
160-
ch_versions = ch_versions.mix(SAMTOOLS_INDEX_DAMAGE_TRIMMED.out.versions.first())
165+
SAMTOOLS_INDEX_DAMAGE_TRIMMED( BAMUTIL_TRIMBAM.out.bam )
166+
ch_versions = ch_versions.mix( SAMTOOLS_INDEX_DAMAGE_TRIMMED.out.versions.first() )
161167
ch_trimmed_index = params.fasta_largeref ? SAMTOOLS_INDEX_DAMAGE_TRIMMED.out.csi : SAMTOOLS_INDEX_DAMAGE_TRIMMED.out.bai
162168

163-
ch_trimmed_bams = BAMUTIL_TRIMBAM.out.bam.join(ch_trimmed_index)
169+
ch_trimmed_bams = BAMUTIL_TRIMBAM.out.bam.join( ch_trimmed_index )
164170
}
165171

166172
emit:
167-
rescaled = ch_rescaled_bams // [ meta, bam, bai ]
168-
filtered = ch_pmd_filtered_bams // [ meta, bam, bai ]
169-
trimmed = ch_trimmed_bams // [ meta, bam, bai ]
173+
rescaled = ch_rescaled_bams // [ meta, bam, bai ]
174+
filtered = ch_pmd_filtered_bams // [ meta, bam, bai ]
175+
trimmed = ch_trimmed_bams // [ meta, bam, bai ]
170176
flagstat = ch_pmd_filtered_flagstat // [ meta, flagstat ]
171177
versions = ch_versions
172178
}

0 commit comments

Comments
 (0)
⚡