text:
Choosing the best number of species for assay
I am doing amplicon sequencing of a virus across many different regions. Lets say I have 20k unique variants of the same virus that I put into my pcr assay and after sequencing and amplifications I am left with, 19k variants that appear. But many of the variants appear at really low read count and I don't know if that read is real or just noise. I have a threshold T that I say if a variant appears above T I'm counting it as real. Out of the 20k variants about only 15% of variants appear to be "real". I have some data on different input species counts (2k, 10k, 20k, 40k), and their corresponding fraction of "real".
My question is how do I determine which input variants count is best? I only looked at 4 different values and maybe the correct value isn't sampled, I obviously can't try everything between 2k and 40k. I think I want a balance between largest fraction of "real" variants and total number of variants.
Is there any research (or better yet code) that answers/discusses this problem?
Thanks in advanced.

@M__ Thanks for the comment adjusted my question.

Added doing amplicon sequencing across different regions

Thanks thats slightly clearer. There is an issue here on the nomenclature. 15+K viral species is 6 times more viruses than designated by the ICTV (International Committee on the Taxonomy of Viruses), i.e. more viruses than the total number of species ever declared. Could you kindly define the term species? What you have stated is not taxonomically possible.

I see, using species not exactly clearly. They are all variants of the same virus. Will adjust

The underlying answer to the question is have you experimentally controlled for in vitro artefacts? Thats taking a cloned virus and performing the same RT-PCR. If not thats makes things very difficult. If there isn't a cloned virus thats a problem. You might need to ask your wet-lab colleagues what old school "cloning a virus" means, or Biology SE will explain this. The reason this is important is given below.
The first step to the question is simply to plot a histogram against copy number. What you will find is a heavily skewed distribution where the copy number of one variant will represent >40% (at a wild guess) of all the viral copy number. Moreover, a handful (less than 100) of variants will comprise 90% of all the copy number. Once that is established it is much easier to assess what is a true variant. What I'm saying is the difference is one log minimum between what you think you are looking at and what is there.
The sources of artefact:

RT-PCR artefacts, particularly if there is an RT step: Taq polymerase alone has an error rate of 1 in 1000 and RT has a very high error rate, worse it's the first step.

What DNA polymerase was used?
The key issue is separating the RT-PCR artefacts from naturally occurring RNA polymerase errors in situ is not possible without a control. If a frequency analysis cut-off strategy is used it needs to be an ultra-conservative threshold, i.e. I mean ridiculously high and this means lot of variants are going to be dumped.
Why? Again it's back to the problem whether the variants were generated in vitro, i.e. artefacts, or in situ, i.e. resulting from RNA polymerase. Generating a cut-off within a skewed distribution is not trivial per se but extremely difficult when the ratio between artefact and genuine variant has not been controlled.

There appears some sort of control. Maximising the threshold: thats a function of the control, so the control needs to accommodate a model of the different numbers of input variables; it will not be a linear function. Basically what input will give me the best representation of non-artifactual variants, that is what the control has to find out. Thus, for every input, you want a control.
Again, if the control is lacking, in any context, then the rule of thumb is to be ultra-conservative, not trying to maximise the number of variants, but minimise them, particularly if the control isn't perfect. Most virological reviewers raise concern with the numbers been reported here on the information given (obviously we don't know the virus).

NGS sequencing is another source of artefact - I should have mentioned there are some famous disasters here. I'm sure you NGS is really good, but a small number of investigators have deposited a disproportionate number of artefact ridden sequences - that everyone knows about.

from comments (again)
A control of 4 is tiny in data science especially a GLM approach. However, other studies would have been performed - possibly lots of them if this was e.g. HIV - and controls would have been run. There are two major strategies (really 3) for biologically cloning a virus. Really you'll need 100+ controls, that may not be feasible but as many within 4 to 100 as are available. GLM would be a good approach IMO. It would need stating that this does not truly control the rate of artefacts, but it's much closer.

Final comments ... don't want to assess thresholds.
I gave this question some thought and I think it would work via machine learning (what I specialise in). However, it may be difficult to find appropriate published controls and/or augmentation. Personally, I'm satisfied this would work if done correctly, but implementation is fairly advance ML. Given the amount of time I've spent on thinking this question through I wouldn't use an alternative approach, because I think it would be valid. Augmentation BTW is a generic term in machine learning and deep learning to increase the size of the training set.

Yes, they have controlled for some in vitro artefacts. I have plotted histograms of frequencies and each of the input variants. I got the expected picture you described, one variant taking a large portion of the reads a small number with a higher read count and a lot of low read count ones. I don't think I want to find a cutoff for the distribution, let's say I have one. I'm looking for a way to choose the number of input variants that maximizes observations over that threahold

Following up on your edit. For each input I have a control but I only have 4 inputs. It is too expensive/time consuming to sequence too many more. Do I just use the controls on the 4 inputs I have? Or is there some way I can simulate behavior for unknown inputs?

ok I think, I'm starting to pick up what you are saying. I have a few different options for what is considered a good observation. A hard threshold, item using some calculus (knee/inflection), as well as some methods based on spike ins (these are ideal probably).
But lets say for each input we can compute the fraction of observations above the threshold. Then fit a GLM to these 4 points? 4 seems a bit small for a GLM right?

Ok I think we are slightly talking about different things. I don't want to evaluate the controls or use them to select a threshold. I want to decide how many unique variants of my virus to sequence at once. Is that what your comment is describing and I'm just missing something?

Ok @TheNumber23 answered in comments above. Upvotes are encouraged.


