text:
How to tell if our ligand-protein docking is good from AutoDock Vina's result
I have perform a ligand-protein docking using Autodock Vina.
The result of the docking looks like this:
WARNING: The search space volume > 27000 Angstrom^3 (See FAQ)
Detected 8 CPUs
Reading input ... done.
Setting up the scoring function ... done.
Analyzing the binding site ... done.
Using random seed: -1553787135
Performing search ... done.
Refining results ... done.

mode |   affinity | dist from best mode
     | (kcal/mol) | rmsd l.b.| rmsd u.b.
-----+------------+----------+----------
   1         -5.9      0.000      0.000.  Pose 1
   2         -5.7     22.945     25.492.  Pose 2
   3         -5.5      1.426      2.046.  Pose 3
   4         -5.5     23.669     25.616
   5         -5.4     25.783     29.152.  .....
   6         -5.3     21.146     23.357
   7         -5.2     20.323     22.545
   8         -5.2     23.864     26.064
   9         -5.1     23.422     26.585.  Pose 9

As far as I understand from these statistics Mode 1(Pose 1) is the best.
However when I actually visualize them in Pymol, Pose 1 has no hydrogen bonding at all
but Pose 2 has.
My question is how can we judge if which of those two Pose is the best to use?
Note in figure below Pose 2 has dashed line (Hydrogen bond).

It is best to contextualise the numbers. -1 kcal/mol is about the potential energy gained from a hydrogen bond —technically described in the r^6 part of the Lenard–Jones term, it is also the average collision energy of a water molecule at 37°C as that is RT ($\frac{k_b\cdot T}{N_A}$, wiki)under a Maxwell–Boltzmann distribution. A salt bridge –2 kcal/mol (Columbic force term).
So your scores are not very low, hence why you are counting two hydrogen bonds. Although you can also see a lovely sulfur-pi interaction which is –1 to –2 kcal/mol, so those bonds alone are probably making a –4 kcal/mol contribution, so I am guessing that some terms may be horrendous, such as repulsion forces etc.
A nice metric is doing a conversion to ligand efficiency, which weeds out affinities driven by size of the molecule... in this case most of the molecule is doing nothing.
Also, most programs have an accuracy of 1 kcal/mol or higher.
So one cannot say what ∆∆G is the best with the data at hand due to noise, but one can say that the ∆∆G is not low enough... Sorry.

Thanks. By  ∆∆G Do you mean  affinity (kcal/mol)?
Then what do you consider ideal score range for  ∆∆G. Since you said that my score above is not low enough?

I did and it generally is fine, but variations of definitions and calculation come into play. Gibbs free energy is ∆G (so one delta to start with), although some tools refer to ∆G_bind as ∆∆G, and generally it's intended as ∆∆G = ∆G_holo - (∆G_apo + ∆G_aq), where ∆G_aq is the free energy in aqueous solvent, but this is not always true as there may be fancy stuff correcting for ∆G_aq (which makes little sense for a high partition coefficient compound), Vina is one of them.

Higher than 10 kcal/mol... If you calculate the ligand efficiency it will be something between –0.1 and –0.2 kcal/mol/heavy. If your ligand is meant to dock with a high affinity, then I would revisit the template —a holo enzyme structure with the substrate removed is best, check the protonation of that aspartate etc.

Ah. I just noticed your box is giving a warning of being over 30 Å x 30 Å x 30 Å. Which is huge.

when you said 10 kcal/mol, how did you come up with that number?
Any reference/publication for that?

It's a ballpark figure. Playing devil's advocate, if you take the example of aspirin (an old small drug that's totally non-specific, that is taken at a 1,000 mg dose, that is charged, that is quickly broken by a liver cytochrome etc), it's binding to it's target is only -7 kcal/mol, which for 13 heavy atoms makes it's LE -0.5 kcal/mol, meaning that threshold doesn't always hold true (with massive caveats)


text:
Programmatically adding hydrogen and remove water to multiple PDB files
I'm trying to dock a ligand to several hundred PDB files (receptors).
I thus need to prepare the those PDB files like removing water and adding hydrogen.
I can do that manually using Autodock or Pymol. But it's too troublesome.
Is there a way to do it programmatically with Python or R or any command line procedure?

If you did it manually in PyMOL, it can also be scripted. But note that adding hydrogens is not that simple. There is more than one way of deciding which hydrogens are present (you may consider pH) and what are their coordinates. So I'd first think what method of adding hydrogens to use.

@marcin In pymol I use h_add command. That is not enough? How would you go about scripting it?

https://chemistry.stackexchange.com/questions/129157/adding-hydrogens-to-multiple-pdb-files-using-openbabel-gui-software here they use openbabel, not scripting but in batch

I can't tell if that's enough for what you're doing. h_add in Python script would become pymol.cmd.h_add('all')

Michael Word's reduce is an established method to do that

GROMACS can add hydrogens. And  water you can filter out with cat filename | grep -v ...

I believe PDB2PQR CLI will do the work wonderfully. Don't let the name trick you: PQR files are organized like PDB ones. Under the hood it runs propka, which is state-of-the-art for predicting a protein residues protonation state. The best part is that you can use PDB2PQR web server and they will give you the corresponding CLI arguments for each option you choose. It also automatically removes water (and other non-protein molecules).
You can also try looking into tutorials on how to setup proteins for dynamics on GROMACS, AMBER or OpenMM and stumble on their own tools for protonating stuff (respectively pdb2gmx (CLI), protonate (CLI), and modeller.addHydrogen (python API)), but as far as I know propka computes residues pKi in a knowledge-based manner and out-performs physics-based (if not just totally automatic, without guess) approaches like these.
For some context: acidic or basic amino-acids in proteins may have different pKi depending on their local environment (i.e. if they are in contact with solvent or not, and residues surrounding them), which is why "adding hydrogen to proteins crystallographic structures at a constant pH" is not as easy as it sounds.

source: Stack Exchange

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

sort: votes


text:
Best way to dock subset of protein PDB and definition of a ligand
I have extracted residue 45 to 88 from a PDB file (2YRQ).
With that I have a new smaller PDB file. Let's call it 45_88.pdb
And I want to dock it to some other full protein.
My question is what's the best way to dock it. I'm considering Autodock Vina.
But I'm not sure if it's appropriate since it is for ligand-protein docking.
What's the definition of a ligand? Can a protein with only 42 residues be considered a ligand?

I would advise using the term peptide more so than ligand, which also includes, and generally is used to refer to, small molecules.

Answer from @pippo1980, converted from comment:
https://pubs.acs.org/doi/10.1021/acs.jctc.9b01208 Comprehensive Evaluation of Fourteen Docking Programs on Protein–Peptide Complexes

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


text:
What does the IGV refseq genes with three different thicknesses of lines mean?
What does the IGV refseq genes with three different thicknesses of lines mean?

The thinkness indicates whether something is an intron, a coding exon, or a UTR. The very thin lines are introns, the full thickness boxes are coding exons and the thinner boxes are UTRs.
UTRs, untranslated regions, are exons, they are included in the final mRNA, but they don't code for protein. So, by convention, they are shown as thinner exons to give the user a visual indication that although this is still an exonic region, it's a special one.
I don't really know where this is documented, it's pretty standard across all genome browsers and bam visualizers I know though so it's a common idiom.


text:
How to decide whether to assume 2 Gene Losses or 1 to find minimum cost trees for Gene Duplication/Expression events?
This paper introduces an approach to estimate a convergence from discord between phylogenetic species trees and gene trees, using Gene Duplication and Expression events.
A Gene Duplication (GD) event refers to when a bifurcation of the phylogenetic gene tree occurs creating two duplicate gene loci. The first image marks them in solid white diamonds.
A Gene Expression (GE) event refers to the loss or Regulatory Mutation (RM) of genes. The types are:

Gene Loss (marked by the cross in the first image) refers to when a locus is lost/deleted along the lineage of an extant species.
Gene Inactivation (marked by the unfilled downward caret in the first image) refers to a Regulatory Mutation where a gene is rendered inactive.
Gene Reactivation (marked by the unfilled upward caret in the first image) refers to a Regulatory Mutation where a previously inactivated gene is reactivated and sampled in extant species.

Now, in the event of a discord between a gene tree and a species tree, we can attempt to establish concordance through assumptions about hypothetical duplication and expression events. The gene tree with a minimum count of (GD + GE) is the most likely assumption for concordance.
In the first image, the right locus resulting from the second bifurcation of the gene tree is sampled in extant species C and D and undergoes a loss event for species E. I would assume that it would also undergo a loss event for species A and B but the paper merges those two loss events into one.
On the other hand, based on the second image from this paper, I've simulated the Gene Duplication/Expression events on the gene tree and found that if such loss events are merged into one, the Duplication/Expression (D/E in the table) count I calculate is lower than what is mentioned in the table.
I'm confused as to whether there is an appropriate approach to this. The author of the first paper mentions that such assumptions regarding Gene Duplication and Expression events are not at all representative of real Gene Duplication and Expression events (established through observed evidence) and serve as a mechanism to approach likely solutions to discordance between species trees and gene trees.
My question is, how do I arrive at the minimum cost tree?


text:
Illumina vs PacBio vs Nanopore price
Can someone help me with comparing the cost per run of Illumina, PacBio and Nanopore sequencing technologies?

Voted to keep open. I honestly think, although it's a wide question, it's essential for planning of any genomics project and it's not like these comparisons don't exist!

Could you please add more context to this, or narrow your question? Searching Google Scholar for illumina nanopore pacbio turns up thousands of comparative results papers, and StackExchange doesn't work well with thousands of answers.

You don't think that the kind of experiment you want to run is at all relevant?

voted to keep open, agree with Kamil S Jaron. This is information that is often intentionally obscured to the consumer by sequencing providers, but is essential for anyone trying to plan. compiling information is important for the community.


text:
GWAS : input files for rrBLUP package in R
I want to do a verification of the result of an RNAseq analysis of human cancer samples by GWAS, I want to use rrBLUP package in R, so I need three input file genotype, phenotype and mapping file, so I am confused about how I can get these files.


text:
Grep peculiar pattern inside a file
I would like to grep everything inside the second square bracket of a file:
>AB1311658.1 unnamed protein product [C.sp. balc.]

>evm.model.pilon.8 [actin][D.r.]

>NP_0232343.1 [fucose mutarotase] [E.l.]

Obtaining:
C.sp. balc.

D.r.

E.l.

How can I solve it?
Marco

Is that everything you have in the file? No other lines? And what operating system are you using?

You can easily install GNU version of the commands into your Mac. Install the GNU Command Line Tools. You can do this using Homebrew, then
brew install coreutils. Default will be MacOS version (unless you intentionally set otherwise), but when you want to use GNU version, just add a "g" before the command. For example, "sort" vs "gsort"; "grep" vs "ggrep". It will make your life easier as a Mac user.

One way using sed:
sed -n '/^>/s/.*\[\([^]]*\)\]$/\1/p' input.fa

Results:
C.sp. balc.
D.r.
E.l.

Explanation:

Turn off automatic printing of the pattern space using -n
Address lines starting with a header using /^>/
Perform a substitution using a capture group to select everything between the last set of brackets on these lines, with s/.*\[\([^]]*\)\]$/\1/
Print the pattern space with p

