# Question
Title: Using SBOL for lineages of recombinant DNA molecules and strains
I wonder about using The Synthetic Biology Open Language (SBOL) for a project involving the development of a web application to track the lineage of recombinant DNA molecules in a laboratory collection. I am wondering whether SBOL can support this use case, and I could really use some help from community members.

### Use case

The ultimate goal is to develop a web application to keep an audit trail of recombinant DNA in a laboratory collection, where sequences are traceable through cloning steps up to their entry point into the collection.

For instance, we can imagine we want to document the generation of a new plasmid by combining fragments of two existing plasmids by restriction and ligation. The web application would provide an interface where you can select the input sequences, select the restriction sites and the assembly, and would generate the child sequence. This information would then be stored in a relational database, and would also be exportable to a text file. Essentially, something like the history of a sequence in SnapGene, but using a relational database instead of the file system, and having an open standard to export the information.

### Implementation (roughly)

This requires the following steps:

1. **Generating a schema to document genetic engineering steps:** The most abstract version of a step would just have inputs and output, but there would have to be a specific implementation for each operation. For the restriction and ligation described above, there would be three steps

* Digestion plasmid A: which enzyme was used, and where it cut.
    * Digestion of plasmid B: which enzyme was used, and where it cut.
    * Ligation of pieces: Order of the pieces of the assembly
```
Plasmid A -[restriction]-> Restriction Fragment A
                                                |-[Ligation]->Final product
Plasmid B -[restriction]-> Restriction Fragment B

```
2. **Building an API that reads the inputs and description of a step and returns the output:** For this, I had thought of using the library pydna, since it seems to cover the initial use cases I could think of.
3. **Building a visual interface to generate this documentation:** A mockup with minimal functionality can be found here
4. **Developing a web application managing a database based on the schema mentioned above**

### Will SBOL ever cover this use case?

As far as I understand it, SBOL currently provides the most abstract layer of the schema described in the previous section, supported by the `<prov:>` tags, so each entity could be related to its ancestors through an `Activity` that represents the ligation, the restriction. etc. However, from looking at the documentation there seems to be no schema for any concrete implementation of `Activity`. For example, to describe a restriction operation. Restriction is relatively simple, but for PCR there would be the need to also define a schema for the alignment of primers, with potential missmatches etc. Provenances should also support importing entries from public resources (such as genome databases). These concrete `Activity` schemas should be readable by an API that would take the inputs and return the output. For cloning steps, this would be performed with pydna or something similar. For sequences imported from public resources, the provenance would have the url of a web request, and the API would request the sequence, and parse it.

From the documentation of `prov:Activity`: *to capture how an entity was derived, it is expected that any additional information needed will be attached as annotations*. From this, I wanted to ask whether specific schemas representing the cloning steps will ever be part of the SBOL, or should be provided as external files.

### How to approach implementation?

My original plan (before realising about `<prov>` tags in SBOL) was to do the schema in json, since it's the easiest way to go for a web application. I suppose that the json schema would be easily convertible into xml to eventually be used in SBOL.

I still think it would be best to start with json, since I don't think for the nature of this task I can leverage any of the features of SBOL beyond the `<prov>` tag. The sequences themselves will be stored in the original import format, and the only new aspect that I will introduce would be the schemas for cloning steps. Once I have a minimal prototype of the schemas and I have tested the application with users, I could come back to these forums and discuss the schema and potential changes that can be made to it, if it is to be included in SBOL.

If you have some thoughts in any of these points, I would really appreciate any feedback.

# Answer
> 1 votes
SBOL 3 is well suited for an application like this, though you may wish to consider some of the specifics (e.g., whether to use a conventional database or an RDF triple store), and has been used for similar (though not exactly the same) workflows in other projects.

To some of your specific questions:

* In addition to tracking the provenance relationships, SubComponent relations can be used to track the inclusion of parts within your final build. The sourceLocation field can be used to express the difference between the original part and a post-assembly trimmed version when there are changes in the reaction.
* The assembly plans and reactions can be represented not just with activities, but using SBOL Interaction objects.
* SBOL Implementation objects can be used for tracking specific aliquots as well as the designs that are intended to be included in them.
* If you still need to enhance Activity records, you can either use the base extension capabilities in pySBOL3 or (better) use SBOL Factory to create extension classes via an ontology.
* With regards to serialization format, since SBOL 3 is pure RDF, you can use any RDF serialization you want, including JSON-LD.

There's a lot to unpack here, so I would suggest that if you want to dig deeper into any of these specific responses that you ask new questions focused on individual aspects.

Here is an example in prototype work on protocols:

https://github.com/SD2E/paml/blob/main/paml/paml.ttl

---
Tags: sbol
---


title: Demovir produces an empty output

text:
# Question
Title: Demovir produces an empty output
I am new to bioinformatics and am working on viral metagenomic (virome) datasets and so would like to use Demovir\* for doing taxonomic annotations on my viral contigs.

First, I installed all prerequisites in my macOS but I found that Usearch 32-bit version could not run on macOS Big Sur, so I installed all and re-ran `demovir.sh my.viral.contig.fasta threads` in ubuntu VM. I could therefore successfully obtained trembl\_ublast.viral.u.contigID.txt as an input for R, but I would always get an empty output after performing the demovir.R command of demovir.sh in ubuntu VM.

Thus, I took these back to my macOS, used conda to create a virtual environment for R-3.5.1 and installed R using `conda install -c r r=3.5.1`, then ran `Rscript demovir.R trembl_ublast.viral.u.contigID.txt`. Again, I got an empty output with only showing the headers of five columns.

I also tried to modify the demovir.sh script by removing all previous parts I've completed and retaining the last command `Rscript $DIR/demovir.R`then ran `modified.demvir.sh my.viral.contigs.fasta threads`again, but still got the same result.

Hopefully anyone can do me a favour to check and resolve this issue.

Many thanks.

Ernie

---

Quoting the GitHub repository:

> "Democratic taxonomic classification of viral contigs to Order and Family level
>
> When performing metagenomic sequencing of Viral-Like Particle (VLP), the majority of returned sequences often bare little to no homology to reference sequences - Viral Dark Matter. Frequently it may be useful to know which viral taxonomic group these novel viruses are likely to belong to as this will give information about nucleic acid type, size and behaviour."

# Answer
> maybe you can open the demovir.R file with Rstudio and execute line by line to check where the failure is. - @zorbax

Thanks for your suggestion. I eventually got a correct final output already.

The key issue was Prodigal would add extra strings after my original contig IDs, so it would cause R not recongnising my viral contigs. Therefore, I added a flag `-trunclabels` in Usearch (after e-value flag) to truncate my contig names. Finally, I re-ran the demovir.sh in ubuntu VM and it has worked well.

Ernie
> 2 votes

---
Tags: sequence-annotation, covid-19, metagenome, taxonomy, virus
---


# Question
Title: BinaryCIF vs MMTF formats, which one to choose?
Both file formats are advertised as more memory/parsing efficent than old PDB format and its successor mmCIF format, but I can't understand the context (I'm a developer, but pretty new to the field of bioinformatics).