text:
Extract data from genome
This question was also asked on Biostars 3 months ago, i.e. a long time - so its perfectly reasonable for the question to be asked here
I started studying bioinformatics and I enjoy it a lot (I'm a software developer with around 15 years of experience but mostly in backend, scripting, architecture, and security) and I want to do, as a project, a database regarding natural compounds. When I saw that on NCBI there are a lot of plant genomes I thought, that, with software, somehow we can "guess" (probabilities, not certainties) what a plant can contain or not.
I'm wondering, if, by bioinformatic means only, it is possible to "extract" genes and, after that proteins and natural compounds from plant genome data that is available on NCBI (please see https://www.ncbi.nlm.nih.gov/datasets/genome/GCA_029618835.1/ as an example)? Example: just by processing the data from the link above we can get to Flavoxanthin or any other substance.
I was wondering if a kind of "de-novo" annotation can be made (if that's the correct term).
What other useful information from a chemical compound perspective can be found from digging on the genome itself?
Thank you!

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

Why do you want to do this? There are lots of different tools; it would help to be more specific in order to make the answers actually useful.

@gringer .. you are totally right, I put full details on another platform when I saw that the interest here was almost none:  https://www.biostars.org/p/9566425/

Thanks for the additional information, but your additional details still describe a generic problem, and are therefore unlikely to get good answers. What is the specific problem that you are trying to solve?

This sounds like you are wanting to perform a machine learning calculation. However, for this question leveraging KEGG, Kyoto Encyclopaedia of Genome and Gene and building an SQL database around it is much more effective. It represents around 60% of all functional genes. KEGG is the best gene annotation out there IMO, particularly for flavonoids etc... Fantastic and all their linking pathways, everything links to everything, just amazing. KEGG represents an enormous volume of work, much of it being manually crunching publications for each gene.
The entire database of KEGG keys is downloadable. Connecting a gene to a class of function is easy via KEGG the relationships are defined its just SQL-ing the whole thing together (querying all the tables). Again, it's totally and completely brilliant. Sure they'll be stuff out there has SQL-ed KEGG key tables, but if your a developer its better to build it independently - this way you're certain whats in there and what its doing.
However, linking raw sequence data to a KEGG gene code is more challenging. lastKOALA and GhostKOALA is ... yeah well an approach no more. There's a main stream approach but with current limitations in implementation. KEGG however isn't limiting, it's linking a gene to a KEGG id is the challenge. Getting hold of the raw sequence for KEGG would be great, fantastic in fact, but I'm yet to find it.

I've used Augustus stand alone (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1160219/) to find genes in a given eukaryotic genome. I agree you can then use things like KofamSCAN to categorize the genes against the KEGG database. If you want to compare against the NCBI database, diamond blast (https://github.com/bbuchfink/diamond) is probably the quickest way to go :)

source: Stack Exchange

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

sort: votes


text:
How to open .weir.fst files witout losing data?
I have .vcf file with over 78 000 variations in 95 experiment subjects. I used vcftools to extract fixation index score for general differentiation between Case and Control groups.
./vcftools --vcf SCH_freebayes_sort_ostE_SnpEff_SnpSift.vcf --weir-fst-pop Case.txt --weir-fst-pop Control.txt --out FstExon_Case_vs_Control

VCFtools - 0.1.17
(C) Adam Auton and Anthony Marcketta 2009

Parameters as interpreted:
    --vcf SCH_freebayes_sort_ostE_SnpEff_SnpSift.vcf
    --weir-fst-pop Case.txt
    --weir-fst-pop Control.txt
    --keep Case.txt
    --keep Control.txt
    --out Fst_TestData_nochange

Warning: Expected at least 2 parts in INFO entry: ID=AF,Number=A,Type=Float,Description="Estimated allele frequency in the range (0,1]">
Warning: Expected at least 2 parts in INFO entry: ID=PRO,Number=1,Type=Float,Description="Reference allele observation count, with partial observations recorded fractionally">
Warning: Expected at least 2 parts in INFO entry: ID=PAO,Number=A,Type=Float,Description="Alternate allele observations, with partial observations recorded fractionally">
Warning: Expected at least 2 parts in INFO entry: ID=SRP,Number=1,Type=Float,Description="Strand balance probability for the reference allele: Phred-scaled upper-bounds estimate of the probability of observing the deviation between SRF and SRR given E(SRF/SRR) ~ 0.5, derived using Hoeffding's inequality">
Warning: Expected at least 2 parts in INFO entry: ID=SAP,Number=A,Type=Float,Description="Strand balance probability for the alternate allele: Phred-scaled upper-bounds estimate of the probability of observing the deviation between SAF and SAR given E(SAF/SAR) ~ 0.5, derived using Hoeffding's inequality">
Warning: Expected at least 2 parts in INFO entry: ID=AB,Number=A,Type=Float,Description="Allele balance at heterozygous sites: a number between 0 and 1 representing the ratio of reads showing the reference allele to all reads, considering only reads from individuals called as heterozygous">
Warning: Expected at least 2 parts in INFO entry: ID=ABP,Number=A,Type=Float,Description="Allele balance probability at heterozygous sites: Phred-scaled upper-bounds estimate of the probability of observing the deviation between ABR and ABA given E(ABR/ABA) ~ 0.5, derived using Hoeffding's inequality">
Warning: Expected at least 2 parts in INFO entry: ID=RPP,Number=A,Type=Float,Description="Read Placement Probability: Phred-scaled upper-bounds estimate of the probability of observing the deviation between RPL and RPR given E(RPL/RPR) ~ 0.5, derived using Hoeffding's inequality">
Warning: Expected at least 2 parts in INFO entry: ID=RPPR,Number=1,Type=Float,Description="Read Placement Probability for reference observations: Phred-scaled upper-bounds estimate of the probability of observing the deviation between RPL and RPR given E(RPL/RPR) ~ 0.5, derived using Hoeffding's inequality">
Warning: Expected at least 2 parts in INFO entry: ID=EPP,Number=A,Type=Float,Description="End Placement Probability: Phred-scaled upper-bounds estimate of the probability of observing the deviation between EL and ER given E(EL/ER) ~ 0.5, derived using Hoeffding's inequality">
Warning: Expected at least 2 parts in INFO entry: ID=EPPR,Number=1,Type=Float,Description="End Placement Probability for reference observations: Phred-scaled upper-bounds estimate of the probability of observing the deviation between EL and ER given E(EL/ER) ~ 0.5, derived using Hoeffding's inequality">
Warning: Expected at least 2 parts in INFO entry: ID=TYPE,Number=A,Type=String,Description="The type of allele, either snp, mnp, ins, del, or complex.">
Warning: Expected at least 2 parts in INFO entry: ID=TYPE,Number=A,Type=String,Description="The type of allele, either snp, mnp, ins, del, or complex.">
Warning: Expected at least 2 parts in INFO entry: ID=TYPE,Number=A,Type=String,Description="The type of allele, either snp, mnp, ins, del, or complex.">
Warning: Expected at least 2 parts in INFO entry: ID=TYPE,Number=A,Type=String,Description="The type of allele, either snp, mnp, ins, del, or complex.">
Warning: Expected at least 2 parts in INFO entry: ID=TYPE,Number=A,Type=String,Description="The type of allele, either snp, mnp, ins, del, or complex.">
Warning: Expected at least 2 parts in INFO entry: ID=CIGAR,Number=A,Type=String,Description="The extended CIGAR representation of each alternate allele, with the exception that '=' is replaced by 'M' to ease VCF parsing.  Note that INDEL alleles do not have the first matched base (which is provided by default, per the spec) referred to by the CIGAR.">
Warning: Expected at least 2 parts in FORMAT entry: ID=GQ,Number=1,Type=Float,Description="Genotype Quality, the Phred-scaled marginal (or unconditional) probability of the called genotype">
Warning: Expected at least 2 parts in FORMAT entry: ID=GL,Number=G,Type=Float,Description="Genotype Likelihood, log10-scaled likelihoods of the data given the called genotype for each possible genotype generated from the reference and alternate alleles given the sample ploidy">
Warning: Expected at least 2 parts in INFO entry: ID=WTD,Number=0,Type=Flag,Description="Is Withdrawn by submitter If one member ss is withdrawn by submitter, then this bit is set.  If all member ss' are withdrawn, then the rs is deleted to SNPHistory">
Warning: Expected at least 2 parts in INFO entry: ID=WTD,Number=0,Type=Flag,Description="Is Withdrawn by submitter If one member ss is withdrawn by submitter, then this bit is set.  If all member ss' are withdrawn, then the rs is deleted to SNPHistory">
Warning: Expected at least 2 parts in INFO entry: ID=TOPMED,Number=.,Type=String,Description="An ordered, comma delimited list of allele frequencies based on TOPMed, starting with the reference allele followed by alternate alleles as ordered in the ALT column. The TOPMed minor allele is the second largest value in the list.">
Warning: Expected at least 2 parts in INFO entry: ID=TOPMED,Number=.,Type=String,Description="An ordered, comma delimited list of allele frequencies based on TOPMed, starting with the reference allele followed by alternate alleles as ordered in the ALT column. The TOPMed minor allele is the second largest value in the list.">
Warning: Expected at least 2 parts in INFO entry: ID=WGT,Number=1,Type=Integer,Description="Weight, 00 - unmapped, 1 - weight 1, 2 - weight 2, 3 - weight 3 or more">
Warning: Expected at least 2 parts in INFO entry: ID=WGT,Number=1,Type=Integer,Description="Weight, 00 - unmapped, 1 - weight 1, 2 - weight 2, 3 - weight 3 or more">
Warning: Expected at least 2 parts in INFO entry: ID=WGT,Number=1,Type=Integer,Description="Weight, 00 - unmapped, 1 - weight 1, 2 - weight 2, 3 - weight 3 or more">
Warning: Expected at least 2 parts in INFO entry: ID=WGT,Number=1,Type=Integer,Description="Weight, 00 - unmapped, 1 - weight 1, 2 - weight 2, 3 - weight 3 or more">
Warning: Expected at least 2 parts in INFO entry: ID=MTP,Number=0,Type=Flag,Description="Microattribution/third-party annotation(TPA:GWAS,PAGE)">
Warning: Expected at least 2 parts in INFO entry: ID=NOC,Number=0,Type=Flag,Description="Contig allele not present in variant allele list. The reference sequence allele at the mapped position is not present in the variant allele list, adjusted for orientation.">
Warning: Expected at least 2 parts in INFO entry: ID=CAF,Number=.,Type=String,Description="An ordered, comma delimited list of allele frequencies based on 1000Genomes, starting with the reference allele followed by alternate alleles as ordered in the ALT column. Where a 1000Genomes alternate allele is not in the dbSNPs alternate allele set, the allele is added to the ALT column. The minor allele is the second largest value in the list, and was previuosly reported in VCF as the GMAF. This is the GMAF reported on the RefSNP and EntrezSNP pages and VariationReporter">
Warning: Expected at least 2 parts in INFO entry: ID=CAF,Number=.,Type=String,Description="An ordered, comma delimited list of allele frequencies based on 1000Genomes, starting with the reference allele followed by alternate alleles as ordered in the ALT column. Where a 1000Genomes alternate allele is not in the dbSNPs alternate allele set, the allele is added to the ALT column. The minor allele is the second largest value in the list, and was previuosly reported in VCF as the GMAF. This is the GMAF reported on the RefSNP and EntrezSNP pages and VariationReporter">
Warning: Expected at least 2 parts in INFO entry: ID=CAF,Number=.,Type=String,Description="An ordered, comma delimited list of allele frequencies based on 1000Genomes, starting with the reference allele followed by alternate alleles as ordered in the ALT column. Where a 1000Genomes alternate allele is not in the dbSNPs alternate allele set, the allele is added to the ALT column. The minor allele is the second largest value in the list, and was previuosly reported in VCF as the GMAF. This is the GMAF reported on the RefSNP and EntrezSNP pages and VariationReporter">
Warning: Expected at least 2 parts in INFO entry: ID=CAF,Number=.,Type=String,Description="An ordered, comma delimited list of allele frequencies based on 1000Genomes, starting with the reference allele followed by alternate alleles as ordered in the ALT column. Where a 1000Genomes alternate allele is not in the dbSNPs alternate allele set, the allele is added to the ALT column. The minor allele is the second largest value in the list, and was previuosly reported in VCF as the GMAF. This is the GMAF reported on the RefSNP and EntrezSNP pages and VariationReporter">
Warning: Expected at least 2 parts in INFO entry: ID=MUT,Number=0,Type=Flag,Description="Is mutation (journal citation, explicit fact): a low frequency variation that is cited in journal and other reputable sources">
Warning: Expected at least 2 parts in INFO entry: ID=SAO,Number=1,Type=Integer,Description="Variant Allele Origin: 0 - unspecified, 1 - Germline, 2 - Somatic, 3 - Both">
Warning: Expected at least 2 parts in INFO entry: ID=SAO,Number=1,Type=Integer,Description="Variant Allele Origin: 0 - unspecified, 1 - Germline, 2 - Somatic, 3 - Both">
Warning: Expected at least 2 parts in INFO entry: ID=SAO,Number=1,Type=Integer,Description="Variant Allele Origin: 0 - unspecified, 1 - Germline, 2 - Somatic, 3 - Both">
Warning: Expected at least 2 parts in INFO entry: ID=SSR,Number=1,Type=Integer,Description="Variant Suspect Reason Codes (may be more than one value added together) 0 - unspecified, 1 - Paralog, 2 - byEST, 4 - oldAlign, 8 - Para_EST, 16 - 1kg_failed, 1024 - other">
Warning: Expected at least 2 parts in INFO entry: ID=SSR,Number=1,Type=Integer,Description="Variant Suspect Reason Codes (may be more than one value added together) 0 - unspecified, 1 - Paralog, 2 - byEST, 4 - oldAlign, 8 - Para_EST, 16 - 1kg_failed, 1024 - other">
Warning: Expected at least 2 parts in INFO entry: ID=SSR,Number=1,Type=Integer,Description="Variant Suspect Reason Codes (may be more than one value added together) 0 - unspecified, 1 - Paralog, 2 - byEST, 4 - oldAlign, 8 - Para_EST, 16 - 1kg_failed, 1024 - other">
Warning: Expected at least 2 parts in INFO entry: ID=SSR,Number=1,Type=Integer,Description="Variant Suspect Reason Codes (may be more than one value added together) 0 - unspecified, 1 - Paralog, 2 - byEST, 4 - oldAlign, 8 - Para_EST, 16 - 1kg_failed, 1024 - other">
Warning: Expected at least 2 parts in INFO entry: ID=SSR,Number=1,Type=Integer,Description="Variant Suspect Reason Codes (may be more than one value added together) 0 - unspecified, 1 - Paralog, 2 - byEST, 4 - oldAlign, 8 - Para_EST, 16 - 1kg_failed, 1024 - other">
Warning: Expected at least 2 parts in INFO entry: ID=SSR,Number=1,Type=Integer,Description="Variant Suspect Reason Codes (may be more than one value added together) 0 - unspecified, 1 - Paralog, 2 - byEST, 4 - oldAlign, 8 - Para_EST, 16 - 1kg_failed, 1024 - other">
Warning: Expected at least 2 parts in INFO entry: ID=PM,Number=0,Type=Flag,Description="Variant is Precious(Clinical,Pubmed Cited)">
Warning: Expected at least 2 parts in INFO entry: ID=AF_popmax,Number=A,Type=Float,Description="Maximum allele frequency across populations (excluding samples of Ashkenazi, Finnish, and indeterminate ancestry)">
Warning: Expected at least 2 parts in INFO entry: ID=AF_popmax,Number=A,Type=Float,Description="Maximum allele frequency across populations (excluding samples of Ashkenazi, Finnish, and indeterminate ancestry)">
Keeping individuals in 'keep' list
After filtering, kept 95 out of 95 Individuals
Outputting Weir and Cockerham Fst estimates.
    Fst: Only using diploid sites.
Weir and Cockerham mean Fst estimate: -0.00031927
Weir and Cockerham weighted Fst estimate: -0.00047509
After filtering, kept 78758 out of a possible 78758 Sites
Run Time = 4.00 seconds

But when I tried to open .weir.fst files with notepad or excel, I lost over 9 000 rows of data. I am sure that there is problem with how I am opening the .weir.fst files.
Please, how are these files opened properly?

Please format your error message/console output properly.

@RamRS I added full log file. Please If you want to help me with formatting my question, be more direct.

I was as direct as can be. What do you think I should do differently?

@RamRS I apologise for my comment. I did not mean to offend you. I was only seriously frustrated and I felt hopeless. This problem took me hours. I do not know, where I made mistake.

I'm not offended, Lukas. I was just wondering what I could have done differently.

Also, don't open files with notepad or Excel, try unix utilities or if you're on Windows, try notepad++

@RamRS Sorry, but I am new contributor. You can tell mi dircetly that I could my log message post fully and edit it little more. But still this is just pointless semantics.

@RamRS Thank you.

Oh this is fairly easy. The file size has exceed the capacity of Excel (50000 rows??). This doesn't affect downstream analysis because the rows are there you just can't visualise them.
BBEdit will open it if this is Mac (free from AppStore - really good BTW). If this is Linux then vi or vim. A Windows version of BBEdit is Notepad++, please not Notepad.
You can request the number of lines (bash Linux) via
wc weir.fst

This will return the number of lines in the file, if it's >50000, i.e. beyond the number of the Excel limit, that's your answer. You just need to find an alternative viewer

Thank you. Now I am really embaresed, that I opened those files in that way.

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


text:
Biopython / ETE3 get subtree of a newick tree
I can not find a good example how to extract a subtree from a phylogenetic tree with either library.
I can load newick files etc, but not sure how to search for a node and extract all children with the right topology. Already thought of parsing the newick file and extract on raw data but there must be a way...
How?
This is a sample tree
(A:0.1,B:0.2,(C:0.3,(D:0.4, (E:0.6, F:0.6))));

and I would like to get everything below C, so D,E,F . I see the ETE3 get_dependants function but it returns a list, not a tree... ? Please correct me if I'm wrong.
I understand and do not expect anyone to write code for me, but just a pointer to the right function ete3 would be enough. I can get all children for C but not sure how to preserve the subtree relationships with that command ete3 : get_dependants().

If we are talking pointers only then the pseudo code I use the ETE3 method
t = Tree((A:0.1,B:0.2,(C:0.3,(D:0.4, (E:0.6, F:0.6))));)
clade = t.get_common_ancestor( .... ) # name the descendants

Thats it, what its doing is cool, but the code is simple

From comments "clade" is the subtree defined by the descendants of the common ancestor

The OPs solution
subtree = node.up.up.up

This is fine if you've a fixed tree, rather than trees that might have a changed topology (for example bootstrap replicates). For example if the paraphyly changes in another "bootstrap replicate', up.up.up may not be the clade of interest. get_common_ancestor ensures the code is generalised. For example if the topology changes the approach above will still identify the correct subtree.
Thus if on one tree you need up.up.up, another tree up.up and a different tree up.up.up.up the above solution solves all scenarios.

so the clade object will be of type Tree ? Need to try in the terminal. Thanks!

Addressed above

Yes, the up.up.up chain has to be adapted to the location of the searched node because it could error out if its hardcoded

Actually ended up doing this to get the neighborhood of a node
tree = Tree(newick_file_path, format=1)
node = tree.search_nodes(name=node_name)[0]
subtree = node.up.up.up

but get_common_ancestor worked too, thanks!

Thanks for your answer. I've responded above about this approach. This can definitely work in some situations, as long as you understand why thats cool.


text:
scRNA: What are good dimensionality reduction/clustering parameters to get biologically plausible groupings?
I've got a moderately large set of PBMCs, over 1M cells. That means I can't easily do a grid search of dimensionality reduction/clustering parameters/methods. Some examples results I'm getting with scanpy default parameters include PCA -- a big, dense grouping with lobes radiating outwards, and UMAP -- "better" groupings in some sense, but many many tiny clusters of cells, more than are biologically plausible (hundreds) as well.
Basically I'm looking for best dimred/clustering practices for such a dataset, other than what's on offer from scanpy/Seurat. My metrics such as they are are that the major cell types should be clearly separated, and that the overall plot have some biological plausibility in terms of number of visually distinct clusters.
I'm aware of the Xiang et al Frontier Genetics 2021 paper, but the number of cells their example datasets use is much smaller AFAIK. They concluded that tSNE is the best choice of method according to their metrics, though UMAP also has some stability advantages. The Chari and Pachter PLOS Comp. Bio. 2023 paper raises some good points (a highly dimensional space can't really accurately have its relations preserved in 2 dimensions, and there are other issues of noise, interpretability, and so on. They recommend clustering from the higher dimensional data, which I do agree with. A different paper rec could be great.
PS: SCVI seems like a good option, I'm playing with that at the moment.

Hi M, thanks for the reply. The first paper I mentioned covers several dimensionality reduction methods and concludes that for the older data they tested, tSNE is best. The second one covers some of the theory and recommends perhaps side stepping dimred entirely.

As for metric, basically I'm looking for separation of major cell types and a plausible overall number of clusters. PCA for instance does separate some PBMC cell types out but not most. UMAP does separate out certain cell types (B cells are well separated for instance) but not others and seems to have lots of spurious tiny clusters.

Yes, I'll edit it in a bit -- will write up a less pithy short summary of the second paper.

I always wonder what 1M cells can answer what 100k cannot. Anyway, what you have to find out based on published datasets of that magnitude is which approaches scale to that number of cells. I never had to do that myself since our datasets are always way smaller. I would consider not to analyse this per-cell but rather try to find some sort of meta-cell approach where you bin cells into metas, to reduce the computational burden. I would not believe that with 1M cells one really benefits from pure single-cell analysis. Maybe some sort of clustering to get "pseudobulks" and then use these.

@M, just SCVI at the moment. @ATpoint, it's 1M+ cells across many individuals, and I do agree that cluster level comparisons are a good way to make the data tractable -- hence the question, really. Point well taken to compare to other similar datasets.


text:
Pathway Enrichment Analysis for one specific pathway?
I have a gene set containing a set of 1500 genes in which I did differential gene expression using limma voom pipeline. After this step I wanted to analyze a specific set of pathways the metabolism set from KEGG database.
hs_kegg_df <- msigdbr(species = "Homo sapiens") %>%
 dplyr::filter(
gs_cat == "C2", # This is to filter only to the C2 curated gene sets
gs_subcat == "CP:KEGG" # This is because we only want KEGG pathways
)

Then since I only want to analyze metabolism I did this (keep only the pathways envolved in metabolism):
hs_kegg_df_with_gs_source <- subset(hs_kegg_df, grepl("^hsa00|^hsa01|^hsa02", gs_exact_source))

My problem is that, for example, for the pathway  KEGG arginine and proline metabolism the set contains 54 genes in which my kit only covers 4, and these for 4 genes are Differentially expressed yet when doing: ( I decided to do ORA because with 1500 genes GSEA is not recommended).
em <- enricher(genes, minGSSize=1, maxGSSize = 20,  universe = names(gene_list), TERM2GENE=hs_kegg_df,pAdjustMethod = "fdr")

This set does not appear with significant p value. I also tried passing to the enricher just that pathway in particular, the p value is 1. Because in hypergeometric test, the background genes are the same as the as DGE genes in all pathways giving this way 1-0=1.
How can I analyze just one pathway in particular given this specific conditions?
These are the 4 genes DGE and included in my kit:
  gs_name                              gene_symbol
  <chr>                                <chr>
1 KEGG_ARGININE_AND_PROLINE_METABOLISM ALDH2
2 KEGG_ARGININE_AND_PROLINE_METABOLISM ARG2
3 KEGG_ARGININE_AND_PROLINE_METABOLISM NOS1
4 KEGG_ARGININE_AND_PROLINE_METABOLISM ODC1

And this are all the 59.
here all the genes dge and all the genes in pathway:
 genes_dge <- c("ALDH2", "ARG2", "NOS1", "ODC1")
 all_genes <- c("ACY1", "AGMAT", "ALDH18A1", "ALDH1B1", "ALDH2", "ALDH3A2", "ALDH4A1", "ALDH7A1", "ALDH9A1", "AMD1", "AOC1", "ARG1", "ARG2", "ASL", "ASS1", "AZIN2", "CKB", "CKM", "CKMT1A", "CKMT1B",s"CKMT2", "CPS1", "DAO", "GAMT", "GATM", "GLS", "GLS2", "GLUD1", "GLUD2", "GLUD2", "GLUL", "GOT1", "GOT2", "LAP3", "MAOA", "MAOB", "NAGS", "NOS1", "NOS2", "NOS3", "OAT", "ODC1", "OTC", "P4HA1", "P4HA2", "P4HA3", "PRODH", "PRODH2", "PYCR1", "PYCR2", "PYCR3", "PYCR3", "SAT1","SAT2", "SMS", "SRM")

Thank you for the answer @M__  its samples from two different subtypes of human cancer. And its "C2".

Two questions

Significance
Lack of consistency between enzymes (log2) within the arginine and proline metabolism pathway.

.  1. We've dealt with significance in a past post on Limma (it was not voom though). I'll try and find the post and link it here.

The answer is straight forward: whilst all four enzymes include Arginine and proline metabolism pathway

only ARG and ODC1 share exactly the same pathways. Thus these should give very similar results.

NOS1 is a key enzyme in 6 different pathways.

ALDH2 is involved in 13 different pathways.

Therefore NOS1 and ALDH2 would have "confounders" (from an experimental perspective, not host metabolism) because they are more complex than singly arginine and proline metabolism pathways.
The other factor is in cancer patients the transcriptomics will get messed up, ALDH2 could go all over the place.
Prediction  ARG = ORD < NOS1 < ALDH2
**ALDH2** aldehyde dehydrogenase 2
Pathways
Glycolysis / Gluconeogenesis
Ascorbate and aldarate metabolism
Fatty acid degradation
Valine, leucine and isoleucine degradation
Lysine degradation
**Arginine and proline metabolism**
Histidine metabolism
Tryptophan metabolism
beta-Alanine metabolism
Glycerolipid metabolism
Pyruvate metabolism
Pantothenate and CoA biosynthesis

**ARG** arginase
Pathways
**Arginine and proline metabolism**
Glutathione metabolism

**NOS1** nitric oxide synthase
Pathways
Arginine biosynthesis
**Arginine and proline metabolism**
Calcium signaling pathway
Phagosome
Apelin signaling pathway

**ODC1** ornithine decarboxylase 1
Pathways
**Arginine and proline metabolism**
Glutathione metabolism


text:
Integrate bulk RNA and ATAC-seq genes
I have up and down-regulated genes from bulk RNA result from DeSeq2. Genes that are differential in chromatin accessibility between two conditions from DiffBind. I would like to find the genes overlap between bulk RNA and ATAC-seq assay as these papers but don't know the code to get it. Would you please have a suggestion? Thank you so much  .

paper 1
paper 2

in silico DEG details
I used intersect function in R which can do the same job. Which I am not sure is the experiment try to find transcription factors that cause the difference in phenotype between 2 conditions. DEG in bulk-RNA from DeSeq2 (log2 fold change >2.5 or log2 fold change < -2.5) intersect with genes differently in chromatin accessibility from DiffBind (log2 fold change >2.5 or log2 fold change < -2.5). I narrowed down to 142 genes.
For atac-seq result, after applying filter log2fold > 2.5 or log2fold < -2.5 and p value < 0.00001, I still got around 12k genes that are different between 2 conditions. Meanwhile bulk-RNA after apply filter has only around 300 genes. The number of genes in human and pig are similar so I am not sure about the number of gene I have in atac-seq
What do you think about this approach because I don't use MEME or TOMTOM?

I share with you what the author did.

" In order to screen new transcription factors by combining with star transcription factors. First, we used MEME software to predict the motif of Duroc pigs' peaks. Second, we compared the transcription factor database with TOMTOM to predict all transcription factors that may bind to the chromatin open regions . Third, we found all possible binding sites of MYOD, MYOG, MYF5, MYF6 and MEF2 families and did a union. A total of 99 sites were found.
"Fourthly, we counted other transcription factors that may bind to these 99 sites and sorted them by the number of binding sites. All of the above statistics were manually counted."

Will read more to give you a good reply. Thank you so much

You requested I make a response. Okay, I'm reluctant to do this firstly, all my  eukaryotes genomics is out of date except for immunology stuff and secondly, this is a hardcore eukaryote gene expression question. I am really not sure about the pig-human side of this experiment.
What you have is an extremely good data set. What you appear to be doing is assessing the complexity of chromatin based gene expression, using very hardcore data. What we don't know is the variation e.g. in experiments leading to the RNA-seq. I would guess the RNA-seq is varying and being compared against atac-seq.
I have not read the papers, read the abstracts you posted.
MEME and TOMTOM don't seem to be impressive in comparison to an absolute biological experiment RNA-seq vs. atac-seq. Predictions are never perfect (the very, very best predictions are 80%) and you have an absolute data set - why even be concerned about predictions based analysis?
I do know DEG and the approach you are using is very conservative. It's good to visualise this on a volcano plot. Personally I just use a probabilistic threshold and then assess the log2 for each gene group. You are using both -/+2.5 log2 and <1E-6 . Personally I would be concerned with the stringency and also look at results with e.g. just <1E-6 or <1E-6 and -/+ 0.5log2. Keep in mind I don't do eukaryotes.
The threshold criteria will be the number 1 issue that affects your results and that is the most important information I can give. So do your approach and also try a more relaxed threshold. Personally I wouldn't relax the p<1E-6, but I would relax the log2. The reasons for this are a bit complicated, but they are statistically sound.
Your intersect approach seems fine, viz.
intersect(object1, object2)

@haci's method is cool, the above does the job too. The question is the genes vs thresholds and the experimental set ups and then the comparison between experimental setups, i.e. which genes are switching under chromatin control.
What I am not clear about is the experimental design and how you are performing DSeq2,

I assume you are using an internal control for RNAseq and a comparable one for atac-seq. I assume you are not comparing atac-seq verse RNAseq in the DEG.
The alternative is atac-seq vs. RNAseq DEG and the comparison is pig versus human. I honestly dunno.

Anway the caveat I would add is an intersect SHOULD NOT ignore whether it is +log2 or -log2 ... you need to do four intersects

+ve verse +ve
-v verse -ve
+ve verse -ve
-ve verse +ve

Again its four separate intersects.
Good luck looks a really interesting experiment and analytics look simple but powerful in context.

Thanks @M__  so much for a very detailed reply! Each RNA-seq and atac-seq has control vs diseased. I can do a volcano plot. So I compared RNA-seq control vs RNA-seq diseased. ATAC seq control vs ATAC seq diseased. Then I try to figure out which genes/transcription factors may be the cause of the difference in the phenotype between two conditions. I see papers do scRNA-seq integrative ATAC-seq but not many papers do bulk RNA-seq integrative ATAC-seq. Maybe I should intersect up-regulated genes vs open chromatin region genes and down-regulated genes vs close chromatin region genes.

Yes, I tried to do it at the beginning but I need two more point to allow upvote. Sorry for that! When I have enough reputation, I will upvote all your answers right away!

I am allowed to upvote now .

I've borrowed one of your smilies, thank you @chris

You are welcome .

If I am getting your question right, all you need is https://bioinfogp.cnb.csic.es/tools/venny/. I assume that your two lists corresponding to RNA-seq and ATAC-seq use the same set of gene identifiers (gene symbol, Entrez ID, Ensembl ID etc).
The Results section is populated with entries that belongs to the corresponding "region" upon clicking on the regions of the Venn diagram.
Personally I have never used it but my wet-lab colleagues use it all the time.


text:
How to get a GISAID account? I registered months ago, still no reply!
Inspired by amateur variant hunters, I would like to join the Pango lineage proposal community and help contribute to variant surveillance.
However, I cannot seem to get access to GISAID, the platform that exclusively hosts a large share of SARS-CoV-2 sequences from around the world.
I applied/registered for an account back in summer 2022 (almost a year ago) using this form: https://gisaid.org/register/
But since then, nothing has happened. I have send an email and filled out the contact form on the website but never heard back.
I am very confused how GISAID can say they are a "public" database, while at the same time not giving accounts to people like me, not even responding to emails.
How does everyone else get access? Is there a secret email you need to contact? Or do you need to pay? I've read here which describes the founder of GISAID in detail.
Any help would be greatly appreciated!

I've got a GISAID account and it's subject to an approval process. It was pretty simple and is awarded quickly. There are definitely no fees needed. I personally didn't see the advantages over NCBI's really very impressive SARS-CoV-2 database system and thats what I'd recommend.
Their resource is here:
https://www.ncbi.nlm.nih.gov/labs/virus/vssi/#/virus?SeqType_s=Nucleotide&VirusLineage_ss=SARS-CoV-2,%20taxid:2697049&ids=MZ039713%20MZ039682%20OV522182%20
Its perfect for interrogation, it has the full spectrum of epidemiological metadata, and further details are in this post here:
How to programatically download SARS-CoV-2 fasta from NCBI/Genbank via API?
The Science journal page seemed a very good write up of GISAID BTW.

If you are wanting GISAID for networking, thats a bit tricky because NCBI is not going to do that.
You can reapply for GISAID access, but at the time I was much more impressed by NCBI. The speed at which NCBI produced an awesome database was truly jaw-dropping and they'd implemented full pangolineages classification at an early stage. Thus I just use NCBI.

Thanks for the edit @gringer, much appreciated. ".. it has all the data" is good interpretation of what I wrote, but my writing wasn't clear. Specifically, "it has all the epidemiological metadata" is what I intended to say.


text:
Parsimony-distance neighborhoods in tree space
I have the following question about the paper A parsimony-based metric for phylogenetic trees by V. Moulton and T. Wu.
In this paper, the authors define a "parsimony-distance" between phylogenetic trees on the same set of taxa. They study how this distance is related to other common distances used in phylogenetics, such as counting tree bisection and reconnection (TBR) moves.
For each notion of distance, we can define the "neighborhood" of a tree as the set of trees which have distance one from the first specified tree.
It is claimed in the paper that the parsimony-distance neighborhood contains the TBR-distance neighborhood, for any tree.
Lemma 6.2 states that only caterpillar trees have unit distance to some other tree. This seems to suggest that for any non-caterpillar tree, the parsimony-neighborhood should be empty.
However, the TBR-neighborhood of any tree should be non-empty, since it is possible to apply a TBR move and obtain a different tree.
Are there some assumptions I am missing about when the parsimony-neighborhood is empty, or when the TBR-neighborhood is contained in the parsimony neighborhood?
Also, when applying the formula in Corollary 6.6 to the non-caterpillar tree on 6 leaves, I am unsure how to apply the formula to give the neighborhood size  $|N_p| = 0$, which should be implied by Lemma 6.2.

Another term for "caterpillar tree" is "ladderized tree", the former seems more common in the math community.

The divergence in language makes it difficult to conceptualise and without that it's difficult to approach algebraic abstraction, albeit I've worked on tree space and its landscapes. I suspect I get point 3, but not much more.

"Caterpillar tree" is a term in graph theory. If that is the definition used here then every parsimony tree would be a caterpillar tree, there would never be a "non-caterpillar tree". This is because each branch length defines one "step": "parsimony's version of a mutation", so each "step" within the branch could be considered to be a "branch-less node", or else bifurcating node.

Ladderized tree I don't understand that. Maybe Np = 0 means a complete polytomy, i.e. a star tree, but could be a branch length = 0, i.e. genetically identical taxa.

The concept for an empty neighborhood ... possibly if all topologies are equally parsimonious this could degenerate into an empty tree neighborhood, again a star tree. What parsimony does is always define a bifurcating tree regardless, over replication all alternative solutions may result collectively in an "empty neighborhood". This might be fine from an algebraic perspective, its weird from algorithm perspective though.

When a TBR-neighborhood is equal to a parsimony neighborhood that actually makes sense. I define that as the TBR-neighborhood being "trapped" within the parsimony neighborhood and thats convergence. Here the "step length" has reached its minimum plateau and that tree is the most parsimonious tree ... that makes complete sense. This also makes sense within the context of the paper (I haven't read it - but I can see where they are heading).

What the authors want to do is pick holes in TBR so they want to define the point at which TBR has reached its most parsimonious tree ... I assume they want to show TBR can make errors and their solution is better. Thus they need to define the "end point". At this "end point" they want to show their measure outperforms TBR to reach a lower step size.
Thus they will want to imply that TBR should be replaced with their derivation. Its a nice methodological way to basic (or "fundamental") phylogentics.


text:
What happened to the PhyLoTA browser?
A paper I'm reading from 2021[1] mentions doing some analysis on data from the PhyLoTA database, with this reference citation:

M. J. Sanderson et. al. 2008 The PhyLoTa Browser: processing GenBank for molecular phylogenetics research, Syst. Biol. 57(3) 335-346.

This paper abstract claims to have an online interface at http://loco.biosci.arizona.edu/pb, but this link is not currently responsive. Does anyone know what happened to this project, and if it is currently accessible at some other location?

[1] S. M. Harrington et. al. 2021 Properties of Markov chain Monte Carlo performance across many empirial alignments MBE

Phylota was implemented in R and is present on GitHub here, described by Bennett et al (2018).
Here's the code...
library(phylotaR)
wd <- '[FILEPATH TO WORKING DIRECTORY]'
ncbi_dr <- '[FILEPATH TO COMPILED BLAST+ TOOLS]'
txid <- 9443  # primates ID
setup(wd = wd, txid = txid, ncbi_dr = ncbi_dr)
run(wd = wd)

Sanderson et al  was 15 years ago, albeit Sys. Biol. is a very reputable journal in phylogenetics. I've deployed solutions published there with outstanding results.


text:
Blast output file only shows 500 lines -outfmt 6
I had created databases of different sets of metagenome datasets - one with 6 runs, other with 48 runs, another with 100 runs, etc using the accession list for each of these datasets and makeblastdb command. Since I wanted to check the abundance of one specific 16S rRNA gene sequence, I ran a blastn command using this 16S sequence (F7) as the query against each of the created databases with a txt file (F7_droughtS.txt, F7_Dy.txt, so on) as the output. Ideally I should have got differently sized files with different number of OTU matches. I did get different files since the SRR ids were different for each database, but the results did not go beyond 500 lines for any database search. Can anyone explain why this is so?
Following is the command I used for blastn of:
sequence.fasta against database firstdb :
blastn -query sequence.fasta -db firstdb -outfmt 6 -out output1.txt

Also ran the following to know the default number of lines :
blastn -help -max_target_seqs
Number of database sequences to show one-line descriptions for

Not applicable for outfmt > 4 Default = 500'     * Incompatible with:  max_target_seqs  -num_alignments <Integer, >=0>    Number of database sequences to show alignments for    Default = 250' * Incompatible with: max_target_seqs

I suggest you have a look at the settings for the blast version you are using. These can usually be accessed with 'blastn -help'. In particular, '-max_target_seqs' may be of interest.

Please [edit] your question and i) remove the images of the files, that isn't telling us anything. You can show us the output of wc -l  file to show the number of lines. ii) show us the exact commands you used. You are probably using a version of blast that shows the top 500 hits by default, but we need to know what blast you use and what exact commands you ran in order to help.

@Laura, Thank you for your suggestions. When I did run the code provided by you, this is what I got :
Number of database sequences to show one-line descriptions for
Not applicable for outfmt > 4
Default = 500'     * Incompatible with:  max_target_seqs  -num_alignments <Integer, >=0>    Number of database sequences to show alignments for    Default = 250'

Incompatible with:  max_target_seqs

How can I change these settings?

Not an expert on 16S, but is BLAST really the correct tool for this(?) Maybe one of these tool suites like QIME

QIIME is an interesting suggestion. Here I'd say BLAST - although there are performance issues - offset parallelising. BLAST is a known approach to scoop all the genetic diversity for a given locus and 16S is a taxonomic goldstandard. This isn't the analysis - the analytics are in the next step, but how you've prepped the data for the next step needs to be understood and be solid. There are a lot of other approaches but blast has flexibility e.g. to shift tblastn .. etc ... and its why I like it.

Personally, IMO "modulisation" custom data pipelines in general, i.e. loads of steps - gives more control and each developer specialises on their component rather than try and do everything at once. The latter is the QIIME approach, but its an interesting debate.

I've reproduced the issue. The solution is to simply set -max_target_seqs to whatever number of hits you want. E.g. -max_target_seqs 1000 gives 1000 and so on ... The problem  is specific to the -outfmt option and not default searches.
blastn version  I'm using v2.12.0, Feb  8 2022

Proof
My database contains around 260000 fasta sequences.
makeblastdb -in mystuff.aln.0.fa   -parse_seqids -blastdb_version 5  -out test.db -title "test" -dbtype nucl
blastn -query query.fas -db test.db -evalue 1e-6 -num_threads 4 -out blasthits.txt

Then wc blasthits.txt

>5000 hits

If we now use -outfmt 6, but remove the stringency of -evalue 1e-6 then ...
blastn -query query.fa -db test.db -outfmt 6 -num_threads 4 -out blasthits.txt

Now perform wc blasthits.txt

500

So the OP is correct, I didn't know -outfmt 6 restricts the output to 500.
If we now introduce the max_target_seqs 1000 parameter, suggested by @Laura then,
blastn -query query.fa -db test.db -outfmt 6 -max_target_seqs 1000 -num_threads 4 -out blasthits.txt

Then wc blasthits.txt

1000

If its 550 wc ... is 550 etc...
Thats it. Note in my blasts I usually parallelise across 4 processes i.e. -num_threads 4.

Thank you very much, it worked!

Thanks @KirtiKulanthaivel glad it worked and much appreciated.


text:
Gene rank scores from LINCS database: Evangelista et al. 2022
The website https://maayanlab.cloud/sigcom-lincs/#/SignatureSearch/UpDown can be used to "Identify reversers and mimickers from over 1 million signatures by entering up and down gene sets .. " I am looking to get these "reversers" and "mimickers", aka chemical perturbagens.
So on this website, I paste my list of a selection of Up-regulated genes, in their respective box, and the same for my list of down-regulated genes. These genes are from an RNA sequencing experiment.
After the website has done its thing, I expand the search for chemical perturbagens (LINCS L1000 Chemical Perturbations (2021)). One of the data tables I can download contain something called a Gene Rank Matrix. This table consist of one row for each gene and columns of the chemical perturbagens found. The values are integers.
A small subsection of the data table would look like this:

My questions are what these integers (12, 144, 411, 641) represent. Are they how much the perturbagens (In this case REP.B025_YAPC_24H_K01_resiquimod_2.22uM) will affect the genes?
The paper describing this website searches (https://academic.oup.com/nar/article/50/W1/W697/6582159) mentions "The expression tables provide the computed CD differential gene expression coefficients for each signature, while the rank matrices provide the integer gene ranks for each signature as determined by the coefficients"
and
"Matrices that represent the gene rankings of the L1000 signatures are stored as HashMaps by the data API in random access memory (RAM). HashMaps are data structures designed to efficiently search for key-value pairs, and this implementation drastically cuts the fetch times of the gene rankings by UUIDs, thus improving the speed of the search"
Does this mean that the scores found wont have any meaning in regards to how well the perturbagens found will affect the genes? Or can some meaning be derived from this?
Thanks in advance!

Summary answer
Although I'm unhappy about how the paper is written and possibly the methodology, the following sentence is likely (almost certainly) to be the answer.
The coefficient (aka weight) is how important that gene is for the classification of that signature. Thats it! :-)
Essentially, there are loads of genes involved in identifying which genes represents which signature. Not all those genes are equally informative. The positive or negative skews in that calculation, called a "classification", are identified through their weights (aka coefficients). Thats simply the answer. The bigger the coefficient the more important that gene is in the classification, i.e. its associated with the signature.
There is no need to be more complicated than that. If you are interested in a critique of that paper, which I had to do to get to the answer and thereby understand the analytics of that paper, this is presented below.

Critique of Evangelista et al. 2022: Nucleic Acids 50:W697–W709
Note 1 I don't have time to read the paper in absolute detail, particularly the appendices, where some of the information needed may have be placed.
Note 2 In the answer below, a "expression coefficient" is what I've called an "expression weight" its the same thing.
Weightings and analytical concerns I do understand weightings in deep learning, particularly machine learning. I read the paper and the ANN (artificial neural network) was very basic. I don't believe it had multiple layers, because it would have said so and it was a very odd description for a deep learning algorithm. It would also be odd to have a single layer in deep learning, very odd. The loss function is also not described at all which is extremely important in deep learning (it will not work otherwise). Essentially, what the authors call "deep learning" appears to a machine learning neural network using an Adams gradient descent.  Adams gradient descent is a standard approach in deep-learning admittedly. Thereafter in the paper the details were very sketchy, maybe they are in the appendix (should be checked). Any deep-learning, machine-learning requires a demonstration of the robustness of the training and this appears absent. The two basic stats are accuracy and AUC-ROC, which were definitely absent. Universal stats are used to quickly understand how good the algorithm is in comparison to other deep-,machine- learning algorithms for other data sets.
Analytical overview Generally I would assume the authors have performed a supervised learning algorithm, but this may not be correct (below). If it is supervised learning they have performed they would have identified their 'signatures' against a very large data set, then trained the different types of signature against that data and tested it. The learning could be terrible - but it's just not stated at all. I'm surprised because Nucleic Acids is at the pinnacle of algorithm publications.
Summary. The authors performed deep-learning (which looked like a neural network machine learning algorithm). The description is at best poor, but further information could be present in the appendix. Therefore its not absolutely clear what a "weight" would represent because there are two possibilities.
Weights possibility: unsupervised learning The one thing that might have happened is they've performed an unsupervised learning algorithm and the weights are a sort of "maximum variance". This very trendy. That would explain the lack of accuracy and AUC-ROC type stuff. IF this is true and I don't believe it's explained - I could be wrong - then the weights represent maximum information. What does that mean? Good question, its context specific, its the gene which is most important in sectioning up the data - e.g. housekeeping gene which is required 'all over the place', not just part of the cell cycle is an example. Basically, in unsupervised weighted output high ranking weights are a mix of really boring and really brilliant stuff. So you just go through it with biological understanding, is this gene boring yeah/neigh, is this one - wait wow! unexpected ... That is actually fashionable analytics :-)
Weights possibility: supervised learning Is already answered at the beginning. This is almost certainly the correct answer - it would be the orthodox approach - and we just accept the robustness checks for the algorithm and possibly the algorithm itself were simply insufficient.
Why isn't the paper clear (at best)? What I assume has happened is the team have brilliant frontend and genomics team backend. The analyst in the middle is unlikely to be a specialist and thats why it's poorly described, or possibly performed. I assume if the reviewers were of similar composition then that's how its got through Nucleic Acids. Thus, I would personally treat the output with caution because deep learning (machine-learning) will be the essential bottleneck between their backend development and your data. Maybe "bottleneck" is not quite an accurate description, because its speeding up the classification, so maybe a "gate-keeper" is a better description.


text:
open and read Fasta file (raw data)
I have a big fasta.dataset file containing half a million proteins (1.0 GB). I have four lines for each protein code:

line 1:the protein code
line 2: protein length in amino acids
line 3: amino acid sequence
line 4: secondary structure

Now, I am trying to open and read it in python (Biopython), and it does not work:
filename = 'pdb.fasta_qual.dataset'
sequences = SeqIO.parse ( filename,'fasta')
for record in sequences:
    example = record
    break
print(example)

How can I read it in python and loop through the file to look at line 3  for each protein to count the sequence length and distribution?
here is the first 5 lines of my file: which my file contains 500,000 proteins for each one has a 4 lines (name ,len of protein in amino acid,the seq represents by letters which what I would to calculate,the secondary structure)
4LGTD
247
M       S       E       K       L       Q       K       V       L       A       R       A       G       H       G      T
.       .       E       E       H       H       H       H       H       H       H       T       T       S       S       .

I want to open and read the file and loop through line 3 for each protein to calculate the length of the sequences and plot a histogram ! to check the distribution.
,,,
The output i am expecting is :
The len for the first seq is  =
The len for the second seq is  =
Until
the len of the last sequence  which is number (500.000)=
And then i can plot a histogram for the len of the sequences
,,,
NOTE: I have opened and read the file's info by Linux, but I could not by python.

Hi and welcome to the site! We need more detail to be able to help you, so please [edit] your question and i) add a few lines of your fasta.dataset file; maybe 3-4 protein entries so we can see what we are working with since what you describe is not fasta; ii) show us the exact output you would want to see from that example input; iii) show us the code you wrote so we don't need to start from 0; iv) explain exactly how it "did not work"; v) tell us what operating system you are working with. Note that since this is not fasta, I doubt there is any reason to use biopython.

