text:
Are the conclusions in "The proximal origin of SARS-CoV-2" legit?
I don't have any background in genetics and bioinformatics, so I ask you if you think that the arguments provided in the article The proximal origin of SARS-CoV-2 by Andersen et al. are convincing. In particular:

While the analyses above suggest that SARS-CoV-2 may bind human ACE2 with high affinity, computational analyses predict that the interaction is not ideal and that the RBD sequence is different from those shown in SARS-CoV to be optimal for receptor binding. Thus, the high-affinity binding of the SARS-CoV-2 spike protein to human ACE2 is most likely the result of natural selection...

I interpret this passage as "if the virus had been engineered, they would have done a better job". Is this interpretation correct? And if it is, does it sound convincing to you as a proof that the virus developed in nature?
The second motivation given in the article is the following:

Furthermore, if genetic manipulation had been performed, one of the several reverse-genetic systems available for betacoronaviruses would probably have been used

Again, do you think that this exclude the possibility of a human intervention in the creation of the virus?
Thank you for your opinions.
Edit: I am a fan of Occam razor, I know that the scenario in which the virus originated in the wild is by far more likely than the human engineered scenario. I just want to know if, given our current knowledge in genetics, would have been possible for some high skilled researchers to engineered COVID19 (without implying that they did this with any bad purpose)?
Edit 2: I share the concerns of the author of this post, that is, that the findings of the article are merely opinions and some of arguments are misleading. Anyone having a solid technical background can comment on this? China owns Nature magazine’s ass – Debunking “The proximal origin of SARS-CoV-2” claiming COVID-19 definitely wasn’t from a lab

Hi @Lorenzo , this is on-topic because of the ligand binding modelling, however could you add citations?

Yes, the title of the article is in the question, you can find it here https://www.nature.com/articles/s41591-020-0820-9?fbclid=IwAR3KyrFty1XSXDsNrYcFfchr8ijK6qZeCria_y5VsAj-Oxdiei--fLOiviY

Please read @Michael's excellent answer here: Is it possible for coronavirus or SARS to be synthetic?. I think that might even be essentially a duplicate of your question. And yes, if we were to actually try and make this, we could do a better job. But that is only one of many reasons why it isn't reasonable to believe that this virus has been manufactured, as explained in the answer I link to.

Does this answer your question? Is it possible for coronavirus or SARS to be synthetic?

@Michael I would like to know if the content of the article excludes possible conspiracy theories with more valid arguments than "human intervention is conspiracy theory". I really want to understand the strenght of the arguments of the author

@terdon, no it doesn't but you did address my questions. "And yes, if we were to actually try and make this, we could do a better job". I assume that if someone can do a "better" job, it is true that the virus could have been created in a lab.

@Lorenzo created? No. But yes, it is possible to modify existing viruses, and the point made in the paragraph you quote is that if this had been modified, if it had been engineered from another virus, then it would have been more efficient. But the main reasons we can safely conclude this is not a man-made virus are nicely outlined in Michael's answer. If your question is generally about whether "man-made" viruses are possible and what, exactly, would "man-made" even mean in this context, then please as a new one without referencing SARS-Cov-19.

sars2 was the result of serial passage with gmo mouse has human ace2 and human immune sys model

In summary, the authors are saying the complete opposite of "human intervention".

While the analyses above suggest that SARS-CoV-2 may bind human ACE2
  with high affinity, computational analyses predict that the
  interaction is not ideal and that the RBD sequence is different from
  those shown in SARS-CoV to be optimal for receptor binding.

The interaction with ACE2 is taking place on the Spike protein (S) in both COVID-19 and SARS-Cov. The S protein is very different between COVID-19 and SARS both in terms of sequence and COVID-19 includes an additional furin cleavage site, meaning structurally it is quite distinct. The "computational analysis predicts" is difficult because there is no crystal structure for COVID-19 S protein therefore modelling ligand binding is essentially guess work, particularly as it comprises an additional (furin) cleaved domain. What the authors appear to conclude is that they cannot map ACE2-Spike binding of COVID-19 using an established SARS S protein structure. I assume they are homology modelling SARS-CoV onto COVID-19 S protein sequence data.  In my personal opinion there is no way to homology model an additional furin protease cleavage site (furin is a host protease) so it would be difficult for "good fit" to occur, but it is just my opinion and the authors may have established methods to overcome this notable limitation, e.g. similar phenomena occur in influenza.

Thus, the high-affinity binding of the SARS-CoV-2 spike protein to
  human ACE2 is most likely the result of natural selection..

Their conclusion of "natural selection" is not strictly accurate and should be "Darwinian positive selection", but its a small issue. What they are saying is Darwinian adaptation has occurred as opposed to "purifying selection, i.e. conservation. Thus they are definitely excluding human intervention because that is not "natural selection"

Furthermore, if genetic manipulation had been performed, one of the
  several reverse-genetic systems available for betacoronaviruses would
  probably have been used

They appear to be saying that if a synthetic virus was constructed it would be created using an established reverse genetics system, and all those produced to date look nothing like COVID-19.
The phylogenetics argument is always a stronger argument than either of the above (although my opinion isn't entirely objective).  Phylogeneticists have had a long history in opposing conspiracy theories, e.g. polio vaccination, HIV and the list goes on, so it is much more natural territory. An example of such an argument is here which @terdon helpfully pointed out.

To address the questions of @Hans
Question 1

Purifying selection = any evolutionary change is deleterious, so the virus is less able to transmit between humans. The vast majority of mutations are deleterious in classical Darwinian thinking*
Adaptation = the amino acid change occurs to ensure the virus is better able to transmit between humans. In phylogenetics this is called 'positive selection' which is detected at a nucleotide level, this is however a very stringest test.

*, There is something called nearly neutral theory, but thats just complicated

Question 2
What we are talking about is a virus infectious clone system. For a 30kb virus this would not be trival at all, but I agree for 10kb viruses it much easier. However, the phylogenetics are opposed to this, because SARS-CoV-2 is 5% divergent from a virus (RaTG13) isolated a long time before the current epidemic. You can't just engineer a virus that is 5% divergent from a bat virus across the entire genome and never been seen before. It would be a work of genius the likes we have never seen using bioinformatics as yet unheard of, because what would happen is you would encounter endless deleterious mutations. Even if you got over all those, how do you know what fitness traits you are aiming for?
Serial passage
This level of change is doable by serial passage of the virus, but this leads to attenuation of virus (less dangerous) not increased fitness. The yellow fever virus vaccine was created in this way. Even if a virus had been serially passaged we could tell, because we know the mutational patterns of other serially passaged viruses.

Thank you @Michael. Meanwhile, I've found the following article, strongly opposing to the conclusions of the paper. I share the view of the author. What do you think about it? (And thank you for your time) https://harvardtothebighouse.com/2020/03/19/china-owns-nature-magazines-ass-debunking-the-proximal-origin-of-sars-cov-2-claiming-covid-19-wasnt-from-a-lab/

@Lorenzo That article is pure, uninformed conspiracy theory codswallop. It presents zero evidence and presents highly implausible thought experiments as likely scenarios.

@Konrad ok, but throwing insults won't help anyone to understand. What do you think about what the author of the article says about the possibility that the virus was the result of induced natural selection "...scientists passed the H5N1 Bird Flu through a series of ferret hosts until it gained ACE2 affinity and then became incredibly virulent". I stress out POSSIBILITY, I'm not talking about "strong evidence".

@Lorenzo You need to understand that there was no insult in my comment, but an entirely accurate, if somewhat exasperated, description of the contents of that website. It’s entirely possible to productively and intelligently discuss the plausibility (or lack thereof) of SARS-nCoV-2 being artificially engineered by that blog post cannot form part of such a discussion. Apart from that I have nothing to add to terdon’s comments, and this answer, which incidentally is written by an expert in the field (in contrast to you, me, and the authors of that blog post).

... well someone looking to publish in this field (not an expert). I am backed by real experts (super experts). What I can say is that the spike protein is particularly weird in comparison to the RNA viruses I know, but thats not 'conspiracy theory', its just important yet difficult to understand.

Thank you @Michael, and good work!

Could you please expand in details "What they are saying is Darwinian adaptation has occurred as opposed to 'purifying' selection, i.e. conservation."? How do you deduce "Thus they are definitely excluding human intervention because that is not 'natural selection' "?
Why would not using an established reverse genetics system necessarily exclude genetic manipulation? Is this not premised on laziness or uncreativity of all possible perpetraters?

@KonradRudolph: But could you please refute in details the claim in that blog post that serial passages through animals would not only generate the efficient ACE2 binding affinity but also create the appearance of natural evolution/selection?
https://harvardtothebighouse.com/2020/03/19/china-owns-nature-magazines-ass-debunking-the-proximal-origin-of-sars- through animal hosts would generate thecov-2-claiming-covid-19-wasnt-from-a-lab/

Hi @Hans I've addressed your point below.

Thank you, Michael. Let me ask you some followup questions. First, on serial passage. Why should serial passage always attenuate the virulence of a virus? Can we not artificially select for the virulence or affinity to say, the ACE2 receptor? Can you explicate the claim "Even if a virus had been serially passaged we could tell, because we know the mutational patterns of other serially passaged viruses." How is serially passed mutation different from that from "natural" evolution/selection?

This is a separatte question.

What do you mean by "a separate question"? Is it not a question on your claim? Or do you mean it is a deeper question requiring more in-depth knowledge to answer? Could you please explicate?

Dear @Hans, okay I agree "a separate question on this site" wouldn't be in a format that this site would permit, because you are starting from a base of limited understanding. So lets just its complicated but doable.

I understand. I will open a new question and invite you to answer. Would you be willing to share your knowledge there?

Hi @Hans it wouldn't work because you would need a base of technical understanding to open a question here. If you look at the number "closed" questions on the site, I think you get the idea. The site is traditionally for programmars/coders. This question is better for an alternative site such as biology stackexchange where it is more discursive. Please don't post it here.

By "base of technical understanding" do you refer to the technical understanding of programming/coding/statistics or (evolutionary) biology? I have quite a reasonable amount of the former but little of the latter. You said the site was for programmer/coder, would it not be fit for me to post the question here? Or perhaps you meant to say my question was more biological related per se than programming/coding related, thus more suitable for biology.stackexchange than here. Which is a more accurate interpretation of your statement?

Dear @Hans, firstly my apologise if you're a coder, so welcome :-). However,  I am recommending your question is posted on an alternative site, because I think the question here would be closed and that would be frustrating for you, and biology stackexchange is in my opinion the correct site. Apologies again for the misunderstanding.

No worries, dear Michael. Thank you for your welcome and suggestion. After browsing the questions on this site, I think you are right that this site deals more with statistics while my question is more concerned with the evolution theory. I have posted my question: https://biology.stackexchange.com/q/93152/8344. Would you care to give an answer there?

Michael, I just want to alert you to my question https://biology.stackexchange.com/q/93152/8344 in case you have not read my last comment. Would you be so kind as to provide an answer there?

Actually, I just recalled that I found a paper before demonstrating virulence gain of H5N1 influenza virus, specifically acquisition of airborne transmission capacity, after 10 series passage through ferrets https://science.sciencemag.org/content/336/6088/1534?fbclid=IwAR3gLlNFHuYbx_8mMgR93CPspMefJSvmOZa5knKKyJRWv4ImN0G_LuEyvjc . This seems to contradict your assertion regarding serial passage.

Here https://www.biorxiv.org/content/10.1101/2020.05.02.073411v1.full.pdf is another paper and result negating your assertion of virulence attenuation after serial passage.

serial passage could be for increase patogenicity... instead of get the virus strain that make live mouse (as selection natural), get the virus strain that kill mouse first or kill more youngs etc this each passage, so you have to cancell last paragraph

they found that sars2 is adapted as well as humans to a particular transgenic mouse with human ace2 gene and human immune system model... end discussion

Very interesting @RosLuP thanks

@Michael  they found that
<<A massive heptapeptide sharing exists between SARS-CoV-2 spike glycoprotein and human proteins. Such a peptide commonality is unexpected and highly improbable from a mathematical point of view, given that, as detailed under the “Methods”>>
in here https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7499017/  and see the graph
is very strange too the same number of peptide proteins in common with sars2 is found in "Mouse (Mus Musculus)"
too... someone has said

Is it true that viruses adapt to the species they invade using also parts of protein that the host uses (e.g. to camouflage)?

if that is true, if exist Mouse "Mus Musculus" GMO with humans genes that are find in sars2... that mouse has the more big set of peptides in common with sars2 more than humans and all mouses...


text:
Feature extraction methods that can handle inconsistent numbers of atoms for molecular dynamics
I want to compare the protein dynamics
) pH 7 versus pH 3, or
) wild type versus mutant
The protein will have slightly different number of atoms at each condition, due to protonation or mutation, which is bad for machine learning. The alternative is to only consider the alpha-carbon atoms, but this will lose significant amount of information (e.g. side chain, salt bridge).
Is there a way in molecular dynamics to extract features from ensembles of different number of atoms?
Note: I have used common methods to decipher the difference between two conditions, including RMSD, RMSF, Rg, secondary structure, salt bridge, SASA, etc. They can explain the difference to some extent. However, I feel the selection of properties is "subjective", and there must be some hidden features not known. So I am seeking a machine learning method, that can "automatically" capture the most essential properties based on the coordinates as inputs.

