# Question
Title: Gene names from components in NMF analysis
*This question was also asked on Biostars*

I am using Rcppml package in R for my NMF analysis, I have a matrix from single-cell analysis. I have cells scores and I wanted to know how to extract genes instead of the components.

i used the following codes.

```
t1.cpm<- LogNormalize(data = as.matrix(t1@assays$RNA@counts)) # Log converting the data form seurat obj

t1.model <- RcppML::nmf(as.matrix(t1.cpm), 30, verbose = F, seed = 1234)

t1.w <- t1.model$w
t1.d <- t1.model$d
t1.h <- t1.model$h

dim(t1.w)
rownames(t1.w)<- rownames(t1.cpm)
colnames(t1.w)<- paste0("component", 1:30)
rownames(t1.h)<- paste0("component", 1:30)

```

and my output is as follows

```
        T_1_TACCTATCAGATGGCA T_1_CCTATTATCAGGCAAG
component1             0.000000000            0.0000000000
component2             0.000000000            0.0023590277
component3             0.001512931            0.0000000000
component4             0.000000000            0.0021574559
component5             0.001175799            0.0020222144
component6             0.000000000            0.0004463346

```

It would be helpful if you can let me know how to extract Gene names instead of components.

Thanks, Dave.

---

Rcppml is a machine learning library.

https://github.com/zdebruine/RcppML

NMF is Non-negative matrix factorization which is used in bioinformatics and AI, which can examine the latent relationships in experimental data sets.

# Answer
> 3 votes
NMF of a gene expression matrix into $k$ components works to decompose the expression matrix into an activity matrix ($H$) and latent factor matrix ($W$).

$W$ describes the $k$ latent factors (often called components) by providing a score for how much each gene contributes to each component. $H$ describes how much each component contributes to the expression of each cell.

To extract the activity and component matrices from your NMF model fit using `RcppML`, use the following code:

```
# cells is a Seurat object
model <- RcppML::nmf(cells@assays$RNA@data, k = 5, verbose = F, seed = 1234)
W <- model$w # Components Matrix
H <- model$h # Activity Matrix

# Name genes/components for latent factor matrix
colnames(W) <- paste0("component", 1:k)
rownames(W) <- rownames(cells) # Gene names

# Name components/cells for activity matrix
colnames(H) <- colnames(cells) # Cell IDs
row names(H) <- paste0("component", 1:k)

```

```
> head(W)
              component1   component2   component3 component4   component5
MIR1302-2HG 0.000000e+00 0.000000e+00 0.000000e+00          0 0.000000e+00
FAM138A     0.000000e+00 0.000000e+00 0.000000e+00          0 0.000000e+00
OR4F5       0.000000e+00 0.000000e+00 0.000000e+00          0 0.000000e+00
AL627309.1  2.927886e-07 9.248074e-07 7.961105e-07          0 2.802758e-06
AL627309.3  3.486284e-08 0.000000e+00 1.916799e-07          0 7.938089e-08
AL627309.2  0.000000e+00 0.000000e+00 0.000000e+00          0 0.000000e+00

> head(H)
                 cell_1       cell_2       cell_3       cell_4       cell_5
component1 0.0001370465 0.0002195795 0.0001628898 0.0002239623 0.0001484157
component2 0.0003604048 0.0002239471 0.0005688147 0.0005242678 0.0004821531
component3 0.0003948121 0.0001485325 0.0002885584 0.0002242272 0.0002865587
component4 0.0001391838 0.0002418495 0.0001361317 0.0001382889 0.0001558255
component5 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000

```

People will often use the top $N$ genes to describe a component. The following code gets the top $N$ genes for each component and stores them in a list

```
# This code assumes the above code block has been already run
require(dplyr) #used for %>% piping
N <- 100
top_genes <- list()
for (i in 1:dim(W)[2]) {
  name <- colnames(W)[i]
  top_genes[[name]] <- sort(W[, i], decreasing = T) %>% head(N) %>% names()
}

```

---
Tags: scrnaseq, machine-learning, nmf
---


# Question
Title: How can I get the duration of the events in Nanopore?
PacBio has the concept of IPD (InterPulse Duration) which is the time between two detected consecutive sequences in the raw signal. I have been trying to extract this value in the Nanopore, but I cannot find them in the .fast5 files. But this info must be there:

This is from the paper Nanopore sequencing technology, bioinformatics, and applications:

''Each current segment contains multiple measurements, and the corresponding mean, variance, and duration of the current measurements together make up the ‘event’ data.'' Directly quoted from the paper, which suggests we should have the information about the duration in the raw signal.

How can I get the duration, I appreciate it if can anybody help?

