text:
Why all values become 1 after dcast?
I am reshaping my dataframe as follows.

I want to set the Wavelength column as the column names with the code as follows.
dcast(unique(df1, by=cols), ID + Condition ~ Wavelength, value.var = 'mean')

But all the values became 1s as follows.

I would appreciate any help.
Yanpeng

I think this should be posted on Stack Overflow - you’d get a better answer there anyway. But please make sure to post data as plain text and not as an image, otherwise it’ll likely get closed :) thanks

That comment is right, it's a general programming question, and text would be much easier for people to help with rather than images.
But, since I ran into this not long ago...
Take a look at the message shown in blue just above the output table.  My guess is you have some repeated ID/Condition/Wavelength combinations and it's defaulting to reporting the number of matching rows in the long-format data frame.  You'll need to either make sure you have all unique combinations, or explicitly define what you expect it to do with multiple values going into a single output row/column position (like, median/mean/max/min/whatever).  The unique() call won't do it since varying mean values would still make for unique input rows.
For example compare this:
> reshape2::dcast(data.frame(
+     Wavelength=c(664, 666, 662),
+     Condition="late",
+     ID=c(6221, 6222, 6223),
+     mean=c(106817, 110357, 112998)),
+   ID + Condition ~ Wavelength, value.var = "mean")
    ID Condition    662    664    666
1 6221      late     NA 106817     NA
2 6222      late     NA     NA 110357
3 6223      late 112998     NA     NA

with this:
> reshape2::dcast(data.frame(
+     Wavelength=c(664, 664, 666, 662),
+     Condition="late",
+     ID=c(6221, 6221, 6222, 6223),
+     mean=c(106800, 106817, 110357, 112998)),
+   ID + Condition ~ Wavelength, value.var = "mean")
Aggregation function missing: defaulting to length
    ID Condition 662 664 666
1 6221      late   0   2   0
2 6222      late   0   0   1
3 6223      late   1   0   0


text:
method to compare the similarity between proteomic profiles
I have a problem finding a methodology to compare proteomic profiles against another. The experiment is the following:

Experiment 1. Ten pdx tumour models with their proteomic profiles were done in one Research Institute with DIA-MS. All of them are the same tumour from different patients. Differential analysis was done with Limma.

Experiment 2. Same tumour, but this time entirely on humans (not pdx) and also controls were included (SWATH-MS in another research institute). Differential analysis was done by t-test.

The biological question is which of the 10 pdx models resembles better with the human one.
I tried different ways of normalizing together the raw data but there is a huge difference in the levels of abundance between both experiments. Also, the pdx ones are really similar between them.
In the end, I did a one vs rest contrast on the pdx conditions to get the unique markers for each model and use them as genesets (positive and negative separated) to try to see if the upregulated and downregulated genes (tumour vs control on experiment 2) are some kind of enriched (pre-ranked GSEA with rank as -log10 p-value*log2FC).
Do you think this can be a solution? Is there a more simple way to do this assessment?
Thank you!!

source: Stack Exchange

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

sort: votes


text:
How to represent a summary of disagreement between two secondary structure predictors given predictions for many proteins?
I have 6000 predicted protein structures.
I am comparing secondary structures predicted by an algorithm with the correct secondary structure assignments.
I need to show which positions in a predicted protein chain have the most differences as to the true protein chain.
How can I do that?
How to represent a summary of disagreement between two secondary structure predictors given predictions for many proteins?
I already prepared confusion matrices. I was looking for something graphical like chart or histogram, etc..

Are you asking:
a. How to represent the disagreement between two secondary structure predictors on one protein
or
b. How to represent a summary of disagreement between two secondary structure predictors given predictions for many proteins
or
c. If there is any common way and simple tool for representing a or b? (kind of yes for a, not to my knowledge for b)

@O.Laprevote, I am asking (b).

(This comment got too long)

positions in a predicted protein chain

You do mean position in SS sequence, not coordinates in tertiary, yes?
From the sounds of it you are talking about SS —not tertiary— and want the analysis to be independent of tertiary.
Therefore, you are comparing for each protein two 1D arrays of sequential elements containing between 3 and 6 states, one for each residue. Technically, you could add another state, a chain break, and have two colossally long 1D arrays. This very much becomes a statistics classification topic.
The properties of these arrays are:

The values are sequential.
Neighbouring elements influence each other as a span of residues makes a SS element. Thus a correct classification will propagate. This will mean that all metrics will be hugely inflated.
Assumption: the divergence is likely on the ends of a span. This could be tested and plotted!

The simplest thing would be collapse the states into a binary loop vs non loop and classify each element of the two arrays and do a typical confusion matrix and its metrics (see Wikipedia). This would have the issue that neighbours influence each other.
For a figure, given N states, one could simply plot a N × N heatmap of the true states vs. predicted states. This would show a hot diagonal of the correct allocations, but will show which states are more prone to misallocation. The frequency of the states may need normalising... but this would still be affected by the length of the SS spans. Therefore, a more sophisticated model may be needed for better normalisation.
For more advanced classifications, more advanced algorithms would be needed, but it still remains a stats problem. It might be worth asking the question again under that discussion framework or look up a relevant multi-state classifier.

I already prepared confusion matrices. I was looking for something graphical like chart or histogram, etc..

If so, given N states, wouldn't a simple N×N heatmap of the true states vs. predicted states suffice? This would show a hot diagonal of the correct allocations, but will show which states are more prone to misallocation.

Please, add your last comment to your answer.

As said, the above is a comment, advising to rephrase or repost the whole Q as a statistical problem as there's a lot of stats know-how here, way way more than compbiochem...

One problem still remains: the proteins are not the same in length. So, how can I align them in the heatmap?

The size of the suggested heatmap is N x N, where N is the number of states (alpha-Helix, Loop, shEet, Pi-helix, Turn etc.). The tally would be on the concatentaed SS sequence as you did for the confusion matrix, but without collapsing the states to binary.

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


text:
Confusing result from Sanger Imputation Service (Eagle v2.4 for phasing, PBWT v3.1 for imputation)
I ran my 23andMe data through the Sanger Imputation Service that uses Eagle v2.4 for phasing, PBWT v3.1 for imputation.
However, some aspects of the results are very confusing to me. Perhaps see #3 (the most confusing) first.
Point 1

There are 62,623 rsIDs for chromosome 3 in my 23andMe file and 2,821,894 SNPs (actually, unique rsIDs) after imputation. However, only 60,144 of the original 62,623 rsIDs are in that ~2.8 million. Where did the other 2,479 go? Here are some examples of missing rsIDs:
rs1516332,
rs12374025,
rs1170695,
rs2731343 and
rs63749864.

EDIT Point 1
Looking at rs1516332, the position is 3:105793, which does have a result in the file, but it doesn't have an rsID:
3    105793    .    C    T    .    PASS    RefPanelAF=0.17878;AN=2;AC=0;INFO=1    GT:ADS:DS:GP    0|0:0,0:0:1,0,0

Looking at rs63749864 (3:13564048) it's missing by rsID and position.
EDIT Point 1 (again)
As described here: https://github.com/richarddurbin/pbwt/issues/51

"... of the 2,479 'missing' SNPs, 1,806 can be 'found' [by position] leaving 673 'anchors' [or rsIDs called by 23andMe] missing from chromosome 3."

Point 2

Some imputed rsIDs appear more than once in the results. How should I interpret that? e.g.

3   64071   rs116059446 G   T   .   PASS    RefPanelAF=0.000369572;AN=2;AC=0;INFO=1 GT:ADS:DS:GP    0|0:0,0:0:1,0,0
3   64071   rs116059446 G   A   .   PASS    RefPanelAF=9.2393e-05;AN=2;AC=0;INFO=1  GT:ADS:DS:GP    0|0:0,0:0:1,0,0

EDIT Point 2
From the answer below, this is a variant with more than two alleles: rs116059446
Seems the format of this VCF file only allows one alt per line, hence the two alts are split over two lines.
See my question here: Is it valid VCF not to 'squash' positions with more than one ALT allele?
And the answer from Durbo here: https://github.com/richarddurbin/pbwt/issues/52#issuecomment-1062951105
In the above case, the genotype is easy to interpret (REF/REF in both lines).
However, I'm now finding examples that I just can't interpret:
3   255395  rs331869    C   G   .   PASS    RefPanelAF=0.385864;AN=2;AC=1;INFO=1    GT:ADS:DS:GP    1|0:1,0:1:0,1,0
3   255395  rs331869    C   T   .   PASS    RefPanelAF=0.980998;AN=2;AC=2;INFO=1    GT:ADS:DS:GP    1|1:1,1:2:0,0,1

Is my genotype CG or TT at this position?
Point 3

Sometimes a rsID appears twice but at different locations:

3   4942430 rs71634747  G   A   .   PASS    RefPanelAF=0.445673;AN=2;AC=1;INFO=1    GT:ADS:DS:GP    1|0:1,0:1:0,1,0
3   4942432 rs71634747  C   G   .   PASS    RefPanelAF=0.248907;AN=2;AC=0;INFO=1    GT:ADS:DS:GP    0|0:0,0:0:1,0,0

Unless I'm really confused, this just seems wrong. If it wasn't for this weirdness, I'd seek concrete answers for #1 and #2, but it could be that they are both just a symptom of this potential error?

Oh, vr-impute@sanger.ac.uk isn't responding to my email. I have asked the above of them.

Don't know. These could be retired or inconsistent SNPs (e.g. shifted locations, or otherwise found to be different from expected). It might be that the imputation suggests that the SNP should be something other than what is observed. I can't see anything obvious from the dbSNP information that I've looked at.

a) In this specific case for rs116059446, that's a SNP that has more than two variants. One way of normalising VCF files is to rearrange them so that each line only refer to a single polymorphism from the reference sequence, so if there are two non-reference variants, the variant will get two lines in the VCF file.
b) For rs331869, this is another SNP with more than two variants recorded in dbSNP. Given that you're digging through the weeds of a huge dataset, the way of thinking about this needs to change from "Is this common" to "Is this possible". Problems are much more likely to exist with SNPs that aren't biallelic, because it indicates that something funky is going on with genetics or testing. I'll propose three possibilities for what could be happening here, in what I consider to be descending order of likelihood:

you have one non-reference allele on one of the sister chromosome copies, and another non-reference allele on another copy
you have a somatic mutation (i.e. a mutation that happened during your lifetime), and some of your cells have one variant, and others have a different variant.
the probes are tagging different regions in your genome from what is expected. Probes should be location-unique, but that's impossible to guarantee for everyone, because genetic variation is a thing that happens.

