https://nf-co.re logo
Join Slack
Powered by
:mega: Calling all maintainers of sequence aligner...
# modules
j
📣 Calling all maintainers of sequence aligner-based modules 📣 Before I go down a 'finding a solution' rabbit hole: A user just reported a case of the 
bowtie2/align
module replacing the context of the reference FASTA file with a SAM file 😱 We've identified the problem being that the 
find
command used to pick up the index directory was finding two sets of 
*.bt2
index files, and unfortunately the prefix of the second set was the same name as the reference FASTA itself. Meaning this:
Copy code
bowtie2 \\
        -x \$INDEX \\
        -1 sample.fastq
ended up resolving to
Copy code
bowtie2 \\
        -x  ./hs37d5/bowtie2 ./hs37d5/hs37d5.fa \\
        -1 sample.fastq
With bowtie2 assuming the first non-argumented file was where the output should go into... So I'm wondering what the next course of action should be... (improve the find pattern? Add a validation check and fail if n > 1 index files) TL;DR has anyone else come across a similar case of an aligner module overwriting an index, and already addressed this with an updated 
find
command I can copy, and/or have a suggestion of a preferred way to counter this?
e
Good catch, and probably relevant to @Adam Talbot 's presentation of yesterday? Would an explicit input parameter make sense, instead of manipulating files with 
find
here?
j
HAven't had a chance to see it yet unfortunately. The problem is bowtie2 (and stuff like bwa) make multiple reference files. If you want to explicitly input parameter it, you will need to add like somewhere between 6-12 different input parameters just for essentially a single index...
e
Hmm I think I know what you mean. And most often not the same number of files. I'm wondering if it makes sense to capture these with an expression instead of a whole directory, and then use them in the next modules input, and if we want just the name, parse the file[0] or something
a
The problem is bowtie2 (and stuff like bwa) make multiple reference files.
It might be easier to provide a tuple input: 
tuple val(name), path(bwt), path(fa), path(bt2)
which is explicit rather than implicit. Sounds like a lot of work. BYEeEEEE 🏃
👍 1
🌟 1
😂 4
"easier"
j
😆
Nextflow pipeline development easier yes
But a good way to have users telling us to 'fuck off'
a
It kinda makes sense with @Ben Sherman's record types, you could 'just' supply a 
BwtIndex
record as an input:
Copy code
input:
    val id
    path forward_fastq
    path reverse_fastq
    BwtIndex bwt_index