This answer is no longer on topic. The question is about what is the best machine learning algorithm to use to analyse atomic coordinates from an MD trajectory in order to infer a novel properties.

Also, it should be preface for any confused reader that unfortunately "feature" is a technical term in both genetics and statistics. For the former a feature is a region of DNA, RNA or protein that has a given property, such as open reading frame or PFAM domain. While in the latter a feature is a preferably independent property that forms part of a feature vector for different types of machine learning.

Heavy atoms are common
Assuming the pH difference only results in protonations, the protons differ between the two conditions, but not the heavy atoms. These will retain all the side chain information.
Merit of backbone only
(Section no longer relevant due to OP edit)
The common metric used to compare two structures is root mean square deviation (RMSD) of the alpha carbons to assess how much does the backbone differ and with analysis of MD simulations it is used to assess whether one is more dynamic than the other (e.g. this Bioinformatics SE question). The reason why the RMSD is done with the alpha carbons and not the minimum common subset of atoms is that backbone changes are large scale, while the side chain changes only affect locally —so their inclusion would needlessly dilute the score. The C prime and backbone nitrogens behave like the C alpha so there is no need for these.
Extraction of coordinates
Assuming by features you mean atomic coordinates, there are loads of ways. Using Pymol module for Python, one could do something like:
import pymol2
with pymol2.PyMOL() as pymol: #can be parallelised
    pymol.cmd.load('trajectory.gro', 'mytraj')
    pymol.cmd.remove('h.') #remove hydrogens
    pymol.cmd.intra_fit("name CA", 1) #align all steps in the trajectory
    for i in range(1, pymol.cmd.count_states('mytraj')):
        for atom in pymol.cmd.get_model(f'state {i}').atom:
            atom.coord #List[int int int]

One could change pymol.cmd.get_model(f'state {i}') to select state {i} and name CA for Cαs only.
Edit: Mutated residues
In the case of mutated residues or residue differences due to the two proteins being close homologues, "mutating" the differing residues to alanine, i.e. by stripping all heavy atom beyond the Cβ (keeping C′, Cα, H, N, O & Cβ) in the MD trajectory —or glycine if one is glycine. If they are homologues with different lengths, this approach would mean skipping residues that match up with gaps in the MSA. Even though not uncommon for statistically analyses, it does pose some issues statistically. It may unbury hidden residues, but length variant generally happens with loops.
In the snippet above, the pymol mutagenesis wizard could be used, but that is overkill, does not work in parallel and the approach does not port to other PDB reading modules in whatever language one wants. So say you have a list to_be_trimmed of integers of PDB residue indices (not pose residue indices):
if atom.resi not in to_be_trimmed or atom.name.strip() in ('C', 'CA', 'H', 'N', 'O', 'CB'):
      do whatever with atom.coord

In acidic environments, you get some funky things happening like exposed aspartate-glycine residue pairs forming isoaspartate or exposed aspartate-proline resulting in a cleavage, etc. For these cases, which are nearly never modelled, ignoring these values is the best option.
Edit: small comment on SASA
The investigated property list is very extensive and some reader may adopt them. So I thought I should mention that per residue SASA is generally converted to relative surface area, which is residue independent, by dividing it by the maximum SASA for that residue (values from Tien et al. 2013: 'A': 121.0, 'R': 265.0, 'N': 187.0, 'D': 187.0, 'C': 148.0, 'E': 214.0, 'Q': 214.0, 'G': 97.0, 'H': 216.0, 'I': 195.0, 'L': 191.0, 'K': 230.0, 'M': 203.0, 'F': 228.0, 'P': 154.0, 'S': 143.0, 'T': 163.0, 'W': 264.0, 'Y': 255.0, 'V': 165.0}), where RSA > 0.2 is an exposed residue.

thanks, I have updated my question.

source: Stack Exchange

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

sort: votes


text:
GWAS, MWAS, EWAS: what are the (in)dependent variables?
I started reading some papers on X-wide association studies, where X can be metabolome, epigenome, etc... The authors usually describe which are the dependent and which are the independent variables of their models. Taking as an example a metabolome-wide association study (MWAS), the authors use as model the following: metabolite ~ condition. That is, each metabolite is considered the dependent variable and the condition (e.g., smoker vs. non-smoker) is considered the independent variable. I was a little bit confused, so I looked into an example study using GWAS: in this case the model can be written as condition ~ SNPs, which is what I would expect (knowing nothing about GWAS but something about Data Science) and quite the opposite from the MWAS described above.
Does the choice of the dependent variable depend on the questions asked in the study, or does it exist some sort of agreement? Moreover, regarding the GWAS model described above, I read that usually researchers fit a model for each SNP: what's wrong with fitting a model which consists of several independent variables (i.e., y ~ SNP_1 + SNP_2 + ... + SNP_N)?

Fitting a model that consists of many independent variables is done, but serves a different purpose. There is nothing inherently wrong with it. A GWAS answers the question "given phenotype and genotype data, which genotypes are associated with the phenotype?" Your proposal answers a question more like "given phenotype and genotype data, how can we best predict phenotype with genotypes?"

See also: https://biology.stackexchange.com/questions/111538/why-would-greml-in-gcta-treat-snp-effects-as-random-effects/111539#111539

To answer the first part of your question, the dependent and independent variable of X-WAS is kind of arbitrary and dependent on the question you asked. But it gradually becomes a convention in the field after the initial name and concept are accepted by the community. For example, GWAS from the very beginning is written as condition ~ SNP, and there is no doubt about it. Later, people did something like condition ~ gene expression and called it TWAS, which is also widely accepted now. In your case, MWAS is indeed metabolite ~ condition, and the only important thing is to know the convention in the research community so people are on the same page when they talk to each other.
As for why running GWAS using single SNP. In a GWAS, it is typical to use 500K to 7M SNPs. So this is a large number. It is obvious that we cannot include all the SNP in one single model, then the question becomes how many SNPs to include in each model, and how to choose those SNPs. This is very hard to decide. Also, due to link disequilibrium, many SNPs are correlated, which means co-linearity, a thing you want to avoid in linear models. It is not clear how to correctly interpret the effect sizes in a model with many SNPs as well. Therefore, due to feasibility, simplicity, and interpretability, it might be best to perform regression on every single SNP.
Now since we are running GWAS in an SNP-by-SNP fashion, there are a lot of work that needs to be done after GWAS in order to find the true causal SNP. Many fine-mapping methods have been developed to tackle this problem.
Interestingly, there are some multi-variate GWASs methods, but they do not come in the flavor of simple linear models. Instead, sometimes complex Bayesian statistics are used, for instance.

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


text:
Software for microbial profiling from 16S rRNA gene sequence
I have hundreds of GBs metagenomic 16S rRNA gene sequence data. I want to do microbiome composition profiling (with relative abundance) from the data. Also after that, I will do functional profiling (gene and metabolome level).
Can you please suggest me the best software(s) for these two tasks?
Thank You

I tried dada2 and is not bad (if you know R).
QIIME2 is also an option.
Many other are available, the choice might also depend on your sample and your exact question.
For functional profiling you may try picrust2

For comprehensive microbial profiling of your 16S rRNA gene sequence data from Illumina MiSeq, I recommend using the web app Microbioma16S (www.microbioma16s.it). This powerful tool offers an end-to-end pipeline that covers everything from preprocessing raw data to conducting alpha and beta diversity analysis. With Microbioma16S, you can explore interactive results of taxonomy associated with features, taxa relative abundances, and perform comparisons between sample groups. It's a user-friendly solution for your microbial profiling needs.

The best thing to do is a phylogenetic tree in my opinion.
The trendy approach and the one deemed acceptable in publication is deep learning for precisely this application, such as the pipeline published here. This approach is attractive in 16S NGS metagenomics classification because it uses advanced approaches in deep learning, notably convolution neural networks (CNN) and deep belief networks (DBN). I don't really understand DBN, but I suspect it is a variant of recurrent neural networks (RNN) such as long-short term memory (LSTM).
Having said all that the accuracy is low, viz.

For instance, at the genus level, both CNN and DBN reached 91.3% of
  accuracy with AMP short-reads, whereas RDP classifier obtained 83.8%
  with the same data.

91% accuarcy for a deep learning neural network is minimal, but it will be quick and will do precisely what you want with little fuss. A phylogenetic tree on the otherhand can take several hundred hours of total processor time. The thing about a 16S tree is that it is 100% accurate, because you can manually identify "grey zones" to assign their classification.
You can retrain CNN models and in future this will likely resolve the accuracy issue. However, I don't know about about DBN and I don't think RSS can retrained, more specifically I don't know how to do it.