You could probably simplify your sed to sed -n '/^>/{s/.*\[//; s/]//p}' file.

The /^>/{ doesn't appear to be necessary.

With GNU grep (the default on Linux):
$ grep -oP '.*\[\K[^]]+' file
C.sp. balc.
D.r.
E.l.

The -o means "print only the matched portion of the line" and -P enables Perl Compatible Regular Expressions which give us \K for "ignore everything matched up to this point". So we then look for the longest possible string ending in a [ and discard it (.*\[\K), and then look for the longest possible stretch of non-] characters after that.
For very large files, you can speed things up a little by focusing only on the header lines:
grep -oP '^>.*\[\K[^]]+' file

For OSX, I don't have a -P option (Perl pattern matching) so the best I can do is finish the one-liner by piping into perl
grep -o '\[.*\]\[[^]]*' input.fa | perl -p -e 's/[][]/ /g'

The -o makes grep match the substring, rather than the entire line.
With GNU grep (the default on Linux) and the -P option, it's a single line of grep where \K or (?=pattern), (?<=pattern) are used for non-capture characters and these are placed around the [
grep -oP '(?<=\[).*(?=\]).*(?<=\[)[^]]*' input.fa

The ?<= is a "positive lookbehind" and the ?= is a "positive lookahead".
It's a lot easier to use a perl one-liner.

You could also pipe to sed, and skip the Perl dependency: grep -o '\[[^]]*\]$' input.fa | sed 's/[][]//g'

@Steve what benefit would that give? You just added sed as a dependency instead. Both sed and perl are installed by default on the vast majority of systems, so unless you are using AIX or a very limited embedded system, perl is just as valid as sed.

@terdon I agree, but nearly all systems with grep will have sed.

@Steve yes, exactly like perl. Unless people are trying to do bioinformatics on toasters or something, I really don't see any point in using sed instead of perl. They are both very powerful, standard, ancient commands that have a long history in the *nix world in general and in bioinformatics in particular. The only full system (as opposed to embedded system) I know of that doesn't install perl by default is AIX and that isn't likely to be used for bioinformatics.

Perl can easily be downloaded @Steve ;), "Dockerised", or "conda-ised". Not many Linux OS have removed Perl, FreeBSD is one. Most of us work with Ubuntu, Debian or Darwin. I think even Alpine Linux has Perl.

@terdon Well, Perl is a sed-killer. For simple substitutions like the above, PCRE are not needed. You can use Perl if you want to, but you could also use sed and not need to escape the brackets.

@M__ Alpine does not come with Perl out of the box

@Steve well yes, but Alpine is mostly for embedded systems, right? Not sure what you mean about escaping, neither tool requires you to escape the brackets. And sure, of course you can use sed. I just don't see any benefit in using one external tool over another here. Nothing wrong with using sed, but nothing wrong with using perl either. Granted, if I bring out perl, I would probably do the whole thing in perl and remove the grep, but still.

@Steve yes technically Alpine Linux holds sed and grep but via Busybox and its grep is basic without -P. I installed perl on Alpine it took a few seconds  maximum. However, the reason is the container is 7.5 MG and Alpine Perl is 55MB. The key to the question is GNU grep.

@terdon You're right, neither tool requires the brackets to be escaped but 's/[[]]/ /g' should be 's/[][]/ /g I think.

@Steve eeeek! You're quite right, thanks!

@Steve here's an article specially written for you ;) https://stackoverflow.blog/2022/07/06/why-perl-is-still-relevant-in-2022/

Thanks @M__ - I hadn't even heard of Perl 7 until now. But then I haven't really touched Perl for nearly a decade now. I just hope you're not confusing me for someone who dislikes Perl - that's certainly not the case. It's just that it's not needed here in my opinion. The Perl dependency also makes pipelining the solution tricky using the standard (i.e. non-mulled) quay.io biocontainers which provide the busybox suite. But that's not likely even a concern here..

@Steve "not confusing me for someone who dislikes Perl " ... just checking ;) Fun as this is the last 3 comments at least need to be deleted, probably more

@terdon Alpine is a "freak" Linux OS, just for Docker containers. I should have checked my information before introducing it here :)

Some of the solutions suggested above utilize GNU grep. If you are a Mac user, you can easily install GNU version of the commands into your Mac. You can do this by installing the GNU Command Line Tools. The easiest way is to use Homebrew, then brew install coreutils.
The default commands will still be MacOS version (unless you intentionally set otherwise), but when you want to use GNU version, you could just add a g before the command. For example, sort vs gsort, grep vs ggrep. You can also check manual of the command the same way, e.g. man gsort. Trust me this will make your life easier as a Mac user.
Also, in recent years default MacOS Terminal shell is zsh. If you desire you can switch it back to bash.

Thanks @Supertech perfect


text:
Help with DIA-Mass Spectrometry data analysis with several conditions (limma?)
I am trying to learn how to analyse normalised DIA-MS data and I am struggling with it ://
The original dataset I got is (6 conditions (2 samples each)) with 3 technical replicates (total: 36 sample columns x 1128 proteins in rows). I removed the proteins with more than 80% of NA, and impute with kNN the rest of missing values.
Then I collapse the replicates doing the average of the 3 replicates (is this correct?)
So now, I have a dataframe with 12 columns x 1072 rows.
The question is... I am not sure how to perform the differential analysis from here (assuming I did well the previous steps lol)
Could someone help me with the tools or code to compare all the 6 conditions between them?
I read something about Limma package, but not sure if it can be used in this DIA-ms data and neither sure of how to use limma for it (I already check the vignette with no succeed)
UPDATE: I tried to use Limma as follows:
fit1 <- lmFit(norm.data_collapsed, design)
fit2 <- contrasts.fit(fit1, contrasts = contrast)
fit3 <- eBayes(fit2)
limma.res.01 <- topTable(fit3, coef = comb.contrast[1], n= Inf)

This returns me a non-sense values for logFC or significance:

Any help would be very welcome,
Best regards

Regarding the big FC values, I just realise the Limma lmfit function requires as input a log2 matrix... so I just transformed my data and now makes more sense. Still I have the doubt regarding the collapsing technical replicates and and if some one can verify the pipeline would be like that.
Thanks!!!
Fer

Please could you mark this answer as accepted? That would be very welcome.


text:
Is there a publication database or search engine offering geneID or UniProtID correlation?
Do we have a publication database or search engine that offers geneID or UniProtID correlation?
For example, I search "IPS cells" -> I want to know what genes or proteins are under active research in IPS cells.

I think this could be close to what you're looking for: http://glad4u.zhang-lab.org/index.php
Paper: GLAD4U: deriving and prioritizing gene lists from PubMed literature
It works with gene names though, so I am not sure it's a 100% fit.

Take a look at GeneWalk, https://churchman.med.harvard.edu/genewalk. It's a very interesting solution to this rather difficult problem. It is well documented and packaged. But It's open source software -- not a service. You'll have to run it on your own system.


text:
How to download data from SRA in Linux systems via the command line?
My workflow for downloading data from SRA has been the following:

Access SRA Run Selector.
Enter the accession number for the project of interest.
Download "Accession List" for the "Total" data. The accession list is basically, as the name says, a .txt file listing all accession numbers for each run in a given bioproject.
For each accession number in the accession list I run the following command:

$ fastq-dump --accession <accesion_number> --outdir .

How to download data from SRA in Linux systems via the command line? What are fast and scalable ways of achieving this? What should one know before downloading data?
Although there are some similar discussions in this StackExchange (eg. Trying to retrieve sra data using sra explorer, Downloading dataset from SRA (SOLiD Platform), Downloading SRA Files from AWS, etc).
The discussion on What's the best way to download data from the SRA? Is it really this slow?
 seems like a good starting point, but I wasn't able to find a straightforward guide on how to download data from SRA, e believe that a well thought and complete answer would help beginners alike. I am preparing an answer to the problem myself, but I thought that the community would probably give some interesting insights (:

I wouldn't say the download itself is slow but in my experience sra to fastq conversion took much more time than the download, and using fasterq-dump instead of fastq-dump did not help much. I would take a look at ENA (https://www.ebi.ac.uk/ena/browser/home) and see if they offer fastq files for the dataset(s) you are interested in.

The dataset I am looking for is actually available at ENA, and indeed, the workflow is way nicer using their interface.

So far, grabseqs is the easiest option, a wrapper for fastq_dump and fasterq-dump, you can install it with conda, I recommend you to use an environment:
conda create -n grabseqs -c louiejtaylor -c bioconda -c conda-forge grabseqs

conda activate grabseqs

grabseqs sra SRRNNNNNNNN

conda deactivate

And download directly projects, runs, and bioprojects from SRA with grabseqs sra SRRNNNNNNNN. Check it out here: https://github.com/louiejtaylor/grabseqs
You can install in a different way just in case the first one doesn't work:
conda create -n grabseqs -yc bioconda parallel-fastq-dump pigz
conda activate grabseqs
pip install grabseqs

grabseqs sra SRRNNNNNNNN

conda deactivate

Nice, I will test this in a couple of hours ^^

I didn't manage to install it easily, it is presenting some kind of conflict between the dependencies. I will try to find something else.

I did a slight modification just in case the official installation steps didn't work.


