# Question
Title: 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?

# Answer
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.
> 5 votes

---
Tags: pdb
---


# Question
Title: How to compare across Microarray datasets for a gene in Mice?
How does one go about comparing expression levels of one gene across multiple microarray datasets pulled from different experiments/papers in mice? In my example, I've downloaded the expression levels for my gene of interest from different microarray datasets for different cell populations of interest (specifically healthy hematopoietic cells). How would I pool the datasets together?

Say I have common cell types from multiple datasets such as B-cells, would I be able to just combine them into one data frame and visualize them as one or would I need to normalize some aspect of that?

I'm really confused because I've had different people tell me different things.

# Answer
> 1 votes
If you are after expression of tissue-specific, but otherwise healthy wild-type mice, Bgee is really a great resource for you. They pull a lot of data, curate it and sort it in a single database that easy to browse. They also have lots of tools for comparative genomics (if you would be wondering about expression of an ortholog in a different species).

If you are just after the methodology of pulling expression data together from multiple sources be because you are interested in pulling data from specific papers, you can check their paper about the database. The pipeline is available on GitHub too.

---
Tags: gene-expression, microarray
---


# Question
Title: How to get the number of complete phage genomes available on ncbi?
I am looking to establish the total number of complete phage genomes available on NCBI.

I am not looking for any specific type, but want to understand the total diversity available.

How can I perform the search?

# Answer
1. Go to the NCBI nuccore database, at https://www.ncbi.nlm.nih.gov/nuccore Copy and paste the search term:

`(bacterial virus[Organism]) AND (3000:110000[sequence length]) NOT bacteria[Organism]`

You'll get 74583 results.

2. Click on the link `send to` and click the checkbox called `file`. It will then download.
3. Check the downloaded file is in `fasta format` with "\>" to start the ID.
4. If fasta format is absent, return to the webpage and where it states `summary` switch that to `fasta text`

I can't remember if point 4 is needed.

**Notes**

* Sequence length, I am guessing that all phages are between 3000 and 105000 in length, this is to screen out partial sequences and to minimise bacteria creeping in to the search. I don't quite understand how a phage can be 105 000bp in length, but you can twiddle these parameters to refine the search.
* @acvill informs us that giant phages greater than \>200 kb have been annoated. Thus a separate search is likely required targeting this genome size.
> 5 votes

# Answer
How about trying this?

1. Go to https://www.ncbi.nlm.nih.gov/nucleotide. Then type the following command in the search box. It resulted in 23971 items.

`(bacterial virus[Organism]) AND (complete genome[ti]) NOT shotgun[ti] NOT plasmid[ti] NOT bacteria[Organism]`
> 2 votes

---
Tags: genome, public-databases, phylogenetics, virus
---


# Question
Title: What are values of FILTER column of vcf files produced by Mutect2
I have called SNVs with Mutect2, now in filter column of a file called filtered vcf I have a lot of things like

```
> length(unique(filtered.vcf$FILTER)) [1] 409 >

> unique(filtered.vcf $FILTER)
  [1] "germline"
  [2] "PASS"
  [3] "clustered_events"
  [4] "clustered_events;panel_of_normals"
  [5] "weak_evidence"
  [6] "haplotype"
  [7] "base_qual;haplotype"
  [8] "clustered_events;germline;haplotype"
  [9] "base_qual;clustered_events;germline;haplotype"
 [10] "clustered_events;haplotype"
 [11] "clustered_events;haplotype;weak_evidence"
 [12] "base_qual;clustered_events;haplotype"
 [13] "haplotype;map_qual"
 [14] "germline;panel_of_normals"
 [15] "base_qual"
 [16] "germline;haplotype"
 [17] "germline;panel_of_normals;position"
 [18] "clustered_events;orientation"
 [19] "clustered_events;orientation;strand_bias"
 [20] "base_qual;haplotype;map_qual"
 [21] "clustered_events;fragment;haplotype"
 [22] "base_qual;clustered_events"
 [23] "germline;haplotype;panel_of_normals"
 [24] "clustered_events;germline;haplotype;panel_of_normals"
 [25] "fragment"
 [26] "clustered_events;haplotype;map_qual"
 [27] "clustered_events;germline"
 [28] "clustered_events;germline;panel_of_normals"
 [29] "germline;weak_evidence"
 [30] "orientation"

```

I tought that the `filtered.vcf` Mutect2 output file contain only the filtered events, but the values I see suggest otherwise. How should I understand the output?

# Answer
You haven't told us what commands you ran, so I am assuming you first ran `Mutect2` and then `FilterMutectCalls`. If so, the `FilterMutectCalls` does not remove any variants from its input, it will simply flag any that fail its builtin filters. Any variants that did not fail a filter are marked as PASS and the rest have a filter added.

Now, whether you want to filter these out or not is up to you. You need to look at the VCF header where each of the filters is defined (admittedly, not very clearly) and make a judgment call on whether you consider this filter to be a deal breaker or not.

For instance, the `germline` flag means:

```
##FILTER=<ID=germline,Description="Evidence indicates this site is germline, not somatic">

```

So that indicates that gatk believes this is more likely to be a germline variant than a somatic one, that it is more likely to be a constituent variant present in the sample's germline cells than a variant that is only found in the sample's tumor cells.

Another clear case is the `base_qual` filter which means:

```
##FILTER=<ID=base_qual,Description="alt median base quality">

```

When a variant is flagged with this filter, it means that the median fastq base quality of the reads supporting the variant was significantly lower than the median fastq base quality of the reads supporting the reference alele at the same position, making it likely that the variant is a false positive.

You need to read the VCF headers and any gatk documentation you can find (warning: these filters are not very well documented at all, in my experience), understand what the filters are and then decide what variants *you* consider real based on what you know about your sample, your experimental design and the question you are trying to answer.
> 5 votes

---
Tags: vcf, filtering, snv, wes, point-mutation
---


title: How to rescale a function?

text:
# Question
Title: 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?

# Answer
> 2 votes
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}$

---
Tags: sequence-analysis
---


title: Python argparse error

text:
# Question
Title: 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

```

# Answer
> 1 votes
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.

---
Tags: python, bash
---


# Question
Title: 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.

# Answer
Small modification,

```
find . -type f -exec md5sum-lite {} \; >> output.txt

```

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

---
Tags: genome, sequence-analysis, illumina
---


title: pathogeniticy prediction tool

text:
# Question
Title: 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

# Answer
> 3 votes
## 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.

---
Tags: genetics, variants, pathogenicity
---


# Question
Title: Correcting for noise in RT-qPCR gene expression data
I have a training set of RT-qPCR gene expression data (not run in triplicate) for a batch of samples with two phenotypes $A$ and $B$ on which I've trained a "logistic regression classifier".

I also have another smaller set of samples which have been run twice as technical replicates on the RT-qPCR machine in order to study the effect of qPCR noise on my classifier.

I have a particular cut-off probability $P$ for my logistic regression model, such that we must be at least $P$ confident to classify a test sample as $A$ (similarly for $B$), otherwise no classification is returned.

I would like to incorporate the effect of RT-qPCR noise into my logistic classifier in order to return a more accurate classification probability for a particular test sample.

My model for the RT-qPCR noise for a particular sample $i$ is given by

$$X\_i\sim N(\mu\_i, \sigma^2I).$$

The logit is then given by

$$t\_i = \beta\_0 + \beta^TX\_i,$$

with distribution

$$t\_i\sim N(\beta\_0+\beta^T\mu\_i, \sigma\_t^2),$$

where $\sigma\_t^2=\\|\beta\\|\_2^2\sigma^2$.

Since $\mu\_i$ is a nuisance parameter, we eliminate it by taking the difference of the two replicates for a particular sample

$$Z\_i=t\_i^{(1)} - t\_i^{(2)} \sim N(0, 2\sigma\_t^2).$$

We can then form the MLE of $\sigma\_t^2$ by computing

$$\hat{\sigma}\_t^2=\frac{1}{2N}\sum\_{i=1}^NZ\_i^2.$$

Therefore, for a particular test sample $x\_i$, we compute the expected value of the logit-normal distribution with location $t\_i$ and squared scale $\hat{\sigma}\_t^2$

$$p\_i=\int\_0^1x\frac{1}{\sqrt{2\pi\hat{\sigma}\_t^2}}\frac{1}{x(1-x)}e^{-\frac{(\text{logit}(x)-t\_i)^2}{2\hat{\sigma}\_t^2}}dx.$$

The resulting $p\_i$ is then our corrected probability of phenotype $A$.

*Is this a sensible way to achieve my stated goal?*

Here is an example plot of the necessary correction when the squared scale is $\sigma^2=2$,

* The black line corresponds to a qPCR assay with no noise.
* The violet line indicates the corrected probability.

For example, if the regression model says that a sample has a 75% chance of being phenotype $A$, then in the corrected model, this actually corresponds to only having about a 65% chance of being phenotype $A$.

# Answer
You're actually not correct about the functional form of the logit. It's nonlinear, but you've written a linear function.

That matters here because once you understand the functional form of logistic regression, you can write it in a way that relates it to your specific generative model for your data. That can tell you what features you ought to use if your model assumptions are correct. You can find the details in PRML section 4.2 (pdf). The essence is:

* If your class labels are generated from a mixture of two classes with prior probabilities $\pi\_1$ and $\pi\_2$
* If your covariates are generated from exponential-family distributions with densities proportional to $\exp(\beta\_1^T\Phi(x))$ and $\exp(\beta\_2^T\Phi(x))$

then the posterior probability that an observation $x$ is in class 1 is

$$\frac{\pi\_1\exp(\beta\_1^T\Phi(x)) }{\pi\_1\exp(\beta\_1^T\Phi(x)) + \pi\_2\exp(\beta\_2^T\Phi(x))$}$$. If you plug in known class labels and mangle this for a while with algebra, you end up with a logistic regression using features $\Phi(x)$. The intercept term of the logistic regression turns into a function of $\pi\_1/\pi\_2$.

If your data are Gaussian, then $\Phi(x)$ would generally include $x$ and $x^2$ and all the cross-terms if $x$ is multivariate. If there are special assumptions about equal covariance between classes, then the squared terms and cross terms may not be needed.
> 1 votes

---
Tags: gene, statistics, pcr, noise, gene-expression
---


# Question
Title: 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?

# Answer
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.
> 3 votes

---
Tags: phylogenetics, proteins, phylogeny
---


title: What happened to sra ftp server?

text:
# Question
Title: 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?

# Answer
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/
> 3 votes

# Answer
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.
> 2 votes

---
Tags: public-databases, data-download, sratoolkit, sra
---


# Question
Title: 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?

# Answer
> 1 votes
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.

# Answer
> 0 votes
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.

---
Tags: rna-seq, gene-expression, microarray
---


# Question
Title: 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?

# Answer
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.
> 3 votes

# Answer
Try

```
p6 + geom_point(alpha = 1/100)