Some info I read was confusing to me (e.g. http://simpsonlab.github.io/2015/04/08/eventalign/ which is nevertheless a great source for understanding what is going on), I knew single measurements' speeds are controlled, however, I didn't know the block stride is fixed and it made me so confused at first.

What I found related to my own question: I guess I was not completely correct in the terminology I was using, which might cause confusion. In nanopore sequencing, *translocation duration* refers to the time a fragment of DNA/RNA passes through the nanopore. Also for better visualization, I made this plot after running *albacore* basecaller:

as @gringer has explained, the *block stride* value is fixed (=5). The figure shows a small portion of albacore Event Table: The x-axis corresponds to the relative time from the start of measurements and the y-axis shows the current level. Each black dot shows a single measurement from an ONT device; the dots are segmented by a set of block strides (=5) used by albacore to call an event. The red lines show the average of the measurements in one segment. The shaded area in light gray represents the event that resulted in calling a base (one to one). Dark gray areas show events that resulted in calling two bases (split events). The white area corresponds to events with no called bases (skipped events). The boxes below the x-axis correspond to bases called in shaded areas and the bold base is the base that is called for that event.

# Answer
Unfortunately ONT no longer produce event data as default in their FAST5 output with the most recent basecallers, and they've been quite vague about how to generate it.

The easiest way to fix that is to use an older basecaller. If you tell Guppy v5.1.15 to output FAST5 files, it should include an event table that includes read start time and move information. This might also work with new Guppy, but ONT has told us that the event tables have gone away.

Processing the move tables to get information about signals per base is a bit complicated. Just using the metadata provided in the FAST5 files it's actually impossible, but you can get close enough to match what the ONT basecallers were trained on (which is probably close enough). To get closer than that, you need to process the raw signal yourself (e.g. see here).

Alternatively, you might be able to get that information out by following the tomboy documentation, which produces output that demonstrates signal-level changes:

https://github.com/nanoporetech/tombo

ONT have had various methods for basecalling over the years, and these methods change how the signal is chunked and processed.

When nanopore sequencing was first commercially available, base change points *were* resolved down to a single-sample level (as described in Jared Simpson's post which you linked). Or, at least, as accurately as they could at the time. Even with that approach, events were occasionally misclassified such that two consecutive events could represent the same base, or two consecutive events could involve a change of more than one base. To account for this, ONT created the *move table*, which described how many bases were shifted between one base and the next base. Initially ONT were using pentamers (i.e. 5 base sequences) for the output of their calling models (this had a physical reason, because 5 bases was about as many as would contribute to the instantaneous current), so the move table could have values of 0 \[indicating two adjacent events described the same base\], 1, 2, 3, 4, or 5. Although 5 was the maximum, because it was a pentamer model, what it actually meant was "there has been a large skip; we have no idea how large."

Unfortunately, sample-level signal parsing is no longer the case, which makes precise event alignment tricky. Most recently, individual signal samples are aggregated into what are called *strides*, which describe the smallest length of time that the basecaller accepts as representing a single base transition.

For this explanation, I'm going to assume that the directory location in the FAST5 file starts from `/baseLoc/`, with the FASTQ sequence found in `/baseLoc/Fastq`. This explanation may become quickly out of date, as ONT change their data and methods frequently.

Raw signal information is found at `/baseLoc/../Raw/Signal`. The stride length information is found in the `block_stride` attribute of `/baseLoc/Summary/basecall_1d_template`. This indicates the number of signal samples allocated for each stride. The move table is found at `/baseLoc/Move`, and represents the number of base shifts in each stride block (usually either 0 or 1). ***It sounds like that move table is the statistic you're interested in***; multiply the number of 0s before a 1 by the number of samples per stride (plus 1 for the base transition), and that's the base duration in samples, which should be divided by the sampling rate to get the number of seconds per base. Sampling rate can be found in the `sampling_rate` attribute of `/???/channel_id` (sorry, don't have a precise location, because I usually just assume 4 kHz... which will probably be changing shortly, from the sound of it).

Some of the signal at the start of the read may have been clipped, and it's necessary to take this into account in order to properly align the move table to the signal. The signal offset is found in the attribute `first_sample_template` for `/baseLoc/../../Segmentation_000/Summary/segmentation`, and represents the number of *samples* (not strides) that should be clipped off the start of the raw signal. There's one more thing that may help for aligning bases to signal, and that's the *base* offset relative to the move table. The basecaller works in chunks of multiple bases (kmers), so the base corresponding to the move table is actually somewhere in the middle of the called/modeled kmer. ONT doesn't provide any insight into this (it's possible they don't know), so I just try some offsets between 1-5 bases and pick the one that looks the best.

I have a FAST5 annotation R script that may help a little bit more in understanding how the alignment works, see here, specifically the code starting from the comment `## Load data from called fast5 file`.

Bear in mind that "number of samples per base" varies in a non-smooth fashion with quite a large dynamic range within an individual read sequence, and could be anything from 1 to over 50 (although I stop believing the basecaller if it reports anything over about 20). It doesn't make much sense calculating it as an aggregate statistic for each sequence, except possibly for detecting really poor basecalling (which is more easily picked up by looking at Q scores in the FASTQ file).
> 4 votes

---
Tags: nanopore
---


title: Convert Mascot mztab file format to XML

text:
# Question
Title: Convert Mascot mztab file format to XML
Cross-posted here

I would like to convert MASCOT mztab file to XML format. I have been looking for converters and I came across PRIDE Converter but this is no longer available for download from here https://github.com/PRIDE-Toolsuite/pride-converter-2.

Could someone please suggest alternate converters, if available?

# Answer
> 1 votes
There are a number of converters to be found at http://www.pastelbioscience.co.uk/resources/databases.html , although not sure that the one you are looking for is among them. May also be better to ask this question on twitter. There is a large #proteomics community that are very helpful and usually come up with an answer :-)

---
Tags: proteomics
---


title: IBM SPSS- 2-years survival of cohort

text:
# Question
Title: IBM SPSS- 2-years survival of cohort
Over the past few months I have been retrospectively analyzing patient files for my master's thesis in medicine. As the median overall survival of my cohort is low, I have decided that I want to analyze the 2-years survival of my population. However, I have not been able to find any literature on how to perform this analysis on SPSS.

Does anyone here have any experience analyzing the median survival over a specific period of time (i.e., 1-year, 2-years or 5-years) using IBM SPSS?

Thank you in advance!

Leon

# Answer
What you want is Cox regression (proportional hazard) models to perform **survival analysis**. It is years since I have used SPSS and would encourage consideration for the use of R (Python is a bit more complicated). There is a commercial go-between which is easier to use than R but more computationally efficient than SPSS, but I cannot remember its name at all.

There is a weblink for using Cox regression in SPSS: https://www.scalestatistics.com/cox-regression.html

I have no idea whether the information is correct, but it looks simple enough to point and click.

**Background** A proportional hazard model is one of numerous survival analysis which model unit risk, e.g. mortality, over time considering the risk e.g. mortality of the baseline population (covariates). Generally a hazard model will comprise a hazard function of the risk of events per unit time against the baseline covariates, such as age, smoking that sort of thing. Thus you wouldn't want to consider the mortality risk of a 25 year old to comparable to that of a 90 year old patient for a given factor, because there's loads of confounders.

At a minimum SPSS could be slow and this is a regression analysis, but you will find out soon enough.

---

It appears there's a forward prediction being sought here. A simple regression may be able to capture this, but the gist of the question appears to be a precise temporal forward prediction. Precise forward prediction are complicated and whilst there are definitely R packages that do this, at a precision level of forecast e.g. 2 vs 5 year this definitely ain't SPSS, a long way from its capabilities. The OP will need to learn R to do this study, it isn't trivial, or approximate from a hazard style regression. For this style of calculation significant epidemiological and clinical information is required. For example, these are often used in cancer studies, but precisely what the background epi/clinical is will be important as to whether a given R library is relevant.

**Summary** The answer to the question is SPSS can't do this and learning R is the way forward, because the OP will need to use an established library and this only place I'm aware that such libraries exist. The only alternative is to approximate using classical survival analysis.
> 1 votes

---
Tags: statistics
---


# Question
Title: Where can I find Single Cell Data with Location "Coordinates"?
Does single cell data typically have the following meta-data: the "coordinates" (e.g. on a tissue, adjacent tissues) saying where each cell in the sample was located relative to other cells? If not, is it possible to reconstruct this with other meta-data on the cells?

With the ultimate goal of working hands on with such location-tagged data, I am hoping for suggestions on the correct terms to search for, and even references to previous studies and/or easy to use public datasets.

For visual reference of the idea I have in mind, consider this picture:

I expect expression of any given gene in the skin cells in Condition 1 to be very different according to the location of the cell. For condition 2, I wouldn't expect the expression of a gene in the skin cells to be overly different, if at all. I want to see if I can formalize this idea using data mining techniques that were created for other purposes, but first I need the proper data.

**edit**: Here is an example of the format I eventually hope to work with.

**Fake table 1**: Each row represents measurements taken for a unique cell, while the last three columns have coordinates for its location. \begin{array}{r|lllllllll} \hline cell & phenotype1 & phenotype2 & \dots & gene1.expr & gene2.expr & \dots & x.loc.coord & y.loc.coord & z.loc.coord \\ \hline 1 & 0.69 & 1.34 & \dots & 1.91 & 0.21 & \dots & 1.12 & 0.05 & 1.09 \\ 2 & 0.34 & 0.92 & \dots & 1.74 & 2.03 & \dots & 0.57 & 0.46 & 0.24 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ n & 1.97 & 1.3 & 0.96 & 0.19 & 0.66 & \dots & 0.25 & 0.02 & 1.27 \\ \hline \end{array}

Alternatively, if the (x,y,z) coordinates of **Fake Table 1** are unavailable or not reconstructable, are there datasets which help us construct an adjacency list of cell pairs according to their location, such as **Fake Table 2**?

**Fake Table 2:** An adjacency table with pairs of cells which were next to each other when the measurement was taken. \begin{array}{r|cc} \hline pair & cell1 & cell2 \\ \hline 1 & 1 & 2 \\ 2 & 1 & 3 \\ 3 & 1 & 4 \\ 4 & 2 & 3 \\ 5 & 2 & 4 \\ \vdots & \vdots & \vdots \\ m-1 & n & n-2 \\ m & n & n-1 \\ \hline \end{array}

I come from a statistics background with a basic understanding of this type of data from talks I've attended, but not hands on experience. I mostly just want to explore this type of data to inform my future work. I am open to any type of scRNA-seq data with cell location "coordinates," if that's the correct thing to ask.

# Answer
> Does single cell data typically have the following meta-data: the "coordinates" (e.g. on a tissue, adjacent tissues) saying where each cell in the sample was located relative to other cells?

No. Typical scRNA-seq is just capturing random cells in a tube with no additional information. The technology you are looking for is spatial transcriptomics where you are measuring RNA levels at a particular location, but that is still not on the single-cell level.

> If not, is it possible to reconstruct this with other meta-data on the cells?

There have been efforts to reconstruct spatial relationship. For example: Satija et al. 2015. There should be more recent approaches that I am not aware of.
> 2 votes

# Answer
Is any of this what you want?

https://support.10xgenomics.com/spatial-gene-expression/datasets
> 1 votes

# Answer
I would suggest searching with terms such as "spatial transcriptomic" or "spatial proteomics". You could also search for platform-specific data such as "CyCIF", "CODEX", "MIBI", "slideseq", "visium" etc..

Here is a website that has curated some spatial data (http://www.spatialomics.org/SpatialDB/). We also recently released about 3TB of spatial imaging data that might be helpful. Unlike transcriptomic data that is often submitted in GEO (https://www.ncbi.nlm.nih.gov/geo/) there is still no one place where all of the spatial data is deposited and so you will find them only as supplemental data in published manuscripts. Hope that helps.
> 1 votes

# Answer
A good solution is to use Human papilloma virus (HPV) PCR (see below) or particularly RNA expression of protein E5, E6 and E7 of the viral genome in skin cells. High expression of E5, E6 or E7 would be diagnostic of skin cancer because of over deregulation of latency of the virus. This is because the viral control regions (L1/L2? - I forget) get disrupted mostly by genome integration leading to E5-E7 disrupting the human cell cycle control.

Over expression would be diagnositic because alot (most? all?) skin cancer is a result of HPV. HPV types 16 and 18 are the key, but most work done in this area focuses on cervical cancer because cervical HPV infection has a high mortality rate. I suspect skin cancer has a comparatively low mortality rate (but is of considerable imporant of course).

Alot of work has been done on the tandem arrays of HPV genomic insertion into the human geneome disrupting the regulatory genes leading to over production of E5-E7 and again is considered important in malignancy. However, its the protein expression of (E5), E6 and E7 that is critical. It is not clear whether human genome integration resulting in a viral tandem arrays that disrupt viral proteins L1/L2 will automatically trigger cancer, or whether it is merely a high risk factor (i.e. tandem array human genome insertions could still remain benign).
> 0 votes

---
Tags: rna-seq, single-cell, public-databases, differential-expression, gene-expression
---


# Question
Title: Is there a software that allows to convert a jpeg graph into dataset?
I have a very old graph as a .jpeg picture that I would like to convert into a dataset to be able to operate with it using R/python. I have no idea if it is even possible.

# Answer
Looks like **Graph Grabber** can do the job.

https://www.quintessa.org/software/downloads-and-demos/graph-grabber-2.0.2#download
> -2 votes

---
Tags: software-recommendation, graphs
---


# Question
Title: Error in seaborn plot " Horizontal orientation requires numeric `x` variable"
I am trying to plot a box plot with seaborn with the following code

```
plot = sns.boxplot(
    y='Min PPI distance',
    x='Synergy_percent',
    color='white',
    orient="h",
    data=synergy_df
)
plot.set(
    xlabel='Synergy (%)',
    ylabel='Min. distance',
);

```

The data `synergy_df` looks like

```
    Combination Min PPI distance    Synergy_percent
0   1003- Camptothecin & 1011- Navitoclax   2.0 22.2
1   1003- Camptothecin & 1017- Olaparib 2.0 6.4
2   1003- Camptothecin & 1022- AZD7762  2.0 40.4
3   1003- Camptothecin & 1025- SB216763 3.0 2.1
4   1003- Camptothecin & 1032- Afatinib 3.0 2.2

```

After running the above code I am getting the error `TypeError: Horizontal orientation requires numeric `x` variable.` I want to add that Synergy\_percent column has many zero values and about 1500 rows; I haven't tried it with a smaller number of rows.

I think the reason for the error is the wide range in the synergy\_percent values. The range is between 0 to 86.7. If I use xlim here what range should I use to divide the x-axis in 4 quarters (0-25-50-100)?

**Background** Seaborn is a Python library for more complex graph plotting than the commonly used `matplotlib` and is often used for boxplots and heatmaps usually using a `pandas` dataframe. The official documents are described here.

# Answer
Very quickly because I'm not around.

I think I know what's going on, but don't have enough time to figure the df. The key thing is `seaborn` needs to know where columns are `x` and which are `y` ... *within the data frame* . This information appears missing and the OP used `y='title', x='title` without reference to the dataframe. I'm pretty certain thats the bug here.

```
import seaborn as sns
import matplotlib.pyplot as plt

plot = sns.boxplot(
    color='white',
    y=synergy_df["Synergy_percent"],
    x=synergy_df["Min PPI distance"]
);
plt.show()

```

Alternative, which seems more what the OP wants ...

```
import seaborn as sns
import matplotlib.pyplot as plt

plot = sns.boxplot(
    color='white',
    # use your orientation command here
    x=synergy_df["Synergy_percent"],
    y=synergy_df["Min PPI distance"]
);
plt.show()

```

**Blurb** Normally the x-axis is the discrete boundary whilst the y-axis is continuous. In you example the x-axis feels like the continuous variable and you are producing a rotated boxplot - thats perfectly fine and seabourn will do that however its good to get it working first before rotating the axes.

**Final** The orientation of the dataframe might possibly be an issue, not sure unless its tried.

---

If there's a formatting issue the only trick I know is

```
 plot = sis.boxplot(
 .... # same old
 dodge=True,
 );

```

I don't know beyond, e.g. imposing percentile ranges onto a boxplot.
> 3 votes

# Answer
I'm not familiar with seaborn, but I'd guess from the error that there's at least one text value inside the `Synergy_percent` field that is encouraging the rest of the values to text. Either that, or spaces in the input file are messing up the parsing of the columns.
> 1 votes

---
Tags: visualization, python, graphs
---


# Question
Title: different results of RMSD for the same trajectory file
I want to analyze the RMSD of my trajectory files using MdTraj.

I used two different codes for this purpose:

the first code:

```
rep = 1
traj_name = "step7_%d.dcd" %rep
rmsd_total = []

while (os.path.isfile(traj_name)):
    t = md.load(traj_name, top = 'step5_input.pdb')
    top = t.topology
    protein_sel = top.select('backbone')
    protein_traj = t.atom_slice(protein_sel)
    rmsd = md.rmsd(protein_traj, protein_traj, 0)
    rmsd_without_zero = [i for i in rmsd if i != 0.0]
    rmsd_total.extend(rmsd_without_zero)

    rep+=1
    traj_name = "step7_%d.dcd" %rep

```

and the second code:

```
traj_1 = md.load('step7_1.dcd', top='step5_input.pdb')
traj_2 = md.load('step7_2.dcd', top='step5_input.pdb')
traj_3 = md.load('step7_3.dcd', top='step5_input.pdb')
traj_4 = md.load('step7_4.dcd', top='step5_input.pdb')
traj_5 = md.load('step7_5.dcd', top='step5_input.pdb')
traj_6 = md.load('step7_6.dcd', top='step5_input.pdb')
traj_7 = md.load('step7_7.dcd', top='step5_input.pdb')
traj_8 = md.load('step7_8.dcd', top='step5_input.pdb')
traj_9 = md.load('step7_9.dcd', top='step5_input.pdb')
traj_10 = md.load('step7_10.dcd', top='step5_input.pdb')
traj_11 = md.load('step7_11.dcd', top='step5_input.pdb')
traj_12 = md.load('step7_12.dcd', top='step5_input.pdb')
traj_13 = md.load('step7_13.dcd', top='step5_input.pdb')
traj_14 = md.load('step7_14.dcd', top='step5_input.pdb')
traj_15 = md.load('step7_15.dcd', top='step5_input.pdb')
traj_16 = md.load('step7_16.dcd', top='step5_input.pdb')
traj_17 = md.load('step7_17.dcd', top='step5_input.pdb')
traj_18 = md.load('step7_18.dcd', top = 'step5_input.pdb')
traj_19 = md.load('step7_19.dcd', top = 'step5_input.pdb')
traj_20 = md.load('step7_20.dcd', top = 'step5_input.pdb')

combine = traj_1 + traj_2 + traj_3 + traj_4 + traj_5 + traj_6 + traj_7 + traj_8 + traj_9 + traj_10 + traj_11 + traj_12 + traj_13 + traj_14 + traj_15 + traj_16 + traj_17 + traj_18 + traj_19 + traj_20

protein_sel = combine.topology.select('backbone')
protein_traj = combine.atom_slice(protein_sel)
rmsd = md.rmsd(protein_traj, protein_traj, 0)

```

I believe these two codes should result in the same RMSD. However, I get two different RMSD diagrams out of these two codes?

Can anyone help me what is the problems with my first code? (since I think the second code is more accurate while is more naive)

# Answer
> 1 votes
It's a reference issue. The second term.

Your first block makes a list of numpy arrays which are the RMSD against frame 0 of each of the stacks, while in your second you make a concatenation of frames, so your reference frame 0 is common.

The first code block could be smartened by splitting (abstracting) the file valid name enumeration and the reading/analysis. Namely, generate a list of valid filenames and *then* analyse. This separation would allow you recycle more easily the code: in the seocnd code block you have the file list problem again, which is instead resolved by a copypaste job —very inflexible. A one liner list comprehension (or, less so, a `functools.reduce`) is a nice option.

---
Tags: python
---


title: Making a bed file for RSeQC

text:
# Question
Title: Making a bed file for RSeQC
I making a bed file for RSeQC, so it can do things like compute the number of reads from exons, introns, 5"UTRs, etc.

I want to use a bed file that corresponds to my GTF file, so I use gtf2bed to make a bed file, like this:

```
awk '{ if ($0 ~ "transcript_id") print $0; else print $0" transcript_id \"\";"; }' Homo_sapiens.GRCh38.gtf | gtf2bed - > $Homo_sapiens.GRCh38.bed

```

This file throws some errors, among other when using it with junction\_annotation.py:

```
Reading reference bed file:  /rst1/2017-0205_illuminaseq/data/references/Reference_Genomes/GRCh38.95/Homo_sapiens.GRCh38.87.bed  ... Traceback (most recent call last):
  File "/rst1/2017-0205_illuminaseq/tools/snakemake_conda_env/f4c16777/bin/junction_annotation.py", line 125, in <module>
    main()
  File "/rst1/2017-0205_illuminaseq/tools/snakemake_conda_env/f4c16777/bin/junction_annotation.py", line 109, in main
    obj.annotate_junction(outfile=options.output_prefix,refgene=options.ref_gene_model,min_intron=options.min_intron, q_cut = options.map_qual)
  File "/rst1/2017-0205_illuminaseq/tools/snakemake_conda_env/f4c16777/lib/python2.7/site-packages/qcmodule/SAM.py", line 3762, in annotate_junction
    exon_starts = map( int, fields[11].rstrip( ',\n' ).split( ',' ) )
ValueError: invalid literal for int() with base 10: 'gene_version'

```

and when I compare this bed file to the RSeQC provide (dated) one these are the differences:

My bed file I created:

```
1       11868   12227   ENSG00000223972 .       +       havana  exon    .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; exon_number "1"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-202"; transcript_source "havana"; transcript_biotype "processed_transcript"; exon_id "ENSE00002234944"; exon_version "1"; tag "basic"; transcript_support_level "1";
1       11868   14409   ENSG00000223972 .       +       havana  gene    .       gene_id "ENSG00000223972"; gene_version "5"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_id "";
1       11868   14409   ENSG00000223972 .       +       havana  transcript      .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-202"; transcript_source "havana"; transcript_biotype "processed_transcript"; tag "basic"; transcript_support_level "1";
1       12009   12057   ENSG00000223972 .       +       havana  exon    .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; exon_number "1"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; exon_id "ENSE00001948541"; exon_version "1"; tag "basic"; transcript_support_level "NA";
1       12009   13670   ENSG00000223972 .       +       havana  transcript      .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; tag "basic"; transcript_support_level "NA";
1       12178   12227   ENSG00000223972 .       +       havana  exon    .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; exon_number "2"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; exon_id "ENSE00001671638"; exon_version "2"; tag "basic"; transcript_support_level "NA";
1       12612   12697   ENSG00000223972 .       +       havana  exon    .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; exon_number "3"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; exon_id "ENSE00001758273"; exon_version "2"; tag "basic"; transcript_support_level "NA";
1       12612   12721   ENSG00000223972 .       +       havana  exon    .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; exon_number "2"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-202"; transcript_source "havana"; transcript_biotype "processed_transcript"; exon_id "ENSE00003582793"; exon_version "1"; tag "basic"; transcript_support_level "1";
1       12974   13052   ENSG00000223972 .       +       havana  exon    .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; exon_number "4"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; exon_id "ENSE00001799933"; exon_version "2"; tag "basic"; transcript_support_level "NA";
1       13220   13374   ENSG00000223972 .       +       havana  exon    .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; exon_number "5"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; exon_id "ENSE00001746346"; exon_version "2"; tag "basic"; transcript_support_level "NA";

```

The RSeQC provided bed file:

```
1       11868   14409   ENST00000456328 0       +       11868   14409   0       3       359,109,1189,   0,744,1352,     DDX11L1
1       12009   13670   ENST00000450305 0       +       12009   13670   0       6       48,49,85,78,154,218,    0,169,603,965,1211,1443,        DDX11L1
1       17368   17436   ENST00000619216 0       -       17368   17436   0       1       68,     0,      MIR6859-1
1       14403   29570   ENST00000488147 0       -       14403   29570   0       11      98,34,152,159,198,136,137,147,99,154,37,        0,601,1392,2203,2454,2829,3202,3511,3864,10334,15130,   WASH7P
1       29553   31097   ENST00000473358 0       +       29553   31097   0       3       486,104,122,    0,1010,1422,    MIR1302-2
1       30266   31109   ENST00000469289 0       +       30266   31109   0       2       401,134,        0,709,  MIR1302-2
1       30365   30503   ENST00000607096 0       +       30365   30503   0       1       138,    0,      MIR1302-2
1       34553   36081   ENST00000417324 0       -       34553   36081   0       3       621,205,361,    0,723,1167,     FAM138A
1       35244   36073   ENST00000461467 0       -       35244   36073   0       2       237,353,        0,476,  FAM138A
1       52472   53312   ENST00000606857 0       +       52472   53312   0       1       840,    0,      OR4G4P

```

Why are the files so different? How would I make a simpler bed file form my GTF file?

# Answer
# UCSC Utilities

Such output results when using the UCSC utilities `gtfToGenePred` and `genePredToBed` in series. Builds for macOS (x86\_64) and Linux are available directly from UCSC. Alternatively, Conda builds are available through the Bioconda channel, under the package names `ucsc-gtftogenepred` and `ucsc-genepredtobed`.

## Example: *Homo sapiens* (Ensembl Release 106)

```
## (optional install through Conda)
#mamba install -c conda-forge -c bioconda ucsc-gtftogenepred ucsc-genepredtobed

## download GTF
wget ftp://ftp.ensembl.org/pub/release-106/gtf/homo_sapiens/Homo_sapiens.GRCh38.106.gtf.gz

gzip -cd Homo_sapiens.GRCh38.106.gtf.gz |\
  gtfToGenePred /dev/stdin /dev/stdout |\
  genePredToBed /dev/stdin /dev/stdout |\
  head

```

**Output BED**

```
1   1211339 1214153 ENST00000379236 0   -   1211554 1214127 0   7   286,129,197,67,102,123,171, 0,364,602,1298,1652,2323,2643,
1   1211339 1214138 ENST00000497869 0   -   1214138 1214138 0   5   493,197,67,794,156, 0,602,1298,1652,2643,
1   1212018 1213498 ENST00000453580 0   -   1213498 1213498 0   4   120,67,102,104, 0,619,973,1376,
1   1203507 1206571 ENST00000328596 0   -   1203590 1206571 0   4   453,88,123,187, 0,891,1862,2877,
1   1203507 1206592 ENST00000379268 0   -   1203843 1206571 0   5   461,203,88,123,208, 0,526,891,1862,2877,
1   1203843 1205680 ENST00000486728 0   -   1203843 1205463 0   4   125,203,88,311, 0,190,555,1526,
1   1203843 1206571 ENST00000379265 0   -   1203843 1206571 0   5   125,182,88,123,187, 0,211,555,1526,2541,
1   1471764 1497848 ENST00000673477 0   +   1471884 1495817 0   16  325,77,102,60,70,166,70,156,57,126,125,52,71,168,109,2364,  0,5509,6879,7284,9102,10373,10780,13251,14017,14345,14779,16098,17439,18492,18798,23720,
1   1478025 1497848 ENST00000472194 0   +   1497848 1497848 0   14  720,60,70,166,70,156,57,126,125,52,71,168,109,2364, 0,1023,2841,4112,4519,6990,7756,8084,8518,9837,11178,12231,12537,17459,
1   1479048 1482662 ENST00000378736 0   +   1482662 1482662 0   4   60,70,166,118,  0,1818,3089,3496,

```
> 3 votes

---
Tags: rna-seq, bed, gtf, rseqc
---


# Question
Title: How can I call one file from a channel that's composed of several files in Nextflow?
In Nextflow I have a process that creates a channel that holds several files in a tuple. I'm trying to call only one of those files as input for a subsequent process.

So basically if I run `channel_name.view()` the channel shows `[/home/nextflow/work/e4/bc47bceb4569cffc01bd76e49245e9/RAxML_bestTree.tree, /home/nextflow/work/e4/bc47bceb4569cffc01bd76e49245e9/RAxML_bipartitions.tree, /home/nextflow/work/e4/bc47bceb4569cffc01bd76e49245e9/RAxML_bipartitionsBranchLabels.tree, /home/nextflow/work/e4/bc47bceb4569cffc01bd76e49245e9/RAxML_bootstrap.tree, /home/nextflow/work/e4/bc47bceb4569cffc01bd76e49245e9/RAxML_info.tree]` . I need just one file so I'm using the full name in the `input` section of the process.

```
input:
file 'RAxML_bipartitionsBranchLabels.tree' from channel_name

```

That does not find the file. On the other hand, if I use `.flatten()`, like this:

```
input:
file 'RAxML_bipartitionsBranchLabels.tree' from raxml_out_ch.flatten()

```

It works but reads some files before reaching to the one that matches and so the process results in error even after completing the task for this file properly and gives several outputs (for each file it looks into), some of them empty and one with the result.

```
[7c/2971f1] Submitted process > tree_visualisation (4)
[8d/dd129b] Submitted process > tree_visualisation (1)
[cb/0f91f5] Submitted process > tree_visualisation (5)
[6c/fb23c8] Submitted process > tree_visualisation (3)
ERROR ~ Error executing process > 'tree_visualisation (5)'

Caused by:
  Process `tree_visualisation (5)` terminated with an error exit status (1)

Command executed:

  Rscript tree.R RAxML_bipartitionsBranchLabels.tree phylogram.png

```

I'm looking for a way to filter this input so that the process only runs for this file and doesn't try with the others.

How can I do that?

# Answer
If, in your input definition, you specify a string to represent your filename and provide multiple file inputs as a collection, a numerical suffix is appended to each filename:

> When a target file name is defined in the input parameter and a collection of files is received by the process, the file name will be appended by a numerical suffix representing its ordinal position in the list.

When you use the flatten operator, the collection is flattened such that each entry becomes a separate entry in the output channel. So for each of the five 'tree' files, five tasks are spawned: one for each 'tree' file. These files are each localized using the filename specified: 'RAxML\_bipartitionsBranchLabels.tree\`. This usually isn't what you want to do.

What you're looking for, I think, is some way to select the file called 'RAxML\_bipartitionsBranchLabels.tree' from the initial collection of files. There are many ways to do this. For example if you wanted to flatten the output channel, it's possible to select the file you want using the first operator and a closure:

```
bipartitionsBranchLabels = raxml_out_ch.flatten().first {
    it.name == 'RAxML_bipartitionsBranchLabels.tree'
}

```

However you should only seldom need to do this. If your pipeline requires just one of these files in a downstream process, my suggestion would be to alter the upstream output declaration to include a separate output channel for this purpose:

```
process RAxML {

    output:
    path 'RAxML_*.tree' into raxml_out_ch
    path 'RAxML_bipartitionsBranchLabels.tree' into raxml_bipartitionsBranchLabels

    """
    touch RAxML_bestTree.tree
    touch RAxML_bipartitions.tree
    touch RAxML_bipartitionsBranchLabels.tree
    touch RAxML_bootstrap.tree
    touch RAxML_info.tree
    """
}

```
> 2 votes

---
Tags: nextflow
---


# Question
Title: What is & how to solve File error: my.xml.state (Remote I/O error)?
I caught the next exception during my phylogeographical analysis in **BEAST 2** with *GEO\_SPHERE*. What could be the reason? & how to evade this in the future?

```
      ...
      856000000     -3662.2647      5969.5577     -9631.8225 42m16s/Msamples
      857000000     -3522.3540      6105.9849     -9628.3389 42m5s/Msamples
      858000000     -3621.7922      6022.9359     -9644.7282 41m54s/Msamples
      859000000     -3463.1050      6107.9766     -9571.0817 41m43s/Msamples
      860000000     -3470.6562      6160.7143     -9631.3706 41m35s/Msamples
File error: VP1_test.xml.state (Remote I/O error)
java.lang.RuntimeException: An error was encounted. Terminating BEAST
        at beast.app.util.ErrorLogHandler.publish(Unknown Source)
        at java.logging/java.util.logging.Logger.log(Logger.java:979)
        at java.logging/java.util.logging.Logger.doLog(Logger.java:1006)
        at java.logging/java.util.logging.Logger.log(Logger.java:1029)
        at java.logging/java.util.logging.Logger.severe(Logger.java:1776)
        at beast.app.beastapp.BeastMain.<init>(Unknown Source)
        at beast.app.beastapp.BeastMain.main(Unknown Source)
        at java.base/jdk.internal.reflect.NativeMethodAccessorImpl.invoke0(Native Method)
        at java.base/jdk.internal.reflect.NativeMethodAccessorImpl.invoke(NativeMethodAccessorImpl.java:62)
        at java.base/jdk.internal.reflect.DelegatingMethodAccessorImpl.invoke(DelegatingMethodAccessorImpl.java:43)
        at java.base/java.lang.reflect.Method.invoke(Method.java:566)
        at beast.app.beastapp.BeastLauncher.run(Unknown Source)
        at beast.app.beastapp.BeastLauncher.main(Unknown Source)

```

Unfortunately, googling this error with **"BEAST"** gave me nothing. If without it, there are a lot of irrelevant problems.

---

**Update**

I have just found out that this error happened simultaneously onto all my nodes (different processors, RAM, & common hards & SSDs). All logs of BEAST 2 analyses have it. What is stranger, some of the analyses pressed on after this & someones were broken. To the present I got more interested not in the broken ones but could this exception spoil those which continued their work?

# Answer
> 2 votes
You can restart this at the place you left 'cause is Beast2.

```
  860000000     -3470.6562      6160.7143     -9631.3706 41m35s/Msamples

```

However, this is almost 1 billion replications and should be sufficient for convergence. The fact it ends on a round number might infer a algorithm limit is in place. If the Beast2 was restarted from this point and it abruptly ended this would be a very likely reason.

The alternative is an error external to the algorithm such as a RAM error - it has been running a long time.

Whatever the reason I would simply perform ESS statistics and convergence diagnostics on what you have now. There's no real reason to continue if it has converged for 10e6 replications maximum and this is much, much bigger than that.

---

The terminated processes have not spoiled the work again you can simply input the terminated chains into the restart for Beast2. This is because Beast2 uses `checkpoints`. Brief information on how to resume from the last checkpoint is highlighted here https://www.beast2.org/2014/04/07/checkpointing-tricks.html. Convergence checks are important.

The situation is only is only catastrophic for Beast1 where you would have needed to start from the beginning.

---
Tags: phylogenetics, phylogeny, beast
---


# Question
Title: What is the different between the I-TASSER, phyre2, SWISS-model in the 3D tertiary structure?
What is the difference between the I-TASSER, phyre2, and SWISS-model in the 3D tertiary structure?

How do they get the results? When I did a prediction, I got a similar result for the highest template in I-TASSER and Swiss, while it was different with Phyre2. What does that mean?

# Answer
> 2 votes
This is a common Q.

Briefly, you have some techniques here that get used separately or combined:

* Ab initio: assemble the protein from scratch, often via a library fragmented peptides (3-mer, 6-mers, 9-mers) based on real structures.
* Template threading: given a protein structure alter it to make it like a give sequence. SWISS-MODEL, Rosetta's SimpleThreadingMover, ModBase, Phyre template threading etc. This is why some of these keep ligands and deal well with multimers. For missing loops in the template these are filled by ab initio methods.
* Hybridisation: give a bunch of models mix and match them for the best energy. I-Tasser, RosettaCM (which uses SimpleThreadingMover), Phyre (which uses Phyre template threading)
* Ab initio with contraints, such as evolutionary covariance. EVFold, AlphaFold2, trRosetta. Often mixed with template based approaches.

Wikipedia has more, but papers are very helpful despite the car-salesmanship.

One thing to note is that protein can adopt different conformations. The AlphaFold2/ColabFold notebooks give multiple models: these may be legit alt conformations, e.g. inactive vs active etc etc

Here is a SwissModel vs AF2: Different structures of same protein from two softwares

Other helpful As:

---
Tags: sequence-alignment, proteins, protein-structure, sequence-analysis, 3d-structure
---


# Question
Title: Reference variant detected as altered one in bam file
I received (from manufacturer) several .bam files and I used four callers (samtools, freebayes, haplotypecaller, deepvariant) to find some sequence variants. In obtained .vcf files, I took a closer look to some calls. I found interesting, homozygous one rs477033 (C/G Ref/Alt) with flag 'COMMON=0' and very low MAF. I also checked the .bam file (in IGV) and this rs itself in ensembl database. It looks like this:

The problem is: in this case ref and alt are switched. The G is reference (ancestral) allele and C is altered. Nevertheless in the .bam file it's interpreted in opposite way. Of course it's false positive variant - there's no mutation there at all. All variant callers returned it as high quality variant. My. bam files are full of them.

**What went wrong in this case? Problem with mapping?**

**How can I fix this thing, to get only novel and rare true positive variants?**

# Answer
The ancestral allele is not always the reference allele. Ref and alt aren't switched, ref is just rare in at least the population you're looking at (it appears to be rare in humans as a whole). So, the C is indeed reference and G is indeed altered, it's just that people homozygous for the G are the overwhelming majority. This is evidenced by the rs477033 variant being named NC\_000001.11:g.165890393C\>G. So, in short, nothing went wrong, all tools functioned correctly.

You'll notice a bunch of these over time, either make a database of them (or find one, someone may have made one) and filter them out or see if maybe the 1000 genomes reference genome has some of these modified (or maybe the T2T genome that was recently published).
> 3 votes

---
Tags: sam, snp, ngs, variant-calling
---


# Question
Title: Minimap2 detects too many 16S sequences in metagenome
I'd like to extract bacterial sequences from a metagenome where ~99% of contigs are from a host insect. My current protocol minimaps every 16S sequences present in the sample, then blasts all those sequences to create a library of full genomes that I minimap against. I'm using the full metagenome as my reference file. I wasn't sure if that makes a difference, but figured I'd mention it. I ran the following command to retrieve 16S sequences that are present in the sample:

```
minimap2 All_reads.fastq 16S_reference_library.fasta > 16S_map.paf

```

Minimap2 retrieves every entry in 16S\_reference\_library.fasta. It's highly improbable some of those sequences are present, or that 25,000 species are present at all. Are there any parameters I could adjust in Minimap2 that might help with this?

Alternatively, I'm trying to implement Krakenuniq and Centrifuge, but I'm concerned that, because I don't have a reference genome for the host insect, they'll have a very high false positive rate.

# Answer
> 3 votes
16S is a very conserved sequence, which is why it's used for targeted phylogenetic analysis; it makes it easy to amplify and analyse. Unfortunately that conservation is an issue with minimap2, which is built around the idea that matching scattered subsequences within a sequence is good enough for identifying matches to the entire sequence. With a highly-conserved sequence, there's a higher chance that subsequences will be identical, meaning that diverged sequences will appear to be the identical from the perspective of the mapper.

Minimap2 can be tweaked to make it less sensitive and more specific, but you'd be better off using a different mapping tool, especially one designed for phylogenetic analysis and conserved sequences. I recommend kraken2.

pre-built 16S indexes for kraken2 can be found here:

https://benlangmead.github.io/aws-indexes/k2

If you really want to map specifically to your own generated metagenome, then you could try LAST (with results filtered through `last-map-probs`), which I have found to be a little bit better in dealing with specificity issues.

Aside: I notice that you've stated "bacterial sequences from a *metagenome*". If these are, in fact, metagenomic reads (and not targeted 16S reads), then you would be better off using one of the larger metagenome bacterial databases (e.g. plusPF or plusPF-8, depending on computer memory) instead of concentrating only on 16S-matching reads. This will provide much better resolution; 16S only really works down to the family level.

# Answer
> 2 votes
You are taking reads as the reference. Given the similarity of 16S, you will find hits for most 16S in other species. You have to filter out poor hits. This is the faster way if you don't have many contigs.

Alternatively, switch the order of reference and reads. But in this case, you need to use a high k-mer occurrence cut off as most of your k-mers will occur \>10,000 times. Something like:

```
minimap2 -c -f30000 -N100 ref-16S.fa reads.fq

```

This will be fairly slow as there are too many hits in the database but it will be easier to process the output.

You can try both and see what works for you.

---
Tags: metagenome, minimap2
---


title: Parsing pre-2007 SMILES string

text:
# Question
Title: Parsing pre-2007 SMILES string
How would one parse the SMILES string

```
 BrC[2]:C[3]:C(:CH:CH:CH:@2):CH:CH:CH:CH:@3

```

I rely on tools like `rdkit` and `OpenBabel` to parse SMILES, but both tools aren't able to parse this string.

More specifically, this SMILES string comes from the supplementary information of a 2007 paper on solubility prediction. After doing some searching, it seems that this paper was received prior to the 2007 standardization of OpenSMILES (this might be why canonical parsers like `rdkit` have trouble with it).

After searching a bit more, I found an older PubChem sketcher that was able to convert the SMILES in question into a more recognizable format

```
 (BrC1=C(C(=C(C2=C1C(=C(C(=C2[H])[H])[H])[H])[H])[H])[H])

```

(parsing this was successful via `rdkit`, see molecule image below) but still need a way to convert similarly-formatted SMILES strings into the OpenSMILES standard/canonical format en mass.

# Answer
I couldn't find a way to directly convert these SMILES strings to a more canonical format. However, thanks to @Matteo Ferla, I was able to get around the problem entirely by looking up the chemical names and resolving the canonical SMILES programmatically.

More specifically, I used the CAS no column given in the supplementary information of the referenced paper to pull the SMILES via Chemical Identifier Resolver. See this answer for an explanation on how to do so.

**EDIT**: Some of the CAS nos couldn't be resolved, so I manually converted them using the PubChem sketcher tool mentioned in the question. There are quite a few (~160) of such molecules, so if anyone needs the full parsed dataset but doesn't want to spend the time manually parsing those molecules, let me know.
> 1 votes

---
Tags: perl, computational-biochemistry, rdkit
---


title: Protein-Ligand in PDB

text:
# Question
Title: Protein-Ligand in PDB
I am a newbie in bioinformatics. In PDB, while searching for proteins, I sometimes get protein with ligand, such as PDB ID 1J1Z. In cases I get pure protein structure, in some cases I get protein-with ligand.

I am having hard time understanding why I can have ligand in PDB structure. I am PDB is supposed to give me pure protein 3D structure. Then how can we have ligand in protein structure??

And sorry if my question is stupid.

# Answer
> 2 votes
In some protein crystallography studies, the protein crystals are grown with the substrate/ligand/inhibitor added, thus the resulting structure contains both the ligand and the protein atom coordinates. This enables studying the protein ligand complex to understand binding, specificity, mechanism of catalysis...

On the other hand in some experiments protein crystals are grown without the ligand, thus the resulting structure contains coordinates just for the protein atoms.

On the third hand, the lucky one, your protein, overexpressed in its natural host \[sometimes even in different one\], co-purify with its own natural ligand, so the resulting structure contains coordinates for both the protein atoms and the sought for ligand.

To understand ligand binding the protein structure without the ligand should be compared to the structure of the protein-ligand complex.

Many published papers report both types of crystals.

For example the crystal structure of pectin methylesterase in complex with hexasaccharide VI: 2NTP is from the publication https://doi.org/10.1038/sj.emboj.7601816 which is the source of six other structures. These structures include the enzyme with different substrates, as well as the active site mutants (D178A) all of which were prepared to understand the molecular basis of substrate specificity. The conclusion after studying all of the mentioned structures is nicely summed in the sentence from the abstract:

> The preferential binding of methylated sugar residues upstream of the catalytic site, and demethylated residues downstream, drives the enzyme along the pectin molecule and accounts for the sequential pattern of demethylation produced by both bacterial and plant PMEs.

---
Tags: protein-structure, pdb
---


# Question
Title: Why must a maximal non-branching path be a contig?
The following is from *Bioinformatics Algorithms*:

> Fortunately, we can derive contigs from the de Bruijn graph. A path in a graph is called non-branching if in(v) = out(v) = 1 for each intermediate node v of this path, i.e., for each node except possibly the starting and ending node of a path. A maximal non-branching path is a non-branching path that cannot be extended into a longer non-branching path. We are interested in these paths because the strings of nucleotides that they spell out must be present in any assembly with a given k-mer composition. For this reason, contigs correspond to strings spelled by maximal non-branching paths in the de Bruijn graph. For example, the de Bruijn graph below, constructed for the 3-mer composition of TAATGCCATGGGATGTT, has nine maximal non-branching paths that spell out the contigs TAAT, TGTT, TGCCAT, ATG, ATG, ATG, TGG, GGG, and GGAT. In practice, biologists have no choice but to break genomes into contigs, even in the case of perfect coverage (like in the figure below), since repeats prevent them from being able to infer a unique Eulerian path.

But this quote doesn't make sense to me:

> We are interested in these paths because the strings of nucleotides that they spell out **must be present in any assembly with a given k-mer composition**.

Is it not possible that we simply have a missing read that would have otherwise branched this path? How do we know that a maximal non-branching path *must* be a contig?

# Answer
The definition of contig is derived from assembly algorithms. A contig is not necessarily present in the genome. PS: i.e. "contig" is algorithmically defined, not biologically defined.

> Is it not possible that we simply have a missing read that would have otherwise branched this path?

Yes, it is possible. Then you are likely to get a misassembled contig.
> 1 votes

---
Tags: sequence-alignment, ngs, assembly
---


# Question
Title: How to properly annotate sequence variants and errors in a GenBank file format and how to keep track of successive versions of a GenBank file
### EDIT

By thinking over the issue, I realized it was mixing two linked but separate questions:

* **Q1. How to handle the situation when there are divergences between a GenBank file and an actual or theoretical molecule:**

+ Should the difference be indicated in the feature table? If so, with which feature key? And if not, how should it be indicated or what should be done?
    + Further, when indicating a divergence using a feature key, which should be considered the reference: the file or the molecule?
* **Q2. How to keep track of changes between versions of a GenBank file when corrections or updates are made?**

I have thus re-titled the question accordingly, and by making further researches into the issues I found what seems to me the appropriate solutions (see my answer below).

### ORIGINAL QUESTION

**Background** I am maintaining a database of sequence files describing a collection of bacterial strains and plasmids that we have in our lab.

I am confronted to two sorts of issues:

* some sequence files contain errors (see comments below for examples of errors)
* other sequence files describe theoretical sequences of molecules that themselves contain variations (mutations, insertion or deletions) that occurred during the molecular cloning process.

**Objective** I want to annotate the sequences to mark those errors and variations using the GenBank format. However, from the feature keys of the GenBank format described on the INSDC website, it is not clear to me which feature key should be used in each case, essentially this is either 'variation' or 'misc\_description'.

**Example** To better explain the confusion, here below are what seems to me the most relevant feature keys, and my comment of the problem I have with them:

```
Feature key           variation
Definition            a related strain contains stable mutations from the same
                      gene (e.g., RFLPs, polymorphisms, etc.) which differ
                      from the presented sequence at this location (and
                      possibly others);
...
Comment               used to describe alleles, RFLP's,and other naturally
                      occurring mutations and  polymorphisms;
                  >>> variability arising as a result of genetic
                      manipulation (e.g. site directed mutagenesis) should
                      be described with the misc_difference feature;

```

According to the line marked `>>>` above, the feature key `misc_difference` should be used...

```
Feature Key           misc_difference
Definition            feature sequence is different from that presented in the
                      entry and cannot be described by any other difference
                      key (old_sequence, variation, or modified_base);
...
Comment               the misc_difference feature key should be used to
                      describe variability that arises as a result of
                      genetic manipulation (e.g. site directed mutagenesis);
                  >>> use /replace="" to annotate deletion, e.g.
                      misc_difference 412..433
                                      /replace=""

```

But now, according to the example marked `>>>` above, it means the variation (here a deletion) is NOT in the considered entry...

**Question** So how are variations (mutation, insertion, deletion, errors) correctly described in Genbank given that are actually present in an entry ?

# Answer
> 1 votes
Your core question is:

> (mutations, insertion or deletions) that occurred during the molecular cloning process.

I asked NCBI directly and their response is:

> For natural genetic variation use 'Variation' feature;

> For 'artificial' genetic variation (lab-produced by genetic modification) then use 'misc\_difference'.

In summary in all scenarios its singly the 'misc\_differences' that are used (see Conclusion).

So in my scenario where I want to screen out targeted mutagenesis its the 'misc\_difference', i.e. the second option of the two you present. This is also true in your second scenario because it is changes possibly unintentional, possibly random mutagenesis / or just unintentional errors during cloning. Finally, this is also true for SOME of the 'errors' you are reporting in comments:

* for the first i, ii and iii that appears to be just correcting stuff - administration, just cross-check with those submitting.
* for 'variations' the second i and ii in comments, that falls into the 'misc\_difference', because a form of mutagenesis, albeit probably unintentional.

---

**Conclusion** The 'variation' feature description isn't clear given the description you reported against NCBI's official Helpdesk response. Its a fair question, its a clear answer, but doesn't explain why the description is confusing.

---

On 'errors' of mis-identification, e.g. correctly identified e.g. via local Blast or Blat, e.g. sequenced the wrong target, I would use an internal marker, e.g. `returned to contributor for comment`. Pass it back to the contributor and ask for their input and depending what they say simply mark the record in the 'Comments section' e.g. `originally defined as XX`. The contributor could have made all sorts of different errors, e.g. annotation - copied and pasted the wrong stuff - experimental anything and you need to give them a chance to respond. Uploading these records on an external database, i.e. Genbank, needs a joint agreement with all parties concerned because it is certainly valid information for someone.

---

Note, the official blurb here, doesn't appear informative:

> Disrupted genes and gene fragments
>
> Sometimes a genome will have adjacent or nearby genes that seem to be only part of a protein. In many cases these indicate a possible problem with the sequence and/or annotation. A related issue is the presence of internal stop codons in the conceptual translation of a CDS that looks like it should be a real CDS. These problems may be due to a variety of reasons, including **mutations** or sequencing or assembly artifacts. They can be annotated in a number of ways:
>
> Annotate the gene with /pseudo to indicate that there is a problem with the gene. Note that this qualifier does NOT mean that the gene is a pseudogene. (see point 2, below, if it is known that the gene IS a pseudogene) If multiple gene fragments were present initially, then add a single gene feature which covers all of the potential coding regions and add the pseudo qualifier. If known, a note qualifier may be added indicating why this gene is disrupted, for example:

```
 1        200     gene
                         gene    phoA
                         gene_desc     alkaline phosphatase
                         locus_tag     OBB_0001
                         pseudo
                         note    nonfunctional due to frameshift

```

# Answer
> 1 votes
### Overview

To properly address the questions, each possible case must be listed:

* **Case 1**: divergence between sequences (sequence in the file vs sequence of the molecule); two cases:

**a**\- (wrong) sequence in the file vs (correct) sequence in the molecule

**b**\- ("correct", theoretical) sequence in the file vs (actual, mutated) sequence in the molecule
* **Case 2**: divergence between a feature annotation in the file vs the correct ("true" or reference) biological feature; three cases:

**a**\- wrong location of the feature

**b**\- wrong feature name

**c**\- variant sequence of the feature

### Short answer

**Case 1a**: The sequence in the GenBank file must be corrected and the *version number* of the file must be incremented (see **Q2**).

**Case 1b**: A new GenBank file consistent with the actual molecule must be created and the *misc\_difference* feature key must be used to annotate the variations where relevant.

**Cases 2a and 2b**: Feature information in the GenBank file must be corrected and the *Modification Date* of the file must be updated (see **Q2**).

**Case 2c**: The *misc\_difference* feature must be used to annotate the variation and the *Modification Date* of the file must be updated if necessary (see **Q2**).

### Long answer

#### Q1. How to handle the situation when there are divergences between a GenBank file and an actual or theoretical molecule?

##### The *misc\_difference* feature key

First, let's analyse the situation with the *misc\_difference* feature key. According to the INSDC documentation:

```
Feature Key           misc_difference

Definition            feature sequence is different from that presented
                      in the entry and cannot be described by any other
                      difference key (old_sequence, variation, or modified_base);

Optional qualifiers   /allele="text"
                      /clone="text"
                      /db_xref="<database>:<identifier>"
                      /experiment="[CATEGORY:]text"
                      /gene="text"
                      /gene_synonym="text"
                      /inference="[CATEGORY:]TYPE[(same species)][:EVIDENCE_BASIS]"
                      /locus_tag="text" (single token)
                      /map="text"
                      /note="text"
                      /old_locus_tag="text" (single token)
                      /phenotype="text"
                      /replace="text"
                      /standard_name="text"

Comment               the misc_difference feature key should be used to
                      describe variability that arises as a result of
                      genetic manipulation (e.g. site directed mutagenesis);
                      use /replace="" to annotate deletion, e.g.
                      misc_difference 412..433
                                      /replace=""

```

From the Definition above, there is still some ambiguity as to what situation is described. Indeed, "*feature sequence is different from that presented in the entry*" can describe both of the following cases:

* the sequence in the file and/or in the molecule is not the one expected for the feature annotated
* the sequence in the molecule is not the one represented in the file

*Note: In fact, as we will see below, there should be no divergence between the sequence in the file and the sequence of the molecule.*

As for the Comment that says that "*the misc\_difference feature key should be used to describe variability that arises as a result of genetic manipulation (e.g. site directed mutagenesis)*", it does not really clarify which of the two cases above is the one intended. Further, the example given in the Comment, which at first sight seems more enlightening, still carries some ambiguity. It says that to annotate a deletion, one should use:

```
misc_difference 412..433
                /replace=""

```

All in all, it implies that the deletion is in the molecule (vs in the represented sequence), that it arose as a result of genetic manipulation and is therefore indicated in the file as explained above.

However, it is very counter-intuitive, and not very practical when searching sequences of variants among sequence files, to have a sequence file that is supposed to describe a molecule and in which the sequence is actually different from the sequence encountered in the molecule it is supposed to describe.

And indeed, this view goes against common observed recent practice where the *misc\_difference* feature key is used to indicate how to go from the sequence represented in the sequence file to another version:

* **Example 1**: https://www.ncbi.nlm.nih.gov/nuccore/D90147.1

```
     misc_difference 308..310
                     /gene="mucA"
                     /note="conflict: deletion"
                     /citation=[1]
                     /replace=""

```

The indicated deletion is not in the sequence but is a way to go from the sequence to a previously reported CDS sequence for the annotated gene.

Translation of the CDS sequence in the GenBank file:

```
MKVDIFESSGASRVHSIPFYLQRISAGFPSPAQGYEKQELNLHEYCVRHPSATYFLRVSGSSMEDGRIHDGDVLVVDRSLTASHGSIVVACIHNEFTVKRLLLRPRPCLMPMNKDFPVYYIDPDNESVEIWGVVTHSLIEHPVCLR*

```

Previously reported CDS (see: https://www.ncbi.nlm.nih.gov/nuccore/M13388; hyphen/dash in the sequence added by me to mark the difference)

```
MKVDIFESSGAS-VHSIPFYLQRISAGFPSPAQGYEKQELNLHEYCVRHPSATYFLRVSGSSMEDGRIHDGDVLVVDRSLTASHGSIVVACIHNEFTVKRLLLRPRPCLMPMNKDFPVYYIDPDNESVEIWGVVTHSLIEHPVCLR*

```

Ref.: https://pubmed.ncbi.nlm.nih.gov/1902211/
* **Example 2**

In each of the GenBank files below, the *misc\_difference* feature key is used to indicate *how to go from the (mutated) sequence in the GenBank file to the sequence of the "wild type"*.

+ *Escherichia coli* *uidA* gene for beta-glucuronidase mutant E413G, clone M2

    https://www.ncbi.nlm.nih.gov/nuccore/AM180563.1

    ```
     misc_difference 1238
                   /gene="uidA"
                   /experiment="experimental evidence, no additional details
                   recorded"
                   /note="site directed mutagenesis"
                   /compare=S69414.1
                   /allele="wild type"
                   /phenotype="Active GUS"
                   /replace="a"

    ```
    + *Escherichia coli* *uidA* gene for beta-glucuronidase mutant C464R, clone M4

    https://www.ncbi.nlm.nih.gov/nuccore/AM180564.1

    ```
    misc_difference 1390
                    /gene="uidA"
                    /experiment="experimental evidence, no additional details
                    recorded"
                    /note="site directed mutagenesis"
                    /compare=S69414.1
                    /allele="wild type"
                    /phenotype="Active GUS"
                    /replace="t"

    ```

See also:
* **Example 3**: https://www.ncbi.nlm.nih.gov/nuccore/KC999035.4

```
 misc_difference 137414..137415

                 /note="TT compared with single T in original ISCfr1
                 sequence that would result in two partial transposase
                 proteins. Confirmed as correct in raw sequence data."

```

Here again the *misc\_difference* feature key is used to describe *how to pass from the represented sequence to a previous sequence version of the considered gene*.

So the *misc\_difference* feature key is appropriate when the sequence in the sequence file *does* match the sequence in the molecule but does ***not** (any more)* match *another reference sequence* (possibly referred to by another gene-related feature) because variations in the sequence have occurred due to some genetic manipulation or, as in Example 3, due to discrepancy with previous sequencing results.

Thus, the *misc\_difference* feature key is appropriate in **Cases 1b and 2c**, and the */replace* qualifier should be used to indicate how to get from the represented sequence to the other reference sequence for the designated feature (in Example 3 above, the */note* qualifier is used instead to describe the difference).

##### Other cases

There are three other cases left:

* **Case 1a**: the features in the GenBank file are consistent with the molecule, but the sequence in the file is wrong
* **Cases 2a and 2b**: the sequence in the GenBank file is consistent with the molecule, but information about some features is wrong

Those three cases can all be classified as a GenBank file containing errors (in contrast with **case 1b** where the "error" is in the obtained molecule and with **case 2c** where the information in the GenBank file is merely incomplete). A GenBank file containing errors is just not appropriate, and the solution here is simply to correct the GenBank file to accurately describe the corresponding molecule.

All that remains is to answer to the second question about how to track changes...

#### Q2. How to keep track of changes between versions of a GenBank file when corrections or updates are made?

##### Versioning and time-stamping GenBank files

GenBank files have two fields that can be used to track their history: the **Modification Date** field and the **VERSION** field.

The **VERSION** field uses an *accession.version* format, where the *accession number* is fixed and the *version number* is incremented only when the sequence data is modified.

See:

*Note that only modifications in the sequence result in a change in the version number. Changes in annotations only affect the timestamp (Modification Date) of the file.*

See:

With this information we can now describe how to properly keep track of the successive versions of a GenBank file in each of the cases described in **Q1**.

##### Proper way to keep track of modifications in GenBank files

**Case 1a**: (wrong) sequence in the file vs (correct) sequence in the molecule

The sequence in the file must be corrected and the version number of the file must be incremented by one.

Moreover the respective file versions (the former and the newer) must refer to each other in the COMMENT field of their respective header.

See:

**Case 1b**: ("correct", theoretical) sequence in the file vs (actual, mutated) sequence in the molecule

A new GenBank file and new molecule name must be created to accurately describe the variant molecule. The *misc\_difference* feature key along with the */replace* qualifier must be used as described above where it is deemed appropriate.

The gene-related features affected should keep the same */gene* qualifier, but an */allele* qualifier might be used to indicate the presence of a difference and the */note* and */phenotype* qualifiers can be used to indicate the reference sequence to which the allele should be compared to and any possible phenotypic consequence of the variation.

The *misc\_difference* feature key can be given the same gene-related qualifier information as the gene where the difference occurs (*/gene*, */allele*, */note*...)

In the special case of variants unintentionally obtained during genetic construction or genetic manipulation, an altered name for the molecule and GenBank file might be used to reflect the fact that the corresponding sequence is different from the theoretically expected sequence, and the *misc\_difference* feature key should be used to indicate any difference that occurred during the construction compared to the initially targetted sequence.

**Cases 2a, 2b or 2c**: wrong or incomplete feature data

The feature annotation in the file must be corrected and the *Modification Date* field must be updated according to the correction date.

### Additional references

**Sequence revision history**

**Examples of revision history of some GenBank files**

### Added for convenience

Here below is reproduced the description and usage of some *feature keys* and *qualifiers* that might be relevant to, or linked with, the question of *how to describe and annotate sequence variants*.

I grouped those definitions there as it might help one to distinguish the different types of variations and divergences and to understand and chose the appropriate term to annotate sequence variants.

#### Selected feature keys

**gene**

```
Feature Key           gene

Definition            region of biological interest identified as a gene
                      and for which a name has been assigned;

Optional qualifiers   /allele="text"
                      /db_xref="<database>:<identifier>"
                      /experiment="[CATEGORY:]text"
                      /function="text"
                      /gene="text"
                      /gene_synonym="text"
                      /inference="[CATEGORY:]TYPE[(same species)][:EVIDENCE_BASIS]"
                      /locus_tag="text" (single token)
                      /map="text"
                      /note="text"
                      /old_locus_tag="text" (single token)
                      /operon="text"
                      /product="text"
                      /pseudo
                      /pseudogene="TYPE"
                      /phenotype="text"
                      /standard_name="text"
                      /trans_splicing

Comment               the gene feature describes the interval of DNA that
                      corresponds to a genetic trait or phenotype; the feature is,
                      by definition, not strictly bound to it's positions at the
                      ends;  it is meant to represent a region where the gene is
                      located.

```

**misc\_difference**

```
Feature Key           misc_difference

Definition            feature sequence is different from that presented
                      in the entry and cannot be described by any other
                      difference key (old_sequence, variation, or modified_base);

Optional qualifiers   /allele="text"
                      /clone="text"
                      /db_xref="<database>:<identifier>"
                      /experiment="[CATEGORY:]text"
                      /gene="text"
                      /gene_synonym="text"
                      /inference="[CATEGORY:]TYPE[(same species)][:EVIDENCE_BASIS]"
                      /locus_tag="text" (single token)
                      /map="text"
                      /note="text"
                      /old_locus_tag="text" (single token)
                      /phenotype="text"
                      /replace="text"
                      /standard_name="text"

Comment               the misc_difference feature key should be used to
                      describe variability that arises as a result of
                      genetic manipulation (e.g. site directed mutagenesis);
                      use /replace="" to annotate deletion, e.g.
                      misc_difference 412..433
                                      /replace=""

```

**modified\_base**

```
Feature Key           modified_base

Definition            the indicated nucleotide is a modified nucleotide and
                      should be substituted for by the indicated molecule
                      (given in the mod_base qualifier value)

Mandatory qualifiers  /mod_base=<modified_base>

Optional qualifiers   /allele="text"
                      /db_xref="<database>:<identifier>"
                      /experiment="[CATEGORY:]text"
                      /frequency="text"
                      /gene="text"
                      /gene_synonym="text"
                      /inference="[CATEGORY:]TYPE[(same species)][:EVIDENCE_BASIS]"
                      /locus_tag="text" (single token)
                      /map="text"
                      /note="text"
                      /old_locus_tag="text" (single token)

Comment               value is limited to the restricted vocabulary for
                      modified base abbreviations;

```

**old\_sequence**

```
Feature Key           old_sequence

Definition            the presented sequence revises a previous version of the
                      sequence at this location;

Mandatory qualifiers  /citation=[number]
                      Or
                      /compare=[accession-number.sequence-version]

Optional qualifiers   /allele="text"
                      /db_xref="<database>:<identifier>"
                      /experiment="[CATEGORY:]text"
                      /gene="text"
                      /gene_synonym="text"
                      /inference="[CATEGORY:]TYPE[(same species)][:EVIDENCE_BASIS]"
                      /locus_tag="text" (single token)
                      /map="text"
                      /note="text"
                      /old_locus_tag="text" (single token)
                      /replace="text"

Comment               /replace="" is used to annotate deletion, e.g.
                      old_sequence 12..15
                      /replace=""
                      NOTE: This feature key is not valid in entries/records
                      created from 15-Oct-2007.

```

**variation**

```
Feature Key           variation

Definition            a related strain contains stable mutations from the same
                      gene (e.g., RFLPs, polymorphisms, etc.) which differ
                      from the presented sequence at this location (and
                      possibly others);

Optional qualifiers   /allele="text"
                      /db_xref="<database>:<identifier>"
                      /experiment="[CATEGORY:]text"
                      /frequency="text"
                      /gene="text"
                      /gene_synonym="text"
                      /inference="[CATEGORY:]TYPE[(same species)][:EVIDENCE_BASIS]"
                      /locus_tag="text" (single token)
                      /map="text"
                      /note="text"
                      /old_locus_tag="text" (single token)
                      /phenotype="text"
                      /product="text"
                      /replace="text"
                      /standard_name="text"

Comment               used to describe alleles, RFLP's,and other naturally
                      occurring mutations and  polymorphisms; variability arising
                      as a result of genetic manipulation (e.g. site directed
                      mutagenesis) should be described with the misc_difference
                      feature; use /replace="" to annotate deletion, e.g.
                      variation   4..5
                                  /replace=""

```

#### Selected qualifiers

**/allele**

```
  Qualifier       /allele=
  Definition      name of the allele for the given gene
  Value format    "text"
  Example         /allele="adh1-1"
  Comment         all gene-related features (exon, CDS etc) for a given
                  gene should share the same /allele qualifier value;
                  the /allele qualifier value must, by definition, be
                  different from the /gene qualifier value; when used with
                  the variation feature key, the allele qualifier value
                  should be that of the variant.

```

**/clone**

```
  Qualifier       /clone=
  Definition      clone from which the sequence was obtained
  Value format    "text"
  Example         /clone="lambda-hIL7.3"
  Comment         not more than one clone should be specified for a given source
                  feature;  to indicate that the sequence was obtained from
                  multiple clones, multiple source features should be given.

```

**/compare**

```
  Qualifier       /compare=
  Definition      Reference details of an existing public INSD entry
                  to which a comparison is made
  Value format    [accession-number.sequence-version]
  Example         /compare=AJ634337.1
  Comment         not to be used for new entries from 14-MAY-2021; This qualifier
                  may be used on the following features: misc_difference, unsure,
                  old_sequence and variation.
                  The feature "old_sequence" must have either a
                  /citation or a /compare qualifier. Multiple /compare
                  qualifiers with different contents are allowed within a
                  single feature.
                  This qualifier is not intended for large-scale annotation
                  of variations, such as SNPs.

```

**/function**

```
  Qualifier       /function=
  Definition      function attributed to a sequence
  Value format    "text"
  Example         function="essential for recognition of cofactor"
  Comment         /function is used when the gene name and/or product name do not
                  convey the function attributable to a sequence.

```

**/gene**

```
  Qualifier       /gene=
  Definition      symbol of the gene corresponding to a sequence region
  Value format    "text"
  Example         /gene="ilvE"

```

**/locus\_tag**

```
  Qualifier       /locus_tag=
  Definition      a submitter-supplied, systematic, stable identifier for a gene
                  and its associated features, used for tracking purposes
  Value Format    "text"(single token)
                  but not "<1-5 letters><5-9 digit integer>[.<integer>]"
  Example         /locus_tag="ABC_0022"
                  /locus_tag="A1C_00001"
  Comment         /locus_tag can be used with any feature that /gene can be used
                  with; identical /locus_tag values may be used within an
                  entry/record, but only if the identical /locus_tag values are
                  associated with the same gene; in all other circumstances the
                  /locus_tag value must be unique within that entry/record.
                  Multiple /locus_tag values are not allowed within one feature for
                  entries created after 15-OCT-2004.
                  If a /locus_tag needs to be re-assigned the /old_locus_tag
                  qualifier should be used to store the old value. The /locus_tag
                  value should not be in a format which resembles INSD accession
                  numbers, accession.version, or /protein_id identifiers.

```

**/map**

```
  Qualifier       /map=
  Definition      genomic map position of feature
  Value format    "text"
  Example         /map="8q12-q13"

```

**/note**

```
  Qualifier       /note=
  Definition      any comment or additional information
  Value format    "text"
  Example         /note="This qualifier is equivalent to a comment."

```

**/operon**

```
  Qualifier       /operon=
  Definition      name of the group of contiguous genes transcribed into a
                  single transcript to which that feature belongs.
  Value format    "text"
  Example         /operon="lac"

```

**/phenotype**

```
  Qualifier       /phenotype=
  Definition      phenotype conferred by the feature, where phenotype is defined as
                  a physical, biochemical or behavioural characteristic or set of
                  characteristics
  Value format    "text"
  Example         /phenotype="erythromycin resistance"

```

**/product**

```
  Qualifier       /product=
  Definition      name of the product associated with the feature, e.g. the mRNA of
                  an mRNA feature, the polypeptide of a CDS, the mature peptide of
                  a mat_peptide, etc.
  Value format    "text"
  Example         /product="trypsinogen" (when qualifier appears in CDS feature)
                  /product="trypsin" (when qualifier appears in mat_peptide
                  feature)
                  /product="XYZ neural-specific transcript" (when qualifier appears
                  in mRNA feature)

```

**/pseudo**

```
  Qualifier       /pseudo
  Definition      indicates that this feature is a non-functional version of the
                  element named by the feature key
  Value format    none
  Example         /pseudo
  Comment         The qualifier /pseudo should be used to describe non-functional
                  genes that are not formally described as pseudogenes, e.g. CDS
                  has no translation due to other reasons than pseudogenisation
                  events.
                  Other reasons may include sequencing or assembly errors.
                  In order to annotate pseudogenes the qualifier /pseudogene= must
                  be used indicating the TYPE which can be taken from the INSDC
                  controlled vocabulary for pseudogenes.

```

**/pseudogene**

```
  Qualifier       /pseudogene=
  Definition      indicates that this feature is a pseudogene of the element
                  named by the feature key
  Value format    "TYPE"
                  where TYPE is one of the following:
                  processed, unprocessed, unitary, allelic, unknown

  Example         /pseudogene="processed"
                  /pseudogene="unprocessed"
                  /pseudogene="unitary"
                  /pseudogene="allelic"
                  /pseudogene="unknown"

  Comment         TYPE is a term taken from the INSDC controlled vocabulary for
                  pseudogenes (http://www.insdc.org/documents/pseudogene-
                  qualifier-vocabulary):

                  processed: the pseudogene has arisen by reverse transcription
                  of a mRNA into cDNA, followed by reintegration into the genome.
                  Therefore, it has lost any intron/exon structure, and it might
                  have a pseudo-polyA-tail.

                  unprocessed: the pseudogene has arisen from a copy of the
                  parent gene by duplication followed by accumulation of random
                  mutations.
                  The changes, compared to their functional homolog, include
                  insertions, deletions, premature stop codons, frameshifts and a
                  higher proportion of non-synonymous versus synonymous
                  substitutions.

                  unitary: the pseudogene has no parent. It is the original gene,
                  which is functional is some species but disrupted in some way
                  (indels, mutation, recombination) in another species or strain.

                  allelic: a (unitary) pseudogene that is stable in the
                  population but importantly it has a functional alternative
                  allele also in the population. i.e., one strain may have the
                  gene, another strain may have the pseudogene.
                  MHC haplotypes have allelic pseudogenes.

                  unknown: the submitter does not know the method of
                  pseudogenisation.

```

**/standard\_name**

```
  Qualifier       /standard_name=
  Definition      accepted standard name for this feature
  Value format    "text"
  Example         /standard_name="dotted"
  Comment         use /standard_name to give full gene name, but use /gene to
                  give gene symbol (in the above example /gene="Dt").

```

---
Tags: file-formats, genbank
---


title: Biopython reads my tree eternally long

text:
# Question
Title: Biopython reads my tree eternally long
I have a nexus tree (1332 taxa) with a lot of additional data. When I tried to read it through `tree = Phylo.read(treepath, "nexus")`, my kernel got eternally loaded. If I abort the process, I get the following message:

```
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Input In [95], in <cell line: 1>()
----> 1 tree = Phylo.read(r"D:\Downloads_hard_disk\Job\NAs_Hepatitis_A\Results\experiments\test_right\test_right_ExpPop_3\test_right_ExpPop_3.biopy.tree", "nexus")

File ~\miniconda3\lib\site-packages\Bio\Phylo\_io.py:60, in read(file, format, **kwargs)
     58 try:
     59     tree_gen = parse(file, format, **kwargs)
---> 60     tree = next(tree_gen)
     61 except StopIteration:
     62     raise ValueError("There are no trees in this file.") from None

File ~\miniconda3\lib\site-packages\Bio\Phylo\_io.py:49, in parse(file, format, **kwargs)
     34 """Parse a file iteratively, and yield each of the trees it contains.
     35
     36 If a file only contains one tree, this still returns an iterable object that
   (...)
     46
     47 """
     48 with File.as_handle(file) as fp:
---> 49     yield from getattr(supported_formats[format], "parse")(fp, **kwargs)

File ~\miniconda3\lib\site-packages\Bio\Phylo\NexusIO.py:40, in parse(handle)
     32 def parse(handle):
     33     """Parse the trees in a Nexus file.
     34
     35     Uses the old Nexus.Trees parser to extract the trees, converts them back to
   (...)
     38     eventually change Nexus to use the new NewickIO parser directly.)
     39     """
---> 40     nex = Nexus.Nexus(handle)
     42     # NB: Once Nexus.Trees is modified to use Tree.Newick objects, do this:
     43     # return iter(nex.trees)
     44     # Until then, convert the Nexus.Trees.Tree object hierarchy:
     45     def node2clade(nxtree, node):

File ~\miniconda3\lib\site-packages\Bio\Nexus\Nexus.py:668, in Nexus.__init__(self, input)
    665 self.options["gapmode"] = "missing"
    667 if input:
--> 668     self.read(input)
    669 else:
    670     self.read(DEFAULTNEXUS)

File ~\miniconda3\lib\site-packages\Bio\Nexus\Nexus.py:718, in Nexus.read(self, input)
    716     break
    717 if title in KNOWN_NEXUS_BLOCKS:
--> 718     self._parse_nexus_block(title, contents)
    719 else:
    720     self._unknown_nexus_block(title, contents)

File ~\miniconda3\lib\site-packages\Bio\Nexus\Nexus.py:759, in Nexus._parse_nexus_block(self, title, contents)
    757 for line in block.commandlines:
    758     try:
--> 759         getattr(self, "_" + line.command)(line.options)
    760     except AttributeError:
    761         raise NexusError("Unknown command: %s " % line.command) from None

File ~\miniconda3\lib\site-packages\Bio\Nexus\Nexus.py:1181, in Nexus._tree(self, options)
   1179     elif special == "W":
   1180         weight = float(value)
-> 1181 tree = Tree(name=name, weight=weight, rooted=rooted, tree=opts.rest().strip())
   1182 # if there's an active translation table, translate
   1183 if self.translate:

File ~\miniconda3\lib\site-packages\Bio\Nexus\Trees.py:82, in Tree.__init__(self, tree, weight, rooted, name, data, values_are_support, max_support)
     80 # there's discrepancy whether newick allows semicolons et the end
     81 tree = tree.rstrip(";")
---> 82 subtree_info, base_info = self._parse(tree)
     83 root.data = self._add_nodedata(root.data, [[], base_info])
     84 self._add_subtree(parent_id=root.id, tree=subtree_info)

File ~\miniconda3\lib\site-packages\Bio\Nexus\Trees.py:123, in Tree._parse(self, tree)
    121 elif not incomment and tree[p] == ")":
    122     plevel -= 1
--> 123 elif tree[p:].startswith(NODECOMMENT_START):
    124     incomment = True
    125 elif incomment and tree[p] == NODECOMMENT_END:

KeyboardInterrupt:

```

I am interested what is the specific reason of this: either the ineffective algorithm or a bug in **Biopython**? Is there a way to read the tree as `Bio.Phylo.Newick.Tree` object? Or should I completely use another package? In the last case, could you recommend effective one?

# Answer
The reason is the RAM has been exceeded for Biopython's capabilities and you are encountering a 'RAM-bottleneck'. I've listed the options from 1 or 3 and recommend option 3c.

The options are use:

1. machine with much greater RAM
2. acceleration coding techniques in Python (such parallelisation)
3. Beast tools written in Java

I assume this calculation is being performed on your laptop because in theory \<1500 isn't much, but agreed there are ~1500 taxa.

**Option 1** I personally don't think you are very far into the bottle-neck. Certainly a Cloud computer (e.g. GCP) could easily perform this calculation and maybe your the headnode of your local cluster. Usually the headnode has the most RAM.

**Option 2** Code acceleration is worth keeping in mind because this achieves \>\> 1 log performance. Parallisation will solve this using an optimal chunk-size. However, large amounts of code can be needed. Code optimisation - Python is very sensitive to this because it uses C-libraries - for example, use `map` and `reduce` and avoid ANY regex in your Python code because IMO its a key weakness of Python. Pandas iterrows is famously slow. If the Biopython library is the problem code efficiency isn't an option (unless you want to rewrite Biopython). Label compression or header recoding techniques work, I'm not sure whether this has been done in Beast output already(?).

**Option 3** The tools for dealing with this which are written and optimised to achieve rapid output are:

1. Tracer \- does lots of analysis on MCMC output. Its really geared towards assessing convergence, but it can do a lot of stuff.
2. TreeAnnotator \- makes a summary tree from the MCMC output you have

3. **TreeStat2** \<- I think this is what you are looking for.
> 1 votes

---
Tags: phylogenetics, python
---