text:
How to calculate switch errors on plink formatted phased genotype data?
everyone!
I'm working with haplotype inference algorithms in pedigree genotype data, and I would like to know if there is a tool or script to calculate switch errors on phased pedigree data (Plink formatted)?
I already know that VCFtools and other scripts do exist, but they only work with VCF format.
Note: I tried to convert my phased pedigree data to VCF but I lost phase, so, this approach does not seem to be helpful, that's why I'm looking for a specific tool to detect switches on phased pedigree data.
A switch errors ocurr when heterozygous positions are misassigned relative to the previous heterozygous position.
See visual example (block 2), True haplotypes up, inferred haplotypes bottom:
https://www.researchgate.net/publication/310766263/figure/fig13/AS:431688409980932@1479933945501/Illustration-of-switch-errors-and-block-breakpoints-In-haplotype-block-2-there-are-two.png
In my case, I got the phased pedigree plink file as follows:
Infered Haplotype, individual 1:
1|2 2|1 1|2 1|2 1|2
Real Haplotype, individual 1:
1|2 1|2 2|1 1|2 1|2
We can see 2 switches, between position 1-2 and 3-4.
Thank you very much in advance!
Regards!

Can you explain what switch errors are?

A switch errors occur when heterozygous positions are misassigned relative to the previous heterozygous position.
See visual example (block 2), True haplotypes up, inferred haplotypes bottom:
https://www.researchgate.net/publication/310766263/figure/fig13/AS:431688409980932@1479933945501/Illustration-of-switch-errors-and-block-breakpoints-In-haplotype-block-2-there-are-two.png
I got the phased pedigree plink file as follows:
Inferred Haplotype, individual 1:
1 1|2 2|1 1|2 1|2 1|2
Real Haplotype, individual 1:
1 1|2 1|2 2|1 1|2 1|2
See 2 switches, between pos. 1-2 and 3-4.


text:
Unable to execute DSSP algorithm from BioPython
I have the following settings in PyCharm:

BioPython is properly installed:

dssp_test.py
from Bio.PDB import PDBParser
from Bio.PDB.DSSP import DSSP
p = PDBParser()
structure = p.get_structure("16PK", "16pk.pdb")
model = structure[0]
dssp = DSSP(model, "16pk.pdb")
a_key = list(dssp.keys())[2]
print(dssp[a_key])