Is BinaryCIF a later iteration of MMTF concepts/ideas (so it's a new competing format) or is the purpose of the two formats different ?

Refs :

# Answer
> 6 votes
So binaryCIF is newer and is heavily inspired by MMTF.

NGL, the most "loved" JS library (IMO) to show protein fetches proteins as MMTF, while Mol\* uses bCIF —Mol\* has the official repo of the bCIF format. Alex Rose, developed the former at the RCSB PDB, and the latter in collaboration with the PDBe. He also worked with Antony Bradley, the main author on the MMTF format, who left the RCSB PDB, some 4+ years ago. So there is definitely no competition there.

* The MMTF download of the RCSB works well: https://mmtf.rcsb.org/download.html
* The MMTF format did not catch on —okay, mmCIF is only slowly replacing PDB, which is still the only accepted format for analyses—, hence the need for a bCIF format
* The libraries to do the conversion aren't without glitches
* In the distant future bCIF will be the sole binary format
* mmCIF \<--\> bCIF conversion is easier

However, I really ought to stress that mmCIF is a better choice to bCIF or mmTF. And even then many tools still hate mmCIFs and use PDB format. If memory footprint is not really a problem, then stick with mmCIF. If you are running a webserver, you may be sending data as gzipped making this a moot point.

Last year for my site, michelanglo, I was faced with a dilemma: Mol\* has dotted lines for gaps and does not crash the browser with a virus or someone wanting to show hyperball representation on all residues. However, it was not documented (and then glitchy) so I firmly stuck with NGL, which is excellent. For file formats, I stick with PDB and a metadata dictionary as I know the format well and is despite having the jankiest syntax and no real metadata support (REMARK... ha!), does not break isopetides, disulfides or for no apparent reason.

# Answer
> 2 votes
Consider **mmJSON**.

Each of the three PDB sites had the same problem: using mmCIF files in web-based viewers is inefficient. To solve this problem three new file formats were introduced:

* MMTF at RCSB,
* BinaryCIF at PDBe,
* mmJSON at PDBj.

**MMTF** has the smallest files – it's aggressively optimized for size. It is based on MessagePack, but on top of it uses a number of custom compression schemes, which make the format rather complex. Usually it's not a problem – MMTF developers made an effort to provide libraries for reading MMTF files in the most popular programming languages. But you need to install these libraries. To make the file smaller, MMTF doesn't contain all the metadata from mmCIF, which may be a blocker in some cases. What is nice is that the MMTF developers reached out to other developers in this field, organized some gathering and they were taking comments and suggestions from the outside. As the result, a few popular programs such as PyMOL and JMol can also read MMTF.

**BinaryCIF** is clearly based on (or at least inspired by) MMTF. So it's also built on top of MessagePack with a number of custom compression schemes. It includes all the mmCIF metadata, which makes the files slightly bigger, but I think it's a good trade-off. Initially there was only TypeScript code for reading it, now I see that there are also implementations in Java in Python.

**mmJSON** is 1:1 equivalent of mmCIF, but uses JSON instead of CIF. There is no extra compression added, but having the data from mmCIF re-ordered column-wise makes gzip compression much more efficient (this was also noted by MMTF developers). If the file is stored as json.gz or downloaded in web browser (with the usual compression between the server and the client) its size is 30-40% of gzipped mmCIF. It's not as small as MMTF and BinaryCIF (which are several times smaller than mmCIF) and this is the only downside.

The advantage is that you can read in it any language like all other JSON files. No need to install/update MessagePack + MMTF/BinaryCIF libraries.

A few years ago I compared reading MMTF and mmJSON in Python – reading mmJSON was order of magnitude faster. It may be very different with libraries in other languages, but generally JSON parsers are fast. The fastest ones parse GBs per second. Reading binary formats can be faster, but decompression schemes can also make it slower. Anyway, all the three formats are fast enough to read.

Sadly, I don't think that mmJSON is as popular as it should be. Perhaps I'm the only person using it outside of PDBj. That's why I'm writing this answer.

The mmJSON files are available at: https://pdbj.org/help/mmjson

---
Tags: file-formats, pdb, 3d-structure
---


# Question
Title: How to optimize the number of amplicons ordered for a PCR wet experiment with several genomic ROIs?
Imagine we have an experiment for which we would like to minimize the number of purchased primers. Let's look at the layout of a single primer below:

Assume we have a precomputed set of all necessary amplicons and we are trying to reduce the costs and find the minimum set of amplicons that tile all regions of interest (ROIs) with some flanking.

The algorithm seems to work, but I am not sure it is ideal, so I would like to know how to assert it finds the minimal covering set of amplicons and how to improve the runtime complexity.

First, we define the underline representation of a primer/amplicon (I use the terms interchangeably), and an auxiliary function for maintaining the same colors among different plotting sessions.

```
from dataclasses import dataclass
import numpy as np

def rand_color():
   return "#" + "%06x" % np.random.randint(0, 0xFFFFFF)

@dataclass(eq=True,frozen=True)
class Primer:
    left_end:  int
    right_end: int
    left_roi: int
    right_roi: int
    color: str
    level: int

```

We then create a random set of primers:

```
N = 25; max_coor = 2_000
min_primer_len, max_primer_len = 200, 500
min_required_flank_len, max_required_flank_len = 10, 50
roi_len = 20
primers = []
for i in range(N):
  left_end = np.random.randint(0,max_coor-max_primer_len)
  primer_length = np.random.randint(min_primer_len,max_primer_len)
  required_roi_flank = np.random.randint(min_required_flank_len,max_required_flank_len)
  left_roi = np.random.randint(left_end+required_roi_flank,left_end+primer_length-required_roi_flank-roi_len-1)
  right_roi = left_roi + roi_len
  primers.append(Primer(left_end, left_end+primer_length, left_roi, right_roi, color=rand_color(), level=i))

```

The algorithm for finding the minimal primer covering set is as follows:

```
primers.sort(key=lambda p:p.left_end)

def does_a_cover_b(a,b, flank=0):
  return a.left_end< b.left_roi-flank and a.right_end> b.right_roi+flank

p_tmp = primers.copy()

cover = set()
for p in p_tmp:
  add_to_cover_set = True
  for c in cover:
    if does_a_cover_b(c,p):
      add_to_cover_set = False

  if add_to_cover_set:
    cover.add(p)

```

Now `cover` is supposed to be a minimal subset of primers that covers all ROIs. To visualize, we use `plot_primer_set` twice:

```
import matplotlib.pyplot as plt
from shapely.geometry import Polygon

def plot_primer_set(primer_set):
  plt.figure(figsize=(4,4), dpi=120)
  for c in primer_set:
    p_all = Polygon([(c.left_end,c.level),(c.left_roi,c.level),(c.right_roi,c.level),(c.right_end,c.level)])
    x,y = p_all.exterior.xy
    plt.plot(x,y, markersize=7, marker='o', linewidth=5, color=c.color)
    plt.yticks(np.arange(0, N, 1))
    axes = plt.axes()
    axes.set_xlim([-100, max_coor+100])

plot_primer_set(primers)
plot_primer_set(cover)

```
1. Is this the optimal and correct minimal set?
2. My complexity is $O(N^2)$, but I expected it to be $O(N \log{N})$, what am I missing?

# Answer
> 1 votes
I ended up implementing the following greedy algorithm that should work fairly good for any binary relation even if it is not transitive nor symmetric.

Assume you have the following binary relation (you may plug in any function you want):

```
def does_a_cover_b(a,b, min_dist=30, max_dist=130):
  return min_dist < b.left_roi - a.left_end < max_dist or min_dist < a.right_end - b.right_roi < max_dist

```

Using the auxiliary functions to define the binary relation, draw the graph, and build the graph:

```
import networkx as nx

def does_a_cover_b(a,b, min_dist=30, max_dist=130):
  return min_dist < b.left_roi - a.left_end < max_dist or min_dist < a.right_end - b.right_roi < max_dist

def draw(G):
  pos = nx.spring_layout(G,scale=2)
  nx.draw(G, pos, with_labels=True, font_weight='bold', node_size=750, width=2)

def build_a_graph(primer_list):
  G = nx.DiGraph()

  for i,p_i in enumerate(primer_list):
    G.add_node(i, data=p_i)

  for i,p_i in enumerate(primer_list):
    for j,p_j in enumerate(primer_list):
      if i == j:
        continue
      if does_a_cover_b(p_i, p_j):
        G.add_edge(i,j)
  return G

```

we may implement the following algorithm:

```
p_set = set(primers)
cover_set = set()
covered_primers = set()
while len(covered_primers)<N and len(p_set)>0:
  G = build_a_graph(p_set)
  nodes =sorted(G.out_degree, key=lambda item: item[1], reverse=True)
  node = nodes[0][0]
  absorber = G.nodes[node]['data']
  p_set.remove(absorber)
  cover_set.add(absorber)
  covered_primers.add(absorber)

  for edge in G.out_edges(node):
    primer = G.nodes[edge[1]]['data']
    p_set.remove(primer)
    covered_primers.add(primer)

print(f"Used a cover set of length {len(cover_set)} primers to cover the total primer set with {len(covered_primers)} primers")
plot_primer_set(primers)
plot_primer_set(cover_set)

```

---
Tags: python, algorithms, pcr
---


# Question
Title: Why do we count reverse compliments when counting kmers in a single DNA strand?
When looking for repeat patterns in a DNA sequence, we can look for the pattern and its reverse compliment with up to `d` mismatches. However, why do we look for reverse compliments if we're analyzing a single DNA strand?

I thought about it and my guess is that this lets us analyze both strands of DNA. For example, if we find the reverse compliment at position `x`, that means that the other strand of DNA has the pattern at position `x - pattern.length`.

Is that correct?

# Answer
> 3 votes
When we sequence genomes, the individual reads come from both strands, and that also applies for assembled sequences\* (contigs / scaffolds). Therefore for most of the genome assemblies out there, strandeness of individual scaffolds is actually an arbitrary property, and reverse complementing any sequence will generate an equally valid representation (with reverted coordinates).

The coordinates of the reverse complementary k-mer matches depend on the notation. Usually, the notations are still in the coordinates of the reference, for example, blast gives for reverse complimentary strand matches the `from` subject coordinate greater than the `to` coordinate. Some others are denoting "the leftmost" position and strand encoded in a different column (e.g. sam).

\*perhaps with the exception of chromosomal level assemblies

---
Tags: k-mer
---


# Question
Title: How to copy only certain counts for genes in tsv file to new file in linux
Hi there I have generated a counts table of samples I need to compare by differential expression analysis.

The layout of the counts table is as follows:

| Gene\_id | Sample\_A\_r1 | Sample\_A\_r2 | Sample\_B\_r1 | Sample\_B\_r2 | Sample\_AB\_r1 | Sample\_AB\_r2 |
| --- | --- | --- | --- | --- | --- | --- |
| KLTH001 | 0 | 0 | 2 | 5 | 8 | 3 |
| KLTH002 | 0 | 0 | 5 | 3 | 9 | 7 |
| YL99812 | 6 | 4 | 0 | 0 | 12 | 14 |
| YL99813 | 9 | 7 | 0 | 0 | 15 | 19 |

So here, for the differential expression analysis, I only want to compare for example the genes starting with "YL-". So in a new file I want the rows including the first row which is "Gene\_id" and the sample names e.g. "Sample\_A\_r1" and then only the rows with genes starting with "YL-" in their names. This will allow me to only compare the changes between Sample\_A alone versus when it was combined with Sample\_B, or in the table Sample\_AB.

If anyone can please give me an indication of how to do this in linux I would really appreciate it. I have seen other similar posts, involving `grep` and `awk` I am just not sure of exactly how to do this. I am still quite new to linux. Any help will be appreciated.

# Answer
Something like this should do job:

```
# creating a new file with the first line of the original
head -1 foo.txt > bar.txt

# ">>" appends to the file on the right hand side
# "^" dictates that the search term should be at the beginning of a line
grep "^YL" foo.txt >> bar.txt

```
> 3 votes

# Answer
One way using :

```
awk 'NR==1 || /^YL/ { print $0 }' file.tsv

```

The default action is to print the line/record, so all you need is the condition:

```
awk 'NR==1 || /^YL/' file.tsv

```
> 1 votes

# Answer
Alternatively, you can use `grep` like so to remove all lines starting with `KLTH`:

```
grep -v '^KLTH' in_file > out_file

```

For example:

```
echo "Gene_id\nKLTH001\nYL99812" | grep -v '^KLTH'

```

Prints:

```
Gene_id
YL99812

```

Here, GNU `grep` uses option:
`-v` : Print lines that do **not** match.

In the regular expression (regex):
`^` : Match the beginning of the line.
> 1 votes

---
Tags: linux, filtering, awk, grep, featurecounts
---


# Question
Title: Given two ADN sequence, can there be multiple optimal alignement? And how to optimally find them?
Given two sequences, I know we can determine optimal alignment of sequence with the Needleman–Wunsch algorithm. But is it possible that two alignments have multiple optimal alignment (I guess it would mean that multiple alignement have the same score)?

And if possible, how would you find them via an algorithm? Just looking for a general idea on how to go about it, I can't see how you could do it beside by checking for all possible combinations.

# Answer
As Needleman-Wunsch is a Dynamic Programming algorithm you do need to complete the entire matrix to find the optimal solution and it is possible to have multiple best alignments.

> The algorithm assigns a score to every possible alignment, and the purpose of the algorithm is to find all possible alignments having the highest score.

https://en.wikipedia.org/wiki/Needleman%E2%80%93Wunsch\_algorithm
> 1 votes

---
Tags: sequence-alignment, algorithms
---


title: ViolinPlot x-axis

text:
# Question
Title: ViolinPlot x-axis
I'm trying to set the x-axis in my Violin Plot to go by each patient rather than it being scattered.

```
VlnPlot(obj, features = c("classical_basal"), pt.size = 0) + NoLegend() + theme(axis.text.x = element_text(angle = 90, hjust = 1))

```

The above code gives me the x-axis but there is no pattern or order. I would like the x-axis to be grouped with the HTB, for example ALL of the HTB2867 first, followed by HTB2876, etc...

# Answer
Seurat does violin plots using ggplot, which orders based on the factor levels. You need to re-factor your grouping category so that the levels matches your preferred order:

```
# default level sort order is alphabetical
obj[["grouping"]] <- factor(unlist(obj[["grouping"]])) # refactor
VlnPlot(obj, group.by = "grouping", features = c("classical_basal"),
        pt.size = 0) +
  NoLegend() + theme(axis.text.x = element_text(angle = 90, hjust = 1))

```
> 2 votes

# Answer
It is difficult to say exactly what you would like to achieve without an example.

But if you would like an individual axis based on a factor like patient maybe you want to use facet\_wrap for the Patient variable like in the example below?

library(ggplot2)

data <- iris

ggplot(data, aes(x=Sepal.Width, y=Sepal.Length, fill=Species)) +
  geom_violin() +
  facet_wrap(~Species, scales="free_x")
```

> 0 votes

---
Tags: r, seurat, single-cell, ggplot2
---


title: NMF clustering not giving results

text:
# Question
Title: NMF clustering not giving results
I'm trying to run the Non-Negative Matrix Factorization(NMF) for my Gene Expression dataset which was originally in matrix form. But it throws errors as given below:

```
expr <- read_excel("./expr.xlsx")
expr<- as.matrix(expr)

#removing the first column containing gene symbols before clustering
expr_rm <- expr[,-1]
class(expr_rm)
[1] "matrix" "array"

```

So converted these expr\_rm to numeric:

```
expr_rm <- as.numeric(expr_rm)

#removing negative values for running NMF
eps <- .Machine$double.eps
expr_rm[expr_rm<=0] = eps

estimate <- NMF::nmf(expr_rm, 3, nrun = 10, method = "lee", seed = 123)

```

> Error: NMF::nmf - 10/10 fit(s) threw an error. Error(s) thrown: run #1: 'x' must be numeric Timing stopped at: 1.95 0.14 4.4\*

```
cluster <- cbind(names(NMF::predict(estimate)), NMF::predict(estimate))

```

> Error in h(simpleError(msg, call)) : error in evaluating the argument 'object' in selecting a method for function 'predict': object 'estimate' not found

```
result <- list(estimate, cluster)

```

> Error: object 'estimate' not found

**For clarity, my expr\_rm has row indices as 1, 2, 3, 4,.. which indicates different genes and column headings as samples.**

Please suggest a solution for this soon!

Additional information (requested from comments):

```
>summary(expr_rm) - before conversion
 Sample_C44_0001.genes.results Sample_C44_0002.genes.results Sample_C44_0003.genes.results Sample_C44_0005.genes.results
 Length:145   Length:145 Length:145  Length:145  Class :character

 > summary(expr_rm) - after numeric type conversion
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
  0.000   1.808   4.264   4.220   6.333  16.246

```

If in expr\_rm is data.matrix() also, it gives the error like:

```
> median_centered_data<- expr_rm - med
Error in expr_rm - med : non-numeric argument to binary operator

```

# Answer
> 0 votes
*Answer from @maximilian-press, converted from comment:*

Sometimes casting to `as.numeric()` will convert from two-dimensional matrix to one-dimensional vector. Based on the output you put into comments (which would be helpful in the question), I wonder if this has occurred and caused an issue. For `data.matrix`, I suggest reading docs for that function- it should have returned only numeric, such that your `med` variable is likely not numeric, or there is a problem with your data. Possibly the tidyverse `tibble` created by `read_excel` might be the issue. I don't know enough tidyverse to say what exactly is going on.

---
Tags: r, genomics, gene-expression, clustering, modelling
---


# Question
Title: Gwas analysis studies using the result of differential analysis
I want to do a GWAS analysis using the differentially expressed genes that I got from a differential analysis. The purpose is the identification of biomarkers.

It's my first time using GWAS, can you please suggest to me a pipeline to use?

# Answer
> 0 votes
I am not an expert on recent packages but I would suggest using Matrix eQTL for purposes of the statistical analysis.

As mentioned in the comments I think that it is worth testing *all* transcripts for eQTL associations. This will additionally give some idea of the false discovery rate. It is possible that testing "candidate gene" differentially expressed genes will lead to a lot of false positives relative to the background level of associations.

---
Tags: gwas
---


# Question
Title: How to find differential expressed genes within pseudotime trajectory with Seurat cluster?
Did anyone know how to find the differential expressed genes within pseudotime trajectory with Seurat cluster? In monocle tutorial, we can use BEAM to find the differential expressed genes between the two branches even with the states difference. But how can I find the differential expressed genes with seurat culster enven with the genotype diffence? Thank you so much!

# Answer
> 1 votes
You probably want approaches such as tradeSeq.

For Seurat to tradeSeq see also this issue.

# Answer
> 0 votes
Not sure if it helps, but there's a tutorial for how to switch between Seurat and monocle here. The basic approach is to convert it to a `cell_data_set` object, then process that object with monocle:

```
erythroid.cds <- as.cell_data_set(erythroid)
erythroid.cds <- cluster_cells(cds = erythroid.cds, reduction_method = "UMAP")
erythroid.cds <- learn_graph(erythroid.cds, use_partition = TRUE)

```

Pseudotime values can be returned back to the seurat object after processing with monocle via AddMetaData:

```
erythroid <- AddMetaData(
  object = erythroid,
  metadata = erythroid.cds@principal_graph_aux@listData$UMAP$pseudotime,
  col.name = "monocle_pseudotime"
)

```

> Thank you for your reply. So I need to do the code you sent from erythroid.cds \<- as.cell\_data\_set(erythroid) ...to ...col.name = "monocle\_pseudotime"), then use the findmarker function of Seurat to find the differential expressed genes? It's that correct?

It's unclear what you're asking, but what you describe here will produce a list of differentially expressed genes based on the cluster definitions of monocle. As far as I'm aware, these clusters don't depend on the pseudotime values.

---
Tags: single-cell, monocle
---


# Question
Title: Statistical or performance test for genomic intervals between 2 samples
Is there an appropriate statistical or performance test to perform when I want to compare genomic intervals between 2 samples. Sample A is the the benchmark of comparing the intervals and sample B is the sample in question. What I am really looking for is not only how many intervals from sample B did indeed overlap with sample A but additionally what intervals from sample B has intervals that don't overlap with sample A and also those intervals in sample A that don't overlap with B and being able to measure the performance of how accurate sample B intervals overlap with sample A. This in some sense, can be translated into the following:

1. True positive: intervals from sample B that overlap with sample A.
2. False positive: intervals from sample B that don't overlap with sample A.
3. False negative: intervals from sample A that don't overlap with sample B.

If I consider the above schema, then I will not be able to do any of the performance and statistical test (ROC curve or accuracy) because it is not possible to know what the true negative and the true negative are needed parameters. Then, what would should I be using to be able to incorporate the above schema and still apply performance and statistical testing.

I highly appreciate any help or pointers in advance.

# Answer
I would recommend looking into BEDtools, assuming that you are talking about genomic intervals, which I *think* you are based on the question. Possibly a diagram would help.

BEDtools has a suite of statistical tools for evaluating overlaps, including Fisher Exact Tests and also approaches based on the relative distance of 2 sets of intervals. So as long as your genomic intervals are defined by BED files and a genome file, this should be quite easy, e.g.

```
bedtools fisher -a A.bed -b B.bed -g genome_file.genome
# will print p-values for left, right, and two-tailed tests

```
> 1 votes

---
Tags: statistics, genomics
---


# Question
Title: Dictionary and index of vcf for base recalibration step
I need suggestions in creating index and dictionary with vcf files. For the base recalibration step, I downloaded Homo\_sapiens\_assembly38.known\_indels.vcf.gz from the given link: https://console.cloud.google.com/storage/browser/genomics-public-data/resources/broad/hg38/v0;tab=objects?prefix=&forceOnObjectsSortingFiltering=false

I am using GATK version 4.2.2.0

When I try to create index and dictoionary of the vcf file, I get the following error: no suitable codecs found. I also tried to check the vcf file using ValidateVariants and I get the same error: no suitable codecs found.

Can anyone suggest which indels vcf file is used for masking sites in the base recalibration step. Has anyone come across this above error and how did they resolve the issue.

Thanks

# Answer
> 1 votes
I have also encountered this error before in cases where there was something wrong with the VCF file I was supplying.

* Double check your command to ensure that you are feeding it the file you believed you were and verify that file is a vcf (gzipped or otherwise)
* If you know how, please compare MD5 or SHA checksums of your copy vs. the original to ensure you have an uncorrupted copy.
* If this is something you haven't dealt with before, a less effective version would be to check and make sure that the file size on the server is the same as the file size on your machine.
* Additionally, you may want to just use the `less` command to look inside the file and see if anything appears to be corrupted (you may need to decompress it for that purpose).
* It may also be worth trying to `gunzip` or `gzip -d` commands to decompress the file (assuming you have the drive space). This will both ensure that decompressing the file on the fly is not causing any grief to GATK (which is unlikely) and that the file decompresses without error (which would also be an indication of corruption).
* Lastly, you could try downloading a fresh version. Please let me know if this helps and you discover the source of the error.

---
Tags: vcf, ngs, gatk, data-preprocessing, variants
---


# Question
Title: If I have data from a full genome sequence, how can I search for a known pathogenic variant of a known gene?
I have a CS background, but am a bioinformatics neophyte

I did a full genome sequence which provided me with ~100GB of files (SNP VCF, Indel VCF, BAM, Indel TBI, SNP TBI, BAM BAI, CNV VCF, CNV TBI, CV VCF, SV TBI, FASTQ R1, and FASTQ R2).

I'd like to search the genome for a pathogenic variant of a known gene. What tools do I need to do this, or where do I start?

# Answer
> 4 votes
I'm going to try and keep this answer from going overboard while still helping you find your feet. The good news is that it sounds like you have started off with some of this work already done for you, and based upon the files you have above, I will try to indicate what you have already and also where you might be able to go with this. The term you should be entering into Google as the generally accepted term for this pipeline is `variant calling` or possibly `variant curation` for the downstream steps.

1. This process starts off with one or more files containing raw observed sequences from a sequencing device (often, but not always, Illumina) This files will usually be in the FASTQ format, but historically I have seen another format called QSEQ. You will likely never see this format and I feel old just bringing it up, but if you ever encounter it, just convert your files to FASTQ immediately. (Files entering here: FASTQ, possibly QSEQ if it was generated well before Carly Rae Jepson suggested that we call her, maybe)
2. Although not a hard requirement, I recommend running something like `fastqc` or `multiqc` or some other form of quality control on your reads and making sure that the reads are of satisfactory quality with no anomalous metrics on them. You may or may not decide that you need to do some clean-up on the reads if you determine that there might be adapter contamination or that some reads are partially/entirely garbage quality. (Files entering here: FASTQ, reference genome in FASTA. Files generated: a quality metrics report and possibly a cleaned up FASTQ)
3. Once you have your raw reads and consider them to be of satisfactory quality, they will need to be aligned to a reference genome. Read alignment is a critical, computationally intensive process that is required for many bioinformatics pipelines, meaning that there are many aligners out there that have been optimized for many different situations. If you are running a variant calling pipeline on Illumina short reads from genomic DNA, a current version of BWA should be a solid choice for you here. (Files entering here: FASTQ. Files generated: SAM and/or BAM files with aligned reads)
4. Now that alignment is complete, we can start getting reads for the real variant calling process by doing any preprocessing needed. The most popular tool(s) for this historically were Picard Tools and GATK, but those two tool sets were used together for so long that they were just begging to be merged. The recent versions of GATK have made it so and Picard tools are now rolled into and called directly from GATK. If you are using an older GATK version, at a minimum you will want to recalibrate the base quality scores to try correcting biases that your sequencer had when assigning confidence scores to individual base calls (they all have this, trust me) and realigning reads around detected insertion/deletion (indel) sites. Newer versions of GATK don't need the indel realignment steps due to the haplotype calling algorithm going even further and reassembling reads around these complex variants. Those are just the minimal preprocessing steps, you may also need to assign read groups if they aren't already in the BAM files and possibly some other stuff (most likely GATK will let you know what's missing as you progress). (Files in: BAM, reference genome and "gold standard" variants in VCF format. Files out: preprocessed BAM)
5. Actual variant calling can happen now, and this is likely going to happen using HaplotypeCaller (mutect2, if you're dealing with cancer/somatic changes). If this is your intention, you can look up the latest and greatest practices at https://gatk.broadinstitute.org/hc/en-us/articles/360035535932-Germline-short-variant-discovery-SNPs-Indels- for the most current guidance. (Files in: preprocessed BAM, reference genome. Files out: variant calls in VCF, possibly an active sites BAM file to visualize reassemblies around active indels)
6. Filtering of variants should now happen. At this stage, you may be very application dependent, but essentially this process should assign some variants to being high confidence and others to being potentially technical artifacts. I have seen people skip this step either due to lack of interest or concern that it will filter out a real variant of interest, but I don't recommend that approach and have seen more people missing real variants due to being inundated with low confidence variants that should have been otherwise filtered out.
7. Variant curation can now happen where one takes their high confidence variants and figures out (or usually has a computer program figure this out) which genes (if any) are being affected by the variant, what the chance will be to the protein coding as well as if the variant might affect splicing or some aspect of gene regulation or structure (such as an out of frame indel). It is often also useful here to compare to some database of mutation frequencies to filter out common SNP variants as potential cause for rare conditions.
8. Identification of candidate mutations based upon the curation results. Identifying the mutations that are in genes that could plausibly control the observed phenotype, are uncommon, and are potentially causing pathogenic changes to the gene or its product (protein).

It looks like you already have VCF files, so at least the initial phase of variant calling has been finished for you already. I would recommend checking out something like Annovar or Varient Effect Predictor (VEP) that can run over your (hopefully filtered) VCF and assign a predicted gene and effect to a variant. Check both the SNP and Indel VCF files, as indels, while more rare, can have profound effects on a gene.

From your description, it sounds like you wanted to check a few genes in particular, and after the variant effect preduction, you should have the needed information to do this.

---
Tags: vcf, software-recommendation, genetics
---


# Question
Title: Meaning of order of alleles in VCF GT field for unphased genotypes
Order of alleles in genotype field (on VCFs) can be used to infer haplotypes

However, what about unphased genotypes? I have a VCF file where I only have unphased genotypes:

```
0/1
1/0
1/1
1/2
2/1

```

However, I have, for example, both `0/1` and `1/0`. Does the order of alleles have some different meaning with unphased genotypes, or genotypes with same alleles (e.g., `0/1` and `1/0`) can be considered synonymous?

# Answer
In the case of unphased VCF `0/1` and `1/0` are exactly equivalent and the order has no meaning.

Throughout the VCF specification, it is always encouraged to use numerical precedence in order to order the genotypes, e.g. `0/0`, `0/1`, `1/1` would be the preferred representations for diploid bi-allelic genotypes.

However it is not strictly enforced in the VCF specification as for now and you could have `1/0`, `2/0`, etc, and bcftools will not complain parsing them. The presence of two representation in the same VCF might indicate that the data comes from multiple variant callers or has been through manual editing, since for example variant callers will always use one representation only (numerical precedence usually)

Also, the order does impact how you read some of the sample information fields, for example the AD (allelic depth), for which the values are read following the order present in the GT field.

I had a look at the VCF 4.2 and 4.3 spec, but indeed is not so precise there on these details.
> 3 votes

---
Tags: vcf, variant-calling, genotyping, phasing
---


# Question
Title: How do you search for new enzymes that are more stable for handling, immobilization?
Noob here. I get that I should deduce what characteristics the ideal new enzyme should meet, and then use tools such as PDB and blast to compare to the old enzymes, and use other tools such as pymol, but I have no idea on how to approach this. Where should I look for tutorials?

# Answer
> 1 votes
From your title, you are after better immobilisation, i.e. your enzyme behaves poorly when bound.

* Poor thermostability: Engineering enzymes for better thermostability, involves predicting mutations that improve its Gibbs free energy (∆G, negative potential = better). Finding mutations that improve a protein by going through the mutational landscape is called scanning. Rosetta `pmut_scan` is a command line tool to do this, while Pross is a really nice webapp that also does it.
* Poor catalysis: if your attachment is close to the active site, that is a design flaw and you need to attach on the other terminus or within a loop. Circular permutation may be an option.
* Aggregation: your enzyme aggregates due to pI close to 7 (see protparam) or has a face that is hydrophobic (see the APBS plugin shipped within PyMOL). Rosetta `supercharge` is great at suggesting mutations to decrease the pI, but if doing it crudely by changing Gln -\> Gln & Asn -\> Asp, do check for surface exposure.

---
Tags: sequence-alignment, blast, pdb, phylogenetics, pymol
---


title: Convert VCF file to mpileup.txt

text:
# Question
Title: Convert VCF file to mpileup.txt
I am working on an iterative analysis that uses orthologous pipelines that require mpileup.txt files as input for a visualization step. This requires me to convert VCF files to mpileup.txt.

This experiment uses error-prone-pcr and follows a somatic variant discovery workflow that was designed by a 3rd party. I have written and implemented several orthologous pipelines with small differences from the original workflow in order to produce a truth set (variants identified consistently by all of the pipelines). The 3rd party used a workflow that calls Samtools to produce mpileup.txt files as the final output and necessitated the development of a custom visualization script. My orthologous pipelines use mutect2 to produce vcf files that need to be compared to the mpileup.txt files.

If I could convert my vcf files to mpileup.txt files without loss of information it would save me significant time cost associated with writing a similar visualization script for vcf files and prevent bias involved in converting the 3rd parties output to vcf format. There are plenty of examples for using samtools and bcftools to produce an intermediate mpileup file that is then converted to a finalized fitered vcf file. But I can't find much direction for converting a vcf file to an mpileup. Can anybody provide direction for this task?

# Answer
Going VCF to mpileup is not really something one does or can do. The mpileup should be generated from a BAM or SAM file or something else that has raw, unfiltered read alignments in it. The VCF just has data on variant sites, and usually just variant sites that passed a certain likelihood threshold at that. With the VCF, the best you can do is generate a very limited pileup based upon the VCF-reported counts of reads at each variant locus with no information in between reported variant loci. This would not be a true mpileup and would most likely cause an error with the downstream process or worse, not cause an error and give a result that is divorced from the reality of your sample.
> 3 votes

---
Tags: rna-seq, samtools, variant-calling, gatk, bcftools
---


# Question
Title: Difference between bioinformatics and cheminformatics
Anyone can answer, and should be in simple way

# Answer
There are four terms to consider, all are computational

* cheminformatics — Big data to explore relationship between small molecules (e.g. enumeration and machine learning stuff)
* computational chemistry ("compchem") - thermodynamics-based methods to explore small molecules (e.g. docking)
* bioinformatics — big data primarily for exploring -omics data (genomics, proteomics, transcriptomics, metabolomics)
* computational biochemistry - thermodynamics-based methods to study structure of protein
> 1 votes

---
Tags: genomics, cheminformatics
---


title: blastn returning an inferior alignment

text:
# Question
Title: blastn returning an inferior alignment
In my work I've found that in circumstances where there is a mismatch near the edges of the query sequence, blastn prefers to return a shorter contiguous alignment, rather than allowing for a mismatch near the 5' or 3' end of the query sequence.

For example, consider the subject sequence (bold indicates alignment with mismatch in brackets):

GACCAACCTG**CTGGATCTGCTGATCGGCG\[G\]T**TTGCGAAATA

and query sequence:

CTGGATCTGCTGATCGGCG\[A\]T

Blasting this will return the 19bp contiguous alignment upstream of the mismatch, rather than returning a 21bp alignment with a single mismatch at the G/A base at position 20 in the query.

I'm using the command line version of blastn, and I've tried messing with the penalties as well as setting `-task blastn-short`, but nothing seems to fix this issue.

I speculate that the various gap and mismatch penalties are working out such that a shorter contiguous alignment is preferred over a slightly longer alignment with a mismatch, but I haven't been able to figure out how to modify the penalties in such a way as to change this preference.

# Answer
> 0 votes
Here's my working solution,

So I found a choice of parameters that resolves this, `-reward 3 -gapopen 5 -gapextend 5`.

---
Tags: sequence-alignment, blast, blastn
---


# Question
Title: Consensus on how to select SNPs for creating a phylogeny when the full dataset is too large
I have got a quite large VCF file (about 20k samples with 160k SNPs after filtering for quality etc.) and I would like to get a phylogeny for it. However, the whole dataset's too large for my computational resources (unless I go for fasttree which would be my last resort). Therefore, I need to subset the SNPs and there's plenty of possible ways to do this (dump all rare SNPs with MAF below a threshold, only use SNPs in coding regions, only use missense SNPs; just to name a few).

Is there a consensus (or at least some comparative work or recent review; I found nothing) on what method works best and gives the tree that's closest to using the full dataset?

# Answer
> 2 votes
I will give an answer that is specific to closely-related microbial genomes. Core-SNPs, SNPs which are within conserved and likely orthologous sequences, have been shown to be relatively informative for constructing phylogenies of large sequence sets\[1\].

Parsnp is the most popular tool for constructing core-genome alignments on microbial genomes, and extracts the core SNPs from the alignments. Parsnp then uses raxml or FastTree in order to produce a phylogeny from the core SNPs. I will say, though, that 20k sequences may be too many for Parsnp to run on (depending on how large the sequences are and your available resources).

Disclaimer: I am the maintainer of Parsnp (:

# Answer
> 1 votes
The subset that gives you the tree closest to the one you obtain with the full dataset does not necessarily is the tree nearest to reality... It could be that you get a more realistic result using only partial data for multiple reasons. Your best chance is in my opinion, as it usually is, making lots of tests and contrasting hypotheses.

---
Tags: phylogenetics, vcf, phylogeny
---


# Question
Title: Removing columns from meta.data in Seurat
I'm trying to remove columns from the meta.data from a Seurat object. Is there a way to do that or would I have to re-create the Seurat object without that column?

# Answer
Set the column you don't want to NULL

```
x$remove.column <- NULL

```
> 2 votes

---
Tags: r, seurat
---


# Question
Title: Snakemake: Accept any of a possible set of inputs (any one being sufficient)
I have a little pipeline that I want to execute using Snakemake.

One of my rules (let's call it `create_somefile`) creates a file `somefile_{n_dims}.json`, and when running the pipeline can specify that I want to run the pipeline for a pre-defined set of values for `n_dims`, which is why I specified this as a wildcard.

Now I want to have another rule, that takes as input **any of** these files. It can take either `somefile_1.json` or `somefile_2.json` or `somefile_32402348.json`, it doesn't matter which one as long as any of these exists. The rule then creates a file that does not depend on the number of dimensions, and thus doesn't want and need the n\_dims in its output.

When I try to have the wildcard only in the input, like this:

```
rule do_stuff:
    input:
        "somefile_{n_dims}.json"
    output:
        "general_file.json"
    shell:
        "echo {input} &&"
        "touch {output}"

```

Snakemake gives me the error `Wildcards in input files cannot be determined from output files: n_dims`. I don't want the output file to contain `n_dims` whatsoever, and unfortunately it's also impossible in snakemake to just create an additional throwaway `temporary(..)`-output that does depend on the n\_dims just so snakemake is happy (error message then: `ot all output, log and benchmark files of rule do_stuff contain the same wildcards. This is crucial though, in order to avoid that two or more jobs write to the same file.`.

Is there any way in snakemake to specify a rule's inputs as "take any of these", without wildcards?

Looking at the documentation, the only thing I found is the concept of `checkpoints` (https://snakemake.readthedocs.io/en/stable/snakefiles/rules.html#data-dependent-conditional-execution), however that doesn't seem to work as well - I tried making the rule `create_somefile` a checkpoint and changing the `do_stuff` rule to

```
rule do_stuff:
    input:
        lambda wildcards: checkpoints.create_somefile.get(**wildcards).output
    output:
        "general_file.json",
    shell:
        "echo {input} &&"
        "touch {output}"

```

however then I get the error `WorkflowError: Missing wildcard values for n_dims`.

Is there any concept in snakemake to allow me to do what I want? To have "wildcards" in inputs, that just stand for "execute this rule as soon as the input for any value of the wildcard is fulfilled, execute only once and be done with it"?

# Answer
I don't think Snakemake can follow a "use any of these" pattern with a unified dependency graph and without additional assumptions, but you can do it if you break the workflow into two pieces.

This will take whatever `somefile_*.json` shows up first, or throw an exception if there aren't any:

```
from pathlib import Path

def input_for_do_stuff(wildcards):
    try:
        return next(Path(".").glob("somefile_*.json"))
    except StopIteration as err:
        raise Exception("No n_dims JSON files found") from err

rule do_stuff:
    input: input_for_do_stuff
    output:
        "general_file.json",
    shell:
        "echo {input} &&"
        "touch {output}"

rule create_somefile:
    output: "somefile_{n_dims}.json"
    shell: "touch {output}"

```

For example if you have `somefile_2.json` already there:

```
$ snakemake -nrp
Building DAG of jobs...
Job stats:
job         count    min threads    max threads
--------  -------  -------------  -------------
do_stuff        1              1              1
total           1              1              1

[Sun Oct 10 12:22:14 2021]
rule do_stuff:
    input: somefile_2.json
    output: general_file.json
    jobid: 0
    reason: Missing output files: general_file.json
    resources: tmpdir=/tmp

echo somefile_2.json &&touch general_file.json
Job stats:
job         count    min threads    max threads
--------  -------  -------------  -------------
do_stuff        1              1              1
total           1              1              1

This was a dry-run (flag -n). The order of jobs does not reflect the order of execution.

```

With that, you'll never be able to run the workflow by just requesting the final output, which kind of breaks Snakemake's way of thinking about things. But at least it'll work like you want if you run it in two pieces. If having a default makes sense for your case, you could also do something like this:

```
def input_for_do_stuff(wildcards):
    try:
        return next(Path(".").glob("somefile_*.json"))
    except StopIteration as err:
        return "somefile_1.json"

```

That way, even if no files exist to start with, Snakemake can come up with a way to make the output you ask for.
> 1 votes

---
Tags: python, snakemake
---


# Question
Title: Pipeline to analyse small RNAseq and tRNAs?
I am a newbie in this world (just finished my master's on bioinformatics. As a practice, a researcher gave me 4 fastq files from a small RNAseq experiment, to see if I am able to reproduce their results but with an updated pipeline and tools (they used bowtie in 2014). I used fastqc for the quality check and I removed the adapters properly.

My first intention was to use Salmon for pseudo-align them, but after do some reading it seems not to be a good idea for such small RNAs.

Can anyone provide me with some light on this? Which alignment tool use and which parameters, and etc. Ideally a proper protocol/pipeline would be amazing.

Thank you in advance...

PD: I was thinking to use the HISAT2 with the full genome (or the pre-builded index from HISAT2 website).

Kind regards, and thank you in advance!

# Answer
> 1 votes
I have heard that the Small RNA Workbench from UEA can be an easy solution for a full small RNA pipeline. It is java based, has a GUI, and is supposed to be easy to install and run. I am not an expert in small RNA analysis, but I think this pipeline will certainly given you an idea about the steps needed for a small RNA analysis, it has been put together by experts.

# Answer
> 1 votes
piPipes is biased for piRNAs but should be a good starting point:

https://github.com/bowhan/piPipes/wiki/smallRNA-seq

Well documented and with many features.

# Answer
> 1 votes
As far as I'm concerned, bowtie 1 is still the method of choice since it rather mature and fast.

And yes, pseudoaligner such as salmon are probably not the best choice For once, they need a reference "transcriptome", which is usually available for RNA-Seq but not so much for small-RNAs. Also the way they handle multiple mapping locations is probably not ideal for small-RNAs that may have excessive number mapping location.

I've also written a simple snakemake-pipeline which simply performs some standard pre-processing and mapping against a reference genome along with QC and coverage tracks. It's still work in progress, but I'd be curious about any feedback.

# Answer
> 1 votes
Check nf-core for all the pipelines they have, including small-RNA seq. To run them, you have to use Nextflow.

---
Tags: trna
---


# Question
Title: Error to running "CIRIquant" function in CIRIquant package at python Language
"CIRIquant" is a Python package for detecting Circular RNA. I installed "CIRIquant" through the below command in terminal:

> pip install CIRIquant

My Linux distribution is Ubuntu 16.04 and my python version is Python 2.7.12. After installing, I got the below message: “Installing collected packages: numpy, scipy, scikit-learn, pysam, PyYAML, argparse, numexpr, CIRIquant Successfully installed CIRIquant-1.1.2 PyYAML-5.1.1 argparse-1.2.1 numexpr-2.6.9 numpy-1.16.4 pysam-0.15.2 scikit-learn-0.20.3 scipy-1.2.2”

So, When I call CIRIquant for running its function in Terminal, I got the below message in Terminal:

```
CIRIquant: command not found

```

Now, I don't know what is my problem? I appreciate it if anybody shares his/her comment with me. Best,

# Answer
> 3 votes
When I try this as a regular user on our Ubuntu 16 LTS system, I end up with `~/.local/bin/CIRIquant` from the pip install, so maybe just add that directory to your path:

```
$ export PATH="$HOME/.local/bin:$PATH"
$ CIRIquant -h
usage: CIRIquant [-h] [--config FILE] [-1 MATE1] [-2 MATE2] [-o DIR]
                 [-p PREFIX] [-t INT] [-a INT] [-l INT] [-v] [--version]
                 [-e LOG] [--bed FILE] [--circ FILE] [--tool TOOL]
                 [--RNaseR FILE] [--bam BAM] [--no-gene]

```

---
Tags: rna-seq, python
---


# Question
Title: Is there a sam flag for all&none of the reads?
a. Is there a SAM flag that specifies all reads? b. Is there a SAM flag that specifies none of the reads?

So that if I run samtools view -f (a) -F (b) the result will be all reads of the file (as if omitting the specification) ?

# Answer
> 3 votes
Instead of doing this kind of tweaking I would write my script in a way that these options are passed as strings that can be empty. Something like:

```
#/ dummy syntax just for illustration:
if(you_want_filtering){
  filter_f='-f something'
  filter_F='-F something_else'
} else {
  filter_f=filter_F=''
}

samtools view (...) ${filter_f} ${filter_F}

```

So if you want all reads then make these strings empty. That is far more generic than your idea because it can well be that the flag values you need do not exist.

---
Tags: samtools, sam
---


title: Different contrasts in DESeq2

text:
# Question
Title: Different contrasts in DESeq2
I have a treatment and control in two time points like this

```
 > design
                                 X condition time
    1       CTRL_24_hrs_replicate1   control   24
    2       CTRL_24_hrs_replicate2   control   24
    3       CTRL_24_hrs_replicate3   control   24
    4  treatment_24_hrs_replicate1        t   24
    5  treatment_24_hrs_replicate2        t   24
    6  treatment_24_hrs_replicate3        t   24
    7       CTRL_48_hrs_replicate1   control   48
    8       CTRL_48_hrs_replicate2   control   48
    9       CTRL_48_hrs_replicate3   control   48
    10 treatment_48_hrs_replicate1        t   48
    11 treatment_48_hrs_replicate2        t   48
    12 treatment_48_hrs_replicate3        t   48
    >

```

I want to test between between treatment and control considering time point 2 hours to 4 hours

I have done like this

```
dds <- DESeqDataSetFromMatrix(countData=a,colData=design, design=~time + condition + time:condition)

ddsTC <- DESeq(dds, test="LRT", reduced = ~ time + time:condition)

> resultsNames(ddsTC)
[1] "Intercept"                "condition_t_vs_control"
[3] "time"                     "conditiont.time"

```

What is the difference of the results of these codes

```
results(ddsTC, name="condition_t_vs_control", test="Wald")

```

**versus**

```
results(ddsTC, name="conditiont.time", test="Wald")

```

Actually I want to know what `condition_t_vs_control` gives and what `conditiont.time` gives

Thank you so much for any intuition

# Answer
`condition_t_vs_control` gives the effect of condition, `conditiont.time` is the interaction of condition and time. If you wanted to test the effect of time, use `name="time"`.
> 3 votes

---
Tags: rna-seq, differential-expression, deseq2
---


# Question
Title: Understanding the different designs in DESeq2
I am using DEseq2 and trying to understand the results obtained using different models.

I have a data design with 2 genotypes and 2 time points.

```
    sample genotype time
1   WT_S1       WT   T1
2   WT_S2       WT   T1
3   WT_S3       WT   T1
4   WT_S4       WT   T2
5   WT_S5       WT   T2
6   WT_S6       WT   T2
7   KO_S1       KO   T1
8   KO_S2       KO   T1
9   KO_S3       KO   T1
10  KO_S4       KO   T2
11  KO_S5       KO   T2
12  KO_S6       KO   T2

```

**I want to know the differences in results obtained from genotype coefficient and time coefficient when using different models.**

Model1) ~ genotype

Since it is comparing the differences in genotype regardless of time , samples 1-6 vs 7-12 are being compared.

Model2) ~ time

This model compares time T2 vs T1 regardless of genotypes. So it is comparing samples (1-3 + 7-9) vs (4-6 +10-12), **is this correct?**

Model3) ~ genotype + time

My understanding is that this model assumes the genotype effect is the same at both time points and so it adds a time effect to both genotypes.

**Does that mean results(obj3, name="genotype\_KO\_vs\_WT") give differences in genotypes by comparing samples 4-6 vs 10-12 ?**

**What samples are being compared in the results obtained from results(obj3, name="time\_T2\_vs\_T1")? How is it different from model 2?**

Model 4) ~ genotype + time + genotype:time

Here I understand results(obj4, name="genotype\_KO\_vs\_WT") gives the differences in genotypes at reference levels ie, samples 1-3 vs 7-9

**results(obj4, name="time\_T2\_vs\_T1"). What samples are compared here? How is it different from ‘time\_T2\_vs\_T1’ results in model 2 or model3?**

The interaction term as I understand is giving the specific effect due to KO at time T2 controlling for the baseline differences in genotypes. results(obj5, name="genotypeKO.timeT2")

**What samples are being compared from the model matrix?**

# Answer
> 0 votes
I strongly suggest that you take your normalized counts, get the averages of different groups, and compare those ratios to the different designs.

That said, you might not be able to replicate the fold changes from ~ genotype + time well in Excel. What that design does is basically say "Lots of the variance within each genotype/time point is due to the time point/genotype. Model that so that those differences are taken into account." Whereas a design of just ~ time, the software will go "Oh, there's a lot of variance within each time point". You won't be able to replicate that well in Excel...but that's why you are using DESeq.

That design is probably better than just time or genotype alone, if you want to look at all the KO versus all the WT

> Does that mean results(obj3, name="genotype\_KO\_vs\_WT") give differences in genotypes by comparing samples 4-6 vs 10-12 ?

No. The simplest way to do that comparison, and comparisons like it is to make a third data column that's genotype\_time, and use contrasts to specify what subgroup to compare to what subgroup.

> Here I understand results(obj4, name="genotype\_KO\_vs\_WT") gives the differences in genotypes at reference levels ie, samples 1-3 vs 7-9

Yes. But since that's super confusing, I'd recommend against it. The way I outlined above is far more readable, and should work the same.

The interaction design is best used when you want to see of the changes caused by time are different in the KO than the WT. In theory, you could just generate the two timepoint fold changes, and subtract one from the other, but you'd have no p-value. Using a design with interactions will get you that fold change difference with a p-value.

---
Tags: r, statistics, differential-expression, deseq2
---


# Question
Title: Different structures of same protein from two softwares
I used AlphaFold and Swiss-Model to get the structure of a protein (KiSS1R) and then compare them. The structures don't match and I was wondering why the softwares produce models that are different in some regions?

Thankful for all help!

# Answer
> 2 votes
## Superposition

> https://michelanglo.sgc.ox.ac.uk/data/d6e427bb-a0cb-4a65-b988-b08237f5054a

## SwissModel

> https://swissmodel.expasy.org/repository/uniprot/Q969F8

SwissModel models are threaded against a template. In this case the homologue has 31% identity, which isn't great. The template is 6TP3 which has a rather high R\_free of 0.32 and it's solved at 3 Å, but it's a membrane structure so I personally would forgive it —others may not. SwissModel has switched primary metric to its QMEANDisCo, which is present as global (0.6 in this case) and per residue (saved in the PDB file as the temperature factor column), which is roughly similar conceptually to AlphaFold2's pLDDT (vide infra), i.e. it tells you what resides to trust. 0.6 not great, not terrible. A nice thing in SwissModel is that, by virtue of being threaded, it's easy to steal ligands and other protein and it shows models of oligomers.

## AlphaFold2

> https://alphafold.ebi.ac.uk/entry/Q969F8

AlphaFold2 is by far the best modelling algorithm as see in the latest CASP competition (#14). However it natively doesn't do oligomers —ColabFold notebook can.

AlphaFold2 gives a per residue score (plDDT). AlphaFold2 models the whole sequence, which results in some nasty "spaghetti loops": these have really bad scores, so need to be ignored. These are not trash, but are likely unstructured in isolation, but bind other protein etc. (more info in this blog post). This structure has N- and C-terminus loops that need to be ignored. The N in particular is the signal that would be cleaved off —to get the membrane embedding dots, one can use the OPM server.

## Choice

The difference in the two models only involves low certainly regions and the support for the SwissModel is problematic. So the AlphaFold2 model is better, but requires removal of the C- and N-terminal loops.

---
Tags: proteins, database, protein-structure, software-quality
---


# Question
Title: Convert Entrez gene ids and HUGO symbols, whose genome assembly is unknown, to GRCh38 Ensembl gene ids
I want to use a file available on a database called `DepMap` (https://depmap.org/portal/download/). This file called `CRISPR_gene_effect.csv` contains gene ids in the format `HUGO (Entrez)` eg. `TP53 (7157)`, without having any complementary information regarding which genome assembly (GRCh37/38) the Entrez ids (and also HUGO symbols) belong to.

I want to convert these gene ids in the `HUGO (Entrez)` format to GRCh38 Ensembl ids (specifically to Ensembl release 100), so that I can map the data to Ensembl's annotations. I would like to get suggestions on this.

I feel that figuring out the genome assembly version (GRCh37/38) of the Entrez ids would be critical in this conversion. Although, I am not aware if Entrez ids change from one genome assembly to other.

# Answer
> 2 votes
Hugo names don't really have an associated assembly. They are just the official name of the gene in human. You will find them in all assemblies (except if they hadn't been annotated/identified in one of the older ones). Entrez entries are also not tied to a specific genome assembly. They usually have mappings to multiple ones. For example, Hugo gene name A1BG, with Entrez ID 1 has:

Now, if you want to get Ensembl gene IDs for a list of Hugo gene names, you can use BioMart to convert. First, get the list of target IDs. You can use either Hugo or Entrez. I found the Entrez simpler to parse out, so I went for those (note that this assumes GNU `grep`):

```
grep -oP '(?<=\()\d+(?=\))' CRISPR_gene_effect.csv > entrez.ids

```

I then went to BioMart and selected "Ensembl genes 104" (note that I could only get Ensembl version 104, I assume you're not adamant about getting version 100 given that it doesn't seem to be current) and "Human genes GRCh38.p13" :

Next, I set the filters to "NCBI gene (formerly Entrezgene) ID(s)" and uploaded my file:

And, finally, selected the Attributes I want to see, namely Ensembl gene ID, HGNC (Hugo) symbol and Entrez ID so I can map everything to each other:

I then clicked on "Results" and got:

\[

You can now download this to a file and have your mapping. You might even be able to use the one I created: https://www.ensembl.org/biomart/martview/352aae8e0b5c774cbbee73d67634d577.

---
Tags: ensembl, identifiers, entrez
---


# Question
Title: Shotgun analysis: Where and how to start?
I am a complete rookie in shotgun analysis, having only performed marker gene analysis before. I have a pipeline that I want to try but I can't wrap my head around how it all works.

Please correct me if I'm wrong:

1. I start with a paired-end Illumina short reads fastq files (.R1,.R2)
2. I will use FASTQC to trim according to quality. I assume that since data has been demultiplexed, the primers and and adaptors are removed as well.

? First question, a lot of people advise to merge all the samples together and then do the co-assembly since apparently it helps with the assembly. Would this be necessary if I also have nanopore reads and only 9 samples? If if I did this, how do I then separate the samples from each other?

3. Because I also have nanopore, I will use a hybrid assembler to get contigs. The output will be a file with contigs.
4. I then plan to use a binning tool (Metawrap) to create bins.

? I believe the output of these will be a file containing bins per sample. How do I estimate the abundance of each bin in the sample? And how do I factor in the depth of sequencing to normalise them? This is my biggest problem, I can't wrap my head around it.

?Can some explain VERY plainly what is sequencing depth and how do I factor it into my analysis? I've seen some people using it for normalisation of the bin/contigs abundance but again, this is idea is not intuitive to me at all.

5. From then on, taxonomy and functional annotation is pretty straightforward.

#Update

So a little background on the project. I have samples sequenced from a bioreactor containing a small community of bacteria (no host removal required). The samples derive from 6 different reactors = 6 samples as well as one reactor will be samples three times (time series) so I have 9 samples in total. The shotgun Illumina will be combined with nanopore. The idea is to get the composition of the community (what kind of bacteria are there?), abundance (what is their relative quantity?) but also functions (what do they do?).

I'm sorry for a very generic question, I would really appreciate somebody answering very plainly without too much jargon. Hopefully it will help me and other people who are just starting out with the shotgun.

Thank you!

# Answer
**Question 1**

For co-assembly vs individual assembly, it is sometimes advisable to do them separately, perhaps if you expect samples to be very different or if you have huge amounts of data. For \<10 species I imagine that you will have quite high coverage, especially with nanopore.

Some metagenomic binners (e.g. MetaBAT2, a standard from a few years ago) can even incorporate information across samples, which might argue for maintaining at least some separation.

"How do I separate the samples from each other" I guess the point is that you **don't**, for assembly at least. You can then go in post hoc with your sequencing data from each sample and try to figure out with abundance analysis which of your genomes are in each sample.

**Question 2**

For abundance analysis/depth, there are strong statistical parallels in RNA-seq analysis, as Lior Pachter is fond of saying, if that is helpful. Any number of metagenomic shotgun abundance computations can be used to estimate abundance of genomes, including Pachter's metakallisto (linked therein). I found metakallisto fairly straightforward to use on contigs resulting from assembly.

Basically, you have some set of DNAs in your sample with some stoichiometry (relative abundance with respect to each other). The idea is to use DNA sequencing to estimate that stoichiometry (absolute abundance is much harder so most people don't bother).

Here is a broad view of the workflow:

1. You come up with some estimate of the DNAs present in the community. This can be a shotgun assembly, or some set of reference genomes, or whatever.
2. You count how many sequencing reads derive from each of those DNAs, usually via alignment or pseudoalignment.
3. YOu do some numerical corrections for systematic covariates such as DNA length. (These can look like RPKMs or TPMs from RNA-seq; metakallisto literally uses TPM calculations)
4. The resulting counts should be abundances, more or less.

You don't have to figure this stuff out yourself. There are various end to end pipelines that do this.

You *will* almost certainly need to compute abundances of contigs if you want to do metagenomic binning.
> 1 votes

---
Tags: ngs, metagenome
---


# Question
Title: Is there any genetic sequence, such that, for all the kinds of mutation, it still has a 100% probability to be passed on?
This question is more on the theoretical-encoding side. Given the known kinds of genetic mutation, transposition, deletion, and such, what can be said about genes that would be left untouched by the mutation, left invariant by the permutation, if you will? Are there any such sequences? Would they be more common in a population, if not detrimental?

# Answer
Not in the way you suggest: a sequence of letters with the property of being immutable by any type of permutation on its constituent bases. I would be interested if anything of the likes exists in mathematics.

Having said that, there are examples of sequences that seem to have such property but only in the light of evolution and the biological processes that come with it:

Ultra-conserved regions have been described as sequences of about 200 base pairs that appear to be exactly the same from humans to mice. They have not experience change in about ~300 million years of evolution (!) and exposure to the various forces of mutation.

The reason for this is not an underlying 'invariant to permutation' property, but rather strong negative selection that prevents individuals that acquire changes on the sequence from being able to pass the changes over to their progeny, thus effectively maintaining only the original version in the population.

The exact mechanisms, functions and properties of these regions are under research. You might want to see if the invariant properties you have in mind have anything to do with these interesting regions.
> 1 votes

---
Tags: genetics, molecular-genetics, phylogenetics
---


# Question
Title: Why does it require 4 vectors to calculate a torsion angle?
The BioPython documentation shows the following:

> 11.6.3 Measuring torsion angles Use the vector representation of the atomic coordinates, and the calc dihedral function from the Vector module:
>
> ```
> >>> vector1 = atom1.get_vector()
> >>> vector2 = atom2.get_vector()
> >>> vector3 = atom3.get_vector()
> >>> vector4 = atom4.get_vector()
> >>> angle = calc_dihedral(vector1, vector2, vector3, vector4)
>
> ```

We know that a torsion angle between N-C-alpha and C-alpha-C planes.

So, the minimum number of points required to form a plane is 3.

Then, Why does BioPython require 4 vectors?

# Answer
> 5 votes
From wikipedia on torsion angles:

> A dihedral angle is the angle between two intersecting planes or half-planes. In chemistry, it is the clockwise angle between half-planes through two sets of three atoms, having two atoms in common.

Two sets of three with two in common makes four total atoms.

You are defining not one plane but two half-planes or planes.

Wikipedia includes the following as an illustrative image:

"Angle between two half-planes (α, β, green) in a third plane (pink) which cuts the line of intersection at right angles"

In this image I interpret the atoms to be at the vertices of the planes, e.g. the two atoms that are in common between the two planes are those at the points of intersection of the planes, and the two other atoms extend the planes $\alpha$ and $\beta$.

# Answer
> 4 votes
To compute dihedral angles you require 4 atoms and hence for vectors to determine two planes (see this figure here) whose angle you want to measure - dihedral angle.

You can also check Fig 55 from the GROMACS documentation

https://manual.gromacs.org/documentation/current/reference-manual/analysis/bond-angle-dihedral.html

Definition of (virtual) Cα backbone bond (κ) and torsion (τ) angles. from https://journals.plos.org/plosone/article/figure?id=10.1371/journal.pone.0215141.g002

---
Tags: python
---


# Question
Title: How to convert species names into common names?
I’m trying to find common names from a list of scientific names (not all will have them though).

I was attempting to use `taxize` in R but it aborts if it doesn’t find an entry in EOL and I don’t know a way around this other than manually editing the list - in which case I might as well check all this way.

Does anyone have a better way of doing this?

# Answer
> 15 votes
As Pierre mentioned, NCBI is a good resource for this kind of transformation.

You can still use `taxize` to perform the conversion:

```
library("taxize")

species <- c('Helianthus annuus', 'Mycobacterium bovis', 'Rattus rattus', 'XX', 'Mus musculus')

uids <- get_uid(species)

# keep only uids which you have in the database
uids.found <- as.uid(uids[!is.na(uids)])
# keep only species names  corresponding to your ids
species.found <- species[!is.na(uids)]

common.names <- sci2comm(uids.found, db = 'ncbi')
names(common.names) <- species.found

common.names

## output:
## $`Helianthus annuus`
## [1] "common sunflower"
##
## $`Mycobacterium bovis`
## character(0)
##
## $`Rattus rattus`
## [1] "black rat"
##
## $`Mus musculus`
## [1] "house mouse"

```

This works correctly if no common name is available, or if there is no species in the database. You can also try a different database, e.g. `db='itis'`.

The advantage of this approach is that you do not have to parse a file yourself. The downside is that it can be slow for large list, since for every species a database request is performed.

# Answer
> 14 votes
use the NCBI taxon dump under ftp://ftp.ncbi.nih.gov/pub/taxonomy

in taxdmp.zip you'll find all the names for a given NCBI taxon

```
$ grep -w ^9606 names.dmp
9606    |   Homo sapiens    |       |   scientific name |
9606    |   Homo sapiens Linnaeus, 1758 |       |   authority   |
9606    |   human   |       |   genbank common name |
9606    |   man |       |   common name |

```

# Answer
> 4 votes
Getting back to the original question about using `taxize` to search specifically in EOL, it looks like the most recent version of the package does not produce this error.

The current version of R I'm using (3.4.2) produces the following output:

> sci2comm( c("Cancer", "Lophopanopeus"), db="eol" )
>
> $Cancer \[1\] "En art taskekrabbe" "Felsenkrabbe" "Atlantic rock crab"
> \[4\] "Jaiba de roca amarilla" "Tourteau poïnclos" "Granciporro atlantico giallo"
> \[7\] "Atlantische rotskrab" "Sapateira de rocha do Atlântico" "Bogenkrabbe"
> \[10\] "red rock crab" "Pacific rock crab" "Jaiba del Pacífico"
> \[13\] "Tourteau du Pacifique" "Granciporro del Pacifico" "Pacifische rotskrab"
> \[16\] "Sapateira de rocha do Pacífico" "Jonahkrabbe" "Jonah crab"
> \[19\] "Jaiba de roca Jonás" "Tourteau jona" "Grancipporo atlantico rosso"
> \[22\] "Jonakrab" "Sapateira boreal" "Bou"
> \[25\] "Almindelig taskekrabbe" "Taskekrabbe" "Taschenkrebs"
> \[28\] "Edible crab" "rock crab" "Buey"
> \[31\] "Buey de mar" "Jaiba de roca masera" "Masera"
> \[34\] "Paguro" "Pato" "Sabago"
> \[37\] "Isotaskurapu" "Crabe dormeur" "Dormeur"
> \[40\] "Tourteau" "crabe-dormeur" "Granchio di mare"
> \[43\] "Granciporro" "Noordzeekrab" "grote krab"
> \[46\] "grote noordzeekrab" "hoofdkrab" "slagkrab"
> \[49\] "steenkrab" "toogkrab" "zeekrab"
> \[52\] "Sapateira"
>
> $Lophopanopeus \[1\] "blackclaw crestleg crab" "knobkneed crestleg crab" "knobknee crestleg crab" "molarless crestleg crab"
>
> $'Pugettia foliata' character(0)

A slightly older version of R (3.3.2) produces an error when one of the taxa was not found:

> sci2comm( c("Cancer", "Lophopanopeus", "Pugettia foliata"), db="eol" )
>
> Error in eol\_pages(taxonconceptID = x, common\_names = TRUE, ...) : Not Found (HTTP 404).

I am very thankful for the package updates!

# Answer
> 1 votes
As Pierre suggested, you can get a dump of such names from NCBI. Then, to query it for the common name for a species using its Latin name, you can do:

```
$ awk -F'\t' -vname="Mus musculus" '($7=="scientific name" && $3==name){
                                      a[$1]=$3
                                    }
                                    ($1 in a && /genbank common name/){
                                      print $3
                                    }' names.dmp
