# Question
Title: Nextflow: Can I parameterize similar processes inside a for loop?
I have set of ~5 processes which perform very similar tasks such that their differences can be parameterized. Rather than putting 5 items through some channel (which seems complicated), or defining 5 separate processes, it appears to be possible in Nextflow to declare the processes in a for loop, and output into an array of channels.

Below is a prototypical example and it seems to work as intended. (I understand this could be done differently with some channel cross product operator but I am just giving this as an example of a process in a for loop.)

Is there anything wrong with doing this? Is there a better way to do it?

Bonus question - the name of the process below is "foo". I would like to have the index reflected in the name i.e. "foo\_$i" but this doesn't seem possible by default - any way to do this??

```
some_inputs = ['prot','dna', 'rna']

some_parameters = [5,3,2,7,1]

// Set up an array of channels to output to
receivers = new ArrayList<Channel>()
for (i = 0; i < 5; i++) {
    receivers.add(Channel.create())
}

// Iterate 5 times to create 5 processes each with different parameters
for (i = 0; i < 5; i++) {

    // WOW - process can be defined inside a for loop!!
    process "foo" {

      // Extract parameters here so correct $i value is used
      p = some_parameters[i]

      input:
      val x from Channel.from(some_inputs)

      output:
      val "$x $p" into receivers[i] // Note that correct index $i seems to be used here

      """
      """
    }
}

// Confirm that each of the processes outputted independently with the correct parameter
receivers[0].view { "Received[0]: $it" }
receivers[1].view { "Received[1]: $it" }
receivers[2].view { "Received[2]: $it" }
receivers[3].view { "Received[3]: $it" }
receivers[4].view { "Received[4]: $it" }

```

Output - note that the parameter is consistent within each input channel i.e. working correctly:

```
N E X T F L O W  ~  version 20.10.0
Launching `play.nf` [fabulous_keller] - revision: c89e2b2281
executor >  local (15)
[10/b477ec] process > foo (3) [100%] 3 of 3 ✔
[16/8282ba] process > foo (2) [100%] 3 of 3 ✔
[db/f7327e] process > foo (1) [100%] 3 of 3 ✔
[fb/886f6f] process > foo (1) [100%] 3 of 3 ✔
[46/c34ed3] process > foo (1) [100%] 3 of 3 ✔
Received[1]: rna 3
Received[1]: prot 3
Received[4]: rna 1
Received[2]: rna 2
Received[0]: dna 5
Received[2]: dna 2
Received[2]: prot 2
Received[0]: prot 5
Received[3]: rna 7
Received[0]: rna 5
Received[1]: dna 3
Received[3]: dna 7
Received[4]: dna 1
Received[3]: prot 7
Received[4]: prot 1

```

# Answer
> 4 votes
Use channels, this is what nextflow's entire design is based on. I'm honestly surprised that spawning 5 processes with the same name in a loop even works.

Turn your parameters into regular channels, use the `combine` operator (https://www.nextflow.io/docs/latest/operator.html#combine) to transform the 2 parameter channels into a channel that has each combination of params/inputs. Feed this channel into the `foo` process. Afterward use `branch` (https://www.nextflow.io/docs/latest/operator.html#branch) to split your `foo` output channel into 5 different channels.

typechan = Channel.from(5,3,2,7,1)
numchan = Channel.from('prot','dna', 'rna')

cartesian = numchan.combine(typechan)

process foo {
  input:
  val numtype from cartesian
  output:
  val numtype into bar
  script:
  """
  echo $numtype
  """
}

bar.branch {
  one: it[1] == 1
  two: it[1] == 2
  three: it[1] == 3
  five: it[1] == 5
  seven : it[1] == 7
}.set {result}
result.one.view()
```
</code></pre>

# Answer
> 3 votes
You don't need to create processes in a loop when you have channels. The each qualifier allows you to repeat the execution of a process for each item in a collection. For example, the contents of 'test.nf':

```
some_inputs = Channel.of( 'prot', 'dna', 'rna' )
some_parameters = Channel.of( 5, 3, 2, 7, 1 )

process foo {

    tag { "task: ${task.index} - ${x},${y}" }

    input:
    val x from some_inputs
    each y from some_parameters

    output:
    tuple val(x), val(y) into receivers

    """
    echo "x: ${x}; y: ${y}"
    """
}

results = receivers.branch { x, y ->

    five: y == 5
        return tuple( x, y )

    three: y == 3
        return tuple( x, y )

    two: y == 2
        return tuple( x, y )

    seven: y == 7
        return tuple( x, y )

    one: y == 1
        return tuple( x, y )
}

results.five.view { "My result for param 5 is: $it" }
results.three.view { "My result for param 3 is: $it" }
results.two.view { "My result for param 2 is: $it" }
results.seven.view { "My result for param 7 is: $it" }
results.one.view { "My result for param 1 is: $it" }

```

Then running like:

```
nextflow run -ansi-log false test.nf

```

Results:

```
N E X T F L O W  ~  version 20.10.0
Launching `test.nf` [fabulous_franklin] - revision: 14c3290d2b
[6e/c12ce7] Submitted process > foo (task: 7 - dna,3)
[93/616f93] Submitted process > foo (task: 4 - prot,7)
[27/0bc3c4] Submitted process > foo (task: 3 - prot,2)
[0d/0d0b51] Submitted process > foo (task: 9 - dna,7)
[d4/f56d18] Submitted process > foo (task: 6 - dna,5)
[b0/d1535f] Submitted process > foo (task: 8 - dna,2)
[13/0a660a] Submitted process > foo (task: 5 - prot,1)
[f0/1c2f57] Submitted process > foo (task: 2 - prot,3)
[3f/3e893d] Submitted process > foo (task: 1 - prot,5)
[19/2921de] Submitted process > foo (task: 10 - dna,1)
My result for param 3 is: [dna, 3]
My result for param 7 is: [prot, 7]
My result for param 7 is: [dna, 7]
[a7/953a2b] Submitted process > foo (task: 11 - rna,3)
[1c/0fd4ef] Submitted process > foo (task: 12 - rna,5)
My result for param 1 is: [prot, 1]
[39/d21adf] Submitted process > foo (task: 14 - rna,7)
My result for param 5 is: [dna, 5]
My result for param 2 is: [prot, 2]
My result for param 2 is: [dna, 2]
[f6/8509a7] Submitted process > foo (task: 15 - rna,1)
My result for param 5 is: [prot, 5]
[6f/cee0ae] Submitted process > foo (task: 13 - rna,2)
My result for param 3 is: [rna, 3]
My result for param 1 is: [dna, 1]
My result for param 2 is: [rna, 2]
My result for param 3 is: [prot, 3]
My result for param 1 is: [rna, 1]
My result for param 5 is: [rna, 5]
My result for param 7 is: [rna, 7]

```

The answer to your bonus question is to use the `task.index` variable in your script. I used this to tag your job, but it can be used anywhere in the process declaration. You likely won't need this variable though.

Usually, using the branch operator is enough for most use cases. However,it's possible to implement your own branching function. The following uses the new DSL, DSL2, to automatically fork the 'receivers' channel:

```
nextflow.enable.dsl=2

process foo {

    tag { "task: ${task.index} - ${x},${y}" }

    input:
    val x
    each y

    output:
    tuple val(x), val(y)

    """
    echo "x: ${x}; y: ${y}"
    """
}

def branchByElement( input_ch, elements, by=0 ) {

    def targets = [:]
    for( element in elements ) {
        def key = element
        targets[ key ] = { input_ch.filter { it[by] == key } }
    }
    return targets
}

workflow {

    some_inputs = [ 'prot', 'dna', 'rna' ]
    some_parameters = [ 5, 3, 2, 7, 1 ]

    inputs_ch = Channel.fromList( some_inputs )
    parameters_ch = Channel.fromList( some_parameters )

    foo( inputs_ch, parameters_ch )

    parameter_results = branchByElement( foo.out, some_parameters, by=1 )
    input_results = branchByElement( foo.out, some_inputs, by=0 )

    parameter_results[7].call().view { "Param 7: $it" }
    input_results.dna().view { "Input dna: $it" }
}

```

Results:

```
N E X T F L O W  ~  version 20.10.0
Launching `test.nf` [distraught_meninsky] - revision: 04df251752
[bb/32325f] Submitted process > foo (task: 4 - prot,7)
[a0/19fdf5] Submitted process > foo (task: 8 - dna,2)
[e9/4c1d09] Submitted process > foo (task: 6 - dna,5)
[a3/dcb002] Submitted process > foo (task: 2 - prot,3)
[db/7e612e] Submitted process > foo (task: 5 - prot,1)
[cd/dfa84e] Submitted process > foo (task: 9 - dna,7)
[d4/6a85f2] Submitted process > foo (task: 7 - dna,3)
[f0/8e8bf4] Submitted process > foo (task: 1 - prot,5)
[89/f20a28] Submitted process > foo (task: 3 - prot,2)
Input dna: [dna, 2]
Param 7: [prot, 7]
Input dna: [dna, 5]
[53/056fe4] Submitted process > foo (task: 10 - dna,1)
[f3/a61abf] Submitted process > foo (task: 12 - rna,3)
Input dna: [dna, 7]
[f9/a568a8] Submitted process > foo (task: 13 - rna,2)
Param 7: [dna, 7]
[06/72bb6c] Submitted process > foo (task: 11 - rna,5)
Input dna: [dna, 3]
[a3/6d80d8] Submitted process > foo (task: 14 - rna,7)
[ac/acdb5a] Submitted process > foo (task: 15 - rna,1)
Param 7: [rna, 7]
Input dna: [dna, 1]

```

---
Tags: nextflow
---


# Question
Title: correlation coefficient versus DEGs analysis: what's the best approach for low expressed genes?
I have 6 folders. Each one contains 7 datasets of a specific type of cancer (RNA-Seq) and 7 datasets of normal tissue (healthy control for that type of cancer). A total of 84 datasets.

I want to investigate my gene of interest in each type of cancer.

I believe the best strategy is to merge the FPKM columns and perform correlation analyses of my gene of interest (GOI) versus all other genes in both cancer and normal tissues. As a result, I will have 14 correlation tables, one for each type of cancer and one for each control tissue. Then, I can investigate the pathways related to the genes most correlated with my GOI in each type of cancer and write my biological interpretation.

Am I in the right direction? Is this a good strategy to answer to my question?

I could also compare the control and cancer samples to investigate the significant DEGs (fold change \> 1) and the pathways activated by these DEGs, but my gene of interest doesn't have a strong pattern of gene expression, so I'm afraid it won't be in the list of statistically significant DEGs, that's why I'm preferring the correlation analysis approach.

Any other idea?

The controls were collected from the same patient and processed at the same time. I have the gene counts and can use EdgeR, but what should I do if my GOI doesn't appear among the statistically significant DEGs (fold change \> 1)? How can I investigate the role of one single gene in a RNA-Seq dataset if it is not differentially expressed? It can have a role in cancer, even if it is not an important DEG. That's why I thought correlation analyses could give me interesting information.

# Answer
First of all, RNASeq is extremely sensitive to batch effects. Are these matched controls processed at the same time as the tumors? IF they are not, lots of your differences will be caused by batch effect, and not cancer.

> I believe the best strategy is to merge the FPKM columns and perform correlation analyses of my gene of interest (GOI) versus all other genes in both cancer and normal tissues.

What you propose isn't publishable. Use the software that people have been using and validating for years. Use either DESeq2 or EdgeR, or limma if you don't have access to raw gene counts.

If the underlying counts are low, the correlation won't mean much, because your real values could be quite off. If your gene has 2 counts in one sample, and 3 counts in another, that does not mean it really correlates with a gene that has 500 reads in the one sample and 750 in the other. I don't know that you can trick your way around it; you need more reads, or more samples for more power.
> 2 votes

---
Tags: r, rna-seq, correlation
---


# Question
Title: Why am I getting empty expression data from GEO?
I am trying to analyze the scRNAseq data from this study. In their Method section they write:

> The accession number for the RNA and DNaseq data reported in this paper is GEO: GSE116237.

When I go ahead and try to pull this data in R using `GEOquery`, I am faced with a data frame with no features (i.e. no genes). The following code

```
library("GEOquery")
gse <- getGEO("GSE116237")
eset1 <- gse[[1]]
eset2 <- gse[[2]]

print(dim(eset1))
print(dim(eset2))

```

yields

```
Features  Samples
       0      184
Features  Samples
       0      681

```

I also didn't manage to download the data manually; I get an error message when trying to unzip the `.zip` files or get redirected to SRA Run Selector from which I also didn't manage to figure out how to download stuff.

Am I missing something or are these files simply corrupted? The paper is published in Cell, so I would assume this data is, to some extend, legit...

# Answer
I am sure the data is legit, you are just approaching it incorrectly. getGEO is an application for microarray data, not for digital count data such as (sc)RNA-seq, therefore what you aim to do is simply not possible by design. Unless you want to start from the raw reads why not taking the at the file named `GSE116237_scRNAseq_expressionMatrix.txt.gz` provided at the bottom of https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE116237?
> 1 votes

---
Tags: scrnaseq, gene-expression, sratoolkit, geoquery, geo
---


# Question
Title: Use conditional in workflow in Nextflow dsl2
How should I do:

```
...
if ($(params.aligner) == "bowtie") {
    align_bowtie2(get_fq_info.out.fq_info)
}
else {
    align_bwa(get_fq_info.out.fq_info)
}
...

```

# Answer
> 2 votes
You can move that into your workflow definition. There's an example in the docs of how best to use workflow parameters:

> A workflow component can access any variable and parameter defined in the outer scope:
>
> ```
> params.data = '/some/data/file'
>
> workflow my_pipeline {
>     if( params.data )
>         bar(params.data)
>     else
>         bar(foo()) }
>
> ```

Your example workflow would look like:

workflow {

    fq_info = get_fq_info.out.fq_info

    if( params.aligner == "bowtie" ) {
        align_bowtie2( fq_info )
    }
    else {
        align_bwa( fq_info )
    }

    ...
}
```
</code></pre>

# Answer
> 2 votes
Is this inside the `main.nf` file or inside the `process.nf` file? Anyway, the `param.anything` variable should already be available in the `nextflow` context when executed, thus, you could simply do:

```
if (params.aligner == "bowtie") {
    align_bowtie2(get_fq_info.out.fq_info)
} else {
    align_bwa(get_fq_info.out.fq_info)
}

```

Check out this `main.nf` of a nextflow pipeline of mine, it contains some working if/else snippets.

> Here it is another StackOverflow question on conditional nextflow pipelines that might help.

Hope it helps.

---
Tags: nextflow
---


# Question
Title: Holistic enzyme activity determination with computation
I'm currently working on a program which will determine the activity difference between an original enzyme and a variant with just 1 or 2 variants. Of course, I'm not talking about the "kinetic activity", but some indications of a probable activity difference should be given. The system has to be sensitive enough to determine the probable activity difference between unthought-of variants.

From the information that I've got from my research, there is two option. I can use one of 3 quantitive enzyme reaction analysis tools (such as EVB, FEB and LRA) or their hybrids.

EVB =\> Empirical Valence Bond

FEB =\> Free Energy Perturbation

LRA =\> Linear Response Approximation

Source: https://www.researchgate.net/post/Can-the-activity-of-proteins-be-studied-using-molecular-dynamic-simulations

Now, my questions are: - Does these application be used for such purpose? - For an accurate prediction of enzyme activity difference, are these programs good enough or require developement? - Which one of them can give the best result? - And lastly, are QM/MM methods a must for an accurate answer?

I'm a high school student so sorry if my questions sound dumb. If you can answer my inquiry's I would appreciate it:)