```

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

---
Tags: scrnaseq, single-cell, visualization, ggplot2
---


# Question
Title: 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.

# Answer
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.
> 2 votes

---
Tags: chip-seq
---


# Question
Title: 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?

# Answer
> 0 votes
> 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.

---
Tags: batch-effects, design
---


# Question
Title: 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 `Channel`s):**

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 `.bam`s 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?

# Answer
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
> 4 votes

# Answer
I think this cheat sheet might help: https://github.com/danrlu/Nextflow\_cheatsheet/blob/main/nextflow\_cheatsheet.pdf
> 2 votes

---
Tags: nextflow
---


# Question
Title: 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?

# Answer
> 1 votes
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.

---
Tags: proteins, protein-structure, pdb, uniprot
---


# Question
Title: 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?

# Answer
> 3 votes
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.

---
Tags: rna-seq, bioconductor
---


# Question
Title: 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 `.fastq`s 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:

1. `samtools sort` each 'mini-bam' and then `samtools merge`?
2. `samtools cat` the 'mini-bams' together and then `samtools sort`?
3. 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:

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

# Answer
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.
> 1 votes

---
Tags: sequence-alignment, sequence-analysis, samtools
---


# Question
Title: 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.

# Answer
> 2 votes
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.

---
Tags: phylogenetics, phylogeny, metagenome, 16rrna, microbiome
---


# Question
Title: Extract sequence context of high-degree nodes in assembly graphs
I often use metaSPAdes to assemble short reads from human microbiomes. My simplified understanding of short-read de Bruijn graph assemblers is that they fail where ambiguous paths cannot be resolved. While it can be said that these points of failure may be due to certain **structural** features of the underlying sequences -- long repetitive sequences or highly conserved sequences found at multiple loci -- I'm interested in the **functional** classification of sequences that are proximal to nodes on assembly graphs. Put another way, what kinds of genes are poorly represented in metagenomes because they are poorly assembled? Naively, I expect mobile elements like transposons and integrons to be in such a list of genes, since they are often flanked by repeated sequences and can be present in many copies across multiple genomic contexts, though I want to take an empirical approach.

I've been dabbling with Bandage, which allows you navigate assembly graphs in a GUI, though I'm having trouble extracting features of the underlying network that Bandage builds from the input (.fastg or .gfa). I'd like to know if there are any tools or workflows that would allow me to, given an assembly graph in GFA format ...

1. filter nodes by degree,
2. extract sequences proximal to high-degree nodes, and
3. filter pulled sequences based on sequencing depth, so as to preclude analyzing node-adjacent sequences in the case where the assembler failed due to low depth.

I'm open to the possibility that I'm taking the wrong approach to this problem, so feel free to answer with a frame challenge.

# Answer
> 4 votes
I don't have a recommendation of specific tools, but we might be able to use the GFA file itself to get at this with just bash commands.

For these examples I use the linked GFA here, which is a tutorial with more useful info.

GFA segment (node) lines start with `S`, link lines start with `L` (spec here):

```
S   3438628 CACAATCTAAGAAATCTGAACTACCTGAAACAGGTGGAGAAGAATCAACAAACAACGGCATGTTGTTCGGCGGATTATTTAGCATTTTAG...
L   3461020 +   3461022 +   55M

```

Links connect the numeric IDs of segments to each other.

## Find high-degree nodes

We can just count the occurrences of node IDs in links, which should give us the degree.

```
egrep "\bL\b" | cut -f4 assembly_graph_with_scaffolds.gfa  | sort | uniq -c | sort | tail

# prints:

      5 3461076
      5 3461246
      5 3461286
      5 591392
      5 606626
      5 634272
      5 721568
      5 766622
      5 822586
      6 3455660