house mouse

```

# Answer
> 0 votes
Here is a code on github that can generate a csv file as well as a python dictionary. https://github.com/rahulnutron/scientific\_name\_to\_common-name\_converter

---
Tags: r, database, taxonomy
---


title: general cheminformatics references

text:
# Question
Title: general cheminformatics references
Can anyone explain what is cheminformatics and how does it help in the future?

I need to write an article that should be understandable for 14-20 yr age group. Please help with references, I don't know anything.

# Answer
> 3 votes

---
Tags: computational-biochemistry, cheminformatics
---


# Question
Title: Long-Covid (LC) GWAS project proposal: should I include a non-LC GWAS?
I am working on a project proposal for an exam, it involves a budget of 1M euros and I got to design a genomic study. I chose a GWAS to assess altered gene expression after covid19 that leads to long-COVID. I thought about 3 cohorts:

Case cohort1: Long-Covid patients (N: 3,000)

Case cohort2: Fully recovered patients (N: 3000)

Control cohort: healthy patients. (N: 6000)

Would it make sense to run two GWAS for the two case cohorts against the control cohort (cohort1/controls and cohort2/controls) then do a miami plot to see if there are some hits that differs from the two cohorts, and pick the ones only presented in the Long-covid for further assessment?

The study will also be supported by a validation step for both gwas. What do you think about this idea? Would it make sense or the fully-recovered patients against control cohorts is totally useless? Thank you in advance

# Answer
> 2 votes
In the end I manage to do 2 cohorts 1: Long Covid Patients (1300N) 2: NON Long Covid Patients that had Covid (1300).

It was the most rational thing to do, and professor liked it. So that's how I manage in th end.

---
Tags: covid-19, gwas
---


# Question
Title: How can I calculate the distance w.r.t. a beta sheet?
Suppose, I have a protein (say, 16PK).

and, I want to calculate

1. the distance between one of its alpha-carbons and the nearest beta-sheet.
2. the distance between two strands in a beta-sheet (i, j).

How should I proceed to achieve this?

# Answer
> 1 votes
The documentation for the PDB module in biopython, which you have tagged, is very clear about how distances are measured between two given Cα atoms.

To find out the secondary structure and interactions between strands of a sheet, you may want to look at the **header** of the PDB. This is the stuff before the `ATOM`/`HETATM` entries. **`SHEET` entries** are what you are after. Please read the PDB format specifications for secondary structure. Biopython PDB does not parse headers well, but reading the file and finding the SHEET entry and positionally splicing out the strings of interest from the entry is beyond straightforward.

Example take from the above link:

```
         1         2         3         4         5         6         7         8
