text:
How do I interpret the function O(f^-2 log n) in Shekhar et al. (2017)?
The results of Shekhar et al. (2017) "indicate that ASTRAL requires
$\mathcal{O}(f^{-2} \log{n})$

gene trees to reconstruct the species tree correctly with high probability where n is the number of species and f is the length of the shortest branch in the species tree.
How do I interpret this function, specifically the big O notation, when applying it to my data to estimate whether I have enough loci?

Please [edit] your question and provide a link or a full citation to the paper you are quoting. There are 86 papers in PubMed matching (Shekhar[Author]) AND ("2017"[Date - Publication] : "2017"[Date - Publication]) . Also clarify what you mean by applying this function. Why would you apply it?

I've updated your post and added a link to the actual paper you mean. A preprint is available here: https://arxiv.org/abs/1704.06831

Big O notation means "of the order of". Thus $\mathcal{O}(f^{-2} \log{n})$ means that the number of trees required goes with $\log{n}$ for number of species (i.e. if you square the number of genes, you should double the number of trees) and with the inverse square of the shortest branch length (i.e. if you double the shortest branch length, you quarter the number of trees needed).
Big O notation generally allows to to compare two input sizes, but cannot really tell you anything about a single point. I.e. it is hard (impossible?) to tell from this function whether you have sufficient loci.


text:
Genomics Statistics Problem
I have used a python script to identify target sequences in a DNA sequence file. There are two classes of sequence: coding and non-coding.
I have identified 728 sequences of interest. 597 of these fall into the coding regions and 131 of these fall into the non-coding regions. This is the equivalent of 18% non-coding. The total non-coding region in the sequence file is 13%.
Is there a statistical tool to demonstrate the python script identified target sequences in a non-random fashion?
If the script identified sequences that were randomly distributed then 13% of them would have been found in the non-coding region, from a total of 728 this seems like it should be reliable.

It sounds like you're fishing for a statistic to prove a point. I wouldn't try to over-analyze your data. You might just consider stating that 18% of the sequences are in non-coding regions which is greater than the expected null hypothesis of 13% of the sequences being in non-coding regions.

One could ask how big are those sequences of interest (SOI) and if there are some of them that span between coding and non-coding regions. Also the distribution of SOI could depend on TADs or other higher structures of your sequence (if it is a crhomosome size sequence). To model the distribution of SOI would require knowing a lot more about the SOI and the sequence you are working with and without a null model it cannot be tested if it is randomly following this model or not

Hi again @conchoecia. Thanks for the feedback. With this point, however, we find ourselves in disagreement. There is a chance that if my results are due to random chance there may well be an uneven distribution. Say if my non-coding hits were 14% it would be obvious this would not be reliable evidence in support of my hypothesis. Furthermore, I have no problem with obtaining a negative result but however I interpret my results it needs to be statistically reliable, there is a chance that 18% is not significantly different to 13% and therefore cannot be used in favour of my hypothesis.

