Hi,
I am conducting virtual screening of some compounds and I am having issues with the docking phase. I wonder why the .cross operator is not generating any output for my ligand vs receptor pairs. I will appreciate any assistance.
Thank you.
You need to provide a code example, otherwise it’s impossible to figure out what’s going on.
workflow {
// Step 0: Load ligand files
ligand_files = Channel.fromPath("${params.ligands}/\*.sdf", checkIfExists: true)
// Step 1: Split multi-ligand SDFs
split_sdf_files = splitLigand(ligand_files).flatten()
// Step 2: Rename by ZINC ID
renamed_sdf = renameSDF(split_sdf_files)
// Step 3: Minimize renamed ligands
minimized_sdf = ligandMinimization(renamed_sdf)
// Step 4: Convert minimized ligands
prepared_ligands = convertLigandtoPDBQT(minimized_sdf)
prepared_ligands_flat = prepared_ligands.flatten()
// Step 5: Load receptor files
receptor_files = Channel.fromPath("${params.receptors}/\*.{pdb,sdf,mol2}", checkIfExists: true)
// Step 6: Convert receptor files to .pdb format
converted_pdbs = convertReceptortoPDB(receptor_files)
// Step 7: Convert .pdb receptor files to .pdbqt
prepared_receptors = convertReceptortoPDBQT(converted_pdbs)
prepared_receptors_flat = prepared_receptors.flatten()
// Call docking process
receptor_centers = Channel.fromPath(‘results/receptor_centers.csv’)
.splitCsv(header:true)
.map { row ->
rec_file = file("results/prepared_pdbqt/${row.receptor}.pdbqt")
\[ rec_file, row.cx.toBigDecimal(),
row.cy.toBigDecimal(),
row.cz.toBigDecimal() \]
}
.filter { it\[0\].exists() } // skip missing files
// Step 9: Cartesian product of receptors and ligands
docking_pairs = prepared_receptors_flat.cross(prepared_ligands_flat)
// Step 10: Run docking
vina_docking(docking_pairs)
The easiest way to debug this is to use the view operator to view the upstream channels. If the cross isn’t producing any output it’s likely because one of the input channels are empty.
Why are you escaping the square brackets in your script (\[)?
The input channels give the required output. .cross is not generating the required output, I wonder if I am using it correctly.
That is probably due to the length of the code. It is without that on my end.
In that case, can you give an example that just creates some fake input channels and calls the cross operator? That way I can run it for myself and see what’s wrong
#!/usr/bin/env nextflow
nextflow.enable.dsl=2
workflow {
ligands = Channel.of("lig1.pdbqt", "lig2.pdbqt", "lig3.pdbqt")
receptors = Channel.of(
["recA.pdbqt", 10.0, 20.0, 30.0],
["recB.pdbqt", 40.0, 50.0, 60.0]
)
pairs = ligands.cross(receptors)
}
Okay I see. The caveat with the cross operator is that it only crosses inputs that have a matching key. Despite it’s name, it does not do a full cross product.
You can specify a closure to define the matching key in terms of an input value. So you could make this closure return a constant value to force a full cross product:
ligands.cross(receptors) { v -> 0 }
Honestly this should probably be the default behavior, but this was before my time 
PAIR: [lig1.pdbqt, recA.pdbqt]
PAIR: [lig1.pdbqt, recB.pdbqt] - it does not ouput all pairs. Is there another way to achieve a cartesian product for 50 x 50 docking?
#!/usr/bin/env nextflow
nextflow.enable.dsl=2
workflow {
ligands = Channel.of("lig1.pdbqt", "lig2.pdbqt", "lig3.pdbqt")
receptors = Channel.of(
["recA.pdbqt", 10.0, 20.0, 30.0],
["recB.pdbqt", 40.0, 50.0, 60.0]
)
docking_pairs = receptors.flatMap { receptor_list ->
ligands.map { ligand ->
[
receptor_list[0], // receptor file
receptor_list[1], // cx
receptor_list[2], // cy
receptor_list[3], // cz
ligand // ligand file
]
}
}
docking_pairs.view { p -> "DOCK_PAIR: ${p}" }
}
>> DOCK_PAIR: DataflowBroadcast around DataflowStream[?]
DOCK_PAIR: DataflowBroadcast around DataflowStream[?]
If you want the outputs to be flat, you can use the combine operator:
ligands.combine(receptors)
1 Like
This topic was automatically closed 7 days after the last reply. New replies are no longer allowed.