Oh, and please use the formatting tools to format your input and output examples and the code you show as code (the {} button).

Thanks for the edit, but please don't post images of text: we can't use those to test our solutions. We need a few lines of your actual file and the output you want to see from those lines. Both as text so we can copy/paste them.

And what output do you want from this? Just the number 247? I ask you again, for the third time, please show us the output you are expecting. And does it have to be in python? I mean, all you want here is awk 'NR%4==2' file, is that enough or do you absolutely need a whole python script?

awk 'NR%4==3' file is calculated just the first sequence while I have half a million of sequences in the same dataset

no, that will print every line in the file whose line number modulo 4 is equal to 3, so it will print the 2nd, 6th, 10th, 14th and so on. But are you now saying you want some string along with the length? Should the solutions print The length for the Nth seq is and then the length?

Sorry before we continue we need to establish why it's not working. Is it the file size is too big and it just freezes or is it throwing a bug? There are bugs in the code, but passing this via Biopython will likely hit a RAM bottleneck

This awk 'NR%5==3{ print length ($0);}' file name  works and gives me the len of characters in the sequences line 3. Now how can I plot a histogram for it? Also calculate the frequency of the same line?

Based on your description, your file is not a FASTA file. I suggest you modify your question. Here is the description of FASTA file> https://www.ncbi.nlm.nih.gov/genbank/fastaformat/. Perhaps that is why it's not working.