12345678901234567890123456789012345678901234567890123456789012345678901234567890
SHEET    1   A 5 THR A 107  ARG A 110  0
SHEET    2   A 5 ILE A  96  THR A  99 -1  N  LYS A  98   O  THR A 107
SHEET    3   A 5 ARG A  87  SER A  91 -1  N  LEU A  89   O  TYR A  97
SHEET    4   A 5 TRP A  71  ASP A  75 -1  N  ALA A  74   O  ILE A  88
SHEET    5   A 5 GLY A  52  PHE A  56 -1  N  PHE A  56   O  TRP A  71
SHEET    1   B 5 THR B 107  ARG B 110  0
SHEET    2   B 5 ILE B  96  THR B  99 -1  N  LYS B  98   O  THR B 107
SHEET    3   B 5 ARG B  87  SER B  91 -1  N  LEU B  89   O  TYR B  97
SHEET    4   B 5 TRP B  71  ASP B  75 -1  N  ALA B  74   O  ILE B  88
SHEET    5   B 5 GLY B  52  ILE B  55 -1  N  ASP B  54   O  GLU B  73

```

---
Tags: python, proteins, protein-structure
---


# Question
Title: curl: (22) The requested URL returned error: 400 Bad Request when using entrez direct
I've just dowloaded entrez direct on my Mac and I was trying to perform some simple queries like:

```
esearch -db pubmed -query "tn3 transposition immunity" | efetch -format medline