however, "16S phylogenetics" and "100% accurate" shouldn't go in the same sentence.

Sorry I mostly disagree @user3479780, but I agree it is better to say "the most accurate". However 1 the sentence is misrepresented, "... 16S tree is that it is 100% accurate, because you can manually identify 'grey zones' to assign their classification". 2 My work (500 cites) and other publications identified it as an excellent taxonomic locus at familial level. 3 Molecular phylogeny is the gold-standard in bacterial taxonomy. 4 Whole genome phylogenetics theory is complicated. 5 At the time many bacterial taxonomic deep-learning studies published had low accuracy (<90%).


text:
Analyzing proteins based on sequence similarity
Using NCBI Gene database, I retrieved list of 7015 protein which are implicated in autism pathogenesis. This list includes tax-ID, organism name, symbol, aliases, map location, chromosome, genomic nucleotide accession version, start position on genomic accession, end position on genomic accession, orientation, exon count and OMIM). Among these 7015, 1144 proteins are involved in neurodevelopmental functions. I want to know that how many of these 1144 proteins are functionally similar to 7015 proteins? Is there a method to query 1144 proteins against 7015 proteins and group them based on structural similarity/ function similarity/ gene family?. Since I have large set of proteins, I am trying to understand if sub-set of proteins share similar function with this large set of proteins.
Thanks
Priya

Please [edit] your question and tell us what kind of data you have. Are these known, annotated proteins? From what species? Do you just have a list of protein fasta sequences? Gene sequences? RNA? Do you instead have a list of protein accessions?

@terdon I have edited the question

Thank you! So these are human sequences? You mention OMIM, but it is better to be sure.

@terdon Yes I have list of protein names (not sequences) from Homo sapiens. Thanks

I suggest you give DAVID a try. Specifically, their Functional Annotation tool. Just enter your list of protein IDs, and it will return groups of proteins where particular GO functions are overrepresented.

This should at least help you start categorizing your list into functional groups.

Thanks. Is there a way to analyze whether if 2 proteins belong to the same protein family or not?

@PriyadarshiniThirunavukkarasu just look at the protein's annotations.


text:
Cuffmerge: EOF marker is absent. The input is probably truncated
I need to merge my all transcripts.gtf from cufflinks output follow this command line :
cuffmerge -o merged_gtf_output -p 15 -s ref.fasta -g anot.gtf assembly.txt

But I have this warning :
[bam_header_read] EOF marker is absent. The input is probably truncated.
[bam_header_read] invalid BAM binary header (this is not a BAM file).
File merge_B3//tmp/mergeSam_fileOaNffV doesn't appear to be a valid BAM file, trying SAM...
[21:40:04] Loading reference annotation.
[21:40:05] Inspecting reads and determining fragment length distribution.
Processed 3921 loci.
> Map Properties:
>       Normalized Map Mass: 33242.00
>       Raw Map Mass: 33242.00
>       Fragment Length Distribution: Truncated Gaussian (default)
>                     Default Mean: 200
>                  Default Std Dev: 80

Can someone tell me what this warning means and how to fix it?

it looks very weird to me, since you are not reading a bam file, some how samtools is called and it throws this warning, see also https://www.biostars.org/p/189139/

Yes it's true that it's weird, but I think it alludes to the files from which the gtf were generated by cufflinks

There is a post about this on seqanswers. The poster says the following about cuffmerge behavior:

I kept the temporary files while running cuffmerge. It builgs a SAM file from transcripts.gtf files. Then it seems first it reads it as BAM since I read samtools take the files as BAM default. Then it finds it is not a BAM so it tries SAM. Then it can read it and resume the task.

