text:
How to display novel genome assemblies or uncommon genome assemblies using the UCSC Genome Browser?
I want to display E.coli BW25113 (GenBank: CP009273.1) strain in UCSC browser. This strain is not listed in http://microbes.ucsc.edu/ browser. How can I display E.coli BW25113 assembly in the browser?

I found two methods:

You want to visit this page for instructions for novel assemblies: http://genomewiki.ucsc.edu/index.php/Assembly_Hubs .
However, before doing that I highly recommend that you check Genome Assembly Hub first:
https://genome-test.gi.ucsc.edu/gbdb/hubs/genbank/ . There is a good chance that your genome or strain is listed there. Here is an example of how you can do that. Example specifically shows how you can display E.coli BW25113 strain in the browser.

Instructions:

Go to NCBI Nucleotide: https://www.ncbi.nlm.nih.gov/nucleotide .
Enter your strain name: E.coli BW25113 (GenBank: CP009273.1)
Find Biosample ID on the page, located under DBLINK: BioSample: SAMN03013572
Go to GenBank assembly hubs: https://genome-test.gi.ucsc.edu/gbdb/hubs/genbank/ .
Select bacteria. This will take you to "Bacterial assembly hub"
at https://genome-test.gi.ucsc.edu/gbdb/hubs/genbank/bacteria/bacteria.html .
Look for E. coli hubs. At the time I looked #15-18 contain E. coli strains.
Click on each hub and search for SAMN03013572 using browser search function (Command+F on macOS). I found my genome in hub #16.
Click on "common name" in the column-3. This will take you to
https://genome-test.gi.ucsc.edu/cgi-bin/hgGateway?hubUrl=http://genome-test.gi.ucsc.edu/gbdb/hubs/genbank/bacteria/16/hub.txt&genome=GCA_000750555.1_ASM75055v1&position=lastDbPos , where FindPosition window will be populated with your genome Assembly.
Click on Go. This should take you to the UCSC Genome Browser with your genome assembly.

Update: The step-9 is no longer working since UCSC genome browser resources described above were experimental prototype proof of concepts
that have since been discontinued according to the Genome Browser Team. Once I find an alternative method, I will post it here.

I am not sure if this qualifies as an official answer. I have contacted UCSC bioinformatics team. They were were very kind and told me that they would add any track that was requested by users. After my request, they kindly added the track for E.coli BW25113. If you want your reference tracks to be added you could submit a request at this website:
Geneme assembly search and request page.
Then, at the upper right corner click on "Can't find your assembly?" button and submit your request. See the screenshot below from the webpage.


text:
How to rescale a function?
When we have a function which takes values in [a,b] and which returns values in [c,d], what transformation can we do so that the function takes values in [0,1] and returns values in [0,1] without the curve changing in both cases?

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer.

Please provide a concrete example that doesn't involve something like "function(x) return x;", which is what you're implying.

Assuming you want a linear transformation (which is what "without the curve changing" sounds like to me), to rescale from range c..d to range 0..1 subtract the minimum then divide by the difference, i.e.
$val_t = (val - c) / (d-c)$
The inverse of this can be applied to work the other way, i.e.
$inp_t = val \times (b-a) + a$
For a function $f(x)$, with ranges as you have stated, the full transformation of input and output would be:
$g(x_t) = \frac{f(x_t \times (b-a) + a) - c}{d-c}$

Very cool response

source: Stack Exchange

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

sort: votes


text:
Are the information stored in the PBD files in the protein data bank complete or incomplete?
I know that there are thousands of PDB files stored in Protein Data Banks.
Are all these files complete in terms of the information they store?
If YES, do the individual files get updated from time to time?
If YES, why/when do they get updated?

A PDB file is the last step in a series of operations.
For a crystal structure in very simple terms, it consists of the coordinates of the atoms placed to best fit the density map, which in turn (with these atoms due to the phase problem) is a conversion from reciprocal space which is the result of transforming the quantified dots from the pretty diffraction pattern picked up from the detector during the experiment.

Are all these files complete in terms of the information they store?

Not at all, a PDB file is not complete. It nearly complete when you count density in there.

If YES, do the individual files get updated from time to time?
If YES, why/when do they get updated?

Only, minor format upgrade, eg. https://www.rcsb.org/versions/1UBQ.
It is however very common to polish structures.
There is a tool called PDB-Redo that is commonly called in many pipelines for example to "upgrade" older files.
Some folk open the files in Coot and let it do its thing.
For most molecular mechanics experiments it is common to energy minimise the structure (freely, restrained by original position or restrained to density map). It is uncommon in several docking tutorials simply to keep it simple.
However, for most uses, if an aromatic ring is not quite planar or a bond is longer than most by a fraction of an Ångström, nothing happens.
Also worth remembering many structures are incorrect and not just imperfect: there may be actually an isopeptide bond in the core of the protein (many did go unnoticed) or an isoaspartate conversion happened or the loop is actually split due to an Asp-Pro pair autocatalysis, waters get mistaken for ligands or sidechains. If these are polished, they may get flagged as problematic, but they may be made worse.

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


text:
What is DSSP algorithm for?

DSSP - Wikipedia

The above article says that DSSP is a hydrogen bond estimation algorithm. Then the article also says that this algorithm is "assigns secondary structure to proteins".
If both are true, how does hydrogen bond estimation helps us to assign secondary structures of proteins?

"Secondary structure is formally defined by the pattern of hydrogen bonds between the amino hydrogen and carboxyl oxygen atoms in the peptide backbone" (source)
To be more specific: secondary structures are formed explicitly because of hydrogen bonds. So when DSSP says that it estimates hydrogen bonds as a proxy for secondary structure classification, this is likely why.


text:
OMIM diseases to UMLS mapping
Is there any public available mapping between Online Mendelian Inheritance in Man (OMIM) diseases and Unified Medical Language System (UMLS) concepts?
OMIM has mapping for genes to entrez or ensembl, but has no mapping for diseases? Does anything exists that would allow me to map OMIM diseases to UMLS?


text:
Python argparse error
Hi I am trying to use this python script (running with bash script) to do my analysis. I am getting the following error:
File "my_python_script.py", line 33
logging.info("Reading input2 file")
^
SyntaxError: invalid syntax
My Python script is as folllows and albeit I omitted the indentation and my full apologies therein, please keep in mind I'm not the coder in this instance.
import itertools
import argparse
import logging
import sys
import numpy as np
import scipy
import pandas as pd

def parse_commandline():
"""Parse command line parameters.
Return type: Object
"""
parser = argparse.ArgumentParser(description="Association of variants")
parser.add_argument("--input1", metavar="FILE1.tsv", help="gene data file", required=True)
parser.add_argument("--input2", metavar="FILE2.tsv", help="mutation data file", required=True)
parser.add_argument("--output", metavar="FILE.tsv", help="Output file name", required=True)
if len(sys.argv) == 1:
    parser.print_help()
    sys.exit(0)
return parser.parse_args()

def main():
params = parse_commandline()

#load genedata
logging.info("Reading input1 file")
try:
gene_data = pd.read_csv(params.input1, sep="\t", index_col ="Model")
sys.exit(0)

#load mutation data
logging.info("Reading input2 file")
try:
mutation_data = pd.read_csv(params.input2, sep="\t",
index_col="Mutation")
sys.exit(0)

##All mutations and gene in data files
all_mutations = list(mutation_data.index)
all_genes = list(gene_data.columns)
cell line Models in mutation_data file
models_mutations = list(mutation_data.columns)

## models in gene_data file
models_genes = list(gene_data.index)
intersection_set = set.intersection(set(models_genes),
set(models_mutations))
models_common = list(intersection_set)
unique_series_gene = gene_data.index.isin(models_common)
end_gene_table = gene_data[unique_series_gene]
transposed_mutation_data = mutation_data.transpose()
unique_series_mutation =
transposed_mutation_data.index.isin(models_common)
end_mutation_table = transposed_mutation_data[unique_series_gene]

##calcuate p-values
all_p_values = []
genes = []
mutations = []
for gene, mutation in itertools.product(all_genes, all_mutations):
gene_table_to_test = np.array(end_gene_table[gene])
mutation_table_to_test = np.array(end_mutation_table[mutation])
p_value = scipy.stats.wilcoxon(gene_table_to_test,
mutation_table_to_test)[1]
genes.append(gene)
mutations.append(mutation)
all_p_values.append(p_value)

## calculating adjusted_pvalues
adjusted_pvalues = p_adjust_bh(all_p_values)
d = {'Gene': genes, 'Mutation': mutations, 'pvalue':
all_p_values,'p_adj':adjusted_pvalues}
df = pd.DataFrame(data=d)

# Print output file
logging.info("Writing into file")
df.to_csv(params.output, sep="\t")

if __name__ == "__main__":
main()

My Bash script
#!/bin/sh

### This script will run python script

python3 my_python_script.py --input1 Gene.tsv --input2
Mutation.tsv --output tmp
diff tmp adjusted_pvalues.tsv
rm tmp

Please correctly format and indent your python code

@Chris_Rands Thank you. I have corrected the format now.

It would be helpful if you could add small, example input files.

Could it be a line break in the line? mutation_data = pd.read_csv(params.input2, sep="\t",
index_col="Mutation")

It is certainly possible the 'error' is the input file 2.
However, both try: routines are weird and if this is the real code it would fall over. Something isn't quite right however, which I'll touch on at the end.
What the code is supposed to say is:
try:
    load stuff
except:
    raise Exception ('Input file 2 is has deviation from csv format and will not load')

raise Exception isn't needed, the following would be a more harmonised approach:
try:
    load stuff