Just read every third line in the file and do sth with it.

I better try and explain the issues a bit better. A dataset of 0.5 million at 1 GB will struggle to load into Biopython. Biopython isn't great at handling large numbers of sequences - for a regular desktop/laptop - and what will happen is the parser will freeze or else take an enormous amount of time (many hours). I call this a "RAM bottleneck". A generator is an easy workaround.
The other issue with the code is there are lots of bugs combining format and code. There is a clarity issue as well it looks like you want to calculate the secondary structure - rather than import the secondary structure alongside the alignment. I could be wrong on this though, certainly secondary structure is present in the ad hoc alignment. There is an alignment format that handles both secondary structure and sequence data. You might think about restructuring the question and your approach into multiple questions.

Reformatting the ad hoc alignment to permit import into Biopython on a pilot data set (not 0.5 million sequences). This is close to @terdon's answer.
Code that will process very large files in Biopython. This is the answer below.
An alignment format via Biopython (probably output - could be input) that integrates linear sequence with secondary structure.
Plotting a histogram for 0.5 million sequences.

Points (Questions) 1. to 3. can be leveraged by Biopython. Point (Question) 4 can be done via Python too - graphing/histograms are pretty good. There seems an underlying issue of project management for coding based work.
If this is a RAM bottleneck then solution is below.