```

but I keep getting this error

```
curl: (22) The requested URL returned error: 400 Bad Request
 ERROR:  curl command failed ( Mar 19 Ott 2021 13:02:26 CEST ) with: 22
https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi -d db=pubmed&id=Unable%2Cto%2Clocate%2Cxtract%2Cexecutable.%2CPlease%2Cexecute%2Cthe%2Cfollowing%2Cnquire%2Cdwn%2Cftp.ncbi.nlm.nih.gov%2Centrez%2Centrezdirect%2Cxtract.Darwin.gz%2Cgunzip%2Cf%2Cxtract.Darwin.gz%2Cchmod%2Cx%2Cxtract.Darwin&rettype=medline&retmode=text&tool=edirect&edirect=16.1&edirect_os=Darwin&email=mariocatalano%40MBP-di-Mario.lan"

```

I'm using bash scripts on my Mac's terminal

How should I fix this?

# Answer
It's not a curl problem. (The 22 error is just passing along the server-side HTTP 400 error, which is because it dumped an error message into a URL.) I think terdon is right and there's something wrong with your edirect installation. You have esearch and efetch on your PATH but apparently not xtract and/or xtract.Darwin.

---

More detail:

Look toward the end of the error message, at the text being passed to curl's `-d` argument:

```
db=pubmed&id=Unable%2Cto%2Clocate%2Cxtract%2Cexecutable.%2CPlease%2Cexecute%2Cthe%2Cfollowing%2Cnquire%2Cdwn%2Cftp.ncbi.nlm.nih.gov%2Centrez%2Centrezdirect%2Cxtract.Darwin.gz%2Cgunzip%2Cf%2Cxtract.Darwin.gz%2Cchmod%2Cx%2Cxtract.Darwin&rettype=medline&retmode=text&tool=edirect&edirect=16.1&edirect_os=Darwin&email=mariocatalano%40MBP-di-Mario.lan"

