-
Notifications
You must be signed in to change notification settings - Fork 21
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Add subsampling #50
Add subsampling #50
Conversation
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thanks a lot for your contribution! The subsampling functionality will be very important to Seqinspector, therefore it is great that you decided to work along this line!
Unfortunately, there is still some room for improvement, but I am optimistic that all of those issues can be resolved during the next two days of the Hackathon.
Generally, I think, that most users would want to select a fraction of reads to subsample rather than deciding for an absolute number of reads. Therefore, your initial idea of using SeqKit sample, which has direct support for probabilities, would make a lot of sense, however, require writing the module first. If going with Seqtk, I think it would be fine merging this PR with absolute numbers only, but if you feel enthusiastic, you could patch the module to include a read counting step and set the sample size according to the number of input reads.
script: | ||
def args = task.ext.args ?: '' | ||
def prefix = task.ext.prefix ?: "${meta.id}" | ||
if (!(args ==~ /.*-s[0-9]+.*/)) { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Do you think it would be useful to expose the random seed of the tool as (possibly hidden) pipeline parameter?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'm not sure about this, is there any example you have in mind?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I am unsure as well. My experience from the RNA-seq pipeline is, that people usually struggle with custom configuration files.
The extra_*
convenience parameters of the pipeline (extra_star_align_args
or extra_ kallistro_quant_args
), which allow modifying the arguments of those tools without a custom config are very popular.
Therefore, I about setting the seed via a parameter, but on a second thought, I agree that the need for a custom config is probably fine.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I can understand the difficulties with custom configuration files, the extra arguments are ok, but on the other side I have the impression that might be easily forgotten when setting the workflow. I think may be better to check if there are paired reads by checking the input channel or just if there is _1 and _2 files, and if present, enables the -s option?
workflows/seqinspector.nf
Outdated
// MODULE: Run Seqkit sample to perform subsampling | ||
// | ||
if (params.sample_size > 0 ) { | ||
ch_sample_sized = SEQTK_SAMPLE(ch_samplesheet.map { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Unfortunately, I don't think the module is capable of processing paired files. The description of the tool suggests that it has to be invoked with the same random seed on each FastQ file separately:
Subsample 10000 read pairs from two large paired FASTQ files (remember to use the same random seed to keep pairing):
seqtk sample -s100 read1.fq 10000 > sub1.fq
seqtk sample -s100 read2.fq 10000 > sub2.fq
and all of the module tests only use a single file.
When you invoke the module, you will have to account for that.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Ok, you know if at least one of the tests have paired reads so I can view how the input channel looks like and also test this functionality?
I have found this post about custom-configuration-files and looking the module, the '-s' option is available using it as ext.args.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I have asked in Barcelona, and apparently there isn't yet a paired-end test, so I has been suggested to open an issue to add it: #55
For now I have added the seed argument in the conf/modules.config
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thanks for the review!
workflows/seqinspector.nf
Outdated
// MODULE: Run Seqkit sample to perform subsampling | ||
// | ||
if (params.sample_size > 0 ) { | ||
ch_sample_sized = SEQTK_SAMPLE(ch_samplesheet.map { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Ok, you know if at least one of the tests have paired reads so I can view how the input channel looks like and also test this functionality?
I have found this post about custom-configuration-files and looking the module, the '-s' option is available using it as ext.args.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think the major issues are addressed, and the other things can be polished later respectively are a matter of taste anyway. Thanks for your contribution!
@@ -18,6 +18,10 @@ process { | |||
saveAs: { filename -> filename.equals('versions.yml') ? null : filename } | |||
] | |||
|
|||
withName: SEQTK_SAMPLE { | |||
ext.args = '-s100' | |||
} |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I have no opinion on whether the pipeline should publish the subsampled files. But if you want to output them in the way you describe in the output.md
, you will probably need to define a corresponding publishDir
directive here? (There is also a new possibility called Workflow Output Schema, that I am however not yet familiar with. Did you use that without me noticing it?)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think is always output because there is a general statement here:
seqinspector/conf/modules.config
Line 15 in a2d589f
publishDir = [ |
Closes #28.
Default 0 means that no subsampling is performed.
The subsampling step is performed at the very beginning of the pipeline.
Regarding this:
SeqKit sample already uses a default according to the documentation:
https://bioinf.shenwei.me/seqkit/usage/#sample see
-s
Snapshots updated after merging #49
PR checklist
nf-core lint
).nf-test test main.nf.test -profile test,docker
).nextflow run . -profile debug,test,docker --outdir <OUTDIR>
).docs/usage.md
is updated.docs/output.md
is updated.CHANGELOG.md
is updated.README.md
is updated (including new tool citations and authors/contributors).