Output
C:\ProgramData\Miniconda3\python.exe C:/Users/pc/source/repos/DSSP_algoritm_test/dssp_test.py
Traceback (most recent call last):
  File "C:\Users\pc\source\repos\DSSP_algoritm_test\dssp_test.py", line 6, in <module>
    dssp = DSSP(model, "16pk.pdb")
  File "C:\ProgramData\Miniconda3\lib\site-packages\Bio\PDB\DSSP.py", line 437, in __init__
    dssp_dict, dssp_keys = dssp_dict_from_pdb_file(in_file, dssp)
  File "C:\ProgramData\Miniconda3\lib\site-packages\Bio\PDB\DSSP.py", line 231, in dssp_dict_from_pdb_file
    p = subprocess.Popen(
  File "C:\ProgramData\Miniconda3\lib\subprocess.py", line 947, in __init__
    self._execute_child(args, executable, preexec_fn, close_fds,
  File "C:\ProgramData\Miniconda3\lib\subprocess.py", line 1416, in _execute_child
    hp, ht, pid, tid = _winapi.CreateProcess(executable, args,
FileNotFoundError: [WinError 2] The system cannot find the file specified

Process finished with exit code 1

Why am I unable to execute DSSP algorithm from BioPython?

Just in the offchange, you install the DSSP binary (https://swift.cmbi.umcn.nl/gv/dssp/) as described in the BioPython documentation, right?

@MatteoFerla, Yes. Check the newly added screenshot.

Installing biopython does not magically install DSSP.

This looks like a python problem and not a “biopython”problem.
Check your file structures. Try giving a full path for 16pk.pdb
These might help. https://stackoverflow.com/questions/17658856/why-am-i-getting-a-filenotfounderror
https://stackoverflow.com/questions/18757127/windowserror-error-2-the-system-cannot-find-the-file-specified-is-not-resol

Up-voted to offset a highly incorrect down-vote.

Ok so this might not provide a solution, but here's a couple things someone might find useful;

First, Biopython does not automatically contain DSSP, which you have to download as mentioned in the comments below (from here swift.cmbi.umcn.nl/gv/dssp)
Secondly, DSSP within Biopython, for whatever reason, accepts the path to the file, not the file itself. I see you have that correctly (dssp = DSSP(model, "16pk.pdb")), but just pointing it ou
Third, DSSP originally had dssp executable, but the newer version has mkdssp executable, so you might need to do dssp = DSSP(model, "16pk.pdb", dssp="mkdssp"), so it knows what to associate it with. This was the problem I had, because without specifying the dssp, it is not able to launch DSSP.


text:
Can the information regarding atomic bonds be interpreted only from atom lines in a PDB file?
The following are the first 14 atoms from the protein 1CRN:
ATOM      1  N   THR A   1      17.047  14.099   3.625  1.00 13.79           N
ATOM      2  CA  THR A   1      16.967  12.784   4.338  1.00 10.80           C
ATOM      3  C   THR A   1      15.685  12.755   5.133  1.00  9.19           C
ATOM      4  O   THR A   1      15.268  13.825   5.594  1.00  9.85           O
ATOM      5  CB  THR A   1      18.170  12.703   5.337  1.00 13.02           C
ATOM      6  OG1 THR A   1      19.334  12.829   4.463  1.00 15.06           O
ATOM      7  CG2 THR A   1      18.150  11.546   6.304  1.00 14.23           C

ATOM      8  N   THR A   2      15.115  11.555   5.265  1.00  7.81           N
ATOM      9  CA  THR A   2      13.856  11.469   6.066  1.00  8.31           C
ATOM     10  C   THR A   2      14.164  10.785   7.379  1.00  5.80           C
ATOM     11  O   THR A   2      14.993   9.862   7.443  1.00  6.94           O
ATOM     12  CB  THR A   2      12.732  10.711   5.261  1.00 10.32           C
ATOM     13  OG1 THR A   2      13.308   9.439   4.926  1.00 12.81           O
ATOM     14  CG2 THR A   2      12.484  11.442   3.895  1.00 11.90           C

These fourteen atoms are members of two THR residues.
Can the information regarding atomic bonds be interpreted only from these fourteen lines?
For example,

the first atom, N, has three bonds attached to CA.
the second atom, CA, has four bonds; three are attached to N, and one is to C.
the third atom, C, has four bonds; one is attached to CA, one to oxygen, and two are empty.
the fourth atom, O, has two bonds; one is attached to the third carbon, and one is to the fifth atom beta-carbon.
the fifth atom, beta-carbon, has four bonds; one is attached to the fourth atom, two are attached to the first branch (OG1), and one is attached to the second branch (CG2).

Or do I have to use my background knowledge of amino acids to find bonds?
Note: I am not talking about hydrogen bonds. I am talking about atomic bonds only in the atoms in the residue.

The official list of bonds in each residue in the PDB is in the Chemical Component Dictionary.
You could also use interatomic distances to infer connectivity, but that's less reliable.

You can use topology files to extract bond information.
Here is an example at CHARMM36


text:
How can I calculate the radius of a residue?
I want to read a PDB file and display/render the atoms using a graphics engine.
If I consider each residue in the protein as a solid sphere, what would be the radius of the sphere?
How can I calculate it? Or, is there any database for that?

You may want to check how existing molecular viewers (PyMOL, Jmol, Chimera, web-based viewers on RCSB, PDBe and PDBj websites) render it.

@marcin, I am developing my own software for learning ourposes.

Van der Waals radii are used to construct a special kind of molecular model called a space-filling model. These models are constructed by drawing each atom as a van der Waals sphere with the atom’s nucleus at the center of the sphere. taken from https://www.reed.edu/chemistry/ROCO/Geometry/space_filling.html not sure how to expand it to a full residue

@pippo1980, thanks! good information. :)

try this paper not sure why 5A radius is choosen https://www.researchgate.net/figure/Each-protein-is-modeled-as-a-union-of-spheres-of-radius-5-A-one-sphere-per-amino-acid_fig2_8614866 the picture seems OK to me tough

ok full article is for free; it states: "....However, it has been suggested that the 5 Å radius is appropriate to model the hydration shell around the protein
[19]. The problem of computing
r for such a model reduces
to computation of the distance to the surface of each sphere,
and taking the minimum over all spheres.... " Worth to notice article 19 its about globular protein, no idead if 5A rules could be extended to other kind of protein too

The protein is
represented as a union of spheres, each of radius 5 Å, and
centered on the C
a atom of each amino acid residue


text:
How can I compute energy from a protein PDB file?
Suppose I have a protein PDB file and want to compute energy using the following formula from the coordinates in the file.

How can I do that?
Sometimes people suggest using the CONECT field for that purpose. However, some PDB files do not have this field.
Example: 4OSK.pdb (Crystal structure of TAL effector reveals the recognition between asparagine and guanine)

https://mattermodeling.stackexchange.com/questions/11678/how-can-i-retrieve-bond-information-froma-protein-pdb-file

https://www.ebi.ac.uk/pdbe-srv/pdbechem/bond/list/CYS Molecule : CYS
Bonds between two atoms of the molecule , then of course in each protein chain you'll have a Nterminal and Cterminal residues and all the ones in between having a peptide bond CO-N between Cter and Nter

http://ligand-expo.rcsb.org/pyapps/ldHandler.py?formid=cc-index-search&target=cys&operation=ccid here with pics http://ligand-expo.rcsb.org/reports/C/CYS/CYS_D3L1.gif

think you'll need to create a dictionary for each residue

see here for better explanation https://chemistry.stackexchange.com/a/162689/118509


text:
Why are DSSP and STRIDE considered reference algorithms, and why do people continue to use them?
I am working on globular protein secondary structure assignments. The purpose of my project is to design an improved algorithm for secondary structure assignment.
I know two tools: DSSP and STRIDE for secondary structure prediction and assignment.
However, their accuracy [the Q3 (Q8) accuracy is defined as the percentage of residues for which the predicted secondary structures are correct] is around 53% and 70%, respectively, which isn't much if we compare them to modern algorithms like PALSSE, etc.
If so, why are DSSP and STRIDE considered reference algorithms, and why do people continue to use them?

This question is not specific enough. Please add more context to your post to describe a specific bioinformatics problem. Generic "What's the best tool" questions are better asked in discussion-type places like Reddit.

Note OP has edited the question.

What does accuracy mean in this context? Comparison to human assignment? 53% for DSSP seems low to me.

@jgreener, The Q3 (Q8) accuracy is defined as the percentage of residues for which the predicted secondary structures are correct.

But what is the source of ground truth? DSSP and STRIDE assign secondary structure from tertiary structure rather than predict it from sequence. I see the number 53% for DSSP in the PALSSE paper, but that is not for Q3 or Q8 accuracy.

@jgreener, then make some correction in the OP, and post an answer.


text:
What does it mean by a Local Coordinate in the case of protein chains?

High-performance transformation of protein structure representation from internal to Cartesian coordinates

In Figure 3a, each residue has its own local coordinates, which are illustrated shown using different colors.

Here are two more papers:

Generative De Novo Protein Design with Global Context

Generating 3D Molecules for Target Protein Binding

By global coordinate, I understand the coordinate system in which we view the entire protein chain.
What does it generally (not only in this particular paper) mean by a Local Coordinate in the case of protein chains?
Why was it necessary to introduce this terminology when they could have just used the term "vector"?
Can you supply me with an example with an explanatory diagram?

For proteins and molecules in general, and within the meaning of this paper, "internal coordinates" are any representation that are not dependent on the molecule rotation or translation in a Cartesian space.
Calling them "internal coordinates" may be slightly confusing, and a more correct term could be "internal representation"
That means, I guess, that internal coordinates could be, as you mention, a list of vector in a Cartesian space, where the Cartesian space "orientation" would be decided by, say, principle component analysis (or any other arbitrary rule depending on the molecule). This however defeats the purpose of most internal representations: dimensionality reduction or faster computations.
Other (more interesting) ways of representing a molecule in internal coordinates would be:

As (Bayati, Lardan & Bardhan 2020) focus on: a set of bond distance, angle, and dihedrals (see Ramachandran plots and rotamers). This is especially useful in the case of proteins. This representation is for example used by Rosetta. Angles and bond distance can be considered as constant (depending on the atomic type), and so one can infer the Cartesian coordinates of multiple atoms with just one dihedral angle, hence providing a huge dimensionality reduction. In this case there is no need for vector. Also, in the case of internal energy computation there is always one energy term depending on dihedrals, so that is already done. Limit of such methods is that in molecular constructs there are energy terms due to long-range interactions, i.e van der waals and electrostatics. There is hence a need in such software to always switch between a Dihedral representation and a Cartesian representation. Another problem is that proteins backbone are then represented as a dihedral "chain", and so the Cartesian coordinates of one atom becomes tightly linked to the previous one, which are linked to the previous one, etc. This is the problem this paper is trying to deal with.

As (Tan et al. 2023) use, a graph representation. Graphein is an interesting python library for converting and analyzing proteins and molecules from cartesian coordinates to graph representation. In this case instead of having a vector, each atoms position is defined as its distance with surrounding atoms. This kind of representation has received more and more attention from the protein structure prediction & design community as it is particularly suited for neural networks. For example I would say the protein representation in AlphaFold is a kind of graph representation (using "triangles" made by the distance between three atoms). In comparison to dihedral representation, graph representation packs less "assumption" about the system and are even worst than Cartesian representation in terms of dimension (if you do a "global" representation of N atoms then you end up with N! values vs only 3*N with Cartesian), but at least things like non-covalent interactions are implicitly represented.

I believe this answer gives enough material to look at. If this still appears unclear, skimming through Rosetta tutorials / papers figures / webinars should provide ample explanations, same for AlphaFold.


text:
How can I use REMO tool to generate seq.dat and seq.txt files?
I am trying to generate seq.dat and seq.txt files myself, so that I can use them with REMO.
I am using mkseq.pl script in REMOv2.
The following was my command:
user_name@server_name:~/spatial/reconstruction_tools/remo_v2$ perl mkseq.pl /home/spatial/input_pdb/pdb1a6j.pdb 1
PSIPRED finished
user_name@server_name:~/spatial/reconstruction_tools/remo_v2$

This command generated four files:

seq.dat
protein.fasta
psitmp27331007f0101.blast
psitmp27331007f0101.fasta

seq.dat is totally empty. psitmp27331007f0101.blast is also practically empty. There is no seq.txt.
What am I doing incorrectly?


text:
How do I do "protein culling"?
I am working on a novel protein secondary structure detection algorithm using machine learning.
As far as I understand, online protein databanks (e.g., RCSB, etc.) have a huge collection of proteins (several millions). All proteins are not useful for research purposes as they contain redundant, duplicate, or incomplete information.
So, we have to choose the minimum number of proteins so that we have enough proteins to cover all the required information for our specific purpose.
This is called "protein culling."
Google search shows that there were some servers like PISCES, ProDaMa, etc. However, as usual, all of them are down.
Now the question is, how do I do protein-culling myself in the context of protein secondary structure detection?
N.B. Proteins with similar folds or secondary structure elements can be considered redundant, even if their sequence identity is low.

How good are your data engineering skills? Removing duplicate, redundant or incomplete stuff isn't hard and a database of several million is honestly pretty small. You can probably do this in pandas2 definitely polars, but after that it gets more specialised.

@M__, I work on secondary and tertiary structures. I need a minimum set of individual PDB files.

@M__, PIESCES server had 5022 PDB files.

It's still applicable to data engineering for the criteria you are considering. I agree that it may not fit inside pandas2. If you were considering genetic distance ... thats more complicated.

This depends entirely on what, exactly, you are trying to do. Please [edit] your question and specify the biological question you are trying to answer.

@terdon, I am working on a novel protein secondary and tertiary structure detection algorithm using machine learning.

OK, please add that to your question but, again, what do you consider not useful here? Presumably, you don't care so much about sequence identity as about tertiary structure similarity, right? So what would you like to cull and what would you like to keep?

@terdon, Proteins with similar folds or secondary or tertiary structure elements can be considered redundant, even if their sequence identity is low.

Have you considered looking at this in terms of protein domains instead of proteins? Say starting with the Pfam data and getting representative examples of each domain? (random shot in the dark here, I have no experience with tertiary structure analysis)

@terdon, I am a beginner at level zero in Bioinformatics.

If you were willing to download the data yourself, you could use a tool like FoldSeek to cluster on various thresholds to remove redundancy. Though I realise this is tertiary structure, rather than secondary structure, but it still could work with some tweaking.


text:
Why am I not being able to detect all the bonds in a residue?
Let us consider the first seven atoms of the  protein 1CRN.pdb:
ATOM      1  N   THR A   1      17.047  14.099   3.625  1.00 13.79           N
ATOM      2  CA  THR A   1      16.967  12.784   4.338  1.00 10.80           C
ATOM      3  C   THR A   1      15.685  12.755   5.133  1.00  9.19           C
ATOM      4  O   THR A   1      15.268  13.825   5.594  1.00  9.85           O
ATOM      5  CB  THR A   1      18.170  12.703   5.337  1.00 13.02           C
ATOM      6  OG1 THR A   1      19.334  12.829   4.463  1.00 15.06           O
ATOM      7  CG2 THR A   1      18.150  11.546   6.304  1.00 14.23           C

This represents an amino acid called Threonine.
The general bond diagram for Threonine is as follows:

So, it has at least six bonds (I am actually not sure!) which are also evident from the seven ATOM lines.
Each atom has its radius. In general, if the distance between two atoms is less than the sum of the radii, there is a bond between these two atoms.
Based on the above information, I tried to identify the bond information, and I could detect only two:
     atom1           atom2        distance       atm1 Radi      atom2 Radi        sum radi        is bond?
         N              CA            1.50            0.65            0.70            1.35           False
         N               C            2.44            0.65            0.70            1.35           False
         N               O            2.67            0.65            0.60            1.25           False
         N              CB            2.48            0.65            0.70            1.35           False
         N             OG1            2.75            0.65            0.60            1.25           False
         N             CG2            3.86            0.65            0.60            1.25           False
        CA               N            1.50            0.70            0.65            1.35           False
        CA               C            1.51            0.70            0.70            1.40           False
        CA               O            2.36            0.70            0.60            1.30           False
        CA              CB            1.57            0.70            0.70            1.40           False
        CA             OG1            2.37            0.70            0.60            1.30           False
        CA             CG2            2.61            0.70            0.60            1.30           False
         C               N            2.44            0.70            0.65            1.35           False
         C              CA            1.51            0.70            0.70            1.40           False
         C               O            1.24            0.70            0.60            1.30            True
         C              CB            2.49            0.70            0.70            1.40           False
         C             OG1            3.71            0.70            0.60            1.30           False
         C             CG2            2.98            0.70            0.60            1.30           False
         O               N            2.67            0.60            0.65            1.25           False
         O              CA            2.36            0.60            0.70            1.30           False
         O               C            1.24            0.60            0.70            1.30            True
         O              CB            3.12            0.60            0.70            1.30           False
         O             OG1            4.34            0.60            0.60            1.20           False
         O             CG2            3.74            0.60            0.60            1.20           False
        CB               N            2.48            0.70            0.65            1.35           False
        CB              CA            1.57            0.70            0.70            1.40           False
        CB               C            2.49            0.70            0.70            1.40           False
        CB               O            3.12            0.70            0.60            1.30           False
        CB             OG1            1.46            0.70            0.60            1.30           False
        CB             CG2            1.51            0.70            0.60            1.30           False
       OG1               N            2.75            0.60            0.65            1.25           False
       OG1              CA            2.37            0.60            0.70            1.30           False
       OG1               C            3.71            0.60            0.70            1.30           False
       OG1               O            4.34            0.60            0.60            1.20           False
       OG1              CB            1.46            0.60            0.70            1.30           False
       OG1             CG2            2.54            0.60            0.60            1.20           False
       CG2               N            3.86            0.60            0.65            1.25           False
       CG2              CA            2.61            0.60            0.70            1.30           False
       CG2               C            2.98            0.60            0.70            1.30           False
       CG2               O            3.74            0.60            0.60            1.20           False
       CG2              CB            1.51            0.60            0.70            1.30           False
       CG2             OG1            2.54            0.60            0.60            1.20           False
2
Done!

How can I detect all bonds?
EDIT: As per pippo1980's request:
public class ADTResidue: Dictionary<int, ADTAtom>, IEquatable<ADTResidue>
{
    public string ProteinName { get; set; }
    public string ChainName { get; set; }
    public int ResidueNo { get; set; }
    public Point3d Coordinate { get; private set; }
    public string Code { get; set; }
    public string Full { get; set; }
    public Color ColorRGB { get; set; }

    public ADTResidue(string code, string full, Color colorRgb)
    {
        Coordinate = new Point3d();

        Code = code;
        Full = full;
        ColorRGB = colorRgb;
    }

    private readonly Dictionary<int, CovalentBond> BondDictionary;

    public List<CovalentBond> GetBonds()
    {
        return new List<CovalentBond>(BondDictionary.Values);
    }

    public ADTResidue()
    {
        BondDictionary = new Dictionary<int, CovalentBond>();
    }

    public void CalculateBonds()
    {
        foreach (ADTAtom atom1 in this.Values)
        {
            foreach (ADTAtom atom2 in this.Values)
            {
                if (atom2.Radius <= 0)
                {
                    throw new Exception("zero radius in " + atom2.Name);
                }

                if (atom1 != atom2)
                {
                    double distance = Point3d.Distance(atom1.Coordinate, atom2.Coordinate);

                    bool exists = ABond.IsBondExists(atom1.Coordinate, atom2.Coordinate, atom1.Radius, atom2.Radius);

                    if (exists)
                    {
                        int hashCode = atom1.GetHashCode(atom2);
                        if (!BondDictionary.ContainsKey(hashCode))
                        {
                            BondDictionary.Add(hashCode, new CovalentBond(atom1, atom2));
                        }
                    }
                }
            }
        }

        int count = BondDictionary.Count;

        Console.WriteLine(count);
    }

    #region equality comparison
    // IEquatable<> implementation...
    #endregion
}

public abstract class ABond: Pair<ADTAtom, ADTAtom>
{
    public ABond():base()
    {
    }

    public ABond(ADTAtom first, ADTAtom second) : base(first, second)
    {
    }

    public static bool IsBondExists(Point3d first, Point3d second, double radius1, double radius2)
    {
        double distance = Point3d.Distance(first, second);

        return (distance <= (radius1 + radius2));
    }
}

public class CovalentBond : ABond
{
    public CovalentBond():base()
    {
    }

    public CovalentBond(ADTAtom first, ADTAtom second) : base(first, second)
    {
    }
}

would you share your code with us ?

@pippo1980, its a part of a bigger framework. I am sharing a smaller part of it. Hope that helps.

unfortunately not, but could help others. Did you try running it on first T of 3NIR.pdb ?

@pippo1980, Yes. My code works now. Previously, I made mistakes with radii values.

is that above C#    ?

@pippo1980, Yes.

is it same as java ?

@pippo1980, Not same, but similar.

Your method doesn't work for two possible reasons;

The radius information may be incorrect.
Your molecule requires some minimization. This option helps to optimize distance between atoms.

You can find some information about the minimization at the following links

Wikipedia
Gromacs


text:
What do these characters mean in a PDB file?
The 33rd atom information in 3nir.pdb is written as follows:
ATOM     33  N  ATHR A   2       4.935 -11.632  15.046  0.74  2.95           N
How should I interpret N  A in the above line?
I mean, what does it mean? Is it Sodium or Nitrogen?
If it is Nitrogen, what does A mean then?

have a look here : https://www.wwpdb.org/documentation/file-format

from https://www.wwpdb.org/documentation/file-format

Protein Data Bank Contents Guide:
Atomic Coordinate Entry Format Description
Version 3.30
Document Published by the wwPDB
This format complies with the PDB Exchange Dictionary (PDBx)
http://mmcif.pdb.org/dictionaries/mmcif_pdbx.dic/Index/index.html.
I believe PDB format is kind of outdated, but not sure..
At page 177 it says about altLoc:

AltLoc is the place holder to indicate alternate conformation. The alternate conformation can be
in the entire polymer chain, or several residues or partial residue (several atoms within one
residue). If an atom is provided in more than one position, then a non-blank alternate location
indicator must be used for each of the atomic positions. Within a residue, all atoms that are
associated with each other in a given conformation are assigned the same alternate position
indicator. There are two ways of representing alternate conformation- either at atom level or at
residue level (see examples).
For atoms that are in alternate sites indicated by the alternate site indicator, sorting of atoms in
the ATOM/HETATM list uses the following general rules:
• In the simple case that involves a few atoms or a few residues with alternate sites, the
coordinates occur one after the other in the entry.
• In the case of a large heterogen groups which are disordered, the atoms for each
conformer are listed together.

your question could be more useful if you post the entire description for residue 2:
ATOM     33  N  ATHR A   2       4.935 -11.632  15.046  0.74  2.95           N
ATOM     34  N  BTHR A   2       5.078 -11.406  15.180  0.31  2.78           N
ATOM     35  CA ATHR A   2       5.757 -11.521  13.850  0.81  3.02           C
ATOM     36  CA BTHR A   2       5.773 -11.153  13.921  0.20  2.67           C
ATOM     37  C  ATHR A   2       7.070 -10.839  14.210  0.74  2.82           C
ATOM     38  C  BTHR A   2       7.155 -10.559  14.193  0.29  1.80           C
ATOM     39  O  ATHR A   2       7.152  -9.941  15.050  0.80  3.31           O
ATOM     40  O  BTHR A   2       7.214  -9.641  15.012  0.25  2.41           O
ATOM     41  CB ATHR A   2       4.976 -10.693  12.813  0.87  5.53           C
ATOM     42  CB BTHR A   2       4.896 -10.354  12.941  0.25 12.07           C
ATOM     43  OG1ATHR A   2       4.611  -9.432  13.388  1.00  6.88           O
ATOM     44  OG1BTHR A   2       3.743 -11.083  12.501  0.25  9.57           O
ATOM     45  CG2ATHR A   2       3.858 -11.584  12.293  0.75 10.03           C
ATOM     46  CG2BTHR A   2       5.683  -9.885  11.726  0.27  5.90           C
ATOM     47  H  ATHR A   2       4.547 -10.814  15.527  0.75  3.44           H
ATOM     48  H  BTHR A   2       5.510 -10.211  15.754  0.25  2.90           H
ATOM     49  HA ATHR A   2       5.962 -12.339  13.548  0.75  3.32           H
ATOM     50  HA BTHR A   2       4.036  -9.929  13.477  0.25  2.86           H
ATOM     51  HB ATHR A   2       5.648 -10.589  11.938  0.75  5.43           H
ATOM     52  HB BTHR A   2       4.644  -9.326  13.574  0.25  5.67           H
ATOM     53  HG1ATHR A   2       5.030  -9.344  14.216  0.75  8.74           H
ATOM     54  HG1BTHR A   2       3.236 -11.198  13.399  0.25 10.21           H
ATOM     55 HG21ATHR A   2       4.096 -12.441  11.924  0.75 10.92           H
ATOM     56 HG21BTHR A   2       6.542  -9.278  12.024  0.25  9.66           H
ATOM     57 HG22ATHR A   2       3.222 -10.974  11.650  0.75 10.92           H
ATOM     58 HG22BTHR A   2       5.039  -9.142  11.179  0.25  9.66           H
ATOM     59 HG23ATHR A   2       3.163 -11.738  13.200  0.75 10.92           H
ATOM     60 HG23BTHR A   2       5.904 -10.639  11.169  0.25  9.66           H


text:
What is protein secondary structure assignment?
I am absolutely new to Structural Bioinformatics (I only started last week). I am working on the secondary structure assignment/prediction (actually I am not sure) of proteins using ML.
When it comes to protein analysis, most of the books or research papers talk about protein structure Prediction:
https://bip.weizmann.ac.il/course/structbioinfo/pdf/SecondaryStructurePrediction_WIS.pdf
Only a handful of books or articles talk about protein structure Assignments:
https://www.biorxiv.org/content/10.1101/2021.02.02.429329v1
Can anyone tell me what is the basic difference between Assignment and Prediction?
Kindly explain with examples.

Do you have any examples for the use of Protein Structure Assignment? This paper uses the term in the title but they don't particularly distinguish the terminology.

To make it crystal clear (and make than pun):

Assignment: a 3D structure is known and the residues are assigned a secondary structure
Prediction: an algorithm predicts from linear (primary) sequence what the secondary structure may be —often incorrectly. This can be homology based, covariance, deep learning etc. But most commonly just the vector norm of window residue composition, e.g. proline is never in a helix and helices and sheets have different periodicity (two for a sheet and three and a bit for a helix).

What do you mean by 3D structure? Is it a Secondary or Tertiary structure?

Tertiary/quaternary structure —which may be of only part of the protein and with gaps due to technical limitations. But I mean specifically a list of atomic coordinates (eg. PDB file) of that chain, possibly as an oligomer, with a ligands or other chains as opposed to say just the name of the fold, e.g. Rossmann fold.

Answer from @jan-wilamowski, converted from comment:
That second paper seems to distinguish between the known data that is used for ML training (i.e. amino acid sequences are assigned some secondary structure) and the attempt to replicate this natural assignment by an "artificial" prediction.


text:
How can I search globular proteins by residue count in the RCSB databse?
I want to find some globular proteins with some of the longest individual chains that have all three types of secondary structures,, namely helix, sheet, and loop.
So, I tried the following search:

This search should have at least returned Titin.
However, it returns nothing.
How can I search for globular proteins by residue count in the RCSB database?

What makes you think that there is a >20,000 amino acid structure for TTN in PDB? The longest residue count for structures referenced from the Wikipedia page is 3,420 monomers (4UOW).

@gringer, see the update.

I'm aware that Titin is "the largest protein chain in your body", but I had trouble finding a PDB model for that complete protein chain. There does seem to be a supercomplex that includes Titin, though: https://www.rcsb.org/structure/8G4L

If I want to find the longest residue count model in PDB that includes Titin, I need to use the statistic "Total Number of Polymer Residues per Deposited Model":

For looking for helices, the statistic seems to be "Count Per Polymer Chain - Protein helices", and for beta sheets "Count Per Polymer Chain - Beta-sheets". These statistics can be combined using the AND operator.
Not quite sure what you mean by "loop" or "Coil". Disordered protein regions have no associated annotation. You can type words into the "Structure Attributes" search box to find statistics that can be searched and filtered.

Not quite sure what you mean by "loop". --- Coils (C).

This post doesn't answer my question. My aim is not to search for Titin. My aim is to find some of the proreins with the longest chains.

@user366312 As you were, I was using Titin as a demonstrative example for showing what I believe to be the appropriate residue statistic. If you don't think this answer is useful in solving your problem, then click on the "This Answer Is Not Useful" down arrow.


text:
What is the cutoff distance in this algorithm?
In the following article:

a fast and efficient Cα-based method for accurately assigning protein secondary structure elements

What are the minimum and maximum cutoff distances in this algorithm?
Also, what does it mean by:

I.e. where does (2x43)x3 come from?

The 43 features are listed in the supplementary materials on page 4. For example,
Feature Type Atoms Involved
1 Distance i, i-5
2 Distance i, i-4
3 Distance i, i-3
4 Distance i, i-2
...
11 Angle i-1, i, i+1
12 Dihedral i, i+1, i+2, i+3

Are some of them; so Feature 11 is the angle between atom(i-1), atom(i), and atom(i+1). Since it says these were calculated for both the C-alphas and pseudo-centers, this gives 2 * 43. Then this is multiplied by three as mentioned in your quoted paragraph.
These features (measurements, almost) are numbers, not booleans - so I do not think that they use distance cutoffs as such, although the random forest method is presumably learning various cutoffs and encodes those in the structure of the forest.


text:
Are the information stored in the PBD files in the protein data bank complete or incomplete?
I know that there are thousands of PDB files stored in Protein Data Banks.
Are all these files complete in terms of the information they store?
If YES, do the individual files get updated from time to time?
If YES, why/when do they get updated?

A PDB file is the last step in a series of operations.
For a crystal structure in very simple terms, it consists of the coordinates of the atoms placed to best fit the density map, which in turn (with these atoms due to the phase problem) is a conversion from reciprocal space which is the result of transforming the quantified dots from the pretty diffraction pattern picked up from the detector during the experiment.

Are all these files complete in terms of the information they store?

Not at all, a PDB file is not complete. It nearly complete when you count density in there.

If YES, do the individual files get updated from time to time?
If YES, why/when do they get updated?

Only, minor format upgrade, eg. https://www.rcsb.org/versions/1UBQ.
It is however very common to polish structures.
There is a tool called PDB-Redo that is commonly called in many pipelines for example to "upgrade" older files.
Some folk open the files in Coot and let it do its thing.
For most molecular mechanics experiments it is common to energy minimise the structure (freely, restrained by original position or restrained to density map). It is uncommon in several docking tutorials simply to keep it simple.
However, for most uses, if an aromatic ring is not quite planar or a bond is longer than most by a fraction of an Ångström, nothing happens.
Also worth remembering many structures are incorrect and not just imperfect: there may be actually an isopeptide bond in the core of the protein (many did go unnoticed) or an isoaspartate conversion happened or the loop is actually split due to an Asp-Pro pair autocatalysis, waters get mistaken for ligands or sidechains. If these are polished, they may get flagged as problematic, but they may be made worse.


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

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

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

(This comment got too long)

positions in a predicted protein chain

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

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

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

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

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

Please, add your last comment to your answer.

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

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

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


text:
How can I change all text in a .tsv file to 1, except the first three fields?
I want to convert a .tsv file table in binary format (1 and 0).  In which all cells containing text would become 1 and blank cells would become 0,  except first three columns.
The given .tsv file has 6 fields
head adrbvbeT_copy.tsv
chr1    1   100
chr1    101 200
chr1    201 300
chr1    301 400
chr1    401 500
chr1    501 600
chr1    601 700
chr1    701 800
chr1    801 900
chr1    901 1000

tail adrbvbeT.tsv
chrY    57226501    57226600
chrY    57226601    57226700    Mmusculus-jaspar2016-Tcf21-MA0832.1
chrY    57226701    57226800    Mmusculus-HOCOMOCOv10-TF7L2_MOUSE.H10MO.A
chrY    57226801    57226900    Mmusculus-HOCOMOCOv10-STAT1_MOUSE.H10MO.A
chrY    57226901    57227000    Mmusculus-HOCOMOCOv10-ZEP2_MOUSE.H10MO.D
chrY    57227001    57227100    NA-SwissRegulon-SMARCC2.SwissRegulon    Mmusculus-HOCOMOCOv10-PBX2_MOUSE.H10MO.C
chrY    57227101    57227200    Mmusculus-HOCOMOCOv10-ZIC3_MOUSE.H10MO.C    Mmusculus-HOCOMOCOv10-TBX5_MOUSE.H10MO.D    Mmusculus-HOCOMOCOv10-PIT1_MOUSE.H10MO.C
chrY    57227201    57227300    Mmusculus-HOCOMOCOv10-MAFB_MOUSE.H10MO.D    Mmusculus-UniPROBE-Tbp.UP00029  Mmusculus-HOCOMOCOv10-MAFB_MOUSE.H10MO.D
chrY    57227301    57227400    Mmusculus-HOCOMOCOv10-PITX2_MOUSE.H10MO.D   Hsapiens-jaspar2016-TGIF2-MA0797.1  Mmusculus-HOCOMOCOv10-PITX2_MOUSE.H10MO.D
chrY    57227401    57227415    Mmusculus-HOCOMOCOv10-TFDP1_MOUSE.H10MO.D   Mmusculus-HOCOMOCOv10-TEAD4_MOUSE.H10MO.D   Mmusculus-HOCOMOCOv10-ZBTB6_MOUSE.H10MO.D

First, I converted all blank cells to 0. Now I'm wondering how I can replace the text in the cells with 1.
awk 'BEGIN { FS = OFS = "\t" } { for(i=1; i<=NF; i++) if($i ~ /^ *$/) $i = 0 }; 1' adrbvbeT_copy.tsv > adrbvbeT_copyN0.tsv

head adrbvbeT_copyN0.tsv
chr1    1   100 0   0   0
chr1    101 200 0   0   0
chr1    201 300 0   0   0
chr1    301 400 0   0   0
chr1    401 500 0   0   0
chr1    501 600 0   0   0
chr1    601 700 0   0   0
chr1    701 800 0   0   0
chr1    801 900 0   0   0
chr1    901 1000    0   0   0

tail adrbvbeT_copyN0.tsv
chrY    57226501    57226600    0   0   0
chrY    57226601    57226700    Mmusculus-jaspar2016-Tcf21-MA0832.1 0   0
chrY    57226701    57226800    Mmusculus-HOCOMOCOv10-TF7L2_MOUSE.H10MO.A   0   0
chrY    57226801    57226900    Mmusculus-HOCOMOCOv10-STAT1_MOUSE.H10MO.A   0   0
chrY    57226901    57227000    Mmusculus-HOCOMOCOv10-ZEP2_MOUSE.H10MO.D    0   0
chrY    57227001    57227100    NA-SwissRegulon-SMARCC2.SwissRegulon    Mmusculus-HOCOMOCOv10-PBX2_MOUSE.H10MO.C    0
chrY    57227101    57227200    Mmusculus-HOCOMOCOv10-ZIC3_MOUSE.H10MO.C    Mmusculus-HOCOMOCOv10-TBX5_MOUSE.H10MO.D    Mmusculus-HOCOMOCOv10-PIT1_MOUSE.H10MO.C
chrY    57227201    57227300    Mmusculus-HOCOMOCOv10-MAFB_MOUSE.H10MO.D    Mmusculus-UniPROBE-Tbp.UP00029  Mmusculus-HOCOMOCOv10-MAFB_MOUSE.H10MO.D
chrY    57227301    57227400    Mmusculus-HOCOMOCOv10-PITX2_MOUSE.H10MO.D   Hsapiens-jaspar2016-TGIF2-MA0797.1  Mmusculus-HOCOMOCOv10-PITX2_MOUSE.H10MO.D
chrY    57227401    57227415    Mmusculus-HOCOMOCOv10-TFDP1_MOUSE.H10MO.D   Mmusculus-HOCOMOCOv10-TEAD4_MOUSE.H10MO.D   Mmusculus-HOCOMOCOv10-ZBTB6_MOUSE.H10MO.D

You can use if (!$i) -- an empty string is false, and 0 is false. And if the field has the value zero, we're setting it to zero anyway.

Why do you need to do this? It seems like an XY problem

Not an expert, but :
awk 'BEGIN { FS = OFS = "\t" } { for(i=4; i<=6; i++) if($i ~ /^ *$/) {$i = 0} else {$i = 1} }; 1'

Since you have a very specific file format, this code starts from field 4, checking whether the fields after it are empty or not. I assume also that know that you file has 6 fields max exactly, so I "hardcoded" this knowledge.

Great it's working... Thanks a lot, @Sam.


text:
How to get coefficient tables from multiple regression model result!
I've done multiple regression modeling in R using the below script. Now I'm wondering how can I make a table for the coefficients of each model of all models from lm, glm, and rf separately.

all_variables <- names(PaC09c_w20_dat3)
response_variables <- all_variables[c(1:3)]
predictors <- all_variables[-c(1:3)]

lm_model <- lapply(
  response_variables,
  function(x) lm(reformulate(termlabels = predictors, response = x), data = PaC09c_w20_dat3)
) |>
  setNames(response_variables)

glm_model <- lapply(
  response_variables,
  function(x) glm(reformulate(termlabels = predictors, response = x), data = PaC09c_w20_dat3)
) |>
  setNames(response_variables)

rf_ranger <- lapply(
  response_variables,
  function(x) ranger(reformulate(termlabels = predictors, response = x), data = PaC09c_w20_dat3, importance="impurity")
) |>
  setNames(response_variables)

save.image(file="PaC09c_w20_Reg.RData")

Input data looks like (Original Table is very large),
  head(PaC09c_w20_dat3)
        EE87865ln1   EE87866ln1   EE87895ln1  blood_vessel_w20 adrenal_gland_w20 bone_element_w20 bronchus_w20
               <dbl>       <dbl>      <dbl>     <dbl>             <dbl>            <dbl>        <dbl>
        1    0.00391     0.00326    0.00332       0                 0                1            0
        2    0.00139     0.00116    0.00132       0                 0                0            0
        3    0.00360     0.00270    0.00469       1                 1                0            1
        4    0.00323     0.00348    0.00339       0                 0                1            0
        5    0.00323     0.00330    0.00382       0                 1                0            0
        6    0.00278     0.00208    0.00214       0                 0                1            0

The expected output would look like this,
                           EE87865_lmcof   EE87866_lmcof
blood_vessel_w20           6.332633e-06    5.916615e-06
adrenal_gland_w20         -2.159470e-06   -2.399176e-06
bone_element_w20           8.578804e-06    8.106153e-06
bronchus_w20              -2.956291e-05   -2.152231e-05

Your code is neither reproducible (one would need the bed file) nor small enough to mentally compute what could be in each variable. How do you expect people on the forum to help you get to a random goal you have?

Well, I found a solution,
 coefs_lm<-as.data.frame(lapply(lm_model,function(x)coef(x)[2:5]))

 coefs_glm<-as.data.frame(lapply(glm_model,function(x)coef(x)[2:5]))

var_imp <- function(object, ...) {
  var_imp_r <- object$variable.importance
  var_imp_r
}

coefs_rf<-as.data.frame(lapply(rf_ranger,function(x)var_imp(x)[1:4]))

2:5 for lm and glm is the number of samples EE.., where 1 is intercept.

What is the 2:28 part? Is it specific to your dataset(s)?

Sorry, it should be 2:5 for this toy dataset in this example.

To get the coefficients alone
lm_model$coefficients
glm_model$coefficients

Everything including coefficients its,
summary(lm_model)
summary(glm_model)

I assume this is the same for summary(rf_ranger), albeit not sure whether coefficients would work here

lm_model$coefficients and glm_model$coefficients produces NULL. As you can see every model object (lm_model, glm_model, and rf_ranger) consists of multiple regression models of samples (e.g., EE...).

Given how OP has made these things lists, we'll probably need an lapply there.


text:
How to convert bed to wig file!
I was wondering how can I convert bed file to wig file.
I have several bed files which look like
head EE88324.shortlong.frag.bed
chr1    10003   10117   1
chr1    10003   10138   1
chr1    10003   10146   1
chr1    10008   10152   1
chr1    10009   10129   1
chr1    10009   10130   1
chr1    10025   10169   1
chr1    10026   10145   1
chr1    10030   10172   1
chr1    10052   10245   1

Now I want convert this bed file to wig file, which should look like this, (this toy bed and wig samples are different)
head -n 50 S9U_ZKDN240002406-1A_H5VMFDSXC_L1_bwaM_readCounts.w1000k.wig
fixedStep chrom=chr1 start=1 step=1000000 span=1000000
2682
8134
7087
6923
5585
5579
7056
5953
6174
7058
6606
6682
5630
4705
5432
6755
9442
6563
5891
6427
6033
6871
6360
6574

What I tried so far,
#bed to wig conversion
chromsizes="~/mapF/GRCh38_chr.txt"
bedGraphToBigWig="~/KentUtils/exe/linux.x86_64/bedGraphToBigWig"
bigWigToWig="~/KentUtils/exe/linux.x86_64/bigWigToWig"

conda activate bedops_env
    bedops --partition EE88324.shortlong.frag.bed |
    bedmap --echo --mean --delim '\t' --window 1000000 --fixed 1000000 - EE88324.shortlong.frag.bed |
    $bedGraphToBigWig - $chromsizes |
    $bigWigToWig - EE88324.shortlong.diff.partition.map.1MB.wig

I do not know bedops, so I assume that this two-step pipe at the end of bedmap produces a regular bedGraph.
The problem is that both bedGraphToBigWig and bigWigToWig like any of the kentUtils tools does not read or write from/to stdin/stdout. Hence, write the bedGraph to disk, then convert to BigWig, then convert to Wig. Tedious, but that's what kentUtils support -- they do not allow piping.

Could you add the relevant command please? As it stands, this is pointing the OP in the right direction, but it isn't actually giving the answer.


text:
How to quantifiy of specific genes from shotgun metagenome?
I have googled a "lot", couldn't find any specific answer to the question. So, I am here seeking for your guidance. My question is similar to this. I have several metagenome (n=30). But for simplicity’s sake, let's assume only one metagenome. I wanted to quantify lasA genes in the metagenome. So far, what I have done (not sure if it's the right approach) is:

Download the lasA genes from uniport.
I then made the multiple sequence alignment (lasA.msa) of the lasA gene downloaded in the previous step.
The protein profile was generated hmmbuild res.hmm lasA.msa
Genome assembly (contigs) were done with megahit megahit -1 r1.fastq.gz -2 r2.fastq.gz --min-contig-len 1000 -m 0.85 -o ASSEMBLY/ -t 20
ORFs were predicted and extracted using getorf from the contigs (from step 4)
The profile (from step 3) was used to search the ORFs using hmmsearch
Several hits were obtained from hmmsearch indicating that the lasA gene is present in the metagenome

Now, my question is, how do I count the number of reads affiliated to the lasA gene? I feel that I need to map the reads to the assembly, but I don't know how to go about that.
NOTE: I have used MG-RAST, Humann3, and eggNOG-mapper, but none of them showed the presence of lasA gene, however the above-mentioned hmmsearch was able to find the gene, what could be the reason? Problem with my approach? Or hmmsearch using the protein profile is more sensitive?
Thanks and regards

The way I would do it is to make an assumption that the reads mapping to your contig == (abundance of the genes in the contig)/contig length.

Map raw reads to the contig containing your hit (original contig, pre-translation)

Normalise it with RPKM:

RPKM = NumberOfReadsThatAligned / ( ContigLength/1000 * totalNumberReads/1,000,000 )

This will give you the relative abundance of the contig which is an APPROXIMATION of the abundance of your gene.

Alternatively, to be more precise you can FragGeneScan your raw reads, extract the sequence of your putative protein and map translated raw reads to it.
For your other question: Yes, HMMs are more sensitive than blast or diamond since its output is probability score that these protein contain motifs extracted from your alignments (given that there is one). For that reason you need to look at the bitscore and E-value. If the bitscore is low and/or E-value is high, you have a high chance of false-positive.
Best wishes,
A


text:
How to copy only certain counts for genes in tsv file to new file in linux
Hi there I have generated a counts table of samples I need to compare by differential expression analysis.
The layout of the counts table is as follows:

Gene_id
Sample_A_r1
Sample_A_r2
Sample_B_r1
Sample_B_r2
Sample_AB_r1
Sample_AB_r2

KLTH001
0
0
2
5
8
3

KLTH002
0
0
5
3
9
7

YL99812
6
4
0
0
12
14

YL99813
9
7
0
0
15
19

So here, for the differential expression analysis, I only want to compare for example the genes starting with "YL-". So in a new file I want the rows including the first row which is "Gene_id" and the sample names e.g. "Sample_A_r1" and then only the rows with genes starting with "YL-" in their names. This will allow me to only compare the changes between Sample_A alone versus when it was combined with Sample_B, or in the table Sample_AB.
If anyone can please give me an indication of how to do this in linux I would really appreciate it. I have seen other similar posts, involving grep and awk I am just not sure of exactly how to do this. I am still quite new to linux. Any help will be appreciated.

Would something like this work?
awk ' $1 YL* ' counts_table.tsv > edited_counts_table.tsv

Would this preserve the header/first line?

Something like this should do job:
# creating a new file with the first line of the original
head -1 foo.txt > bar.txt

# ">>" appends to the file on the right hand side
# "^" dictates that the search term should be at the beginning of a line
grep "^YL" foo.txt >> bar.txt

Thanks very much, will give it a try now

Alternatively, you can use grep like so to remove all lines starting with KLTH:
grep -v '^KLTH' in_file > out_file

For example:
echo "Gene_id\nKLTH001\nYL99812" | grep -v '^KLTH'

Prints:
Gene_id
YL99812

Here, GNU grep uses option:
-v : Print lines that do not match.
In the regular expression (regex):
^ : Match the beginning of the line.

One way using awk:
awk 'NR==1 || /^YL/ { print $0 }' file.tsv

The default action is to print the line/record, so all you need is the condition:
awk 'NR==1 || /^YL/' file.tsv

Thanks I will try this out, it seems the previous command didn't copy all of the occurrences of "YL-" across. Could you maybe suggest how I could just remove any gene names starting with "KLTH-" from my original counts file as well? I would just like to be able to check that I got all records from the original file in the separated files I create.

@Justin1609 You can just negate the pattern match expression: awk 'NR==1 || !/^KLTH/' file. This will just print the first line, or lines that don't start with KLTH.

okay awesome thanks @Steve I appreciate it, will give it a try


text:
Slow bwa shm "preloading" into memory
Coming from devops background, working on automating usual genomics pipelines, I am benchmarking existing pipelines in order to choose cloud resources properly.
Past week most of work was on testing bwa mem.
I've noticed that bwa takes about 10 minutes to start processing data. When monitoring bwa resource usage I've noticed that RAM is slowly building up to about 6GB. That is amount of reference indexed data.
After bwa loads the reference first time, any subsequent run is fast. But why this preloading is happening so slow? My reference data set is already indexed with bwa index.
When I monitor disk activity with iotop I noticed that disk read is only 5MB/s, when I run bwa mem alignment (only first time), which is very slow. Performance of the drive is lot more than that and a simple copy test confirms that.
Is it a "limitation" of bwa, is this slow preloading expected to be like that and I should just accept it and move on forward or there is something I don't understand and maybe someone can help me understand what could I do?
The same "slow preload" happens when I run bwa shm /path/to/ref.
Then bwa mem starts immediately. But again slow preload is something that is bothering me and I am trying to understand what can I do about it.

Issue was with EBS snapshot which holds reference indexed data. Issue is that when you run EC2 from snapshot, EBS storage performance is expected to be low until you optimize the volume. My assumptions and questions for two weeks, all "goes to the water".
Unfortunately I was not aware of this hence I asked the question out from ignorance.

Can you please add this as an actual answer and accept it? This will help the StackExchange algorithms to identify this as a question with and answer, and stop promoting it in the future to fish for more answers.

I will. Didn not want to do that imediatelly since it is a little bit awkward to answer own question = )

Yeah, I find it awkward as well, but you did already answer the question (as a comment to your own question). Might as well get a bit more credit for that. Also, it's worth thinking about comments as transient things that are there only to help people improve their work; they're not meant to be used for anything substantial that should hang around for a long time (beyond gratitude).

I'll add an answer to my own question, which was asked out of lack of experience with cloud (AWS) EC2 service.
The way I was testing bwa is from a custom AMI image created for our pipelines.
The image has "preloaded" reference data that we are using for alignement so that we spare some time downloading reference data from S3 and spare some S3 api calls since we have a lot of files used as reference data.
Simply our image is single 200GB EBS volume used with AWS ECS AMI where we stored our reference data and created a custom image for Batch out of it. Then we use that AMI for AWS Batch jobs.
Once I used iotop (after few days of testing) to monitor disk activity I've noticed that bwa loads index into memory very slow, cca 5MB/s. Syncing files from S3 was faster then using this snapshot approach.
For benchmarking I've used the same AMI manually launched. In the AWS documentation there is part which explains why and how you should optimize EBS when AMI is launched from a snapshot.
https://docs.aws.amazon.com/AWSEC2/latest/UserGuide/ebs-initialize.html
After spining up fresh EC2, syncing reference data and running bwa everything was clear. My EBS volume was slow and had to be initialized first.
That is all. My lack of EBS storage knowledge was the issue.
This helped me understand that I should not use custom AMI with ref data stored and will probably use EFS or simply syncing the data from S3 inside Docker entrypoint.

Answer from @user438383, converted from comment:
BWA has slow indexing due to the indexing size; this can be improved by using bwa-mem2:
https://github.com/bwa-mem2/bwa-mem2

We are happy to announce that the index size on disk is down by 8 times and in memory by 4 times due to moving to only one type of FM-index (2bit.64 instead of 2bit.64 and 8bit.32) and 8x compression of suffix array. For example, for human genome, index size on disk is down to ~10GB from ~80GB and memory footprint is down to ~10GB from ~40GB. There is a substantial reduction in index IO time due to the reduction and hardly any performance impact on read mapping.

Yes but it is not yet production ready. Talking from my very poor amount of reading on the bwa-mem2 topics.


text:
DESEQ2 Question about results()
I am currently learning to perform Differential Analysis via DESEQ2 R Package, and I believe I've made progress, able to format the data correctly [maybe] for DDS(). When I run the results function to see the output, the data seems fine, as can be seen below:
log2 fold change (MLE): Group RA vs C  Wald test p-value: Group RA vs
C  DataFrame with 6 rows and 6 columns
       baseMean log2FoldChange     lfcSE      stat    pvalue      padj
      <numeric>      <numeric> <numeric> <numeric> <numeric> <numeric> Gnai3   95.1978      0.3458366  0.275366  1.255916  0.209146  0.999383
Cdc45   21.9942     -0.0545177  0.346199 -0.157475  0.874870  0.999383
H19     29.0133     -1.4954220  1.516390 -0.986173  0.324048  0.999383
Narf    16.6967      0.2571059  0.370880  0.693232  0.488164  0.999383
Cav2    63.6356      0.5060019  0.354969  1.425480  0.154018  0.999383
Klf6    65.8469      0.2851175  0.265755  1.072858  0.283335  0.999383

However, when I run summary() or any other function, they do not seem to behave correctly, not counting any stats. Same with MAplot not coloring genes that should be marked red:
out of 13184 with nonzero total read count adjusted p-value < 0.1 LFC
> 0 (up)       : 0, 0% LFC < 0 (down)     : 0, 0% outliers [1]       : 150, 1.1% low counts [2]     : 0, 0% (mean count < 1)

I'm confused as to why this is occurring and was wondering if anyone had ideas.
EDIT: Added More Code for context
sergey <- (read.csv("cleanmatrix25+.csv", sep = "\t", header = T,
    check.names = F, row.names = "Gene_names"))

meta <-(read.csv("metadata25+.csv", sep = "\t", row.names = 1))

dds <- DESeqDataSetFromMatrix(countData = sergey, colData = meta,
    design = ~Group)

keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]

dds$Group <- relevel(dds$Group, ref = "C")

dds <- DESeq(dds)

res <- results(dds)

head(res)

summary(res)

plotMA(res, ylim=c(-2,2))

EDIT2:
I sorted result by p-value, so all of these should be red on the following MAplot due to having a p-value < 0.1


text:
Translate all reads in a .fastq into protein sequence from deep mutational scanning experiment
I'm working with paired-end NGS reads from an Illumina platform. The sequences I have are all of the same gene, but have one or more substitution mutations each.
Here is a rough workflow for generating the reads:

Expression of a library of sequence variants for a heterogeneous protein of interest in phage.
Selection of high-fitness variants according to a custom assay.
Extraction of the variant pool via amplification of the variable region of the gene of interest.
Preparation of Illumina sequencing library from the amplicon pool.
Deep sequencing of the amplicon pool.

I've already combined the paired-end reads with USearch and manually removed sequences with subpar quality scores, but I now need to translate all of the remaining DNA sequences into amino acid sequences. I'm uncertain of where to begin addressing this issue, and I was hoping someone could point me in the direction of an application or python package that could assist me.
My goal: Count the abundances of each sequence variant within the Illumina data to identify which of the protein sequences are the most fit.

Won't you need to map to a reference? You will need to know what frame to use to translate the reads. Unclear to me what the experimental design is here. Where do these reads come from? Is it an amplicon library? (If it is amplicon it is possible that you don't need to map anywhere.)

Reads are coming from a mixed population of phages which express multiple sequence variants of a protein of interest on their surface. The protein of interest is then put under selection. after selection the gene sequences of phages were collected and quantified to identify copy number of each variant in the population. This allows us to determine the fitness of each variant by seeing how many of that sequence there is.
I suppose I might need a reference map, but I am pretty inexperienced to know which program would be best.

"gene sequences of phages were collected" how? Presumably by sequencing, but what was the sequencing library preparation procedure? Was it whole genome shotgun e.g. nextera prep, locus-specific amplification (i.e. your protein of interest), barcode amplicon, something else? Presumably you limited sequencing to your protein of interest somehow, otherwise this would be a somewhat inefficient workflow.

Thank you for these questions! This is not my wet experiment, I am simply trying to use the data. But I can tell you the read libraries were generated with locus specific PCR.

Ok, one further question: Do all of your forward reads start at the same coordinate? E.g. are most reads mostly the same sequence? If paired-end, same question for reverse reads.

I already accepted your answer, but yes this is the case! Thank you very much!

I might edit your question slightly to include additional relevant info.

Based on back and forth in the comments, I think that you might be interested in a deep mutational scanning workflow, for example dms_tools2 or enrich2.
These are built specifically for the kind of experiment that you are interested in, if I understand correctly. See attached figure from dms_tools2:

Thank you! I think I will be using enrich2 for this data processing, it matches very well with my needs!


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

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

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

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

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

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

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

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

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


text:
Funcotator reference file error in GATK4
I am using a function called Funcotator (here) in GATK4. The code I am trying to run is:
gatk Funcotator \
--variant chr21.vcf \
--reference hg38.fa.fai \
--ref-version hg38 \
--data-sources-path funcotator_dataSources.v1.2.20180329 \
--output chr21.funcotated.maf \
--output-file-format MAF

However, it gives me an error message like this:
java.lang.IllegalArgumentException: File is not a supported reference file type: /Users/data/Bam/hg38.fa.fai
at htsjdk.samtools.reference.ReferenceSequenceFileFactory.lambda$getFastaExtension$0(ReferenceSequenceFileFactory.java:253)
at java.util.Optional.orElseGet(Optional.java:267)
at htsjdk.samtools.reference.ReferenceSequenceFileFactory.getFastaExtension(ReferenceSequenceFileFactory.java:253)
at htsjdk.samtools.reference.ReferenceSequenceFileFactory.getDefaultDictionaryForReferenceSequence(ReferenceSequenceFileFactory.java:223)
at org.broadinstitute.hellbender.utils.fasta.CachingIndexedFastaSequenceFile.checkFastaPath(CachingIndexedFastaSequenceFile.java:180)
at org.broadinstitute.hellbender.utils.fasta.CachingIndexedFastaSequenceFile.<init>(CachingIndexedFastaSequenceFile.java:143)
at org.broadinstitute.hellbender.utils.fasta.CachingIndexedFastaSequenceFile.<init>(CachingIndexedFastaSequenceFile.java:125)
at org.broadinstitute.hellbender.utils.fasta.CachingIndexedFastaSequenceFile.<init>(CachingIndexedFastaSequenceFile.java:110)
at org.broadinstitute.hellbender.engine.ReferenceFileSource.<init>(ReferenceFileSource.java:35)
at org.broadinstitute.hellbender.engine.ReferenceDataSource.of(ReferenceDataSource.java:27)
at org.broadinstitute.hellbender.engine.GATKTool.initializeReference(GATKTool.java:439)
at org.broadinstitute.hellbender.engine.GATKTool.onStartup(GATKTool.java:722)
at org.broadinstitute.hellbender.engine.VariantWalker.onStartup(VariantWalker.java:45)
at org.broadinstitute.hellbender.cmdline.CommandLineProgram.runTool(CommandLineProgram.java:138)
at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMainPostParseArgs(CommandLineProgram.java:192)
at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:211)
at org.broadinstitute.hellbender.Main.runCommandLineProgram(Main.java:160)
at org.broadinstitute.hellbender.Main.mainEntry(Main.java:203)
at org.broadinstitute.hellbender.Main.main(Main.java:289)