The easiest approach is to use a generator (please see 'format').
example = (record for record in SeqIO.parse(filename,'fasta'))
print (next(example)) # example.id ?

A generator will deal with the sequences one by one instead of trying to load the whole thing to memory. It may not be suitable for certain complex calculations, but it is by far the easiest approach to the problem (if the problem is the OP is struggling to load a million sequences into memory).
Note I always used yield commands but @Steve introduced the () format and I definitely prefer it. Also note that manipulating the output of a generator can be different to a list, tuple or dictionary (calculation depending).

Format
The format supplied isn't fasta so the SeqIO will not work in fasta mode and should throw an error/exception. The format you need to leverage Biopython is very different. I would recommend a standard 'open' command and passing this to a custom dictionary. Also the output of print(next(record)) should produce an object symbol, not the actual data e.g. record.id would be needed to see data: it does however demonstrate the generator is working.

It sounds like all you really want is to print every 2nd line of a text file. If so, you don't need Python, let alone BioPython, you can do it with basic *nix tools:
$ awk 'NR%4==2' pdb.fasta_qual.dataset
247

NR is the current line number, and % is the modulo operator. Therefore, NR % 4 will equal 2 on the second line of each group of 4 lines in the file. In awk, the default action when something evaluates to true is to print the current line, so this will print out all length lines from your input file (assuming you don't have blank lines anywhere, which we don't know since you only showed the first entry).
If you want to add some sort of message to be printed, you could do:
$ awk 'NR%4==2{ print "The length of protein",++c,"is:",$0}' pdb.fasta_qual.dataset
The length of protein 1 is: 247

