text:
Why is bulk RNA sequencing reflecting AVERAGE expression but not TOTAL expression of all cells?
When I am reading papers that compares bulk RNA sequencing and single-cell RNA sequencing, we often see papers describe bulk RNA seq measures the average cell expression.
For example, in this paper Single-Cell RNA-Seq Technologies and Related Computational Data Analysis

bulk RNA-seq mainly reflects the averaged gene expression across thousands of cells

And this paper A practical guide to single-cell RNA-sequencing for biomedical research and clinical applications

the averaging that occurs in pooling large numbers of cells does not allow detailed assessment of the fundamental biological unit—the cell—or the individual nuclei that package the genome

And this paper Differential expression analyses for single-cell RNA-Seq: old questions on new data

Standard RNA-Seq experiments need millions of cells for sequencing [2,3], and therefore can only get averaged measurements of gene expressions of the cells sequenced.

My question is, when we perform the bulk RNA seq and calculates parameters like TPM, we don't actually divide the total number of mRNA transcripts with the number of cells (I think sometimes we don't even have an accurate number on the number of cells that we used in bulk RNA seq), we simply normalise it with the transcript length and count.
Therefore, if we DID NOT divide the number of transcripts with the total cell number, how are we measuring the average of expression, but not the total expression?
Or in other words, what is the "average expression" referring to? Average to what?

What these statements refer to is that bulk RNA-seq will overlook the different transcriptional cell types, hence the value of single cell RNA-seq in unraveling heterogeneous cell populations/tissues.
I think "total" is also avoided because total RNA has another meaning related to capturing all the RNA transcripts in a cell (e.g. including non polyA transcripts), which is still not the case for most RNA-seq experiments.
Also note that RNA-seq experiments generally (without spikins) are only relative and not giving absolute quantifications of the amount of RNA in a cell or groups of cells.

Thank you for your answer! I understand what you said, but I'm asking about total expression, not total RNA. You said "bulk RNA-seq will overlook the different transcriptional cell types", yes I agree, but then it also means that it is interpreting the total transcriptional actives within the mixture of cells. I'm still not sure why the word "average" is used in these publications. Or in other words, what is the "average expression" referring to? Average to what?

@benson23 As I said, these numbers are generally relative and not absolute so these statements are not implying an actual mathematical calculation of mean expression (or indeed summed total expression). 'Average' also means 'typical' and that is the essence here, but this is semantics

Thank you for explaining again! I agree that these numbers would be relative, however, I do not quite agree with you that the use of the word "averaged" is semantic. Based on the first paper I quoted, it said "averaged gene expression across thousands of cells", in the second paper "the averaging that occurs in pooling large numbers of cells" and in the third paper "only get averaged measurements of gene expressions of the cells sequenced", which to me, certainly implies some mathematical operations...

Moreover, in academia, the words used should be precise. If the meaning here really refers to "typical", I believe these authors would not use such an ambiguous word of "average". If it really doesn't have a mathematical meaning, the phrase "relative expression" or as you said, "typical expression" would be much better. Also, if you try to replace the word "averaged" with "typical" in the examples that I quoted, it'd be kind of weird.

I believe using the average is fair for this purpose. The moment you break thousands or millions of cells you are already averaging relative mRNA content; you are making a soup of it. Keep in mind that when you sequence RNA (RNA-seq), you are not sequencing every single molecule in the sample. In other words, you do not know how many molecules of each mRNA you have.

What you are doing is comparing relative level of certain mRNAs between two samples in a giving sequencing depth. It's like comparing salinity iin two different seas: Dead Sea has salinity of 337 g/kg and Black Sea has salinity of 20 g/kg. You can consider the salinity a single mRNA. You do not need to know the exact size of Dead Sea to say that salt concentration of Dead Sea is ~16-fold higher than Black Sea.


text:
Different ClinVar classifications/SNPs are linked to the same rsid
I am using R to extract ClinVar classification of SNPs using rsid. Here's an example of my workflow:
Patient data:

Sample
Func
Gene
Codon_Change
dbSNPv151_GRCh38
Chr
Start
End
Ref
Obs

Mary
missense_variant
SPINK5
Cga/Gga
rs199757347
chr5
148091214
148091214
C
G

With this rsid rs199757347, two ClinVar classifications were found with different codon changes, one is C > G (Uncertain significance, missense variant) and the other one C > T (Pathogenic, nonsense variant).
With the rentrez package in R, I can extract both ClinVar classifications:
library(XML)
library(rentrez)

extract_from_esummary(
  entrez_summary(db = "clinvar",
                 id = entrez_search(db = "clinvar", term = "rs199757347")$ids,
                 rettype='xml',
                 parsed = F),
  "germline_classification"
  )

               968176
description    "Uncertain significance"
last_evaluated "2022/01/01 00:00"
review_status  "criteria provided, single submitter"
trait_set      data.frame,2
               351517
description    "Pathogenic"
last_evaluated "2019/01/01 00:00"
review_status  "criteria provided, multiple submitters, no conflicts"
trait_set      data.frame,2

However, only the classification "Uncertain significance" is relevant to my sample, how can I incorporate the alternative alleles in querying ClinVar with rsid? Or do you have a better way to query ClinVar with my sample data?

Do you absolutely need to use rsIDs? The problem isn't clinvar here, it is that the rsID describes two variants.

@terdon, Yes I understand, that's why I wish to "incorporate" the info of ALT allele in the search. And based on my sample data, I think the most useful info for querying ClinVar is the rsID. Do you have other recommendations on querying ClinVar if I don't use rsID? Thank you!

I am afraid not. I am spoiled, the company I work for has a really good API and that's what I use, but it isn't a free product so any answer I post using it would only be helpful to our clients :/ However, my guess is that you should be able to get the nt change if you don't restrict the output to germline_classification only, so then you should be able to parse that to get what you need. Alternatively, if you can convert the rsIDs to HGVS notation (i.e. NM_006846.4(SPINK5):c.652C>G in your case), you can query using that.

source: Stack Exchange

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

sort: votes