In this specific case for rs71634747, it looks like that SNP has been retired, and has merged with rs724135. It's possible that this is a multi-nucleotide variant where variants are only observed in tandem. This is an intronic variant, so seeing a 3-nucleotide variation is a little unexpected, but given that you're picking out exceptions from a very large dataset, not completely unexpected.

I think the main take-home from this is that you shouldn't ever expect results to be a perfect representation of reality. It can be useful to look at raw data to help understand what's going on, but there's too much weird stuff going on with biology, chemistry and physics for anything to be 100% accurate.
If you did exactly the same test at a later date, some of your 23andMe results will be different. If you did imputation at a later date, with a different program version, or different imputation reference dataset, your results are likely to be different.
Accuracy can be substantially improved by using different technologies to attempt to discover the same thing, and only concentrating on consistently-reported differences. You could try an Infinium Omni2.5 SNPchip, for example, or whole-genome sequencing. Unfortunately, this will lead to another problem: what should be done with the inconsistent results? The answer to that is situation dependent, and is one of the hardest challenges of bioinformatics.

Thanks for these answers. I'm not expecting data to reflect reality (I've been doing bioinformatics for too long for that ;-) I just expect data to be sane. This is just a basic sanity test, flagging things that looked wrong. I think point 1 and 3 are now tidied up (I don't think the software is behaving correctly, but I can at least understand how the software is behaving, and it's not bonkers.

I still don't understand point 2. Lets look at your potential sources for this data sanity check failure... 2b.1) makes perfect sense, but I would expect the genotype to be 0/1 and 0/1. That would be easy to interpret. That's not what I see (1/0 and 1/1). 2b.2) This may be true, but the data could never show that (imputation). 2b.3) Same, not possible in this data (imputed). Increasingly I think this call is just wrong. Not biologically wrong, just fails basic sanity checks.

Like imagine I showed you a vcf that had 3 3333 . A G . . . GT 4/4... Well, you could say, ah ha, this is a sign of a quadruplex DNA at this location in this individual... well maybe. I'd just say the VCF is broken. It's just wrong / impossible to interpret (imho).

The answer is to relax and throw away the 0.001% of stuff that doesn't make sense, lower your standards and just get on with finishing the job you set out to do without actually caring about consistency, logic, meaning or formats that are actually defined or interpretable.
Nobody cares about what these SNPs mean, how to interpret them programmatically, why there are version issues in the output of the imputation server or why it throws away ground truth in the output, leaving the user to struggle with all these issues. I mean... if nobody cares, why should I?
I've written up as detailed a description of the problems as I can here:
https://github.com/richarddurbin/pbwt/issues/52 and here: https://github.com/richarddurbin/pbwt/issues/51
and written a script to throw away all the junk here:
https://github.com/Geromics/covcheck/blob/wip/report-v2/Research/debug_imputation_results.py


text:
How can I convert protein distance maps and sequences to a pdb files?
How can I convert protein distance maps and sequences to pdb files?
For the same problem, in the case of a model predicting the structure of a protein, such as AlphaFold, how do we convert the generated output tensor to a PDB file and visualize it?

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

It would be helpful to see what your outputs look like now. PDB files are supposedly provided as part of alphafold predictions (see github page: https://github.com/deepmind/alphafold#alphafold-output). PDB format described here: https://en.wikipedia.org/wiki/Protein_Data_Bank_(file_format)


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 best way to compare differences in DEG expression between two populations?
I am interested in differential response to a drug in the "RED" group versus the "GRN" group. I mean, I need to compare differences in DEG expression between two populations (RED vs GRN)
I have 3 pairs of control vs treated samples. 3 pairs for each population.
Here is an example of results for a gene named "gen1".
Gname   RED_ctl1  RED_drg1  RED_ctl2  RED_drg2  RED_ctl3  RED_drg3
gen1           4        10         5        10         6         9
Gname   GRN_ctl1  GRN_drg1  GRN_ctl2  GRN_drg2  GRN_ctl3  GRN_drg3
gen1           4         9         5        10         4         9

In the end, I expect to plot a boxplot that would show me (or not) the difference between RED and GRN groups in terms of DEG expression. I expect to find a group (either RED or GRN) with DEGs that respond to drugs much more (or less) in terms of expression.
Since this is my first time doing this analysis, I would like to know what is the classical way of analyzing that kind of data? What are the steps?
A related question was presented here on Bioinformatics SE

A t-test or eyeballing a box plot might be okay for something very informal, but that's not the statistically right way to do it.
If this is based on RNASeq data, the answer is clear.  You need DESeq2 or EdgeR.  Read their tutorials, and pick one to use.
If this is microarray, limma is probably best.
If it's something else, you'll need to read some papers to see what people use.

It doesn't matter what data you are dealing with classical statistics can handle this data and non-parametric are very versatile, but there is one notable limitation.
If this is strictly paired data (which what I suspect):
GRN 1 | GRN 2 | GRN 3
RED 1 | RED 2 | RED 3

This would be a two-way ANOVA and is a very powerful test it is fulfils the criteria (below). Interpreting the output needs a bit of background reading.
However, if it is simply GRN vs RED then it is one-way ANOVA (or T-test) and Kruskal-Wallis test (non-parametric). In addition, if you are comparing
RED1
RED2
RED3
GRN1
GRN2
GRN3

Thus 6 data sets a T-test doesn't work and its one-way ANOVA or KW test.
ANOVA The Fmax test is needed for ANOVA analysis to assess the homogeneity of variance.

Additional ANOVA criteria

The test works best when the number of samples drawn from each population is the same (check in this case). This is not an essential requirement but if the sampling between groups is heavily that's not good.
The underlying populations are normally distributed (same is true of the T-test). If Fmax says 'yes' skip this bit.

I honestly don't know if Fmax is available in Python, it is definitely not numpy.fmax. It is available in R:
hartleyTest(x, ...)

where x is the vector. You can do it manually.
For 1-way ANOVA
from scipy.stats import f_oneway

f_oneway(GRN123, RED123)
f_oneway(GRN1,GRN2,GRN3,RED1,RED2,RED3) # This makes it like the second example

The data is within a list e.g. GRN1 = [1,2,3,1,1,14,1]
Easy.
Two-way ANOVA is complicated
import statsmodels.api as sm
from statsmodels.formula.api import old

model = ols('treatment ~ C(RED1) + C(GRN1) ... : C(GRN1)', data=df).fit() # see note
sm.stats.anova_lm(model, typ=2)

You should be able to circumvent the ols - because its complicated. It is however useful, particularly is Fmax fails.
Here for details.
You might substitute your one-way ANOVA for a T-test in this instance, but it depends how you structure you comparison.
If Fmax fails and you were wanting to perform a one-way ANOVA you shift to the non-parametric Kruskal-Wallis test. In python,
from scipy import stats

stats.kruskal(GRN123, RED123)

Work around if its paired data ...
stats.kruskal (GRN1, RED1)
stats.kruskal (GRN2, RED2)
....

Its a bit cheeky but you can cover yourself with Bonferoni correction and strictly Wilcoxan is better here.

So its easy.
HOWEVER if your data is paired and fails an Fmax test, strictly you go to the Friedman test, this is the non-parametric version of the 2-way ANOVA (by Milton Friedman). Now for some reason no-one ever does this, but  I don't remember what the reason why. Thus your options are:

in that situation you could just use the ols regression model fitted into the ANOVA - but you need to think hard about his model because it needs to present nicely paired data to ANOVA.
You could of course do the Friedman test and its up to the reviewer to understand the theoretical criticism of this (I don't remember it and generally the quantitive standards in molecular peer-review isn't great)
You could use multiple KW

Alternative the 'correct' way to get ANOVA working when Fmax says 'no' is transformation. This gets complicated however and is one step before a custom model.

Final points

Exception If this is microarray data - there is an issue about standardisation and that is complicated (comparing different microarray slides against experimental variation).
What algorithms such as DSeq2 do will likely do is apply a fixed transformation to the data, but I don't know about standardisation. Non-parametric statistics will circumvent the requirement for transformation ... their whole point is you don't transform and must not transform - must be raw data.
Again non-parametrics become difficult in replicating 2-way ANOVA ... there are several approaches here and of course there is nothing to stop you using it and leave to someone else explain the limitation.

Note the Wilcoxon is the non-parametric equivalent of the T-test, which I think we covered in a previous post (these are mice right?).

What would I do? Personally, if a non-parametric fits I use it every time ... its simple, everyone understand it and a significant result is not in dispute. The reality is you are probably into a skewed distribution. Furthermore, no-one gets worried about a few KW tests (type 2 error) and you can always deploy bonferoni's correction.
I should understand the limitation of Friedman test and I might be working under a false assumption. If a simple transformation results in good Fmax I'd use ANOVA (if raw data isn't homogeneous). Beyond its statistical modelling and that's complicated.


text:
Problems when using vcf2maf
I am trying to use a commonly recommended package named vcf2mafto convert vcf file to maf file and is reliant on Ensembl's Variant Effect Predictor (VEP) here. However, if I follow the guidelines shown on the github, I would get:
perl vcf2maf.pl --ref-fasta hg38.fa --input-vcf chr21.vcf --output-maf chr21.maf
dyld[33981]: Library not loaded: @rpath/libcrypto.1.0.0.dylib
  Referenced from: /usr/local/anaconda3/bin/samtools
  Reason: tried: '/usr/local/anaconda3/bin/../lib/libcrypto.1.0.0.dylib' (no such file), '/usr/local/anaconda3/bin/../lib/libcrypto.1.0.0.dylib' (no such file), '/usr/local/lib/libcrypto.1.0.0.dylib' (no such file), '/usr/lib/libcrypto.1.0.0.dylib' (no such file)
dyld[33982]: Library not loaded: @rpath/libcrypto.1.0.0.dylib
  Referenced from: /usr/local/anaconda3/bin/samtools
  Reason: tried: '/usr/local/anaconda3/bin/../lib/libcrypto.1.0.0.dylib' (no such file), '/usr/local/anaconda3/bin/../lib/libcrypto.1.0.0.dylib' (no such file), '/usr/local/lib/libcrypto.1.0.0.dylib' (no such file), '/usr/lib/libcrypto.1.0.0.dylib' (no such file)