```

It's putting this in the `id=` part:

```
Unable%2Cto%2Clocate%2Cxtract%2Cexecutable.%2CPlease%2Cexecute%2Cthe%2Cfollowing%2Cnquire%2Cdwn%2Cftp.ncbi.nlm.nih.gov%2Centrez%2Centrezdirect%2Cxtract.Darwin.gz%2Cgunzip%2Cf%2Cxtract.Darwin.gz%2Cchmod%2Cx%2Cxtract.Darwin

```

which is maybe something like (trying to guess which characters got turned into `%2C`):

> Unable to locate xtract executable. Please execute the following: nquire -dwn ftp.ncbi.nlm.nih.gov entrez/entrezdirect xtract.Darwin.gz; gunzip -f xtract.Darwin.gz; chmod +x xtract.Darwin

But you should already have xtract and xtract.Darwin from the same edirect install that gave you esearch and efetch.

For what it's worth the official Entrez Direct install on Mac OS looks a little goofy. When I try it myself it ends with:

```
Entrez Direct has been successfully downloaded and installed.

-e In order to complete the configuration process, please execute the following:

  echo "export PATH=\\\${PATH}:" >> $HOME/.bashrc
  echo "export PATH=\\\${PATH}:" >> $HOME/.zshrc

or manually edit the PATH variable assignments in your .bashrc and .zshrc files.