If you really must do this in Python, you could do something like this:
#!/usr/bin/env python3

import sys

#my_file= sys.argv[1]
line_number = 0
with open(sys.argv[1]) as f:
    for line in f:
        line_number += 1
        if (line_number % 4 == 2):
            print("The length of protein %d is: %d" % (line_number,int(line.strip())))

And then you would run it like this (assuming you named the script foo.py):
$ python3 foo.py pdb.fasta_qual.dataset
The length of protein 2 is: 247

Out of curiosity, I created a test file by repeating your example 10000000 times, which left me with a 2.4G file:
$ perl -e 'open(A,"pdb.fasta_qual.dataset"); @lines=<A>; print "@lines" x 10000000' > testFile
$ ls -lh testFile
-rw-r--r-- 1 terdon terdon 2.4G Feb  7 17:06 testFile

I then tested both approaches for speed. Unsurprisingly, awk was  much faster:
$ time awk 'NR%4==2' testFile > /dev/null

real    0m6.447s
user    0m5.917s
sys     0m0.492s

$ time foo.py testFile > /dev/null

real    0m10.856s
user    0m10.111s
sys     0m0.523s


text:
Why might parabricks not be writing output to the output file I provided?
This feels like a silly question, but I am testing parabricks bammetrics and nothing is being written to my output file. This is the command I am using:
docker run \
  --rm \
  --volume ~/workdir:/workdir \
  --volume ~/ouputdir:/outputdir \
  -w /workdir nvcr.io/nvidia/clara/clara-parabricks:4.0.1-1 \
  pbrun \
  bammetrics \
  --ref /workdir/reference.fasta \
  --bam /workdir/cram.cram \
  --interval-file /workdir/regions.bed \
  --out-metrics-file /outputdir/metrics.txt