text:
How to compare across Microarray datasets for a gene in Mice?
How does one go about comparing expression levels of one gene across multiple microarray datasets pulled from different experiments/papers in mice? In my example, I've downloaded the expression levels for my gene of interest from different microarray datasets for different cell populations of interest (specifically healthy hematopoietic cells). How would I pool the datasets together?
Say I have common cell types from multiple datasets such as B-cells, would I be able to just combine them into one data frame and visualize them as one or would I need to normalize some aspect of that?
I'm really confused because I've had different people tell me different things.

If you are after expression of tissue-specific, but otherwise healthy wild-type mice, Bgee is really a great resource for you. They pull a lot of data, curate it and sort it in a single database that easy to browse. They also have lots of tools for comparative genomics (if you would be wondering about expression of an ortholog in a different species).
If you are just after the methodology of pulling expression data together from multiple sources be because you are interested in pulling data from specific papers, you can check their paper about the database. The pipeline is available on GitHub too.

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


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

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

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

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


text:
Using ticks on the x axis or a scale bar when visualisizing a phylogram?
When visuliazing a phylogram is it preferable to show ticks on the x axis or a scale bar?
Given an example tree with branch lengths in newick format (A:0.01,(B:0.02,C:0.05):0.02);, the two options are illustrated in the plot below:

Both will work. The convention is the scale-bar and the rationale is below (which is a bit complicated).
By convention:

the right-hand tree where there is a scale-bar and is the traditional diagram for phylograms.
the left-hand tree scale-bar ("x-axis") is used for molecular clock representations, i.e. not for phylograms, but dendrograms.

In terms of clarity the "x-axis" (or extended scale bar) appears clearer. However, firstly, for a seasoned phylogenecist it makes no difference and simplicity is preferred. Secondly, there's coalescent theory.
The left-hand tree ... the thing is in a molecular clock representation the scale is a reverse axes starting at the right-hand tips. Whether thats time or genetic diversity its a good representation (well for time its essential) because coalescent theory means that there is zero genetic diversity for that taxa at time 0 (the idea is you are extrapolating, or reconstructing, the common ancestor). The representation you presented confuses that convention because 0 starts at the assumed node in the outgroup.
Thus the scale-bar is better.


text:
LD Score Regression Derivation hard to follow
I am trying to understand the derivations from Sullivan et al. (2015) in the Supplementary Material. There, it is mentioned in the first page that the least squares estimate of the j-th SNP effect, considering the polygenicity linear model $φ=Xβ + ε$, is $\hat{β_j} := X^T_j φ/N$ . Normally, the least squares solution of $\hat{\beta}$ is $(X^TX)^{-1}X^T\phi$. By comparing the two equations with each other, this means that, in this case, $X^TX=N I_{M\times M}$. However, LD is assumed, therefore $E(X_{ij}X_{ik})\neq 0$, so the covariance matrix cannot be diagonal. Could someone explain to me how the derivation for the beta estimate was performed, please? What am I missing? I checked that the same thing is mentioned in more than one sources such as here.

I still don't understand what you're asking. If you want to try and rephrase your question I can give it a shot, but I'm not sure where the confusion lies right now.

I edited the question, I hope now it is clearer what is the mathematical issue I am having trouble understanding.

Not entirely clear on what you're asking but the beta estimates in LDSC come from GWAS summary statistics. As I understand it, they aren't "derived" so much as just set equivalent to what they linear model defines them as, which is their additive effects on the phenotype as implied by the model. So for each SNP j, it's effect size (beta) is whatever difference one copy of the allele makes on the phenotype (phi), divided by N, the number of SNPs.

Sorry for being late, but what do you mean by the phrasing "equivalent to what the linear model defines them as"? If only the linear model solution was considered, then the estimated beta-s would be the pseudo-inverse of X times X' times phi, namely the least squares solution.

Accepting your answer, I think I misunderstood the phrasing in the supplementary note, it is indeed only a way of defining a weighted effect on the explained output, nothing to do with least squares solution...