dyld[33983]: Library not loaded: @rpath/libcrypto.1.0.0.dylib
  Referenced from: /usr/local/anaconda3/bin/samtools
  Reason: tried: '/usr/local/anaconda3/bin/../lib/libcrypto.1.0.0.dylib' (no such file), '/usr/local/anaconda3/bin/../lib/libcrypto.1.0.0.dylib' (no such file), '/usr/local/lib/libcrypto.1.0.0.dylib' (no such file), '/usr/lib/libcrypto.1.0.0.dylib' (no such file)
dyld[33984]: Library not loaded: @rpath/libcrypto.1.0.0.dylib
  Referenced from: /usr/local/anaconda3/bin/samtools
  Reason: tried: '/usr/local/anaconda3/bin/../lib/libcrypto.1.0.0.dylib' (no such file), '/usr/local/anaconda3/bin/../lib/libcrypto.1.0.0.dylib' (no such file), '/usr/local/lib/libcrypto.1.0.0.dylib' (no such file), '/usr/lib/libcrypto.1.0.0.dylib' (no such file)
ERROR: You're either using an outdated samtools, or --ref-fasta is not the same genome build as your --input-vcf. at vcf2maf.pl line 410.

I think the reference is the correct one, and I guess the error might be because of the VEP installation (as I came across a lot of bugs when following their gist).
As a result, I tried to use the Conda installation and use the codes:
conda install -c bioconda vcf2maf

However, after the package was successfully installed, I could not call it in the command line.
(base) user-MacBook-Pro:~ scottxu$ conda install -c bioconda vcf2maf
Collecting package metadata (current_repodata.json): done
Solving environment: done

# All requested packages already installed.

(base) user-MacBook-Pro:~ user$ vcf2maf
-bash: vcf2maf: command not found

Or, if I want to use the Perl script:
I will get:
(base) user-MacBook-Pro:~ user$ perl vcf2maf.pl
Can't open Perl script "vcf2maf.pl": No such file or directory

Could anyone give me suggestions on how to install and use vcf2maf correctly? Thank you very much in advance~~
I am working on macOS Monterey 12.2.1 with an Intel Core i5 processor.

I've just installed it via conda and it works ... I have answered below . Geee wished I'd installed on a custom environment because the uninstallation is a bit messy

Install it via Anaconda/conda.
conda install -c bioconda vcf2maf
conda install -c bioconda/label/cf201901 vcf2maf

I personally like Anaconda-navigator for OSX but you would likely go for 'command line version' in part because it will place it in your home directory. Its here https://www.anaconda.com/products/individual and just use the either the above and then its done. The conda environments can be considered once vcf2maf is working.
It is never cool to run Perl on the system command line.

Once installed ensure the conda activated
conda activate <env name> # if its base I don't think the <> is needed

Then simply type
vcf2maf.pl

Output
> Usage:
>      perl maf2vcf.pl --help
>      perl maf2vcf.pl --input-maf test.maf --output-dir vcfs

Its a Perl script hence it will not recognise maf2vcf without the suffix pl. Anaconda doesn't need perl <perlscript>.
Oh just to mention I just conda install vcf2maf without the channel but I think I've got bioconda on my repo list.

I've just uninstalled it. There appears to be an error with the help docs. The help docs (probably PODs) may not have been uploaded and in any case it will run without them. What conda gives you is the full set of dependencies for the script (makes uninstallation a headache) -  however you don't need to worry hence I would just use it maf2vcf.pl --input-maf test.maf --output-dir vcfs that should work.

Note of explanation: why Anaconda works
I better explain, the script run by itself is complaining it can't see the dependencies it requires. It doesn't just need them, it needs to see them. Anaconda/conda takes care of all that by downloading the dependencies and ensuring they are all visible. The downside is uninstalling everything afterwards, but that's my problem not the OPs.

Yeah, but how could I use it after it is installed? Running these codes I will only get the outcomes that the package is already installed. However, I could not call vcf2maf in the command line.

Thank you very much for your kind help. I could get the output if I type vcf2maf.pl; however, if I try to use perl vcf2maf.pl --help, I would get an error message: Can't open perl script "vcf2maf.pl": No such file or directory. Do you know why this could happen? Thanks a lot~~~~~

I've just uninstalled .. it was a bit messy.Answered above, anyway its been fun.


text:
Bedtools count split reads spanning whole feature
[Note: this is for bedtools v2.30]
I want to count the number of RNA-Seq reads that fully cover a given splice junction. For that, I thought of defining a BED feature around the junction, then using bedtools coverage.
Here is a simplified situation:
1-BASED POS 1|      10|       20|       30|
CHROMOSOME   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
BED FILE              **
BAM READ FWD >>>>>>>>>>
BAM READ REV                  <<<

We can generate the corresponding example files:
echo -e "I\t9\t11\tmy_position\t.\t+" > single_place.bed

echo -e "@HD\tVN:1.4\tSO:coordinate
@SQ\tSN:I\tLN:30
first_read\t99\tI\t1\t255\t10M\t=\t20\t20\tAATTCCGGAA\tDDDDDDDDDD\tNH:i:1\tHI:i:1\tAS:i:127\tnM:i:0
first_mate\t147\tI\t20\t255\t3M\t=\t1\t20\tATC\tDDD\tNH:i:1\tHI:i:1\tAS:i:127\tnM:i:0" | samtools view -b > fwd_only.bam
samtools index fwd_only.bam

cat single_place.bed
I       9       11      my_position     .       +
samtools view fwd_only.bam
first_read      99      I       1       255     10M     =       20      20      AATTCCGGAA      DDDDDDDDDD      NH:i:1  HI:i:1  AS:i:127        nM:i:0
first_mate      147     I       20      255     3M      =       1       20      ATC     DDD     NH:i:1  HI:i:1  AS:i:127        nM:i:0

The per-base count (with the -d option) gives the expected result:
coverageBed -d -a single_place.bed -b fwd_only.bam
I       9       11      my_position     .       +       1       1
I       9       11      my_position     .       +       2       0

And using a filter -f 1.0 to only keep reads that cover the entire feature indeed removes our read:
coverageBed -a single_place.bed -b fwd_only.bam
I       9       11      my_position     .       +       1       1       2       0.5000000
coverageBed -f 1.0 -a single_place.bed -b fwd_only.bam
I       9       11      my_position     .       +       0       0       2       0.0000000

However, in practice, I have a lot of spliced reads, so I wish to use the -split option. In that case, the filtering doesn't work anymore:
coverageBed -split -f 1.0 -a single_place.bed -b fwd_only.bam
I       9       11      my_position     .       +       1       1       2       0.5000000

even though the per-base count is still unchanged:
coverageBed -d -split -a single_place.bed -b fwd_only.bam
I       9       11      my_position     .       +       1       1
I       9       11      my_position     .       +       2       0

Why does -split interfere with -f 1.0?
With spliced reads
Note, in practice I'm more interested in cases where the reads can be spliced, such as this:
1-BASED POS 1|      10|       20|       30|
CHROMOSOME   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
BED FILE              **
BAM READ FWD >>>>>>>>>>__>>
BAM READ REV                  <<<

where _ is a spliced read. The corresponding BAM can be generated with:
echo -e "@HD\tVN:1.4\tSO:coordinate
@SQ\tSN:I\tLN:30
first_read\t99\tI\t1\t255\t10M2N2M\t=\t20\t20\tAAAAAAAAAACC\tDDDDDDDDDDDD\tNH:i:1\tHI:i:1\tAS:i:127\tnM:i:0
first_mate\t147\tI\t20\t255\t3M\t=\t1\t20\tATC\tDDD\tNH:i:1\tHI:i:1\tAS:i:127\tnM:i:0" | samtools view -b > fwd_only.bam
samtools index fwd_only.bam

In that situation, it seems I need to use -split and -f 1.0. I thought of filtering a posteriori with e.g. awk '$10>0.5, but that doesn't seem to work either: when not using -d, the output of coverageBed does not separate the reads that cover the feature partially vs fully.

It appears this is an issue that has been raised on Github, and that I had unfortunately missed.
This issue happens in bedtools coverage (my examples above), but also in bedtools intersect. So, I believe this is actually the same issue as #928. Indeed, switching to bedtools v2.27 solves my issue described above (with v2.30).
Note that this particular issue likely appeared as a consequence of solving issues #773, #750, #673.
My workaround is to use bedtools intersect INVERTING the bed and bam, so that I just filter in the reads that match the feature, then count the resulting number of rows.
intersectBed -split -bed -wb -F 1.0 \
             -a fwd_only.bam -b single_place.bed | \
  cut -f16 | \
  sort | \
  uniq -c

Not quite sure if this is the issue you're having, but I've encountered a similar situation with coverage estimation when trying to map nanopore cDNA reads to the genome. I have a full workflow for mapping oriented long reads here, but here's a summary of the split issue:

The default options for BEDTools treat sequence deletions (which happen frequently in nanopore reads) as a drop in coverage, which can make exon hunting and coverage calculation more difficult. I submitted a pull request to add an additional ignoreD parameter to the command line to allow cDNA reads with split coverage across introns to ignore deletions when considering coverage; this request has now been incorporated into the main BEDtools repository (as of v2.30.0).

Depending on how your mapping is set up, adding -ignoreDto your command line may allow spliced reads to be properly counted for coverage with no deletion gaps:
$ bedtools genomecov -bga -strand '+' -split -ignoreD -ibam mapped_BC01.sam
chrX    0       7384244 0
chrX    7384244 7384413 1
chrX    7384413 171031299       0

$ bedtools genomecov -bga -strand '+' -split -ignoreD -ibam mapped_BC01_exon.sam
chr11   0       62551499        0
chr11   62551499        62551595        1
chr11   62551595        62552139        0
chr11   62552139        62553212        1
chr11   62553212        122082543       0

I do not think this solves it, as the split is to handle the Ns in the CIGAR string, while these are short reads with no deletion. But by pointing to your PR, you made me find a relevant Github issue that I had somehow missed, thank you! I'll post it as an answer for future reference.


