https://nf-co.re logo
Join Slack
Powered by
Hi all, years ago I asked for some operators to ex...
# nextflow-plugins
p
Hi all, years ago I asked for some operators to extract data (samples, chromosome) from common hts files (bam, vcf, fasta) https://github.com/nextflow-io/nextflow/pull/1266 . @p suggested to implement this idea in a plugin and I'm happy to release my first plugin using the htsjdk library ( the Java API for high-throughput sequencing data (HTS) formats )
👏🏻 1
👏 1
wow 7
my plugin is available at https://github.com/lindenb/nf-htsjdk
🚀 1
A first operator 
dictionary
extracts chromosomes from BAM/SAM/CRAM/BCF/VCF/FA files.
Copy code
local_files_ch = Channel.fromPath(["data/*.bam","data/*.sam","data/*.bcf","data/*.vcf","data/*.vcf.gz","data/*.fai"])

remote_files_ch = Channel.of(
	"<http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa>",
	"<http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase1/data/HG00101/exome_alignment/HG00101.mapped.illumina.mosaik.GBR.exome.20110411.bam>",
	"<https://storage.googleapis.com/gcp-public-data--gnomad/release/4.0/vcf/exomes/gnomad.exomes.v4.0.sites.chr1.vcf.bgz>",
	"<https://www.encodeproject.org/files/ENCFF656TNN/@@download/ENCFF656TNN.bam>",
	"<https://github.com/lindenb/jvarkit/raw/master/src/test/resources/toy.cram>"
	)

htsfiles_ch= local_files_ch.mix(remote_files_ch)

htsfiles_ch.dictionary(withLength:true).
	view{"DICTIONARY: ${it}"}