But metrics.txt doesn't have anything written in it when I open it after the command is run. Am I doing something wrong? I am very new to parabricks, so I assume I am missing something simple
I apologize if this is the wrong place to ask this question.

EDIT Summary
For future users who encounter the problem the solution is,
docker run --platform linux/amd64 --rm -it -u 0 --entrypoint /bin/sh  nvcr.io/nvidia/clara/clara-parabricks:4.0.1-1

I am unhappy with the practice of limiting Docker access via blocking standard runtime: there are statements of intent on SE as a strategy to attempt to restrict Docker access. Parabricks isn't an isolated incident of deploying this setup within bioinformatics. However, I checked their forum and there is no statement of intent that it's deliberate. Parabrick's copyrights are here. There would be issues if there was ever an official statement to this extent. They would likely argue it was for other reasons to automate stuff.

Original answer ...
I don't know parabricks, but I know Docker. What you are attempting is to run Docker from the outside ... I would get on the inside to view the error log. Whats happening is the error is not being passed outside the container. Getting on the inside requires
docker run  --rm -it --name mydocker nvcr.io/nvidia/clara/clara-parabricks:4.0.1-1

This should give you access to the container. If it doesn't post back and I'll forward a different command that will. It's impossible to block access to a Docker container albeit some people try.
For simplicity is worth trying
docker run  --rm -it --name mydocker clara-parabricks:4.0.1-1

Anyway, if you are using this approach you'll need,
docker cp *outside container* *insider container* nvcr.io/nvidia/clara/clara-parabricks:4.0.1-1 and vice versa

... this is an amazing command and allows you to move data in and out of the container. You've got to get the data on the inside of the docker and docker cp does that very well. I would then look to running the whole thing from the inside of the container. Thats exactly how I run Dockers, then I exit and return when it's done. There might be a bit of clean via this approach, I'm not sure --rm might take care of it.
You'll also be able to look at the programs and copy it out onto hard drive. This last bit may be not advisable - the whole point of Docker is to protect your system from their algorithms ... but you can definitely do it and if you get it paired up with NextFlow thats cool.
I don't particularly want to pull the Docker container because they're often huge and I've no interest in the algorithm albeit I've spotted its at nvcr.io/nvidia/clara/
Docker, to be honest it's lazy by the developers. Docker is - techie packaging. So it's good for reviewers and behind web firewalls. It's useful but like you've just experienced  if it flops you're stuck. The one positive thing about Docker (apart from DevOps) is it's generally tolerated to package someone else's algorithm as part of a pipeline and let a third party 'pull it'. This is because they'll be a Dockerfile which can reconstruct the whole thing without copyright infringement. So if that first developer objects ... the Docker developer would just issue the Dockerfile. conda and stuff max out on that concept.

There's a weird error as the OP described below.
Right Docker uses by default a tag latest that does not work for parabricks. The tag must be 4.0.1-1. Could you retry the above code please?

Update For good measure try the following:
docker run  --rm -it --name mystuff  nvcr.io/nvidia/clara/clara-parabricks:4.0.1-1 --entrypoint /bin/bash

Might not be bash, thus
docker run  --rm -it --name mystuff  nvcr.io/nvidia/clara/clara-parabricks:4.0.1-1 --entrypoint /bin/sh

They shouldn't throw an error, but you might still find you are outside the container after running it. Lets see.

Update
The error looks okay, because a runtime error is exactly what would be expected if they are attempting to prevent access. Therefore please try:
docker run --platform linux/amd64 --rm -it -u 0 --entrypoint /bin/sh  nvcr.io/nvidia/clara/clara-parabricks:4.0.1-1

This, or a variant of it, should work.

Thank you for the response! When I ran that command it didn't work, I get this error:
edited
I've read that NVIDIA does not want to give access to the parabricks container

THANK YOU! I was able to get inside the container and run the command. For some reason, running it from within the container caused the output to be written to the metrics file... I'm not sure why it would be different, but I'm glad that it worked nonetheless. Thank you so much again for your help!

Just want to clarify that I read that nvidia doesn't want people to access the container somewhere on the dev forums, not in the documentation (and I was paraphrasing, the original dev response was "the best way to run parabricks is outside the container").  edited

Okay thanks @womy. Checked the forum, its fine.  Note some of the comments in this thread have been deleted or edited to comply with extended discussion SE stuff


text:
Phylogeny building in R from FASTA files:
Im working on building a phylogeny from scratch with downloaded FASTA sequences from GeneBank. I think Im doing alright up until the multi sequence alignment in the msa package where I get an error Error in convertAlnRows(result$msa, type) : There is an invalid aln file! As always thanks in advance for taking a look: Here is the reproducible example:
library (ape)
library (rentrez)
library (msa)

#Download Species Data
# B_terrestris
B_terrestris <- entrez_fetch(db = "nucleotide",
                                  id = "NC_045179.1",
                                  rettype = "fasta")

# B_hypocrita
B_hypocrita <- entrez_fetch(db = "nucleotide",
                             id = "NC_011923.1",
                             rettype = "fasta")

#Vespa
Vespa <- entrez_fetch(db = "nucleotide",
                      id = " MT137096.1",
                      rettype = "fasta")

seq<- c(B_terrestris, B_hypocrita, Vespa) # gotta clean it

#FASTA Cleaner Function
fasta_cleaner <- function(fasta_object, parse = TRUE)
{
  fasta_object <- sub("^(>)(.*?)(\\n)(.*)(\\n\\n)","\\4", fasta_object)
  fasta_object <- gsub("\n", "", fasta_object)

  if(parse == TRUE){
    fasta_object <- stringr::str_split(fasta_object,
                                       pattern = "",
                                       simplify = FALSE)
  }

  return(fasta_object[[1]])
}

#This should be ready to go!
clean.seq<-fasta_cleaner (seq, parse = FALSE)

str (clean.seq)

msa (clean.seq,  type="dna", method="Muscle", cluster="neighborjoining")

@M__ Thank you! What sequence alignment method would you reccomend instead?

The formal question I've been given is what aligner would I use? There is an issue between forming a data pipeline for ape and the latest and greatest. Its trade and the compromise would be msaClustalOmega(), but the rationale is complicated.
I strictly use muscle5 or specifically muscle -super5 option as a standalone. I am certain R msa has not introduced muscle5 and I don't even know if it's published yet. It's on bioRxiv.org. The super5 option must be specified and once done its amongst the most incredible algorithms per se notably in speed - no other aligner comes close - and personal benchmarks in accuracy. To give some idea, I was going to write an aligner as part of a package (its makes stuff simpler) and thought - whats the point with muscle5?
Its competitor was MAFFT and published benchmarks demonstrated it was the strongest, before muscle5 was released.
The dilemma is simply that ape is cool and if you are wanting to have a seamless pipeline then greatest aligners don't fit with that. MAFFT isn't in msa and seriously doubt muscle5 will be. HOWEVER Clustal Omega is in msa and that is the strongest (most recent) of the options it has. Omega does have published benchmarks and its neither quite as fast nor quite as accurate as MAFFT.
If you've got 50000 sequences - it's muscle5 no contest. If you've 50 sequences then keep the pipeline together via msaClustalOmega.

The original question and code appeared to involve converting fasta inside R onto msa. Thats how subprocessing works outside R (or Python) but it's not the OO approach.
My original answer was in comments and remains an unofficial answer because usually I run the code I'm recommending and I don't fancy subprocessing muscle via R.
The error was the msa algorithm was attempting to shift the object into a format that will subprocess via muscle and it doesn't understand fasta format (it wants the raw object). Thus, you appear to be wanting to pass "official fasta" into msa. It doesn't need it, it simply needs a sequence object. It's imported as fasta but then it should be a sequence object (it's no longer fasta). That should go straight into msa. Thus seq should load straight into msa with "fasta" option removed. Finally could you remove msa and used msaMuscle then simply msaMuscle(seq).