```

So we can see that there are a lot of links of segment `3455660`. We can confirm this in the GFA:

```
grep 3455660 assembly_graph_with_scaffolds.gfa
# prints:
S       3455660 CTGATAAACGGCGGGAACAGAACCAACACTACGCGTTGCTCCGATCTCAACAACTTCTAAG   KC:i:2936
L       3455660 +       56698   -       55M
L       3455660 +       3455666 +       55M
L       3455660 +       3456290 -       55M
L       3455660 -       3455652 -       55M
L       3455660 -       3456398 -       55M
P       NODE_6_length_300278_cov_129.624456_1   3448698+,3456292-,3448662+,3448670+,3451896-,3460544-,3460388+,3434104-,3434124-,56698+,3455660-,3455652-,3460374+,3460806-,3460824-,3460806-,3460824-,3460806-,3460530-,3449022-,3454210+,3459744-,3459742-    *
P       NODE_8_length_282610_cov_136.723533_1   3460440-,3434818+,3276488+,3461076-,3461070-,3461062-,3461054-,3461046-,3461038-,3460426-,3461020+,3461022+,3461020+,3460516-,3455644+,3455652+,3455660+,3456290-,3448662+,3460546-,3275518-,3387950+,3456304+,3367628+,523826-,3418048-,3456388+,3460376-,3327686+,3455382-,3460586+,655086-,3275436+,3460670+,3278936+,3460218+,3460346+,3460498+,3459664+,3460496+,3459664+,3459200+,3459424-,3458356+,3455426-,3455420-     *
P       NODE_39_length_23557_cov_173.180283_1   3456398+,3455660+,3455666+,3418048-,3273044-,3274776+,3433484-,3433478-,3433470-,3433462-,3456812-,3460134+,3455826-,3319248+,3454700+,3421874+,3453038+,3459348-,3453046-,3234810-,3460596+,3460980+,3460654+  *

```

## Extract sequences proximal to a node

Having extracted the high-degree node, we can use the output of `grep` above to find links to other sequences:

```
grep 3455660 assembly_graph_with_scaffolds.gfa | egrep "\bL"
# prints:
L   3455660 +   56698   -   55M
L   3455660 +   3455666 +   55M
L   3455660 +   3456290 -   55M
L   3455660 -   3455652 -   55M
L   3455660 -   3456398 -   55M

```

We can then go lookup `S` lines for those segments:

```
egrep "\bS\W56698" assembly_graph_with_scaffolds.gfa
# prints:
S   56698   GCATGTCGAGAAAACACCTTAGAAGTTGTTGAGATCGGAGCAACGCGTAGTGTTGGTTCTGTTCCCGCCGTT    KC:i:2067

```

## Coverage filtering

Read count (`RC` tag) is an optional field for GFA, which does not exist in the linked file.

However, we do have k-mer count (`KC`), which you could use. Presumably high k-mer count would correlate well to coverage. I don't have a suggestion of a nifty tool that filters on this, I'd use a custom script.

# Summary

Probably there exists some command line tool that does this stuff, but I don't know about it. It has to all exist inside assemblers, but it might be pretty hard to pull out.

If I had a little more time, I might write a simple tool that uses knowledge of the GFA structure to perform these kinds of data slicing operations, it seems pretty useful. It really shouldn't be hard to do this stuff!

A related approach, that I believe has a somewhat different focus, can be found here.

# Answer
> 3 votes
**Edit**: I've written a short blog post describing my use of nodeSeqs.sh to investigate the functions of sequences at the nodes of short-read assembly graphs for a set of human oral microbiome samples.

https://albertvill.com/posts/nodeseqs\_oral/

TL;DR -- transposases and recombinases are enriched, as expected.

---

Based on Maximilian's excellent answer, I have created a script to address my own question.

https://github.com/acvill/nodeSeqs

Here are the details on how it works:

> Link lines (beginning with `L`) are extracted from the GFA file, and segment occurrences are counted from both the `From` and `To` fields. To ensure that the final sequences capture the overlapping (linked) portions of segments, links are processed as segment-orientation pairs. Segment IDs with at least the specified minimum degree (`-d`) are written to `temp_XXXXX/segments.txt`, and the complementary segment lines (beginning with `S`) are pulled. *k*-mer counts (`KC` tag) are used to compute coverage for each segment using the approximation from Gonnella & Kurtz 2016; rewritten:
>
> where *K**S* is the *k*-mer count for a segment, *L**S* is the length of a segment, and *k* is the *k*-mer size used in de Bruijn graph construction (`-k`). For (meta)SPAdes, graphs are built iteratively using an increasing *k*-mer size, so *k* is equal to the largest *k*-mer size used. Sequences are extracted corresponding to segments whose coverage is at least the specified minimum coverage (`-m`). Sequences longer than the specified context (`-c`) are truncated to the context length with respect to the link orientation. Including the `-t` flag will automatically remove the temporary folder when the program is finished, which holds intermediate files and `log.txt`.

The code is indeed "shell stuff", though hopefully it is not too "kludgey". ;-)

---
Tags: assembly, metagenome, graphs, de-bruijn-graphs
---


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

1. Generate `.fastq` files from `.bams` (while preserving the `@RG` group from the original bam)
2. Split the `.fastqs`
3. Align
4. 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:

1. How do I synchronize the emitted `@RG` file (`sample_RG_file`) with the split `.fastqs`?
2. Is there a better way?
3. **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:

1. 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.
2. 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.*

# Answer
> 4 votes
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.

> 1. 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.

> 2. 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.

> 3. 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.

# Answer
> 2 votes
# 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).

<pre><code>    /*
     * 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)
```
</code></pre>