Thanks

# Answer
> 1 votes
Your title says `holistic`. This is a tad problematic as there's layers upon layers. Say, post-translation regulation, inhibiting metabolites, interacting protein etc. That is why when talking of an enzyme inhibitor, at the biochemical level one speaks in terms of k\_i (inhibition constant), while at the cellular level ("holistic") one talks of IC50.

This could be a biocatalysis/synthetic biology problem (make an enzyme better) or translational medicine problem (why is this mutation pathogenic).

In terms of *in silico* enzyme design, one can have all sorts of issues that need to be addressed:

* catalysis
* side reactions
* stability
* aggregation
* solubility

Also there's the expression problems, which may be due to subpar promoter, insufficient cofactors and poor codon usage.

For clinical variant assessment, a variant may:

* destabilise leaving to degradation
* remove a post-translational modification
* disrupt the active site, affecting catalysis
* lose a binding partner
* result in incorrect localisation
* result in solubility issues

So force-field scores are not the sole solution. Furthermore, a forcefield method gives one a ∆∆G (difference in relative Gibbs free energy), not a K\_M (Michaelis constant) and a k\_cat (catalytic rate), although it is somewhat possible to calculate these, or a relative variant of these.

> a variant with just 1 or 2 variants

A variant protein can have one or two *mutations* (e.g. G10R) in its sequence. A *mutant* and a *variant* are near synonyms, but variant is generally intended to be a mutant in a population.

> Does these application be used for such purpose?

For enzyme design MD simulations are not generally the first step. Rosetta pmut\_scan for example is a nice starting step as the mutational landscape of a protein is N^19, while a MD simulation takes weeks and 6+ core per single variant. This calculates only an increase in stability (relative difference in Gibbs folding energy), not reactivity. For that you can use a transition state, whose results are inaccurate, but majorly affected by the set-up, i.e. junk-in junk-out plays a big part. There's a saying 'a week in the lab saves you an hour at the library': this is backwards for sarcasm as it happens too often, namely that is background research is always a lifesaver. An QM-MD simulation on the wild type reaction definitely can help in elucidating things if it is not clear.

## Online tools

Generally, however, there are resources out there to do many stability calculations already.

**SynBio** If you are interested in picking a residue or two to make an enzyme more stable, without resorting to a thermophilic enzyme, PROSS is a very nice tool.

**Clinical** (Caveat: If your variant is a variant that affects you, a friend or a relative, you/they should speak to a genetic councillor and not do research on it yourself) There are several tools that can describe a variant at the protein level. Do remember that ∆∆G is a potential so more negative is more stable —except in SDM where the score is inverted.

---
Tags: proteins, protein-structure, proteomics, structural-biology
---