@llrs Hi. The sequences of interest are 29 bp (promotors) and I am searching chromosome FASTA. There are a some overlaps, but I have not counted them as non-coding for simplicity. The non-coding regions were defined by a different study (https://doi.org/10.1186/s13068-017-0742-z) Table 1.

Thanks everyone, question answered on cross-validated.
https://stats.stackexchange.com/questions/386173/genomics-stats-problem/386177#386177

If your null hypothesis is that the event of being selected by your program is statistically independent from the event of being non-coding, you can use Fisher's exact test to reject the null hypothesis.
Suppose you have 10,000 (I'm making that up but presumably you have the real number) total sequences, 8,700 coding and 1,300 noncoding. Given that your program selects 728 sequences, you would expect around 13%*728 ~= 95 noncoding sequences. But your program picked 131 noncoding sequences and you want to know if 131 is far enough away from 95 to conclude that the difference is not just due to chance.
The probability of selecting 131 or more of these noncoding sequences (given that we are choosing 728 total sequences) is given by the so-called hypergeometric distribution. The two-sided test will also incorporate the probably of getting an equally unlikely low number (in this case, somewhere around 57) or less.
In R for example, you can calculate this using the 2-by-2 contingency table as follows:
fisher.test(rbind(c(597, 131), c(8103, 1169)))

which gives:
p-value = 7.399e-05

So in this situation (and assuming the probability model described), there's around a 0.0074% chance of observing a ratio as "extreme" as the one you observed, which indicates it's unlikely that the two events are independent.
Note: As others have mentioned, a binomial test could be appropriate too. If your number of total sequences is very large, then the two results will converge. But if it isn't, then Fisher's exact test will be more appropriate (essentially since you're sampling without replacement).

A Poisson or binomial would suffice.
I'm not at all sure the difference will significantly deviate from randomness.
Perl: 95% confidence interval of a non-parametric bootstrap for e.g. 1000 replicates is an alternative approach. I'd be more comfortable coding that solution in Perl than I would in Python, so is unlikely to be of much use here. Anyway, I would use a 2-dimensional array for your total data-set, the second dimension being coding/non-coding, I would then use the inbuilt random-number generator and multiply by the size of the array to randomly pluck out each element. 95% CI of the mean works by ranking the data and removing the top and bottom 2.5% and taking the interval for a single replicate (which is repeated 999 times more).
There will be a deep learning approach to this data, which would be easier to implement in Python, but I'd need to think about this because it is not simple classification.
Again my warning is that given the sample size, the result at best might be borderline. There are likely to be lots of interesting in your data, but coding/non-coding may not be one of them.

Hi, Michael. Thank you. Yes, I am trying binomial analysis now, don't worry about whether this results in a null hypothesis. I am simply trying to rule out a possible explanation for some results I obtained in the wet-lab. Whatever result I get is fine by me I am just happy to show I have properly investigated possible explanations.

Good to know. I'm definitely advocating DL for molecular solutions, they are very, very cool.

Haha, I'm interested in DL however my data science, coding and statistics skills need a lot of work before I can use such a tool. But I read a lot about the applications and I'm keen to get to the point where I can use such tools.

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
Power law distribution alpha values
I did a Pearson co-expression analysis for generating networks for my tissue-specific (chondrocytes) RNA-seq data. I used R package poweRlaw to check for the power-law distribution. We got the following values:
 Xmin = 14
 goodness of fit = 0.035
 α = 4.76

My question is what does alpha (α) of more than 4 indicates?
Is it is okay to have such values?

Thank you

It means the distribution is heavily skewed and has a 'long-tail', thus the variance is much greater than the mean. In other words most values are around zero, with a small number of values having extreme values.
[
You need to check the goodness of fit because that is saying the deviation from a power-laws distribution is just significant. So it ain't quite power-laws, but the general gist of a 'long-tailed distribution' is correct.
Its a good result and a nice analysis, but you are not quite there at present (needs to be >0.05). Check the negative binomial distribution.

I've seen you data, it looks very good and the power series is not a bad explanation of the data. There are two possibilities for the goodness of fit being 0.035

The goodness of fit is low because of 'outliers' at the beginning of the distribution
An alternative distribution, such as a negative binomial distribution fits better
Check it is 0.035 not 0.35 because that changes everything.

The best way to do this is to plot the power series expected distribution onto the graph and look where the fit 'fails'.
If it is not fitting at the beginning of the distribution thats is about outliers and easily solvable. If it is not really capturing the shape then its probably a slightly different distribution.
If its point 1 then just remove the very first point in the graph and try again and check the p-value. Or take the mean of the clump of three points next to the first point and try again (or just remove those three point). Don't do any more than that and the technical word is "outlier" analysis. In your write-up you must explain the issue around any 'outliers' and that is perfectly acceptable. If it is a power-series distrupted by an outlier(s) there is formal analysis to demonstrate an 'outlier' using outlier analysis if any one ever asks you.
If its point 2. then use e.g. negative bionomial or the gamma distribution (alpha <1) and check the goodness of fit.
You want a goodness of fit > 0.05 and preferably >0.1 to be sure its the right distribution. There are formal statistics to do this, but that is the main one.
My guess is there are outliers at the beginning because your sample size is not large enough to identify the variation correctly, hence the fit is low. In 'long-tails' its easy to have a 'low sample size' because of the way they work. Thus in this case mapping a distribution is a perfect way to compensate for the sampling.

Thank you for your reply. I have now attached my result graph and I would really appreciate it if you could have a look at it

Hi @Priya I've modificed the answer above

Thank you very much for explaining it so nicely.

@Michael, it is not clear to me whether alpha here represents skewness or something else. Similarly, goodness of fit interpretation does not seem to help.

The power laws distribution is skewed, but alpha will be 'more skewed' as the value increases. The value you have will be optimal for your the data set you plotted. I'm not sure what you mean by goodness of fit "doesn't seem to help".

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
How to combine multiple files into one file?
I have multiple files (n=86000) with one column each and I want to combine them all into one file with 86000 columns.
I tried the following command
Paste *.txt > mean.txt

And as a for loop
for i in *.txt, do paste $i > mean.txt; done

This did not appear to work because there are 86000 to combine, I then divided my files into multiple folders with each folder not containing more than 15000 files but this did not working
Format of my single file is as follows (It has 16000 lines)
0.047202943506056
0.0383549703366668
0.0159162908608477
0.0238048058550539
0.0694728928685859
0.0389052658410832

What is the bug in the code?

Thanks. Do you want one long list? In this case the command is 'cat' not 'paste'. Just check you are not  exceeding your disk space.

What is the error message you get?

Equivalently answered here assuming OP wants to concatenate vertically, not horizontally, though the approach could be modified anyway: https://stackoverflow.com/questions/18695105/how-to-concatenate-huge-number-of-files

Re-reading the question: 16000 lines x 86000 files ... this is not alot. I think the error report is important, or else perhaps it just hangs?

@ Fabio Marroni There is no error. When I check the output file it has only one column.

@Michael: No don't want one long list. I want an output file with 86000 columns.

I'm voting to close this question as off-topic because it is not about bioinformatics but file manipulation. In any case it has been answered here.

@Priya The fact that you wanted a file with 86000 columns is not clear from your original post. You simply said you wanted to combine them into a single file. When posting questions, please take time to clearly and precisely describe the issue.

@DanielStandage Apologies for not being clear but I have not got my answer yet.

@Priya Please explain what is not clear from the thread I linked in my first comment.

@DanielStandage Sorry, I didn't click the link in your post. Thank you

@DanielStandage such questions are on topic since manipulating and parsing text files is a very large part of a bioinformatician's work. This question is really not clear, mind you, but it is on topic.

@Priya  tried to clarify the question a bit and I also suspected a typo in the code you tried. I hope the edits I have done are appropriate, if not please roll it back or [edit] it by your own :-).

The bug in your loop appears be "> filename" per column, this over-rides the last iteration resulting in a file with a single column, which has been over-written 85999 times. Thus the file output you observe should be the single column present in the last file of your array. If you replaced "> filename" with ">> filename" within the loop you would then append to the file a new column with each iteration. If the loop fell over at some point you could then simply use one of the outputs of wc filename to assess how many columns you had (you might have to transponse the data for wc to work (lots of scripts for this in Perl), although there will almost certainly be an easy solution for counting columns in awk.
Personally, I would use a dataframe for this operation (possibly R, definately Python pandas), because 86000 columns is a lot of columns in any other walk of life, even if it is just post-processing. Generally, if you just want bash, awk is a better solution in this case (I don't code with it, but some here certainly do). However, I do get the idea that paste will work if each column is an identical length.
BTW paste must be lower case, you've got Paste in the first line of your code, but its okay in the second bit of code.
BTW BTW, I did not know shell would take *.txt directly as an array, with Perl/Python you need to "glob" it and assign it as an array/list.


text:
What is outliers and how to determine outliers in the following phylogenetic tree?

I an confused with the term 'outlier' and 'outgroup' in phylogeny. Could someone please explain with the example phylogenetic tree above which is the correct term? Thank you very much !!!

Welcome to the group.  Can you explain your question more.  What do you need explained about the tree?

An "outlier" is a term used in general statistics to infer that a minority of the data set deviates substantially from the majority. An example would be the number of radical substitutions present in a house-keeping gene. An outgroup is described below and singularly refers to phylogeny.

Distance- and ML-based algorithms using reversible models can't find the root of trees. A classical method to root a tree is to use an outgroup (not outlier), which is a species/sequence known to directly descend from the root of the tree. In your case, it is relatively easy: add a chicken or fish ortholog to your dataset and put the root on the chicken/fish branch.
Outgroup works well for single-copy genes used in species tree construction (e.g. ADH). However, it doesn't always work. The culprit is the "known" part. Say a gene has two copies A and B in vertebrates. Copy A was lost in rodent and copy B lost in primate. If you choose a chicken A gene as the outgroup, the correct tree should be ((primary-A,chicken-A),rodent-B). Without knowing the true history, you may forcefully put the root at the chicken-A branch and build a wrong tree ((primate-A,rodent-B),chicken-A).
There are a few other tree rooting methods. An easy approach is to put the root at the longest branch in the tree, assuming the presence of molecular clock. When the species tree is known, you can root a gene tree by minimizing the number gene duplication/loss events in the history. I generally prefer the latter approach when the relevant information is available.

How to differentiate between outgroup and outlier? Thank you very much :)
@user172818

@ZhengKeongNg In phylogenetics, "outlier" is not clearly defined. Avoid that word.

you may use the OD-seq algorithm to filter outliers in a multiple alignment.
https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-015-0702-1


text:
Does the ".full.aln" file produced by snippy-core contain all bases of my input sequences aligned to the reference genome?
I have a number of sequences and a reference genome.
I used snippy to align each individual sequence with the reference genome. I then used snippy-core on the *.aln files from snippy which produced a new *.aln file, and a *.full.aln.
It is my understanding that the new *.aln file contains a core SNP alignment with only bases shared by all of my sequences, and that the *.full.aln file contains the core genome alignment with all bases of each sequence aligned to the reference genome - regardless of variation between sequences. Is my understanding correct? If not where am I going wrong?
Further, the program snp-sites uses both the new *.aln and *.full.aln file to produce a file which can be run through FastTree to produce a phylogenetic tree. Does this just process the files to give the phylogenetic distance between each of my sequences? I don't really understand what this step achieves at all...

snippy-core is a tool that processes the snippy_outdir_{1..n} folders and produces the following files:

core.aln – a fasta format file that contains only the polymorphic sites (of {'A', 'C', 'T' and 'G'}) in your input set, whose length tends to decrease asymptotically as more isolates are added
core.full.aln – a fasta format file that contains all the core sites and base calls at any site where the reference has a site. The reference is not padded but the sequences mapped to it may be padded with gaps and 'N's to match the length of the reference.  The alignment length does not change as more isolates are added.

snippy{-core} was written by Torsten Seemann et al.  snp-sites was written by Andrew Page et al.   snp-sites will reduce any alignment to the core polymorphic sites (columns) and will still output a fasta.  Therefore the output of a core.aln processed through snp-sites should be the identical the input.  All fasta alignments can then be processed downstream using your tree inference package of choice (e.g., FastTree, IQ-Tree, MrBayes, BEAST ...).  If you want a distance matrix output from a fasta alignment, use something like snp-dists, also written by Torsten Seemann, or the APE package in R.  Be sure to consider how gaps and Ns are treated during this step.

Thank you so much for this. Was breaking my head for a long time. Its working now.


text:
Prediction of prokaryotic origins of replication (ORI)
I want to predict origins of replication (ORI) on hundreds of prokaryotic genomes. The most straight-forward solution would be to use most commonly used tool, Ori-Finder.
It uses integrated gene prediction, analysis of base composition asymmetry, distribution of DnaA boxes, occurrence of genes frequently close to oriC regions and phylogenetic relationships to predict ORIs.
However, I came across a paper by Parikh et al which describes a single measure (correlated entropy measure / CEM) with which they achieved 100% ORI-identification on 500 genomes.

I like the CEM approach more than a close-source Windows-only black-box.
Parikh tested their algorithm on NCBI genomes. How are ORI annotated on these reference genomes? I couldn't find such an annotation. Does the sequence simply begin at the ORI in leading strand direction by convention? Were the ORI of the reference genomes identified by experiment?

Hi and welcome to the site! I removed your second question since that is not something we can give a definitive answer to and this site only deals with questions that can be answered without relying on opinion. Sorry about that!

I have not seen the papers but I know the technique SVM. This is a machine learning approach its pretty good, but has been around this the 1980s.
Its quite simple known ORIs are given with the associated genetic trait, e.g. GC skew, SVM will then regress the GC skew of a known ORI against the GC content of the genome, there will be a training, parameter optimization and test split in the data, typically 80:20:20. De novo data will then be feed into the trained algorithm, in the same way the test data was, and a prediction will be generated. There are weaknesses, the authors have used only one ML technique, typically numerous will be used e.g. SVM, random forest, lasso regression etc... In addition, the prediction accuracy for GC is low (75%), in addition 500 genomes isn't great for ML, you need 1000s.


text:
Are phylogenetic tree construction algorithms any different than general clustering algorithms?
Are phylogenetic tree construction algorithms any different from clustering algorhithms? I suspect the answer is no.
Of course phylogenetic tree construction uses biological knowledge, e.g special distance metrics, but does it brings anything new to clustering algorithm (e.g hierarchical, neighbor-joining etc).

Yes, they are very different. If no one else responds I'll provide a full theoretical perspective about what is going on.

@Michael G. Eagerly waiting for your answer. :)

Most phylogenetics theory is based on the probability of reversion mutations and this drastically distinguishes it from cluster algorithms based on a simple distance matrix

Well, distance-based tree building algorithms are based on a simple distance matrix.

The goal of a phylogeny is to estimate the "expected" number of mutations between all taxa in the analysis and their hypothetical common ancestors. A cluster-analysis will only identify the "observed" mutations and "expected" and "observed" mutations can be majorly different due to the major artefact of reversion mutation. This is particularly true of nucleotide phylogenies.
The key difference between clustering algorithms based on a "uncorrected" distance matrix and phylogeny is the latter is based on an explicit model to accommodate reversion mutations. The real problem is that there are only 4 bases, so randomly there is 1/4 chance that a mutation at a given position will mutate back to the original, e.g. A->C->A. Observed mutational differences = 0, expected (real) mutational differences = 2. What phylogeny is concerned about is reconstructing that "2". The problem is significant because almost every gene has regions of rapid mutations and regions of low mutations.
The primary way of doing this is via the Jukes-Cantor correction and is foundational in all nucleotide trees expect "p-distances". If nucleotide divergence is lower than 75% then it is possible to estimate the expected number of mutations using the observed via the JC correction. Moreover, when combined with a method of estimating rate variation within a gene, (usually the discrete gamma distribution), JC correction is very effective at recovering the "true tree". This is because fast evolving homologous sites are grouped together - big correction via JC, slow evolving sites are grouped together - small correction via JC. Other approaches to improving the JC correction are via identifying the bias between purine to purine mutations and pyrimidine to pyrimidine and purine to pyrimidine mutations.
The importance of JC was demonstrated by simulation studies of 4 taxa ((a,b), (c,d)), if two strains evolved very quickly when there sister lineages evolve slowly a clustering algorithm will report the fast strains are a sister group, i.e. ((b, c),(a, d)). If the JC method is implemented via maximum likelihood (or Bayesian) then it correctly recovers the true tree ((a,b), (c,d)). The artefact is known as long branch attraction.
Clustering algorithms based on a distance matrix that implements the JC correction, tend to perform poorly for long branch-attraction artefacts. This doesn't mean the correction is useless, just not particularly powerful. The issue is that distance matrix methods don't "stick" to the model and the clustering will introduce a layer of inaccuracy. Normally presentation of a "corrected" distance matrix combined with a clustering algorithm will require bootstrapping (resampling with replacement) to assess whether a given cluster is supported. Parameterised distance matrices, using neighbour-joining clustering in in-conjunction with a bootstrap are considered okay.
The @user172818 mentioned parsimony and this method is considered less reliable because it can't implement a JC correction. IMO it is possible weighted parsimony could make a "come-back", but it would be really complicated to implement a biological weighting method and would require extensive, independent calculations.

A great question, though a little ambiguous. I don't know what "general clustering algorithms" refer to. For biological sequences, building a tree can be thought as a way of clustering. Anyway...
There are different tree building algorithms. Max parsimony (MP), max likelihood (ML) and bayesian algorithms directly take sequences as input. They are distinct from distance based clustering.
Then there is a class of distance based algorithms in phylogenetics. They start from an all-pair distance matrix and aim to find a tree that is most compatible with the matrix. Among them, UPGMA is largely hierarchical clustering. Neighbor-joining is somewhat similar to hierarchical clustering, but it builds unrooted trees. FastME uses a very different approach. Briefly, for each topology it tries to find the best branch lengths by weighted least square (for details, see its paper); the best topology minimizes total branch lengths. FastME finds the best topology via nearest neighbor interchange (NNI), closer to some ML algorithms (e.g. PhyML).

does it brings anything new at the level of clustering algorithm?

Personally, I think the distance based tree building methods are more sophisticated than naive hierarchical clustering, provided that the real relationship follows a tree topology.


text:
Snakemake: Cannot find first rule?
Back again with another snakemake query. This time I decided to port my read cleaning and alignment pipeline to snakemake. Repeating the steps from previous question and trying not to make any typo here, I'm now receiving this "error" while trying a dry run using snakemake snakemake -nr -s alignMake:
Building DAG of jobs...
MissingRuleException:
No rule to produce snakemake (if you use input functions make sure that they don't raise unexpected exceptions).

Below is my snakefile:
import os
shell.executable("/bin/bash")

workdir: '/proj/uppstore2018136/publication/jupyterNotebooks'

pathToPEReads = 'raw_data/reads/BL6/totalRNA'
pathToSMReads = 'raw_data/reads/BL6/smallRNA'

pathToHISAT2index = 'indexes/hisat2/reference'
pathToBOWTIE2index = 'indexes/bowtie2/reference'

pathToTemp = 'temp/'
pathToBAM = 'processed_data/BAM'

readsPE, = glob_wildcards('raw_data/reads/BL6/totalRNA/{sample}_1.fastq.gz')
readsSE, = glob_wildcards('raw_data/reads/BL6/smallRNA/{sample}.fq.gz')

samToolsThreads = 6

rule all:
    input:
        expand('processed_data/BAM/{sample}.bam', sample=readsPE),
        expand('processed_data/BAM/{sample}.bam', sample=readsSE)

rule trimgalorePE:
    input:
        fastq1 = pathToPEReads + '{sample}_1.fastq.gz',
        fastq2 = pathToPEReads + '{sample}_2.fastq.gz'
    output:
        output1 = pathToTemp + '{sample}_1_val_1.fq.gz',
        output2 = pathToTemp + '{sample}_2_val_2.fq.gz'
    shell:
        '''
        module load bioinfo-tools
        module load cutadapt
        module load TrimGalore
        trim_galore --quality 20 --paired --length 80 --illumina --stringency 2 --clip_R1 11 --clip_R2 11 -o {pathToTemp} --paired {input.fastq1} {input.fastq2}
        '''

rule trimgaloreSEandDust:
    input:
        fastq = pathToSMReads + '{sample}.fq.gz',
    output:
        output = pathToTemp + '{sample}_trimmed.fq.gz',
    shell:
        '''
        module load bioinfo-tools
        module load cutadapt
        module load TrimGalore
        trim_galore --quality 20 --length 16 --max_length 32 --small_rna --fastqc --stringency 2 -o {pathToTemp} {input.fastq}
        '''

rule align_hisat2:
    input:
        fastq1 = pathToTemp + '{sample}_1_val_1.fq.gz',
        fastq2 = pathToTemp + '{sample}_2_val_2.fq.gz'
    output:
        output = 'processed_data/BAM/{sample}.bam'
    threads: 10
    shell:
        '''
        module load bioinfo-tools
        module load HISAT2
        module load samtools
        hisat2 -q -p {threads} --no-discordant --no-mixed --rna-strandness RF --dta -x {pathToHISAT2index} -1 {input.fastq1} -2 {input.fastq2} | samtools sort -@ {samToolsThreads} -o {output}
        '''

rule align_bowtie2:
    input:
        fastq = pathToTemp + '{sample}_trimmed.fq.gz'
    output:
        output = 'processed_data/BAM/{sample}.bam'
    threads: 10
    shell:
        '''
        module load bioinfo-tools
        module load bowtie2
        module load samtools
        bowtie2 --no-unal --end-to-end -D 20 -R 3 -N 1 -L 20 -i S,1,0.50 -p {threads} -x {pathToBOWTIE2index} -U {fastq} | samtools sort -@ {samToolsThreads} -o {output}
        '''

rule clean:
    shell:
        '''
        rm temp/*
        '''

I believe it is due snakemake being unable to find the first rule ? Maybe there is something wrong with parsing ? I tested the glob_wildcards and expand function independently and they appear to work properly.
Edit: Here is the directory tree:
.
├── citations
├── doc
├── genome
│   ├── genome.fa
│   └── Mus_musculus.GRCm38.94.gtf
├── indexes
│   ├── bowtie2
│   │   ├── reference.1.bt2
│   │   ├── reference.2.bt2
│   │   ├── reference.3.bt2
│   │   ├── reference.4.bt2
│   │   ├── reference.rev.1.bt2
│   │   └── reference.rev.2.bt2
│   └── hisat2
│       ├── reference.1.ht2
│       ├── reference.2.ht2
│       ├── reference.3.ht2
│       ├── reference.4.ht2
│       ├── reference.5.ht2
│       ├── reference.6.ht2
│       ├── reference.7.ht2
│       └── reference.8.ht2
├── LICENSE
├── notebooks
├── plots
├── processed_data
│   └── BAM
├── raw_data
│   └── reads
│       └── BL6
│           ├── smallRNA
│           │   ├── smBL6_1.fq.gz
│           │   ├── smBL6_2.fq.gz
│           │   └── smBL6_3.fq.gz
│           └── totalRNA
│               ├── EKIxxxx012_1.fastq.gz
│               ├── EKIxxxx012_2.fastq.gz
│               ├── EKIxxxx013_1.fastq.gz
│               ├── EKIxxxx013_2.fastq.gz
│               ├── EKIxxxx014_1.fastq.gz
│               ├── EKIxxxx014_2.fastq.gz
│               ├── ESSxxxx077_1.fastq.gz
│               └── ESSxxxx077_2.fastq.gz
├── README.md
├── recipes
│   ├── alignMake
│   ├── cluster.json
│   ├── indexMake
│   ├── indexMakeScript.out
│   └── indexMakeScript.sh
├── results
├── scripts
├── slurmLogs
│   ├── slurm-6575541.out
│   ├── slurm-6575543.out
│   ├── slurm-6575815.out
│   └── slurm-6580867.out
├── source
└── temp

Loosely related, I think your rule clean is never executed unless your run snakemake as snakemake [opts] clean. If you want to clean up stuff on completion, maybe look at the onsuccess and similar handlers. Also consider the temp flag...

@dariober - I agree. I didn't knew about onsuccess. I will include such handlers in future scripts. But the reason why I didn't want to automatically use clean was to manually check if things happened properly in background.

you are starting snakemake the wrong way.
snakemake snakemake

means "Running the pipeline in Snakefile until you reach a rule named snakemake or a file name snakemake is created."
If your snakefile is called snakemake and this is the pipeline you want to start, the correct syntax would be:
snakemake -s snakemake

Thanks - another one of the stupid mistake - since dry runs are also automated; I didn't notice that the driver was using using snakemake snakemake - should have noticed it while posting the question. Although I discovered it earlier, I guess this question should still stay online.


text:
Nextflow: No signature of Method: chr() is applicable for argument types; unknown variable
I was trying to use plink in nextflow. My first process will split the input (single bed|bim|fam file) into 22 chromosomes (each with bed|bim|fam files) which would be used in subsequent processes. My code and error log is below. Any help?
params.dir='data/'
params.publish="output2/"
Channel
.fromFilePairs("${params.dir}*.{bed,bim,fam}",size:3)
.set {data}

Channel
.of(1..22)
.set {chr_names}

process step1a {
publishDir params.publish
input:
tuple val(sample_name), path(bfiles) from data
output:
tuple val("${prefix}"), path("${prefix}.step1a.{bed,bim,fam}") into step1a_results
tuple val("${prefix}"), path("${prefix}.step1b.{bed,bim,fam}") into step1a_out
tuple val("${prefix}"), path("${prefix}.step1b.bim") into step1a_bims
path("${prefix}.log") into step1a_logs
path("${prefix}.step1a.dupvar") into step1a_duplicates
path("${prefix}.step1a.log") into step1a_log
path("${prefix}.step1b.log") into step1_log
script:
prefix=chr"${chr_names}"
"""
plink --bfile "${sample_name}" --chr "${chr_names}" --out "${prefix}" --make-bed
plink --bfile "${prefix}" --list-duplicate-vars ids-only suppress-first --out "${prefix}.step1a"
plink --bfile "${prefix}" --exclude "${prefix}.step1a.dupvar" --make-bed --out "${prefix}.step1a"
plink --bfile "${prefix}.step1a" --snps-only just-acgt --make-bed --out "${prefix}.step1b"
"""
}
process step1b {
publishDir params.publish
input:
tuple val(sample_name), path(bfiles) from step1a_out
tuple val(bims), path("${prefix}.step1b.bim") from step1a_bims
output:
tuple val("${sample_name}"), path("${sample_name}.step1c.{bed,bim,fam}") into step1b_results
path("${sample_name}.step1c.log") into step1b_logs
path("${sample_name}.step1b.snplist.txt") into step1b_snplist
script:
"""
awk '{ if ((\\\$5=="T" && \\\$6=="A")||(\\\$5=="A" && \\\$6=="T")||(\\\$5=="C" && \\\$6=="G")||(\\\$5=="G" && \\\$6=="C")) print \\\$2, "ambig" ; else print \\\$2 ;}' "${prefix}.step1b.bim" | grep ambig > "${prefix}.step1b.snplist.txt"
plink --bfile "${prefix}.step1b --exclude "${prefix}.step1b.snplist.txt" --make-bed --out "${prefix}.step1c"
"""
}

The error message is:
N E X T F L O W  ~  version 20.10.0
Launching `t4.hf` [zen_picasso] - revision: b0ff8db413
[-        ] process > step1a -
[-        ] process > step1b -
Error executing process > 'step1a (1)'

Caused by:
[-        ] process > step1a -
[-        ] process > step1b -
Error executing process > 'step1a (1)'

Caused by:
  No signature of method: chr() is applicable for argument types: (org.codehaus.groovy.runtime.GStringImpl) values: [DataflowQueue(queue=[DataflowVariable(value=1), DataflowVariable(value=2), DataflowVariable(value=3), DataflowVariable(value=4), DataflowVariable(value=5), DataflowVariable(value=6), DataflowVariable(value=7), DataflowVariable(value=8), DataflowVariable(value=9), DataflowVariable(value=10), DataflowVariable(value=11), DataflowVariable(value=12), DataflowVariable(value=13), DataflowVariable(value=14), DataflowVariable(value=15), DataflowVariable(value=16), DataflowVariable(value=17), DataflowVariable(value=18), DataflowVariable(value=19), DataflowVariable(value=20), DataflowVariable(value=21), DataflowVariable(value=22), DataflowVariable(value=groovyx.gpars.dataflow.operator.PoisonPill@28c22659)])]
Possible solutions: run(), run(), any(), is(java.lang.Object), use([Ljava.lang.Object;), tap(groovy.lang.Closure)

Source block:
  prefix=chr"${chr_names}"
  """
  plink --bfile "${sample_name}" --chr "${chr_names}" --out "${prefix}" --make-bed
  plink --bfile "${prefix}" --list-duplicate-vars ids-only suppress-first --out "${prefix}.step1a"
  plink --bfile "${prefix}" --exclude "${prefix}.step1a.dupvar" --make-bed --out "${prefix}.step1a"
  plink --bfile "${prefix}.step1a" --snps-only just-acgt --make-bed --out "${prefix}.step1b"
  """

Tip: when you have fixed the problem you can continue the execution adding the option `-resume` to the run command line

If I use "prefix=chr"${chr_names}" " inside the code block (""") it shows: unknown variable 'prefix'

The documation re conditional scripts is useful here:

Process scripts can contain conditional statements by simply prefixing
the script block with the keyword script:. By doing that the
interpreter will evaluate all the following statements as a code block
that must return the script string to be executed.

But the statements don't necessarily need to be conditional. I.e. they can just be declarations like in your example. You get the error 'No signature of Method: chr()' because Nextflow interprets 'chr' as a method call when you use:
prefix = chr"${chr_name}"

If your intention is to simply add the 'chr' prefix to each chromosome name, you'd need:
prefix = "chr${chr_name}"

Notice that I changed your variable name. This is because in your code, you still need to feed in the channel of chromosome names. The easiest way to do this is to use the each qualifier:
process myprocess {

    input:
    tuple val(sample_name), path(bfiles) from plink_data
    each chr_name from chr_names

    script:
    prefix = "chr${chr_name}"

    """
    echo "My sample name is: ${sample_name}"
    echo "My prefixed chromosome name is: ${prefix}"
    """
}

Unless you need chromosome names without the 'chr' prefix, another solution might be to simply add the prefix to each one initially and simply avoid having to prefix them later:
Channel
    .of(1..22)
    .map { "chr${it}" }
    .set { chr_names }

It seems process 1 and 2 running in random order. Is there any way to specify that after finishing the first process, the second process will start. Because the outputs of the first process will be the input of the second process. In my case my second process shows error because it does not find an input, which is an output of the first process. Did I clarify the matter?

I suspect the problem arises because you've specified two input and output channels ('step1a_out' and 'step1a_bims'). It's not clear to me why you need the two channels, since the 'step1a_bims' looks like a subset of 'step1a_out'? Please see understanding how multiple input channels work

Ohh! Thank you. Just one thing, How could I declare prefix of the previous process output as an input of the current process. I have chr(1-22).step1b.bed|bim|fam as output of process step1a; Now I want to declare only chr(1-22) as input in the step1b process. I want different suffix inside process step1b. I am not sure if I could clarify the issue.

Does my answer to your most recent question answer your question?

Yes. Thank you. Still wondering how could I declare prefix (characters before first .) as input in the next process which is the output of the previous process.

No worries at all @zillurrahman. If I understand correctly, the easy way would be to include it in your output tuple as an extra value for your prefix. But if it's already part of (i.e. the prefix) of another value then you could again use the map operator to get it back. Perhaps something like: map { some_name, some_thing -> tuple( some_name.substring(0, some_name.indexOf('.')), some_name, some_thing ) } ? A concrete example would be great if the above doesn't help.

Thanks a lot @Steve. I tried different ways but failed, maybe because I am new in nextflow. My files are chr1-22.step1b.bed|bim|fam in out_ch1. I need a way to find chr1 to chr22 as prefix for the next input so I can modify the suffixes. for example input chr1-22.step1b.bed|bim|fam --some command --out chr1-22.step1c.bed|bim|fam into out_ch3.
Also how can I use subset of the outputs channel. From above example I need only chr1-22.step1b.bim as a $bim --some command to run -out chr1-22.step1b.snp.txt. So far I tried "filter" but failed. Thanks again for helping all the way.

@Steve, I tried to filter as: tuple val(bim), path(bims) from step1a_bim.filter{bim -> bims=~ /*bim/} But I am getting error unexpected input: Also I am still not getting the prefix correctly. Perhaps any example would be helpful. Do you know any training module to run plink, gwas in nextflow?

Sorry @zillurrahman - I'm not sure I can get you a useful answer without a minimal example of exactly what you're trying to do. On file-like objects you can call getSimpleName() to get the prefix, see check file attributes. The problem is that this won't work on a tuple of files. I guess I'm confused as to why you can't just include the chromosome name in the output tuple of the first process? This way you'd avoid having to try and find the prefix in the first place.

Using my answer above as an example, your output declaration could look like: tuple val(sample_name), val(chr_name), path("${sample_name}.${chr_name}.step1b.{bed,bim,fam}") into step1a_out. Also, I'm not aware of any good plink/gwas pipelines but you could try having a look here: https://github.com/nextflow-io/awesome-nextflow. HTH.

Thanks a lot. In that case how can I declare inputs in the next process? If I use step1a_out two time (first for bfiles, and then for bim) I got error Channel 'step1a_out' has been used twice as an input by process 'step1b'

If I understand correctly, the input declaration for the next process will be basically the same as the first processes' output declaration. If that's the case, you'd need: tuple val(sample_name), val(chr_name), path(bfiles). Since the bfiles tuple will contain the bim  file as the second element, you should be able to access it directly with ${bfiles[1]}. Alternatively, in your script block, you could unpack the bfiles: (bed, bim, fam) = bfiles, then just use ${bim} as your variable. If these variables don't need to be global, consider instead: def (bed, bim, fam) = bfiles. HTH.

Thanks a lot @Steve . It works. Just one thing. How could I specify output directory of my results? I tried params.publish='output2' then in each process publishDir params.publish. I got only link of the output files in output2. How could I get actual result files in my desired directory? Again thanks a lot for helping me all the way.

Ohh! I got it. I just need to mention mode: 'copy' in publishDir declaration. Thanks!


text:
nextflow: Filter outputs of a process
How could I filter outputs of a process in the input of the next process? Filtering works fine in channel, but if I try to filter outputs I got compilation error. I tried as follows:
Channel
.fromFilePairs("${params.dir}*.{bed,bim,fam}",size:3)
.set {data}
    process step1a {
    input:
    tuple val(sample_name), path(bfiles) from data
    output:
    tuple val("${prefix}"), path("${prefix}.{bed,bim,fam}") into (step1a_out,step1a_bim)
    script:
    """
    plink --bfile "${sample_name}" --chr 1-22 --out "${prefix}.step1b" --make-bed
    """
    }
    process step1b {
    input:
    tuple val(bim), path(bims) from step1a_bim.filter{bim, files -> name =~/*bim/}
    tuple val(sample_name), path(bfiles) from step1a_out
    output:
    path("${prefix}.step1b.snplist.txt")
    tuple val("${sample_name}"), path("${sample_name}.step1c.{bed,bim,fam}") into step1b_results
    script:
    """
    awk '{ if ((\\\$5=="T" && \\\$6=="A")||(\\\$5=="A" && \\\$6=="T")||(\\\$5=="C" && \\\$6=="G")||(\\\$5=="G" && \\\$6=="C")) print \\\$2, "ambig" ; else print \\\$2 ;}' "${prefix}.step1b.bim" | grep ambig > "${prefix}.step1b.snplist.txt"
    plink --bfile "${prefix}.step1b" --exclude "${prefix}.step1b.snplist.txt" --make-bed --out "${prefix}.step1c"
"""

It gives me
Compilation error
-cause: Unexpected input: '{'
process step1b {

Step1a gives me chr1.step1b.bed|bim|fam, chr2.step1b.bed|bim|fam.......chr22.step1b.bed|bim|fam (22 plink files)
I need to access prefix of those outputs to run plink and other tools to generate chr1.step1c.bed|bim|fam ...chr22.step1c.bed|bim|fam (another 22 plink files). run awk to generate chr1.step1b.snplist.txt .... chr22.step1b.snplist.txt (22 txt files).
Any help to filter the inputs so that I could use only *bim files from step1a_bim? Also how to select prefix of files in the input of processes that are from step1a_out to run plink? Does it make sense now?
Best Reagrds
Zillur

In case of a compilation error it is more useful if you post your actual code, properly formatted, and your error message.

@Pallie, I edited my question above. Does it clarify my issues on outputs and inputs?

Instead of outputting all files and filtering for only bim, output only bims:
output:
tuple val("${prefix}"), path("${prefix}.{bed,bim,fam}") into (step1a_out,step1a_bim)

Do
output:
tuple val("${prefix}"), path("${prefix}.{bed,bim,fam}") into step1a_out
tuple val("${prefix}"), path("${prefix}.bim") into step1a_bim

So you can skip your complicated filter:
process step1b {
    input:
    tuple val(bim), path(bimfile) from step1a_bim
    output:
    path("${prefix}.step1b.snplist.txt")
tuple val("${sample_name}"), path("${sample_name}.step1c.{bed,bim,fam}") into
step1b_results
    script:
    """
    awk '{ if ((\\\\\$5=="T" && \\\\\$6=="A")||(\\\\\$5=="A" && \\\\\$6=="T")||(\\\\\$5=="C" && \\\\\$6=="G")||(\\\\\$5=="G" && \\\\\$6=="C")) print \\\\\$2, "ambig" ; else print \\\\\$2 ;}' "$bimfile" | grep ambig > "${prefix}.step1b.snplist.txt"
    plink --bfile "${prefix}.step1b" --exclude "${prefix}.step1b.snplist.txt" --make-bed --out "${prefix}.step1c"

"""

```

Thanks @Pallie . I tried your method. In that case, I got Process 'step1b' input file name collision -- There are multiple input files for each of the following file names: chr5.step1b.bim. Any suggestion?

The step1a_bim is a subset of step1a_out which produce bed|bim|fam files.

I could divide process step1b into two processes and could use step1a_out. But how could I use only *.bim files as input? Because step1a_out produce all *.bim, *.bed and *.fam files. For one process I need all three (bed|bim|fam) and for another I need only bim files. I am not sure how to subset the elements of a channel.


text:
parse csv/text file to create a channel with array in nextflow
[![example csv file][1]][1]
First I want to start by parsing and doing some operation on a text/csv file. How could I read a csv/text file in nextflow?
Suppose, I have above csv file with 3 columns. Lets say I want to create a channel which emits an array from minimum of column 2 (7.9e5) to maximum of column 2 (8.8e5) with an increment of 1e4. How could I create the following channel and later access it in the script? Any help? The
array=([1]="7.9e5 8.0e5" [2]="8.0e5 8.1e5" [3]="8.1e5 8.2e5" [4]="8.2e5 8.3e5" [5]="8.3e5 8.4e5" [6]="8.4e5 8.5e5" [7]="8.5e5 8.6e5" [8]="8.6e5 8.7e5" [9]="8.7e5 8.8e5" [10]="8.8e5 8.9e5")

I need to do some operations with this array in the script which will generate multiple files for each element in the array. These multiple files need to be merged in a single file.
Sorry, I could not provide any nextflow script that I tried, because I not sure how create an array with splitCsv
I have tried in python:
df1=pd.read_csv('data1/omni-ec-eur-fin-pass-qc.bim',sep='\t',header=None)
arr1=np.arange(df1.iloc[:,3].min()-1,df1.iloc[:,3].max()+1,3000000)
diffs = []
for x in zip(arr1[0::],arr1[1::]):
    diffs.append(x)
arr2=np.arange(1,len(diffs)+1,1)
dict(zip(arr2,diffs))

Which gave me:
{1: (12343, 3012343),
 2: (3012343, 6012343),
 3: (6012343, 9012343),
 4: (9012343, 12012343),...

But I need something like this:
([1]="12343 3012343" [2]="3012343 6012343" [3]="6012343   9012343" [4]="9012343 12012343"....)

It could be arrays (22 arrays for 22 chromosome) or a csv file containing arrays for each chromosome which could be used as an input channel. There are 22 bim files (one for each chromosome) stored in the output channel which would be used to create the above mentioned array.
The original data inside the bim looks like this:
    0   1   2   3   4   5
0   1   rs11240777  0   798959  A   G
1   1   rs4475691   0   846808  A   G
2   1   rs7537756   0   854250  G   A
3   1   rs13302982  0   861808  A   G
4   1   rs1110052   0   873558  C   A

The 4th column would be used to create the arrays.
[1]: https://i.sstatic.net/cxqyl.png

PLease [edit] your question and give us an example file and minimal nextflow script to work with.

Is there any book to learn nextflow operations? I was trying to find any book to learn nextflow but failed. At least if there is any example of creating array in nextflow would be great.

Have you read the nextflow documentation? I don't know much about it, but the whole idea is to use it for control flow; it runs other programs and they can have arrays. But next flow isn't a programing language per se. You write your script in Python or peel or whatever and then run the pipeline if scripts with nextflow.

@zillurrahman Do you really need an array? Or do you mean a channel? Your last paragraph suggests the latter would be a better fit. Please edit your question to clarify exactly what is needed.

@zillurrahman Generally, it's better to off-load all work to one or more processes. Writing a process to create a range of values, like in your example, is trivial in Perl/Python. You can use your favorite language in a process' script block. Please see: https://www.nextflow.io/docs/latest/process.html#scripts-a-la-carte

@terdon It's an extension of the Groovy programming language :-)

@steve sure, but it isn't designed to be used as a regular programing language but as a process control system to build and manage pipelines. It seems far more reasonable to manage the data in another script and then have nextflow manage that script instead of trying to use nextflow as your scripting language directly. Which is exacy what you suggested too.

Thanks @Steve and @terdon, I edited the above question. Hope, it would clarify my necessity.  However, I am also trying to create such array in bash . But, I think probably a nextflow channel would be more appropriate in this case.

@terdon wasn't disagreeing with you on that...

@zillurrahman Is there some reason for choosing Bash over Python for example? Bash isn't a good choice given your requirements.

@Steve No reason to choose bash. I tried in python and got an array. np.arange(df1.iloc[:,3].min()-1,df1.iloc[:,3].max()+1,3000000) which gives me: array([    12343,   3012343,   6012343,   9012343,  12012343, .... . How could I get something like: ([1]="12343 3012343" [2]="3012343 6012343" [3]="6012343   9012343" [4]="9012343 12012343"....)

So far, I have done: diffs = [] for x in zip(arr1[0::],arr1[1::]):     diffs.append(x) diffs=np.array(diffs) arr2=np.arange(1,len(diffs)+1,1) arr3=np.column_stack((arr2,diffs)) arr4=arr3.tolist() which gave me [[1, 12343, 3012343],  [2, 3012343, 6012343],  [3, 6012343, 9012343], . Should I convert it to string and use regex to get ([1]="12343 3012343" [2]="3012343 6012343" [3]="6012343   9012343" [4]="9012343 12012343"....) or there are better alternatives?

@zillurrahman Do you really need to create the array string? If all you eventually want to do is create a channel of tuples, your first process should write a simple CSV which can then be split using the splitCsv operator.

@zillurrahman Also, your comments above should be edits to your question. Small snippets of unformatted code isn't especially helpful. Please include a small minimal example and the code that reproduces the issue.

Thanks @Steve. I have edited my above question. Hope it makes sense now.

@zillurrahman Like terdon mentioned earlier, please include some example input (i.e. cat data1/omni-ec-eur-fin-pass-qc.bim). Please format it in a code block.

@Steve . The first few rows of the input data is added to the above question.

Hi, any help/suggestion for the above problem? Did I clarify enough the issue?


text:
Merge multiple channels each emits multiple files in nextflow
I need to merge multiple channels each emits multiple file. So far I did something as below:
process prephasing {
publishDir params.out, mode:'copy'
input:
tuple val(prefix), path(bfiles) from plink_data
each chrom from chroms
path hrc_sites_file
output:
tuple val("${chrom}"),path("chr${chrom}*.{bed,bim,fam}") into qc1_chr
"""
PRE_PASHING_QC2.sh "${hrc_sites_file}" "${prefix}" "${chrom}"
"""
}
channel.fromPath('/some/path/genetic_map_chr*_combined_b37.txt').into {shape_map;i
mp_map}
channel.fromPath('/some/path/*hrc.r1-1.ega.grch*.hap.gz').into {shape_ref;
imp_ref}
channel.fromPath('/some/path/*hrc.r1-1.ega.grch37.chr*.legend.gz').into {s
hape_legend;imp_legend}
channel.fromPath('/some/path/*hrc.r1-1.ega.grch37.chr*.samples').into {sha
pe_sample;imp_sample}
qc1_chr.filter{name, bfiles -> bfiles.collect {it.getBaseName()}.every {it.endsWith("step10")}}.into {shape_in;imp_in}

shape_para=shape_in.merge(shape_map).merge(shape_ref).merge(shape_legend).merge(shape_sample)

process phasing {
publishDir params.out, mode:'copy'
input:
each file(penta) from shape_para
output:
path("*step10*") into phased
script:
def (in, map, ref, legend, sample) = penta
"""
shapeit2 -T 32 --input-bed "${in}" --input-map "${map}" --input-ref "${ref}" "${legend}" "${sample}" --duohmm -O "${in}.phased" --seed 54321
"""
}

The "prephasing" process works fine and produce expected output. But the "phasing" process does not produce any output nor any error. Any help?
I have tried according to @Steve as below:
channel.fromFilePairs('/cluster/projects/p697/projects/moba_qc_imputation/resources/1000Genomes/genetic_map_chr*_combined_b37.txt',size:1).map {group_key, file_list -> tuple(group_key.replaceFirst(/^genetic_map_/,""),file_list.first())}.set {genetic_map}

channel.fromFilePairs('/cluster/projects/p697/projects/moba_qc_imputation/resources/HRC/haplegendsample/*hrc.r1-1.ega.grch37.chr*.haplotypes.{hap.gz,legend.gz,samples}',size:3).map{group_key,file_list -> tuple(group_key.replaceFirst(/^*HRC\.r1-1\.ega\.grch37\./,""),file_list)}.set {ref_panel}

qc1_chr.filter{chrom, bfiles -> bfiles.collect {it.getBaseName()}.every {it.endsWith("step10")}}.join(ref_panel).join(genetic_map).into {shape_in;imp_in}

process phasing {
publishDir params.out, mode:'copy'
input:
tuple val(chrom), path(input_bed), path(intput_ref), path(input_map) from shape_in
output:
tuple val("${chrom}"),path("${chrom}.{phased.haps,phased.sample,log}") into phased
script:
def (bed, bim, fam) = input_bed
def (haplotype, legend, sample) = input_ref
"""
shapeit2 --thread "${task.cpus}" --input-bed "${bed}" "${bim}" "${fam}" --input-map "${input_map}" --input-ref "${haplotype}" "${legend}" "${sample}" --duohmm --output-max "${chrom}.phased" --seed 54321
"""
}

But getting no output so far, (attached image)
I have tried to check individual channel separately using channel.view(), identified the problems, modified the code and checked again. Using the code below I can see individual channel emits 22 lines [chr1-22,[file_list]]. After merging I can see [chr1-22,[joined_file_list].
channel
    .fromFilePairs( "output3/*.{bed,bim,fam}", size:3 )
    .set { plink_data }
plink_data.filter{chrom, bfiles -> bfiles.collect {it.getBaseName()}.every {it.endsWith("step10")}}.map{group_key,file_list -> tuple (group_key.
replaceFirst(/^*\.step10/,""), file_list)}.set {in_pl}

channel.fromFilePairs('/cluster/projects/p697/projects/moba_qc_imputation/resources/1000Genomes/genetic_map_chr*_combined_b37.txt',size:1).map {
group_key, file_list -> tuple(group_key.replaceFirst(/^*genetic_map_/,""),file_list.first())}.set {genetic_map}
//genetic_map.view()

channel.fromFilePairs('/cluster/projects/p697/projects/moba_qc_imputation/resources/HRC/haplegendsample/*hrc.r1-1.ega.grch37.chr*.haplotypes.{ha
p.gz,legend.gz,samples}',size:3).map {group_key,file_list -> tuple(group_key.replaceFirst(/^.*.hrc\.r1-1\.ega\.grch37\./,""),file_list)}.set {re
f_panel}
//ref_panel.view()

in_pl.join(ref_panel).join(genetic_map).into {shape_in;imp_in}
process phasing {
//publishDir params.out, mode:'copy'
input:
tuple val(chrom), path(input_bed), path(intput_ref), path(input_map) from shape_in
output:
tuple val("${chrom}"),path("${chrom}.{phased.haps,phased.sample,log}") in
to phased
script:
def (bed, bim, fam) = input_bed
def (haplotype, legend, sample) = input_ref
"""
//shapeit="/cluster/projects/p697/projects/moba_qc_imputation/software/shap
eit"
shapeit --thread "${task.cpus}" --input-bed "${bed}" "${bim}" "${fam}" --input-map "${input_map}" --input-ref "${haplotype}" "${legend}" "${samp
le}" --duohmm --output-max "${chrom}.phased" --seed 54321
"""
}

But now I get error:
No such variable: chrom

But I decalred "chrom" in the input
The input bed|bim|fam files from channel from FilePairs works but from a output channel of previous process does not. If I run phasing solo and all the inputs from channel from FilePairs it works. But if I run prephasing and phasing together, the prephasing process works but for phasing process not output and no error; Not sure what I am doing wrong.
hrc_sites_file = file(params.hrc_sites_file)

Channel.fromFilePairs( "${params.input_dir}/*.{bed,bim,fam}", size:3 ).set { plink_data }

Channel.of(1..22).set { chroms }
params.out='output/'

process prephasing {
publishDir params.out, mode:'copy'
input:
tuple val(prefix), path(bfiles) from plink_data
each chrom from chroms
path hrc_sites_file
output:
tuple val("${chrom}"),path("chr${chrom}*.{bed,bim,fam}") into qc1_chr
tuple val("${chrom}"),path("chr${chrom}*log") into qc1_log
tuple val("${chrom}"),path("*-chr${chrom}.step4-HRC.txt") into qc1_hrc
path("chr${chrom}.step1b.snplist.txt") into qc1_snps
path("chr${chrom}.step1a.dupvar") into qc1_dups
"""
PRE_PASHING_QC2.sh "${hrc_sites_file}" "${prefix}" "${chrom}"
"""
}
channel.fromFilePairs('/cluster/projects/p697/projects/moba_qc_imputation/resources/1000Genomes/genetic_map_chr*_combined_b37.txt',size:1).map {group_key, file_list -> tuple(group_key.replaceFirst(/^*gen
etic_map_/,""),file_list.first())}.set {genetic_map}

channel.fromFilePairs('/cluster/projects/p697/projects/moba_qc_imputation/resources/HRC/haplegendsample/*hrc.r1-1.ega.grch37.chr*.haplotypes.{hap.gz,legend.gz,samples}',size:3).map{group_key,file_list ->
 tuple(group_key.replaceFirst(/^.*.hrc\.r1-1\.ega\.grch37\./,""),file_list)}.set {ref_panel}

qc1_chr.filter{chrom, bfiles -> bfiles.collect {it.getBaseName()}.every {it.endsWith("step10")}}.map{group_key,file_list -> tuple (group_key.replaceFirst(/^*\.step10/,""), file_list)}.join(ref_panel).joi
n(genetic_map).into {shape_in;imp_in}
imp_in.view()
process phasing {
publishDir params.out, mode:'copy'
input:
tuple val(chrom), path(input_bed), path(input_ref), path(input_map) from shape_in
output:
tuple val("${chrom}"),path("${chrom}.*") into phased
script:
def (bed, bim, fam) = input_bed
def (haplotype, legend, sample) = input_ref
"""
shapeit --thread "${task.cpus}" --input-bed "${bed}" "${bim}" "${fam}" --input-map "${input_map}" --input-ref "${haplotype}" "${legend}" "${sample}" --duohmm --output-max "${chrom}.step10" --seed 54321
"""
}

The merge operator is deprecated and will be removed from future Nextflow releases. You'll need instead the join operator to join by chromosome. You can use the fromFilePairs operator to quickly get the group key, formatting it as required:
channel
    .fromFilePairs("/some/path/HRC.r1-1.EGA.GRCh37.chr*.{hap.gz,legend.gz,samples}", size: 3)
    .map { group_key, file_list ->
        tuple( group_key.replaceFirst(/^HRC\.r1-1\.EGA\.GRCh37\./, ""), file_list )
    }
    .set { ref_panel_files }

channel
    .fromFilePairs("/some/path/genetic_map_chr*_combined_b37.txt", size: 1)
    .map { group_key, file_list ->
        tuple( group_key.replaceFirst(/^genetic_map_/, ""), file_list.first() )
    }
    .set { genetic_map_files }

qc1_chr
    .map { chrom, bfiles ->
        tuple( chrom, bfiles.findAll { it.baseName.endsWith("step10") } )
    }
    .join( ref_panel_files )
    .join( genetic_map_files )
    .set { phasing_inputs }

process phasing {

    input:
    tuple val(chrom), path(input_bed), path(input_ref), path(input_map) from phasing_inputs

    script:
    def (bed, bim, fam) = input_bed
    def (haplotypes, legend, sample) = input_ref

    """
    shapeit2 \\
        --thread "${task.cpus}" \\
        --input-bed "${bed}" "${bim}" "${fam}" \\
        --input-ref "${haplotypes}" "${legend}" "${sample}" \\
        --input-map "${input_map}" \\
        --duohmm \\
        --output-max "${chrom}.phased" \\
        --seed 54321
    """
}

Thanks @Steve. Only one thing, the prefixes of hap.gz, legend.gz and samples are not same for each chromosome. That's why in my code there is a * at the beginning and another * at the chromosome label.

Hi @Steve I edited my question adding picture of running script and code snippet.

@zillurrahman: I think, though, you understand the general gist of the answer - that is, to create a tuple with a key as the first element which can then be used to join (multiple) channels. The above is just an example that uses the chromosome name as the join key. You should not expect a drop-in replacement for your code. I don't know what files your glob patterns match, and still don't. Your patterns could match multiple files, but I assumed you had one (hap, legend, sample) per chromosome. If the reference panel files have different prefixes, you may prefer separate channels for these.

Sorry for the inconveniences. I rechecked the code, there were some type. Modified them and can see joined items. But getting a new error. No such variable: chrom But I declared the chrom variable in the input.

The new code snippet was added to my above question.

@zillurrahman I can't reproduce the error you get using your new code snippet. It should have failed earlier because of this typo: path(intput_ref) should be path(input_ref). Do you know what line the error is being produced on?

Thanks @Steve. The error is occurring at line 20 where the output is declared.

In the .nextflow.log file, it is showing that the error is occurring at line 20. And line 20 is tuple val("${chrom}"),path("${chrom}.{phased.haps,phased.sample,log}") into phased where the output was declared.

@Steve you are right. There was a type in my code. The error is Unknown variable 'input_ref. The typo is fixed and the code is running now! Thanks a lot!

If I run phasing separately, it works. But in combination with prephasing  it does not work. I mean bed|bim|fam files from channel from FilePairs works but from the output channel of previous process does not. Even I can not see the joined channel.view() in .command.out in the work directory. No error also. not sure what goes wrong. I added the updated code snippet in the above question.

@zillurrahman: I think your ouput declaration for qc1_chr needs to be: tuple val("chr${chrom}"),path("chr${chrom}*.{bed,bim,fam}") into qc1_chr. HTH.

@Steve Thanks. I updated the output declaration of process prephasing. Still the process phasing does not work nor any error message. In the last screenshot of the above question, you can see that there is a work directory (at the left) assigned for prephasing process but nothing for phasing process. Not sure what I am doing wrong. The same phasing process works if I run it in a separate script where input bed|bim|fam  are assigned from FilePairs.

@zillurrahman: I suspect one of your input channels to join on isn't working as expected. Could you please post the output of qc1_chr.view(), ref_panel.view() and  genetic_map.view()? Thanks!

Thanks @Steve. ref_panel.view() gives [chr,[respective hap.gz,legend.gz,sample.gz]], genetic_map.view() gives [chr,[respective maps]]. The problem with filtering qc1_chr. qc1_chr.filter{chrom, bfiles -> bfiles.collect {it.getBaseName()}.every {it.endsWith("step10")}}.view()  shows nothing. Interestingly, same filtering method works when the bed|bim|fam files were from FilePairs. I also tried qc1_chr.filter { it[1].name.contains('step10') }.view() but nothing. qc1_chr.view() gives [chr1-22,[respective all plink files]]

My mistake, I think you'll need to replace the filter statement with .map { chrom, bfiles -> tuple( chrom, bfiles.findAll { it.baseName.endsWith("step10") } ) }. No filtering is required, we just need to select out the required step10 files.

Thank you very much @Steve. It works now. as there are multiple patterns, I used channel.fromFilePairs('/cluster/projects/p697/projects/moba_qc_imputation/resources/HRC/haplegendsample/*hrc.r1-1.ega.grch37.chr*.{hap.gz,legend.gz,samples}',size:3).map {group_key,file_list -> tuple(gro up_key.replaceAll(/^.*.hrc\.r1-1\.ega\.grch37\.|\.haplotypes\.noibd|\.haplotypes/,""),file_list)}.set {ref_panel}. Do you think replaceAll is ok instead of replaceFirst . I got expected output using using replaceAll. Could I use multiple pattern using replaceFirst?

@zillurrahman: I don't think you can since the first looks like a prefix and the last two look like suffixes. I think you'd need two replaceFirst calls. I'd be tempted to just use the first string that looks like a chromosome name: group_key.find(~/chr\d{1,2}/). HTH.

Thanks @Steve I just tested the find option and it works. Thanks a lot! It seems fewer code and easier to write!

Hi @Steve, Wondering if you get a chance to look at https://bioinformatics.stackexchange.com/questions/15698/parse-csv-text-file-to-create-arrays-in-nextflow . I just need a little help to create an array from a csv file in nextflow.


text:
Metagenomics pipeline recommendations for short-read data
de novo metagenomics on viral NGS data is a hot-topic. On this site alone at least 4 specific algorithms have been used to identify multi-strain/multi-species for a given data set, however these do not necessarily involve short-read data, and are listed below.
Questions
What approach is recommended for traditional Illumina multiple infection data both for:

for intraspecific (intra-species)  data; and
for different species co-infecting the patient, but within the same virus family?

Considerations:

Benchmarking using a common standard (for example)
Error rates, not miscalling a SNP but producing mosaic artefacts (an artificially mixed genome)
Suitability of short-read data for intra-species metagenomics.

Programs discussed on Bioinformatics SE:
SAVAGE https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5411778/ Baaijens et al. De novo assembly of viral quasispecies using overlap graphs Genome Research 2017 May; 27: 835–848.
S-aligner (@juanjo75es)
https://www.biorxiv.org/content/biorxiv/early/2021/02/02/2021.02.02.429443.full.pdf Juanjo Bermúdez s-aligner: a greedy algorithm for non-greedy de novo genome assembly 2021 bioRχiv February: 1-7
Rey http://denovoassembler.sourceforge.net Sébastien Boisvert et al Genome Biology. 2012 13:R122, December
@LRJoshi here used the Broad Viral NGS pipeline for de novo assembly.
https://viral-ngs.readthedocs.io/en/latest/assembly.html

These are very small genomes < 30 kb. Quasi-species is a very real phenomenon (loads of mutations), however mixed infections are also well-known. In a mixed infection, producing one genome which is a mix of two genetically distinct viruses is a key concern.

Under my experience, the best ones in completeness and N50 are s-aligner, Megahit and rnaSPAdes, but I admit I didn't test the other ones for different reasons. Broad indeed seems to use internally SPAdes and Trinity. I compared the results from these three with reported results in benchmark sets from the other assemblers and the results were better for these three.
From these three I get the best results with s-aligner but it has some issues: it actually requires some practice using it and it is slightly slower for this kind of data (<30kbp genomes). In the other hand, it usually has lower requirements in memory than other assemblers except maybe Megahit.
Here there are some comparative results (NGA50) using public benchmark sets:

Obviously, I am kinda biased as I am the author of s-aligner but the data doesn't lie.
Getting the best results also may not be the first requirement in some cases, like if you want it fast and with a fast learning curve, and in this case Megahit I think is usually quite good. But in these cases, you also can design a script that will automate s-aligner for your specific needs. You can contact me if you need help.
BTW, s-aligner also gets better results with other kinds of data like bacterial genomes, transcriptomes, or small eukaryotes but always at the cost of requiring more running time or more computational resources (running it in AWS, for example).

Many thanks @juanjo75es. The key take home message thus far is the consistency of the benchmarks used and s-align does resonate strongly on this index. I'm considering this from a conservative point of view. The variable performance of Spade/MetaSpades is weird (everyone is aware of Spades). Run time is not the issue as you've pointed out. I've some reading to do.

Great. Indeed with the last version t actually gets better results than what it's reported in the document and does not require Flye even for larger viruses. If you want to test it at some point, let me know and I will provide some advice and the last version.


text:
Caching and parallelization
Currently algorithmically developing a calculation for "co-dependent" mutation (quotation appropriate). The algorithm uses chunked biopython alignment objects which are parallelized.
Aim
One optization is caching a database within parallelization, but I ain't sure whether the parallelization duplicates the cache.
Method
Each alignment chunk (biopython object) is passed onto a separate CPU and refers back to the original database to obtain its metadata in the final stage of the calculation (dragging metadata around this algorithm makes no sense).
The idea is to cache the database to prevent pulling it in from the harddrive:
from multiprocessor import Pool
from myOOPlibrary import mystuff

Facade design, the parallelised index to the alignment chunks enter into the myOOPlibrary:
from cachetools import cached, TTLCache

cache = TTLCache(maxsize=1, ttl=500)

Class mystuff (stuff2):
... # init omitted

    def mymethod(self, ..):
    # to be updated

Class stuff2 ():
... # init omitted

    @cached(cache)
    def getDB(self, ...):
.... #  pull it in (harddrive) and parse to a dictionary
    return metadataDB

Question
I'm wondering if caching is the right approach in parallelization, specifically whether this will end up with multiplying the database held in RAM...
Alternatives
The alternativees to caching are:

IO to JSON (not cool)
IO via pickle (binary IO): a JSON will readilypickle

The possible advantage of the alternativees are the database is written once and accessed by all cores for certain.
Underlying question

When caching is the database  cached for each core? If so the RAM overhead becomes an issue because the database is duplicated total CPUs - 1 times
Alternatively is the cache shared between CPUs centrally and each CPU refers to this location?

As currently written, this might be more appropriate for Stack Overflow

Possibly but parallelization of biopython objects must be fairly common in algorithm/pipeline development. Alot work with biopython here

What are you using for the parallelization part?  If it's just the built-in threading module, I don't think anything will be duplicated there, you just have to make sure you don't step on your own toes with any read/write aspects (e.g. use something like Queue objects if needed).  But in that case, why using a caching library at all, versus just a single shared DB object you load once?  (More of the rest of the code might help clarify.)

Updated, @Jesse comment is cool and closes in on an easy solution.

If python code is the bottleneck, rewriting in C/C++ may get you a 50-fold speedup.

Thanks, full Cython coding is the next development (C within Python and compiling into C). You're 100% right @user172818 what I've not described is there's extensive looping and that will be speeded up via C. For numpy there's no gain in Python (which I'm sure you know). Coding in raw C ... lots core-dumps and not quite so enthusiastic.

Summarising

@Jesse's comment I think is the right answer, open the database outside parallelisation. What I didn't thoroughly describe (and what was essentially requested) is the architecture, the alignments entering into an entire OOP library and the whole library is parallelised, i.e. every new alignment entering library is on its own CPU. The answer is obvious: instantiate the OOP library with the entire database and thats a key purpose of the design.
The technical answer is joblib. This is built to cache across parallelised archetechture as the joblib team point out here. The criticism is that there's parallelisation within parallelisation. The discussion* is as much in the comments as in the response.
@user172818 response is correct too, the algorithm would benefit from C. Its a lot of code to rewrite in C, but sidestepping via Cython is a future development (its a hybrid language connecting C and Python).

*, https://stackoverflow.com/questions/25033631/multiple-processes-sharing-a-single-joblib-cache

There was a comment hinting this is Stack Overflow, that's difficult. Ultimately the technical answer (which is Stack Overflow) isn't the correct answer as far as I understand OOP.

Feedback
Just to wrap this, point 1. worked fine, thank you. It means a database sits on each processor in parallelisation, but is imported  once. In truth, it worked nicely and provided cleaner code. Re-writing to instantiate once would affect other stuff so would be more costly than a cost-saving.
Python can be very fast (e.g. numpy), but also really slow (hence it has critiques), extensively parallelising helps bridge the gap.
Python parallelisation has critiques which I strongly, strongly believe are unfounded: I parallelise loads via multiprocessor library no problems.
Caching here, i.e. parallelising across a parallelisation could stretch things too far.
Last point
Point 3, re-writing in Cython/C, is a greater priority than the testing required for parallelised caching: can't parallelise all data and that code is better as C.


text:
HDF5 and BioSQL solutions
I'm looking at better database/storage solutions for NCBI virus data, with all attributes particularly year and country of isolation, together with structural data, possible antibody data, T-cell data and bioinformatics data such as %GC etc... The amount of data is comparatively small, with a maximum size of 3 x ([6000 < no. viruses < 20000] for 10 genes) and each entry will contain nucleotide, amino acid and one or more pieces of alignment data (i.e. with indels).
The data will be used for numerous different calculations from structural biology, phylogenetics, and deep learning. Some of the calculations will require parallel processing. A typical calculation would be to identify a clade/similarity/difference then screen the database to look at the isolates information.
Previously I've used lots of flat files which are loaded at initiation and are held within a data pipeline either as an object (Python) or especially a multi-dimensional hash (Perl, preferred).
The two solutions that I think look suitable are:

BioSQL - The official Genbank solution, however my SQL commands are very basic and the documentation appears minimal (for BioSQL)
HDF5 (hierarchical data format) - this looks good because it uses complex hashes (dictionaries), and allows layering of multiple information onto a single Genbank sequence code (key). It has strong documentation. However, I understand this isn't a "database" (long term storage), but thats how I intend to use it.

JSON/Pickling I had a look at JSON, which to me looks okay, because I'm fond of complex hashes / dictionaries, but it appears to be more of a Web solution. Finally, "pickling" (Python) is cool, very cool, but I understand this not a database solution, because it risk compatibility issues, e.g. with future Python releases.
I am interested in knowing what archival systems board members are using, and a critique whether my assessment is ok for moving forward.

If you are querying repeatedly and have multiple features on which you are applying search criteria, use a proper database with indexed features/columns/fields.

Just to conclude on the solution to this ancestral question. HDF5 is very useful, but there are better solutions for this specific question if it is focused on pandas dataframes. If you've lots of data the pythonic way to do this is pandas and therefore storing as a pandas dataframe makes sense.
For long term storage, which was the original question, the storage of choice for pandas is parquet. Described here

Apache Parquet provides a partitioned binary columnar serialization
for data frames. It is designed to make reading and writing data
frames efficient, and to make sharing data across data analysis
languages easy. Parquet can use a variety of compression techniques to
shrink the file size as much as possible while still maintaining good
read performance.

In the following example, two methods parquet methods are used for a pandas dataframe proteinDF which is first written and then imported.
from fastparquet import write
write('outfile.parq', proteinDF)

and
import pandas as pd
pd.read_parquet('proteindata.parquet', engine='fastparquet')

The above approach appears to have received most current developmental attention. However,
import pyarrow as pa
import pyarrow.parquet as pq
arrow_table = pa.Table.from_pandas(proteinDF)
pq.write_table(arrow_table, 'myoutpath', use_dictionary=False,
                   compression='SNAPPY')

and
import pandas as pd
pd.read_parquet('example_pa.parquet', engine='pyarrow')

The above approach is compatible with numba. 'Multithreading' is further described here.
The differences between the methods appear slight. Parallelisation is a big deal and depends on the problem.
A general discussion on the methods is here
HDF5 vs parquet
The backdrop for parquet is the method is used for batch I/O not record by record writing to disk and that is this question and for that is pretty fast (at a guess bit better than pickle, but a lot more compression).  HDF5 is perfect for rapid record by record writing and that is what it is designed for, so 'benchmarks' comparing I/O speed are not realistic unless they compare batch vs. record by record timings, where HDF5 excels. parquet does achieve significantly better data compression and is considered a long term storage solution. In my hands HDF5 for fasta writing results in around half the size of the flat file (not much).

Just to conclude, @user172818 and @AlexReynolds are collectively saying either learn SQL or stick with CSV.

Update,
I can provide an answer. It's HDF5 with pandas dataframes (Python).
Simply
df.to_hdf('path/file.h5', 'key', format='table')

This isn't SQL but it works:
What it gives,

it's fast, that's it's purpose (not my main reason though)
it perfectly preserves all pandas dataframe information, this is  important
The 'key' ( object) allows you to store multiple dataframes into the same file, you just provide different key names

The last point is very useful, because you can pull out a dataframe of choice and the SQL commands are mimicked within pandas
I ran into trouble using CSV to store/recover complex pandas dataframe formats.
So it's good enough and much faster than SQL. I could have problems with format changes in time (unlike SQL), but for the here and now it works well.

Just to add a caveat, for dictionaries the json dump works well, but it is just shunting dictionaries around, so isn't relevant.

Hashes are memory intensive. Perhaps keep an eye on your memory usage, if you go this route. This is why indexing and sorting approaches are used.

Thanks, I appreciate critical feedback. In defence of Perl v5.20 now has postfix dereferencing ->@* syntax, which a step forward for hashes. The data loads involved are less than a single higher eukaryotic genome. Obviously not doing dataframe style calculations in a hash ... its doable done (higher order programming), but possibly surpassed.


text:
eggNOG: DNA-DNA or protein-protein orthologue identification?
I'm trying to work out is how reliable/accurate my eggNOG metagenomic outputs are. Critical for this is how its identifying orthologues.
Background eggNOG is commonly used in metagenomics gene identification. Take a sample, randomly sequence all the DNA therein, then "annotate" via eggNOG. It identifies the putative gene function against lots of databases e.g. NCBI, KEGG, Brite, Go ...
eggNOG uses is own database for gene identification via 'seed orthology' each seed has its own KEGG ID.
Question Does the comparison between the query (metagenomics data) and the seeds take place using DNA-DNA comparisons? Alternatively, is it performing protein-protein comparisons? Algorithmically, by the easiest is DNA-DNA, however DNA-DNA would result in limitations of interpretation ('cause of complex evolutionary stuff). It could be using an optimised a protein database instead and returning the output as DNA.

I can answer this now. I've been 'under the hood' of eggNOG and its using "hmm" (hidden Markov model) very likely, almost certain, the Eddie (sp?) group's hmmer to cluster proteins associated. It represents the cluster members strictly via seed orthologue code, which then are aligned to form a 'raw' and 'trimmed' alignment. Each alignment is assigned a code which corresponds to its know function. This is a database of alignments. Once it gets a protein - protein hit against one of its orthologue (I suspect it will be using a nested system) it will then assign the query to that alignment - each of which is assigned to a protein function.
It representative database not an all encompassing database, i.e. sequences are being omitted.
Taxonomically it is simply declaring a the familial name. There will be an internal database mapping the seed orthologue onto the species, but that is not being publicly disclosed.
I make a better explanation than that if anyone wants it.

There is some description in the original paper


text:
Remote NCBI's Blast Perl API: maximising hits returned via command line
Question and Background I'm using NCBI's Blast Perl API to interrogate Genbank for example blastp. I don't get the throughput using their point and click web stuff.
I tried Biopython's Bio.Blast import NCBIWWW but its under access restrictions from the host, notably attempting to use it during the day (works sporadically). The API however gets through and runs as fast as a web search ... that's good.
Here's an example of the command I'm running:
CMD=Put&PROGRAM=blastp&DATABASE=nr&HITLIST_SIZE=2500&ALIGNMENTS=2500&QUERY=%3Emytest%0AMSILANKDTR%0A

The above query gives default Blast output. What I'm trying to do is increase the number of hit recovered to the maximum (5000). The options for this are HITLIST_SIZE, ALIGNMENTS, DESCRIPTIONS. In the above example I'm attempting a hit_size of 2500 ... but the same output is returned as for default. More generally  I can't get any of the "number of hits/descriptions/alignments returned" options to work, the default is returned each time.
Has anybody solved this?
Typical options for the API are given below
Typical options
Parameter   Definition  Type    CMD Allowed values
QUERY       Search query    String  Put *   Accession, GI, or FASTA.
DATABASE    BLAST database  String  Put *   Database from appendix 2 or one uploaded to blastdb_custom (see appendix 4)
PROGRAM     BLAST program   String  Put *   One of blastn, blastp, blastx, tblastn, tblastx. To enable megablast, use PROGRAM=blastn&MEGABLAST=on.
FILTER      Low complexity filtering    String  Put F to disable. T or L to enable. Prepend “m” for mask at lookup (e.g., mL)
FORMAT_TYPE Report type String  Put, Get    HTML, Text, XML, XML2, JSON2, or Tabular. HTML is the default.
EXPECT      Expect value    Double  Put Number greater than zero.
HITLIST_SIZE    Number of databases sequences to keep   Integer Put,Get Integer greater than zero.
DESCRIPTIONS    Number of descriptions to print (applies to HTML and Text)  Integer Put,Get Integer greater than zero.
ALIGNMENTS      Number of alignments to print (applies to HTML and Text)    Integer Put,Get Integer greater than zero.

Just to feedback on the answer and comment

Using -remote by locally installing blast as recommended by @Scott. This is a valid answer so definitely upvoted and recommended its further upvoted. I had tried this but not mentioned it in the question (really sorry), but it times out due to traffic issues. The API overcomes traffic issues. Again, @Scott's is a valid solution (considering I'd not explained I'd tried it).

More details ...

The Perl API is here https://blast.ncbi.nlm.nih.gov/docs/web_blast.pl
Line 95 below ...
$args = "CMD=Put&PROGRAM=$program&DATABASE=$database&QUERY=" .
$encoded_query;

This is changed to,
$args = "CMD=Put&PROGRAM=$program&DATABASE=$database&HITLIST_SIZE=$maxhits&QUERY=" . $encoded_query;

Also need to add an additional line at line 73
$maxhits = shift;

From the command line,
perl ./web-blast.pl blastp nr 2500 ./example.fa

(Note, The variable order is switched in the program).
example.fa - e.g.
>rleguminosarum
WGFDGSSTQQAEGRSSDCVLKPVAIYPDPARTNGALVMCEVMMPDGVTPHASNARATILDDEDAWFGFEQ
EYFFYQNGRPLGFPEQGYPAPQGPYYTGVGYSNVGDVAREIVEEHLDLCLAAGINHEGINAEVAKGQWEF
QIFGKGSKKAADQIWMARYLLQRLTEKYGIDI

I'm proficient in Perl5, i.e. there's no error in changes to the code in the API, albeit my expertise is object oriented and high performance Python.

Limit on the request ... I'm testing that now.
First results eeek @terdon might be right ... :-(

Can you give us a minimal example of the actual code that we can test on our local machines, please? I wonder if it could be a limit in individual post requests and if you might need to set the downloading of the results in some kind of loop to get N at a time.

While I haven't used your exact method above, using the BLAST command line tools return the desired number of target sequences
blastp  -query query.fa \
        -task blastp-short \
        -db nr \
        -max_target_seqs 5000 \
        -outfmt 6 \
        -out query.out \
        -remote

query.fa
>mytest
MSILANKDTR

The resulting query.out file does contain 5000 unique hits.

I'll mark this as "accepted" because this is a good solution given the original question I asked (I tried this approach, but didn't state it in the question). The full explanation is given below and the above answer should be read within the context of my explanation for future reference. Note unless the other Mods over-rule my "accept" please keep this comment here to provide future readers with the full solution and thanks to @scot

Right here's the response from NCBI,

Please send your future blast related questions to blast-help@ncbi.nlm.nih.gov directly. The gb-admin alias is for submitters with questions related to their sequence deposition/submission.

For your issue, you should look at the following two parameters:
DESCRIPTIONS  Number of descriptions to print (applies to HTML and Text)  Integer Put,Get Integer greater than zero.
ALIGNMENTS    Number of alignments to print (applies to HTML and Text)    Integer Put,Get Integer greater than zero.

Do please note that asking for this level of hits saved/returned could be problematic for the search. They could make the search exceeds the cpu time allocation and lead to error.

I tried this ... it was the same as "HITLIST_SIZE".
I figured the answer ...
Why its happening? The API uses a default parameter for the *maximum number of hits. It's imposing a limit on the number of hits returned, but this wasn't stated in the documentation of the parameters.
Whats the limit? It's >70 but <= 100. I suspect the limit is 100 hits. Thus anything above this its set at 100.
Summary The reason why the API is working and other approaches simply time out is NCBI is imposing an upper limit on the number of hits for the API and the Blast servers are not permitting requests above this threshold at peak time. Thus because other approaches are not imposing the restrictions on the NCBI server, the API gets through, whilst the other methods (blastp directly and Biopython) are timing out.
NCBI are recommending the API is because it is automating a 100 hit cut-off and the servers are deprioritising limits beyond this
I have filed a note of advice to NCBI because it should be clear in the documentation there's an upper limit on the HITLIST_SIZE parameter. It took me at least a full half day to figure this out (which isn't cool for limited documentation).

Proof of concept (untested) Thus, using @Scott's solution
blastp  -query query.fa \
        -task blastp-short \
        -db nr \
        -max_target_seqs 100 \
        -outfmt 6 \
        -out query.out \
        -remote

This will work with the -max_target_seqs 100 and thats all the API is doing.


text:
Parameter optimisation of minimap2
I'm leveraging minimap2 to map select genes from short-read microbial fastq metagenomics zip files.
minimap2 --heap-sort=yes --splice-flank=no metagenome.fastq.gz fastaDB.fas

What I noticed is the low tolerance of mismatch when I need to maximise mismatch because microbes have loads of genetic diversity. Thus the collection of references sequences are never really a reference, in fact the "microbe reference genome" is a misnomer, because in situ diversity of microbes is huge.
Question How do you relax the mapping stringency criteria for minimap2? If you are familiar with Blast parameters this would help alot here. I know blast but I've never parameterised it.
The options are this are:
Alignment options

-A Matching score [2]

-B INT    Mismatching penalty [4]

-O INT1[,INT2] Gap open penalty [4,24].

-s INT    Minimal peak DP alignment score to output [40].

--score-N INT
Score of a mismatch involving ambiguous bases 1

Anything else that is relevant.
Note, -f is a famous option, not sure it helps here.
I'm going to code up a hyperparameterisation unit-test, but understanding the parameters first would really, really help.
Alternative approaches welcome as well.

Plugging minimap2. I code high performance (hadoop-style, back-tracking, graph, traverse tall buildings in a stride  etc ...). The coder that wrote minimap has knowledge beyond, i.e. very fast run times directly from gzip file. I assume they are indexing then chunking the raw zip file then unzipping each 'chunk' JIT (just in time).

Thanks @MaximilianPress great answer. I'll fill in the backstory.
I'm trapping total genetic diversity across Genbank and the metagenomes. My protocol to do both:

Blasting out Genbank and then
piping the output across metagenomes via in this case minimap2.

Hence speed is key. So I ain't replacing Blast with minimap2 but leveraging minimap2 via the output of blast.
However on reflection of @MaximilianPress's discussion the answer is to generate a gold standard of diversity to assess alternative approaches.

First of all, use short read option
The short read preset (-x sr) for minimap2 is heavily recommended for doing such mapping:
minimap2 -x sr
Heng Li on minimap2 for related genomes
Heng Li (minimap2 author) has some suggestions on a somewhat related github issue (though it is now 4 years old):

Minimap2 is not really intended for this purpose. For short reads he recommends stampy, but it is clearly quite old (python2.7) and does not seem to have much recent developer history.
Failing that, he suggests an option set such as follows:

If you want to give minimap2 a try anyway, you can use the following setting:  minimap2 -axsr -A2 -B4 -O4,24 -E2,1
If minimap2 is not sensitive enough, reduce -k and -w at the cost of mapping speed. Reducing -w should have relatively lower impact on speed, but reducing -k should be more sensitive.

I would suggest engaging with Heng, he is quite responsive in general and may have more suggestions.
Other aligners for metagenomics
I would suggest that you look for more possibilities for alignment, or possibly not using alignment in the first place. However, this is very sensitive to your ultimate application, which I'm not sure that I understand from your post.
I was pleased to see that others have reported good results using minimap2 for metagenomics, which I was skeptical about, having been out of the field for a bit. Usually in metagenomics there is an expected tradeoff between sensitivity and specificity for the reasons that you state (this paper is old but it uses BLAST, so it is a sort of best case for sensitivity on unrelated genomes at the level of divergence that they test).
There are of course new shotgun metagenomic aligners coming out all the time, here is one recent example that I know nothing about but which may be interesting as a comparison.
If you are really looking for a rich parameter set that allows low-stringency mappings, I would suggest going back to BLAST, though it will be much slower. minimap2 and related algorithms are to some extent dependent on the ability to find short seed matches, so there will necessarily be less than perfect sensitivity in your application.

Many thanks. Cool answer. I've addressed the issues you've below the final "-------"

Just to follow-up on the performance of minimap2: the solution @MaximillianPress described works. Thousands of hits are obtained for comparatively obscure genes (some will be duplicates). Parameterisation gives a 2-log10 increase (astonishing).
The issue is that its:

relatively slow - precisely as described by @MP;
very, very sensitive to parameterisation, just a single integer changes the output by a couple of thousand hits.

I can see why there was concern about using it - such notable sensitive is an issue - and alternatively its clear why other investigators have had great success with it - hit the "sweet-spot" and its a bonanza. For my purposes it's fine - it's not a quantitation of complete genetic diversity, as long as I capture a reasonable proportion of diversity that is fine. Clearly, they'll be times the "sweet-spot" is missed but thats still okay (famous last words :-) ).