I tied this bit with the set of sequences  msaClustalW(c(B_terrestris, Vespa, B_hypocrita), type="dna", cluster="nj") but now I get this error: "Error in .Call2("new_XStringSet_from_CHARACTER", class(x0), elementType(x0),  :
key 69 (char 'E') not in lookup table"- Any further suggestions?

At least it can see the object. If this was Python I could debug it (well I could write the code to do it straight up), R will be more challenging. You need to print the seq object and check its integrity.  It seems to be complaining that E isn't DNA, suggesting either this is protein (unlikely) or an E from the header has got caught in the sequence. There are alternative ways to construct the sequence object, so as a last resort could do that. The idea here is to split the object attributes into sequence and ids then reassemble it.


text:
Why is COI metabarcoding not used for prokaryotes?
I'm a beginner in bioinformatics working with environmental metabarcoding, and while reading up on different marker genes, I understand that COI is preferably used for animals thanks to its higher variability over 18S (although it depends on which taxonomic resolution fits best and which reference data is available). What I don't understand though, and can't find an answer to, is why COI metabarcoding isn't used for bacteria and archaea as well. Why is 16S preferred, although bacteria and archaea clearly have the COI gene? Is it not variable enough? Or are the primer binding sites too variable? Something else?

Hi @Joel welcome to the site. This is a reasonable question and I'm happy to answer it, however could you state the taxonomic groups in your study?

The exact answer to your question is bacteria do not have a mitochondrion, therefore cox1 is not a universal gene combining both Archaebacteria and bacteria. In contrast 16S is universal for these two kingdoms AND it is also single copy. Moreover, bacteria have diverse tropisms beyond the Krebs or Calvin cycles. Extensive studies have been performed on 16S and it provides an excellent taxonomic target in bacteria.
18S is a problematic target in eukaryotes because it can paralogous and that is not an ideal target in taxonomic studies, but has nicely conserved primer site: thus single copy genes are required/preferred. In addition, 18S is prone to indels (insertions/deletions) and these are difficult to model, align and understand the underlying molecular mechanism or selection advantage.
Every animal on earth has a mitochondrion and therefore cox1. Mitochondrion genes are guaranteed to be single copy in higher eukaryotes and cox1 has been extensively sequenced throughout metazoans. Thereby barcoding is appropriate because there's an exact reference to identify a given species is usually available. It has been extremely successful in insects which are the most biodiverse and species rich taxonomic Class on earth: there are more insects than all other higher eukaryotic species combined.

In answer to the comments ...

Whether mitochondrion were originally bacteria is a very difficult question and remains a theory (compartmentalisation theory - very old). I have not seen any rigorous phylogenetic support for this.

cox1 (COI) or cox2 (COI) is not universal to bacteria, e.g. here. Its primary use is for aerobic organisms (oxygen consuming). Thus cox1 has a strong Krebs cycle association but there are a multitude of anaerobic pathways.

Anaerobic bacteria don't necessarily need cox1/cox2 and there are a lot of obligate anaerobic bacteria in addition to non-obligatory anaerobic bacteria (which can do both).
Marine metazoan are obligate aerobic organisms.

There is no reason to suppose bacteria  are strictly  single copy cox1. mtDNA is tiny and hence copy number heavily constrained, bacterial genomes by contrast are much larger and could easily accommodate duplication events.

There isn't much work on bacterial cox1, it has been pursued to attempt to augment antibiotic therapy. Degli Esposti (2020) is the closest I know to a comprehensive discussion on the bacterial cox genes.

Thanks for your answer. A few questions:
• I know prokaryotes don't have mitochondria, but that shouldn't make a difference since mitochondria are descendants of bacteria, right? As far as I understand, both eukaryotes, bacteria and archaea have the COI gene.
• What do you mean by "diverse tropisms beyond the Krebs or Calvin cycles" and how does it matter?
• I had not considered paralogy, but now that you bring it up I understand that it must be problematic to have different variants of the same marker gene in the same individual. But do paralogues of COI really exist in bacteria/archaea?

I'm primarily interested in marine metazoan diversity, so COI is what I'm being told to use. But I just fancied what a lovely story it would be if those universal COI primers could be designed to work on prokaryotes as well, and then it got me trying to find out why not. If you know any paper that discusses this, I'd be happy for it!

@Joel addressed in comments.

Thank you @M__. You bring back parts of my biology education that I'd forgotten, and now it all makes sense. If I'll do prokaryotic metabarcoding, I can now rest assured that 16S is the correct choice. But I still hope for someone to find a universal marker gene that does work for eukaryotes and prokaryotes alike :)


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:
Generate VCF from different .bam files with different chromosome names
I have two resources of .bam files. One is generated by our lab (1 sample = 1 bam). One is downloaded online (again 1 sample = 1 bam). For the downloaded samples the chromosomes are labelled: chr1, chr2, chr3 etc For our lab samples, the chromosomes are labelled: 1, 2, 3 etc. I want to generate a single VCF file of variants across all samples.
I'm using bcftools:
bcftools mpileup -Ov -f ref.fasta -b samples.txt | bcftools call -mv -o bamMge.vcf

However, I get no calls and the repeated error: [E::faidx_adjust_position] The sequence "1" was not found
I have two questions:

Is my strategy correct (i.e. bcftools mpileup)? Do I need to incorporate extra steps (NB I also have matching gVCF files)
How can I either (i) alter the chromosome labeling of one of the subsets of bam files, or, (ii) use a mapping file to match chromosome labels during the mpileup run?

Thanks

The main issue here is that your bam files have different chromosome labels. (i.e. 1,2,3 vs. chr1,chr2,chr3) as you mentioned.
This suggests that the data was mapped to a slightly different version of the reference genome. As a first layer, you might want to double check that these reference chromosomes are cross compatible (e.g. both correspond to hg19 for instance). To find more information on the reference genome used, you can check the header information using "samtools view -H "
Assuming that the reference genomes are cross compatible, you can try converting it to sam, adding the missing "chr" label in a pipe, and converting it back to a bam file. However, this will also require modification of the header file which can be challenging in itself.
The easiest (and more rigorous) approach is probably just to convert the batch of bam files with the wrong chromosome labels back into fastq files. You will then want to remap them using the exact same reference genome as the other batch of bam files. After which, you should be able to perform variant calling on the combined batch of bam files.

To add to this very correct answer, here are a couple of documentation articles that explain the key points in greater detail: https://gatk.broadinstitute.org/hc/en-us/articles/360035891071-Reference-genome and
https://gatk.broadinstitute.org/hc/en-us/articles/360035531652-FASTA-Reference-genome-format

I am not familiar with getting variants from bam files, so I will just answer the second part. From the error message, it seems that it is necessary to have consistent chromosome labels across datasets. It is also a good idea to do so because this keeps your datasets consistent and makes things a lot easier. To change the chromosomes labels, you can write your own scripts (I use Python, for example). You can also find several solutions on BioStars (https://www.biostars.org/p/98582/).

Editing BAMs or VCFs to change chromosome names is really dangerous if you're not 100% sure you know what you're doing. I've seen scripts that screwed up the file contents in subtle enough ways to mess up downstream analyses. Even if you avoid that, you can end up combining data aligned to different genome reference builds, which also ends up affecting analysis results. It's much safer to realign from scratch like KT8 describes.


text:
Possibility to save output blastn table in memory using biopython
Is there a possibility, using biopython to save output table of blastn in memory, not in file on hard drive, for process it using pandas for example and then delete it from memory?
I mean, is there a possibility to redirect output of "out" argument of this function into memory? In pandas DataFrame for example.
NcbiblastnCommandline(query="fasta1.fasta", subject="fasta2.fasta", out="output_table.csv", outfmt=6)

I always saved output tables in file before, but now I need to align one set of small fasta files against another set of small fasta files and I think that the best way is avoid writing of many small tables in hard drive because after recieving one table I have to read and process this table immediately and after that I don't need this table anymore.

You can specify send to stdout using out='-' in the Biopython wrapper.
from Bio.Blast.Applications import NcbiblastnCommandline
import pandas as pd

cline = NcbiblastnCommandline(query='seq.fna', subject='seq2.fna', outfmt=6, out='-')
output = cline()[0].strip()

rows = [line.split() for line in output.splitlines()]
cols = ['qseqid', 'sseqid', 'pident', 'length',
        'mismatch', 'gapopen', 'qstart', 'qend',
        'sstart', 'send', 'evalue', 'bitscore']

data_types = {'pident': float, 'length': int, 'mismatch': int,
              'gapopen': int, 'qstart': int, 'qend': int,
              'sstart': int, 'send': int, 'evalue': float,
              'bitscore': float}

df = pd.DataFrame(rows, columns=cols).astype(data_types)


text:
Analyzing a blast result of NGS data
I am analyzing BLAST results from total RNA NGS, pair ends. This BLAST results was performed against the ssRNA virus database download from NCBI.
I filtered those read data whose paired ends uniquely matched a single virus and I have both (removing single reads).
Moreover, based on the sstart and send values from a unique virus in each library I calculated the difference (ssdif) of the lowest (sstart) of the highest (send) value that I got in the BLAST results,
for example: If I got 4 reads that match against one virus in one library (2 reads with pair ends) I got the lowest sstart and the highest send for these four hits.
Then, I calculated the coverage using the ssdif; coverage=(ssdif/genome length for that virus) * 100
After I calculate the coverage for each virus in each library based of the BLAST results I calculated the mean coverage for each genus and for each family of all viruses detected in my analysis.
Questions

Does it make sense the way that I calculated the coverage? I would like to filter my data for those viruses that have similar or higher coverage to the genus or family (previously calculated)

Based on what taxonomic level does it make more sense filter them?

Thanks in advance