Would you like to do that automatically now? [y/N]

```

Those commands don't look right (what's getting added to the path? an empty string?) and if I say y, I get this at the end of .bashrc and .zshrc which definitely isn't right:

```
export PATH=${PATH}:

```

It seems like an empty string at the end of PATH behaves like `.` (current directory), so I can go into ~/edirect and run your command and it *does* actually work. So this isn't exactly your problem... but still, that seems wrong.
> 2 votes

---
Tags: bash, entrez, curl
---


# Question
Title: Convert R RNA-seq data object to a Python object
I have done some work in R and would like to try a Python tool.

What is a good way to **import** the data (and its annotations etc) **as a Python object?**

I am particularly interested in converting a `Seurat` object into an AnnData object. (Either directly or as a Python object which can be converted into an `AnnData`.)

# Answer
> 10 votes
A simple solution for converting Seurat objects into AnnData, as described in this vignette:

```
library(reticulate)
seuratobject_ad <- Convert(from=seuratobject, to="anndata", filename="seuratobject.h5ad")

```

---

Alternatively, one can use Loom, *"a file format designed to store and work with single-cell RNA-seq data".*

In R, install LoomR:

```
devtools::install_github(repo = "mojaveazure/loomR")

```

Convert from Seurat object to loom:

```
pfile <- Convert(from = pbmc_small, to = "loom", filename = "pbmc_small.loom")
pfile$close() # close loom objects when finished using them.