# Question
Title: Does SBOL support timing and threshold value parameters?
Like the propagation delay and threshold value parameters in electronic circuits, does SBOL format support any such parameters/glyph to represent the **propagation delay** (time to trigger the output protein once the input is provided) and **threshold value** (minimum input concentration required to trigger the model's output) of genetic models?

If currently not supported, can these parameters be incorporated in the SBOL visual representation which would give reader an idea (just by looking at the schematic diagram) about the timing behavior of a model?

# Answer
> 3 votes
*Note: I will use SBOL3 in my descriptions; the answer can be equivalently translated to SBOL2*

At present, SBOL doesn't have any explicit recommendations for how to represent these types of timing parameters. This is because there is not yet a clear scientific consensus on whether or not such parameters are a good way to represent genetic circuits.

That said, one of the advantages of using SBOL is that its extensible RDF-based model allows additional information like this to be added into a model. Two good ways of doing so are:

1. If the parameters are simple and can stand alone as measurements, then they can be added to designs as a `Measure` object with appropriate units and type. For example, if you have biological devices with a 90 minute propagation delay, then on the `Component` representing the device you can add a `Measure` with numerical value 90, units of `om:minute`. You'd likely want two types, one of `SBO:0000225 (Delay)`, and a custom type for your own notion of propagation delay, since I don't think that's included in any ontology yet.
2. If the parameters are more complex and/or inter-related, e.g., a full data-sheet, then you'd likely want to instead create one or more extension classes that are used to represent the information. A library that makes this easy is sbol\_factory, though one can do it with any of the SBOL libraries.

This becomes more complicated if you want to talk about behavior in different environmental contexts, but the basic principles are still the same.

Finally, as for SBOL visual: again, there's no explicit representation, but the model is extensible, and this information can be added as any sort of text or custom glyph as long as it doesn't conflict with ones already in the standard.

---
Tags: sbol
---


# Question
Title: how to sort and join two files based on First Column Id
I have got two files large files like this in Tab-delimited format, trying to merge in R, do we need to sort the file before merging?

```
 otu <- read.table("asv.txt", sep="\t", header=T, blank.lines.skip=F, check.names=F)
> head(otu)
                ID                   Soil.9.S42 Soil-27.S33.L001 Soil-45.S54.L001
1 d2ec9f3b77975c0f457e4b7413b217ff         34               51               84
2 3147790f0d5a78316fb9dd64f53b9473         28               32               52
3 97aecc1f35cc1f50db507ad71dd22367         33                0                7
4 bfad6370d28182cc6304844e9bec7fb6        545              237              127
5 5fa2a987221a1d9ca416148570c18086          0                0                0
6 8e60d301122d7aa359eb6b0b00f37f62        124              143              128

> taxa <- read.table("taxonomy.txt", sep="\t", header=F ,stringsAsFactors=F,quote="")
> head(taxa)
                  ID                     V2             V3                  V4
d2ec9f3b77975c0f457e4b7413b217ff Bacteria  Cyanobacteria    Oxyphotobacteria
3147790f0d5a78316fb9dd64f53b9473 Bacteria Proteobacteria Alphaproteobacteria
97aecc1f35cc1f50db507ad71dd22367 Bacteria Proteobacteria Gammaproteobacteria
bfad6370d28182cc6304844e9bec7fb6 Bacteria Proteobacteria Alphaproteobacteria
5fa2a987221a1d9ca416148570c18086 Bacteria Proteobacteria Alphaproteobacteria
8e60d301122d7aa359eb6b0b00f37f62 Bacteria Proteobacteria Deltaproteobacteria

```

Now having issue with sorting file but merge command work perfectly;

```
newotu <- otu[order(ID),]
Error in order(ID) : object 'ID' not found
> total <- merge(otu,taxa,by="ID")

```

Thanks for your help.

# Answer
You should check `left_join()` method in `dplyr` package. In this way you do not need to sort the data frames, as the function automatically take care of this. For example: `merged_df <- left_join(otu, taxa, by = "ID")`.
> 3 votes

---
Tags: r, bash, linux
---


# Question
Title: Unexpected error using "stdin" in bedtools intersect as a piped command
From the bedtools intersect man (for version 2.30.0):

```
-a BAM/BED/GFF/VCF file “A”. Each feature in A is compared to B in search of overlaps.
   Use “stdin” if passing A with a UNIX pipe.

```

The stdin part, however, does not apply to version 2.27.1. The help on this version gives no advice for piping.

```
bam mem ref fq1.fastq fq2.fastq | samtools view -b | bedtools intersect -a stdin -b blacklist.bed -v > align_black.sam
Error: unable to open file or unable to determine types for file stdin

- Please ensure that your file is TAB delimited (e.g., cat -t FILE).
- Also ensure that your file has integer chromosome coordinates in the
  expected columns (e.g., cols 2 and 3 for BED).

```

I also tried:

```
(map to bam) | bedtools intersect -a "stdin" -b blacklist.bed -v > align_black.sam
(map to bam) | bedtools intersect -a -b blacklist.bed -v > align_black.sam
(map to bam) | bedtools intersect -b blacklist.bed -v > align_black.sam

```

But these don't work either, either for the same or a different reason. Here is a site that tells me to use the above syntax: https://bedtools.readthedocs.io/en/latest/content/example-usage.html

What is the solution to this if I want to pipe everything together for my users?

# Answer
> 3 votes
With v2.27.1 (or any version I ever used) instead of `stdin` use `-`

```
(...) | bedtools intersect -a - -b blacklist.bed -v > out.bam

```

Be aware that the output of this is a BAM, not SAM file. It would be easier though to simply make a list of regions you want to keep, e.g. the complement of the backlist file with the entire genome (bedtools complement), and then use `samtools view -L` option to only keep overlapping reads. That saves you the run through `bedtools`.

# Answer
> 1 votes
Here's a way to skip conversion to BAM and pipe between tools, which performs the same set operations:

```
$ bam mem ref fq1.fastq fq2.fastq | sam2bed - | bedops -e 1 - blacklist.bed > blacklisted_reads.bed

```

With this toolkit, the hyphen (`-`) follows a Unix convention of serving as a placeholder for stdin.

# Answer
> 1 votes
While I have chosen ATpoint as the official answer, another possible syntax to use is

```
bedtools intersect -a stdin -b blacklist.bed -v < `(map to bam)` > align_black.bam

```

The main issue with this syntax is that it obfuscates the direction of flow for the pipeline, piping step 2 after step 1 is better in the official answer.

---
Tags: bedtools
---


# Question
Title: How to extract / convert gff3 CDS sequences to multifasta
I would like to extract all the CDS from a batch of genomes. I have found a perl script from BioStars but this does not seem to work for me. I would preferably like to have a script/ method which will use the locus tag as a header.

**Example**

My gff files are from prokka output

```
ATCC0000    Prodigal:2.6    CDS 243 434 .   +   0   ID=FKKLIMLP_00001;Parent=FKKLIMLP_00001_gene;inference=ab initio prediction:Prodigal:2.6;locus_tag=FKKLIMLP_00001;product=hypothetical protein
ATCC0000    prokka  gene    243 434 .   +   .   ID=FKKLIMLP_00001_gene;locus_tag=FKKLIMLP_00001
ATCC0000    Prodigal:2.6    CDS 1727    2131    .   -   0   ID=FKKLIMLP_00002;Parent=FKKLIMLP_00002_gene;inference=ab initio prediction:Prodigal:2.6;locus_tag=FKKLIMLP_00002;product=hypothetical protein

```

# Answer
I like `bedtools getfasta`. My typical option set is `bedtools getfasta -fi <reference> -bed <gff_file> -name -s`. Be aware of the `-s` to make sure you are pulling the correct strand. I like bedtools because it is a versatile tool overall for handling bed, gff and vcf file manipulations.

```
# bedtools getfasta

Tool:    bedtools getfasta (aka fastaFromBed)
Version: v2.26.0
Summary: Extract DNA sequences from a fasta file based on feature coordinates.

Usage:   bedtools getfasta [OPTIONS] -fi <fasta> -bed <bed/gff/vcf>

Options:
    -fi Input FASTA file
    -bed    BED/GFF/VCF file of ranges to extract from -fi
    -name   Use the name field for the FASTA header
    -split  given BED12 fmt., extract and concatenate the sequencesfrom the BED "blocks" (e.g., exons)
    -tab    Write output in TAB delimited format.
        - Default is FASTA format.

    -s  Force strandedness. If the feature occupies the antisense,
        strand, the sequence will be reverse complemented.
        - By default, strand information is ignored.

    -fullHeader Use full fasta header.
        - By default, only the word before the first space or tab is used.

```
> 6 votes

# Answer
The gffread utility in Cufflinks package might be interesting for you. To generate a multi-fasta file with nucleotide sequences from your GFF file, then you can try:

```
gffread -w output_transcripts.fasta -g reference_genome.fa input_transcripts.gff

```
> 4 votes

# Answer
Using python and this GFF parser that mimics Biopython's SeqIO parsers:

```
from BCBio import GFF

# Read the gff
for seq in GFF.parse('my_file.gff'):
    # only focus on the CDSs
    for feat in filter(lambda x: x.type == 'CDS',
                       seq.features):
        # extract the locus tag
        locus_tag = feat.qualifiers.get('locus_tag',
                                        ['unspecified'])[0]
        # extract the sequence
        dna_seq = seq = str(feat.extract(seq).seq)
        # simply print the sequence in fasta format
        print('>%s\n%s' % (locus_tag, dna_seq))

```
> 3 votes

# Answer
The `xtractore` program from the AEGeAn Toolkit was designed for this type of use case. Just set `--type=CDS`.

```
$ xtractore -h

xtractore: extract sequences corresponding to annotated features from the
           given sequence file

Usage: xtractore [options] features.gff3 sequences.fasta
  Options:
    -d|--debug            print debugging output
    -h|--help             print this help message and exit
    -i|--idfile: FILE     file containing a list of feature IDs (1 per line
                          with no spaces); if provided, only features with
                          IDs in this file will be extracted
    -o|--outfile: FILE    file to which output sequences will be written;
                          default is terminal (stdout)
    -t|--type: STRING     feature type to extract; can be used multiple
                          times to extract features of multiple types
    -v|--version          print version number and exit
    -V|--verbose          print verbose warning and error messages
    -w|--width: INT       width of each line of sequence in the Fasta
                          output; default is 80; set to 0 for no
                          formatting

```

It looks like you're processing prokaryotic genomes, but this program works also on eukaryotic genomes where CDSs often have introns interspersed.
> 2 votes

# Answer
You can use `agat_sp_extract_sequences.pl` from AGAT. See here for all option: https://agat.readthedocs.io/en/latest/tools/agat\_sp\_extract\_sequences.html.

Below example of use:
> 0 votes

---
Tags: fasta, software-recommendation, format-conversion, gff3
---


title: Optitype for Singularity

text:
# Question
Title: Optitype for Singularity
I am trying to use `optitype` on an HPC without docker support or sudo privileges. A search brought me https://github.com/researchapps/optitype and https://github.com/nf-core/hlatyping. I tried bootstrapping the singularity image as follows:

```
singularity build optitype.simg docker://fred2/optitype

```

However, this does not seem to run as expected. Could someone with experience bootstrapping docker files onto singularity suggest how I could change this? How do the commands change when doing a bootstrap? I would appreciate your help!

# Answer
> 1 votes
Bit of an old post sorry, but I think you might be getting a little confused here. The nf-core pipelines come with their own software already, so you don't need to install or bootstrap optitype yourself. If it's used in the pipeline, it comes in the pipeline container(s) (or conda environment or whatever packaging system you use). So you should never have to touch `docker://fred2/optitype`.

When you run the pipeline with the `-profile singularity` argument, it will automatically pull the image from Docker Hub and convert to Singularity for you *(some of the newer pipelines pull pre-build Singularity images directly)*. You can also use the `nf-core download hlatyping --singularity` command from the helper tools to pull the image ahead of time.

I hope this helps! Come and join the nf-core Slack organisation if you need more help. We have a channel for each pipeline, including `#hlatyping`.

---
Tags: docker
---


title: Find common SNPs in multiple VCF files

text:
# Question
Title: Find common SNPs in multiple VCF files
I have 3 VCF files.

```
A1.vcf
A2.vcf
A3.vcf

```

I want to get the common SNPs that are present in all these three files. And output must be in vcf format.

output:

```
common_A.vcf

```

# Answer
A potential solution is to compress and index the files using bgzip and tabix, then use bcftools isec to find SNPs common to all three files, e.g.

for files in A1.vcf A2.vcf A3.vcf
  do
  bgzip "$files"
  tabix "$files"
done

bcftools isec -n=3 -c all -o all_common_SNPs.vcf A1.vcf.gz A2.vcf.gz A3.vcf.gz
```

> 3 votes

---
Tags: vcf, snp, variants
---


title: Percentage of each cluster in Seurat

text:
# Question
Title: Percentage of each cluster in Seurat
I am using Seurat to analyze my single cell data. I have 2 conditions, treated and untreated. I am trying to create a stacked bar graph in order to show the differences in cell types for each condition but need to collect the percentages of each cluster for the specific cell types. How to I put together a sheet that contains the percentages of each cluster and the make a stacked bar graph of comparison of each condition.

# Answer
> 2 votes
```
pt <- table(Idents(myobj), myobj$orig.ident)
pt <- as.data.frame(pt)
pt$Var1 <- as.character(pt$Var1)

ggplot(pt, aes(x = Var2, y = Freq, fill = Var1)) +
  theme_bw(base_size = 15) +
  geom_col(position = "fill", width = 0.5) +
  xlab("Sample") +
  ylab("Proportion") +
  scale_fill_manual(values = brewer.pal(12, "Paired")) +
  theme(legend.title = element_blank())

```

This shall work. `myobj` is your seurat object. This is my old code so if I remember correct, `orig.ident` is your treated and untreated conditions. `Idents` should give you cell clusters. If this does not work, please let me know

# Answer
> 1 votes
You can access these data from a Seurat object as columns of `object@meta.data` and tabulate them.

There are functions such as “propellor” which will also implement this:

https://github.com/Oshlack/speckle

Note that a Seurat vignette also explores this topic:

https://satijalab.org/seurat/archive/v3.0/interaction\_vignette.html

You need to be careful when doing this to ensure different proportions are biological as single cell data is prone to technical errors.

---
Tags: r, scrnaseq, seurat
---


title: Clustering with Seurat

text:
# Question
Title: Clustering with Seurat
My data consists of two groups control and stim. I'm concerned that my clustering in the stim group is dominated by the unregulated genes. I can remove these genes and recluster but is there a way to then project the removed genes onto the new umap space so see if they are indeed all clustered together. They are no longer in the seurat object so I can't think of how to do this. I hope this makes sense.

# Answer
> 1 votes
Yes, this can be done, particularly if your concern is with the clustering / UMAP rather than the transformation process; I've done this recently with some Rhapsody single cell data that we generated. The main trick is adding *'return.only.var.genes = FALSE'*, as well as calculating a variable set that is a bit larger than the desired set. Depending on the number of cells, this will balloon the size of the resulting Seurat object pretty significantly.

Note that UMAP is a technique to place *cells*, rather than genes. If you have the information on where the cells should go, then plotting gene expression on top of that via FeaturePlot should work without any reprojection needed.

Here's what I did (with experiment-specific variables excluded):

```
sct.data <- SCTransform(sct.data, vars.to.regress="percent.mt",
                        variable.features.n = 3000 + length(genesToExclude),
                        return.only.var.genes = FALSE);

sct.data <- RunPCA(sct.data, verbose = TRUE,
                   features=setdiff(VariableFeatures(sct.data), genesToExclude));
sct.data <- RunUMAP(sct.data);
sct.data <- FindNeighbors(sct.data);
sct.data <- FindClusters(sct.data);

```

All transformed genes should be available for querying in the results.

---
Tags: r, seurat, clustering
---


# Question
Title: Can an alternate allele be more common than a reference allele?
I'm trying to analyze data from a genotype-phenotype association study. The genotype data is stored in csv files that record for each variant site and each subject, the alleles that the subject has at that site. A quick online search reveals that 0 means the reference allele, 1 means the most common alternate allele, 2 means the 2nd most common alternate allele, etc.

But I noticed that at some variant sites, most (\>90%) subjects have 1/1, i.e. two copies of the alternate allele. I asked around and some suggested I just swap the roles of reference and alternate alleles whenever I find that the alternate allele frequency \> 0.5. My questions are:

1. Can there be a good reason to classify a less common allele as the reference allele? I thought that the reference allele was defined as the most common one in the population.
2. Is swapping the roles of reference and alternate alleles whenever I find that the alternate allele frequency \> 0.5 a standard practice?

Thanks

# Answer
> 1 votes
1. You are mixing up a bit the terminology here. One thing is a reference allele and another is a major allele. The first refers to the genotype present at a given position in the "reference" genome (e.g. if you work with human data this can be the reference version hg19 or hg38, etc). The second refers to the most common allele in a population. This definition can lead to problems since what is the most common allele can differ among different populations. In this respect, it is perfectly fine to have a reference allele with a lower frequency than another alternate allele.
2. No for the reasons above. If you mean whether swapping major to minor is standard practice whenever you find out the major allele is actually minor as per the definition (less frequent), then I guess you can but this will lead to inconsistencies whenever you try to merge your data with other data sets, or other people merge their data with your data. Since this for sure happens and issues about strandness are well known, one has to make sure that you make explicit how you treat your data once you release it, and also one has to make sure that major/minor are defined the same way when merging data, and do the fix whenever needed.

# Answer
> 1 votes
1. Yes, depending on how the reference allele was determined. If, for example, you want to refer to hg38 as the reference genome, then it's important to use the hg38 allele as reference.
2. Without stating the reference allele convention? Not that I'm aware of. However many people refer to minor/major alleles rather than non-reference alleles, in which case it is assumed that the major allele is the most common allele.

As long as you state what your reference genome and/or population is, then it's fine to choose your own convention about reference alleles. But please be consistent and explicit about what you've done. One of my biggest struggles in GWAS has been working out the correct strand orientation for stated variants (important when merging datasets), and it doesn't help for complementary variants (i.e. A/T, C/G) when the ref/alt alleles are arbitrarily chosen.

---
Tags: variant-calling, gwas, allele-frequency
---


# Question
Title: DESeqDataSetFromTximport all(lengths > 0) is not TRUE
I am suddenly running into an error when running the DESeqDataSetFromTximport

```
txi.rsem <- tximport(files, type = "rsem", txIn = FALSE, txOut = FALSE)
dds <- DESeqDataSetFromTximport(txi = txi.rsem,
                                colData = SampleFile,
                                design = ~ Compound)`

using counts and average transcript lengths from tximport
Error in DESeqDataSetFromTximport(txi = txi.rsem, colData = SampleFile,  :
  all(lengths > 0) is not TRUE

```

any suggestions what causes this error?

# Answer
I found this same error as of today using this code:

```
pad.gene.data <-
  tximport(files = file.list.for.gene.analysis,
           type = "rsem",
           txIn = FALSE,
           txOut = FALSE)

dds.gene.level <-
  DESeqDataSetFromTximport(txi = pad.gene.data,
                           colData = pad.characteristic.data,
                           design = ~ condition)

```

I had done this same differential gene expression analysis using transcript counts aggregated to the gene level, on a mac about 5 weeks ago, but as of this past week switched to a windows and re-installed R 3.6.3 and whatever the current version of tximport and DESeq2 is.

I was able to replicate the results I had gotten on the mac 5 weeks ago removing all rows where all gene lengths are greater than 0:

```
pad.gene.data$abundance <-
  pad.gene.data$abundance[apply(pad.gene.data$length,
                                1,
                                function(row) all(row !=0 )),]

pad.gene.data$counts <-
  pad.gene.data$counts[apply(pad.gene.data$length,
                             1,
                             function(row) all(row !=0 )),]

pad.gene.data$length <-
  pad.gene.data$length[apply(pad.gene.data$length,
                             1,
                             function(row) all(row !=0 )),]

```

So first removing any gene where there is a length count of 0 from the abundance and counts matrix, then removing the genes from the length matrix so that they are all the same size. It was able to run that way.

I am definitely skeptical of why I would need to do that as it doesn't seem that it was advertised, if this is a new version, that we need to preprocess the data beforehand. I hope someone who might know more of the code can comment on this.
> 1 votes

# Answer
I had this issue as well. It looks like all zero-length genes also have zero expression across all samples. The code below makes sure only to exclude genes that have both zero length and zero expression.

txi = tximport(file_paths, type = "rsem", txIn = TRUE, txOut = TRUE)

zero_length_and_unexpressed = (apply(txi$abundance, 1, max) == 0) &
                              (apply(txi$length, 1, min) == 0)

txi$length = txi$length[!zero_length_and_unexpressed,]
txi$abundance = txi$abundance[!zero_length_and_unexpressed,]
txi$counts = txi$counts[!zero_length_and_unexpressed,]

dds = DESeqDataSetFromTximport(txi, sampleTable, ~ 1)
```

> 0 votes

---
Tags: deseq2, tximport
---


title: Sequence allignment with suffix array?

text:
# Question
Title: Sequence allignment with suffix array?
I'm writing the code from scratch and not using libraries for the actual indexing/search since it's a project for school. Any blackbox advice is fairly useless.

I'm planning on aligning some short reads with the *D. melanogaster* genome but, since I was thrown into this project with very little info/ experience I'm stuck on a few things.

Firstly, for preprocessing the fasta file:

* What should I do with Ns (ambiguous bases), I've seen some people removing them or replacing them by gaps but this does not seem like the same information? Should they be left in?
* As far as I understand, the lowercase letters stand for potential repeats. However, If I just ignore this, it is possible my sequences might not align to the repeat because it is not actually stated as a repeat. E.g: "ATATATAT" will not align to 'at'.

As far as I understand the difference between using suffix array and suffix tree: suffix tree will search in linear time (but is memory intensive), while suffix array will be space efficient (but searches through ALL sequences each time). I'm leaning towards suffix array because I'm working on my laptop and memory is kind of limited. I'm also thinking that, since you sort the array lexicographically, you really only have to check the sequences if they start with the same first base as the read you're searching for?

Another question on the arrays: if my genome has 100M base pairs, the array will have 100M entries, some of them being close to 100M bases long? Should I be using the seed and extend method?

On the other hand, since the fruit fly genome is one of the smaller genomes, maybe it would be feasible to use a tree?

Is there any way to validate the alignment quickly at the end? How much error should I be allowing?

This were a lot of question, I've been trying to read up on things but I can mostly only find applications using biopython etc.

Any help would be hugely appreciated!

# Answer
> 2 votes
> * What should I do with Ns (ambiguous bases), I've seen some people removing them or replacing them by gaps but this does not seem like the same information? Should they be left in?
> * As far as I understand, the lowercase letters stand for potential repeats. However, If I just ignore this, it is possible my sequences might not align to the repeat because it is not actually stated as a repeat. E.g: "ATATATAT" will not align to 'at'.

Good questions. I think this depends on the scope of your class project. Frankly, I think it would be best to just start with unambiguous sequences (only ATGC) and not worry about potential repeats. It's also a good software engineering exercise to solve the easy case first in a way that's flexible and allows you to expand to the additional complexities above. For example, the original Bowtie software considered references with ambiguous characters invalid. You can see for yourself here how Bowtie2 handles it.

> As far as I understand the difference between using suffix array and suffix tree: suffix tree will search in linear time (but is memory intensive), while suffix array will be space efficient (but searches through ALL sequences each time). I'm leaning towards suffix array because I'm working on my laptop and memory is kind of limited. I'm also thinking that, since you sort the array lexicographically, you really only have to check the sequences if they start with the same first base as the read you're searching for?

Yes, the most naive implementation of a search with a suffix array will be $O(m\log n)$ time, where $m$ is the read length and $n$ is the genome length. However, with enhanced suffix arrays you can also obtain the same $O(n)$ search complexity of suffix trees. There is also the use of suffix arrays to construct Burrows-Wheeler transforms and even further compression with the FM index. There are open source implementations of all of these data structures that may be interesting to you.

> Another question on the arrays: if my genome has 100M base pairs, the array will have 100M entries, some of them being close to 100M bases long? Should I be using the seed and extend method?

Suffix arrays store the indices of the start of each suffix in order to save space, so each entry is a single integer.

> On the other hand, since the fruit fly genome is one of the smaller genomes, maybe it would be feasible to use a tree?

Its worth trying to compute this yourself (: You can determine the number of bytes required per base pair in your input and see how that scales with your available memory. For example, if your suffix tree implementation required $20n$ bytes, you'd end up w/ roughly 2Gb needed for a fruit fly genome.

> Is there any way to validate the alignment quickly at the end? How much error should I be allowing?

Again, this depends on the scope of the project. Allowing for errors in mapping reads is an active area of research. See the supplemental of the Bowtie2 paper for a dynamic programming algorithm that handles a certain number of mismatches.

> This were a lot of question, I've been trying to read up on things but I can mostly only find applications using biopython etc.

All good questions, and I think many answers depend on the scope of the project. Reimplementing a read mapper is no small feat, especially when you start considering all of the different complexities such as sequencing error and ambiguous bases. I would strongly consider starting with looking at just the pure string problem of read mapping: Only consider ATCG in both reference and reads, and don't even consider reads who do not have exact matches to the reference. From there, time permitting, you can start adding on to your solution. As far as reading goes, there's a deep rabbit whole of these topics, from tools like bwa-mem, bowtie2, mini-map, to the more theory-oriented papers discussing optimized versions of the data structures. I can attach some here if you would like. Hope these answers helped!

# Answer
> 2 votes
Suffix arrays vs Suffix tree is a problem quite often discussed in bioinformatics classes, but it's not that much used in practice. Most of people use already optimized mappers bowtie2 or bwa-mem.

I would advise you against your intention of mapping the reads using your own implementation of anything. You don't need to reinvent the wheel, just use an existing solution.

By the way, the Ns in the reference are either hard-masked nucleotides (that would be repeats of some sorts) or (more often) simply unknown nucleotides. So there is no way to map your reads there anyway (and the mentioned mappers will handle the Ns for you).

---
Tags: sequence-alignment, fasta, homework
---


# Question
Title: Median string problem & multiple sequence alignment
I read about the median string problem as an introduction to the multiple sequence alignment, however none of the MSA algorithms used seems to be using the idea of finding the median string.

To my understanding from the consensus (median) string, it would be possible to retrieve the actual alignment by simply doing the pairwise alignment of the consensus string with each of the strings.

Why is this approach not used? Are there any algorithms that use it?

And related to this, what kind of tasks in bioinformatics/biology are median/consensus string algorithms used for?

# Answer
> Why is this approach not used? Are there any algorithms that use it?

As shown in this paper, the MEDIAN STRING problem is NP-Complete on alphabets of carnality greater than 3. While the problem isn't used directly in approaches due to the complexity, the idea of trying to create a consensus sequence between pairs of input and using that as building blocks is common. These are considered "progressive alignment methods".
> 1 votes

---
Tags: sequence-alignment, multiple-sequence-alignment
---


title: Create GFF3 feature exporter

text:
# Question
Title: Create GFF3 feature exporter
Problem Statement

> The GFF3 format is a commonly-used one in bioinformatics for representing sequence annotation. You can find the specification here:
>
> http://www.sequenceontology.org/gff3.shtml
>
> The genome and annotation for Saccharomyces cerevisiae S288C is on the class server here:
>
> ```
> /home/jorvis1/Saccharomyces_cerevisiae_S288C.annotation.gff
>
> ```
>
> Note that this same file has both the annotation feature table and the FASTA sequence for the molecules referenced. (See the '##FASTA' directive in the specification.) Within the feature table another column of note is the 9 th , where we can store any key=value pairs relevant to that row's feature such as ID, Ontology\_term or Note. Your task is to write a GFF3 feature exporter. A user should be able to run your script like this:
>
> ```
> $ export_gff3_feature.py --source_gff=/path/to/some.gff3 --type=gene --attribute=ID --value=YAR003W
>
> ```
>
> There are 4 arguments here that correspond to values in the GFF3 columns. In this case, your script should read the path to a GFF3 file, find any gene (column 3) which has an ID=YAR003W (column 9). When it finds this, it should use the coordinates for that feature (columns 4, 5 and 7) and the FASTA sequence at the end of the document to return its FASTA sequence. Your script should work regardless of the parameter values passed, warning the user if no features were found that matched their query. (It should also check and warn if more than one feature matches the query.)
>
> The output should just be printed on STDOUT (no writing to a file is necessary.) It should have a header which matches their query, like this:
>
> > gene:ID:YAR003W .... sequence here ...
>
> As an extra challenge, you can format the sequence portion of the FASTA output as 60-characters per line, which follows the standard.

So far I have this below. I'm not even sure if this is correct to start with:

```
#!/usr/bin/env python3

for line in open("Saccharomyces_cerevisiae_S288C.annotation.gff"):
   line1 = line.rstrip()
   if line.startswith("#"):
       continue
   column = line1.split("\t")
   if len(column) != 9:
       continue
   id = column[8]
   type = column[2]

```

Can anyone help me with this? I'm totally lost. Thank you!

# Answer
Your answer looks like a good start. You'll probably want to use sys.argv, argparse, or another argument parsing tool to read in options from the command line, then compare them to the fields (e.g. the "type" argument gets compared with the type column).

Given that attributes can be arranged in an arbitrary order, you'll probably need to read all of them from each line and store in a dict in order to determine whether the target is present. A regular expression search is another way to do that.

There are a few libraries available (e.g. from biopython) for reading and subsetting FASTA records.
> 0 votes

---
Tags: gff3, homework
---


# Question
Title: RNAseq DE comparison with samples of different read length
I have RNA samples but with different read lengths (Eg, HiSeq 2x125 and NovaSeq 2x150bp data). I would like to do DE analysis on these samples. What do you recommend? Do you recommend to trim the reads to 125bp for the NovaSeq to compare with the other batches from HiSeq? Or any other suggestions.

# Answer
> 0 votes
You don't need to do anything, this is taken care of by the linear model - if you don't have any confounding effects.

If you do have confounding effects, e.g. all control samples on HiSeq and all treated on NovaSeq there's not a lot you can do anyway, other than be honest when you present your results.

Taken care of by the model means adding an extra 2 level factor for sequencer by sample, you need to do this anyway to take care of the batch effect between sequencer runs.

# Answer
> 0 votes
Yes, I would simply trim the longer one to the shorter read length, then map the data as usual. I doubt though that it makes a big difference on the global scale but the processing step is simple and it makes sure you have no mappability bias. If you do not have any other confounding effects then you can do a DE analysis without any extra covariates, I think it is always desirable to keep things as simple as possible and trimming is the simplest solution.

---
Tags: rna-seq
---


# Question
Title: What are the units for the force output by nstfout, the velocity output by nstvout, and the coordinates output by nstxout in gromacs?
What are the units of force, velocity, and coordinates output by gromacs trajectory? In `.mdp`, they are specified by `nstfout`,`nstvout`,`nstxout`.

# Answer
> 3 votes
`nstxout`, `nstvout` and `nstfout` are the number of **steps** between coordinates, speed and forces written in the output trajectory (Output control documentation).

As for the units present in your trajectory files (say a `.trr`): coordinates are in **nm**, speeds in **nm.μs⁻¹** (so in **mm.m⁻¹**) and forces in **kJ.mol⁻¹.nm⁻¹**. These should be consistent through all GROMACS files (Definitions and Units documentation).

Note that `.xtc` files multiply coordinates by 1000 to compress them, so they are stored in **pm** (see GROMACS xtc documentation). Given MDAnalysis documentation on the subject, I believe `.xtc` files can *only* store coordinates, no speeds or forces.

I would be wary of softwares used to analyze your trajectory results (MDAnalysis, VMD, etc), as they may automatically do the conversion of these units, and for example show you coordinates in Å.

---
Tags: simulated-data, gromacs
---


# Question
Title: How to input the opposite strand of DNA?
I am trying to write a program that will count the repeats of a given DNA strand:

I am using a Counter library, something like this code shown below

```
from collections import Counter

        s="XXXXXAGCCTGCCAAGCAAACTTCACTGGAGXXXXXTGTGCGTAGCATGCTAGTAACTGCAXXXXXTCTGAATCTTTCAGCTGCTTGTXXXXXTGGGCCTCTCACAAGGCAGAGTGTCTTCATGGGACTTTGATATTTATTTTTGTACAACCTAAGAGGAACAAATCCTTTGACACTGACAAATTGGCTTCCATATTTTATACCTTAATCATCTCCATGTTGAATTCATTGATCAACAGTTTAAGAAAAAAAGATGTAAAAATGCTTTTAGAAAGAGAGGCAAAGTTATGCACAATAACTTCTCATGAAGTCACAGTTTGTTAAAAGTTGCCTTAGTTCACAATAAATAATTATGTATGCTCTATAATTTCAGTGA"

def chunk_string(s, n):
    return [s[i:i+n] for i in range(len(s)-n+1)]

counter = Counter(chunk_string(s, 5))

```

Repeat supposed to be "declared" as an identical sequence on the same strand (like `XXXXX` here that occurs 4 times). However, I can also have this repeat (for ex. `XXXXX` in this particular test case) on the OPPOSITE strand.

How am I suppose to intake the opposite strand in case of using `.fasta` files? As I know for today - fasta is presented only in a single strand way. Or am I completely mixed the terms?

# Answer
Yes, you're going to want to take the reverse complement of your string in order to analyze the other strand, that is the standard.

Also, it is worth pointing out that you may want to use a generator expression instead of a list expression for your `chunk_string` function:

```
def chunk_string(s, n):
    return (s[i:i+n] for i in range(len(s)-n+1))

```

This will result in a generator that yields kmers as the `Counter` constructor iterates through them, which means that you don't need to first construct and store the entire list of kmers first and then pass to `Counter` but instead can feed them directly to it. Not essential for small strings, of course, but a good practice and useful for larger strings. See this SO post for a better explanation.

You can also use Biopython's `reverse_complement` function in order to easily obtain the reverse complement.
> 1 votes

# Answer
What's going to be simplest is to generate the reverse complement of your short search sequence, and search for that too.
> 0 votes

---
Tags: python, fasta
---


# Question
Title: Using Entrez.efetch() to retrive .fasta file from any NCBI database?
I tried to use these python codes to retrieve a .fasta file from gene database, but the output were the parameters of the report

```
from Bio import Entrez
Entrez.email = "send_helps@earth.com"
handle = Entrez.efetch(db = "gene",id = "682",rettype="fasta", retmode="xml")
print(Entrez.read(handle))

```

i also read the documentation on this site about efetch() and seems like this function can only retrieve certain kinds of data for each database. So i can't take fasta file from gene database !?

I don't want to use the NCBI web UI. What's wrong with my code? Is there a way i can guarantee to retreive a .fasta file through python from NCBI ? Am i misundertanding something severly ?

Sorry for bad English. thank you.

# Answer
> 2 votes
In order to download the FASTA sequence of transcript variants using Entrez eutilities, you will have to use `Entrez.elink` to get a list of transcript accessions for the gene of interest and then use `Entrez.efetch` to fetch the sequence(s). Just to give you an idea, you can use Entrez Direct for this as follows:

```
elink -db gene -id 682 -target nuccore -name gene_nuccore_refseqrna \
  | efetch -format fasta

```

You can find more Entrez Link descriptions here. Note, there is a `gene_fasta` option for `efetch` but it still returns only the transcript sequence (i.e., no intron sequence).

If you need the full gene sequence (including the intronic regions), you can use the commands mentioned here.

Alternatively, you can use NCBI Datasets for this. You can use the web interface, command line tool or the python API. This page with the how-to articles has the information you seek.

---
Tags: fasta, python, phylogenetics, entrez, efetch
---


# Question
Title: How do I add a colour annotation bar to the heatmap generated by DoHeatmap function of Seurat v2?
I am using **Seurat v2** for professional reasons (I am aware of the availablity of **Seurat v3**). I am clustering and analysing single cell RNA seq data. How do I add a coloured annotation bar to the heatmap generated by the **DoHeatmap** function from Seurat v2? I want to be able to demarcate my cluster numbers on the heatmap over a coloured annotation bar. I want to have the an image like this:

(taken from Satija V3 tutorial).

My original code:

```
library(Seurat)
DoHeatmap(object = object, genes.use = top20_pos$gene)

```

produces the following image which is not what I want:

I am aware that **DoHeatmap** function returns a **ggplot2 object** but after searching the internet for a long time, I have not been able to find which **geom** layer I can add to achive my objective.

On the other hand, I have tried to generate a separate coloured annotation bar from another package called **ComplexHeatmap** with the intention of cropping it to fit the cluster demarcations on my heatmap generated with **DoHeatmap** but this has proven to be very time-consuming and inaccurrate:

```
library(ComplexHeatmap)

df = data.frame(type = c(rep("a", 18), rep("b", 3), rep("c", 2), rep("d", 1), rep("e", 1)))
ha = HeatmapAnnotation(df = df, col = list(type = c("a" =  "blue", "b" = "red", "c" = "green", "d" = "cyan", "e" = "magenta")))
draw(ha, 1:25)

```

which produces:

I will be very grateful for any hint on how to achieve my goal automatically and accurately.

Thank you in advance for your help.

# Answer
> 1 votes
I am not sure if it's a v3 or also v2 thing, but have you tried by setting `group.by` and `label`?

From `?Seurat::DoHeatmap`:

```
group.by   A vector of variables to group cells by; pass 'ident' to group by cell identity classes
label      Label the cell identies above the color bar

```

If you prefer external tools, you can use `pheatmap` and setting `annotation_col` e `annotation_colors`.

First define a data.frame (df) containing the data you want in the bar (clusters in your case)

```
Description Condition
liver       Tumour
lung        Normal

```

Then you have to assign the colours

```
anno <- as.data.frame(df)
anno_colori <- list(
   Condition = c(Normal = '#FFC000', Tumour = '#0070C0'),
   Description = c(liver = '#D55E00', lung = '#66CC99')
)

pheatmap(mat,
     annotation_col = anno,
     annotation_colors = anno_colori,
)

```

The code above will add 2 bars on the top of your heatmap, one for `Condition` and another for `Description`.

# Answer
> 0 votes
There was a discussion about this in Seurat developer's github and someone actually made a package to easily add several color bars using some DoHeamap-like function. Check it here.

---
Tags: r, scrnaseq, seurat, ggplot2, heatmap
---


# Question
Title: Why models of stochastic gene expression predict that intrinsic noise should increase as the amount of transcript decrease
I am reading Elowitz et al. (2002). When talking about intrinsic noise, i.e. noise due to microscopic events that govern which reactions occur and in what order, it is stated that:

> Models of stochastic gene expression predict that intrinsic noise should increase as the amount of transcript decreases.

Noise, defined in the article as the standard deviation divided by the mean, can be divided in intrinsic (due to microscopic events during transcription) and extrinsic (due to localization, state and concentration of other cellular components, e.g. polymerases).

I think the intrinsic noise should increase at low level of transcript amount because the average (mean) amount will be lower but deviation decrease is smaller than mean decrease. I don't know, however, if this is the actual reason (and why deviation decrease should be smaller than mean decrease) or if deviation decreases with the mean and there's another more likely explanation.

1. Elowitz, M. B., Levine, A. J., Siggia, E. D. & Swain, P. S. Stochastic gene expression in a single cell. Science 297, 1183–1186 (2002).

# Answer
> 2 votes
Reading Kaufmann and van Oudenaarden (2007), it seems to validate the first alternative (using results from the Central Limit Theorem):

> Although biochemical fluctuations influence all stages of gene expression, those involving molecules in extremely low abundance are expected from a statistical standpoint to be larger in magnitude and therefore to contribute disproportionately to the overall variation between cells.

1. Kaufmann, B. B. & van Oudenaarden, A. Stochastic gene expression: from single molecules to the proteome. Curr. Opin. Genet. Dev. 17, 107–112 (2007).

# Answer
> 0 votes
The inverse relationship between variance and mean is a consequence of the stochastic nature of gene expression, which can be modeled as a discrete birth-death process. One of the simplest models assumes that the gene is expressed constitutively. Then, using the chemical master equation, mRNA production can be shown to evolve as a discrete Poisson process. The steady-state distribution is Poisson. For the Poisson distribution, coefficient of variation ($\sigma^2/\mu^2$) is equal to $1/\mu$. Further, if protein is translated in bursts from mRNA (translation burst), protein noise can be written as

\begin{equation}\label{eq1}\tag{1} \eta^2\_P = \frac{\sigma\_P^2}{\mu\_P^2} = \frac{B}{\mu\_P}, \end{equation}

where $\eta\_P^2, \sigma\_P^2,\mu\_P,$ and B are the coefficient of variation, variance, mean and burst size for protein. From \eqref{eq1}, the inverse relationship is evident.

Gene expression models can be made more detailed by assuming different states for the promoter. The promoter switches between these states, and each of these states can have different rates of transcription. For these models, $\eta\_P^2$ can again be shown to follow an inverse relationship with $\mu\_P^2$.

Reference:

Paulsson, Johan. "Models of stochastic gene expression." Physics of life reviews 2.2 (2005): 157-175.

---
Tags: noise, gene-expression
---


# Question
Title: How to increase the size of the shapes that denotes effect sizes of different datasets and metaanalysis?
I am doing metaanalysis with forestplot from ggforestplot package. I have 7 different datasets from 7 different countries. After metaanalysis, I have 7 Effect sizes and a Meta-analysis effect size, i.e., total of 8 effect sizes. Now, I have generated a forest plot using this code:

```
# Forestplot
my_plot <- forestplot(
  df = final,
  estimate = estimate,
  logodds = FALSE,
  colour = country=="metaanalysis",
  shape = country,
  title = "Associations to disease accross populations",
  xlab = "Effect Size",
  #xlim = c(-1.5, 1.5),
  xtickbreaks = c(-1.0, -0.5, -0.45, -0.4, 0, 0.4, 0.45,0.5, 1.0)
)+
  # You may also want to add a manual shape scale to mark meta-analysis with a
  # diamond shape
  ggplot2::scale_shape_manual(
    values = c(5L, 0L, 1L, 2L, 8L, 19L, 4L, 3L),

  )

```

But, the shapes denoting the countries and the metanalysis seem very small. How can I make them larger for better visualization? Can anyone please help me?

**Representative data for replication:**

```
final <- data.frame(
  stringsAsFactors = FALSE,
              name = c("RIBOSYN2-PWY__flavin_biosynthesis_I_(bacteria_and_plants)",
                       "PWY-6385__peptidoglycan_biosynthesis_III_(mycobacteria)",
                       "PWY-6151__S-adenosyl-L-methionine_cycle_I",
                       "PWY-6387__UDP-N-acetylmuramoyl-pentapeptide_biosynthesis_I_(meso-diaminopimelate_containing)",
                       "PWY-6700__queuosine_biosynthesis",
                       "PEPTIDOGLYCANSYN-PWY__peptidoglycan_biosynthesis_I_(meso-diaminopimelate_containing)",
                       "PWY-6386__UDP-N-acetylmuramoyl-pentapeptide_biosynthesis_II_(lysine-containing)",
                       "PWY-5667__CDP-diacylglycerol_biosynthesis_I",
                       "PWY0-1319__CDP-diacylglycerol_biosynthesis_II",
                       "PWY-6122__5-aminoimidazole_ribonucleotide_biosynthesis_II",
                       "PWY-6277__superpathway_of_5-aminoimidazole_ribonucleotide_biosynthesis",
                       "PWY-6163__chorismate_biosynthesis_from_3-dehydroquinate",
                       "TRPSYN-PWY__L-tryptophan_biosynthesis",
                       "PANTOSYN-PWY__pantothenate_and_coenzyme_A_biosynthesis_I",
                       "GLYCOCAT-PWY__glycogen_degradation_I_(bacterial)",
                       "P161-PWY__acetylene_degradation",
                       "FUC-RHAMCAT-PWY__superpathway_of_fucose_and_rhamnose_degradation",
                       "P108-PWY__pyruvate_fermentation_to_propanoate_I",
                       "PWY-7328__superpathway_of_UDP-glucose-derived_O-antigen_building_blocks_biosynthesis",
                       "TCA__TCA_cycle_I_(prokaryotic)",
                       "PWY-6901__superpathway_of_glucose_and_xylose_degradation","PWY66-400__glycolysis_VI_(metazoan)",
                       "FUCCAT-PWY__fucose_degradation",
                       "PRPP-PWY__superpathway_of_histidine,_purine,_and_pyrimidine_biosynthesis",
                       "PWY-7234__inosine-5'-phosphate_biosynthesis_III",
                       "TRPSYN-PWY__L-tryptophan_biosynthesis",
                       "PWY-6122__5-aminoimidazole_ribonucleotide_biosynthesis_II",
                       "PWY-6277__superpathway_of_5-aminoimidazole_ribonucleotide_biosynthesis",
                       "RIBOSYN2-PWY__flavin_biosynthesis_I_(bacteria_and_plants)",
                       "PWY-6163__chorismate_biosynthesis_from_3-dehydroquinate",
                       "PWY-5667__CDP-diacylglycerol_biosynthesis_I",
                       "PWY0-1319__CDP-diacylglycerol_biosynthesis_II",
                       "PANTOSYN-PWY__pantothenate_and_coenzyme_A_biosynthesis_I",
                       "RIBOSYN2-PWY__flavin_biosynthesis_I_(bacteria_and_plants)",
                       "PWY-6386__UDP-N-acetylmuramoyl-pentapeptide_biosynthesis_II_(lysine-containing)",
                       "RIBOSYN2-PWY__flavin_biosynthesis_I_(bacteria_and_plants)",
                       "PANTOSYN-PWY__pantothenate_and_coenzyme_A_biosynthesis_I",
                       "PWY-6385__peptidoglycan_biosynthesis_III_(mycobacteria)",
                       "PWY-6387__UDP-N-acetylmuramoyl-pentapeptide_biosynthesis_I_(meso-diaminopimelate_containing)","PWY-6700__queuosine_biosynthesis",
                       "PWY-6151__S-adenosyl-L-methionine_cycle_I",
                       "PWY0-1319__CDP-diacylglycerol_biosynthesis_II",
                       "PWY-5667__CDP-diacylglycerol_biosynthesis_I",
                       "PEPTIDOGLYCANSYN-PWY__peptidoglycan_biosynthesis_I_(meso-diaminopimelate_containing)",
                       "PWY-6387__UDP-N-acetylmuramoyl-pentapeptide_biosynthesis_I_(meso-diaminopimelate_containing)",
                       "PWY-6385__peptidoglycan_biosynthesis_III_(mycobacteria)",
                       "PANTOSYN-PWY__pantothenate_and_coenzyme_A_biosynthesis_I",
                       "PEPTIDOGLYCANSYN-PWY__peptidoglycan_biosynthesis_I_(meso-diaminopimelate_containing)","PWY-5667__CDP-diacylglycerol_biosynthesis_I",
                       "PWY0-1319__CDP-diacylglycerol_biosynthesis_II"),
          estimate = c(0.492903365186722,
                       0.490116714405755,0.470681520902659,0.465726944992238,
                       0.461786439528495,0.455469435286562,0.449174326308655,
                       0.433755737833203,0.433670791448977,0.412009484469875,
                       0.412009484469875,0.409849212159568,0.407818381473059,
                       0.406187518931948,-0.405573831841341,-0.407806403622317,
                       -0.418351108232386,-0.419107645158175,-0.420781760133428,
                       -0.441930843508614,-0.444723954813865,-0.451224756237618,
                       -0.504642907784859,-0.515843696383127,-0.523823484265933,
                       1.01203740887784,0.860911808112254,0.860911808112254,
                       0.782824395343453,0.779871038531055,0.763412247455786,
                       0.763412247455786,0.744965434413444,0.735086809796288,
                       0.733132943355509,0.731855318166865,0.730523378560539,
                       0.728356885227315,0.694487466356475,0.694335496605963,
                       0.68469489116557,0.680753784497535,0.680732580641012,
                       0.668909750729078,0.665107737526883,0.662352062210176,
                       0.643767077562848,0.640204416487105,0.639460661647805,
                       0.637911093993277),
                se = c(0.10858902830706,
                       0.101375498952004,0.101193730594653,0.101213265518358,
                       0.101165178275411,0.101149604926044,0.101114528367103,
                       0.101177617463096,0.101176196339295,0.10102421433221,
                       0.10102421433221,0.100910703206783,0.101046445558438,
                       0.12265509469116,0.10486856029419,0.103267024648877,0.100867566762125,
                       0.11467836793049,0.10790947678297,0.100997454994136,
                       0.101684644439116,0.103275688295495,0.101358559584742,
                       0.101485340026678,0.131050876787633,0.12672527384678,
                       0.0911987727585976,0.0911987727585976,0.0473284627336248,
                       0.120947206065826,0.120594420271743,0.120594420271743,
                       0.120207965256539,0.12000489747134,0.119965054342127,
                       0.0890575546896603,0.0468935746714197,
                       0.0890043416222531,0.119198789457309,0.0466102635462743,
                       0.119011209638725,0.0883056120686039,0.0883053113523379,
                       0.0881390968551948,0.0880862640207957,0.118593197976584,
                       0.0877953030586674,0.118192523540133,0.0462081105569776,
                       0.0461972348759043),
           country = c("metaanalysis","metaanalysis",
                       "metaanalysis","metaanalysis","metaanalysis",
                       "metaanalysis","metaanalysis","metaanalysis","metaanalysis",
                       "metaanalysis","metaanalysis","metaanalysis",
                       "metaanalysis","metaanalysis","metaanalysis","metaanalysis",
                       "metaanalysis","metaanalysis","metaanalysis",
                       "metaanalysis","metaanalysis","metaanalysis","metaanalysis",
                       "metaanalysis","metaanalysis","Swedish","Kazakh","Kazakh",
                       "Japanese","Swedish","Swedish","Swedish","Swedish",
                       "Swedish","Swedish","Kazakh","Japanese","Kazakh",
                       "Swedish","Japanese","Swedish","Kazakh","Kazakh",
                       "Kazakh","Kazakh","Swedish","Kazakh","Swedish","Japanese",
                       "Japanese")
)

```

**Output:**

Thanks

# Answer
A quick look at the code for forestplot shows that geom\_effect is what needs to be modified. If you add a value for fatten changes the size of the shapes in the the figure (but not in the legend), or a value for size will change the shape, line and legend sizes (add to not geom\_effect not to the aes layer)

This may be easier said than done as forestplot also seems to call an unexposed function (function quo\_is\_null). A couple of simple hacks are possible; copy and make your own forestplot in the default environment and remove the call to the function (it seems to just check if pvalue is set), or change the environment of the "new" forestplot to ggforestplot

```
envorionment(forestplot) <- envornment(ggforestplot::forestplot)

```

There may be other less awkward methods, but this will work. Fatten changed to size 6 which is massive.

```
forestplot <-
function (df, name = name, estimate = estimate, se = se, pvalue = NULL,
          colour = NULL, shape = NULL, logodds = FALSE, psignif = 0.05,
          ci = 0.95, ...)
{
  stopifnot(is.data.frame(df))
  stopifnot(is.logical(logodds))
  name <- enquo(name)
  estimate <- enquo(estimate)
  se <- enquo(se)
  pvalue <- enquo(pvalue)
  colour <- enquo(colour)
  shape <- enquo(shape)
  args <- list(...)
  const <- stats::qnorm(1 - (1 - ci)/2)
  df <- df %>% dplyr::mutate(`:=`(!!name, factor(!!name, levels = !!name %>%
                                                   unique() %>% rev(), ordered = TRUE)), .xmin = !!estimate -
                               const * !!se, .xmax = !!estimate + const * !!se, .filled = TRUE,
                             .label = sprintf("%.2f", !!estimate))
  if (logodds) {
    df <- df %>% mutate(.xmin = exp(.data$.xmin), .xmax = exp(.data$.xmax),
                        `:=`(!!estimate, exp(!!estimate)))
  }
  if (!quo_is_null(pvalue)) {
    df <- df %>% dplyr::mutate(.filled = !!pvalue < !!psignif)
  }
  g <- ggplot2::ggplot(df, aes(x = !!estimate, y = !!name))
  if (logodds) {
    if ("xtickbreaks" %in% names(args)) {
      g <- g + scale_x_continuous(trans = "log10", breaks = args$xtickbreaks)
    }
    else {
      g <- g + scale_x_continuous(trans = "log10", breaks = scales::log_breaks(n = 7))
    }
  }
  g <- g + theme_forest() + scale_colour_ng_d() + scale_fill_ng_d() +
    geom_stripes() + geom_vline(xintercept = ifelse(test = logodds,
                                                    yes = 1, no = 0), linetype = "solid", size = 0.4, colour = "black")
  g <- g + geom_effect(ggplot2::aes(xmin = .data$.xmin, xmax = .data$.xmax,
                                    colour = !!colour, shape = !!shape, filled = .data$.filled),
                       position = ggstance::position_dodgev(height = 0.5),
                       fatten=6) +
    ggplot2::scale_shape_manual(values = c(21L, 22L, 23L,
                                           24L, 25L)) + guides(colour = guide_legend(reverse = TRUE),
                                                               shape = guide_legend(reverse = TRUE))
  if ("title" %in% names(args)) {
    g <- g + labs(title = args$title)
  }
  if ("subtitle" %in% names(args)) {
    g <- g + labs(subtitle = args$subtitle)
  }
  if ("caption" %in% names(args)) {
    g <- g + labs(caption = args$caption)
  }
  if ("xlab" %in% names(args)) {
    g <- g + labs(x = args$xlab)
  }
  if (!"ylab" %in% names(args)) {
    args$ylab <- ""
  }
  g <- g + labs(y = args$ylab)
  if ("xlim" %in% names(args)) {
    g <- g + coord_cartesian(xlim = args$xlim)
  }
  if ("ylim" %in% names(args)) {
    g <- g + ylim(args$ylim)
  }
  if ("xtickbreaks" %in% names(args) & !logodds) {
    g <- g + scale_x_continuous(breaks = args$xtickbreaks)
  }
  g
}
environment(forestplot) <- environment(ggforestplot::forestplot)

```
> 1 votes

---
Tags: r
---


# Question
Title: What are common ways to calculate gene length for TPM calculations?
I am aware of this similar question. But the accepted answer there answers how to calculate TPM given a mapping from gene name to gene length.

My question is, given an annotation file of genes (for example from Ensembl here), how do you calculate gene length for the purpose of calculating TPM on a gene count matrix coming from RNA seq data. That is, I understand the actualy calculation of TPM given gene lengths, but I do not understand what the gene lengths are.

From what I understand, to calculate gene lengths from annotation file for RNA seq data, you only consider the exons of each gene, but even so, many exons overlap. What is the standard method of dealing with this?

For example, I have read that using the length of the union of all exons of the gene is acceptable, but from what I understand TPM is meant to represent the "expected number of reads coming from gene X when randomly selecting 1e6 transcripts". So taking the union of all exons will be an estimate of this but it might under-estimate this.

EDIT: From what I can tell this question is not a duplicate of RNASeq: Normalization, stabilization, gene length and rlog That question does discuss how to perform normalization that takes gene length into account, but it does not ask "what is gene length? How does one determine the length of each gene?" Possible answers could be like what I already mentioned: length of the union of all exons in the gene, or another possible answer is one from the linked question "length of the longest transcript isoform" but I feel that it is not really discussed in that answer/question.

# Answer
I think the most common way is to use something like Kallisto which figures out what proportion of reads are likely assigned to each transcript variant, and calculated TPM based on that.
> 1 votes

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


# Question
Title: How to know if ID correspond to nucleotide sequence or protein with a script
I have a text file that contains a list of IDs (314 sequences):

```
AVP78031.1
AVP78042.1
ATO98108.1
ATO98120.1
ATO98132.1
...

```

My goal is to make a script (maybe using Python or Perl) to check in the list if all the IDs are nucleotide or protein sequences. For example:

```
AVP78031.1 -> protein (this is a nucleotide sequence, I change nucleotide for protein to show an example).
AVP78042.1 -> nucleotide
ATO98108.1 -> nucleotide
ATO98120.1 -> nucleotide
ATO98132.1 -> nucleotide

```

Any idea to do a script?

Thank everybody!

# Answer
> 2 votes
If these are all GenBank or RefSeq accessions, you can use Entrez Direct for this as shown below:

```
$ cat accs.txt
ATO98108.1
ATO98120.1
ATO98132.1
AVP78031.1
AVP78042.1
$ cat accs.txt | epost -db nuccore | efetch -format acc
## no output because none of them are nucleotide accessions
$ cat accs.txt | epost -db protein -format acc | efetch -format acc
AVP78042.1
AVP78031.1
ATO98132.1
ATO98120.1
ATO98108.1

```

**NOTE**: This will work only if the accessions are currently live because `epost` does not find any suppressed accessions. For example:

```
$ cat accs.txt
NM_002826.3
NM_002826.4
NM_002826.5
$ cat accs.txt | epost -db nuccore -format acc | efetch -format acc
NM_002826.5

```

Here, all three accessions are valid nucleotide accessions but only the last one, NM\_002826.5, is alive.

An alternate way is to use the accession prefixes defined here and come up with an appropriate regular expression query.

---
Tags: python, perl, sequence
---


title: mRNA cap identification by sequencing?

text:
# Question
Title: mRNA cap identification by sequencing?
If you sequence a positive-sense RNA virus genome using NGS via Nanopore, Illumina or even the old Sanger method, is it possible to tell if the virus has 5' (five prime) cap at the 5' end of the RNA genome?

# Answer
It's possible to do 5' cap enrichment sequencing (e.g. see here), and presumably you could look at differences between an enriched vs non-enriched sequencing run to determine whether or not there's a substantial number of sequences with the cap.

Unfortunately, because the cap is right at the end of the sequence, it zips through a nanopore very quickly, so may be quite difficult to distinguish from any other base at the end. There are also no official models that have been trained to specifically detect the cap, so you'd be initially restricted to looking for signal deviations. If you could find a way to extend the sequence for 20 or so based beyond the cap site, it might be possible to detect it at a single-molecule level.

Despite all that, this is something that UCSC is working on; see their repository here:

https://github.com/mitenjain/dRNA\_capping\_analysis
> 1 votes

# Answer
In short the answer is no, however the 'car-jumping method' described here it is possible, but it is about the enzymic MMLV reverse transcription in response to the capping structure. I assume that an alternative reverse transcriptase doesn't do 'car jumping' therefore any sequencing method will pick this up, but the difference between MMLV verse the heat-tolerant style reverse transcriptases. I would assume the secondary structure associated with viral capping prevents a low temperature reverse transcriptase, i.e. MMLV from correctly reproducing the vRNA into cDNA.

It is important to note that a cap, isn't just a virus thing, its mRNA eukaryotes per se. The cap is a 7-methylguanosine (m7G) linked to the 5′ end of the transcript by a 5′-5′ triphosphate bridge \[m7G(5′)ppp(5′)N- or m7G(5′)ppp(5′)Nm in which Nm is a 2′-O methylated nucleoside\].
> 0 votes

---
Tags: sequence-alignment, covid-19, sars-cov-2
---


# Question
Title: Normalize RNA seq data from multiple runs for expression analysis
I have RNA samples sequenced with TruSeq Stranded Total RNA kit protocol in Illumina HiSeq (2x125bp) and NovaSeq platforms (2x150bp) - almost 100 samples altogether. I have to use the samples data for expression analysis. Question is how do we normalize samples from different runs and of read length? What is best way (with available packages) we can normalize the count (to overcome the bias) and do a DE analysis? ANy suggestion would be helpful.

# Answer
> 1 votes
There are a number of methods. If you're doing DE, you have ComBatSeq, SVAseq, RUVseq, BUSseq. You could also try Z-score normalization ($(x - \overline{x}) / \sigma$) or even quantile transformation. For the latter two, make sure you work on each batch individually, not the whole dataset at once. See more on that here.

To visually compare if transformation works, plotting PCA/UMAP/t-SNE on raw and transformed data can perhaps be of some insight.

# Answer
> 0 votes
I don't think RNASEq is that sensitive to running different things on different instruments. Differing read lengths will change things. You could always trim the reads down so they all match. Or include read length as a factor in the design when looking for DE genes.

# Answer
> 0 votes
Normalization and batch-effect correction are two different things.

You always have to do normalization of RNAseq data in order to be able to compare the runs and do DE analysis, for instance. Plenty of methods and you have to search for them as normalization methods indeed.

What you are asking and concerned about is "batch-effects" in your data as a function of having used two different Illumina machines. This is a valid concern, but before you try to correct for it or mitigate it, you first have to assess whether it is actually there. Maybe is not at all. To figure out if you have batch-effects you have to first do data exploration, this is why the use of PCA/UMAP/t-SNE has been already suggested. In this plots you will be trying to see whether your data cluster together by machine type instead of by the biology it represents. Read about how to assess the presence of batch effects.

If you do observe it or anyway want to correct for its presence, then indeed you can try one of the methods also already mentioned: ComBatSeq, SVAseq, RUVseq, BUSseq.

Rather than the two different machines, the batch-effect might actually be due to different people doing the lib preparation, experimental settings if experiments were done at different moments, etc.

# Answer
> 0 votes
I agree with swbarnes that you don't need to explicitly model for different read length, but you need to normalize for different production of sequencing runs. Among the several available packages, DEseq2 is one of the most widely used. The normalization method computes the geometric mean of expression of each gene across samples; then, expression of the gene in each sample is divided by the geometric mean. Once this has been don for all genes, in each sample it computes the median of values, and this is the "size factor". Ok, it looks convoluted, but it is indeed easy, and it is implemented in DESeq2 and in other software packages. You can find methodological details here. Finally, DESeq2 also offer the possibility to account for additional experimental details (e.g. read length) and this should reduce the risk of spurious results.

# Answer
> 0 votes
1. If you have different read length then trim all samples to the shortest length across batches, that will avoid mappability bias. Including it into the DE design is possible but as you can completely eliminate this bias in silico by trimming this is what I'd do.
2. Sequencing on different machines is usually not much of a confounding factor, but the fact that two machines were involved suggest different batches of experiments. I hope that batch is not confounded by group so each batch should contain samples of all groups to avoid a general batch effect that cannot be corrected for. Library prep should (in a well-designed experiment) be fully identical.
3. Use any of the established methods, be it RLE or TMM (they usually perform ery similar), and then perform PCA based on the top variable log2-transformed counts to explore whether there are batch effects that are worth correcting, so whether you see clustering or differences based on batch rather than experimental design in the early PCs. If so (given it is not confounded by experimental groups) you can include it into the DE design as a covariate.

As said by others, normalization and batch correction are independent processes. Also as I said in on top, you can eliminate read length bias (even though I doubt it is a major factor here) in silico which is what I'd recommend so you do not have to bother with it downstream.

---
Tags: rna-seq, differential-expression, gene-expression, normalization, featurecounts
---


# Question
Title: Why Do I Get Different Results Each Time From AutoDock Vina?
I just tried to use Vina for calculating binding affinity. I used the same protein, the same ligand, the same options, but there are different results. Why do I get different results? Is it the limit Point of docking tools? Then why we use it?

# Answer
> 2 votes
Vina uses a non-deterministic docking algorithm. See the manual.

> The docking algorithm is non-deterministic. Even though with this receptor-ligand pair, the minimum of the scoring function corresponds to the correct conformation, the docking algorithm sometimes fails to find it. Try several times and see for yourself. Note that the probability of failing to find the mininum may be different with a different system.

The manual states that you can set the random seed in order to obtain reproducibility.

> I changed something, and now the docking results are different. Why? Firstly, had you not changed anything, some results could have been different anyway, due to the non-deterministic nature of the search algorithm. Exact reproducibility can be assured by supplying the same random seed to both calculations, but only if all other inputs and parameters are the same as well. Even minor changes to the input can have an effect similar to a new random seed. What does make sense discussing are the statistical properties of the calculations: e.g. "with the new protonation state, Vina is much less likely to find the correct docked conformation".

# Answer
> 1 votes
The answer by @Throckmorton is 100% correct in that it's not deterministic. I just wanted to address a few things —so this is a long comment.

All too often the box is set as the whole protein, this is bound to make it trip up on local minima. Sure, with native ligands, binding pockets close to the active site increases the dissociation constant of the complex, but docking is not the way to calculate these. Generally this is done by tethered MD, where the ligand is pulled out.

Docking returns generally either a **binding potential** (∆∆G\_bind) as a kcal/mol (Autodock) or kJ/mol, where negative is good, or a rank. This is not the typically binding affinity measurements, e.g. an IC50 or dissociation constant (K\_D) —The Eyring equation converted in a relative form (ratio makes constants disappear) allows you do some brute conversions to a k-on for the site, but this is highly not recommended.

Docking is not perfect and, worse, gets abused a lot. However I'd flip the OP's statement around and say the limitation by far is wet lab side, where doing an ITC properly is hard and crystallising a protein with a ligand is even worse...

---
Tags: docking
---


# Question
Title: How do I plot clusters based on DNA sequence alignment?
I have a group of about 1050 sequences that have been aligned using blastn. I am now trying to create a cluster plot of the sequences in order to identify similar sequences, something like this:

However, I am having trouble finding a tool to do this. My professor recommended using ClustalW, and so I have spent plenty of time searching for a tool that can take a sequence alignment and plot them such as in the linked image. However, I have not been successful in finding what I need and am beginning to think I am missing something or something is going over my head. If anyone could assist with finding the tool I need I would very much appreciate it. Thanks so much!

# Answer
> 2 votes
I think one of the reasons you struggle is that clustering DNA sequences is not a clearly defined task. In general intention of clustering is to reconstruct, or approximate the relatidness of the DNA seqeunces.

If all these sequences are homologous you can do a multiple sequence alignment (using MAFFT or Clustal) and build a phylogenetic tree (using something like IQ tree). If your prof suggested ClustalW, this is probably what they meant. Evolutionary speaking, that's the optimal method for getting relations of sequences, but they must be "alignable". Also with this many sequences, it will be computationally quite heavy.

For within species relatedness, phylogenetic reconstruction is sometimes impossible, because recombination breaks down the linkage of individual alleles. However, if all of the sequences are long enough, you can do a principal component analysis on SNPs, for example using something like plink. The plot you showed actually looks a lot like PCA (it even has the explained variability on axes).

Finally, regardless of their relationship, you can always do a dendrogram based on a similarity matrix, but for that, you need to define similarity. That will be usually something kmer-based. I am not following this field very closely lately, but there are some recent papers about kmer-based clusterings, like this one.

When looking for the paper I also run into this one, which proposes genome clustering using "Inter-nucleotide Covariance". No idea how useful/sound it is. Just found the abstract quite interesting.

Long story short, it really depends on what you are after and what are your sequences! If you can build a reliable tree, it is the most powerful way to determine the relationships among sequences.

---
Tags: visualization, clustering, multiple-sequence-alignment
---


# Question
Title: Modeling number of reads mapped to a gene
I am looking for a probability distribution of a number of reads mapped to a particular gene in metagenomic sequencing (NGS, shotgun, likely illumina).

Naively one could model it via a binomial (or Poisson) distribution as: $$ P\_n = {N\choose n}\pi^n(1-\pi)^{N-n}, $$ where

* $n$ is the number of reads mapped to the gene
* $N$ is the sequencing depth (total number of reads)
* $\pi\propto l$ is the probability of read mapping to the particular gene, proportioal to the gene length (possibly adjusted for the read length)

This is however complicated by such intricacies of sequencing as:

* finite length of the fragments, which can be longer or shorter than the gene length
* paired-end sequencing, which means one has to be careful when counting evens where only one or both reads map to a gene
* something else?

Is there an established model for this case? Could you recommend reading on the subject?

# Answer
> 2 votes
This is quite a similar question to this one, but I notice that you're specifically asking about a *probability distribution* rather than differential expression. That sounds like gene dispersion and/or shot noise; a good reference for me on that has been the documentation for DESeq2:

http://www.bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html

The gene dispersion issue is discussed in more detail in their 2014 paper:

https://doi.org/10.1186/s13059-014-0550-8

> We model read counts K ij as following a negative binomial distribution (sometimes also called a gamma-Poisson distribution) with mean μ ij and dispersion α i . The mean is taken as a quantity q ij, proportional to the concentration of cDNA fragments from the gene in the sample, scaled by a normalization factor s ij , i.e., μ ij =s ij q ij . For many applications, the same constant s j can be used for all genes in a sample, which then accounts for differences in sequencing depth between samples. To estimate these size factors, the DESeq2 package offers the median-of-ratios method already used in DESeq. However, it can be advantageous to calculate gene-specific normalization factors s ij to account for further sources of technical biases such as differing dependence on GC content, gene length or the like, using published methods, and these can be supplied instead.

---
Tags: ngs, read-mapping, reads, illumina, genome-sequencing
---