output
Copy code
DICTIONARY: [10, 135534747, <http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase1/data/HG00101/exome_alignment/HG00101.mapped.illumina.mosaik.GBR.exome.20110411.bam>]
DICTIONARY: [9, 141213431, <http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase1/data/HG00101/exome_alignment/HG00101.mapped.illumina.mosaik.GBR.exome.20110411.bam>]
DICTIONARY: [chr10, 133797422, <http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa>]
DICTIONARY: [chr10, 133797422, <https://storage.googleapis.com/gcp-public-data--gnomad/release/4.0/vcf/exomes/gnomad.exomes.v4.0.sites.chr1.vcf.bgz>]
DICTIONARY: [chr10, 133797422, <https://www.encodeproject.org/files/ENCFF656TNN/@@download/ENCFF656TNN.bam>]
DICTIONARY: [chr10_GL383545v1_alt, 179254, <http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa>]
DICTIONARY: [chr11, 135086622, <http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa>]
DICTIONARY: [chr11, 135086622, <https://storage.googleapis.com/gcp-public-data--gnomad/release/4.0/vcf/exomes/gnomad.exomes.v4.0.sites.chr1.vcf.bgz>]
DICTIONARY: [chr11, 135086622, <https://www.encodeproject.org/files/ENCFF656TNN/@@download/ENCFF656TNN.bam>]
DICTIONARY: [chr11_KI270721v1_random, 100316, <http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa>]
DICTIONARY: [chr11_KI270721v1_random, 100316, <https://www.encodeproject.org/files/ENCFF656TNN/@@download/ENCFF656TNN.bam>]
DICTIONARY: [chr11_KI270902v1_alt, 106711, <http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa>]
DICTIONARY: [chr11_KI270903v1_alt, 214625, <http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa>]
(...)
😍 1
a second operator 
samples
extracts the samples from SAM/CRAM/BAM/VCF/BCF files:
Copy code
htsfiles_ch.
	filter{!(it.toString().endsWith(".fai") || it.toString().endsWith(".fa"))}.
	samples(defaultName:"NO_SAMPLE").
	view{"SAMPLE: ${it}"}
output:
Copy code
SAMPLE: [HG00101, <http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase1/data/HG00101/exome_alignment/HG00101.mapped.illumina.mosaik.GBR.exome.20110411.bam>]
SAMPLE: [NO_SAMPLE, <https://storage.googleapis.com/gcp-public-data--gnomad/release/4.0/vcf/exomes/gnomad.exomes.v4.0.sites.chr1.vcf.bgz>]
SAMPLE: [NO_SAMPLE, <https://www.encodeproject.org/files/ENCFF656TNN/@@download/ENCFF656TNN.bam>]
SAMPLE: [S1, <https://github.com/lindenb/jvarkit/raw/master/src/test/resources/toy.cram>]
SAMPLE: [S1, /home/lindenb/src/nf-htsjdk/data/S1.rota.bam]
SAMPLE: [NO_SAMPLE, /home/lindenb/src/nf-htsjdk/data/rotavirus_rf.noRG.sam]
SAMPLE: [S1, /home/lindenb/src/nf-htsjdk/data/rotavirus_rf.bcf]
SAMPLE: [S2, /home/lindenb/src/nf-htsjdk/data/rotavirus_rf.bcf]
SAMPLE: [S3, /home/lindenb/src/nf-htsjdk/data/rotavirus_rf.bcf]
SAMPLE: [S4, /home/lindenb/src/nf-htsjdk/data/rotavirus_rf.bcf]
SAMPLE: [S5, /home/lindenb/src/nf-htsjdk/data/rotavirus_rf.bcf]
SAMPLE: [NO_SAMPLE, /home/lindenb/src/nf-htsjdk/data/vcf_no_genotype.vcf]
SAMPLE: [S1, /home/lindenb/src/nf-htsjdk/data/rotavirus_rf.vcf.gz]
SAMPLE: [S2, /home/lindenb/src/nf-htsjdk/data/rotavirus_rf.vcf.gz]
SAMPLE: [S3, /home/lindenb/src/nf-htsjdk/data/rotavirus_rf.vcf.gz]
SAMPLE: [S4, /home/lindenb/src/nf-htsjdk/data/rotavirus_rf.vcf.gz]
SAMPLE: [S5, /home/lindenb/src/nf-htsjdk/data/rotavirus_rf.vcf.gz]
🤯 1
and the 
build
operator adds the name of the build (defined in ``./plugins/nf-htsjdk/build/resources/main/META-INF/builds.xml`)`:
Copy code
htsfiles_ch.build().
        view{"BUILD: ${it}"}
output
Copy code
BUILD: [GRCh38, <http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa>]
BUILD: [GRCh37, <http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase1/data/HG00101/exome_alignment/HG00101.mapped.illumina.mosaik.GBR.exome.20110411.bam>]
BUILD: [GRCh38, <https://storage.googleapis.com/gcp-public-data--gnomad/release/4.0/vcf/exomes/gnomad.exomes.v4.0.sites.chr1.vcf.bgz>]
BUILD: [GRCh38, <https://www.encodeproject.org/files/ENCFF656TNN/@@download/ENCFF656TNN.bam>]
BUILD: [samtools, <https://github.com/lindenb/jvarkit/raw/master/src/test/resources/toy.cram>]
BUILD: [rotavirus, /home/lindenb/src/nf-htsjdk/data/S1.rota.bam]
BUILD: [rotavirus, /home/lindenb/src/nf-htsjdk/data/rotavirus_rf.noRG.sam]
BUILD: [rotavirus, /home/lindenb/src/nf-htsjdk/data/rotavirus_rf.bcf]
BUILD: [rotavirus, /home/lindenb/src/nf-htsjdk/data/vcf_no_genotype.vcf]
BUILD: [rotavirus, /home/lindenb/src/nf-htsjdk/data/rotavirus_rf.vcf.gz]
BUILD: [rotavirus, /home/lindenb/src/nf-htsjdk/data/rotavirus_rf.fa.fai]
wow 1
that's it. 😛
@Ben Sherman I just saw your comment : https://github.com/nextflow-io/nextflow/pull/1266#issuecomment-2231376536
Those operators should just be functions. Then you can use them in a map operator or a process or wherever you want
Thanks for the suggestion. what do you mean here ? Can I use a simple function to generate multiple output from a channel (chr1, chr2, chr3, etc...) or do you think about using something like 
flatMap
?
b
Instead of implementing them as operators (right now they are essentially specializations of the 
map
operator), just implement them as functions which process a single input/output. That gives you more flexibility and you don't have to deal with channels and operators in your plugin code. I will try to rewrite your example to show what I mean:
Copy code
htsfiles_ch
  .map { hts -> dictionary(hts, withLength:true) }
  .view { "DICTIONARY: ${it}" }

htsfiles_ch
  .filter { hts -> !(hts.endsWith(".fai") || hts.endsWith(".fa")) }
  .map { hts -> samples(hts, defaultName:"NO_SAMPLE") }
  .view { "SAMPLE: ${it}" }

htsfiles_ch
  .map { hts -> build(hts) }
  .view { "BUILD: ${it}" }
✅ 1
This way it's also clear that you're just mapping over some collection of inputs This happens all the time in plugins (e.g. nf-validation) where people implement custom operators that turn out to just be specialized 
map
operators. Unfortunately we started the idea with operators like 
splitCsv
in Nextflow, but we are trying to course-correct
p
While it’s possible apply any function as an operator via map or flatMap. The “native” approach is better at least of reasons: • easier for the user • above all, scale better with a large data, since the operation is non-blocking.
the latter applies in particular to this use-case, because a function would need to read all the bam file, return a BIG list that then is converted to a channel
Using a native (operator) approach, the channel can emit a value as soon as it is parsed from the bam file.
same for 
slitCsv
this is why it does not need to be corrected
p
@p I'm not sure I understand; What is the "native approach" here ? a "simple" function as @Ben Sherman suggested 🤔 ?
b
scale better with a large data, since the operation is non-blocking
this is not true. there is no difference under the hood between Pierre's operators and the function equivalent that I proposed
p
that’s pierre problem
❓ 1
b
?
same thing for 
splitCsv
. whether it is a function or operator, you are still mapping over some channel of values. the only difference is whether the 
map
operator is explicit in the user's code
p
the operator allows you to implement a stream
look at splitter code
it can emit each entry while reading the file
p
well, most functions in htsjdk reads all the data at once. For example the function extracting a dictionary from a VCF/BAM/FAI reads ALL the data in memory and returns a 
List
(e.g: https://samtools.github.io/htsjdk/javadoc/htsjdk/htsjdk/variant/vcf/VCFHeader.html#getSequenceDictionary-- ), there is no code to emit the the data on the fly. So if you tell me that a better approach is to move the code to a function instead of an operator, I can do this.
b
Paolo is right about 
splitCsv
because it is really a flatMap rather than a map, so it is useful to stream each collection. Though I would argue that 
flatMap
should be able to flatten a stream just like a list, so that you can just do 
flatMap { v -> f(v) }
instead of creating a new operator every time you want to use this pattern. But that is another topic 😄 The hts operators seem to be a map rather than a flatMap so I think my advice stands
✅ 1
p
Thank you all, I've update my code to move the operators as a set of functions. This is indeed more easy as object are returned and the user can call the getters to compose the content of a channel. https://github.com/lindenb/nf-htsjdk e.g:
Copy code
htsfiles_ch.                
    map{[dictionary(it),it]}.
    flatMap{row->row[0].getSequences().collect{dict->[dict.getContig(),row[1]]}}.
    filter{it[0].matches("((chr)?[12]|RF.*)")}.
    view{"DICTIONARY: ${it}"}


DICTIONARY: [chr1, <http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa>]
DICTIONARY: [chr2, <http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa>]
DICTIONARY: [1, <http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase1/data/HG00101/exome_alignment/HG00101.mapped.illumina.mosaik.GBR.exome.20110411.bam>]
DICTIONARY: [2, <http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase1/data/HG00101/exome_alignment/HG00101.mapped.illumina.mosaik.GBR.exome.20110411.bam>]
DICTIONARY: [chr1, <https://storage.googleapis.com/gcp-public-data--gnomad/release/4.0/vcf/exomes/gnomad.exomes.v4.0.sites.chr1.vcf.bgz>]
DICTIONARY: [chr2, <https://storage.googleapis.com/gcp-public-data--gnomad/release/4.0/vcf/exomes/gnomad.exomes.v4.0.sites.chr1.vcf.bgz>]
DICTIONARY: [chr1, <https://www.encodeproject.org/files/ENCFF656TNN/@@download/ENCFF656TNN.bam>]
DICTIONARY: [chr2, <https://www.encodeproject.org/files/ENCFF656TNN/@@download/ENCFF656TNN.bam>]
DICTIONARY: [RF01, /home/lindenb/src/nf-htsjdk/data/S1.rota.bam]
(...)
🙌 1
Open in Slack
PreviousNext
Powered by