Does anyone know how to solve this problem? I tried to use the reference file downloaded from the GATK cloud, but it still doesn't work and has the same error message.
Also, I wonder if the reference files must be the version provided by the GATK team? The one I am using now is not from the GATK team, and as they force me to provide --data-sources-path, I just copied the one they used on their post. Could that cause any problem in the further analysis?
Thanks a lot ~~~
Here is an update:
I first tried hg38.fa as my reference and got this error:
A USER ERROR has occurred: Fasta dict file file:///Users/data/Bam/hg38.dict for reference file:///Users/data/Bam/hg38.fa does not exist. Please see http://gatkforums.broadinstitute.org/discussion/1601/how-can-i-prepare-a-fasta-file-to-use-as-reference for help creating it.

So I generated the hg38.fa.fai, then I got the error message shown above.
Then, I tried the reference file downloaded from the GATK cloud, resources-broad-hg38-v0-Homo_sapiens_assembly38.fasta.fai, but the same error appeared again.
Another update:
I tried to use the hg38.dict as the reference and used the same code:
    gatk Funcotator \
    --variant chr21.vcf \
    --reference hg38.dict \
    --ref-version hg38 \
    --data-sources-path funcotator_dataSources.v1.2.20180329 \
    --output chr21.funcotated.maf \
    --output-file-format MAF

However, I got the same error message again:
java.lang.IllegalArgumentException: File is not a supported reference file type: /Users/data/hg38.dict

Are all the files .fai, .dict and .fa in the same directory and have the same names?

Yeah, I think so.

You need to pass the fasta file instead of fa.fai or dict file as reference
    gatk Funcotator \
    --variant chr21.vcf \
    --reference hg38.fa \
    --ref-version hg38 \
    --data-sources-path funcotator_dataSources.v1.2.20180329 \
    --output chr21.funcotated.maf \
    --output-file-format MAF

The directory where you have hg38.fa must contain also the hg38.fa.fai and hg38.dict
Hope this helps
Zhan