j
But how does 
bwt_index
get picked up up though is the issue
OR you mean it does the validation for you?
a
Oh I've done that
😆 1
j
(the record type expects X number and if ti doesn't match error out?)
a
Sorry secret repo I can't show you but I can copy+paste the code
❤️ 1
Copy code
// Detect BWA index files (auto-discover based on FASTA path)
bwa_pattern = bwa ?: fasta + ".{bwt,pac,ann,amb,sa}"

bwa_files = files(bwa_pattern, checkIfExists: true)
ch_bwa = Channel.of(bwa_files)
                .map { bwa_file_list -> validateBwaIndexFiles(bwa_file_list) }
                .map { bwa_file_list -> tuple([id: "bwa"], bwa_file_list) }
you could simplify it to:
Copy code
map { bwa_file_list -> createBwtIndex(bwa_file_list) }
j
and validateBwaIndexFiles includes the record type, or is this a separate function?
But I'm also going to throw in a spanner in the works: with bowtie2, you don't need to have the prefix the same as the FASTA file name 😬 🎉
a
Well that's OK, you just ignore the fasta bit:
Copy code
bwa_pattern = params.bowtie2_index + ".{bwt,pac,ann,amb,sa}"
It's just constructing a string which you pass to 
file()
Going back to the original problem. How did two bt2 index files end up in the process? Can we fix that and be more specific about which files are in that index input?
j
OK, good we reached my exact question! (but that was a useful tangent, to be clear!)
tux sorry 1
So that's the thing, currently we use 
find
Copy code
INDEX=`find -L ./ -name "*.rev.1.bt2" | sed "s/\\.rev.1.bt2\$//"`
    [ -z "\$INDEX" ] && INDEX=`find -L ./ -name "*.rev.1.bt2l" | sed "s/\\.rev.1.bt2l\$//"`
    [ -z "\$INDEX" ] && echo "Bowtie2 index files not found" 1>&2 && exit 1
🤮 1
And that pattern is used a lot in other modules too 😉
With the (reasonable) assumption you would only have one set of files
a
Can we do this:
Add a validation check and fail if n > 1 index files
In Nextflow?
In the channel?
j
That's the the thing: can we have 
nextflow
discover the files already? Or do we rely on the existing 
find
pattern
Given that what goes into the module is a directory (that should contain the index files)
a
Nextflow can do it...I'm assuming the input is something like 
bowtie_index/*
?
j
Could be named anything
Ideally we want to do this in the module
rather than pipeline level code (as a caveat)
a
errrr in the `script`block? 
def bowtie_index = files.anyOf { groovy goes here }
?
🏴‍☠️ 1
j
Maybe? I've always had mixed experiecne when it got into the groovy bit 😅
When trying to get it to list files from a directory for example
and then getting it to reliably count
(but that is my inexperience)
a
From my reading of the docs, you have a list of files like so: 
test.1.bt2  test.2.bt2  test.3.bt2  test.4.bt2  test.rev.1.bt2  test.rev.2.bt2
The name of the index is completely arbitrary when you run bowtie2-build, in ☝️ case it's 
test
In the 
find
, we detect the 
.rev.1.bt2
, then remove that suffix. If that fails, we find the 
.rev.1.bt2l
. If that fails we error. When people are using the BOWTIE2_ALIGN module are they generally supplying a directory or list of files?
👍 1
j
Directory
AFAIK
(but anecdotal from the few pipelines I know are using it)
👍 1
a
It makes it messier but we can support both
👍 1
it ain't pretty:
Copy code
script:
// Handle both directory and file collection inputs
def index_files = index instanceof List ? index : file(index).listFiles()?.toList() ?: [index]

def idx_file = index_files.find { it.name.endsWith('.rev.1.bt2') || it.name.endsWith('.rev.1.bt2l') }
if (!idx_file) error "Bowtie2 index files not found in: ${index}"
def index_base = idx_file.name.replaceFirst(/\.rev\.1\.bt2l?$/, '')

// If input was a directory, prepend the directory name
def index_path = index instanceof List ? index_base : "${index.name}/${index_base}"

def args = task.ext.args ?: ''
def prefix = task.ext.prefix ?: "${meta.id}"

"""
bowtie2 \\
    -x ${index_path} \\
    -p ${task.cpus} \\
    ${args} \\
    ...
"""
but in that method, the fasta file would never appear in the 
${index_path}
, because it would only ever use a single value.
but it might pick the wrong one 🤔
j
wrong one as in if you had multiple sets?
of index files?
a
The first problem had two files:
Copy code
./hs37d5/hs37d5.rev.1.bt2
./hs37d5/hs37d5.fa.rev.1.bt2 // Nextflow might pick this one!
j
Yeah exactly OK
a
Which is better but not correct
j
But could we then on the 
idx_file
definition have a check that the result should be <= 1
?
find{}
-> 
findAll{}
(if that's a thing?)
a
Claude wrote this:
Copy code
script:
def index_files = index instanceof List ? index : file(index).listFiles()?.toList() ?: [index]

// Find ALL potential index bases, not just the first
def idx_matches = index_files.findAll { it.name.endsWith('.rev.1.bt2') || it.name.endsWith('.rev.1.bt2l') }

if (idx_matches.size() == 0) {
    error "Bowtie2 index files not found in: ${index}"
}
if (idx_matches.size() > 1) {
    error "Multiple Bowtie2 index sets found: ${idx_matches.collect { it.name.replaceFirst(/\.rev\.1\.bt2l?$/, '') }.join(', ')}. Please provide only one index."
}

def index_base = idx_matches[0].name.replaceFirst(/\.rev\.1\.bt2l?$/, '')
def index_path = index instanceof List ? index_base : "${index.name}/${index_base}"

// ... rest of script
j
HA! SEEI DONT NEED NO AI
I think that's nice...
Actually this may fix another issue we had with someone mounting files on a samba based NAS and running on WSL2 and not working with 
find
😅
a
The pipeline developer has the same problem, but at least it's an explicit error. I would suggest cleaning up the inputs to the process if possible, something is being greedy and putting too many files into the process.
👍 1
There's probably a stray glob munching up files
j
Yeah in retrospect this would be better but I think we need to make the bowtie2 module more robust given it is already 'in production' so to say
This has been very helpful and insightful for me, thank you very much @Adam Talbot
I will try out the suggestions above nexgt time I'm able to do some nextflow coding and see if this fixes the issue
And then maybe can roll out to other modules
a
The best method might be to go upstream to BOWTIE2_INDEX:
Copy code
output:
tuple val(meta), val(index_base), path("${index_base}*.bt2{,l}"), emit: index
Which can be directly plumbed into `BOWTIE2_ALIGN`:
Copy code
input:
tuple val(meta), path(reads)
tuple val(meta2), val(index_base), path(index_files)

script:
"""
bowtie2 \\
    -x ${index_base} \\
this 2
e
Yep, that was my thought, and would be my preference :D
j
Doesn't help when someone is supplying there own - big waste of compute to regenerate on each run
e
Hmm, is this mutually exclusive? If someone provides their own index then the 
BOWTIE2_INDEX
does not have to run. But then the downstream module would definitely need the extra checking step to make sure it's the proper format. I see your point
a
you would calculate it like the def script, but in the workflow (e.g. 
map
)
j
AFAICR listFiles() had issues when we tried it before: https://github.com/nextflow-io/nextflow/issues/6488 Looks like the replacement is listDirectory() but it will only be in newer NF versions
🚫 1
s
Checking that the prefix that is found is only one fine and erroring if not makes a lot of sense. Could even be done in bash rather than Groovy
m
@Yumi Sims found that 
listFiles()
doesn’t work on symlinks. Need to call 
toRealPath()
first to resolve the symlink.
👍 1
a
Symbolic links are platform specific so always cause headaches facepalm Does the listFiles() not work when used on the symbolic link or when symbolic links are within the directory. e.g. which one of these:
Copy code
symlink.listFiles() -> [ file1, file2, ... ]
realfile.listFiles() -> [ symlink1, symlink2, ...]
using 
.toRealPath()
won't work on remote filesystems
m
listFiles()
doesn’t work when used on the symbolic link itself
a
j
a
my fix there addresses the underlying method (
listFiles0()
) so it should help both 
listFiles
and 
listDirectory
🤞
🎉 3
👍 2
j
So we need to wait until the next release and we can clean up a lot of stuf then!?!
Super fast, thank you @Adam Talbot!
a
it aint merged yet!
j
tenor_gif7183277751414770962.gif
a
Alright the change got rolled back: https://github.com/nextflow-io/nextflow/pull/6808#issuecomment-3867729501 but you should be able to use 
listDirectory()
in the near future: https://github.com/nextflow-io/nextflow/pull/6821#event-22676287462
❤️ 1
Open in Slack
PreviousNext
Powered by