except:
    print ("Input file 2 is has deviation from csv format and will not load')
finally:
    sys.exit('Abandon ship!')

The idea is that if the CSV file didn't load, it would stop the code and tell you approximately where to start looking. You might then ask, can it actually sees the file and then peer into the loaded CSV. If it loads the code continues.
You could just stay:
try:
    load stuff
except:
    sys.exit('Oh the file did not to load')

This isn't cool, but who am I (who sticks Exceptions in except statements), but it will work and it could be that is the actual code.
What the code actually says is
try:
    load stuff
    sys.exit()

or
try:
    load stuff
sys.exit()

Same outcome the code stop because sys.exit stops everything ... nothing is left. What the code should really say is not just perform a try: except: but perform a 'getter' and 'setter' on the input file to check the CSV format, this is because your the user of the code not the developer and the coder has to put safe guards in because stuff happens.
sys.exit is really used in development, for example in debugging.
However, somehow you got past the first try:, sys.exit and I don't know how you did it, because that should not have been possible.  I would carefully check the original code.
Just to point out as @ChrisR mentioned, if this was your actual python code it would fall over because you've omitted all the identation, which in python is essential. Its not like shell where indentitation matters not a jot.
Finally, as @supertech has mentioned, there are numerous line breaks in the code and whether you are in shell or python this can cause problems. I think the line break(s) are within an extra layer of identation, which Python can read through (I never understood how).
Personally, I'd go back to the coder and point out the try: sys.exit() and just point out any beginners webpage on try: except:, this is not advanced coding. If anything sys.exit() risks weird behaviour for a user particularly as no distinguishing information is provided.


text:
How to run MD5 check for multiple fastq files in different subdirectories?
I have received Illumina sequencing reads for 100 samples. I have 8 R1.fastq.gz and 8 R2.fastq.gz files for each sample in each subfolder. I want to run an MD5 check for all the fastq files in each subfolder.
My folder structure looks like:
/mydata/sequencing/clean
/mydata/sequencing/clean/Sample1
/mydata/sequencing/clean/Sample2
......
/mydata/sequencing/clean/Sample100

I am using the following code being on the /mydata/sequencing/clean directory.
find . -type f -exec md5sum-lite {} \;

This code is printing the md5 result on Terminal.
How do I get these results printed to a .txt file? Besides, is there a more elegant way to do this MD5 check?
Thank you.

You posted the same thing in Biostars.org :) https://www.biostars.org/p/9509386/