So it sounds like the BAM file is generated from GTF as a kludge inside cuffmerge, not the other way around.
Did you change any of the intermediate files or delete anything? Is there maybe a different cuffmerge run that it might be colliding with?
Did you get any output? it sounds like this error is frequently observed when actually the command completes fine, e.g. https://www.biostars.org/p/189139/. See other posts by googling the problem maybe, e.g. https://www.biostars.org/p/344439/
(I know this isn't a proper answer but it was too involved for a comment.)

Thank you for your answer, No I did not delete anything and despite the warning the command was executed well.


text:
Portals besides the National Cancer Institute GDC Data Portal for downloading histology slides
I am downloading sample slides (in svs format) from the GDC Portal for a convolutional neural network model, but I also need normal tissue slides to perform classification.
There are some normal tissue slides that are adjacent to the tumors, but I'm looking for other easy to use portals on the internet where I can download different types of normal tissue (bone, connective, etc).
Can anyone recommend similar type portals where I can download svs files (or equivalent) of normal tissue samples?


text:
Phasing partially phased genomic data
I am trying to create an integrated callset with coding regions from phased exomes and noncoding regions from a genotyping array.  I would like to merge the data together in the correct phase.  Since the exomes are already phased with read informed phasing, I don't want to rephase them because some of the phasing changes after doing this.  Is there a way to statistically phase genetic data that has already been partially phased?
My data is in VCF format.
Thanks

Answer from @user438383, converted from comment:
I would recommend checking out shapeit4 as it is able to integrate statistical phasing and read-based phasing like you are after. The manual page for WhatsHap might also be useful for you:
https://whatshap.readthedocs.io/en/latest/index.html


text:
How can I download the SNP genotyping information from GTEX?
I am new to the field of Bioinformatics. GTEx database is great; I know how to check whether a SNP is eQtl with the online tools of GTEx.
But my question is how to download the right data to draw a boxplot on my own. like this paper:

Please help me . I have been stuck on this for many days.

I agree that GTEx is great, however for SNP variation however, I recommend using gnomAD due to sheer number —and the high possibility that a reviewer will ask the same.

In recent years violin plots have taken over bloxplots as the common way to show such data —assuming few dataseries. Yes, five years ago showing a violin plot always required explaining but now the broader scientific community understands them.

@MatteoFerla THX for your note. I will go to learn how to use gnomAD. Would you please give me some hints that how to download the genotype data of a SNP (eg, sample 1 is AA ; sample 2 is AG, etc) from GTEx . I really want to know .

@MatteoFerla I want to perform eQTL analysis from GTEx. Help me please and that will save my life!

Cross-posted: https://www.biostars.org/p/429533/

@ATpoint It‘s deleted. Thank you for your advice.

There are two types of data you can get from GTEx. The first is summary statistics, basically the slopes and p-values you see in the boxplot. You can also get expression levels for all genes measured in GTEx.
However, in order to plot boxplot by yourself, you also need genotype data. This is not publicly available because of privacy issues. In order to get it, you have to apply to it on dbGaP. I think you can find relevant information on the GTEx website.


text:
How did researchers derive the Ramachandran "validation" contours?
I'm a beginner to structural biology and for fun, calculated the tortional angles of some 100, 000 proteins. Here is my Ramachandran plot:
When I went to look for "canonical" Ramachandran plots, I discovered most researchers overlay a "contour" plot of allowed torsional angles (see below). My two questions are:

How are the distributions (ie. countours) of protein secondary structures calculated/derived? (ie. How do I "cluster" these into beta-sheets, alpha helices?)
How can I derive them myself?

You have already calculated you ψ and φ angles, hence your plot and you simply wish to do a contour map, correct? A contour map is a 2D histogram basically. There are many plotting programs. I use plotly, eg. plotly/contour.
Your scatter plot looks like it's in MS Excel blue, but with Matlab style ticks set outwards. For the former, I have zero idea, while for the latter there is a function called contour. So could you specify the program?

Alternative, a commonly done thing is setting the alpha of the markers to something like 0.2 so that the overlapping ones are shown better. And also worth mentioning, that hexbin plots are becoming increasing popular this year for 2D histograms (via ggplots or py matplotlib)

Hi @MatteoFerla, thanks for your response. To clarify, the plot was generated in matplotlib. I have a better appreciation for how the plot is generated, but my question is: how do I know which cluster belongs to alpha helix, beta pleated sheets, etc.

It is derived from the crystal structure and is a known variable. It can be spotted by eye by plotting the scatter plot (or hexbin plot...) with the data filtered by secondary structure. There is no fancy k-mer clustering of similar needed. Alternatively it can be argued a priori: residues on a β sheet go up-down so the ψ must be about 180°, while the turn on a helix is 3 and a bit, so must be about 60 °...

To expand on Matteo Ferla's comment as an answer, there are two approaches to getting these "allowed" contours:

From probabilities found in observed crystal structures. You could imagine plotting a Ramachandran plot for the whole PDB and colouring the points by assigned secondary structure from DSSP. The contours would then be regions that encompass e.g. 99% of alpha-helical residues, or 99% of all residues. In doing this you would notice for example that glycine can adopt more regions of the Ramachandran space that other amino acids. PROCHECK is software that scores a protein based partly on the likelihood of the phi/psi angles of each residue.
From reasoning about the chemistry of the secondary structure of proteins. Only certain phi/psi combinations allow the i -> i + 4 pattern in an alpha-helix or the alternating up/down sidechain pattern in a beta-strand.

Thanks @jgreener! The plot above is a portion of PDB, but captures something on the order of 10^5 proteins. I can already make out what appears to match the secondary structures. I'm curious how I can rigorously derive the probability distributions from crystal structures, as an exercise!


text:
Merging based on Chromosome coordinates
I have a simple task which seems very complicated now
Two data-frame I have One contains the Chromosome coordinates of distal enhancer elements and other one rlog values of accessibility which I want to map.
The distal coordinates data
dput(head(a))
structure(list(Chr = c("chr6", "chr11", "chr20", "chr20", "chr10",
"chrX"), Start = c(154986726, 70314780, 15895696, 40706992, 29132014,
129525939), End = c(154987062, 70315059, 15895987, 40707264,
29132266, 129526265)), class = c("tbl_df", "tbl", "data.frame"
), row.names = c(NA, -6L))

The Rlog_ATAC_location data
Chr    Start    End  CMP1  CMP2  CMP3  CMP4  CMP5  CMP6  GMP1  GMP2  GMP3  GMP4  GMP5  GMP6  HSC1  HSC2  HSC3  HSC4  HSC5  HSC6
  <chr>  <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
 chr1   16103  16348 2.70  2.99  2.87  1.88   2.06  1.94 2.63   2.37  2.61  1.81  1.31  1.21 1.89  1.67  1.95  1.09  2.17  0.861
 chr1   96493  96796 0.976 0.325 0.943 0.216  1.33  1.77 1.66   1.33  1.47  2.39  1.76  2.37 0.939 0.198 1.37  0.439 0.868 0.756
 chr1  271114 271432 1.95  2.62  2.04  2.29   2.49  2.86 2.12   2.89  2.30  2.25  3.16  3.50 0.695 1.03  0.491 0.785 1.15  1.00
 chr1  273026 273333 0.740 1.88  1.48  1.71   2.50  2.88 2.81   2.66  2.89  2.07  3.02  3.61 1.24  1.40  0.984 1.35  0.651 1.53
 chr1  274265 274599 0.954 1.59  2.06  0.227  2.05  1.19 0.906  1.60  1.62  1.99  2.05  2.47 0.920 0.211 1.06  0.426 1.99  0.529
 chr1  402478 402771 2.98  3.26  3.01  3.07   2.51  2.76 3.29   2.65  2.34  3.03  3.00  2.99 2.54  3.04  2.88  2.96  2.94  2.68

I tried inner join ,merge etc but i can map the chromosome cooridnates but the expression values comes empty
tmp <- merge(a,Rlog_ATAC_location, by=c("Chr","Start" ,"End"), all.x=TRUE)

Chr  Start    End CMP1 CMP2 CMP3 CMP4 CMP5 CMP6 GMP1 GMP2 GMP3 GMP4 GMP5 GMP6 HSC1 HSC2 HSC3 HSC4 HSC5 HSC6 Mono1 Mono2 Mono3
 chr1  16104  16348   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
 chr1  96494  96796   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
 chr1 271115 271432   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
 chr1 273027 273333   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
 chr1 274266 274599   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
 chr1 402479 402771   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA

Tried left join as well i get same result.
I want to map the coordinates as well as the samples columns which contains the expression
Any suggestion or help would be really appreciated

Do your two data frames have the same coordinations? For instance, does chr1   16103  16348 exist in another data frame as well? If not, then you cannot use left_join or similar functions. If this is the case, maybe a more precise description of your problem is to annotate distal elements with your ATAC location?

I have updated my question with the data actul data Im using if i do left join I can map all the coordinates but the values are not coming "
Do your two data frames have the same coordinations?" yes one is subset of others

in the Example you've given us, there is no overlap in Start positions between the two data frames, so the behaviour of something like join is to return NA. You should go back and check with any(a$Start  %in%  Rlog$Start) if there is any overlap between the to data frames

i think i have messed up somewhere i checked its not matching

In your ATAC data frame, you have chr1   16103  16348. But in your merged data frame, you have chr1   16104  16348. Is it possible that you a is 1-based whereas your ATAC data frame is 0-based? This can happen in bioinformatics and can be annoying. This might be why you have a lot of NA in your merged data frame. You should subtract 1 from the start position of your a data frame while keeping the end position unchanged. Then try to merge again.
You should always be aware of your data format to know whether it is 1-based or 0-based.

i will definitely  do that


text:
Understanding some of the computational bottlenecks of Covid-19 research
I am a researcher in high-performance computing (with very little bioinformatics background), and I am trying to understand what are the current biggest computational bottlenecks of software used for research on the topic of Covid-19 (testing, vaccine development, drug identification, etc...).
I know this is a very broad question, but I am looking for the biggest common denominators. Said differently, if you could pick a few operations / algorithms / data structures to speed up Covid-19 research and they would magically run 1000 times faster, what would they be?
Note: links to books, papers, and scientific literature are very welcomed.

I can only speak of drug design (and even then I am terrible at turning down the jargon).
In the case of drug design, this is pretty much plan C. Namely, none of compounds that entered clinical trial at the start of the year work (let's call this plan A) and none of the vaccines that are entering now clinical trial work (let's call this plan B although vaccination prevents so is a better solution). Clinical trials take a lot time as they need to be thorough (cf. the infamous thalidomide from the days with more lax testing)
There are several bottlenecks in drug design. These are technical based on our imperfect understanding of things, which results in less than ideal models. CPU time is not really an issue. Empirical validation is required at many steps along the way. This is greatly increased by automation. The XChem team at Diamond light source have done a fragment-screening against the SARS-CoV-2 main protease at an insane pace: things are changing fast!
Drug screen
There are many ways to make a drug. In one way, a test system was set up (a reporter assay) that gives you a measurable read-out, say mm of clearance of bacteria on a Petri dish (viability) around a cellulose dish soaked in different separated natural extracts —and you might find that Penicillium chrysogenum has a compound of interest. What you are doing is drug screen. Screening a large panel of compounds (of varying purity) does not necessarily require knowing the protein target.
Rational drug design
Fast forward to this decade and it becomes possible to predict the free energy release upon binding of a small molecule against a protein —mostly.
This requires a solved 3D structure of the protein by X-ray crystallography. Crystallising a protein is a tricky business, as a pure protein sample has to come out of solution in a folded form to form a neat lattice as opposed to just making unfolded gunk. What salt conditions work? Anyone's guess. Some have done deep learning on it, but the best solution is to have a robot try a very large range of combinations.
Once you have a protein structure you can predict what fits into the part that you want to target, say an active site or an interface. This is not done by simple cartesian geometric shape fitting like LEGO bricks as there are many molecular forces involved. So a mathematical model of these forces is made, called a force field. Many have been made, a lot of study has gone into these, but they are not perfect. It is all made worse by the fact that we cannot really model water properly —most models use a proxy (implicit solvent) or use non-polarisable water.
This is all while made worse by the fact that a protein is not static but wobbles (following the same mathematical models). So a ligand does not slot in really like a key in a lock, but has an induced fit. Modelling all this is accurately is problematic, but one can still try: calculating where exactly a ligand binds to a protein is called for which one gets energy difference upon binding. And one can try an inordinate amount of small molecules —Zinc DB contains 230 million purchasable compounds. This is called a virtual screen. One can even start with the natural substrate and "weaponise it" by adding chemical groups that do a chemical reaction that breaks the protein. However, a metabolite is very different from a drug (cf. Lipinski's rule of five) as the latter must cross membranes etc. But generally the accuracy is not great without multiple steps of empirical validation.
To improve the results from a virtual screen the most common approach is to do a fragment screen to guide it. This is to automatedly soak protein crystals with a large library of small compounds and see which crystal contains a bound hit. The hits can be expanded into bigger fragments in a virtual screens and the short list of compounds tested (these are not cheap) by soaking and crystallising again.
In the case of Covid19 virus, the cool thing is that this whole process is open in the aforementioned current campaign (Diamond is a government owned synchrotron and not a big pharma). What is really nice about this is that for the second cycle comp-chemists from everywhere have done their own guided virtual screen and have submit compounds (which you can see publicly).
It should be said that drug screens that are not rational design are still a thing. The recent antibiotic-discovering deep-learning paper in Cell]6 is such an example, although it should be said that it was a feat not only in AI but also in automated empirical testing.
Future
Automation via better and open software will definitely improve. With a small drive possibly from the dissemination of easy to use and to code off-the-shelf electronics, such as Rasperry Pis or Arduinos. But mainly with more complex machines such as the Opentron system (a python-controlled robot).
But indubitably where the greatest improvement can come is better mathematical models. Possibly even with deep learning. After all researches were amused last year-ish (2018 competition) when DeepMind "won" the protein 3D structure (fold) modelling competition (CASP13) and there have definitely been some really nice results in the field of docking with hybrid models. So deep-learning may still surprise us yet...
Lastly, open science will be a strong driver: using premade packages or tools on GitHub and code snippets from StackOverflow etc. will accelerate these.

Drug screening
I understand the idea behind this, particularly with deep learning playing such an important role in computation (see GPU). Moreover, molecular dynamics is famous for sumping computer power, so yeah this is a bottleneck. However, this is not trivial because the wet-lab work is such a long way off. The first crystal/cyro/EM/NMR structure that will be resolved is Spike (S) and this would not be a 'traditional' drug target, its a vaccine candidate (okay influenza has surface antigen drug target, Tamiflu). Expressing S is a priority (vaccines) and solving the structure is a by-product and its not that big, plus the majority of it is water-soluable so it will cystalise nicely (its a surface antigen). Obtaining the structure for a more 'traditional drug target' which would be ORFa and ORFb (RNA dependent RNA polymerase ... and other replication associated proteins), will be difficult, because they are huge (ORFa + ORFb are 20 Kb of RNA, which for an RNA virus is twice as big as the genomes for most RNA viruses) and for crystals ... solubility(?). Ultimately, how do you cystalise a protein that 4kb of amino acids - well you'd have to do it in subdomains, you'd probably need NMR/ really hot EM.
General
HPC is a big field and it highlights the importance of one area most struggle with which is GPU (graphic processing unit) for deep learning (DL) applications. HPC's role molecular dynamics/drug discovery) and a mixture of DL and MD.
Traditional CPU 'farms' most of use have access to, is less critical because whether it takes 2 days per calculation or 2 hours isn't a huge issue (but its very seriously annoying indeed)..
The most pressing need for computational power will be epidemiological modelling, particularly if someone is pursing a microsimulation. The firepower would enable them to explore a large number of scenarios and parameter explorations. Apart from a vaccines then epi modelling is the hot-topic of the moment.


text:
Is there a aminoacid mutation that forms similar protein?
Is there an amino-acid shift that will 99% of the time end up resulting the same biochemical function/structure and implications? Like Alanine replaced by Leucine will code in 99% of the cases(at any position that first has Alanine) towards proteins that will function the same without big percentual changes in how good it binds or catalyses.
I do have a big amount of mutations. I'd like to filter them out on chance of biological/phenotypical impact, and since some amino acids look like good substitutes by eye (if I look at the molecular structure)

I am always surprised that people assume a one-liner question without any details on the context will lead to any sophisticated answers. Please invest some effort and add details. Since this is not really bioinformatics I suggest you post this in a more suitable community such as SE Biology. Be sure to provide more details there.

Maybe. For a more detailed answer, please add more details to your question.

What do you mean by "shift", if you mean frameshift the answer is no

A mutation is most often likely to be either neutral (no effect) or destabiling (protein misfolds). Changes in activity require changes at the active site or the entranceway to it —which are often only a handful of residues, so are very uncommon.
So the answer to your question is most mutations keep the activity, although a large fraction break it. While to get a new activity the process is more fluid and is not boolean.
NB. Not all proteins are enzymes, but for simplicity I will concentrate on enzymes.
Neutrality and core vs. surface
The outer residues of a protein, the surface residues, are very tolerant to mutations. Unless they are part of an binding surface (interface) or the change results in a charge change that is. The core residues on the other hand are more sensitive to mutations.
Destabilisation
One can do a scan in silico or in vitro or in vivo to determine the mutational landscape of a protein, which can take various forms.
This being Bioinformatics stack overflow I'll give an in silico example. The following for example is a heatmap I did some time back with Rosetta pmut_scan, which is a classic tool for this and shows in the y axis the residue it was mutated to and on the x axis the position. A larger change in Gibbs free energy of folding) relative to wild type (∆∆G) means it is destabiling, but anything less than 5 kcal/mol is neutral. As you can see many are neutral, but some half are destabiling —with a caveat for proline, which require a backbone change.

In the in vitro setting one speaks of melting temperature (DSC assays) or turnover number over Michaelis constant (enzyme specific assay). In the in vivo setting one speaks of a fitness cost of a mutation, i.e. the protein is less active so the organism grows less.
Drug resistance
Drugs inhibit enzymes by binding the active site. Therefore resistance is often costly. There are some nice papers that look into the fitness cost and mutational landscape of HIV protease or β-lactamase. In this situation a strong epistasis is seen. That is, a first mutation that confers resistance to the drug greatly destabilises the protein, while a second compensates for the first mutation.
Change in function
An enzyme can catalyse a reaction well on a given substrate, because it was evolved this way (technically the encoding gene evolved, but shh). This is the physiological activity. However, enzymes can also catalyse side activities on similar substrates with a lower turnover and higher Michaelis constant. This is especially true for substrates it was not selected against. This is called "enzyme promiscuity". The tolerance of the physiological activity to mutation is called "robustness", while the ability to increase promiscuous activities is "plasticity". Different enzymes have a different balance of the two.
So a mutation may change the strength of the physiological activity, but may result in a new promiscuous activity that can be used as a stepping stone for selection to increase this activity. Therefore, as you can see it is not a binary switch, except for a few very specific case.
Best amino acid for changes?
Proline is a very strange amino acid as it is a secondary amino acid. It lacks a hydrogen on the backbone nitrogen so it cannot form the required hydrogen bond for an α-helix. And can adopt a cis-conformation with minimal penalty (whereas it is about 5 kcal/mol for other chiral amino acids). So changes to this amino acid are often very destabilising.
Protein cores tend to be hydrophobic so changes to the hydrophobic residues are often deleterious, especially when going from a small residue to a larger one.
The remaining residues are frequently found on the surface so may be neutral. However, some are found in the active site and may be able to shift the parameters of the ligand space for that enzyme. A change in charge may result in an increased specificity for a different set of substrates for example. A size change from a larger amino acid to a smaller one may result in a broader range of promiscuous activities —mutations to alanine for example.
So ironically, polar and charged residues as they are frequently found on surfaces are the residues which are most often neutral, but in a small amount of cases they may change activity.


text:
truncate fastq.gz file without error in Snakemake
All,
As the first step in my Snakemake pipeline, I generate small fastq files from all input fastq.gz files. This is to run a quick round of "infer_experiment.py" (from RSeQC, to get the kit strandedness). However, Snakemake says this step fails. I think I know why, it is because a command like this
zcat fastq_raw/20190011-043_S43_R1_001.fastq.gz | head -n 400000 > infer_strandedness/20190011-043_S43_subset_R1.fastq

Ends with an error: gzip: stdout: Broken pipe and Snakemake uses Bash strict mode. (Snakemake reports: (one of the commands exited with non-zero exit code; note that snakemake uses bash strict mode!)). So my question is:
How to avoid this, or, are there other truncating methods that do finish cleanly? Preferebly they are also efficient (seqtk offers subsampling but this parses the entire fastq file!).

You're better off using seqtk even if it's slower since you're then not performing a highly biased sampling from the beginning of the file, where the quality tends to be worse.

But the library is randomized during prep and flowing into the chip, why would the first part of the fastq file (first scanned tile of the chip) be biased? A bit worse quality would not be a problem if I only want to detect the strandedness...

The beginning of the files represent the edge of the flow cell where the quality is lower.

Snakemake is using bash strict mode to ensure best practice error reporting in shell commands. Unfortunately some tools violate bash strict mode.
More details of the error, here. Alternatively, you can use if you don't want disable the strict mode.
zcat file.fastq.gz  | awk '(NR<=400000)' > subset.fastq

This solves my problem exactly (the awk solution), it is fast and does what my original command does but without pipefail error. Thank you.

@Freek another option would be to add || true to give a 0 exit status


text:
nextflow `fromPath` not passing all files with *
main.nf is
#!/usr/bin/env nextflow

params.input = ''
tchan = Channel.fromFile(params.input)

tchan.println()

I have a.bam, b.bam and c.bam
I run:
nextflow run main.nf --input *bam

But it only prints a.bam
According to the documentation, wildcard character it is interpreted as a glob path matcher
What am I doing wrong?
Thank you

Shell expands wildcards, try ls * vs  ls "*" .

got it!
I should run:
nextflow run main.nf --input "*bam"

instead of
nextflow run main.nf --input *bam


text:
Running a script in python
I am trying to run a script in python but I getting this error
loan-mac-13:vcf2maf fi1d18$ python
Python 2.7.10 (default, Feb 22 2019, 21:55:15)
[GCC 4.2.1 Compatible Apple LLVM 10.0.1 (clang-1001.0.37.14)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import os
>>> os.chdir("/Users/fi1d18/Downloads/vcf2maf-master/vcf2maf")
>>>
>>>
>>> os.getcwd()
'/Users/fi1d18/Downloads/vcf2maf-master/vcf2maf'
>>> vcf2maf.py -h
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
NameError: name 'vcf2maf' is not defined
>>> python vcf2maf.py -h
  File "<stdin>", line 1
    python vcf2maf.py -h
                 ^
SyntaxError: invalid syntax
>>>

Where I am doing wrong?
Thanks for any help
By your comments I run
Sorry says that Last login: Wed Apr  1 11:50:14 on ttys000
pc-133-235:~ fi1d18$ pwd
/Users/fi1d18
pc-133-235:~ fi1d18$
pc-133-235:~ fi1d18$ cd Downloads
pc-133-235:Downloads fi1d18$ python vcf2maf.py -h
Traceback (most recent call last):
  File "vcf2maf.py", line 19, in <module>
    import vcf
ImportError: No module named vcf
pc-133-235:Downloads fi1d18$

As I said, that's a completely separate issue. Please ask a new question about it. First, however, just copy/paste the error into your favorite search engine. That will give you the answer.

I am sorry to say, but this is a bioinformatics equivalent of "What is a cell" in biology. It appears to me that you really should go through some python tutorials to get over the basics before you proceed with your research. I heavily recommend https://software-carpentry.org/ as a starting point.

You %100 right @KamilSJaron, I am nothing in the first line of python

You are trying to run a script file in the python console. That doesn't make sense. The console is for running commands manually, interactively. A script is a collection of commands. You can't just run that by name in the console. If you want to use it, copy the contents of your script and try to run them in the console.
The way you run a python script file is simply:
python vcf2maf.py -h

Or, assuming you have made the file executable and it has the right shebang line, just:
./vcf2maf.py -h

The specific error you are getting is because foo.bar in Python means "run the bar method from the object foo. So vcf2maf.py means run py from vcf2maf. Since there is no object defined with the name vcf2maf, you get the name 'vcf2maf' is not defined.

Sorry says that Last login: Wed Apr  1 11:50:14 on ttys000
pc-133-235:~ fi1d18$ pwd
/Users/fi1d18
pc-133-235:~ fi1d18$
pc-133-235:~ fi1d18$ cd Downloads
pc-133-235:Downloads fi1d18$ python vcf2maf.py -h
Traceback (most recent call last):
File "vcf2maf.py", line 19, in
import vcf
ImportError: No module named vcf
pc-133-235:Downloads fi1d18$

@Angel I can't read that. Please [edit] your question to add extra information, don't put it in the comments. And tell me who/what "says" and when. I can't guess what commands you are running. If that's the output of running the script, then it is a completely different issue, you are missing the vcf module. So please ask a new question about that.


text:
Allele count is zero
I ran freebayes on bacterial data and got a vcf file. Some of the SNPs found have AC=0, but in the bam file I can see that there are alternative alleles, so in my understanding it supposed to be AC>0.
example:
in vcf file:
CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  55509   55508   55506   55507
scaffold1       78      .       A       C       9.27683e-13     .       AB=0;ABP=0;AC=0;AF=0;AN=3;AO=466;CIGAR=1X;DP=2254;DPB=2254;DPRA=0;EPP=20.9226;EPPR=29.3274;GTI=0;LEN=1;MEANALT=2;MQM=50.3412;MQMR=45.5379;NS=3;NUMALT=1;ODDS=2080.75;PAIRED=0.899142;PAIREDR=0.533932;PAO=0;PQA=0;PQR=0;PRO=0;QA=16347;QR=63107;RO=1783;RPL=203;RPP=19.7856;RPPR=291.853;RPR=263;RUN=1;SAF=465;SAP=1006.25;SAR=1;SRF=894;SRP=3.04075;SRR=889;TYPE=snp       GT:GQ:DP:AD:RO:QR:AO:QA:GL      .       0:126.704:776:626,149:626:22118:149:5255:0,-1217.7      0:126.704:834:663,170:663:23671:170:6004:0,-1224.79     0:126.704:644:494,147:494:17318:147:5088:0,-770.853

but with samtools tview:

Can someone tell me if this makes sense and why it happens?

The AC is the number of alternate alleles from called genotypes, you can see this in the vcf header:
##INFO=<ID=AC,Number=A,Type=Integer,Description="Total number of
alternate alleles in called genotypes

Assuming diploid, if we have 3 individuals, 0/0 , 0/1, 0/0, AC will be 1, because it detects 1 count of alternate allele.
If you look at your genotype calls, all are 0, so AC = 0

Oh, now I get it. freebayes chose 1 genotype per samples, due to majority of the reads. thank you for your answer!

you're welcome :) happy genotyping


text:
qvalue and p.adjust functions
Between the qvalue package and the p.adjust function, which is more appropriate to use when trying to calculate the q-values of a dataset? According to the manual for the q-value package, the q-value calculated is not an adjusted p-value which is what the p.adjust function would return. There seems to be an equal amount of people that use one or the other, but I wanted to see other opinions that could help clarify what would be most appropriate to use.

"then q-value is equal to FDR" Does it mean if I use the parameter"pi0=1" in qvalue(), the outcome will be equal to adjust.P(pvalues,"BH")?
Could you explain it by citing an example?
Any answer is appreciated!Thanks a lot!

p.adjust can calculate adjusted p-values with different methods. I think you can choose Bonferroni correction and BH correction, for instance. BH correction gives you False Discovery Rate (FDR), which is one of the most commonly used I guess.
Q-value is different and cannot be calculated by p.adjust. The idea is that when you have many p-values to adjust, they form a distribution. If all p-values come from the null hypothesis, the p-value distribution is going to be uniform on [0, 1]. But usually, our p-values are a mixture of the null distribution and alternative distribution. The parameter $\pi_0$ describes the proportion of p-values (or tests) that comes from the null hypothesis. To calculated q-values, the algorithm first estimates $\pi_0$, then calculates the q-value.
Note that if you set $\pi_0=1$, which means that you believe all the p-values come from the null distribution, then q-value is equal to FDR.
I think most people do not know about the meaning and details of q-values and it might be sufficient to use BH procedure in most cases.


text:
Calculating the number of probes for a given genomic range
I have a big data frame like this for copy number (exome seq)
    > head(df)
  Chromosome    Start      End
1 1             64613  5707515
2 1           5712940  5732322
3 1           5732399 16383682
4 1          16383742 16389288
5 1          16390813 16830026
6 1          16830201 17278112

I also have coordinate of human exams like this
Exone   Choromosome     Start_of_exone  End_of_exone
uc001aaa.3  chr1    11873   12227
uc001aaa.3  chr1    12612   12721
uc001aaa.3  chr1    13220   14409
uc010nxr.1  chr1    11873   12227
uc010nxr.1  chr1    12645   12697
uc010nxr.1  chr1    13220   14409
uc010nxq.1  chr1    11873   12227
uc010nxq.1  chr1    12594   12721
uc010nxq.1  chr1    13402   14409

Actually I need to interset my copy number coordinates and human exome coordinates to see there are how many exomes in these coordinates
I have used this code but all returns 0
snp_table <- expanse coordinate file
genomic_ranges <- df

how_many_in_range <- function(coords){
    coords = as.vector(coords)
    sum(snp_table$Chromosome == coords[1] & snp_table$Position > as.numeric(coords[2]) & snp_table$Position < as.numeric(coords[3]))
}

genomic_ranges$number_of_snps <- apply(genomic_ranges, 1, how_many_in_range)

Can you help me?

I am wondering how you got that calculation from google about the X2 thing. Can you describe in detail what is in your data frame?

I'm sure that's not what @PhoenixMu is asking you - you're literally just telling us what the header already tells us. Some real information would give us context on the data contained in the data frame. Are they probes? exons? random regions? And "google says" is not a valid citation. What is the source of your information on probes to exons ratio?

This is called by "merge by overlap" see for example this (my preferred) answer using data.table::foverlap()

Merge then aggregate.

This command gave me the number of axons
I hope this is correct
Please correct me if I am wrong
how_many_in_range <- function(coords){
    coords = as.vector(coords)
    sum(human_exome_coordinates$Choromosome == coords[1] & human_exome_coordinates$Position > as.numeric(coords[2]) & human_exome_coordinates$Position < as.numeric(coords[3]))
}

genomic_ranges$number_of_exon <- apply(segment_mean_table, 1, how_many_in_range)

Regardless of what your specific data is, the goal of your analysis is really to find overlaps between two data frames of genomic coordinations. My favorite package for doing this is GenomicRanges. It is very fast and can take care of different strands (if your data contains strand information).

Sorry I am too stupid in calculation, how I can do that?

Just a quick example, suppose I have two data frames. df1.gr=GRanges(Rle(df1$chr), IRanges(df1$start, df1$end)), df2.gr=GRanges(Rle(df2$chr), IRanges(df2$start, df2$end)). And then ov = findOverlaps(df1.gr, df2.gr). You can also customize some arguments in the findOverlaps function. You can also find more on the GenomicRanges webpage.


text:
co-localisation error code using co-loc
Can anyone offer advice on what the below error code relates to. I am using co-loc R package with 2 GWAS datasets. Using coloc.abf function with pvalues, MAF, sample size and case/control cohort. Thank you
  results <- coloc.abf(  dataset1=list(snp=as.character(1_subset[,2]),pvalues=1_subset[,11],N=1_data_subset$sample_size, type="cc",s=0.327,MAF=1_data_subset[,12]),   dataset2=list(snp=as.character(2_subset[,2]),pvalues=2_subset[,8],N=2_subset$sample_size,s=0.4319,MAF=1_subset[,12], type="cc"),p12=1e-06)

error code: Error in data.frame(V, z, r, lABF) :
  arguments imply differing number of rows: 1994, 1992
In addition: Warning messages:
1: In 2 * N * f :
  longer object length is not a multiple of shorter object length
2: In 2 * N * f * (1 - f) :
  longer object length is not a multiple of shorter object length
3: In r * z^2 :
  longer object length is not a multiple of shorter object length

Hi, I can't replicate your error message. Can you double check that the line of code you inputted is the same as that which gave the error message.

The only thing different is that I put the code in a loop to make code quicker if(length(tmp)>0){
results = coloc.abf(  dataset1=list(snp=as.character(1_subset[,2]),pvalues=1_subset[,11],N=1_subset$sample_size, type="cc",s=0.327,MAF=1_subset[,12]),
dataset2=list(snp=as.character(2_subset[,2]),pvalues=2_subset[,8],N=2_subset$sample_size,s=CCRATIO,MAF=1_subset[,12], type="cc"),p12=1e-06)
}

In your dataset2 = list(), the MAF argument should be MAF=2_subset[,12] instead of MAF=1_subset[,12]. Since your two datasets have different numbers of lines, this is what causes the error.
I am not sure whether this is a typo in your code, or it is because you do not have MAF information in the second dataset. What I like to do before COLOC is to filter two datasets so they have the same number of lines, and then feed those filtered datasets to the coloc function.
(Also just out of curiosity, is 1_subset a data frame for your GWAS results? I think it is not recommended to have variable names starting with a number.)

Thank you for your reply. That is really kind. I dont have MAF for dataset two so used dataset 1 for both. I will try the filter function and then maybe input MAF detail first onto the second dataset using 1st dataset. I did change the names, didnt want to use the disease names on forum. Thank you for double checking! Have a safe day, kind regards!


text:
umap highlighting two different models
I'm trying to create a umap for single cell data from human samples and ptx samples. I can get the umap to where it shows the umap with the different clusters but I want to show where the ptx samples are and where the human samples are.
My code is as follows:
#create the Seurat object

OD_10K_HUMAN <- CreateSeuratObject(counts = HUMAN_OD_10K.data, min.cells = 0, project = "human")
SD_5K_HUMAN <- CreateSeuratObject(counts = HUMAN_SD_5K.data, min.cells = 0, project = "human")
BNL.5K <- CreateSeuratObject(counts = SD_BNL_5K.data, min.cells = 0, project = "ptx")
BNL.10K <- CreateSeuratObject(counts = OD_BNL_10K.data, min.cells = 0, project = "ptx")
BNM.10K <- CreateSeuratObject(counts = OD_BNM_10K.data, min.cells = 0, project = "ptx")
BNM.5K <- CreateSeuratObject(counts = SD_BNM_5K.data, min.cells = 0, project = "ptx")

#merge data
scData <- merge(BNL.10K, y = c(BNL.5K, BNM.10K, BNM.5K, SD_5K_HUMAN, OD_10K_HUMAN), add.cell.ids = c("A", "B", "C", "D", "E", "F"), project = "HTB2876")

mark the mito genes
mito.genes <- grep(pattern = "^MT-", x = rownames(x = scData), value = TRUE)
length(mito.genes)
scData[["percent.mt"]] <- PercentageFeatureSet(scData, pattern = "^MT-")

scData[["log_nCount_RNA"]] <- log2(scData[["nCount_RNA"]]+1)

# remove cells with <200 RNA molecules, or >6000 molecules, or >30% mito
scData <- subset(scData, subset = nFeature_RNA > 200 & nFeature_RNA < 6000 & percent.mt < 30)
scData <- NormalizeData(scData, normalization.method = "LogNormalize", scale.factor = 10000)
scData <- FindVariableFeatures(scData, selection.method = "vst", nfeatures = 2000)

all.genes <- rownames(scData)
scData <- ScaleData(scData, features = VariableFeatures(object = scData), vars.to.regress = c("nCount_RNA"))
scData <- RunPCA(scData, features = VariableFeatures(object = scData))
#DimPlot(scData, reduction = "pca")

numPC = 20
scData <- FindNeighbors(scData, dims = 1:numPC)
scData <- FindClusters(scData, resolution = 0.4)
scData <- RunUMAP(scData, dims = 1:numPC)

DimPlot(scData, reduction = "umap", label = TRUE)

By default, DimPlot colors the cells by their "identity class", stored in Ident(scData) in your case. This indeed contains the clustering information.
However, you can you use the group.by argument for DimPlot to pass another column from scData@meta.data to color cells by. In your case you would need a metadata column that stores the ptx/human information.
The merge function automatically stores the original "project" information of the individual Seurat objects that you combine under orig.ident in the metadata. Since you set the project names to either "ptx" or "human", the orig.ident column should contain the information you need.
Something like this should work:
DimPlot(scData, reduction = "umap", group.by = "orig.ident", label = TRUE)


text:
Is it possible to determine the genomic context of a gene in a Whole Genome Shotgun project?
I performed several BLAST searches and identified several genes of interest in a Whole Genome Shotgun (wgs) project. I know that gene X is located on scaffold 1234, from nucleotide 1 - 2250, while gene Y is located on scaffold 2571, from nucleotide 3380 - 5790.
Given this information is it possible to determine if these genes are close to each other in the genome?
Many thanks for your help,

Leah

great question.. i guess if you cannot assemble the contigs, chances are you are unlikely to get it from the current dataset.. which organism is this?

Thanks @StupidWolf! It is a wgs from an amoebozoan protist.

@michael, I would very much appreciate your thoughts on this question!


text:
merging two/or more bed file into one bed file
I am trying to merge two bed files (more in future) to one. my bed files are something like :
.
I need to merge them in a way to have the shared chromosome location.
Is there a way to do that ?

In my opinion, the easiest way to merge bed files is to use bedtools merge.

Convert each of your Excel spreadsheets to tab-delimited text files (how-to).
Install and run Cygwin, if using Windows. If you are using OS X, open up the Terminal app.
Install the BEDOPS toolkit to get sorting and set operation tools (if you're not averse to using them).
Strip out the Microsoft line endings, strip the header, and sort each of the N files:
$ tr -d '\r' < fileN.tsv | tail -n+2 - | sort-bed - > fileN.bed

Merge the files:
$ bedops --merge file1.bed file2.bed ... fileN.bed > answer.bed

The file answer.bed will contain all merged regions from files file1.bed through fileN.bed. Merging will be done on intervals for each unique chromosome name.

Thanks that worked


text:
Is there any value in scaffolding the output contigs of MEGAHIT assembler given a metagenomic dataset?
As far as I understood, for most assembly programs, the scaffolding step takes into consideration paired-end information in order to get from contigs (contiguous sequences) to scaffolds (longer sequences that might have some N-filled gaps).
My assembly software of choice, MEGAHIT, uses paired-end information to build the contigs, but it does not output a standard scaffold. So I am wondering the following, related things:
-is it meaningful to run a scaffold program on the output of MEGAHIT? I imagine there might be some instances in which paired-end information could span a gap.
-Which software would you recommend for it? (I've tried soapdenovo2 and SSPACE, but they appear not to be actively maintained so I have the issue of 'it doesn't work and I can't do anything about it')
-Could relevant information (eg. two contigs being connected by pair-end information) be recovered by alternative and perhaps more user-friendly means, such as exploring the assembly graph with Bandage?
Thank you for your time!

Update 2: It looks like your approach has actually been suggested here as one way to use the PE information. I guess megahit may not be really using the PE information anyways. I still believe that it's kind of weird but other people do suggest it, so maybe it's worth trying what Torsten suggests.
I think that scaffolding with PE information from the reads used as input to Megahit is somewhat sketchy. You could still try it, but I'd be worried about artifacts, just because de novo assemblers are so heuristic.
However, I think that it is perfectly ok to scaffold using orthogonal data. Here are some examples of orthogonal data:

A reference assembly (weird for metagenomics but still possibly ok).
long reads, e.g. ONT or PacBio (low coverage probably fine).
Proximity ligation data, e.g. Hi-C reads.
Mate-pair reads (also lead to artifacts at a somewhat high rate, but technically valid).
Linked-reads (e.g. 10X).
Probably other methods as well.

The tools used in each case would be somewhat different. For a little more information about how you might use this, here is a recent review of tech.
Full disclosure: I work for a company that sells Hi-C kits for such applications.
Update: Realized that I missed one part of the Q. I think that visually exploring the assembly using e.g. bandage is always a good idea. Quite possibly you can make some scaffolding decisions that way, but to me it sounds somewhat painful to do, especially in a metagenome where there are going to be a lot of multiple-branching collapsed regions.

You can use SOAPdenovo-Fusion to scaffold contigs produced by MEGAHIT as suggested by one of the developers: https://github.com/aquaskyline/SOAPdenovo2


text:
Cannot blast against specific NCBI databases
I am having issues with some prokaryote reference genome databases (exact names : ref_prok_rep_genomes.*), that I downloaded from the NCBI website : https://ftp.ncbi.nlm.nih.gov/blast/db/. Files in this database go from ref_prok_rep_genomes.00 to ref_prok_rep_genomes.12.
In Linux command line, I was able to launch successfully the bash command  blastn -query (my_query) -db ref_prok_re_.genomes.* -out (my_output) with indexes from 00 to 05, but this does not work from indexes 06 to 12 (with the same query and output files). The exact same thing occurs with blastdbcmd -info -db ref_prok_rep_genomes.*.
Each time the error displayed is the same:
BLAST Database error: Input db vol does not match lmdb vol.
All the ref_prok_rep_genomes files are in the same folder. I downloaded them using wget and extracted them using tar -zxvf. All ref_prok_rep_genomes.* databases (indexes 00 to 12) have nhr, nin, nnd, nni, nog and nsq files. The folder also contains ref_prok_rep_genomes files (without numbers), with  nal, ndb, nos, not, ntf and nto extensions respectively, as well as taxdb.btd and taxdb.bti files.
Do you guys have an idea of what's going on?
Thank you very much for your time :)

Did you checked the integrity of the downloaded files? wget ftp://ftp.ncbi.nlm.nih.gov/blast/db/ref_prok_rep_genomes.??.tar.gz.md5 && md5sum -c *md5

I just did. With any of the md5 files, it gave me
ref_prok_rep_genomes.01.tar.gz: OK
ref_prok_rep_genomes.02.tar.gz: OK
ref_prok_rep_genomes.03.tar.gz: OK
ref_prok_rep_genomes.04.tar.gz: OK
ref_prok_rep_genomes.05.tar.gz: OK
ref_prok_rep_genomes.06.tar.gz: OK
ref_prok_rep_genomes.07.tar.gz: OK
ref_prok_rep_genomes.08.tar.gz: OK
ref_prok_rep_genomes.09.tar.gz: OK
ref_prok_rep_genomes.10.tar.gz: OK
ref_prok_rep_genomes.11.tar.gz: OK
ref_prok_rep_genomes.12.tar.gz: OK```

Actually I understood what happened. NCBI's new version 5 databases use virtual LMDB memory. They give issues with the previous versions of ncbi-blast+.
I managed to fix the problem by (1) updating ncbi-blast+ (currently it is the 2.10.0 version), (2) deleting my previous databases, (3) downloading them again using update_blastdb.pl --decompress <database_name> instead of wget


text:
Compressing the human genome to few megabytes
Multiple sources (see for instance this or this) discuss how genetic data will have scalability problems, given the huge file size of the human genome. The most straightforward encoding (see here) of the human genome requires about 700 Mb.
I came across this paper claiming to be able to store the human genome in about 4 Mb, having a reference genome and exploiting the fact that all human genomes are mostly equal. This paper is much older than the the other references discussing the scalability problems. Why is this technique not widely used?

https://bioinformatics.stackexchange.com/questions/12824/theoretical-limit-of-human-genome-compression?noredirect=1&lq=1

Answer from @alex-reynolds, converted from comments:
As with any compression scheme, if you have to store a reference genome, along with the 4 Mb "diff", is it really compressed?
http://mattmahoney.net/dc/dce.html
For sequencing data, you might also look at the CRAM format, which also uses diffs against a reference genome to store reduced datasets. Labs are switching from BAM to CRAM to manage storage of increased numbers of sequencing datasets, with associated changes to pipelines and toolkits to use this format. For performance reasons, BAM may be preferred for analysis, and CRAM for long-term storage.
https://gatk.broadinstitute.org/hc/en-us/articles/360035890791-SAM-or-BAM-or-CRAM-Mapped-sequence-data-formats
https://www.ga4gh.org/news/cram-compression-for-genomics/

I saw this paper a few weeks back that claims to offer 4x more compression than typical formats for both fastqs and BAM/CRAM files.
Genozip - A Universal Extensible Genomic Data Compressor
https://academic.oup.com/bioinformatics/advance-article/doi/10.1093/bioinformatics/btab102/6135077


text:
run a tool with python>=3.6 on galaxy-project instance that runs on python2.7
I am developing a wrapper for a bioinformatics tool writen in python3.
The native Galaxy-project runs on python2.7.
How can I the tool that requires python3 in a galaxy-project instance?

Each Galaxy tool uses its own (typically conda-based) environment, which will normally not use the same python as Galaxy itself (note that Galaxy can run under python 3). An example of this is deepTools, which is python 3 only but will happily run under a python 2 Galaxy since conda is being used to create a python 3 environment for it. It's best to distribute your tool via bioconda and then Galaxy will be easily able to install it.

I hope I did it the right way: 1- conda create ..., 2-conda install ..., etc...

Yes, though it's best to allow galaxy to make the env for you and put a copy of the wrapper in the toolshed.


text:
A box plot of qualitative variables
I have 2 groups of patients; responders to chemotherapy and non-responders
I have calculated cancer cell fraction (CCF) of a set of driver genes (by variant allele frequency) for each group and I have
> head(dat)
    Response       CCF
1 Responders 1.0000000
2 Responders 0.5413323
3 Responders 1.0000000
4 Responders 1.0000000
5 Responders 1.0000000
6 Responders 1.0000000

unique(dat$Response)

[1] Responders    Nonresponders
Levels: Nonresponders Responders
>

If CCF > 0.95 , the mutation is clonal otherwise sub clonal
I want to show how many clone and subclone are in these two groups by box plot or something similar like below

I have tried this which was nonsense
ggplot(dat, aes(x=Response, y=CCF)) +
    geom_boxplot()

Can you help me?

I'm not sure that boxplot will be the more appropriate representation, as you will end with two numbers (count of Clonal and SubClonal) per groups of patients.
One solution will be to first create a new categorical variable based on CCF values using for example ifelse statement for example by writing:
df$Clonal <- ifelse(df$CCF > 0.95,"Clone","SubClonal")

Then you can get the count of each Clone and subclone for each responders and non responders by using table:
DF <- as.data.frame(table(df[,c("Response","Clonal")]))

You can finally convert these new values as percentage by grouping according to the Response. For example, using dplyr, you can do something like that:
library(dplyr)

DF %>% group_by(Response) %>% mutate(Freq_percent = Freq / sum(Freq))

And finally, you can plot it as a barchart or single point in ggplot2.
Hope it helps you to figure it out how to deal with your data.


text:
Theoretical limit of human genome compression
How small can a compressed file containing the human genome be?
I'm aware that this question cannot actually be answered, since it is asking for the Kolmogorov complexity of the human genome, which is not computable. So reformulating:
What is the largest known lower bound that we can put on the size of a file containing the human genome?
An example answer could be: on average we know that the human genome contain at least this amount X of randomness, therefore it is not possible to do better than Y bytes.
I'm agnostic on which human genome, but I would say without polymorphisms to make it simpler.
The compressed sequence would need to be stand alone, without using a reference (I asked a question on that topic here: Compressing the human genome to few megabytes).
There should be a program able to compress and decompress the genome in a lossless way. Apart from that, any kind of modelling and coding is allowed.

For the most modern CHM13 assembly (v2.0; Complete T2T reconstruction of a human genome; https://github.com/marbl/CHM13), the best known compression without any use of a reference is 567 MB (1.45 bits/bp; 1.80 bits/distinct canonical 31-mer).
Upper estimates of genome compressibility can be made using standard general or specialized compressors (see, e.g., this comparison). In my experience, xz can achieve great compression ratios on genomic sequences, especially when used with the maximum compression level (-9) and a single thread (-T1). The compression performance can be further boosted by using a large dictionary (-e --lzma2=preset=9,dict=512MiB). Importantly, the original FASTA file needs to be cleaned – linebreaks hamper compression as they act as noise; the reformatting can be done by seqtk seq -U (which also converts all nucleotides to the upper case). It the case of the CHM13v2 assembly, xz compresses the human genome to 634 MB.
However, even better performance can be achieved with specialized compressors, especially if they have particular support for reverse complements and mutational models. Among the best ones is GeCo3, which has been shown to compress CHM13v2, with a particular combination of parameters, to 567 MB 3.
Lower estimates are impossible in general (we can't reject the existence of a program that would create the whole sequence from small amount of starting conditions), but can be done for particular compression techniques, e.g., particular variants of dictionary compressors. However, this leads to intractable optimization problems (see e.g. the problem of minimal string attractors).
A useful auxiliary measure for estimating the compressibility of single genomes is the number of distinct k-mers for a "reasonably" selected k (e.g., k=31). The lower bound on the compressibility of k-mer sets can then be roughly approximated by 2 bits per distinct k-mer, however, this is limited by the embedded structural assumptions, e.g., that single genomes are in principle random strings with occasional (possibly distant) repeats. For instance, the CHM13v2 human genome assembly contains 2,512,390,070 distinct canonical 31-mers, which would then correspond to 628 MB assuming 2bits/kmers, while the best available specialized compressors can compress it down to 567 MB, i.e., by approx 10% better.
Updates

May 27, 2022 – updated the upper bound based on the latest results with CHM13v2 (xz & GeCo3), changed the lower bound wording to an auxiliary estimate as the assumptions have shown not to hold

There's a benchmark that compared dozens of compressors on DNA data, including human genome: http://kirr.dyndns.org/sequence-compression-benchmark/
In particular, this table shows file size of compressed human genome, produced by the strongest setting of each compressor: http://kirr.dyndns.org/sequence-compression-benchmark/?d=Homo+sapiens+GCA_000001405.28+%283.31+GB%29&cs=1&cg=1&com=yes&src=all&nt=4&only-best=1&bn=1&bm=size&sm=same&tm0=name&tm1=size&button=Show+table
The best result is currently 657 MB on that particular genome build (GCA_000001405.28), when compressed with fastqz.
By the way, xz (mentioned in another answer) is currently number 20 in that table. (Though the genome is un-modified straight from the database, not the "cleaned" form.)
Disclosure: I worked on that benchmark.
Note: What is being compressed here is a FASTA file containing a particular assembly of a human genome. Compared to an actual genome it has this extra content: Sequence names and brakes between contigs, end of line characters, N characters representing unknown bases, and soft-masked repeats (represented as lowercase). On the other hand, it also misses some information from the actual genome: 1. It's an unphased haploid sequence, while the real genome is diploid (two copies of each chromosome). 2. It includes no epigenetic information (DNA methylation, histone modifications), nucleosome positioning, 3D structure of packed DNA, etc. So, it's just a digital model of a genome with a certain degree of accuracy, that is useful,  simple, and relatively affordable to obtain. It's a format commonly used for representing genomes in databases.

I quickly tested xz on this build (GCA_000001405.28) in the same way as in my reply (seqtk seq -U; xz -9 -T1) and got 667 MB, which is very close to fastqz-slow in the referenced table (657 MB). With general compressors, it's important to remove eols in the sequences, as well as to convert all bases to the upper case. Note that none of these steps modify the genome in any way – fastas can always be converted to a fixed width, and the upper/lower case info are usually just masked repeats.

As the table above says 770 MB (i.e., +103 MB) for xz with similar parameters (-e9), I assume the file probably wasn't in a form that would facilitate compression using general compressors, and the observed gain of specialized compressors came mainly from removing the eols/upper case conversion.

An additional experiment: I verified if xz could outcompete even fastqz-slow, and it indeed can. GCA_000001405.28 can be compressed down to 652 MB by 'xz -T1 -e --lzma2=preset=9,dict=512MiB' (again, after 'seqtk seq -U').

@Kamil S Jaron, cleaning is one reason, another is that it's just different genome build. Though as you see Karel obtain similar number with his method.

@Karel Břinda, by the way, it's good idea to tweak xz's dict (dict=512MiB), I should benchmark this setting too. Thanks.

Thanks both, although you are right Kirr, that in real-life situation you don't want to lose information about softmasked bases, the OP asked about the number of bites to store the human genome, therefore Karel's approach is completely valid (human genome does not have really softmasked bases, that's just our way how to add there one more layer of information).
Both of you, would you like to update the answers? I think this discussion is really interesting.

@Kamil S Jaron, I updated the answer. You're right that it's just one way of representing a genome, and also many genomes don't include soft-mask (e.g. prokaryote genomes).


text:
Create a genotype matrix from VCF file using bcftools
I have a vcf file that contains SNPs from chromosome 2 for 2000 samples. I want to create a genotype matrix from this vcf file.
The file is in zipped format as well so that I can use bcftools. I want to filter snips with

MAF>= 1%
more than 2 alleles.

I'm not quite sure why I have to filter for more than 2 alleles, is it common to do this if I'm trying to perform PCA analysis?
The instruction I'm using right now is:
bcftools query -e'AF<0.01' -e'AF>0.99' chr2.vcf.gz | \
  sed 's/0|0/0/g' | sed 's/0|1/1/g' | sed 's/1|0/1/g' | sed 's/1|1/2/g' | \
  grep -v "|" > extract.tab

How do I only extract the multi-allelic cases?

I have a way to do it using bcftools view -m2 chr2.vcf.gz but the results seem fishy.

How do I correct this?

What do you mean 'the results seem fishy'? Can you post a small sample input and output.


text:
plot correlation between several bam files
I have several mab files from medip-seq. they supposed to be replicates and I have to draw their correlation. First I used multiBamSummary from deeptools to merge all the bam files based on bi/bed. Now when I am using the plotCorrelation I do not get any file. the argument is :
plotCorrelation --corData myMultiBamSumary.npz --corMethod spearman --whatToPlot scatterplot --skipZeros --removeOutliers --plotFile Myplot --plotFileFormat pdf

Is that the correct strategy to pot the correlation? Is there any suggestion?
Thanks alot

That's a completely normal way to plot correlation between them. Spearman correlation was a wise choice on your part, since it's rank based and you don't have to care as much about normalization.
Update: It turns out that matplotlib is broken and I need to have a new deepTools release to work around that. That should be available sometime Sunday.

here is the file.  Thanks much   https://drive.google.com/open?id=1wKuHg3IfJWBCUwrInk4UZIA7mO80Pl7U

I suspect this is a matplotlib bug, I'll investigate.

There's a matplotlib bug that I need to program around. The fix is very simple and I'll release a new deepTools versions as soon as possible (sometime Sunday).

Thanks alot looking forward to it

3.4.2 is now on pypi and will be available via conda sometime today.

Much appreciate.


text:
Attempting to generate a heat map like in this figure
I have a data frame
Response TO CHEMOTHERAPY    Gene    CCF
Responders                  POLQ    1
Responders                 CDKN2A   1
Responders                  TP53    1
Responders              CCDC102B    0.906729075
Non-responders           CCDC102B   1
Non-responders           CCDC102B   1
Non-responders           CCDC102B   1
Non-responders             TSHZ3    1
Non-responders             PIK3CA   1
Responders                 SLIT2    1
Non-responders              SLIT2   1

I have two groups of patients; Responders to chemotherapy and non-responders. For each group I have calculated cancer cell faction (CCF) for a fixed set of genes so that an individual gene has several reads in each group
If CCF = 1 the gene clonal and if CCF < 1 gene is sub clonal, for a given gene for each group, I therefore would have the number of clonal and sub clonal for this gene, for example 2 clonal 7 sub clonal which would be 2/7 something like the below heat map

By your kindly code I produced this

Which is not well informative
This is another potential example heat map

Which show the clonal frequency for each group
By the way thanks for your kindness
I don't know how to use dput() but this is my complete data
https://www.dropbox.com/s/ol7m4nlkuyuy85a/c1.txt?dl=0

What have you tried? You could use geom_rect after setting x and y aes to Sample and Gene.

There is multiple examples of how to do heatmap with ggplot2 on SO, for example, you can start by taking a look here: https://stackoverflow.com/questions/8406394/how-to-produce-a-heatmap-with-ggplot2/8406774#8406774 or https://stackoverflow.com/questions/14290364/heatmap-with-values-ggplot2 or https://stackoverflow.com/questions/60162198/customizing-heatmap-appearance-using-ggplot-in-r/60163171#60163171 ... Try to reproduce those examples and then you will know how to deal with your data

Given the lack of effort from the OP, probably I shouldn't have wasted my time on this question with a lot of problems:
i) The question is not reproducible, the OP should have supplied their data with dput().
ii) The data as pasted above does not quite make sense to me, why would you have three readings from the same gene and even in such a case how should these replicates be represented in a heatmap?
iii) The question is not specific enough and should not go like "how can I generate a heatmap" and should be more like "in an attempt to generate a heatmap like in this figure I have tried this piece of code but ..."
Having said these, and after downvoting the question, here is a solution leveraging the ComplexHeatmap package:
library(data.table)
library(ComplexHeatmap)

my_data <-
"
Sample  Gene    CCF
Sample1 POLQ    1
Sample1 CDKN2A  1
Sample1 TP53    1
Sample1 CCDC102B    0.906729075
Sample1 CCDC102B    1
Sample1 CCDC102B    1
Sample1 CCDC102B    1
Sample2 TSHZ3   1
Sample2 PIK3CA  1
Sample2 SLIT2   1
Sample2 SLIT2   1
"

my_data <- fread(my_data)

my_data <- dcast.data.table(my_data,
                            formula = Gene ~ Sample,
                            value.var = "CCF",
                            fun=mean)

# dcast generates NA values as the provided data does not have values for each gene
my_data[is.na(my_data)] <- 0

# Heatmap() expects a matrix input
my_matrix <- as.matrix(my_data[,c("Sample1", "Sample2")])
rownames(my_matrix) <- my_data$Gene

Heatmap(my_matrix,
        cell_fun = function(j, i, x, y, width, height, fill) {
          grid.text(my_matrix[i, j], x, y, gp = gpar(fontsize = 10))})

This is how the data look after dcast, in which multiple entries per gene (my point ii above) is summarized with mean():
> head(my_matrix)
           Sample1 Sample2
CCDC102B 0.9766823       0
CDKN2A   1.0000000       0
PIK3CA   0.0000000       1
POLQ     1.0000000       0
SLIT2    0.0000000       1
TP53     1.0000000       0

Here is the resulting heatmap:

Sorry please don't close my question, I have edited that and I am certain I know what I want and what is the final goal


text:
Getting the number of weeks from a time span
I have a very big clinical data like
Date.Death  Date Of Operation
02/08/2015  27/11/2012
02/08/2015  27/11/2012

I want to know how many weeks each patient has lived after operation
I then tried
> t1 <- as.POSIXct("27/11/2012")
> t2 <- as.POSIXct("02/08/2015")
> as.double(difftime(t2,t1,unit="weeks"))
[1] -1317.571

Or
> span <- interval(as.POSIXct("27/11/2012"), as.POSIXct("02/08/2015"))
> span
[1] 0027-11-20 LMT--0002-08-20 LMT
> t <- as.period(span, unit="day")
> t
[1] "-9223d 0H 0M 0S"
> t / as.period(dweeks(1))
estimate only: convert to intervals for accuracy
[1] -1317.571

Both give wrong number of weeks
Manually calculating one by one is painful and error prone
You know anyway to get the number of weeks?
Thanks

I'm voting to close this question as off-topic because it has nothing to do with bioinformatics.

https://www.statmethods.net/input/dates.html

I think your error is located in your use of as.POSIXct which gives you wrong date formats because you did not provide a format in your function:
> as.POSIXct("27/11/2012")
[1] "27-11-20 LMT"

> as.POSIXct("27/11/2012", format = "%d/%m/%Y")
[1] "2012-11-27 PST"

So, you need to first convert your date in the appropriate date format and then you will be able to calculate the time interval.
Personally, I prefer to manipulate dates with lubridate package because the name of their function is a little bit more intuitive to me. For example, if you have a date to the format "Day/Month/Year", you can use dmy function.
With your example, you can do:
library(lubridate)
library(dplyr)

df <- data.frame(Date.Death = c("02/08/2015","02/05/2015"),
                 Date.Operation = c("27/11/2012","27/11/2012"))

df %>% mutate_all(~dmy(.)) %>%
  mutate(Diff_Weeks = interval(Date.Operation,Date.Death) / weeks(1),
                                      Diff_days = Date.Death - Date.Operation)

  Date.Death Date.Operation Diff_Weeks Diff_days
1 2015-08-02     2012-11-27   139.7143  978 days
2 2015-05-02     2012-11-27   126.5714  886 days
```

Sorry why inviting me? :( I googled a lot but I did not find the solution