text:
Querying The Cancer Genome Atlas (TCGA) for gene expression
The UCI gene expression cancer RNA-Seq dataset (https://archive.ics.uci.edu/ml/datasets/gene+expression+cancer+RNA-Seq) consists of 5 types of cancers (BRCA, KIRC, LUAD, COAD, PRAD) and gene expression levels of 20,531 genes for 801 samples. The dataset looks like this:

How to query in The Cancer Genome Atlas (TCGA) to create an updated version of this dataset?
I've done the following query.

Is the given query and the UCI dataset equivalent? If not, how to query properly?
P.S.: I'm a CS student; do not have any biology background.

Not quite what I was hoping for the query rather than the raw link ... Anyway I guess I wasn't clear

To start with, your query above is selecting for breast cancer samples (TCGA-BRCA) only. The original dataset seems to contain 4 other cancer types.
The link you provided has a link to https://www.synapse.org/#!Synapse:syn4301332 and by looking at filenames of the 14 files there, I can tell that there are different types of data: miRNA expression, gene expression, methylation, ... And your query is limited to gene expression only. I am not even sure if you can use just a single query to fetch all of the 14 data files, you might need to use multiple queries.
If you are familiar with R, I would recommend TCGABiolinks, it will abstract away a lot of otherwise annoying manual steps. Apparently you can also use Python retrieve TCGA data:https://docs.gdc.cancer.gov/API/Users_Guide/Python_Examples/

I've added an image of the UCI dataset. I just need to perform a query that returns similar kind of gene expression level tabular data.
As the query only returns TCGA-BRCA data only, will replicating the query for other four types create similar (biologically equivalent and meaningful) kind of dataset?

Yes, you can change the cancer type and re-use the query but my main concern is the type of data that this query would return. Your query fetches count data via the HTSeq - Counts filter but in the posted image, there are floats instead of integers (gene counts), this can be normalized data. If I remember correctly, HTSeq-Counts should be returning actual counts, so integers.


text:
Missing ',' in line when Biopython reads a nexus tree
I want to edit a tree that I got from BEAST2 treeannotator in nexus-format.
Usually I use the module Phylo from Biopython for such work but Phylo.read(r"filename.tree", "nexus") gave me the next exception:
---------------------------------------------------------------------------
NexusError                                Traceback (most recent call last)
Input In [29], in <cell line: 1>()
----> 1 Phylo.read(r"filename.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:1144, in Nexus._translate(self, options)
   1142         break
   1143     elif c != ",":
-> 1144         raise NexusError("Missing ',' in line %s." % options)
   1145 except NexusError:
   1146     raise

NexusError: Missing ',' in line 1 AB298157.1_2015_-7.9133750332192605_114.8086828279248, 2 AB298158.1_2007_-8.41698974207…

Using Nexus.read(Nexus(), input=r"filename.tree") gave the same result. Please could anyone help with this? I cannot understand the reason of this error because nexus file looks correct.

Nexus format
A NEXUS tree file is a rigid but very useful format which can be used to store a phylogenetic tree, along with associated metadata, for example, the alignment.
It will always comprise a:

TAXA block

The TAXA block contains information about taxa.
... can comprise a,

DATA block

The DATA block contains the data matrix (e.g. sequence alignment).
... in this context importantly will have a,

TREES block

The TREES block contains trees described via the format, e.g. ((A,B),C);:
Each block will open with a BEGIN statement and terminate with an END statement. Again its versatility to store specific information within a block is valuable property of the format.

NB! I encountered some problems with my following work. So I would like to warn readers about them.
It is allowed in the newick format to surround the label with quotes, & some programmes or scripts add them to those names that have ambiguous characters. But it can lead to exceptions during the following phylogenetic analysis, for instance, in BEAST. I wish you would be careful with this.

That's a complicated output.
Two solutions, the easiest is simply to use Andrew Rambaut's figtree http://tree.bio.ed.ac.uk/software/figtree/ . Figtree is a very, very good tree drawing program and it will work; it works on Mac, Windows and Linux.
You can then export the file from Figure - from memory in numerous formats, e.g. phylip format, and import it into Biopython if that is where it is needed. If you are parsing Beast output, Beast has A LOT of internal tools and you might be attempting to replace those via Biopython. I would avoid it. Also there is the package for processing Beast output (which I forget), I think its called Tracer.
With Biopython try parse, it might help if there are multiple trees:
from Bio import Phylo
trees = Phylo.parse("examples.nex", "nexus")
for tree in trees:
    print(tree.name)

I would avoid the 'r'. I agree it makes no difference, but really you'd need to be using unusual file names for 'r' to be needed. Its possible you are using it to get around \s, i.e. spaces in the file name. You could just think about simply not including \ and space type characters in file names. In regex it is used for certain but that's because its 'inescapable' (no pun intended, but I thought it was funny).

Just to wrap up and clarify, I recommended:

To open the tree in Figtree;
Find the export menu and output as either 'newick' or preferably 'nexus' or both;
Re-run Biopython importing as 'nexus' or 'newick' (if nexus fails).

The OP reported ...

I have just compared them (nexus->FigTree->nexus).

Good so it opened up in FigTree and exported, this was expected.

Fig Tree compose all this actually in newick & stash it in nexus block decor. While in the previous case newick tree contains IDs, & there is a seclude block with ID-name pairs

I missed this bit, "newick stashed in a nexus block". That is potentially really cool because I suspect that is within the TREES block. If that's correct then that is nexus format. This should now open Biopython.
Why the bug?
One of the issues with Nexus is there are two variants regarding how the TAXA block relates to the TREES block. In one format they recode the taxa - I definitely understand why. Recoding is possibly (very likely) be being used in Beast to increase algorithm efficiency - it makes no sense slushing a load of taxa names around a convergence algorithm. Writing all that stuff out is not cool either.
However, when you are looking at a tree that's not what you want at all - efficiency is an irrelevance and readability is important. Thus recoding the taxa labels is definitely not wanted. This is a possible explanation for the read error.
Whatever the reason for the bug I am certain Biopython has a work around because in phylip alignment format this is certainly true. There is a custom module that handles the specific issues the format generates.
It was fun but time to move on.

Could I use Figtree algorithms through command line? If it has to be converted to another format, I would like to do this automatically.

My advice is to get it working, even if point and click, - it will work then compare the input formats to understand what the issue is. You can then script the parsing.

So do you suggest to convert nexus to nexus & compare the difference?

I have just compared them (nexus->FigTree->nexus). Fig Tree compose all this actually in newick & stash it in nexus block decor. While in the previous case newick tree contains IDs, & there is a seclude block with ID-name pairs.

Good, should work. I've rewritten the answer to accommodate this.


text:
How to use GEO_SPHERE?
I have some problems with geosphere package which extends/applies the Beast Bayesian MCMC package for calculating phylogenetic trees. Please, could anyone give any hints about at least one of them?
Firstly, I cannot reproduce the tutorial. I followed the instructions step-by-step but my results are too different (Fig 1) from specified in it. What could I do wrong?
Secondly, it is oft occurred that my mean and standard deviation increase too high (Fig 2). I am aware that this is bad but what is the right way not to allow such situations? Just put an upper border? In this case, what is the reasonable border?
Thirdly, why is a site model of the geodata partition set as the first option in a list of  models for amino acids or nucleotides substitution, & I cannot pick any other (Fig 3 & 4)? What do precision and fast checkbox mean? Is there any way to get under the hood of model to vary more parameters?
Finally, are there any documentation, guides, tutorials &c. about geosphere besides the aforementioned tutorial? I did not find any on its GitHub repository.
Fig 1

Fig 2

Fig 3

Fig 4

My trees:

As it is turned out from the author's words:

This is a normal result for the tutorial.
It is recommended to border the deviation of geomodel to 3 powers.
There is set up a homogenous diffusion model but BEAUTi does not show this correctly. Precision is defy τ as bt^-1 where t – a branch length, b – precision. It can be an arbitrary number that the quotient of it divided by the potential tree height would much lower than 1. Otherwise, the next approximation has some deterioration:

Threshold is an angle (in degrees 0...180) below which there is no contribution to the probability density. The fast check box turns up approximation of arccos of angles close to zero (cos(x) > 0.9). The error is at the most 1e-10. More hidden parameters can be found through XML-editors in prepared files.
By the way, the author, Remco Bouckaert, has released new version recently where we can specify regions for clades explicitly in BEAUTi.
I hope that this information will be useful for anyone like me who is stuck in GEO_SPHERE.

Thanks @Vovin it is good to know the underlying spatial model. I'm aware of  recent upgrades, its part of a very active Bayesian group

This may not be the exact answer but its a pretty good guess and suggests an absence of getters and setters for this particular issue. Thus I accept the error does not look informative.
I know 100% what is wrong with the MCMC tree model you performed on HBV, hepatitis B virus.
The ESS are mostly <<200, that is where they are in red. The MCMC needs to run for longer or a simpler model (probably both). For example, the prior is 4 which is way too low.
In particular sprecision.geo is - and that does not look good for the calculation you are performing because it suggests it (the geo object is not present and looks remarkably like the error you are receiving). Simply this will be written in Java and that value it will see as an object. Whether this is singly due to insufficient runtime, I don't know, but the ESS values are a pretty good guess at diagnosing the situation.
What might be happening is that if the values in red are too low and insufficient to calculate sprecision.geo hence its - and the absence of this value is preventing the diffusionModel.s:geo from being calculated.
Do consider Beast is a hugely impressive but complicated population dynamics Bayesian simulation which is further complicated because you are stacking a further level of parameters onto it.

I have responded to the comments below here:
HBV is known to be a difficult data set in phylogeny, it has a very weird genome and if further iterations are resulting in errors with the basic MCMC something is definitely not right with the core calculation. 'Too much MCMC' isn't a valid error in a Bayesian calculation. I would simply switch data set and try again. The ESS is simply too low for any further analysis and that's why its in red.
One of the photos provided doesn't seem to stack up.

sprecision.geo is "-" because this is a fixed precision parameter (Fig 3) so it does not vary. This was run 1 million iterations as it was specified in the tutorial. As you can check by the link in the question their results are not ideal but certainly better. Moreover, when I run more iterations (10 & 100 millions), I get problems with rate mean and deviation. By the way, the trees were gotten the same topologically & geographically (I am going to apply them to question).

I've responded above

I am trying to reproduce the tutorial and want to distinct the difference between my actions and the author's ones. I made a crack to process the whole sampling (129 sequences) but understand that I miss something about geo clock priors. My own research is dedicated to HBA. I have about a half of a thousand our exclusive sequences and about a thousand GenBank ones. So I need to study this instrument in detail before full-blown using.

The spatial inference depends on and accurate tree.


text:
What is the meaning of "Filtering allele frequency" in the gnomAD database?
Recent gnomAD versions include a "filtering allele frequency" which tells you when a variant can be safely adjudged not to be disease-causing. Unfortunately, I'm having trouble making sense of the definition, which is

"Technically, this is the highest disease-specific maximum credible population AF for which the observed AC is not compatible with pathogenicity. More practically, If the filter allele frequency of a variant is above the maximum credible population AF for a condition of interest, then that variant should be filtered (ie not considered a candidate causative variant). See http://cardiodb.org/allelefrequencyapp/ and Whiffin et al. 2017 for additional information."

I interpret this statement as "If

gnomAD_WG_AF_{population} > gnomAD_WG_FAF95_{population},

then the variant should be filtered."
(Here "FAF95" stands for "filtering allele frequency at 95% confidence".) However, the above inequality seems to hold for every variant in gnomAD, so I'm clearly missing something. On the other hand, perhaps the statement "maximum credible population AF for a condition of interest" means "proportion of people in the population who have the condition," but in that case the use of the term AF doesn't make sense. Does anyone have a better understanding of this FAF95 number than I do?

My understanding is from Whiffin et al. 2017.
Imagine a dominant mendelian disease A with complete penetrance, caused only by one allele, which is currently unknown. In your search for that unknown allele, a "credible" candidate variants is one whose population allele frequency is lower than the prevalence of that disease in the population. In that sense,
$$\text{max. credible. AF}_\text{A} = \text{prevalence}$$
In a more realistic scenario, a similar disease B can be caused by more than one allele. A credible candidate variant would therefore have a lower expected allele frequency, because only some proportion of diseases in the population should be attributable to that variant (the "allelic contribution" of that variant),
$$\text{max. credible AF}_\text{B} = \text{prevalence} \times \text{max. allelic contrib.}$$
Finally, consider a similar disease C with incomplete penetrance. Intuitively, a credible candidate variant would be allowed to have higher population allele frequency than the disease prevalence, due to the presence of unaffected carriers,
$$\text{max. credible AF}_\text{C} = \text{prevalence} \times \text{max. allelic contrib.} \times 1/\text{penetrance}$$
Thus far, we have talked about the maximum credible population allele frequency. Population reference databases like gnomAD only have a sample of the population. Whiffin and colleagues observe that one can set a value for the maximum credible allele frequency, and subsequently compute by modelling a Poisson distribution, the one-tailed 95% confidence interval for the expected allele count in a sample of a given size.
Conversely, they suggest, you could work backwards from an allele count and find the highest possible value for the maximum credible population allele frequency that gives the allele count as the 95% confidence interval. I believe this is the FAF95 recorded in gnomAD, or at least its definition should be similar.
In other words: find a reasonable estimate for $\text{max. credible AF}_\text{C}$, and filter FAF95 against that.


text:
Mismatch between mutant gene and reference gene outside the site of mutation
I am tasked with designing primers for a particular mutant (Target Gene Locus: At1g28490 i.e. SYP61) of Arabidopsis, obtained from Gabi-Kat:

The NCBI reference sequence for the genome goes:
GGCTTTATAAGATTCTCCAGTAGCCCAAAAATAAACTTTCAGCGAATCTCAAAGCCAACCGCAGGTGATGCGTTTTGCTATGTAAGAGAGACGAATCTCACGCCGGAATCTCTTTTACTTCACCGTCGTTTGTTCTCCGTCGCGTCACCTTTCGATTAATTCCTTCACAGGGTAGCAAAATTTTGTTGATTTCGGG$\color{red}{\textrm{[}}$ATGTCTTCAGCTCAAGATCCATTCTACATTGTTAAAGAAGAGATTCAGGATTCTATTGATAAGTTGCAATCTACATTCCACAAATGGGAGCGTATCTCTCCCGACATGGGTGATCAAGCACATGTGGCGAAGGAGCTTGTTGCTACTTGTGGAAGCATTGAGTGGCAGGTA͟G͟A͟T͟G͟A͟G͟C͟T͟G͟G͟A͟A͟A͟A͟A͟G͟C͟G͟A͟T͟T͟A͟C͟T͟G͟T͟T͟G͟C͟A͟G͟C͟A͟A͟A͟A͟G͟A͟T͟C͟C͟T͟T͟C͟ATGGTATGGCATTGATGAAGCTGAGCTTGAAAAACGGAGGAGATGGACTAGTAATGCTAGGACACAGGTGCGGAATGTGAAGTCTGGTGTTTTAGCTGGAAAGGTTAGTTCTGGAGCAGGTCATGCAAGTGAAGTGCGTAGAGAACTAATGAGAATGCCAAACTCGGGTGAAGCAAGTAGATATGATCAGTATGGAGGGAGAGACGATGATGGTTTTGTACAGTCAGAATCAGATAGGCAAATGCTGTTGATAAAGCAACAAGACGAAGAATTGGATGAGCTTAGTAAAAGTGTACAGAGAATTGGAGGAGTGGGGCTCACTATACACGACGAACTCGTTGCACAGGAGAGAATAATAGATGAATTGGATACGGAGATGGACAGCACAAAGAACCGTCTTGAATTTGTTCAGAAAAAAGTGGGAATGGTAATGAAGAAAGCTGGAGCAAAAGGTCAAATGATGATGATATGTTTCTTGCTCGTCTTGTTCATCATCCTATTCGTTCTTGTCTTCTTGACCTAA$\color{red}{\textrm{]}}$AGAACAGAGCCGCGGGATCTTATAGATGGGAGTTTAATGTCTTACATTATCATGTTTGTGGTGAAACAAATGGAATATTTACATAAGGTTGTCTGTTTTTTCTTCTCCTATCCAATCCCACATATAATTTTGAACTTTGGCTGTGTGATCTCTAGAGAACCTGAAATTAAATACTTTCATAGCTTGATCCCATGGCTAAGTTCCT
Gabi-Kat provides the following sequence from the mutant as a confirmatory match with the genome locus:
AGCGGATCTCTCCCGACATGGCTGATCAAGCACATGTGGCGAAGGAGCTTGTTGCTACTTGTGGAAGCATTGAGTGGCAGGTC͟C͟T͟C͟C͟T͟A͟T͟T͟C͟A͟C͟C͟C͟A͟C͟G͟G͟T͟T͟C͟T͟A͟T͟T͟A͟A͟T͟C͟C͟C͟C͟T͟C͟T͟T͟T͟A͟T͟A͟T͟A͟T͟A͟A͟

I have bold-faced the common region to both sequences. What explains the complete mismatch of the trailing part (u͟n͟d͟e͟r͟l͟i͟n͟e͟d͟) of the Gabi-Kat sequence with the ref genome?
It is not the case that the mutation happened there; we are in the coding region (demarcated by red square brackets in the NCBI sequence) while Gabi-Kat's mutation is in an untranslated region.
Update 1:
Running Ensembl Blast on the GabiKat seq returns a 92.9% identicality with SYP61. Their sequence for the gene (taken from TAIR10.1) goes:
AGCGGATCTCTCCCGACATGGCTGATCAAGCACATGTGGCGAAGGAGCTTGTTGCTACTTGTGGAAGCATTGAGTGGCAGGTA͟C͟T͟C͟G͟T͟A͟T͟T͟T͟A͟G͟T͟C͟A͟C͟G͟G͟T͟T͟C͟T͟A͟T͟T͟A͟A͟T͟C͟T͟C͟G͟T͟C͟T͟T͟T͟A͟T͟A͟T͟A͟T͟A͟A͟
Note the difference in the trailing part from NCBI ref; this almost coincides with the Gabi-Kat seq!
So, is the NCBI ref-genome and the data for the particular TAIR genome inaccurate? Maybe, TAIR updated their whole genome set but forgot to update data for individual genes?
Update 2:
My proposed explanation for the anomaly makes no sense. TAIR 7 (2004) has the same sequence. It looks like I am missing out on something more fundamental.

Conclusion The evidence suggests the region you are looking at Gabi-kat appears either 27045 downstream or 192250 upstream of syntaxin SNARE (SYP61), but is definitely not within the gene. It is chance homology (e.g. a pseudo-gene) that appears to occur within syntaxin (SYP61). The intergenic region could be much closer to the SYP61 gene if the same reference genome was used for both analyses Genbank LR782542.1 versus Genbank NC_003070.9 (different reference genomes in my analysis). Hence, the University is probably correct.
What I suggest is the genome of A. thaliana is downloaded and Blast analysis is performed locally. This way differences in the reference genome are removed.
The analysis is as follows:
The target protein in question is involved in vesicle trafficking (here) called SNARE, comprising coiled-coil domains and are highly conserved, for example show 100% conservation between Arabidopsis thaliana, Arabidopsis suecica and   Camelina sativa

MSSAQDPFYIVKEEIQDSIDKLQSTFHKWERISPDMGDQAHVAKELVATCGSIEWQVDELEKAITVAAKDP

If Gabi-kat was a protein its sequence would be ...

RISPDMADQAHVAKELVATCGSIEWQVLLFTHGSINPLFIY

Here's a blast of the whole regions of the Gabi-kat. The region in green is the SNARE proteins ...

So the issue is what is the terminal 3' protein and Blast that, i.e. the final bit that matches no known plant and the locus you flagged as being anomalous ..

This is basically zero result via Blast, the protein with low homology are lower metazoans (animals) and a few bacteria ... a long way from plants. In other words this isn't a known protein of natural origin at least on NCBI.
The locus is literally like no protein on earth.
However ... Blasting the nucleotide ...

Reference Arabidopsis thaliana = syntaxin (SNARE; SYP61)

Gabi-Kat Arabidopsis thaliana top two Blast hits are ...

Chromosome 1 at 10039082 - 10039208, or 9819787 to 981991 or
Both value    1e-42
Both dentity is 93%

So Gabi-Kat is a real part of the Arabidopsis thaliana, Gabi-kat is a non-coding region. Thats the major conclusion.
The question where is the location of SYP61? And that is

Chromosome 1 at 10012037 to 10020311

Described here from NCBI's full report

Thus the mutant is 27045bp downstream or 192250 bp upstream IF the same reference genome is used for the blasting (Genome LR782542.1) and the Genbank report NC_003070.9. They could easily be much closer to each other is the same reference genome was used for both sets of analysis.

To answer the comments

I [the OP] do not follow your [my] blast analysis
My [the OP] blasts gives ... ['hits within the gene']
etc ....

It appears to me there is a psuedo-gene phenomenon occurring which is making you think a mutant is located within a gene. However, if you disagree thats perfectly okay and the original version of this answer I had assumed the Gaba-kat team had got this wrong (I always from the basis of supporting the OP).
NCBI was the database repository in the remote analysis, but there's an issue with multiple references ... so to reiterate I suggest you perform this analysis as a local blast to remove multiple reference phenomena.
To conclude I personally don't think there is a basis to believe the University responsible for Gaba-kat has got this wrong https://www.gabi-kat.de. However, if you feel otherwise then I would recommend contacting them with your analysis and any supportive analysis for example from this thread or elsewhere. Please note, I don't have a connection with this university.

I do not follow your BLAST of the Gabi-Kat Nucleotide Seq. My blast gives a 92.9% alignment for the entire stretch (127 bp long) with SYP61 (1:10016601 to 1:10016721). E = 4.6e-45

Indeed, SYP61 is located from 1:10015856 to 10018047 in NC_003070.9 (GCA_000001735.2; TAIR v. 10.1). Which is inclusive of 1:10016601 to 1:10016721. You have got the specifics wrong.

Running a nucleotide blast on the underlined part of the NCBI ref-seq (that is the part, which diverged from Gabi-Kat seq; 45 bp) gives 100% match to SYP61 (1:10016871 to 1:10016915). E = 3.5e-18.

To conclude, the 45 bp stretch is non-contiguous to the 127 bp stretch in NC_003070.9 (GCA_000001735.2; TAIR v. 10.1).


text:
How to explain disparities in Flavonoid Docking results between Vina and Schrodinger?
I am planning to dock a protein to Selleck Chem's Flavonoid Library. I have designed two workflows using Autodock Vina and Schrodinger as follows:
Schrodinger

Use their proprietary tools for protein preparation followed by grid preparation under default settings.
Run LigPrep on the SelleckChem library to generate one ligand (at most) per compound under default settings.

I executed another run (let's call it 'special run'; see below) generating eight ligands (at most) per compound under default settings.

Run Glide in SP mode under default settings + add a reward for H-bonds.

Autodock Vina

Run every module of UCSF Chimera's Dock Prep on the Protein under default settings.

Except use Reduce in place of AddH - suggested by Chimera as a better alternative.

Calculate grid-center to be the centroid of the default ligand and use the same grid size as used in Schrodinger workflow.
Run prepare_receptor4.py but in non-default settings to only remove lone pairs and non-standard residues.

We do not remove waters or add charges, as done under default settings; they have been already executed by Chimera.

Download 3D files of each compound in the Selleckchem library from Pubchem in sdf format.

Use Obabel to convert them to pdb format.
Run prepare_ligand4.py under default settings to get output in pdbqt format.
Use Obabel to replace Gasteigger charges with eem2015ba.
Use Obabel to minimize energy using MMFF94 forcefield under default settings.

Run Vina 1.2 to dock at exhaustiveness of 32.

The results from the two workflows exhibit vast disparity with the correlation coefficient of the rank being about 0.2. For example, the best result of Schrodinger is ranked 150 by Vina (the library has about 200 molecules in total) and so on. If I take the best 20 results from Vina and Schrodinger each, there is a single common compound.
In my 'special run', I saw that Schrodinger assigns vastly different overall docking scores to various conformations of the same ligand. This might answer the above disparity unless I am doing something acutely wrong. But that raises the question about how to choose the initial conformation of the ligand if I am constrained to use open-source software and starting conformations influence the final result so severely?

Please show the first few lines of results, or whatever other appropriate way can demonstrate this desparity.

@gringer I have added a couple of lines. Does this suffice?

Yes, thanks. I think that's enough.


text:
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 ?

This is a good question. One thing I'm puzzled by are 'some sequence files contain errors' if its an ambiguity that is simply dealt with using the degenerate code, e.g. R for A or G ambiguity, N for total ambiguity. If its formatting errors within the file, I'd personally remove them.

@M__ Thank you for the comment. By "errors" I don't mean "ambiguities", I mean sequence files that are not accurate in the sense that the sequence they contain is different from the sequence carried by the molecule they are supposed to describe as checked by sequencing. One could then just correct or discard the files, but for some of them, for the sake of archiving and keeping track of the developments, I would like to be able to annotate those errors.

@M__ Following up your (temporary ?) answer below: * Examples of errors: (i) a feature is properly annotated, but the sequence given is incorrect, (ii) a feature is correctly mentioned, but the sequence span is incorrect (iii) the name of a feature is incorrect (wrong variant for example). * Examples of "variations": (i) during the cloning process a T mutated to a G (variation = 'mutation') (ii) during the cloning process some bases got inserted/deleted (variation = 'insertion' or 'deletion' or 'indel'). Let me know if more details can help.

@M__ That's true, I can see that... Yes, I think too that the most important aspects are mutations (variations), seconded by errors in the sequence data (over annotation errors). Thank you very much again for looking into this.

okay post-updated, official response below.

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

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:

Escherichia coli uidA gene for beta-glucuronidase mutant Y468C, clone M6
https://www.ncbi.nlm.nih.gov/nuccore/AM180565.1

Escherichia coli uidA gene for beta-glucuronidase mutant W549R, clone M12
https://www.ncbi.nlm.nih.gov/nuccore/AM180566.1

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:

https://www.ncbi.nlm.nih.gov/genbank/samplerecord/
https://www.ncbi.nlm.nih.gov/genbank/samplerecord/#ModificationDateB
https://www.ncbi.nlm.nih.gov/genbank/samplerecord/#VersionB

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:

https://www.ncbi.nlm.nih.gov/genbank/sequenceids/

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:

https://www.ncbi.nlm.nih.gov/genbank/release/current/ --> 3.4.2  Entry Organization

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

https://www.ncbi.nlm.nih.gov/genbank/sequencerevisionhistory/

Examples of revision history of some GenBank files

https://www.ncbi.nlm.nih.gov/nuccore/U46667?report=girevhist
https://www.ncbi.nlm.nih.gov/nuccore/NM_029716?report=girevhist
https://www.ncbi.nlm.nih.gov/nuccore/NM_033116?report=girevhist

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").


text:
Mapping to assemblies in subdirectories
I am performing a bash loop one liner:.
for k in */Assembly/*/*.fastq; do minimap2 -ax map-ont assembly.fasta $k > mapping.sam; echo $k; done

The file hierarchy is as follows:
\Drain_4
    \Assembly
        \A4_assembly
           assembly.fasta
           A4.fastq
        \D3_assembly
        \I8_assembly
        \U7_assembly
\Drain_5
\Drain_7

Where \ denotes a directory. This is the error i am getting:
[ERROR] failed to open file '.fasta': No such file or directory A4.fastq

Both input files are in the same directory.
The fastq is my current directory. I'm guessing i will need to provide the complete pathname for both my input files?
I should mention that in all folders the assembly.fasta file is named the same but the fastq file is different (has a different prefix depending on which folder it is in.

your file hierarchy doesn't indicate any file called NC_D2/NC_D2_trim_filt.fastq; where is this file from?

The most likely issue is that the assembly.fasta file is not found in the current directory. Taking @m's comments into consideration, a minimal rewriting of your script to fix this would be as follows:
for k in */Assembly/*/*.fastq; do minimap2 -ax map-ont $(dirname ${k})/assembly.fasta $k > $(dirname ${k})/mapping.sam; echo $k; done

But I'd recommend at least adding additional quote and brace protection, just in case file names contained spaces:
for k in */Assembly/*/*.fastq;
  fDir = $(dirname ${k});
  do minimap2 -ax map-ont "${fDir}/assembly.fasta" "${k}" > "${fDir}/mapping.sam";
  echo ${k};
done

edit: fixed up variable capitalisation

The solution is quite simple you issue the one liner either from within A4_assembly directory or else NC_D2 (see below). It will then find assembly.fasta and the error will disappear. That's because the current working directory is assigned to the directory in which the one liner is issued. Directory A4_assembly is where assembly.fasta is situated, thus is will work.
If will fail to find the fastq file thus you will need to rewrite the loop, simple,
for k in *.fastq; do minimap2 -ax map-ont assembly.fasta $k > mapping.sam; echo $k; done

... will work, executed from within Assembly_4
The one liner suggests you are executing the command from the directory above the Drain_4 directory.
If this is how you intended the one liner to operate I don't think its the solution to the one liner, because there is a risk you will encounter logical errors (not bugs) ..
*/Assembly/*/*.fastq

If you did execute this commend from above Drain_4 it will looking all over the place here in every Drain and every assembly directory therein. If you executed this anywhere else it will fail to find *.fastq because the start of the wildcards is the current working directory in a bash for relative paths.
However, when you dump the file it goes into one one place:
 > mapping.sam

So there is a risk are performing loads and loads of minimaps but they are simply overwriting each of other, i.e. they are replacing each other.
The other issue is
NC_D2/

This directory is not represented in your directory file.
Fundamentally the problem is that assembly.fasta is always in one location which appears to be in the directory above the directory structure you described, viz. */Assembly/*/*.fastq ... according to the one liner, but in the directory structure its in A4_assembly it is also likely to be in \D3_assembly amongst loads of others and other Drains (master directories).
Thus even if it 'finds' the reference genome (if that is what it is?), it will see that for everything ... and I mean everything.
There are dynamic solutions to what I think the code using bash is looking to achieve because you can issue regex from within the do ... done, or simply do it as a script.

I suggest to simply execute the first script I wrote in every Assembly directory then everything will simply work and there's no risk of logical errors.


text:
How to use EPAM Cloud Pipelines?
I am searching for an easy-to-use GUI management tool to handle my WDL workflows, I found EPAM Cloud Pipelines. Unfortunately, it uses cloud services such as AWS, GCP, etc. but I wanted to use my own backend instead, thus, I attempted to see if it would be possible to configure it to use my backend, however, there seems to be little to no instructions on how to install it. If anyone has any experience using this, can you please provide a description? Or, if you use or know of another GUI tool similar to EPAM Cloud Pipelines that can help manage and build WDL workflows, please share some information on that aswell.


text:
How do I use bash for loops to loop through a WDL Array?
Within the example below a WDL array of strings, test is defined as "a","b","c","d" the value total is a variable that signifies the length of the WDL array. Using the bash for loop, which will iterate from 1 up to length of the array, I seek to print out the values at each indicy in the array and output it to another array. For example when the loop is on 1 the value at indicy 1 within the array will be printed out in this case "a" and stored in the array out, if the loop is on 2 then the value "b" will be printed out and stored in the array out, etc. I tried the below example, however, I get an error when running. Can anyone help me with this problem?
How I tried to do this in WDL:
task A {
    input {
        Array[String] test = ["a","b","c","d"]
        Int total = length(test)
    }

    command {
       for (( c = 1; c <= ${total}; c++ ))
       do
           echo ${test[$c]}
       done
    }
    output {
        Array[String] out = read_lines(stdout())
    }
}

I suspect the issue here is that you are trying to use a bash variable inside of a WDL placeholder. The ${} are interpreted by WDL as places where some value should should be inserted into the command. So when you say echo ${test[$c]}, test[$c] is seen a WDL expression. This isn't a valid WDL expression, causing a parsing error.
The choice to use ${} as placeholder syntax in WDL was a rather inconvenient decision by the WDL team, because it clashes with the bash syntax. There is an alternative syntax available for the command section in WDL 1.0, which is preferred for this very reason. This alternative uses ~{} as placeholder syntax, instead of ${}, and wraps the command in <<< >>> (referred to as the "heredoc-like syntax"), rather than the usual { }. In this heredoc-like syntax the ${} isn't seen as a WDL placeholder and can, thus, be used in bash commands like normal. You can read about it in the WDL SPEC.
Using this syntax your WDL would look something like this:
version 1.0

workflow test {
    call A
}

task A {
    input {
        Array[String] test = ["a","b","c","d"]
        Int total = length(test)
    }

    command <<<
       ARRAY=(~{sep=" " test}) # Load array into bash variable
       for (( c = 0; c < ~{total}; c++ )) # bash array are 0-indexed ;)
       do
           echo ${ARRAY[$c]}
       done
    >>>

    output {
        Array[String] out = read_lines(stdout())
    }
}

Note that if the values in your array are more complex and might include characters with special meaning in bash (like spaces), it would probably be better to wrap the values in single quotes like so:
ARRAY=('~{sep="' '" test}')
And just in case you weren't aware, bash for-loops can also loop over arrays directly:
command <<<
       for c in ~{sep=" " test}
       do
           echo $c
       done
>>>

This is a comprehensive answer, giving a clear insight into WDL and its integration with bash.

Thank you for this. In the first answer you gave, is it possible to load the array into a bash variable i.e. ARRAY=(~{sep=" " test}) within the workflow definition so it would like the bellow example? Can I then take that as an input into my task definition?
workflow test {
Array[String] test = ["a","b","c","d"]
Int total = length(test)
ARRAY=(~{sep=" " test})
call A  {
      inputs:  testarray =  ARRAY
 }

}
task A {
inputs {
Array[String] = testarray
}
}

Currently, WDL doesn't have a way of defining environment variables in WDL itself. You have to do so in the command section. There has been some discussion on introducing such a feature though: https://github.com/openwdl/wdl/discussions/458


text:
What is the contributing factor in single measurements of nanopore sequencing? One base, a k-mer or difference between the base left and entered?
My question is about nanopore sequencing and specifically about the current that is measured by the device in each measurement. The question is:
In each measurement in nanopore sequencing, the change in the current between two sides of the membrane at the sensing region is recorded by the sequencer. Is it a single nucleotide or a k-mer or the difference between the first and last base in the kmer that contributes to one single measurement?
I am citing two papers here that I looked at:

Nanopore sequencing technology, bioinformatics and applications
(Nature 2021)
Three decades of nanopore sequencing (Nature 2016)

Based on my finding in these papers, I realized that all the nucleotide molecules in the k-mer in the sensing region contribute to the measurement, and not only one base, yet I am not sure if it is a good assumption.
I know some models like HMM and NN sequential models can computationally infer the individual bases (for example some basecallers). But my question is more about what is physically measured?
In the second cited paper, it is stated that, in more recent nanopore designs, DNA strand can be fixed such that one individual nucleotide is positioned in the sensing region. Still, the diameter of the sensing region is not exactly one base.
Moreover, since the current change is recorded not the current itself, the contributor could be the difference between the base which just passed from the sensing region and the one base that just entered?
So I restate my question; what is the contributing factor in single measurements of nanopore sequencing?
In case it is the k-mer that is measured, how could we find k? Is it based on the properties of the device and the kit information?

The sequencing channels measure electric current, i.e. the flow of charged particles through the nanopore. The circuits are not measuring the bases directly; they're measuring the way the flow of ions is changed when different things are introduced into the pores.
Older R9 pores had a length that could fit about five bases in the pore at any one time, the newer R10 pores have a double-size pore that fits about ten bases. Any of these bases can potentially influence the flow of current, and that degree of influence is increased at the narrower pinch points of the Nanopore. Inconsistent unwinding of double-stranded DNA (or hybrid DNA/RNA) can also affect the current profile over time, so there may also be some effects on the current that transfer from where the unwinding occurs, 20 or so bases prior to the pore introduction.
The 5/6 kmer models are (or were, with R9) a reasonable sequencing model, because that's the number of bases that would fit in the pore, so those are the bases that have the greatest influence on the current.
I'll explain using a bigger physical model of hydroelectricity:
The mental model that I have in my head of a flow cell is a dam with a lake / reservoir, and holes in the dam that slowly drain the lake, spinning turbines that produce electricity (i.e. current). The interesting bits in the lake are the things going through the holes that aren’t water, but the water flow is the only thing that can be measured: there are no cameras attached to the holes, only current meters attached to the turbines.
It just so happens that those non-water things have a strange enough shape that they interrupt the water flow in large and predictable ways, which means that observing the water flow can be used as a good proxy for observing the things themselves.
If the holes get blocked, nothing can move through the holes, so there is no current flow: sequencing cannot happen.
Eventually that water reservoir will get depleted. If the reservoir is fully drained [ionic balance of the carrier ions in the flow cell is equalised], then it doesn’t matter how clean the holes are; the dam will produce no more electricity [sequencing cannot happen].
As an aside, this model helps me to understand why keeping the pores occupied with long DNA prolongs the sequencing life of the flow cell. The flow of carrier ions is reduced when a pore is actively sequencing, so the carrier ion reservoir drains slower, and can keep sequencing for longer.


text:
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.

I was under the impression that PacBio was one sequence per ZMW. Did you mean two detected consecutive bases?

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

Thank you again! This is by far the best answer and explanation I received. Now the procedure is much clearer to me. After running guppy, now I have the move table, block stride and sampling rate information. Tombo also added some other info, but I think I will move forward with guppy for now.