text:
Handling failed blocks within scattered WDL workflow
I'm trying to adapt an existing workflow to make it more robust to failures. Below is a simplified example of the workflow.
workflow main {
   File filelist="list_of_input_files.txt"
   Array[String] files=read_lines(filelist)
   scatter(file in files) {
      call task1 { input: in=file }
      if(defined(task1.out)) { call task2 { input: in=task1.out } }
      if(defined(task2.out)) { call task2 { input: in=task2.out }
   }
}

task task1 {
   File in
   String out1 = basename(in,'.txt') + '_1.txt'
   command { if(...) { ./script1.sh ${in} > ${out1} } }
   output { File? out = out1 }
   runtime { ... }
}

task task2 {
   File in
   String out2 = basename(in,'.txt') + '_2.txt'
   command { if(...) {./script2.sh ${in} > ${out2} } }
   output { File? out = out2 }
   runtime { ... }
}

task task3 {
   File in
   String out3 = basename(in,'.txt') + '_3.txt'
   command { if(...) {./script3.sh ${in} > ${out3} } }
   output { File? out = out3 }
   runtime {...}
}

with my list_of_input_files.txt looking like
firstfile.txt
secondfile.txt
thirdfile.txt
...

Obviously this is oversimplified, but I am trying to catch failed blocks by making the output of each task optional (using File?) and then having an if(defined(taskX.out)) condition to call the next task. My reasoning is that if something fails (represented by the if(...) within the command block of each task), then the out file will never be created in that task and so won't exist as input for the next task and thus if(defined(taskX.out)) will resolve as false and skip the next task (and therefore all the subsequent tasks).
With current testing, the first file in the list satisfies the if(...) conditions of all the tasks and succeeds in giving me the final output file firstfile_1_2_3.txt. However, I am getting a Required file output 'secondfile_1.txt' does not exist' error message because it fails the if(...) condition in task1 and so the output file is never created.
From what I've read in the WDL v1.1 spec, optional File outputs are allowed and are undefined if the files they point to don't exist. However, after much confused googling, I've realised it's probably because I'm using Cromwell, which is stuck with WDL v1.0.
My question is whether this should be working in WDL v1.0 and, if not, is my logic correct and should this (psuedo-)code work with WDL v1.1? Is there a way to achieve what I'm trying to do using pure WDL v1.0?

I stayed clear of WDL because seemed to over complicate things. Can re-write this as a bash script but would need to study the pipeline a bit more.

I think the thing to do is contact Broad Institute directly because this appears to be a version issue pertaining the upgrade from 1.0 to 1.1

Your WDL code looks more like draft-2 WDL than 1.0 or 1.1 WDL. In WDL 1.0/1.1 there is a separate input section in tasks and workflows, and the version must be specified on the first line. In any case, optional file outputs should be supported in all versions.

@M__ Thanks - my confusion is based around the different specs. It seems that the features I'm trying to use (optional files) don't seem to be supported in the spec that Cromwell uses.

@DavyCats Thanks for the comment! I think I've made the mistake to assume that Cromwell incorporates the latest drafts, but it actually uses WDL 1.0 and so these features are not supported. Despite the WDL 1.1 spec being > 1 year old, I can't see any Cromwell implementation based around it yet.

Optional outputs should work in WDL 1.0. Here is an example that runs without issue through cromwell (version 74).
Note that for the second call in the scatter the input file needs to be converted to from an optional to an non optional value. Since this only gets run if this value is defined, this conversion can safely be done using the select_first function (which returns the first defined value in a given list of optionals).
version 1.0

workflow test {
    input {
        File filelist = "files.txt"
    }
    Array[String] files = read_lines(filelist)

    scatter (file in  files) {
         call blah as b1 {input: in=file}
         if (defined(b1.out)) {
             call blah as b2 {input: in=select_first([b1.out])}
         }
    }
}

task blah {
    input {
        File in
    }
    String outPath = basename(in,'.txt') + '_blah.txt'

    command {
        if $(grep -q "yay" ~{in})
          then
             echo yay > ~{outPath}
        fi
    }

    output {
        File? out = outPath
    }
}

files.txt looks like this:
test1.txt
test2.txt

With test1.txt containing yay and test2.txt containing nay.


text:
Potential side effects of replacing read group tags in BAM file
I have a set of BAM files where the read group tags have some (default?) values, i.e.:
@RG ID:RG0  LB:LB0  PU:PU0  SM:SM0

This creates issues in my downstream analyses, where multiple BAM files with the same SM tag are used.
Samtools provides a command to replace the read group tag. However, I am not sure if there are possible side-effects that I should be aware of, and I might need to remap the BAM file after this replacement. I want to replace ID and SM with the sample id, which is unique for each BAM file.
Do I need to remap and/or run additional steps, or replacing the RG tags should be sufficient to update the BAM files in a consistent way?

potential (not exactly?) duplicate of https://bioinformatics.stackexchange.com/q/908/57. Let me know if those answers fix your problem too.

Thanks @KamilSJaron, as indicated above, my question was mainly concerned on potential side effects of using samtools to modify the header - and if remapping was needed. Therefore, even if the question you linked might be interesting, it is not exactly a duplicate of that. Thanks for your help, though.

You do not need to remap the files, replacing the read group information with samtools is sufficient to deal with this. When you update your pipeline, have it use the sample name to construct the read group information. That way you won't run into this problem again.

I assume you are using GATK software tools for downstream analysis. If that's the case, then you will most likely encounter issues due to missing information on read groups. I believe that GATK can't differentiate between read groups if they are all similar.
In all cases, you can find more information about SAM/BAM and related specifications on the 1.3 The header section as they describe what is mandatory.
So, if you have access to FastQ files, you can add them directly during the alignment process with bwa mem:

bwa mem

-R STR    Complete read group header line. ’\t’ can be used in STR and will be converted to a TAB in the output SAM. The read group ID will
be attached to every read in the output. An example is
’@RG\tID:foo\tSM:bar’. [null]

An example command line would be:
$bwa mem \
-M \                 // Mark shorter split hits as secondary (for Picard compatibility).
-t "$cpus" \         // number of threads
"$genome" \          // genome reference
"$fastq_r1" \        // R1
"$fastq_r2" \        // provide this line if you have R2
-R "@RG\tID:1\tLB:lib1\tPL:illumina\tSM:$sampleName\tPU:runBarcode" \
> "$sampleName".sam  // Output

bwa2 mem

Bwa-mem2 is the next version of the bwa-mem algorithm in bwa. It
produces alignment identical to bwa and is ~1.3-3.1x faster depending
on the use-case, dataset and the running machine.

The command line is exactly similar to the command line above, except that you would be using bwa2.

Samtools

As you have mentioned, you can completely overwrite (default mode) the read groups already existing by doing:
samtools addreplacerg \
-r "@RG\tID:1\tLB:lib1\tPL:illumina\tSM:$sampleName\tPU:runBarcode" \
-m overwrite_all \         // mode to replace all
-@ 5 \                     // number of cpus
-O b \                     // output bam
-o "$sampleName"_RG.bam \
"$sampleName".bam

Picard AddOrReplaceReadGroups

java -jar $picard AddOrReplaceReadGroups \     // Also available in GATK tools
I="$sampleName".bam \
OUTPUT="$sampleName"_RG.bam \
RGID=1 \
RGLB=library \
RGPL=illumina \
RGPU=runBarcode \
RGSM=${sampleName} \
SORT_ORDER=coordinate

Thank you for your detailed answer. My question was actually more limited to: do I need to remap or update the bam files after having modified the read group with samtools? Therefore, I decided to use the other answer. I will however award the bounty to both once it is possible for me to do it.

@gc5 you can only award a bounty to a single answer.

Hi @terdon, yes. Although I meant to award a second bounty for the help. I will be awarding it in 23 hours when it is possible for me to do it.


text:
Error in as.character.default(new("BSgenome", pkgname = "BSgenome.Hsapiens.UCSC.hg19", : no method for coercing this S4 class to a vector
I am trying to use the following function in R:
DBR1 <- readDNAStringSet(file = 'DBR1.fasta')
DBR1.string <- as.character(DBR1)
RepeatMatchOutput <- motif_detection(seqName = DBR1, chrs = '', start.position = NA,
                                     end.position = NA, motif = 'GAAGAA', nr.mismatch = 0,
                                     reverse.comp = F, print.status =F)

The error message is:

Error in as.vector(x) : no method for coercing this S4 class to a
vector

I must be missing something.
-- output of sessionInfo():
R version 4.1.2 (2021-11-01
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19042)

Matrix products: default

locale:
[1] LC_COLLATE=English_United Kingdom.1252  LC_CTYPE=English_United Kingdom.1252
[3] LC_MONETARY=English_United Kingdom.1252 LC_NUMERIC=C
[5] LC_TIME=English_United Kingdom.1252

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
 [1] reshape2_1.4.4                    ggplot2_3.3.5                     gquad_2.1-1
 [4] ape_5.6-2                         seqinr_4.2-8                      BiocManager_1.30.16
 [7] stringr_1.4.0                     STRAH_1.0                         BSgenome.Hsapiens.UCSC.hg19_1.4.3
[10] BSgenome_1.62.0                   rtracklayer_1.54.0                GenomicRanges_1.46.1
[13] Biostrings_2.62.0                 GenomeInfoDb_1.30.1               XVector_0.34.0
[16] IRanges_2.28.0                    S4Vectors_0.32.3                  BiocGenerics_0.40.0

loaded via a namespace (and not attached):
 [1] SummarizedExperiment_1.24.0 lattice_0.20-45             colorspace_2.0-3            vctrs_0.3.8
 [5] yaml_2.3.5                  utf8_1.2.2                  XML_3.99-0.8                rlang_1.0.2
 [9] pillar_1.7.0                withr_2.5.0                 glue_1.6.2                  BiocParallel_1.28.3
[13] plyr_1.8.6                  matrixStats_0.61.0          GenomeInfoDbData_1.2.7      lifecycle_1.0.1
[17] zlibbioc_1.40.0             MatrixGenerics_1.6.0        munsell_0.5.0               gtable_0.3.0
[21] restfulr_0.0.13             Biobase_2.54.0              parallel_4.1.2              fansi_1.0.2
[25] Rcpp_1.0.8                  scales_1.1.1                DelayedArray_0.20.0         Rsamtools_2.10.0
[29] rjson_0.2.21                stringi_1.7.6               BiocIO_1.4.0                grid_4.1.2
[33] ade4_1.7-18                 cli_3.2.0                   tools_4.1.2                 bitops_1.0-7
[37] magrittr_2.0.2              RCurl_1.98-1.6              tibble_3.1.6                crayon_1.5.0
[41] pkgconfig_2.0.3             MASS_7.3-54                 ellipsis_0.3.2              Matrix_1.3-4
[45] R6_2.5.1                    GenomicAlignments_1.30.0    nlme_3.1-153                compiler_4.1.2
```

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

what library does motif_detection come from? It looks like it might be this but hard to say from what you're doing. What type is DBR1? Is it what seqName is supposed to be? Is chrs the expected type?


text:
Metagenomics pipeline recommendations for short-read data
de novo metagenomics on viral NGS data is a hot-topic. On this site alone at least 4 specific algorithms have been used to identify multi-strain/multi-species for a given data set, however these do not necessarily involve short-read data, and are listed below.
Questions
What approach is recommended for traditional Illumina multiple infection data both for:

for intraspecific (intra-species)  data; and
for different species co-infecting the patient, but within the same virus family?

Considerations:

Benchmarking using a common standard (for example)
Error rates, not miscalling a SNP but producing mosaic artefacts (an artificially mixed genome)
Suitability of short-read data for intra-species metagenomics.

Programs discussed on Bioinformatics SE:
SAVAGE https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5411778/ Baaijens et al. De novo assembly of viral quasispecies using overlap graphs Genome Research 2017 May; 27: 835–848.
S-aligner (@juanjo75es)
https://www.biorxiv.org/content/biorxiv/early/2021/02/02/2021.02.02.429443.full.pdf Juanjo Bermúdez s-aligner: a greedy algorithm for non-greedy de novo genome assembly 2021 bioRχiv February: 1-7
Rey http://denovoassembler.sourceforge.net Sébastien Boisvert et al Genome Biology. 2012 13:R122, December
@LRJoshi here used the Broad Viral NGS pipeline for de novo assembly.
https://viral-ngs.readthedocs.io/en/latest/assembly.html

These are very small genomes < 30 kb. Quasi-species is a very real phenomenon (loads of mutations), however mixed infections are also well-known. In a mixed infection, producing one genome which is a mix of two genetically distinct viruses is a key concern.

Under my experience, the best ones in completeness and N50 are s-aligner, Megahit and rnaSPAdes, but I admit I didn't test the other ones for different reasons. Broad indeed seems to use internally SPAdes and Trinity. I compared the results from these three with reported results in benchmark sets from the other assemblers and the results were better for these three.
From these three I get the best results with s-aligner but it has some issues: it actually requires some practice using it and it is slightly slower for this kind of data (<30kbp genomes). In the other hand, it usually has lower requirements in memory than other assemblers except maybe Megahit.
Here there are some comparative results (NGA50) using public benchmark sets:

Obviously, I am kinda biased as I am the author of s-aligner but the data doesn't lie.
Getting the best results also may not be the first requirement in some cases, like if you want it fast and with a fast learning curve, and in this case Megahit I think is usually quite good. But in these cases, you also can design a script that will automate s-aligner for your specific needs. You can contact me if you need help.
BTW, s-aligner also gets better results with other kinds of data like bacterial genomes, transcriptomes, or small eukaryotes but always at the cost of requiring more running time or more computational resources (running it in AWS, for example).

Many thanks @juanjo75es. The key take home message thus far is the consistency of the benchmarks used and s-align does resonate strongly on this index. I'm considering this from a conservative point of view. The variable performance of Spade/MetaSpades is weird (everyone is aware of Spades). Run time is not the issue as you've pointed out. I've some reading to do.

Great. Indeed with the last version t actually gets better results than what it's reported in the document and does not require Flye even for larger viruses. If you want to test it at some point, let me know and I will provide some advice and the last version.


text:
Summary statistics version problem
I am learning GWAS study and PRS for predicting disease. I was aligning my sequencing data with GRCH38 reference panel. However, in the PRS analysis step, I found that most of the summary statistics are developed in GRCH37. It was a bit inconvenient for me because I have to convert vcf/snp file into GRCH37 which may end up losing some snps. I was checking this database link.
Is there any database of summary statistics in GRCH38?

It's a good question, could you please consider upvoting/accepting the answer by user438383. It's definitely useful IMO.

This is a common issue. Instead of lifting over your genomes, it might be computationally and conceptually easier just to lift over the summary statistics to GRCh37. Alternatively (and preferably), this is done for you for every GWAS published in the GWAS catalog. Simply click on the "harmonized" folder after clicking on "download" and you can download whatever build aligns with your genotype data.

I am sticking with generating vcf in GRCH37. Thanks anyway.

@ShafayetRahat You'll also need to ensure the scorefile is b37 too

ok, thanks for the tips

You might want to check out this tool which contains some software to liftover your coordinates from one build to another.


text:
Taking too long to download metadata file (estimated size 30MB) from Covid-19 Data Portal
I am working on developing a tool to classify covid variants. So I am trying to download SARS-CoV-2 sequence data along with lineage annotation from Covid-19 Data Portal (I could not find a better source that provides both sequence and lineage annotated data.)
Since it is a very large data (given 4,050,433 samples), my idea is to first download the metadata (.tsv file) with accession ids and respective lineages from the portal, filter them to keep only the variants I am interested in, and pass the accession_ids as a text file that I would like to download to CDP File downloader as suggested in the portal.
However, downloading metadata .tsv file from the portal is taking forever (download speed is in bytes sometimes). I have tried downloading directly as a file, using wget and also using curl. I am not sure if there is a way to download it quicker since it is losing connection in between and restarting the download. Just as a side note, the sequence data is downloading pretty fast at 4+Mbps.
I hope my question is clear. Any suggestions/tips on how can I make this happen is very much appreciated. Thanks in advance.

30MB isn't a lot of data at all, given the issue is the EBI servers and not your institutes provision - you would have noticed this before on other projects if it was an institute issue. It therefore suggests someone else is hogging the EBI capacity and there is no provision to bypass this.
The best solution is contact them directly

ecovid19@ebi.ac.uk

... and state the download size and download speed.
There could be a bot sucking up the database (that is what it sounds like), which would be easy to write given FTP is permitted. They will have the ability to block an IP and NCBI have done this before albeit a long time ago.
Finally filling the user experience form would help EBI which is here. However the issue you described suggests escalation. NCBI are very responsive and I would be surprised if EBI did not respond quickly.
Good luck with the project, sounds interesting.
Generic response below if its just traditional over-demand for the service.

I suspect that the beautiful database EBI has provided (and it is amazing) has not anticipated the level of demand for the database size. They are far from being alone.
Modern web architectures, such as Spotify, dedicate enormous effort to matching demand with service provision, notably using Kubernetes. When the server is overloaded they will simply duplicate their core database onto additional machines, which they will be buying in from AWS/GCP. It can even be done dynamically without intervention. NCBI are certainly heading in that direction.
If the EBI wrote this under Kubernetes (it could be a monolith however), the limiting factor is the budget in computer resource allocation (again you can increase Kubernetes provision dynamically). However, if its a monolith this will not be easy and again ultimately someone has pay, i.e. EBI.
Regardless, the fact you can use FTP means EBI are hosting their own service, GCP doesn't permit this and AWS is almost certainly similar. This would mean they have to physically provided additional in-house machines. I visited EBI some time ago and at that time the then director went to some lengths describing how focused EBI were on matching their in-house storage capacity to exponential rise in sequence data. Sudden popularity of a webserver is not trivial issue if it hasn't been anticipated.

Alternatives
Using the NCBI service COVID service is the alternative if there is an urgent need to get the project moving. You may also wish to use GISAID here.
If you want to use the NCBI service, I provided the detailed instructions here: How to programatically download SARS-CoV-2 fasta from NCBI/Genbank via API?
It is however important to understand the possible underlying causes from EBI. I do recall the then EBI director stating their biggest issue was a user simply switching to Genbank, which from their point of view isn't cool given the amount effort going into service provision.

Thank you for the detailed explanation and alternatives. I have found an alternative to EBI from the answer you have suggested. It was quick to download from NCBI.
I will report/inform EBI about the data speed issues and update their response here.

You are welcome. It is a pity for EBI because I think they try hard, however it is also a learning curve for them

I have gotten a response from EBI saying, "They are aware of these downloading issues with few files and they have allocated my ticket to a colleague and help me resolve it asap"

That's good, fingers crossed. They will be sensitive about you switching to Genbank.


text:
Is it valid VCF not to 'squash' positions with more than one ALT allele?
I'm seeing output from PBWT that looks like this:
3   112841  rs78923776  C   G   .   PASS    RefPanelAF=0.000153988;AN=2;AC=0;INFO=1 GT:ADS:DS:GP    0|0:0,0:0:1,0,0
3   112841  rs78923776  C   T   .   PASS    RefPanelAF=0.078534;AN=2;AC=1;INFO=1    GT:ADS:DS:GP    0|1:0,1:1:0,1,0

The first line says that my genome is CC at this position, but the second line (for the second alt allele) says that my genome is CT at this position. OK, the first line couldn't say this, so it calls me as REF/REF, but this call has to be interpreted in the context of the second call, and can't be taken at face value.
That's why I wonder if this is valid VCF or if it's actually breaking the specification.
The line should be 'squashed' down to the following (I think):
3   112841  rs78923776  C   G,T .   PASS    RefPanelAF=0.000153988,0.078534;AN=3;AC=0,1;INFO=1  GT:ADS:DS:GP    0|2:0,0:0:1,0,0

Many pieces of software can only work with biallelic diploid calls, so they have to split up multiallelic genotypes into separate lines.

Thanks. Is there a standard way to work out the genotype call?

It is a valid multi-allelic entry, and most variant callers produce such entries. However, comparison between two such gets more complicated. With biallelic entries, 2 CHR-POS-REF-ALT entries either match or they don't. With multi-allelic entries, they can match completely, partially or not at all, which means you also need to be able to specify how much of a match is a partial match etc. This is a hindrance in annotation processes that annotate using the CHR-POS-REF-ALT fields.
This also introduces cases where an indel can become a multi-allelic entry with a nearby SNV like so:
chr1    100    GA    G
chr1    101    A    AA

can now be represented as
chr1    100    GA    G,AA

which can get confusing. Operations are always more predictable when VCF has been decomposed and left-aligned.

Good to know this is valid. I'm having trouble interpreting some of the genotypes I'm seeing. Is there a standard way to do this?

Can you give me an example of the difficulty you're facing?

Hi @Ram, there are several genotype call combinations that I don't know how to 'decode'. I've written up a report here:
https://github.com/richarddurbin/pbwt/issues/52#issuecomment-1064125959
The simplest way I can describe it is to say that, if the call is split across two lines (for a tri-allelic variant) the REF part of the call should be the same on both lines, no?

I glanced through your thread there and I think you're in the gray area between biology and bioinformatics. VCF entries assume reference sequence is reality, so yes, a multi-allelic site, when split, will have the same REF bases. Assuming that the mutated site is still diploid, a tri-allelic variant will have GT 1/2 where 1 is ALT-1 and 2-ALT-2 (so no REF base)

Right, I'd put that in VCF as something like ./1 and ./1 (assuming I had to split alternative alleles out into separate bi-allelic lines for some reason).

Yes, that's what decomposers do - they give you 1/. GT entries.


text:
A two color bar plot in R
I have raw read counts of miRNAs in two conditions like below
> head(long1)
             miRNA Read count Condition
1:  hsa-miR-539-3p         33    living
2: hsa-miR-199b-5p        102    living
3: hsa-miR-146a-3p         27    living
4:   hsa-let-7a-3p        186    living
5:    hsa-miR-1275         54    living
6:  hsa-miR-29a-3p       6464    living
>

For each individual miRNAs, I want to show the difference of raw read counts across conditions by
p=ggplot(data = long1, aes(x=long1$miRNA, y=long1$`Read count`)) + geom_boxplot(aes(fill=Condition))
p <- p + guides(fill=guide_legend(title="Condition"))
p <- p + xlab("miRNAs") + ylab("Raw mapped reads") + ggtitle("s vs living with cancer")
p + facet_wrap( ~ miRNA, scales="free")+theme_bw() add_pval(p, pairs = list(c("s", "living")), test='wilcox.test')

Which returns this
[![enter image description here][1]][1]
This plot is not good
How I can a plot like below in which miRNAs are on Y axis and two colours shows the conditions like this

Using Google well is one of the best skillsets you can possibly learn in programming-allied fields. The style of plot you're looking for is a "grouped bar plot". I googled "r grouped bar plot" and found a nice link: https://www.r-graph-gallery.com/48-grouped-barplot-with-ggplot2.html
The key is to use position=dodge alongside color=Condition, x="Read count", y=miRNA )not sure ow you'd quote the name "Read count", might be better off renaming that column to readcount).


text:
Convert Arlequin ("DNA" markers) to Genepop (PGDspider or otherwise)
I have Arlequin output files (*.arp) from fastSimcoal2 that I'm trying to convert to genepop files (to read into adegenet). The Arlequin files are using the "DNA" marker (50 bp long), and each includes ~2,500 samples.
When trying to use the program PGDspider to convert these (either the GUI or the CLI; v2.1.1.5), things proceed smoothly until the 1,000th allele is reached, where the program generates an error:
INFO  07-03 11:18:10,630 (GenepopWriter.java:parseAlignedSampleData:596)  -Allele "TGGCATCAGCTTCTTCTAAGCAGGCCCGCTTTGTCTAGCGACAGATCTCG" converted to "999"
INFO  07-03 11:18:10,630 (GenepopWriter.java:parseAlignedSampleData:596)  -Allele "ATCTTGTGGTCATGGAGCACCCAGAGCCAAACTCCTAAAGAAGACCCCGC" converted to "1000"
ERROR 07-03 11:18:10,630 (GenepopWriter.java:parseAlignedSampleData:611)  -data "1000" in ind "2_375" (population: Sample 2) is too long!!!
INFO  07-03 11:18:10,631 (GenepopWriter.java:parseAlignedSampleData:618)  -In GENEPOP, alleles cannot be coded with more than 3 digits.

Is there a way around this error? I've processed RADseq data with many 1,000s of alleles using adegenet before, so I know these types of files can be written into genpop--I just don't know how to specify this to PGDspider. There's a field in the .spid file that I might be able to use to fix this (see below), but I don't know what to populate that field with, and there aren't any suggestions from the (limited) documentation.
# Specify the locus/locus combination you want to write to the GENEPOP file:
GENEPOP_WRITER_LOCUS_COMBINATION_QUESTION=

Alternatively, if there are suggestions for running this conversion without using PGDspider, I'm all ears.
Thanks for your help!

Spooky bug, it looks like a simple bug in the code. However, the author my have been concerned with RAM capacity and hence the ceiling but a more explicit error code should be generated. You can split the data set into two (more?) lots of 1 to 999 and pipe these through PGDspider.  You could write  a new PGDspider or ask the author for a fix. However 1000 inputs isn't much data.

Given this is just above the threshold I would simply break the file into 3? parts.

I'm hesitant to break the files apart, only because the error is a result of PGDspider not separating alleles by sample. During the conversion PGDspider writes the alleles as 3 digit numbers, but it isn't restarting the numbering for new individuals in the Arelquin file. Hence, allele 625 in individual 1 is 625, and allele 1 in individual 2 is 626. I don't know how to fix that behavior, though...

Right okay its a bug, you would have to go into the code. If it was counting in a higher base that could circumvent the issue. You could write an external database which collected (appended) the output before the clock resets and that's probably the easiest way to do it.

Ultimately, I achieved this using the R package strataG. The steps taken (in R) are shown below.
arp.object <- arlequinRead(arp.path)
gtype.object <- arp2gtypes(arp.object)
genind.object <- gtypes2genind(gtype.object)

I found that fastSimcoal2 shows an odd behavior when simulating DNA markers, in which sequences are concatenated together (regardless of marker blocks specified; relevant post here). Running the fastSimcoal2 simulations from strataG addressed this bug (and is more straightforward/reproducible, anyways).


text:
BioPython internal_coords module returns different dihedral angles for the (seemingly) same protein structure
I am using the internal_coords module from BioPython in order to compare dihedral angles of two different conformations of the same protein. The conformations are processed by me in order to make sure that their residues are aligned and they always contain exactly the same atoms, so that in the end I expect to have the same number of dihedral angles in the two structures.
However, despite the match in the residues and atoms, I obtain a different number of dihedral angles for the two conformations.
You can find a reproducible code sample here.
from Bio.PDB import *

def extract_torsion_angles(structure):
    ic_chain = internal_coords.IC_Chain(structure)
    ic_chain.atom_to_internal_coordinates()

    d = ic_chain.dihedra

    """
    Update the keys of a dihedra dictionary to make searching easier.

    The keys of the dihedra dictionary (as obtained from
    extract_torsion_angles()) are 4-tuples where each element is
    an AtomKey containing 5 fields (respos, icode, resname, atm,
    altloc, occ). We use the residue position, residue name and
    atom name in order to create a new representation for each
    AtomKey. The final dihedra dictionary will have 4-tuple keys
    where each element of the tuple is one such string.
    """

    new_d = {}
    for key, value in sorted(d.items(), key=lambda x: x[0]):
        new_key = []
        for atom_key in key:
            # Keep: residue position + resiude name + atom name.
            atom_key_str = (atom_key.akl[0] + '_' +
                atom_key.akl[2] + '_' + atom_key.akl[3])

            new_key.append(atom_key_str)

        new_d[tuple(new_key)] = value

    return new_d

f1 = 'protein1.pdb'
f2 = 'protein2.pdb'

parser = PDBParser(PERMISSIVE=1, QUIET=0)
s1 = parser.get_structure('', f1)
s2 = parser.get_structure('', f2)

c1 = list(s1.get_chains())
c2 = list(s2.get_chains())

print('Chains: ', len(c1), ' ', len(c2))

r1 = list(c1[0].get_residues())
r2 = list(c2[0].get_residues())

print('Residues: ', len(r1), ' ', len(r2))

for i in range(len(r1)):
    resname1 = r1[i].get_resname()
    resname2 = r2[i].get_resname()

    if resname1 != resname2:
        print('Residues not matching')
        break

    atoms1 = list(r1[i].get_atoms())
    atoms2 = list(r2[i].get_atoms())

    atom_names1 = sorted([atom.get_id() for atom in atoms1])
    atom_names2 = sorted([atom.get_id() for atom in atoms2])

    if atom_names1 != atom_names2:
        print('Different atoms!')

d1 = extract_torsion_angles(s1)
d2 = extract_torsion_angles(s2)

print('Torsion angles: ', len(d1), ' ', len(d2))

diff_keys = set(d2.keys()) - set(d1.keys())
print(diff_keys)

These are the two PDB files protein1.pdb and protein2.pdb used.
As you can see, the number of chains, the number of residues, the residue names and the atom names in each residue are the same, yet there are 4 extra dihedral angles in the second protein.
My question is, am I missing something from my comparison? Is it unreasonable to expect the same dihedral angles for these two structures? (And if yes, why?). Is this the expected behaviour of the atom_to_internal_coordinates() function?
And extra: do you have any recommendations for other tools that can easily compute all dihedral angles (including the ones in side chains) for a protein?

By different you mean wildly different? Do the residues look like different when you superpose the protein in PyMOL (align) and look at the residues by eye —PyMOL does have the command get_dihedral, but the laziest way would be to click the four atoms (in edit mode) and it will show the dihedral. Do the values differ?

I should also say, the two conformers should have some slight differences. Especially, surface sidechains and especially if the two protein are crystallised under different conditions...

@MatteoFerla Thanks for your input! The problem in this case was not that the angles are different (which was indeed expected as the proteins were coming from different structures), but rather that BioPython just didn't output certain angles at all for one of the proteins.

The problem got solved after posting on the official BioPython repo.
Because of the processing (where some residues were removed in order to align the proteins), some consecutive residues were too far away from each other and thus the bond length was too large. The solution was to increase the maximum peptide bond allowed, as follows:
internal_coords.IC_Chain.MaxPeptideBond = 5.0
ic_chain = internal_coords.IC_Chain(structure, verbose=True)

Thanks for posting back @cubehead. If you mark the post as accepted that would be cool

Just out of curiosity: were these naughty bonds Asp-Pro or Glu-Pro in loops (self-cleavage and isoglutamate)? Or the more classic case of bad density fitting of loops?

OK , visualizing both proteins in PyMOL , using set cartoon_gap_cutoff = 0  ,
I get :

Apparently the two conformations , have more than one difference in various loops.
I modify the code above , thanks to your hint :
internal_coords.IC_Chain.MaxPeptideBond = 5.0
ic_chain = internal_coords.IC_Chain(structure, verbose=True)

My code :
import Bio

print('########## Biopython  VERSION ##########################################')
print('\n                BIOPYTHON_VERSION : ', Bio.__version__)
print('\n########################################################################\n\n')

from Bio.PDB import *

def extract_torsion_angles(structure):

    print('\nSTRUCTURE ---> ' , structure)

    ic_chain = internal_coords.IC_Chain(structure, verbose=True)

    ic_chain.atom_to_internal_coordinates()

    d = ic_chain.dihedra

    """
    Update the keys of a dihedra dictionary to make searching easier.

    The keys of the dihedra dictionary (as obtained from
    extract_torsion_angles()) are 4-tuples where each element is
    an AtomKey containing 5 fields (respos, icode, resname, atm,
    altloc, occ). We use the residue position, residue name and
    atom name in order to create a new representation for each
    AtomKey. The final dihedra dictionary will have 4-tuple keys
    where each element of the tuple is one such string.
    """

    new_d = {}
    for key, value in sorted(d.items(), key=lambda x: x[0]):
        new_key = []
        for atom_key in key:
            # Keep: residue position + resiude name + atom name.
            atom_key_str = (atom_key.akl[0] + '_' +
                atom_key.akl[2] + '_' + atom_key.akl[3])

            new_key.append(atom_key_str)

        # new_d[tuple(new_key)] = value

        new_d[tuple(new_key)] = tuple((value , structure.get_full_id()))

    return new_d

f1 = '18760-BIOINFORMATICS-protein1.pdb'
f2 = '18760-BIOINFORMATICS-protein2.pdb'

parser = PDBParser(PERMISSIVE=1, QUIET=0)
s1 = parser.get_structure('protein1', f1)
s2 = parser.get_structure('protein2', f2)

c1 = list(s1.get_chains())
c2 = list(s2.get_chains())

print('Chains: ', len(c1), ' ', len(c2))

r1 = list(c1[0].get_residues())
r2 = list(c2[0].get_residues())

print('Residues: ', len(r1), ' ', len(r2))

for i in range(len(r1)):
    resname1 = r1[i].get_resname()
    resname2 = r2[i].get_resname()

    if resname1 != resname2:
        print('Residues not matching')
        break

    atoms1 = list(r1[i].get_atoms())
    atoms2 = list(r2[i].get_atoms())

    atom_names1 = sorted([atom.get_id() for atom in atoms1])
    atom_names2 = sorted([atom.get_id() for atom in atoms2])

    if atom_names1 != atom_names2:
        print('Different atoms!')

def get_no_tors_diff(structure1 , structure2) :

    pept_bond = 1.4

    while pept_bond < 10.0 :

        pept_bond = round(pept_bond , 2)

        print('\n\n Peptide Bond Lenght = ' , pept_bond)

        internal_coords.IC_Chain.MaxPeptideBond = pept_bond

        d1 = extract_torsion_angles(s1)
        d2 = extract_torsion_angles(s2)

        print('Torsion angles: ', len(d1), ' ', len(d2))

        diff_keys = sorted((set(d1.keys()) ^ set(d2.keys())) , key = lambda tup : tup[0].split('_')[0])

        print('diff_keys lenght : ', len(diff_keys))

        print(diff_keys , type(diff_keys))

        for i in diff_keys :

            print(i , i[0].split('_')[0] , type(i))

        diff_keys_values = {}

        for i in diff_keys :

            if i in d1.keys() :

                diff_keys_values[i] = d1[i]

            if i in d2.keys()  :

                diff_keys_values[i] = d2[i]

            else:

                pass

        print('diff_keys_values lenght : ', len(diff_keys_values))

        for i in diff_keys_values :

            print(i , diff_keys_values[i][1])

        print('\n\n')

        if len(diff_keys) == 0 :

            print('\n\n Peptide Bond Lenght without differences = ' , pept_bond)

            return pept_bond

        pept_bond += 0.1

    print('\n\n Cannot get rid of differences up to 10.0A in Peptide Bond ')

    return None

print('get_no_tors_diff(structure1 , structure2)' , get_no_tors_diff(s1 , s2))

First and last bit of output :
Chains:  1   1
Residues:  313   313

 Peptide Bond Lenght =  1.4

STRUCTURE --->  <Structure id=protein1>
chain break at GLY  6  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at LEU  37  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at ARG  38  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at LEU  75  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at LEU  84  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at VAL  92  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at PHE  160  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at PRO  163  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at VAL  164  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at HIS  198  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at PHE  224  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at PHE  245  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at LYS  246  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at LEU  251  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at LYS  279  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at ALA  285  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at ARG  296  due to MaxPeptideBond (1.4 angstroms) exceeded

STRUCTURE --->  <Structure id=protein2>
chain break at PHE  1  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at GLY  6  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at LEU  37  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at ARG  38  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at LEU  75  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at LEU  84  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at VAL  92  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at PHE  160  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at PRO  163  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at VAL  164  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at HIS  198  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at PHE  224  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at PHE  245  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at LYS  246  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at LEU  251  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at ALA  278  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at ALA  285  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at ARG  296  due to MaxPeptideBond (1.4 angstroms) exceeded
Torsion angles:  2504   2501
diff_keys lenght :  9
[('0_N_C', '1_F_N', '1_F_CA', '1_F_C'), ('0_N_CA', '0_N_C', '1_F_N', '1_F_CA'), ('0_N_N', '0_N_CA', '0_N_C', '1_F_N'), ('277_D_C', '278_A_N', '278_A_CA', '278_A_C'), ('277_D_N', '277_D_CA', '277_D_C', '278_A_N'), ('277_D_CA', '277_D_C', '278_A_N', '278_A_CA'), ('278_A_N', '278_A_CA', '278_A_C', '279_K_N'), ('278_A_C', '279_K_N', '279_K_CA', '279_K_C'), ('278_A_CA', '278_A_C', '279_K_N', '279_K_CA')] <class 'list'>
('0_N_C', '1_F_N', '1_F_CA', '1_F_C') 0 <class 'tuple'>
('0_N_CA', '0_N_C', '1_F_N', '1_F_CA') 0 <class 'tuple'>
('0_N_N', '0_N_CA', '0_N_C', '1_F_N') 0 <class 'tuple'>
('277_D_C', '278_A_N', '278_A_CA', '278_A_C') 277 <class 'tuple'>
('277_D_N', '277_D_CA', '277_D_C', '278_A_N') 277 <class 'tuple'>
('277_D_CA', '277_D_C', '278_A_N', '278_A_CA') 277 <class 'tuple'>
('278_A_N', '278_A_CA', '278_A_C', '279_K_N') 278 <class 'tuple'>
('278_A_C', '279_K_N', '279_K_CA', '279_K_C') 278 <class 'tuple'>
('278_A_CA', '278_A_C', '279_K_N', '279_K_CA') 278 <class 'tuple'>
diff_keys_values lenght :  9
('0_N_C', '1_F_N', '1_F_CA', '1_F_C') ('protein1',)
('0_N_CA', '0_N_C', '1_F_N', '1_F_CA') ('protein1',)
('0_N_N', '0_N_CA', '0_N_C', '1_F_N') ('protein1',)
('277_D_C', '278_A_N', '278_A_CA', '278_A_C') ('protein1',)
('277_D_N', '277_D_CA', '277_D_C', '278_A_N') ('protein1',)
('277_D_CA', '277_D_C', '278_A_N', '278_A_CA') ('protein1',)
('278_A_N', '278_A_CA', '278_A_C', '279_K_N') ('protein2',)
('278_A_C', '279_K_N', '279_K_CA', '279_K_C') ('protein2',)
('278_A_CA', '278_A_C', '279_K_N', '279_K_CA') ('protein2',)

.........................
.................
..........
......
 Peptide Bond Lenght =  4.6

STRUCTURE --->  <Structure id=protein1>
chain break at PHE  160  due to MaxPeptideBond (4.6 angstroms) exceeded
chain break at PHE  245  due to MaxPeptideBond (4.6 angstroms) exceeded
chain break at ARG  296  due to MaxPeptideBond (4.6 angstroms) exceeded

STRUCTURE --->  <Structure id=protein2>
chain break at LEU  84  due to MaxPeptideBond (4.6 angstroms) exceeded
chain break at PHE  160  due to MaxPeptideBond (4.6 angstroms) exceeded
chain break at PHE  245  due to MaxPeptideBond (4.6 angstroms) exceeded
chain break at ARG  296  due to MaxPeptideBond (4.6 angstroms) exceeded
Torsion angles:  2546   2543
diff_keys lenght :  3
[('83_K_CA', '83_K_C', '84_L_N', '84_L_CA'), ('83_K_N', '83_K_CA', '83_K_C', '84_L_N'), ('83_K_C', '84_L_N', '84_L_CA', '84_L_C')] <class 'list'>
('83_K_CA', '83_K_C', '84_L_N', '84_L_CA') 83 <class 'tuple'>
('83_K_N', '83_K_CA', '83_K_C', '84_L_N') 83 <class 'tuple'>
('83_K_C', '84_L_N', '84_L_CA', '84_L_C') 83 <class 'tuple'>
diff_keys_values lenght :  3
('83_K_CA', '83_K_C', '84_L_N', '84_L_CA') ('protein1',)
('83_K_N', '83_K_CA', '83_K_C', '84_L_N') ('protein1',)
('83_K_C', '84_L_N', '84_L_CA', '84_L_C') ('protein1',)

 Peptide Bond Lenght =  4.7

STRUCTURE --->  <Structure id=protein1>
chain break at PHE  160  due to MaxPeptideBond (4.7 angstroms) exceeded
chain break at PHE  245  due to MaxPeptideBond (4.7 angstroms) exceeded
chain break at ARG  296  due to MaxPeptideBond (4.7 angstroms) exceeded

STRUCTURE --->  <Structure id=protein2>
chain break at PHE  160  due to MaxPeptideBond (4.7 angstroms) exceeded
chain break at PHE  245  due to MaxPeptideBond (4.7 angstroms) exceeded
chain break at ARG  296  due to MaxPeptideBond (4.7 angstroms) exceeded
Torsion angles:  2546   2546
diff_keys lenght :  0
[] <class 'list'>
diff_keys_values lenght :  0

 Peptide Bond Lenght without differences =  4.7
get_no_tors_diff(structure1 , structure2) 4.7

Using internal_coords.IC_Chain.MaxPeptideBond = 4.7 both models have same Chain Breaks at resi 160, 245 e 296 and are described by 2546 dihedral angles.
My question now is how, if feasible , to force PyMOL to paint petide bonds of such lenght both in cartoons and lines/sticks mode.
Second question if Biopython has a tool to detect chain breaks in Structures other than just printing it out with ic_chain = internal_coords.IC_Chain(structure, verbose=True) , asking this because setting internal_coords.IC_Chain.MaxPeptideBond = 50 gets rid of all the chainbreaks but don't have a way to monitor their disappearing. Pics below show chain breaks at 4.7.


text:
Why does GATK produce both 0/1 and 1/0 genotypes in the same file? Are the two not equivalent?
I have always thought that 1/0 and 0/1 in VCF genotype fields are equivalent. And yet, GATK uses both. For example, these are two variants called in the same sample and the same run of GATK 4.1.4.0:
chr7    117120317   .   ATTCATTGTTTTGAAAGAAAGATGGAAGAATGAACTGAAG    A   748.97  .   AC=1;AF=0.5;AN=2;DP=64;ExcessHet=3.0103;FS=0;MLEAC=1;MLEAF=0.5;MQ=60;QD=11.89;SOR=7.223 GT:AD:DP:GQ:PL:SB   1/0:0,36:63:99:2294,1042,933:0,0,0,36
chr7    117120306   .   GACGCGCT    G   516.6   .   AC=1;AF=0.5;AN=2;BaseQRankSum=0.826;DP=70;ExcessHet=3.0103;FS=0;MLEAC=1;MLEAF=0.5;MQ=60;MQRankSum=0;QD=8.2;ReadPosRankSum=-5.786;SOR=0.976  GT:AD:DP:GQ:PL:SB   0/1:37,26:63:99:524,0,1203:0,20,0,26

And here are the same lines redacted for clarity:
chr7    117120317   .   AT..AG  A   [...]   GT:AD:DP:GQ:PL:SB   1/0:0,36:63:99:2294,1042,933:0,0,0,36
chr7    117120306   .   GA..CT  G   [...]   GT:AD:DP:GQ:PL:SB   0/1:37,26:63:99:524,0,1203:0,20,0,26

As you can see, one is genotyped as 1/0 and the other as 0/1. The VCF specs say:

GT (String): Genotype, encoded as allele values separated by either of / or |. The allele values are 0 for the
reference allele (what is in the REF field), 1 for the first allele listed in ALT, 2 for the second allele list in
ALT and so on. For diploid calls examples could be 0/1, 1 | 0, or 1/2, etc. Haploid calls, e.g. on Y, male
non-pseudoautosomal X, or mitochondrion, are indicated by having only one allele value. A triploid call might
look like 0/0/1. If a call cannot be made for a sample at a given locus, ‘.’ must be specified for each missing
allele in the GT field (for example ‘./.’ for a diploid genotype and ‘.’ for haploid genotype). The meanings of
the separators are as follows (see the PS field below for more details on incorporating phasing information into
the genotypes):

/ : genotype unphased
| : genotype phased

This suggests that although 1|0 is not the same as 0|1, the unphased 0/1 and 1/0 genotypes are completely equivalent. Which makes perfect sense given that the N/N genotypes are unphased, so n/N must be equivalent to N/n.
So, my questions are: why does gatk produce some lines with 0/1 and others with 1/0? Is it a bug? Am I missing something? I can't imagine a developer writing code to pick either 0/1 or 1/0 randomly, so presumably they indicate some difference in what gatk sees about a variant. So are 0/1 and 1/0 not 100% equivalent somehow?

 I have seen How is the GT field in a VCF file defined? which suggests the two are equivalent, but this seems to be contradicted by GATK's own output which is why I am asking this again to get a clear and preferably referenced answer.

@M__ this was from a human sample. But you can't know dominant or recessive unless you also know what gene, if any, is affected by the variant and what that gene's mode of inheritance is with respect to a specific disease, so I don't think it can be about AD vs AR.

Okay finally got it. The convention is it should be in ascending ordering, i.e. 0/1. The question is, why the convention broken?

I just realized that's almost certainly on me. I just took a couple of variants from a file randomly, and I hadn't really noticed the positions. I'm pretty sure they were in the right (ascending) order in the original file. Which lends more credence to your suggestion, @M__.

Convention
The convention is it should be in ascending ordering, i.e. 0/1. The question is, why the convention broken?
Phasing I think, at a guess, the answer is it depends what reference genome you are using. I would suggest that this result would not occur if the reference genome is mum (or dad) and the alleles are occurring in the child. The output supports this idea because it is marked with '/', i.e. they are not phased, thus the reference genome is neither mum nor dad.
Example provided The thing I would like to point out is in the example given 117120306 and 117120317 are not randomly scattered across the genome but appear to be 11 bp apart. For a human genome those are very strongly linked.
Linkage hypothesis Thus I would suggest the answer is linkage in this example, given the reference genome isn't a parent. The mutations are occurring on different strands therefore better represented as 0/1 and then 1/0, i.e. one strand comprises 'A' allele and the opposite strand is the 'G' allele.
Why is unphased linkage important In this scenario the interpretation could be hugely important because it would be stating that each allele has come from different parents. This would be more important than maintaining convention. In an inherited disorder there is a risk the child has a dominant phenotype, when each parent could have been heterozygous. Clearly it depends whether the mutation is forcing an amino acid mutation, at 11bp it is certainly possible, i.e. mutation in the first codon upstream and mutation in the second codon downstream.
More generally in any population analysis, linkage is an enormously important issue and central to any calculation. The point is simply that most analysis infers independent segregation of each allele and in this example that remains a good assumption because the allele (I am suggesting) are on opposite strands. If they were on the same strand it would create a bias away from panmixia.
The hypothesis is formed on this example alone and if the alleles were 10 megabp apart it would be impossible to conclude they were linked, therefore the convention of 0/1 should hold. That's my guess, its not a bad guess.

This is an interesting suggestion. I just picked two random variants from the file, so I will have to investigate a bit to see if what you suggest is actually the case more generally. I'd love to see this backed up by a gatk source (post, or dev, or source file or whatever). Certainly plausible, anyway, thanks!

Hope it helps @terdon. Its a hypothesis to work on. Geraldine@Broad will know of course.

I'm honored by your confidence in me, @M__ — But in this case I don't know for sure  I don't recall seeing something like this before. And if it was a phasing thing I would expect there to be phase related annotations. Are you seeing this in more than one record? And is this single or multi sample, have you done any subsetting, etc?
FWIW you might get a definitive answer faster on the GATK forum.

Alas, but regardless my confidence remains @Geraldine_VdAuwera  Thanks.

Thanks, @Geraldine_VdAuwera, I'll try there.


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:
Detect mutation context in a read of a sam file
After sequencing (Illumina) of some DNA, I generated a sam file through alignment of a fastq file (using Bowtie2). Instead of using variant calling programs, I want to know the specific variances, compared to the reference genome, in each specific read.
For example, I'm only interested in C to A mutations. I would like to include reads that contain MORE than 5 C to A mutations. Afterwards I would like to analyze the specific sequence context of that mutation and for doing that I need to know the 3 nucleotides before and after the C to A mutation, which gives a context of NNNCNNN.
I looked at samtools mpileup or bcftools, but I can't find the exact thing that I want to analyze. Any ideas?
Thank you!

That kind of custom fiddling with reads and variants is very cumbersome, non-standard and also error-prone. Do a standard variant calling pipeline and then filter for the mutations that you want. Then extract the variant position (so the coordinates) and get the variant context from the reference genome. Using individual reads for this is not advisable as only a properly called variant is reliable. The power of NGS is the validation of an observation by redundancy (=many reads covering a variant). With individual reads your variant can simply be noise, there is hardly a way to tell. Alignment will confidently map the read to the "context" as you call it, so in my opinion there is not really an advantage of calling the "context" directly from individual reads. You could better extract every confident C>A mutation, then get the three bases up/downstream and then make some kind of position frequency matrix to create a motif.
Example, just for illustration, imagine these were for variants (A in the middle) with three bases up/downstream:
CTAATTG
CTCATAG
CTTATAG
CCTAGTG

The motif would then be this, you can easily see whether there is some conservation in the "context" of your variant.ad

The variant calling you could do with an end-to-end pipeline such as nf-core sarek, this will eventually give you a VCF file(s) with the variant coordinates. From there simply filter for the variants you want, get the coordinate, add three bases up/downstream to start/end and then use something like bedtools getfasta to get the sequences from the reference genome.

If you want to look at a few, specific reads manually, you can do this on VarSome. We have recently released a new feature that allows you to enter single NGS or Sanger reads and get the list of variants detected in that read. For example, here:

In this example, there was one variant identified when comparing to the human hg38 genome assembly:

CYP24A1(NM_000782.5):c.1186C>T which is at chr20:54158136 G⇒A (note that the CYP24A1 gene is on the - strand which is why the C to T variant in HGVS notation is a G to A in the genome). You can then click on the variant link(s) to investigate the variant's effect and annotation.

Obviously, this isn't practical for many reads, but it can be handy if you just want to look at a few manually. More importantly, NGS is very error prone (which is why we sequence at high depth), and you should absolutely not assume that just because a variant is present in a single read it was also actually present in the source DNA you sequenced.

Disclaimer: I work for the company that builds VarSome.


text:
BPGA error: there was no protein data in file
I'm trying to do pan-genome analysis using BPGA and I face this error:

This error didn't show up in the other data set where all the .gbk is downloaded from NCBI. This .gbk file is our own data and has not been published in NCBI yet. There's of course protein data in it, there's this CDS qualifier "translation", but somehow BPGA didn't detect this.
The difference between our data and NCBI data is that most of the NCBI .gbk files are annotated with PGAP. On the other hand, our .gbk is annotated with Prokka (have tried it in DFAST and KBase). NCBI .gbk files that are annotated with PGAP have this CDS qualifier "protein_id" in them. Besides that, I don't think there's any difference.
Does the BPGA script counts the protein number based on the "protein_id" qualifier? How should I solve this? I can't annotate the .fasta file with PGAP because it's so resource-intensive. I believe that, if Prokka is installed locally, there's this "--compliant" option that can integrate the "protein_id" qualifier into the output .gbk file. I have tried to install Prokka locally, but it didn't run well too, I encountered another error for that and all the solutions didn't work for me.
Any opinion on this issue? Thanks for reading my question.
EDIT:
I tried to insert one dummy "protein_id" on my Prokka-annotated .gbk files and it run perfectly without errors.

But when I tried to insert dummy protein_id on my DFAST-annotated .gbk files, it won't run and faced another error. What is this "protein_id" really defined by the script?

@raysteven,
If you follow the documnetation of BPGA at https://sourceforge.net/p/bpgatool/wiki/Home/, you can find the example of SUPPORTED INPUTS for all types of input files including gbk or gbff.
Unfortunately, Prokka generated files are not supported due to different or missing tags than NCBI.
One simple solution you may try is to edit your genBank file to include the tags such as /protein_id="Any_unique_protein_id_inside_quotes" and everything else shown in the gbk example.
I hope this helps.