---
Tags: nextflow
---


# Question
Title: Displaying BIC and AICs from MSA nucleotide multisequence fasta file
I am initiating a project reproducing MEGA style trees using R to find the best DNA protein models under AIC/BIC. I have read the R AICs and BIC documentation but I'm struggling to implement this.

Could I implement AICs and BICs without having to specify the number of sequences in the fasta file?

My initial code that works the same as MEGA for aligning and building a phylogenetic tree where I want to view AICs and BIC (as per reference picture):

```
library(seqinr)
library(adegenet)
library(ape)
library(ggtree)
library(DECIPHER)
library(Biostrings)
library(viridis)
library(ggplot2)
library(msa)

mySequenceFile <- system.file ("sequences", "allSequences.fasta", package ="msa")
mySequence <- readDNAStringSet(mySequenceFile)

mySequence

myAlignment <- msa(mySequence, "Muscle")

myAlignment

alignment <- msaConvert(myAlignment, "seqinr::alignment")
distMatrix <- dist.alignment(alignnment, "similarity")
clustering <- hclust(distMatrix)

plot(clustering)

dendrogram <- as.dendrogram(clustering)
phylotree <- as.phylo(clustering)

plot(phylotree, type="radial")

```

Please ignore some libraries as I previously ran into some "-" deletion issues when I was writing the code

# Answer
> 1 votes
Having looked carefully at the question, I understand the 'gap in the market' but all this code will do is generate a provisional tree. I would suggest looking at protein models first because you need models to perform an AIC/BIC. So far you haven't generated used or described a protein model and without that there's nothing to compare via AIC.

To answer the direct question you should seek to avoid removing the sample size from the AIC and base the work on the equation below. If you have large samples sizes n can be removed (which is what you are wanting).

The work you should read is jModeltest which is here, https://github.com/ddarriba/jmodeltest2. This originated from Modeltest by David Pasado.

Their code on GitHub will tell you exactly how to implement AIC and BIC. I recommend ditching R in this instance. If done via subprocessing (below) its not hard.

Their work concerns nucleotides and never really hit the mainstream, but was used by some.

The key criticism is like your approach that the tree generated was based on neighbor-joining. The tree will influence the matrix so there is a danger that the ML or Bayesian phylogeny resulting from the model will be biased by neighbour-joining. Secondly, it was not computationally expensive to do this entirely by ML or Bayesian approaches because it's a 4\*4 matrix.

Protein phylogeny is very different, obviously the matrix is bigger 20\*20 the other issue is that it's a 'static matrix'. There are approaches where they "combine" matrices but I don't know the algorithm. Thus we know exactly how the protein matrices are generated, "combined" needs a clearer explanation and how this is obtained could be a good application for AIC. Thus could provide useful clarity. There is in addition.the potential to generate new matrices and assess them, but the concept of introducing meaningful dynamism into the matrix does have potential and this is absent (except for cryptic combined options). AIC would then assess whether this is helpful. It might not be hugely helpful in traditional static matrices, but would compare them against further hypotheses of mutation and innovative approaches. If used deployed correctly could be cool because it allows you to experiment and ask the question what precisely does this protein matrix mean for my data?

More generally stretching the utility of AIC type metrics is very interesting in 'bleeding edge' machine learning (admittedly its a bit futuristic).

---

I suspect (could be wrong) jModelTest is coded in Java, an alternative would be to subprocess via Python using relevant compiled code. This was how the original Modeltest was written, essentially driven through the leading maximum likelihood package of the time known as PAUP. Python does implement AIC/BIC via the statsmodel library, `scikit-learn` (machine learning) or use the equation directly.

k= number of independent variables to build model L= maximum likelihood estimate of model .. for small sample sizes (n)

R is not known for subprocessing and I am not aware of any ML/Bayes libraries for phylogeny.

In summary, this is a good project and it is worth investing the time because it has a clear, meaningful output and potential to extend into more innovative experimental approaches. I for one would use established code (i.e. written by someone else) that implemented AIC/BIC for proteins matrices and would be fairly straight forward (for you) to write. How you maximise its use in deploying it to help understand protein matrices needs creative thinking (remember its parameters verse likelihood).

Good luck.

---
Tags: r, fasta
---


# Question
Title: How to visualize electron density of a specific part of the protein residue?
Hi everyone (mostly who are interested in structural bioinformatics). I face the problem that I cannot plot the electron density of a specific part of my protein. For the whole protein it is not a problem:

```
from Bio.PDB import *
import nglview as nv
import ipywidgets

pdbl = PDBList()
path = pdbl.retrieve_pdb_file('1FAT')
parser = MMCIFParser()
struct = parser.get_structure("1fat", path)
view = nv.show_biopython(struct)
view

```

But when I tried to pick up residues from 10 till 109 (`res_list = Selection.unfold_entities(struct, "R")[10:109]`) and plot I can't do that, because for `nglview` I need to have a special Object. As I'm a newbie, I'd like to address this question to specialists.

At the moment I'm trying to visualise model (structure), not electron density. But the next step is trying to find a package to visualize electron density also. I've just started becoming familiar with some open packages.

Thank You in Advance!

# Answer
> 2 votes
NGLview runs off the JavaScript libray NGL, which is documented as a series of very helpful examples: in the gallery. The commands are not one-to-one (slightly frustratingly), but it is helpful as you can figure it out easily with a few `dir` and `help` commands. In this case whatever example of density being loaded it will reveals that it's a different component —this is the same everywhere as density files do not contain atomic coordinates. Unfortunately, there's not a direct way to generate a new `ComponentViewer` from different sources to an `NGLWidget` and you cannot merge two of them: the various `.show_*` functions return an `NGLWidget` but aren't part of it.

First, I am not sure why you are going through a convoluted route for loading a PDB code:

```
pdb_code = '1fat'

import nglview as nv

view: nv.NGLWidget = nv.show_pdbid(pdb_code)

```

Now, on loading the component a representation is added to it. These can hide the representation you then add and things will not make sense, until this detail is recalled! I fall in that pit nearly every time. But we can change it with update command:

```
# I don't like wonky sheets
view.update_cartoon(smoothSheets=True, color='salmon', component=0)

```

Now let's add the density map the circuitous way (no `.show_url`).

```
import requests, io
response: requests.Response = requests.get('https://www.ebi.ac.uk/pdbe/coordinates/files/1fat.ccp4')
response.raise_for_status()
fh = io.BytesIO(response.content)
comp: nv.component.ComponentViewer = view.add_component(fh, ext='ccp4')

```

Now we can add a surface representation, but first let's clear the default. For the arguments I had a gander at the example:

```
view.clear_representations(component=1)
view.add_representation('surface',
                        color='grey',
                        wireframe=True,
                        isolevel=1.5,
                        component=1)
view.center('100-110:A')
view

```

---
Tags: protein-structure, 3d-structure, structural-biology
---


# Question
Title: Generating contrast matrix for limma in loop
I am trying to generate contrast matrix for different conditions using **for** loop.

```
cluster_len <- 3
req_contrasts <- c()

for (m in 1:cluster_len) {
  for (n in 1:cluster_len) {
    if (n > m) {
      req_contrasts <- c(req_contrasts, (paste0(paste0("Group", n), "-", paste0("Group", m))))

    }
  }
}

```

req\_contrasts returns `"Group2-Group1" "Group3-Group1" "Group3-Group2"`

The contrast that needs to be passed is `"Group2-Group1","Group3-Group1","Group3-Group2"`

I have tried the following code to combine each string and form a single vector

```
# To separate each condition by a double quote and comma
temp_1 <- paste0(req_contrasts, collapse='","', sep="")
temp_1
[1] "Group2-Group1\",\"Group3-Group1\",\"Group3-Group2"
# To remove the escape character
temp_2 <- gsub("\\", "", temp_1, fixed = TRUE)
temp_2
[1] "Group2-Group1\",\"Group3-Group1\",\"Group3-Group2"

```

This does not work

```
contrast.matrix <- makeContrasts(temp_2, levels=design_mat)
Error in parse(text = ej) : <text>:1:14: unexpected string constant
1: Group2-Group1","
^

```

Limma throws the following error in presence of escape character. I am assuming it is because it considers that the string is yet to end.

Following input works

```
contrast.matrix_2 <- makeContrasts("Group2-Group1", "Group3-Group1", "Group3-Group2",                 levels=design_mat)

```

In my analysis the number of clusters to be analysed is 16. Hence, the contrast conditions cannot be generated manually.

# Answer
> 1 votes
I wrote this snippet a while back to generate all pairwise contrasts:

```

make_all_contrasts <- function (group, delim="_vs_"){

  suppressMessages(require(limma))

  #/ ensure that group levels are unique
  group <- sort(unique(as.character(group)))

  #/ make all combinations
  cb   <- combn(group, 2, FUN = function(x){paste0(x[1], "-", x[2])})

  #/ make contrasts
  contrasts<- limma::makeContrasts(contrasts=cb, levels=group)
  colnames(contrasts) <- gsub("-", delim, colnames(contrasts))

  return(contrasts)
}

```

Assuming three contrasts this will create a limma-like contrast matrix with all possible pairs:

```
make_all_contrasts(group=c("groupA", "groupB", "groupC"), delim="_vs_")

        Contrasts
Levels   groupA_vs_groupB groupA_vs_groupC groupB_vs_groupC
  groupA                1                1                0
  groupB               -1                0                1
  groupC                0               -1               -1

```

---
Tags: limma, loop
---


# Question
Title: Alternatives to CAP-miRSeq for microRNA processing
I am involved in a Drosophila transcriptomics project in which we have RNAseq data for small RNAs (for microRNAs) in addition to total RNA (for mRNAs).

* The total RNA reads have been analysed satisfactorily using the Tuxedo pipeline using TopHat, Cufflinks, Cuffquant and Cuffdiff. (I run this locally on a fast desktop Mac.)
* Initially we applied the Tuxedo pipeline to the short RNA reads, but abandoned this because of an incidence (1 in 10) of anomalously low or high replicates that was inconsistent with the bigwig plots from the bam stage, displayed in the UCSC browser.
* We there switched to analysing the short RNA reads using the CAP-miRSeq pipeline on Unix, which initially — though very slow — gave consistent results. However recent attempts to process later samples by this method failed as the pipeline can no longer find some dependencies it requires.

I see no point in spending further time on either of these methods.

**Question**

Has anyone *experience* of a *more modern* analysis method for microRNA RNAseq data, preferably one that would run on MacOS?

I have no trouble in running the Tuxedo pipeline and get good results with total RNA. The problem there is not my computer. Unless someone who has analysed microRNA using this package understands the problem I have encountered, there is no rationale basis for trying this again. As regards Cap-miRSeq, it does not run on the Mac, only on Unix. A colleague set it up on a virtual box. We are our support. Is there anyone working on RNAseq of microRNA?

# Answer
*Answer from @supertech, converted from comment:*

This can be frustrating I understand. In my experience it's easier to read the paper and install the modern versions of the tools compatible with your computer (paper was published in 2014). Even more convenient you could open a free account at https://galaxyproject.org/support/account-quotas/ (gives 250GB free space), and do your analysis there. Tophat aligner can be very slow. STAR is better from mapping accuracy to speed; I would recommend that you use that if possible. If you are doing analysis with large genomes your laptop computer is unlikely to have necessary RAM to run the aligners.

* Answer from @m, converted from comment:\*

Why not then simply use Bioconductor microRNA, which is under active/recent development?

Looks like STAR has been used for microRNA mapping for ENCODE:

https://www.encodeproject.org/microrna/microrna-seq/

*From the ENCODE pipeline developer (@user3274289):*

The ENCODE microRNA pipeline uses STAR for mapping and quantification and should be runnable on mac (you might need one on the beefier side though to have enough memory but i have run some workflows on a new macbook pro). We do not have drosophila STAR index available so you need to generate that using the workflow that is provided. See instructions on github.
> 1 votes

---
Tags: rna-seq, transcriptome
---


# Question
Title: Where can we download a curated dataframe of the full archive of UniProt?
I want to download a curated dataframe of the full archive of UniProt.

Any format as long as it is a table (`.csv`, `.tsv`, `SQL` etc.).

Solution 1: Download the .fasta and do the formating myself -\> very slow + information incomplete.

Solution 2: Use API like BioServices -\> still very slow but more information.

Do we have an already-curated dataframe containing as complete as possible the entry + metadata of UniProt that can be downloaded right away?

# Answer
*Answer from @matteo-ferla, converted from comment:*

There are many tables in the downloads page such tables enumerating each protein in a species proteome with basic stats, but not all the data. The uniprot records are very structured so would convert poorly to a table (even filling cells with lists/arrays and dictionaries/objects)… so a fuller table would be a mess.
> 2 votes

---
Tags: uniprot
---


# Question
Title: Normalization methods for single cell RNA sequencing that take read count into account
I have two RNA-seq datasets. One was sequenced at an average read count of 1.5 million per cell the other at 43K average reads per cell. For the first I also have the meta data from reads alligned and unaliagnable for each cell, but I may have to do some digging for the other dataset. I am looking for the presence of a few very informative genes being expressed at similar levels in the second dataset to confirmed examples of these cells in the first. The first dataset is much smaller.

My question is this: Are there normalization methods which I could use given the data I have to allow me to accurately say these cells are in the 90th or 80th percentile of a a particular gene being expressed?

# Answer
> 4 votes
The default normalisation methods of DESeq2 have a quantile-like normalisation process (based on raw read counts), and this default process should be able to handle datasets with different read depths.

If that's not good enough, DESeq2 allows you to provide alternative information (e.g. cell counts) to bypass the usual normalisation process using the `sizeFactors` function:

> The sizeFactors vector assigns to each column of the count matrix a value, the size factor, such that count values in the columns can be brought to a common scale by dividing by the corresponding size factor (as performed by counts(dds, normalized=TRUE)). See DESeq for a description of the use of size factors. If gene-specific normalization is desired for each sample, use normalizationFactors.

More details here:

https://www.rdocumentation.org/packages/DESeq2/versions/1.12.3/topics/sizeFactors