Haha yes. Because sometimes I do not get a reply to my questions (may be they are  too naive some times) :(

Glad you got the answer.

thats because its on the wrong stack

Anyways it is always great to receive feedback on bioinformatics problems. Nice to see that many roads lead to Rome.

Small modification,
find . -type f -exec md5sum-lite {} \; >> output.txt

If you use > the output will almost certainly contain a single record.

Thank you for the prompt reply. Is there any other code that you know of for looping the md5 command over different subfolders?

Yes of course ... there are loads of ways to do this. A bash script is good (separate question IMO). However, the > and >> will still be central. In terms of elegance, your one liner looks cool. BTW >> is very likely the solution because of the find approach. One liners are cool BTW.

Thank you very much.

Your welcome. Your  second question (alternative code) is worth asking because of the variety of ways it can answered, e.g. a bash for loop, Python subprocessing, Perl one liner are the  approaches I can think of.

Yeah. totally agreed.


text:
TGCA LAML rsem data
How do I test the gene expression differences between multiple groups?
I have normalized TCGA RSEM data, this is the dataset
Human__TCGA_LAML__UNC__RNAseq__HiSeq_RNA__01_28_2016__BI__Gene__Firehose_RSEM_log2.cct.
I have different groups of data which were defined in terms of blast percentage as well as karyotype. Here I dont have control sample as such but I do have normal karyotype. Now since these are all already normalized data what would be ideal test to find the differences.
What methods can I use if I have to compare between these different sub types of data? So far I have always used raw counts and used DESeq2 to perform differential and other downstream analysis how in this case should I go ahead?
UPDATE
I'm trying to find the important or significant genes that defines a particular subtypes of cancer in this AML.
So far I have tried clustering and PCA on it I do see some clusters and patterns. Now I want to take those cluster for downstream analysis such as pathways plus we want to validate some of the key genes that might come out.
So far I have done as such I get raw counts which i use with deseq2 to run differential analysis between multiple groups or subypes.
How to do the same here when I don't have raw counts to start?


text:
pathogeniticy prediction tool
I have a question related to the evaluation of variant pathogenicity in the context of human genetics. There are many predictors available, however, is there a specific predictor that takes into account the 3D protein structure as well?
Imagine there is a mutational hotspot, i.e. a sequence of amino acids of a protein region that are pathogenic (assessed by functional studies, for example). And then, let's say several hundred residues far from this region there is a variant whose pathogenicity status is unclear. Now, what can happen here is that this variant, although far from the mutational hotspot (considering the protein secondary structure) can actually be quite close to the mutational hotspot in the 3D space. I mean, this variant may actually form a part of a functional site, together with residues from the mutational hotspot. This could indicate that this single variant is pathogenic as well, although there is no direct evidence.
And so, is there is a predictor that takes 3D structure into account?
In other words, I'd like to ask a question: Are there any pathogenic variants (defined by clinvar) within a sphere with a radius of 100A centered on a particular residue of interest?
Or, are there, for example, any glycine residues within such a sphere?
It's not really a prediction, it would be just a search algorithm, taking 3 parameters:

protein 3D shape
search sphere [center residue; radius]
needle - the element I'm searching for, such as AA residue type, for example

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer.

Venus
I am the developer of venus.cmd.ox.ac.uk, which is a site that given a mutation makes the variant (and gives the ∆∆G) and shows what residues are nearby and whether these appear in clinVar or gnomAD (and scores them too) or have some annotation such as disulfide or post-translational modifications. The aim is to allow the user to see and judge for themselves what is going on as there are multiple possible mechanisms which may make a variant different (and then saving and sharing the results via Michelanglo).  In other words it is not a predictor but a data aggregator of soughts.
It shows, in addition to the structure in the NGL viewport, a table where the nearby residues are listed (and if clicked focuses on in the protein)
Limitations
An issue according to non-users is that it does not give a statistical relevance score for the features, but I really did not want that as these can be completely off and easy to manipulate and just encourage users to blindly believe the calculations.
An issue according to me, is that it is based on empirical data: if nothing is known to bind or where it binds, it cannot say much, unlike MutPred2, a linear sequence deep-learning tool that assigns a prediction of whether there's a binding site etc. without prior empirical evident.
Very very few users provide a PDB file and this aspect is limited. Venus is that it is centred around canonical Uniprot numbering, so providing PDB files with differing numbering will not work. Furthermore, the app is coded to "give up easily" on a structure and pick the next option, so does not try to parameterise unusual ligands in a provided structure (I notice now a job failure due to the ACE C-terminal cap, for example).
Other options
There are other tools that show the variant in a structural context, such as Missense3D (worth a gander), but do not show other variants or annotations.
However, for a specific query, such as the closest proline to the residue of interest in a structure with no relevant Uniprot, it may be best to simply use a viewer, such as Mol* or a local one, such PyMOL.


text:
Given a 'random' rsID, how can I find the closest (LD) SNP in a given set of SNPs?
I'm assuming that this is such a common question that I don't have to do my own imputation...
I have a SNP from a paper that is associated with a specific trait... Actually I have 8 of them:
https://www.ebi.ac.uk/gwas/publications/33307546
However, when I search my SNP data from 23andMe (v3 NOT v5), only one of them is in the set of SNPs used by that chip, namely rs2109069.
I'm wondering how I find the closest SNP in my v3 data for the remaining 7 SNPs missing in that data.
Is there a quick way to find a 'tag' SNP in v3 for a given 'random' SNP?
I'm aware that a tag SNP will only be as good (on average) as the LD between it and the given SNP.

What is a "closest SNP"? You mean in distance? Why would that be relevant? Don't you want to look for pathogenic or at least LP or VUS variants? Why would the nearby SNP be of interest?

Sorry, thought it was clear. SNP x is associated with a trait but isn't on the SNP-Chip. SNP y is in LD with X and is on the chip.

So you want to try and guess at the variants your sample carries based on their linkage disequilibrium values? And is this using the raw, unfiltered data of 23&me? (if the answer is yes, please don't do this, their raw data are not supposed to be used that way and are not reliable at all).

I don't understand... you mean I cant run 23andMe v3 data through imputation?

No, I mean you shouldn't trust 23&me data beyond the filtered results they show you. If you have downloaded the full, unfiltered dataset, that will contain loads of false positive calls that should not be considered. Imputation won't be a problem, I just find it a bit odd that you would use what is essentially a statistical approach (using LD values as a proxy to guess what variants are really present) when looking at a single, specific sample. Is getting proper sequencing data so you can actually see what you have in the loci of interest not an option?

Yes it's not an option. Please can you provide a reference for the fact that 23andMe 'unfiltered' data contains a lot of false positives?
I had understood that SNP-Chip accuracy was very high when validated by WGS techniques.

I don't have a reference, but I have looked at the unfiltered results of four different (healthy) samples and all of them had variants that, if true, would have made these healthy individuals very, very sick or, more likely, dead. Variant calls need to be filtered and the 23&me unfiltered results are just a list of variants with no information on call quality, coverage, fastq quality or anything that can be used to sort false positives from true calls. That's why they only show the filtered results.

You're talking about sequencing based variation calls. I'm talking about SNP-Chip based calls. As far as I know, 23andMe isn't doing any sequencing commercially.

All I can tell you is that in all 4 cases I have checked, the unfiltered data, the dataset you can download but is not shown in their web interface, was riddled with obviously unreliable calls.

That makes no sense. For example, I ran my own 23andMe data through GenomePrep (results here: https://supfam.mrc-lmb.cam.ac.uk/GenomePrep/results.html?id=48318966999) and you can see that 80% of the calls match the reference. Sure 20% don't, but I wouldn't be me if they didn't.

As you wish. I can only share my experience based on what I have seen and my ~7 years of working in NGS variant calling in a clinical setting. I would never trust unfiltered 23&me results. Their filtered stuff is fine, but unfiltered variant calls are generally not reliable.

That's true when it comes to GBS data, but as I said, 23andMe DON'T do GBS, they use SNP-Chips, which are generally regarded as very reliable.

And I'm sure their filtered data are indeed reliable. Not so their unfiltered data, at least not in my experience. There's a reason they don't show the unfiltered results, after all. Here's a paper that seems to confirm my experience (I only glanced at the abstract): https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7879796/: "SNP chips are extremely unreliable for genotyping very rare pathogenic variants and should not be used to guide health decisions without validation."

This is weird... You don't seem to be understanding what I'm saying, but you clearly are very good in the field, so I'm confused. Perhaps we should have a call, because we seem to be talking over each other. e.g. of course a SNP chip won't see a rare variant... it's totally beside the point... Somehow I think we're talking about different things.

Heh, I'm getting the same impression :)  Feel free to ping me in the site's chat room. Might be easier there.

Also, while I do have experience with NGS in general and NGS in a human clinical setting in particular, I do not have any to speak of with SNP chips. I just know that the unfiltered 23&me data I have seen have consistently been crap. I may very well be wrong or may have misunderstood what you're doing entirely.


text:
What could explain the difference in RNA-seq and Microarray expression levels?
What are some possible reasons why some trends observed in Microarray expression levels is not observed in RNA-seq. Example the difference between 2 cell types for a gene of interest show major differences when microarray is used with one cell type clearly having higher expression but rna-seq for the same cell types shows equal expression levels for the gene of interest?
What are possible explanations for the discrepancy between such results ? Shouldn't they have relatively identical trends otherwise gene expression in one of them is unreliable?

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer.

I think you need to supply examples to support your observations, you are making a very sweeping statement there. For a start are you considering a relative analysis?

That is a good question. I have never done any Microarray experiments or data analysis myself. But I am very familiar with RNA-seq library workflows and data analysis. Library prep methods for RNA-seq involves many steps including use of several enzymes, which introduce certain level of bias. For example, reverse transcriptase used for making cDNA may copy certain mRNA sequences better than others. Other steps that involve enzymes: second strand synthesis, end-repair, ligation, and PCR amplification. Non-enzymatic steps could also introduce biases. I imagine Microarray and probe designs have also significant biases. There could be additional reasons I have not mentioned above, such as improper data analysis. If global gene expression levels are significantly different, I would suspect that it could be experimental reasons e.g. bad lab practices, unreliable reagents etc. I hope this helps.

And the biological state of the cell :-)

Pallie:  I assumed that question is asked about the differences between using two different methods using the same RNA sample. If that is not the case, I agree with you it's not even fair comparison. And question should be edited to prevent confusion.

It depends on the microarray and the cDNA sequencing sample prep.
One difference that I know about is that microarrays can be designed to target specific isoforms (e.g. spanning a isoform-specific splice boundary), whereas the shotgun sequencing approach of short-read cDNA sequencing makes this difficult to achieve (although in many cases isoform counting can be done to some degree by observing differences in coverage across transcripts).
However, microarrays are very heavily dependent on a good annotation. They can only detect existing known sequences, so it's possible that multiple isoforms are covered by a sequence that was previously believed to be isoform-unique, or that a sequence picked up by shotgun cDNA sequencing is not detected at all by a microarray.


text:
Where do I find the primary evidence for regulatory relationships included in RegulonDB?
When I look up high-throughout (ChIP) evidence for binding of the ferric uptake regulator Fur to sites in the genome of E. coli, RegulonDB lists a whopping 4214 targets. (To repeat this search, go to the downloads page, click on "High-throughput datasets",  enter Fur in the right-most search box, and click on the fourth link, for HTRI00000827.) Some of these have the same gene (alternate promoters?), but I think over 2,000 unique target genes are given. When I click through to the supporting citation, Seo et al. say "Collectively, a total of 98 Fur-binding sites were newly identified in this study[.]" Clearly I'm missing something here. Why does RegulonDB list so many more target genes than the original paper? What's the evidence, if not ChIP? Thanks!
Edit: I've tried to reach out to RegulonDB using the contact form on their website, but I am unable to pass their captcha test! I am uncertain what the problem could be. My protocols are completely up to date. I will update the question if I manage to contact them.
Edit: the paper describing incorporation of ChIP data into RegulonDB is here. It is obviously well-informed, including senior authors from labs behind multiple ChIP datasets. Tables 4 and 5 describe the quantity of new interactions added and are closer to the original reference than to what I am seeing from the website.

I've looked into this a little and I can't find an answer. If I had to guess, the 4214 targets represent sequences with any number of ChIP hits, and the final 98 Fur-binding sites are those with ChIP peaks above the described threshold with paired RNA-seq showing the effect of Fur on expression at those loci. It's a bit frustrating that the "Effect" column of the dataset is "not determined" for all positions.

I reached out to the maintainers, and these results are derived from re-analysis of the raw reads obtained from GEO. This is described in the regulonDB 10.5 paper sections "HT dataset processing" and "Reconstructing regulatory interactions". The processing pipeline, described in fig S1, is fairly typical -- finding motif matches in MACS2 peaks. The cutoff used was simply lower than the source paper, yielding more hits.
Edit: This pipeline description also explains why the tables you download contain duplicates (same regulator, same target). Duplicates can arise from different binding sites upstream of the same gene.


text:
What happened to sra ftp server?
My whole (professional) life I am downloading sequences from SRA and ENA databases. Long time ago, I used prefetch from sra-tools, but later I switched to simply using wget from ftp servers of the two respective servers (used to work for both SRA and ENA). But it seems that SRA ftp server does not have the same data in there anymore and I kind of worry about reproducibility of my old pipelines - so what happened to sra ftp server?
The paths that used to work were formatted as follows
ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/${ACCESION::3}/${ACCESION::6}/"$ACCESION"/"$ACCESION".sra

but it seems sra-instant directory is not accessible anymore. Did the files move? Or were the access been just restricted?
The alternative workaround is to dl the files using prefetch instead of wgetting. Prefetch is part of sra-tools and available on conda and it's not such a bad solution. I guess I am mostly curious what happened to the ftp server and if anyone knows what is coming next, is ENA going to restrict the access too some time soon?

SRA moved all of the data to the cloud per https://ncbiinsights.ncbi.nlm.nih.gov/2020/02/24/sra-cloud/
You can still download data for free using sratoolkit as described in the official documentation here: https://www.ncbi.nlm.nih.gov/sra/docs/sradownload/
You can download the data using sratoolkit without having to pay cloud providers. There are some user-pays data on the cloud but for FASTQ files you don't have to pay.
For documentation on how to download data from the cloud, see https://www.ncbi.nlm.nih.gov/sra/docs/sra-aws-download/

Cool and thanks

This is certainly useful and relevant, but does not truly answers the question. The article does not say anything about discontinuing the ftp servers, rather offering a new way to access the data.

I don't think downloading from the FTP was never an officially advertised method. The entire transition from NCBI FTP to cloud is not a concern for the users of sratoolkit, the preferred method to download data.

NCBI appeared at least in part to shift onto Google Cloud buckets, one previous OP raised a question about access via GCP (I think the question was closed). Fundamentally, traditional FTP is not sustainable against data growtch. GCP is accessed via gsutil. It would be good if you could contact NCBI and post their official response here. They will response when their database provision is in a state of flux.
sratoolkit appears to be a front-end to GCPs gsutil. Just to point out, it would be very unusual for NCBI to charge the user for e.g. direct gsutil download, because usually the provider will charge e.g. NCBI on a per download basis. This will likely work out cheaper overall for NCBI.
Just to formalise my bit, if SRA are hosting on a GCP bucket, Google will not permit ftp access which is at least part the reason behind your question.
They do have their own ftp equivalent called gsutil which will permit direct download if you know the bucket name. gsutil is part of the ClouD SDK package and can be downloaded via conda/anaconda conda install gcloud-sdk. This will give you command line access if you know the bucket names. There is a past post on this BTW.

You don't really need to use gsutil. For example, if you run srapath SRR13397412 you will see the path as https://storage.googleapis.com/nih-sequence-read-archive/run/SRR13397412/SRR13397412 which can be downloaded using wget.

update:  conda install -c conda-forge google-cloud-sdk


text:
Protein-coding gene evolution age and species divergence time: are they correlated?
Is it possible to infer protein-coding gene evolution age by species divergence time?
For example, species A diverged 300 million years (MYA) ago from its common ancestor with species C, whilst species B diverged from C 10 MYA ago.
Assume they share a protein with a similar sequence, i.e. they are  assumed to be homologous. Can I say this homologous gene in species A has an older evolution age than the homologous gene in species B?

Yes this is exactly true, proteins are ideal for deep molecular dating. DNA can suffer from compositional bias. The thing you need is a reliable tree, with strong bootstrap support. You will then calibrate it with your externally established dates.
Programs such as Beast allow you to do the two simultaneously. Generally, whatever approach you take you will be required or else requested to use Beast. It is however impressive.
What are you dating? Ten million years is not 'deep', however 300 million is very deep indeed for animals (Metazoa), especially so for vertebrates.

To answer your question in general terms molecular dating tends to estimate earlier divergence against the fossil record, or against the perceived view of divergence. Molecular phylogenetics argue this is because the fossil record is not sufficiently well sampled, ie the earliest fossil is missing. Palaeontology historically had a different viewpoint.
The plethora of molecular dating methods has continued, an impressive Beast based approach was recently published and I'm not clear how these changes/ improve the dating estimate. My personal view is that generally each new technique improves the accuracy of molecular dating.

Answering the second part of your question, reconstructing the ancestral sequence at a given node is formal analysis and can be performed using ML, Bayes or parsimoney. Thus it attempts to estimate the ancestral sequence and you can then compare the mutations against extant taxa. It is perhaps better to say it establishes a robust hypothesis regarding ancestral mutations.

Thank you for the detailed answer! I am now learning what Beast can do for me.  Btw, then I am curious: is it plausible that a species didn't evolve while some of its protein sequences dramatically evolved? If it is plausible, how could we trust the correlation between divergence time and molecular age?

Answer to first question, absolutely, however modern molecular clocks take this into account (you are thinking of a global clock which very rarely occur) - you need to check but I think Beast has a global clock option. To answer your second question, you should seek compatible external evidence to assess a given clocks accuracy, for example the fossil record which was not used for the molecular date estimate.


text:
Arrange ggplot Figure for scRNA-seq data
I have generated a ggplot for 8 single-cell libraries, with the purpose of visualizing the tSNE facet plot by sample, colored by cell type -- with percentages. The best I could get to is this -

however, it looks too crowded, and I also want the cell types to be labeled for each cluster, which makes it look even worse.
How can I improve the data visualization here?
my code:
p6 <- ggplot(tsnePlot, aes(x=tSNE_1, y=tSNE_2, color=RefinedCellTypes)) +
  scale_color_manual(values=refinedCellTypeColors) +
  geom_point(data = tsnePlotAll, mapping = aes(x=tSNE_1, y=tSNE_2), size = 1, alpha = 1, color = "grey75") +
  geom_point(size = 1, alpha = 1) +
  geom_text(data = tsneClusterLabelsRefined, mapping = aes(x = tSNE_1, y = tSNE_2, label = CellType),color="black",size=0.75,hjust=0) +
  geom_label(data = tsneCellTypeFractions, mapping = aes(x = tSNE_1, y = tSNE_2, label = scales::percent(Freq,accuracy = 0.1)),color="black",fill = "white", alpha = 0.5, label.size = 0) +
  facet_wrap(~Sample) +
  theme_bw(base_size = 12) +
  theme(legend.position="left",
        axis.text=element_text(size=12, color = "black"),
        axis.title=element_text(size=10))

Is there any way I can control the number of cols/rows for the figure?

On top of M__'s suggestion to reduce the point alpha, here are a few other things you can implement to make your figures cleaner:
Remove grid lines from the panel backgrounds
theme(panel.grid.major = element_blank(),
      panel.grid.minor = element_blank())

I think a figure legend like the one you have now is a better approach than labeling every cell type in every facet. However, if you still want to go that route, I suggest implementing repulsive labels with the ggrepel package. geom_label_repel() can also be applied to your percentages so they are not overlapping.

Agreed, BTW I like the user icon

Phage are phantastic!

Except for bacteria :-/

Try
p6 + geom_point(alpha = 1/100)

And variations of 1/10 to 1/1000, this will give a better idea of point density.


text:
Is it possible to do cell composition for scATAC-seq data without scRNA-seq data?
Background My understanding was that if I did scATAC-seq and I have some clusters of cell groups, the only way I can label it is by correlating those groups with scRNA-seq data.
My lab ordered some scATAC-seq data and the vendor did some complimentary analysis and provided cell composition plots indicating which portion of the cells of type X, Y, Z, so on.
Question How is this possible to identity cell composition identifying given cell types? Alternatively is it possible they just used some default RNA-seq data from 10x genomics?

Definitions
scATAC-seq stands for "Single-cell sequencing assay for transposase-accessible chromatin"  This attempts to consider genome-wide regulatory landscapes in single cells. scATAC-seq data can be sparse and noisy presenting a challenge for analysis.


text:
References request: deciding what duplicate samples to include to avoid batch effects with previous studies
Suppose there is cohort with many samples to assay. In phase I of the project, subjects with a certain disease incidence and healthy controls were assayed. Now, in phase II of the project, we are going to assay samples with another disease and analyze them with healthy controls from phase I of the project. From past experience, we know that we have to assay healthy controls from phase I (whose bio samples are still available for assay) to avoid batch effects between phase I and phase II of the project. However, I don't know how much duplicate samples I should include and how the duplicates should be distributed.
My question is certainly related to experiment design, and perhaps there have been extensive discussions on this kind of design, but I don't know the exact term of this kind of experiment design such that I couldn't find anything relevant to my case in the literature. Can anyone point me to references related to this kind of experiment design?

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

at least six biological replicates should be used, rising to at least 12 when it is important to identify SDE genes for all fold changes.

See this paper for an explanation of why six replicates is good for differential expression analysis:
https://dx.doi.org/10.1261/rna.053959.115
When looking at conditions, it's necessary to think about what exactly needs to be compared (ideally prior to doing any experiments or sequencing, because that reduces the cost and frustration of re-doing things).
I recommend you consult with a statistician to find out the appropriate mixture of samples for your experiment, because there can be subtle differences between studies that mean a particular design does not work.
One example from our own investigations is that we were doing cDNA sequencing of cell lines with and without mitochondrial genomes, and it was only after doing two years worth of sequencing and looking at the data that we realised we'd never done any sequencing runs that only involved cell lines without mitochondrial DNA - our rationale for doing that was to reduce batch effects. This meant we couldn't perfectly distinguish between barcode spillover and a very small proportion of mitochondrial DNA in our "mitochondria-free" samples.


text:
Nextflow: how to create many-to-many tuple for process
Using Nextflow, I need to submit x * y jobs where x is the number of input .bam files and y is the number of genome intervals (e.g., <chrom>:<start>-<end>). i.e., for every .bam file, submit a single job for every defined interval.
I have tried a few variations with no luck (see below).
Ultimate question: For each .bam file, how do I submit a single job for every defined genome interval using NextFlow?
Attempts
Attempt 1 (submit both as individual Channels):
Only submits one .bam for all intervals.
sample_bams = Channel.fromPath(sample_input_path + "*.bam")
intervals = Channel.from(['1:10000-20000', '5:55555-77777'])
PROC(sample_bams, intervals)

Attempt 2 (submit as tuple of mapped Channels):
Submits only 10 jobs and fails because the interval Channel is passed as DataflowBroadcast around DataflowStream[?] object. The PROC was changed to receive a tuple rather than individual arguments.
sample_interval_tuples = Channel.fromPath(sample_input_path + "*.bam")
                         .map { sample_file -> tuple(sample_file, align_to_ref, DRF_jar, Channel.from(['1:10000-20000', '5:55555-77777'])) }

Attempt 3 (submit as tuple of mapped .bams to non-Channel of intervals):
Submits 10 jobs and fails because the intervals are passed as a single list. The PROC was changed to receive a tuple rather than individual arguments.
sample_interval_tuples = Channel.fromPath(sample_input_path + "*.bam")
                         .map { sample_file -> tuple(sample_file, align_to_ref, DRF_jar, ['1:10000-20000', '5:55555-77777']) }

Really appreciate your help!
related query: for DSL v2, we wouldn't include from combined_inputs in the input, yeah?

The combine operator can be used to produce the Cartesian product:
sample_bams = Channel.fromPath( './path/to/bams/*.bam' )
intervals = Channel.of( '1:10000-20000', '5:55555-77777' )

sample_bams
    .combine(intervals)
    .set { combined_inputs }

process test {

    echo true

    input:
    tuple path(bam), val(interval) from combined_inputs

    """
    echo -n "${interval} ${bam}"
    """
}

Regarding the related query: Correct. The from and into 'bind' declarations can be omitted

we should crown Steve as King Of Nextflow Answers

I think this cheat sheet might help:
https://github.com/danrlu/Nextflow_cheatsheet/blob/main/nextflow_cheatsheet.pdf

Love the cheat sheet!


text:
Samtools sort: most efficient memory and thread settings for many samples on a cluster
We're preparing to analyze thousands of .bam files beginning with re-alignment, sorting, etc.
Complicated question: Has anyone investigated the optimal thread and memory settings to sort hundreds to thousands of .bam files using samtools sort on a cluster?
By 'optimal', I'm looking to sort as many samples as possible in the shortest time frame. I'm mostly interested in optimizing samtools sort because the other steps (e.g., alignment) are much more obvious to fine-tune.
For a single sample, throwing large amounts of memory + threads is trivial, but that isn't scalable for thousands of samples. There is also likely a diminishing return at a certain point, even for a single sample.
Ideally, I'd like to see plots identifying the 'sweet spot' for threads + memory for a single sample (e.g., x == threads, y == memory, z == time), and then plots maximizing throughput for more samples. I would expect the data to follow a -log function, where we could identify the optimal settings at the knee in the curve.
There are obviously many other factors (e.g., cluster size, administrative scheduler settings, etc.), but I suspect findings on a given cluster will generalize to other clusters.
Really appreciate any feedback.

Update 2
Ok, for what I expect will be my final update for this answer, I compared the following:

samtools 1.14 + zlib
samtools 1.15 + libdeflate (this is a different version, but shouldn't have a major effect. They released a new version just before I asked cluster admins to compile against libdeflate and didn't notice until later.)
sambamba v0.8.2

Note: This time I used the following settings:

Only two samples instead of three
Max Nextflow queue size of 30 to avoid too many threads reading from the same two files.
I tested CPU requests from 1 to 7 (step size 1) and then from 9 to 17 (step size 2).
For jobs where I allotted few CPUs and low memory, I provided a 200% buffer to prevent the jobs from failing with OUT OF MEMORY errors. i.e., I gave the SLURM job 2x the amount of memory than I told the tools they could use. sambamba was the worst offender for exceeding memory in these conditions, but I gave them all 200% buffer to keep comparisons equal.

Seems pretty clear to me that sambamba is the fastest and uses requested CPU resources most efficiently (i.e., is able to use all of the CPUs requested up until ~8 or 9 CPUs). For reference, the difference between sambamba and samtools + libdeflate with 9 CPUs is ~10 minutes on average.
We will certainly be using sambamba going forward, probably with ~9 CPUs and ~9GB per thread since I had to provide a 200% buffer in lower CPU/Mem jobs to avoid job failure, anyway.
Again, all code and results are available on my GitHub repository (Samtools sort optimization test).
(Angle 1) Comparison between samtools + zlib, samtools + libdeflate (Lsamtools), and sambamba
Here, I'm showing the general trend that real time decreases with increased CPUs for all three tools, but Mem Per CPU doesn't have a huge effect. sambamba out performs samtools + zlib regardless of settings. samtools + libdeflate out performs samtools + zlib until ~11 CPUs, where they perform the same (plot 3 below explains why).
Interesting finding: samtools + libdeflate performs better than sambamba with a single thread, but sambamba takes over from there as CPUs increase (again, plot #3 explains why).

(Angle 2) Comparison between samtools + zlib, samtools + libdeflate (Lsamtools), and sambamba
Here, I'm showing that samtools + zlib and sambamba performance flattens out at ~11 CPUs (again, plot 3 below explans why).

Requested CPUs vs. CPU utilization for samtools + zlib, samtools + libdeflate (Lsamtools), and sambamba
Here, I'm showing how well each tool is able to utilize the CPUs allotted to it. It was fascinating to see that samtools + libdeflate flattens out very early in it's CPU utilization (i.e., ~3 CPUS) and struggles to cross the 400% CPU utilization with even >10 CPUs. This is why samtools + libdeflate performance flattens out quickly.
samtools + zlib does a much better job utilizing allotted CPUs, but still lags behind sambamba throughout (and is slower, anyway).

Update 1

My first plots (below) were showing duration on the z-axis, which is the time from job submission to completion rather than job start to completion. That's why the larger memory requests were much longer in the original plots. Duration is still a useful measure, but I'm now including real time, which is from start to completion.
Judging solely by real time (above), there is a clear gain between 5 & 7 CPUs (~10 minutes, on average when compared to the same mem-per-thread).
More data coming, including samtools compiled against libdeflate and sambamba.

Original post
I made time to semi-formally look at this. I'm glad I did because the results surprised me. I don't answer all of my questions directly, but these results are good enough to get me going in the right direction.
My code and results are in my GitHub repository (Samtools_sort_optimization_test) where you will also find an interactive plot you can download to see the results more clearly (see SAMTOOLS_SORT/samtools_sort_CPU_and_memory_comparison.tar.gz in the repository). It's hard to truly appreciate the 3D plotly plot until you start rotating it.
NOTES:

Even though this is from only three samples, I think my results are accurate. I welcome any feedback in case there's a bug.
In my experience, samtools often exceeds its allotted memory (at least on a SLURM scheduler). i.e., the job gets killed by SLURM. To avoid this, I requested an additional 20% of memory to provide a buffer for samtools. I don't know if this discrepancy is because of SLURM or samtools, though.
These results are using samtools compiled against zlib. We are preparing a version of samtools compiled against libdeflate, per @user172818's suggestion. I will add those results once completed. I may also compare against sambamba, per @acvill's suggestion.
These results are related to my other question about the best sorting/merging method. I will post those results in that question and provide a brief summary here.

tl;dr (too long, didn't read)
Primary conclusions based on the plots below

Optimal CPUs (results as expected): There's a significant gain with some additional CPUs, but it levels off around 7 CPUs (total, not 'additional'; i.e., samtools sort -@ 6). I thought we might still see some nice gains above 7, but this is probably what most would expect (e.g., similar to what @user172818 suggested).
Optimal memory per thread (results NOT as expected): There is very little to gain with more than 1GB of memory (per thread); that surprised me. There are some performance gains up until ~7GB per thread (~2-4 minutes faster). What really surprised me, however, is that performance dramatically decreases as memory per thread goes >7GB. The updated plot (above) shows there is some gain with >1GB mem per thread. When running lots of samples at once (hundreds to 1000), mem per thread seems to play a much larger role (data coming).
Effect from sorting status (results as expected): Largely in an attempt to do this as quickly as possible (laziness?), I used three samples, where one was from an aligned and unsorted .bam file, and the other two were from .bam files that had been aligned and sorted, but I then collated them (samtools collate), rather than re-aligning them. Unsurprisingly, the collated samples (not the average use case) sorted slightly faster than the
unsorted sample (average use case).

Initial opinion

"Small" sample set: For a small sample set, I would probably use 7 total CPUs and 7GB per thread.
"Large" sample set: I will probably use 7 CPUs (maybe 5) with only 1GB memory per thread. That will allow the scheduler to keep as many jobs going as possible.

samtools sort: duration vs. CPUs and mem per thread
As mentioned above, an interactive version of this plot is available on my GitHub.
Angle 1

Angle 2

Additional details
Analyses were run on University of Kentucky's MCC Cluster, which uses AMD 7702P processors, using samtools v.1.14

I usually use samtools sort -@4 -m4g. Samtools will be faster with more threads and more RAM, but I feel the gain is not significant. I haven't done a proper evaluation, but I have run samtools sort with different settings from time to time on clusters when I have to consider speed-resource tradeoffs.
When you use samtools, make sure you compile it with libdeflate, which is about twice as fast as zlib on compression. I don't know if the bioconda version built with libdeflate. You may consider the precompiled binary here.
Samtools sort, especially with libdeflate, is fairly fast these days. Performance metrics in the other answer are outdated.

Super helpful info about libdeflate. I was not aware of that. I really want to do a comparison, but it's difficult to make the time. I may do it soon, and will try to include a comparison between zlib and libdeflate.

Side note: if you want to see if the samtools binary you're using is using zlib or libdeflate, you can see which libraries it was compiled against using ldd /path/to/<samtools_binary>. zlib will show as something similar to libz.so.1.


text:
Samtools sort: most efficient approach to sort a single sample after aligning split .fastq files
Related to my other question (Samtools sort: most efficient memory and thread settings for many samples on a cluster), we need to optimize samtools sort as we prepare to analyze thousands of .bam files, beginning with re-alignment, sorting, etc.
Ultimate question: what is the fastest way to sort a single sample after aligning the sample across many separate cluster jobs? i.e., after splitting .fastqs for a single sample and aligning as distinct jobs. I'm looking for recent and clear time/thread/memory comparisons.
Unless I'm missing something, samtools merge was designed exactly for the purpose of merging pre-sorted bams. An excerpt from the manpage (http://www.htslib.org/doc/samtools-merge.html): "Merge multiple sorted alignment files, producing a single sorted output file that contains all the input records and maintains the existing sort order."
Is it fastest to:

samtools sort each 'mini-bam' and then samtools merge?
samtools cat the 'mini-bams' together and then samtools sort?
Another option I haven't considered?

Note: I have seen the following similar questions on other forums, but they are outdated (~8 years ago) and do not answer this question directly:

Best Way To Merge A Many Thousand Small Bam Files Into One Big Bam File?
Merging BAM files before sorting

As part of its operations, samtools sort will split the input up into multiple temporary BAM files (which are each individually sorted) if the reads can't fit entirely within memory. Those files are merged for the final output. From the manpage:

The  sorted output is written to standard output by default, or to the specified file (out.bam) when -o is used. This command will also create temporary files tmpprefix.%d.bam as needed when the entire alignment data cannot fit into  memory (as controlled via the -m option).

Given this, I would expect that a concatenated, unsorted BAM will take more time to process than merging sorted BAM files, because in addition to creating the concatenated BAM file, the splitting would need to be done again (assuming the reads don't fit in memory) prior to the final merge carried out by samtools sort.
I appreciate that this is not a quantified answer; it's just my guess based on what I have observed to happen when doing my own analyses.

It's funny. Now that you mention that, it seems so obvious! Great point. I may still quantify it just to see the difference.


text:
Relationship between entries for same UniProt Id
I was playing around, looking at some AlphaFold predictions for certain proteins and I realized that sometimes the sequences of aminoacid look very different for the same UniProt id. For example, 15C8 from PDB has 2 chains:
DIVLTQSPAIMSASLGERVTMTCTASSSVSSSNLHWYQQKPGSSPKLWIYSTSNLA
SGVPARFSGSGSGTSYSLTISSMEAEDAATYYCHQYHRSPYTFGGGTKLEIKRADA
APTVSIFPPSSEQLTSGGASVVCFLNNFYPKDINVKWKIDGSERQNGVLNSWTDQD
SKDSTYSMSSTLTLTKDEYERHNSYTCEATHKTSTSPIVKSFNRN

and
EVQLQQSGAELVKPGASVKLSCTASGFNIKDTYMHWVKQKPEQGLEWIAQIDPANG
NTKYDPKFQGKATITADTSSNTAYLHLSSLTSEDSAVYYCAADPPYYGHGDYWGQG
TTLTVSSAKTTPPSVYPLAPGSAAQTNSMVTLGCLVKGYFPEPVTVTWNSGSLSSG
VHTFPAVLQSDLYTLSSSVTVPSSTWPSETVTCNVAHPASSTKVDKKIV

At the same time, when looking at the UniProt database (the corresponding id is P01869), there's a list of multiple structures and their 3D shapes (here: https://www.uniprot.org/uniprot/P01869). However, looking at the AlphaFold entry, I see that the amino acid sequence looks different, eg. for the case above it's:
AKTTPPSVYPLAPGSAAQTNSMVTLGCLVKGYFPEPVTVTWNSGSLSSGVHTFPAV
LQSDLYTLSSSVTVPSSPRPSETVTCNVAHPASSTKVDKKIVPRDCGCKPCICTVP
EVSSVFIFPPKPKDVLTITLTPKVTCVVVDISKDDPEVQFSWFVDDVEVHTAQTQP
REEQFNSTFRSVSELPIMHQDWLNGKEFKCRVNSAAFPAPIEKTISKTKGRPKAPQ
VYTIPPPKEQMAKDKVSLTCMITDFFPEDITVEWQWNGQPAENYKNTQPIMNTNGS
YFVYSKLNVQKSNWEAGNTFTCSVLHEGLHNHHTEKSLSHSPGLQLDETCAEAQDG
ELDGLWTTITIFISLFLLSVCYSAAVTLFKVKWIFSSVVELKQTLVPEYKNMIGQA
P

My question (as someone who's just getting started with proteins) would be: what is the relationship between all these entries corresponding to the same UniProt Id?

PDB structure 15C8 corresponds to the protein product of the mouse Igkc gene, with the UniProt accession P01837 as you can see in the linked page:

Your P01869 UniProt accession is for the protein product of the mouse Ighg1 gene instead, so why would you expect these to be the same sequence? They are different proteins coded for by different genes.

Thanks for your response! I got the P01869 after looking it up in this tool: https://www.uniprot.org/uploadlists/ (from PDB to UniProtKB). This returns both P01837 and P01869 -- do you happen to understand why?

And in any case, my confusion regarding the amino acid sequence still holds using the correct UniProt id you mentioned, as AlphaFold has the following sequence: RADAAPTVSIFPPSSEQLTSGGASVVCFLNNFYPKDINVKWKIDGSERQNGVLNSWTDQDSKDSTYSMSSTLTLTKDEYERHNSYTCEATHKTSTSPIVKSFNRNEC which is a substring of the first chain in my original post + "EC" at the end

@Andreea that's... a different question but yes, the structure predictions of AlphaFold or SwissFold or any other such collection tend to be of only parts of a protein and not the entire thing. As for why you got both, that's because P01837 and P01869 are similar proteins and the specific sub-sequence whose structure has been predicted occurs in both of them.


text:
Co-occurrence networks in Metagenomics studies
I have recently acquired some 16S metagenomics data, and was wondering if anyone can speak of the potential limitations, challenges as well as advantages to conducting a network-based study on metagenomics, such as what is done w.r.t co-occurrence networks.
I thought it would be interesting to deviate from the "classic" study of taxonomic abundancies present in different timepoints / conditions. However, it seems as though for such studies they are usually dealing with a relatively large number of samples - please find paper below:
https://animalmicrobiome.biomedcentral.com/articles/10.1186/s42523-020-00054-4
Could such analyses be done with much smaller sample sizes such as those generated via mouse models? What are the potential limitations of such an approach, and how could they be overcome?
Thanks.

Do you mean covariance, ie this is the measurement of co-occurance?

No. There are many possible ways to generate the networks. Based on my understanding, they are generated primarily using correlation metrics, in order to capture 'pairwise affinities', between different taxonomic classifications - i.e. bacteria i, and bacteria j.

I do agree, its just 'co-occurance' isn't a statistical method. Pearson's or Spearman's correlations for example? Pairwise is a statistical method but not a metric and might even suggest distance method like P-distances :-/. Without knowing the metric you run the risk that the assumptions for one, may not be applicable. At present its not clear whether the analysis is applied directly to sequence data - which would likely be a mistake - or their taxonomic designation.

In one of the studies I am currently reading, they use a 'treshold-Kendall correlation'. My understanding is that this is to be applied to processed sequencing data, after taxonomic classification, rather than simply the sequences.

Okay, that makes sense. Would be good to update your question. Its Spearman type stuff which is part of SciPy stats. I dunno enough about R, but there are some very capable R-coders here. Really this is population dyanamics of a microbial community, lactobaccilus (if its gut) and all that. Its definitely not my thing and I wouldn't approach the question in this way.

Would you care expanding on what you don't like about this particular approach? Doesn't have to be a long answer, I just want to get an idea of the limitations of this approach ..

To answer your question more formally. I strongly believe you are treating the mice with something that will shift its microbial flora. Clearly the experiment is not described.
Thus what I simply feel is you essential have a discrete time series comprising before, during and after treatment. 'Treatment' could be anything drug, vaccine, antibiotic, diet who knows.
I think that is a very different time-series to a one where classic correlation statistics will be informative, I'm just guessing in the absence of a hypothesis, thus I could be wrong.
I also kinda suspect lactobacillus are strong players.
Anyway the stat you are looking at is a special case of the general correlation coefficient in much the same way Spearman's classic statistic is. What I can say is that the statistic is not a specialist approach for large sample sizes, which I believe is your question.
'New fangled cool', might be an approach of data mining and that might be okay if you've loads of data you don't know how to analyse, like it goes this way, that way. Ultimately, statistics are used to assess an exact hypothesis. If e.g. lacks power (very common) or the assumptions are violated then look deeper. I do agree that in population biology generating a hypothesis is often helped by data mining, but what I think is you already have an exact hypothesis and I'm far from sure this would nail it.


text:
Best way to match genes in two RNA-seq count matrices from two geo accessions
I have downloaded two datasets, and I am trying to remove the genes they do not share in common.  Is there a biocmanager package with a function that can compare their gene ids, match them and reorder them?

Assuming the datasets are arranged in matrices with one gene per row, this can be done by intersecting row names:
commonNames <- sort(intersect(rownames(data1.mat), rownames(data2.mat)));
filtered.data1.mat <- data1.mat[commonNames,];
filtered.data2.mat <- data2.mat[commonNames,];

It would be helpful if you could provide more information about your datasets. There are many other ways to do this, and the way that works best for you will depend on your particular situation.


text:
Nextflow (DSL v2): how to best synchronize multiple outputs from a single process
I have a workflow that needs to:

Generate .fastq files from .bams (while preserving the @RG group from the original bam)
Split the .fastqs
Align
etc.

The only way I could think to preserve the @RG group was to print it to a file during step 1 and emit the file. I can't figure out how to synchronize the emitted @RG file with the .fastqs, however. I've tried several variations without luck.
Questions:

How do I synchronize the emitted @RG file (sample_RG_file) with the split .fastqs?
Is there a better way?
Bonus question (may be worth a separate post): does Nextflow not stage files until the actual .bash script starts? I tried to operate on the sample_RG_file within the script using Groovy, but got an error that the file didn't exist. (The code was working before I found I needed to synchronize sample_RG_file with the

Other variations:

The code was 'working' before I found I needed to synchronize sample_RG_file with the split .fastqs. I was originally passing GENERATE_FASTQS_PROC.out.sample_RG as a separate parameter to ALIGN_FASTQ_BWA_PROC, but it was always passing the same (first) file.
I also tried GENERATE_FASTQS_PROC.out.sample_RG.flatten() in the .map operator, but that had the same error.

Current error:
The code below fails with the following error:

Error executing process > 'REALIGN_SAMPLES_WF:ALIGN_FASTQ_BWA_PROC (1)'

Caused by:
Not a valid path value type: groovyx.gpars.dataflow.DataflowBroadcast (DataflowBroadcast around DataflowStream[?])

Semi-pseudo code
Workflow
workflow my_WF {
    take:
        sample_input_file_ch
        original_ref

    main:
        GENERATE_FASTQS_PROC(sample_input_file_ch, original_ref)

        reads_per_job = 1_000_000
        fq_split_reads_and_RG_ch = GENERATE_FASTQS_PROC.out.fastqs
                        .splitFastq(by: reads_per_job, pe: true, file: true)
                        .map {
                            fq1, fq2 ->
                                tuple( fq1, fq2, GENERATE_FASTQS_PROC.out.sample_RG )
                        }
                        .view()

        ALIGN_FASTQ_BWA_PROC(fq_split_reads_and_RG_ch, ...)
}

GENERATE_FASTQS_PROC
process GENERATE_FASTQS_PROC {

    input:
        path(sample_input_file)
        val(original_ref)

    output:
        tuple path("**/*_R1.fastq"), path("**/*_R2.fastq"), emit: fastqs
        path("**/*_RG.txt"), emit: sample_RG

    script:
    """
    generate_fastqs.sh $sample_input_file $original_ref $task.cpus $mem_per_thread
    """
}

ALIGN_FASTQ_BWA_PROC
process ALIGN_FASTQ_BWA_PROC {

    /*
     * Publish results. 'mode: copy' will copy the files into the publishDir
     * rather than only making links.
     */
    // publishDir("${params.results_dir}/01-REALIGN_SAMPLES", mode: 'copy')

    label 'REALIGN_SAMPLES'

    input:
        tuple path(fq1), path(fq2), path(sample_RG_file)
        val(align_to_ref)
        val(align_to_ref_tag)
        val(original_ref)
        val(output_format)

    output:
        path '*.{bam,cram}', emit: final_alignment

    script:

    /*
     * Tried to read the RG file using Groovy, but
     * got error that the file didn't exist.
     *
     * Does Nextflow not stage files until the actual .bash script starts?
     */
//    println "sample RG file: $sample_RG_file"
//    file_lines = sample_RG_file.readLines()
//    sample_name = file_lines[0]
//    sample_RG = file_lines[1]

    """
    lines=()
    while read line
    do
        lines+=("\\\$line")
    done < $sample_RG_file

    sample_name="\\\${lines[0]}"
    sample_RG="\\\${lines[1]}"

    bash realign_bwa-split.sh $fq1 $fq2 \\\$sample_name \"\\\$sample_RG\" ${align_to_ref} ${output_format} $task.cpus
    """
}

Update
Per @Steve's suggestion, I'm trying to have GENERATE_FASTQS_PROC create the split .fastq files using Unix split, but I'm still struggling with various concepts.
Supposed I have this emitted from my process:
[sample_A_RG.txt, [s1_R1.split_00.fq, s1_R2.split_00.fq, ..., s1_R1.split_XX.fq, s1_R2.split_XX.fq]]`

In my mind, I want to pull out s1_RG.txt file from the Channel and extract the @RG. Then I want to call:
[<all_fastqs_for_sample_A>]
    .map{
        // pull out split ID numbers
    }
    .groupTuple()

I believe this would group _R1 and _R2 for each split set, and then call
[[00.fq1, 00.fq2], ..., [XX.fq1, XX.fq2]]
    .combine(sample_RG_var)

to get [[@RG, 00.fq1, 00.fq2], ..., [@RG, XX.fq1, XX.fq2]].
Seems like it would be so much easier if I could use Channel factories and operators (e.g., .fromFilePairs().combine( rg )) in the output:.
Getting files into a process is relatively straightforward, but getting them out in a reasonable manner is less so.

Nextflow channels guarantee that items are delivered in the same order as they are sent. So if a process declares two (or more) output channels, the items that are emitted should already be synchronized. This is true unless of course one of the channels sets the optional output attribute.

How do I synchronize the emitted @RG file (sample_RG_file) with the
split .fastqs?

Since the splitFastq operator can create n outputs for each input, you need some way to combine the @RG files using some key. This key, for example, could be the basename of the sample input file. This, I think, would work nicely given the input file should just have a '.bam' extension. Another option, would be to supply some other value in the input declaration and just use that in your output declaration. If we went with the former, your output declaration might look like:
output:
tuple val(sample_input_file.baseName), path("sample_R1.fastq"), path("sample_R2.fastq"), \
    emit: fastqs
tuple val(sample_input_file.baseName), path("sample_RG.txt"), \
    emit: sample_RG

Then, the body of your workflow might look like:
workflow {

    reads = GENERATE_FASTQS_PROC(sample_input_file_ch, original_ref)

    reads.fastqs
        .splitFastq( by: params.reads_per_job, pe: true, file: true )
        .combine( reads.sample_RG, by: 0 )
        .map { sample, fq1, fq2, rg_txt -> tuple(fq1, fq2, rg_txt) }
        .set { fq_split_reads_and_RG_ch }

    ...
}

You'll notice that I removed your fuzzy glob patterns from your output declaration above. Your workflow definition suggests that the generate_fastqs.sh script produces a single pair of FASTQ files and a single readgroup file for a given set of inputs. However, as you know, glob patterns can match one or more output files. When declaring output files, it's always best to name these guys explicitly. If you do use a glob pattern and for some reason your process was to produce two files that match a given pattern, you'll instead get back a tuple, which would break the pipeline.

Is there a better way?

Since you need the @RG tags with the split FASTQ files, I don't think there is. However, I generally avoid using the splitting operators on large files. My preference is always to pass "work" off (including IO) to a separate process. A separate process can also provide increased flexibility if/when requirements change.

Does Nextflow not stage files until the actual .bash script starts? I tried to operate
on the sample_RG_file within the script using Groovy, but got an error
that the file didn't exist.

Nextflow evaluates the code in the script definition before returning the script string to be executed. This is why it's possible to implement a conditional script. Basically, you want to avoid all IO operations here. The script string to be executed has not even been decided at this point.

Update
I'm not sure how you plan to operate per sample. IIUC, [<all_fastqs_for_sample_A>] is basically a Channel of lists, where each list contains all split FASTQ files for a given sample. Using map here to pull out the split ID numbers followed by the groupTuple operator might not give the expected results. You might also run into the problems with calling groupTuple, since, without a size attribute, it will basically wait for all outputs before emitting the grouped result:

You should always specify the number of expected elements in each
tuple using the size attribute to allow the groupTuple operator to
stream the collected values as soon as possible. However there are use
cases in which each tuple has a different size depending on the
grouping key. In this case use the built-in function groupKey that
allows you to create a special grouping key object to which it’s
possible to associate the group size for a given key.

Part of the problem is that we don't know what the generate_fastqs.sh script does and how the other scripts work. If you are aligning using BWA MEM and have samtools available in your $PATH, it might even be easier to produce interleaved output (since BWA MEM can read interleaved FASTQs) using the samtools fastq command. This would save some work having to group each list of FASTQs into pairs, ensuring they are ordered as expected. This would be reasonably involved.
Not sure what BAMs you need to process, or if you're looking for a more generic solution, but I thought I'd throw in an alternative solution you might find works well for you or that you may not have yet considered. Since BAMs often contain multiple readgroups, I find that if I split by readgroup (using samtools split, flatten the output Channel, then samtools collate followed by samtools fastq, I can achieve somewhat uniform mapping walltimes (assuming the input BAMs are similar) without needing to chunk the reads and then merge the alignments back together. If your input BAMs are not similar at all, then you will likely need to chunk to generate uniform alignment jobs.
Below is an example, that shows how to create a groupKey, followed by the transpose, join/combine and groupTuple pattern. There are some caveats and there is some room for optimization:
nextflow.enable.dsl=2

params.original_ref = './ref.fasta'
params.reads_per_job = 10000

process samtools_view_header {

    tag { "${sample_name}:${bam.name}" }

    cpus 1
    memory 1.GB

    input:
    tuple val(sample_name), path(bam)

    output:
    tuple val(sample_name), path("${bam.baseName}.header.txt")

    """
    samtools view \\
        -H \\
        -o "${bam.baseName}.header.txt" \\
        "${bam}"
    """
}

process samtools_name_sort {

    tag { "${sample_name}:${bam.name}" }

    cpus 4
    memory 8.GB

    input:
    tuple val(sample_name), path(bam)

    output:
    tuple val(sample_name), path("${bam.baseName}.nsorted.bam")

    script:
    def avail_mem = task.memory ? task.memory.toGiga().intdiv(task.cpus) : 0
    def mem_per_thread = avail_mem ? "-m ${avail_mem}G" : ''

    """
    samtools sort \\
        -@ "${task.cpus - 1}" \\
        ${mem_per_thread} \\
        -n \\
        -o "${bam.baseName}.nsorted.bam" \\
        -T "${bam.baseName}.nsorted" \\
        "${bam}"
    """
}

process samtools_fastq {

    tag { "${sample_name}:${bam.name}" }

    cpus 1
    memory 2.GB

    input:
    tuple val(sample_name), path(bam)

    output:
    tuple val(sample_name), path("${bam.baseName}.fastq")

    """
    samtools fastq \\
        -O \\
        -T RG,BC \\
        -0 /dev/null \\
        "${bam}" \\
        > \\
        "${bam.baseName}.fastq"
    """
}

process split_fastq {

    tag { "${sample_name}:${fastq.name}" }

    cpus 1
    memory 1.GB

    input:
    tuple val(sample_name), path(fastq)

    output:
    tuple val(sample_name), path("${fastq.baseName}.${/[0-9]/*5}.fastq")

    """
    split \\
        --suffix-length=5 \\
        --additional-suffix=".fastq" \\
        -d \\
        -l "${params.reads_per_job*4}" \\
        "${fastq}" \\
        "${fastq.baseName}."
    """
}

process bwa_index {

    tag { fasta.name }

    cpus 1
    memory 12.GB

    input:
    path fasta

    output:
    tuple val(fasta.name), path("${fasta}.{amb,ann,bwt,pac,sa}")

    """
    bwa index "${fasta}"
    """
}

process bwa_mem {

    tag { "${sample_name}:${fastq.name}" }

    cpus 8
    memory 12.GB

    input:
    tuple val(sample_name), path(fastq), path(header)
    tuple val(idxbase), path(bwa_index)

    output:
    tuple val(sample_name), path("${fastq.baseName}.bam")

    script:
    def task_cpus = task.cpus > 1 ? task.cpus - 1 : task.cpus

    """
    bwa mem \\
        -p \\
        -t ${task_cpus} \\
        -C \\
        -H <(grep "^@RG" "${header}") \\
        "${idxbase}" \\
        "${fastq}" |
    samtools view \\
        -1 \\
        -o "${fastq.baseName}.bam" \\
        -
    """
}

process samtools_coordinate_sort {

    tag { "${sample_name}:${bam.name}" }

    cpus 4
    memory 8.GB

    input:
    tuple val(sample_name), path(bam)

    output:
    tuple val(sample_name), path("${bam.baseName}.csorted.bam")

    script:
    def avail_mem = task.memory ? task.memory.toGiga().intdiv(task.cpus) : 0
    def mem_per_thread = avail_mem ? "-m ${avail_mem}G" : ''

    """
    samtools sort \\
        -@ "${task.cpus - 1}" \\
        ${mem_per_thread} \\
        -o "${bam.baseName}.csorted.bam" \\
        -T "${bam.baseName}.csorted" \\
        "${bam}"
    """
}

process samtools_merge {

    tag { sample_name }

    cpus 1
    memory 4.GB

    input:
    tuple val(sample_name), path(bams)

    output:
    tuple val(sample_name), path("${sample_name}.bam{,.bai}")

    """
    samtools merge \\
        -c \\
        -p \\
        "${sample_name}.bam" \\
        ${bams}
    samtools index \\
        "${sample_name}.bam"
    """
}

workflow {

    original_ref = file( params.original_ref )
    bwa_index( original_ref )

    Channel.fromPath('./data/*.bam')
        | map { tuple( it.baseName, it ) }
        | set { sample_bam_files }

    samtools_view_header( sample_bam_files )

    sample_bam_files
        | samtools_name_sort
        | samtools_fastq
        | split_fastq
        | map { sample_name, fastq_files ->
            tuple( groupKey(sample_name, fastq_files.size()), fastq_files )
        }
        | transpose()
        | combine( samtools_view_header.out, by: 0 )
        | set { realignment_inputs }

    bwa_mem( realignment_inputs, bwa_index.out )
        | samtools_coordinate_sort
        | groupTuple()
        | map { sample_name, bam_files ->
            tuple( sample_name.toString(), bam_files )
        }
        | samtools_merge
        | view()
}

There are some caveats however, in no particular order:

It unravels your shell scripts. I think this is a better solution though because it gives you control over each of the components.

No CRAM support. But this should be easy to add if needed.

A single process to read a BAM header is wasteful, especially if the input BAM is large. More time will be spent localizing the file than actually doing work.

Replace samtools sort -n with samtools collate. There's really no need to name sort. You may prefer to pipe the output of samtools collate into samtools fastq for avoid writing an extra file. This might be useful if you need to operate on an extremely large number of large BAMs.

The 'samtools_fastq' step should write gzipped output with low compression. This means replacing coreutils split with something else. A solution using Biopython could be useful here to read and write gzipped output in a single pass.

The last 'map' command in the workflow is unnecessary. I just thought I'd show that calling .toString() on the GroupKey gives you back an actual String object.

Also, it's not necessary to pass around a tuple of 'sample name' and paths as process inputs and outputs. There's no harm in doing this though. When you are more confident with Nextflow Channels, refactoring these processes should be easy.

My samtools sort processes ensure that they do not oversubscribe the nodes that they land on. They use approximately 320% - 380% CPU, depending on the size of the BAM. The resource requests for the other processes are not optimized at all.

With newer versions of samtools, BAM indexing can be done at the same time as merging, with the --write-index option.

Thanks, @Steve. You've made several important points. Question about your suggested workflow code: wouldn't .splitFastq fail since reads.fastqs contains a tuple including the sample_file_input.baseName? I've been battling how to work with the tuples of tuples in various scenarios (e.g., in this situation and after a .groupTuple()).

Also, per your suggestion, I've been working to have GENERATE_FASTQS_PROC split the .fastqs itself and emit those. That's causing me a headache. The level of abstraction in Nextflow has really made these details quite elusive for me.

I added an update to my question explaining what I feel like I should be able to do and what I ultimately want to get.

Thanks, @Steve. I think the interleaved option might be the simplest. Great suggestion. Still, I'm shocked by how difficult it is to get files out of a process without mixing them between samples. What kills me is, I know all files for a single sample are safely contained in a unique directory, but I'm forced to mix them up with files from other samples in the same Channel, so I'm left with convoluted methods to try to distinguish them from each other. If I emit the path to the directories and then use Channel.fromPath(), it would be sooooo easy. But that's not possible, AFAIK.

generate_fastqs.sh takes a single .bam -> sorts by name -> converts to .fastq (currently w/ bedtools), and splits them with Unix split. The files are stored safe and sound, all tidy-like in a unique directory, where I have clear control over what's generated and mixed together (unlike Nextflow Channels). Forgive the frustration. :-)

Sorry @MarkEbbert, I was meant to let you know yesterday that I had updated my answer but got pulled in another direction before I could finish the update. WRT your first question above, splitFastq will work fine with a tuple that includes a sample identifier. For example, it will work with the output of fromFilePairs with flat: true.

@MarkEbbert I struggled with the same thing when I first started with NF. I found that outputting tuples of files, prefixed with some identifier (usually a sample name) helps considerably when working with multiple channels.

This was a huge help, @Steve. You've also made several additional helpful suggestions. It had never dawned on me that samtools was oversubscribing because -@ is for additional threads, though that should have been completely obvious. Really appreciate your help. I'm still not a fan of NF output Channels, but maybe I'll get there. :-)

Update 3
While the code in Update 2 works, @Steve pointed out that my approach may not be ideal in cloud environments (see comments). He also provided a Nextflow solution that works. I'm still not a fan of how Nextflow output: Channels work, but that may be from my inexperience.
Update 2
Pretty sure I got it this time. The following code works. I'm getting all of the tuples I was expecting.
Here's the updated code for my workflow to create the tuples from each line in the sample_tuple_file I created from a bash for loop. I spent quite a while trying to get a combination of .splitText() and .splitCSV() to work, but for some reason, you cannot chain as follows: it.splitText().splitCSV().
     GENERATE_FASTQS_PROC.out.sample_tuple_file
         .map {
            it.splitText()
         }
         .flatten()
         .map {
             toks = it.trim().split(',')
             sample_name = toks[0]
             fq = toks[1]
             rg = toks[2]
             tuple(sample_name, fq, rg)
             }
         .view()
        .set { sample_fq_tuples_ch }

Update 1
I thought I had it, but I'm actually only getting ONE of the tuples from my GENERATE_FASTQS_PROC.out.sample_tuple_file.map { } statement (the last one). And I have no idea why...

Original posted answer
Despite @Steve's exceptional help and patience with me, I have really battled to figure out how to work with the Nextflow output: Channels, but simply could not figure out any solution using them.
I finally figured out a super simple method that is much easier than dealing with Nextflow Channels, which I share here in case it will help others. I personally will probably use this method from now on, unless someone can convince me that's a bad idea. In my view, using my approach gives me complete control over the way outputs and inputs are handled, and I don't have to worry whether Nextflow is doing what I hope it's doing.
generate_fastqs.sh
generate_fastqs.sh takes a single .bam -> sorts by name -> converts to .fastq (now w/ samtools fastq, thanks to @Steve), and splits them with Unix split. Then, at the end of the script, I create a single file with all of the tuples that I want to pass to my next process.
Here are the key lines for these steps (including exit codes to ensure Nextflow receives error notifications). The last step is the key (for me, anyway):
...

#####################
# Sort .bam by name #
#####################
if ! time samtools sort -@ $n_threads -n -m ${mem_per_thread}G $sample_input > $sample_sorted_by_name; then
    echo "ERROR: sort failed"
    exit 1
fi

...

#############################
# Create interleaved .fastq #
#############################
if ! time samtools fastq -o $fq -s /dev/null -0 /dev/null --reference $original_ref $sample_sorted_by_name; then
    echo "ERROR: bamtofastq failed"
    exit 1
fi

...

################
# Split fastqs #
################
if ! time split --additional-suffix=".fastq" -dl $n_lines $fq "${split_dir}/${sample_name}.interleaved_R1_R2.split."; then
    echo "ERROR: split failed"
    exit 1
fi

...

###########################################################
# Store sample name, .fastq, and RG to file for alignment #
###########################################################
RG_file="${sample_name}.tuples.txt"

# Create empty file
> $RG_file
for fq_file in $PWD/${split_dir}/*
do
echo "${sample_name},$fq_file,$RG" >> $RG_file    # <-- Creating future tuples in file, here (comma separated).
done

GENERATE_FASTQS_PROC
In GENERATE_FASTQS_PROC, emit one simple file--the tuple file I created myself with a simple bash loop (so easy and feels so good):
    emit:
        path("*.tuples.txt"), emit: sample_tuple_file

Workflow
In my workflow, I call GENERATE_FASTQS_PROC, access the file from the Channel, read each line, and create a Nextflow tuple for each line (again, the lines I so easily created in a bash for loop).
    /*
     * Generate fastq files from the original .(cr|b)am file and split into
     * mini .fastq files with `n_reads_per_run` in each mini .fastq
     */
    GENERATE_FASTQS_PROC(sample_input_file_ch, n_reads_per_run, original_ref)

    /*
     * Get sample name, fastq, and read group (RG) from the tuple file using Nextflow
     */

     GENERATE_FASTQS_PROC.out.sample_tuple_file
        .map{
            sample_tuple_file ->
                sample_tuple_file.withReader {          // Open file
                    it.eachLine {                       // Loop over each line in file
                        line ->
                            toks = line.split(',')      // Split by comma
                            sample_name = toks[0]       // Get each item
                            fq = toks[1]
                            rg = toks[2]

                            tuple(sample_name, fq, rg)  // Create the simple tuple that has eluded me for two days!
                    }
                }
        }
        .set { sample_fq_tuples_ch }                    // Assign it to a new Channel name and move along!

    /*
     * Align fastq file pairs
     */

    ALIGN_FASTQ_BWA_PROC(sample_fq_tuples_ch, align_to_ref,
                             align_to_ref_tag, original_ref, output_format)
```

Ugh. I thought I had it, but I confused myself with my own print statements. I'm actually only getting ONE of the tuples from my GENERATE_FASTQS_PROC.out.sample_tuple_file.map { } statement. And I have no idea why...

Pretty sure I got it working with update #2. Hopefully others can learn from my mistakes and final victory.

I can see you have put in quite a bit of effort here, but this is not something I would recommend to others. This is a solution that will work on a local filesystem only, where all files exist under a common root. This is what makes it possible to pick them up later on if you have a pointer to them (your tuples.txt contains file pointers). Since the required files are not declared in the output block, Nextflow will not attempt to "unstage" them (upload them to a remote object store, like a bucket) when you deploy your workflow to the cloud. They will just be lost.

Thanks for including your generate_fastqs.sh. I have a better picture now of what you are trying to do. I will try to update my answer with a working solution soon. I just need to grab a couple of test BAMs. I should have a couple somewhere..

True. I hadn’t considered cloud file systems. I was worried you’d burst my bubble. :-)

Couldn’t I still emit the files you mention, but use my method to generate tuples, etc? Or, couldn’t I publish results and access them that way? I realize this breaks convention, but it appears that, while Nextflow has many advantages,  Nextflow’s output Channels are fundamentally convoluted, flawed, and maybe even dangerous. I’d love to be proven wrong. I don’t know Nextflow’s internals, but seems like this could be easily fixed. e.g., if every process instance assigned an ID to the process that we could filter by upon return, I could trust that NF is giving me what I think I’m getting.

You could perhaps. But that would not be a pleasant experience at all and would produce a convoluted solution. I don't have any real problem with Nextflow Channels. I think they're great and they make writing complex workflows possible. I've updated my answer above. I think it's ok as an example but there are some caveats listed. You may or may not need to address these. But it will likely need some tinkering before you can run with it.


text:
Imputing 23andMe 3v data to v4 and v5?
The 23andMe version 3 data (v3) is now quite 'old' and new data is being collected on the updated v5 Illumina SNP-Chip (details here).
Is there an online service for imputing 23andMe v3 files to v4 and v5?
Seems this would be quite handy!

I am seriously not sure why is this getting downvoted, I find the post-edit question a clear cut reasonable

Thanks Kamil, I was surprised... Good to know it's not just me :-)

The different V numbers are just different number of SNPs as far as I'm aware. So there's not going to be a service which just imputes the variants to make e.g. V3 into V5.
The Michigan Imputation Server uses the Haplotype Reference Consortium to impute genotypes and allows input in 23andme format. You'll end up with many more SNPs than in the V5 chip, but filtering out the other ones should be easy.

Thanks, I see what you mean. I'm now curious what coverage of v5 you get from v3 :-)