```

Then import loom object in Python using loompy, or directly as AnnData:

`scanpy.api.read_loom`

---

Alternatively, see `feather`.

---

Or export as some text format (csv, json) then import into Python.

# Answer
> 0 votes
The accepted solution is probably the best for older objects of type `seurat` created with Seurat package v2. For newer `Seurat Objects`, there is a new tool designed specifically for this purpose, called `SeuratDisk`.

There is a nicely documented vignette about the Seurat \<-\> AnnData conversion. In short:

1. In R, save the Seurat object as an h5Seurat file:

```
library(Seurat)
library(SeuratDisk)

SaveH5Seurat(SeuratObject, filename = "file.h5Seurat")
Convert("file.h5Seurat", dest = "h5ad")

```
2. In python, simply load the file:

```
import scanpy

adata = scanpy.read_h5ad("file.h5Seurat")

```

# Answer
> 0 votes
Since both Seurat and Scanpy are evolving rapidly, the safest way is to save your Seurat object into MatrixMarket format + metadata csv, and import them into Scanpy.

---
Tags: rna-seq, r, python
---


# Question
Title: Given a multiple sequence alignment, how do I output all of the non-consensus characters by location?
I have an MSA from MAFFT (as a FASTA file; there are several sequences in the MSA), and I would like to quantify the "variants" between these sequences.

What is the standard tool to do something like this?

I would prefer to know the non-consensus variants in the the entire MSA, but that may be too tricky, so I would like to figure out the non-consensus variants between those select sequences.

I'm familiar with calculating an MSA, and then visualizing it---but I would like an output file with all of the different bases.

# Answer
> 2 votes
A very basic solution in python:

```
from Bio import SeqIO
seq_dict = {i : rec.seq for i,rec in enumerate(SeqIO.parse("your_msa.aln", "fasta"))}
for pos in range(len(seq_dict[0])):
    pos_bases = ([seqio[i] for seqio in [*seq_dict.values()]])
    # If position contains more than one unique base
    if len([*set(pos_bases)])>1:
        print(pos, pos_bases)        # prints the position and alternative bases there

```

If you need something more useful, just leave a comment here and I'll develop this further.

---
Tags: sequence-alignment
---


# Question
Title: Creating Experimental Components on SBOL
I am developing software to store experiments or studies using SBOL3. In the regular synthetic biology workflow we will find the experimental components Strain and Media. A safe way to implement them is using `tyto.SBO.functional_entity` as type for both.

Other attempts to implement media (https://github.com/Bioprotocols/paml/blob/main/examples/growth\_curve.py)looks like this:

```
SC_media = sbol3.Component('SC_Media', 'TBD', name='Synthetic Complete Media')

```

Therefore there is a need for standardise these types on SBOL3 to provide more accurate representations. The Experimental Factor Ontology (EFO) have ontology terms for Strain http://www.ebi.ac.uk/efo/EFO\_0005135 and Media http://www.ebi.ac.uk/efo/EFO\_0000579 that can be used as type. Then draft SBOL components looks as follow:

```
strain = sbol3.Component(fj_strain.name(), http://www.ebi.ac.uk/efo/EFO_0005135)
strain.name =fj_strain.name
strain.description = fj_strain.description()

media = sbol3.Component(fj_media.name(), http://www.ebi.ac.uk/efo/EFO_0000579)
media.name =fj_media.name
media.description = fj_media.description()

```

The problem with this is the lack of readability on the code given the obscure ontology URIs. This problem has been solved for other ontologies using tyto, but EFO is not supported out-of-the-box and the configuration of new ontologies is not working. On tyto repository there is an example:

```
from tyto import EBIOntologyLookupService, Ontology
KISAO = Ontology(uri='http://www.biomodels.net/kisao/KISAO_FULL#', endpoints=[EBIOntologyLookupService])
KISAO.Gillespie_direct_algorithm
'http://www.biomodels.net/kisao/KISAO#KISAO_0000029'

```

But I copied and runed that code but it raise this error: \[ERROR\] tyto.py 70: Ontology http://www.biomodels.net/kisao/KISAO\_FULL# is not available at EBI Ontology Lookup Service

What I would like to do is something like this:

```
from tyto import EBIOntologyLookupService, Ontology
efo = Ontology(uri='http://www.ebi.ac.uk/efo', endpoints=[EBIOntologyLookupService])
efo.media

```

But this raises the same error. How can I add EFO to tyto? What is the URI that should be provided? Should I raise an issue in tyto to ask for EFO out-of-the-box support?

Do you agree with the ontology terms for Strain and Media? Any feedback?

# Answer
Section 7.9 of the SBOL3 specification already provides recommended terms for experimental elements. Among these, Section 7.9.2 specifically addresses strain and media, reading in pat:

> The roles of reagents may vary in context: for example, arabinose may serve as an inducer or as a media carbon source. As such, contextual role SHOULD be indicated by an NCI Thesaurus (NCIT) term in a role property of the Feature. Examples include:
>
> * NCIT:C64356: PositiveControl
> * NCIT:C48694: Cell
> * NCIT:C85504: Media
> * NCIT:C14419: Strain
> * NCIT:C120268: Inducer

NCIT is one of the ontologies that is built-in for tyto:

```
tyto.NCIT.Strain
Out[4]: 'https://identifiers.org/ncit:C95501'
tyto.NCIT.Media
Out[5]: 'https://identifiers.org/ncit:C19442'

```
> 0 votes

---
Tags: sbol, ontology
---