---
Tags: rna-seq, normalization
---


title: Can Primer-BLAST be run locally?

text:
# Question
Title: Can Primer-BLAST be run locally?
Primer-BLAST is a tool for designing target-specific primers, available through the Primer-BLAST web interface.

Is it possible to run the tool locally instead of accessing it via the internet?

---

Primer-BLAST obviously uses NCBI Blast which we all know about, but it is a pipeline running through Primer3. This is an extremely popular and reliable primer design software, but only does single stranded DNA/RNA. There is a prettier version, Primer3plus here.

# Answer
To strictly answer the question its not possible because it is a pipeline based on Primer3 which is not written by NCBI, therefore they would not be in a position to distribute it. Nor would Primer3 be able to distribute Blast because that is NCBI. It is simply a web-server.

NCBI and Primer3 would need to come to a joint agreement, which they must have done to present the webserver (at least it would be polite). You can of course replicate the pipeline yourself.

Brief details are here, https://www.ncbi.nlm.nih.gov/tools/primer-blast/primerinfo.html

---

It can be done via Docker. This exact pipeline hasn't been distributed via Docker but Blast is available as a Docker image.

---

**New Answer** Just to follow this up, primers4clades is a Docker image aimed at very similar objectives to Primer-Blast. I have not used it, but its worth considering if standalone is essential, e.g. for ease of building a pipeline.

Here's the blurb

> Primers4clades provides a fully automated pipeline for the design of PCR primers for cross-species amplification of novel sequences from metagenomic DNA or from uncharacterized organisms belonging to user-specified phylogenetic lineages. It implements an extended CODEHOP strategy based on both DNA and protein multiple alignments of coding genes and evaluates thermodynamic properties of the oligonucleotide pairs, as well as the phylogenetic information content of predicted amplicons, computed from the branch support values of maximum likelihood phylogenies.

I don't know the user, but Docker is a safe running environment that is isolated from the system. There's no harm in checking.

Finally both Blast and Primer3 are available via `conda`, but not Primer-Blast. In theory, I believe Primer-Blast could be distributed via `conda` under the rules of `conda`. In fact, I think, but could be wrong, that anyone could rebuild this pipeline and post it back to `conda`. @merv will know exactly whether this is 'permitable' and is very likely to correct me if I'm wrong.
> 0 votes

---
Tags: primer
---


title: Error in installing vcf2maf

text:
# Question
Title: Error in installing vcf2maf
I tried to install `vcf2maf` following this tutorial from GitHub: https://github.com/mskcc/vcf2maf.

I get a repeat error when I try to test out vcf2maf using the code:

```
perl vcf2maf.pl --input-vcf tests/test.vcf --output-maf tests/test.vep.maf

```

The error message looks like this:

```
'''STATUS: Running VEP and writing to: tests/test.vep.vcf
ERROR: Cannot find VEP script under: /Users/user/miniconda3/bin'''

```

Basically I followed every steps provided on the tutorial. Does anyone could help me find out what is going wrong?

I am working on macOS Monterey 12.2.1 with an Intel Core i5 processor.

# Answer
> 1 votes
Install it via anaconda/conda either,

```
conda install -c bioconda vcf2maf
conda install -c bioconda/label/cf201901 vcf2maf

```

I would be surprised if that installation fails. A possible issue is Perl version and system version.

---
Tags: vcf, maf
---


# Question
Title: Chosing an imputation panel for SNP-Chip data?
I have about 1,000 SNP-Chip data (samples) that I'd like to impute over (for the purpose of having more rsIDs to match against GWAS data).

However, I don't know the ancestry of each sample / the ancestry hasn't been recorded in a reliable way.

Is there a 'quick and dirty' method to decide which imputation panel to use based on the genotype data itself? e.g.

* sample 1 'looks' British in England and Scotland (GBR),
* sample 2 'looks' Colombian in Medellin, Colombia (CLM),
* sample 3 'looks' Bengali in Bangladesh (BEB), etc.

Once I know which group each sample matches most closely I can then perform imputation with the appropriate imputation panel.

# Answer
If you're using an hmm-based method such as for example impute2 or later, then the imputation method will perform the haplotype matching for you. No need to cut down the reference panel. Try to get a big reference panel such as 1000 genomes phase 3 or the haplotype reference consortium (disclaimer: I am an author on both those papers).

But do be careful about relative power: a sample with haplotypes that are better represented in the ref panel than another sample will have better-imputed genotypes than the other sample.
> 1 votes

# Answer
If I were you, I would use ADMIXTURE in supervised mode `--supervised` and use some of the 1000 genomes populations as reference populations. This is a fast and accurate way to obtain broad-scale ancestry proportion estimates for each of your individuals.
> 1 votes

---
Tags: imputation, 1000genomes, snp-chip
---


title: Imputing 23andMe 3v data to v4 and v5?

text:
# Question
Title: 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!

# Answer
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.
> 2 votes

---
Tags: vcf, imputation
---