text:
How can I programmatically add a Hydrogen 'Atom' to a 'Residue' object?
I know the algorithm for creating a Hydrogen atom and adding to a residue:
Point3d  create_hydrogen(Point3d C, Point3d N, Point3d CA, Point3d H)
{
  H.set(N);
  H -= C;
  H.norm();

  Point3d tmp2(N);
  tmp2 -= CA;
  tmp2.norm();

  H += tmp2;
  H.norm(N_H_BOND_LENGTH);
  H += N;

  return H;
}

Point3d  create_hydrogen(Residue prev_residue, Residue residue, Point3d H)
{
  return create_hydrogen(prev_residue.find_atom(" C  "),
                              residue.find_atom(" N  "),
                              residue.find_atom(" CA "),
                                   H);
}

Atom create_hydrogen(Residue prev_residue, Residue residue)
{
  Atom H = make_shared_atom(0, " H  ", 1);
  create_hydrogen(prev_residue, residue, H);
  residue.add(H);
  return H;
}

I didn't find this implementation in BioPython.
Is this already implemented in BioPython? If YES, in which module?
If NO, what would be the implementation like in BioPython?
Added image of peptide bond too try to depict the problem at hand, I believe the hydrogen in cause is the one on the peptide bond plane(pale yellow) indentifiedin purple

Search for Bio.PDB.Structure.Structure object in biopython code

I would suggest again using already existing tools to add hydrogens. You said your PI had concerns with this, but I must say that it is not a random process like letting kids adding baubles to a Christmas tree: for sp3 heavy atoms bound by two heavy atoms and for terminal heavy atom in sp2 the hydrogens have a defined geometric position. Terminal amines, hydroxyls and methyls will require more advanced calculations, but bar for ST turns these will not contribute to SS.

@MatteoFerla, I understand. I talked to him, and he doesn't agree. This is his project and his decision is final. :) What can I do?

@pippo1980, Search for Bio.PDB.Structure.Structure object in biopython code --- there is nothing useful in that module.

Switch to openbabel ? Look how the Arpeggio guys handled the visualization of ligand protein interaction of pdb : https://github.com/harryjubb/arpeggio/blob/master/arpeggio.py

@pippo1980, How about Biotite?

Or http://thegrantlab.org/bio3d/

@pippo1980, Bio3D is offered in R. That's not useful .

You could refactor it in Python ;-)

Meaning translate/port

@pippo1980, That will take months. I don't have that much time to invest right now.

Is it stille related to Calculate Hydrogen bonds between each two chains Biopython : like https://stackoverflow.com/questions/43030557/calculate-hydrogen-bonds-between-each-two-chains-biopython

The code that you posted is a trivial function that adds H atom to N atom in a given distance. Do you want to reproduce just this?

@marcin, For me it is non-trivial. I am finding to similarity in BioPython and this code.

@marcin, Do you want to reproduce just this? --- For now, YES.

First you need to understand how this code work. It's a geometrical exercise. C, N, CA and H are in the same plane. The code seems to assume that angles H-N-C and H-N-CA are equal (which is approximately true). It calculates two vectors, normalizes them, adds them and sets the length of the resulting vector. So you'd need to do the same in Python. I don't know why you want to use BioPython (and what's the point of the whole exercise), but you can implement the same in Python using, for example, Vector from Bio.PDB and then add H atom to Bio.PDB.Residue.

@marcin, what framework do u recommend other than BioPython?

Hard to tell without knowing more about this project. I see that many people use PyMOL for simple scripting (it's molecular viewer, but supports scripts in Python). If you'd like to add all hydrogen atoms, the most popular program is REDUCE from the Richardson lab.

@marcin, I need a python-based Bioinformatics protein manipulation framework that is Open source and easy to learn. E.g. SeqAn has great documentation. They also published a book which is very convenient for me.

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2724740/ HAAD: A Quick Algorithm for Accurate Prediction of Hydrogen Atoms in Protein Structures it compares and cites other alghoritm too, reduce seems to be the older one maybe still the more robust accurate

Not sure but here https://biopython.org/docs/1.75/api/Bio.PDB.Residue.html   Bio.PDB.Residue module Your can find :
add(self, atom) :
      Add an Atom object.

Checks for adding duplicate atoms, and raises a PDBConstructionException if so.
Can’t try it just now but who knows it may work
Googling around I found out that if you plan to add Hydrogens to PDB files [I believe nowadays PDBx is the repository grade standard] you could do thar following PDB format rules. I am copying them from UCSF Chimera website (https://www.cgl.ucsf.edu/chimera/docs/UsersGuide/tutorials/pdbintro.html#hydrogens):

Hydrogen Atoms

In brief, conventions for hydrogen atoms in version 3.0 PDB format are as follows:

Hydrogen atom records follow the records of all other atoms of a particular residue.
A hydrogen atom name starts with H. The next part of the name is based on the name of the connected nonhydrogen atom. For example, in amino acid residues, H is followed by the remoteness indicator (if any) of the connected atom, followed by the branch indicator (if any) of the connected atom; if more than one hydrogen is connected to the same atom, an additional digit is appended so that each hydrogen atom will have a unique name. Hydrogen atoms in standard nucleotides and amino acids (other than the rarely seen HXT) are named according to the IUPAC recommendations (Markley et al., Pure Appl Chem 70:117 (1998)). Names of hydrogen atoms in HETATM residues are determined in a similar fashion.
If the name of a hydrogen has four characters, it is left-justified starting in column 13; if it has fewer than four characters, it is left-justified starting in column 14.

In the following excerpt from entry 1vm3, atom H is attached to atom N. Atom HA is attached to atom CA; the remoteness indicator A is the same for these atoms. Two hydrogen atoms are connected to CB, one is connected to CG, three are connected to CD1, and three are connected to CD2.

ATOM     10  N   LEU A   2       4.595   6.365   3.756  1.00  0.00           N
ATOM     11  CA  LEU A   2       4.471   5.443   2.633  1.00  0.00           C
ATOM     12  C   LEU A   2       5.841   5.176   2.015  1.00  0.00           C
ATOM     13  O   LEU A   2       6.205   4.029   1.755  1.00  0.00           O
ATOM     14  CB  LEU A   2       3.526   6.037   1.578  1.00  0.00           C
ATOM     15  CG  LEU A   2       2.790   4.919   0.823  1.00  0.00           C
ATOM     16  CD1 LEU A   2       3.803   3.916   0.262  1.00  0.00           C
ATOM     17  CD2 LEU A   2       1.817   4.196   1.769  1.00  0.00           C
ATOM     18  H   LEU A   2       4.169   7.246   3.704  1.00  0.00           H
ATOM     19  HA  LEU A   2       4.063   4.514   2.992  1.00  0.00           H
ATOM     20  HB2 LEU A   2       2.804   6.675   2.065  1.00  0.00           H
ATOM     21  HB3 LEU A   2       4.099   6.623   0.873  1.00  0.00           H
ATOM     22  HG  LEU A   2       2.234   5.353   0.004  1.00  0.00           H
ATOM     23 HD11 LEU A   2       4.648   4.447  -0.148  1.00  0.00           H
ATOM     24 HD12 LEU A   2       3.334   3.331  -0.516  1.00  0.00           H
ATOM     25 HD13 LEU A   2       4.137   3.260   1.052  1.00  0.00           H
ATOM     26 HD21 LEU A   2       0.941   3.892   1.216  1.00  0.00           H
ATOM     27 HD22 LEU A   2       1.522   4.860   2.568  1.00  0.00           H
ATOM     28 HD23 LEU A   2       2.296   3.323   2.188  1.00  0.00           H

```

As you can see from the algorithm, adding Hydrogen is not that straightforward. :D If it were that simple, I wouldn't have posted this question!!

Are you talking about the geometry of the missing H or about the protonation state of a lateral chain or about polar hyrogens ?

see here for H atom names for Phenilalanine , bond lenght and angle have to be calculated from the atom coords https://www.ebi.ac.uk/pdbe-srv/pdbechem/atom/list/PHE or tabulated here https://www.ebi.ac.uk/pdbe-srv/pdbechem/bond/list/PHE

The other answer is correct. But I thought I'd give a few pointers in how one can understand how to do something with a given Python module —teach a man to fish kind of thing...
In a Jupyter Notebook or IPython shell (don't), you can do:

help(function) print out the docstring
dir(object) shows all the public and magic attributes and methods —note a class will not show the instance methods
type(object) or object.__class__ gives you what the class is, but object.__class__.__mro__ tells you what it inherits.
in the standard library module inspect is the function inspect.getsource(function), which for a notebook is messy due to newlines so needs to  be print_code = lambda fun: print(inspect.getsource(fun)) and print_code(fun).

In your case you read a file...
from Bio import PDB
parser = PDB.PDBParser()
structure : PDB.Structure.Structure = parser.get_structure("tryp-cage", "1L2Y.first.pdb")

So parser.get_structure "magically" makes a structure, eh? This means it can add stuff! As mentioned print_code(parser.get_structure) will reveal its secrets, or within it call some self._private_fun(), which can be inspected (self = parser) all the way down the rabbit hole. Until you find a nice constructor that declares nicely an atom instance or whatever you ever fancy searching for —that is the great thing of open source code.
In the specific case of PDB module, everything inherits Entity, which has the attribute child_list, which is what the iterators (__iter__ is the topmost) read without making a copy. The end of the above rabbit hole may just add a child Entity-inheriting class instance to this list, but it also do other things that is worth keeping an eye out.
In the case of residue.add(atom), creating the Bio.PDB.Atom.Atom instance is slightly problematic as there's no typehinting or :type xxx: declarations in the docstring. But doing type(atom.bfactor) etc. on the attributes of an existing atom will reveal how to do it. One the serial number is a PDB thing, which is sequential for the model.
last_serial = list(structure.get_residues())[-1].child_list[-1].serial_number + 1
residue.add(PDB.Atom.Atom(name=' H  ', coord=np.array([1,1,1]), bfactor=0., occupancy=1., altloc=' ', fullname=' H  ', serial_number=last_serial+1,element='H'))

The coordinates will not be [1,1,1], but the point that you have calculate. Say for H, this would be coplanar with C, CA and N, that is 1 and smidge Å away from N at an angle not quite in the middle as found in amide bonds.

If the other answer is correct, what do you think about the algorithm I posted? This algorithm is taken from the c++ source code I have in my hand.

As you would realize, I can write Python and C++ program, but I am just confused about hydrogens. :)

Biopython documentation is not that rich especially for ones like me that use google as primary/unique docs search

I have expanded the answer with an example, showing that what I wrote is how one finds the correct solution without documentation. In terms of how to do the maths, you have to do that in numpy as its simple maths —I do not know if the C++ code is 100% correct as it is projecting on the line between C-CA midpoint and N, which I am not sure is the correct angle. Totally easy to verify though

I would be curious about how Bio.PDB.PDBIO module (https://biopython.org/docs/1.74/api/Bio.PDB.PDBIO.html) would save such pdb with added H atoms, like : https://www.cgl.ucsf.edu/chimera/docs/UsersGuide/tutorials/pdbintro.html#hydrogens or something else ?

@MatteoFerla Hi, sorry to bother, would you be so kind to point me to a PDB (x-ray) into https://www.rcsb.org/ , that contains H atoms. Don't know how to search for them, or some other PDB compliant file with Hs, to play with biopython I am curious about it, I know I should learn PDBx by now but is far more complicated to me

dont worry got it from http://biophysics.cs.vt.edu/ even if it marks all H as atom 0, without renumbering everything from scratch

Neutron diffraction structures have protons experimentally. Say PDB:5MOP —A beauty as it has a hydronium or two. For small molecules x-ray crystallography can detect protons but the database for that is the CCDC

Do you have code for the correct naming of H atoms ? Does PDBx follow same rules ?

All formats that support atom names use the same atom names and should match the definition in the PBD. e.g. https://files.rcsb.org/ligands/view/ALA.cif Ligands are not always numbered by canonical index order, but amino acids follows a precise lettering scheme (greek) for heavy atoms —the annoying Val <-> Ile superposition issue is common to all tools. However, for protons there are two schemes:  HB1 vs. 1HB , the former is like the heavy atoms while the latter is more sensible as the number prefix is compatible with heavy atoms who have a number suffix, cf. ARG with HH11 vs. 1HH1

—the annoying Val <-> Ile superposition issue is common to all tools ... what is that ?

The CB atom is not chiral in valine, but the atom names make it "name chiral" (not a real word, but just like isotopic chirality): looking at valine with the HB atom forewards, you have CG2 on the side of the N (left) and CG1 on the C (right) in a predictable counterclockwise fashion. Isoleucine is chiral on the CB atom and the CG1 is on the left as it is more complex due to CD1. As a result they do not superpose well by atom names, but do if CG2 is mapped to CG1 and CG1: CG2.

@MatteoFerla  AttributeError: 'Structure' object has no attribute 'get_residue' --> last_serial = list(structure.get_atoms())[-1].serial_number ?

OK my attempt, input pdb created in pymol with fab AVGLWPT , AVGLWPT , hydro=0 and saved as avgwpt.pdb :
ATOM      1  N   ALA     1      67.846  51.727  15.161  1.00  0.00           N
ATOM      2  CA  ALA     1      68.522  53.021  15.161  1.00  0.00           C
ATOM      3  C   ALA     1      70.022  52.847  15.161  1.00  0.00           C
ATOM      4  O   ALA     1      70.553  51.730  15.150  1.00  0.00           O
ATOM      5  CB  ALA     1      68.014  53.813  16.378  1.00  0.00           C
ATOM      6  N   VAL     2      70.773  53.896  15.173  1.00  0.00           N
ATOM      7  CA  VAL     2      72.223  53.728  15.173  1.00  0.00           C
ATOM      8  C   VAL     2      72.922  55.066  15.188  1.00  0.00           C
ATOM      9  O   VAL     2      72.296  56.134  15.198  1.00  0.00           O
ATOM     10  CB  VAL     2      72.665  52.877  13.916  1.00  0.00           C
ATOM     11  CG1 VAL     2      74.193  52.762  13.664  1.00  0.00           C
ATOM     12  CG2 VAL     2      72.143  51.420  13.932  1.00  0.00           C
ATOM     13  N   GLY     3      74.212  55.098  15.190  1.00  0.00           N
ATOM     14  CA  GLY     3      74.888  56.392  15.204  1.00  0.00           C
ATOM     15  C   GLY     3      76.388  56.218  15.203  1.00  0.00           C
ATOM     16  O   GLY     3      76.917  55.104  15.192  1.00  0.00           O
ATOM     17  N   LEU     4      77.138  57.267  15.215  1.00  0.00           N
ATOM     18  CA  LEU     4      78.589  57.099  15.215  1.00  0.00           C
ATOM     19  C   LEU     4      79.288  58.437  15.230  1.00  0.00           C
ATOM     20  O   LEU     4      78.660  59.507  15.241  1.00  0.00           O
ATOM     21  CB  LEU     4      79.040  56.251  13.992  1.00  0.00           C
ATOM     22  CG  LEU     4      78.445  54.825  13.832  1.00  0.00           C
ATOM     23  CD1 LEU     4      79.094  54.113  12.637  1.00  0.00           C
ATOM     24  CD2 LEU     4      78.610  53.968  15.098  1.00  0.00           C
ATOM     25  N   TRP     5      80.577  58.469  15.232  1.00  0.00           N
ATOM     26  CA  TRP     5      81.254  59.763  15.246  1.00  0.00           C
ATOM     27  C   TRP     5      82.754  59.589  15.245  1.00  0.00           C
ATOM     28  O   TRP     5      83.281  58.478  15.233  1.00  0.00           O
ATOM     29  CB  TRP     5      80.817  60.551  16.512  1.00  0.00           C
ATOM     30  CG  TRP     5      79.341  60.958  16.536  1.00  0.00           C
ATOM     31  CD1 TRP     5      78.304  60.258  17.187  1.00  0.00           C
ATOM     32  CD2 TRP     5      78.731  62.003  15.872  1.00  0.00           C
ATOM     33  CE2 TRP     5      77.340  61.929  16.130  1.00  0.00           C
ATOM     34  CE3 TRP     5      79.257  63.017  15.031  1.00  0.00           C
ATOM     35  NE1 TRP     5      77.049  60.848  16.950  1.00  0.00           N
ATOM     36  CZ2 TRP     5      76.464  62.876  15.556  1.00  0.00           C
ATOM     37  CZ3 TRP     5      78.372  63.946  14.482  1.00  0.00           C
ATOM     38  CH2 TRP     5      76.996  63.879  14.742  1.00  0.00           C
ATOM     39  N   PRO     6      83.505  60.638  15.257  1.00  0.00           N
ATOM     40  CA  PRO     6      84.959  60.560  15.257  1.00  0.00           C
ATOM     41  C   PRO     6      85.498  60.036  16.599  1.00  0.00           C
ATOM     42  O   PRO     6      84.724  59.735  17.507  1.00  0.00           O
ATOM     43  CB  PRO     6      85.354  62.018  14.978  1.00  0.00           C
ATOM     44  CG  PRO     6      84.201  62.829  15.588  1.00  0.00           C
ATOM     45  CD  PRO     6      82.967  61.980  15.271  1.00  0.00           C
ATOM     46  N   THR     7      86.770  59.913  16.774  1.00  0.00           N
ATOM     47  CA  THR     7      87.263  59.414  18.054  1.00  0.00           C
ATOM     48  C   THR     7      88.771  59.342  18.062  1.00  0.00           C
ATOM     49  O   THR     7      89.436  59.668  17.100  1.00  0.00           O
ATOM     50  CB  THR     7      86.639  58.014  18.376  1.00  0.00           C
ATOM     51  CG2 THR     7      87.080  57.358  19.702  1.00  0.00           C
ATOM     52  OG1 THR     7      85.224  58.110  18.479  1.00  0.00           O
TER
END

pic :

code :
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
molecule create in pymol with fab AVGLWPT , AVGLWPT , hydro=0

"""

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

from Bio.PDB import (
                    PDBParser,
                    Atom,
                    PDBIO
                    )

structure_code = 'avgwpt'

structure_name = structure_code

parser_pdb= PDBParser(QUIET=False)    ### shows warnings

structure_whole = parser_pdb.get_structure(structure_name, structure_name+".pdb")

last_serial = list(structure_whole.get_atoms())[-1].serial_number

print(" last_serial : ",last_serial)

for model in structure_whole:

    for chain in model:

        for resi in chain:

            print('1 : ', resi)

            print('2 : ',resi.id)

            try :

                print('3 : ',chain[resi.id[1]-1])

                doable = True

            except :

                print('3 : ','residue -1 missing !!!!!!!!!!!!!!!!!!!!!!1')

                doable = False

            print('4 : ',[i for i in resi.get_atoms()])

            if doable:

                print('5 : ','ok')

            else:

                print('5 : ','NO')

            if doable:

                print(chain[resi.id[1]-1]['C'])

                print(chain[resi.id[1]-1]['C'].get_vector(), type(chain[resi.id[1]-1]['C'].get_vector()))

                print(chain[resi.id[1]]['N'].get_vector())

                c = chain[resi.id[1]-1]['C'].get_vector()

                n = chain[resi.id[1]]['N'].get_vector()

                ca = chain[resi.id[1]]['CA'].get_vector()

                nc = n-c

                nca = n-ca

                print('vector BN : ', nc)

                nc.normalize()

                print('vector BN normalize : ', nc)

                print('vector B2B : ', nca)

                nca.normalize()

                print('vector B2B normalize : ', nca)

                print('sum of nc + nca : ' , nsum := nc + nca)

                nsum.normalize()

                nsum**1.023

                nsum += n

                print(nsum)

                print(nsum.get_array(), type(nsum.get_array()))

                resi.add(Atom.Atom(name=' H  ', coord = nsum.get_array(), bfactor=0., occupancy=1., altloc=' ', fullname=' H  ', serial_number = last_serial ,element='H'))

                last_serial += 1

io=PDBIO()

io.set_structure(structure_whole )

io.save(structure_name+"_H_added.pdb" , preserve_atom_numbering = True)

output , saved as avgwpt_H_added.pdb:
ATOM      1  N   ALA     1      67.846  51.727  15.161  1.00  0.00           N
ATOM      2  CA  ALA     1      68.522  53.021  15.161  1.00  0.00           C
ATOM      3  C   ALA     1      70.022  52.847  15.161  1.00  0.00           C
ATOM      4  O   ALA     1      70.553  51.730  15.150  1.00  0.00           O
ATOM      5  CB  ALA     1      68.014  53.813  16.378  1.00  0.00           C
ATOM      6  N   VAL     2      70.773  53.896  15.173  1.00  0.00           N
ATOM      7  CA  VAL     2      72.223  53.728  15.173  1.00  0.00           C
ATOM      8  C   VAL     2      72.922  55.066  15.188  1.00  0.00           C
ATOM      9  O   VAL     2      72.296  56.134  15.198  1.00  0.00           O
ATOM     10  CB  VAL     2      72.665  52.877  13.916  1.00  0.00           C
ATOM     11  CG1 VAL     2      74.193  52.762  13.664  1.00  0.00           C
ATOM     12  CG2 VAL     2      72.143  51.420  13.932  1.00  0.00           C
ATOM     52  H   VAL     2      70.368  54.810  15.182  1.00  0.00           H
ATOM     13  N   GLY     3      74.212  55.098  15.190  1.00  0.00           N
ATOM     14  CA  GLY     3      74.888  56.392  15.204  1.00  0.00           C
ATOM     15  C   GLY     3      76.388  56.218  15.203  1.00  0.00           C
ATOM     16  O   GLY     3      76.917  55.104  15.192  1.00  0.00           O
ATOM     53  H   GLY     3      74.741  54.249  15.182  1.00  0.00           H
ATOM     17  N   LEU     4      77.138  57.267  15.215  1.00  0.00           N
ATOM     18  CA  LEU     4      78.589  57.099  15.215  1.00  0.00           C
ATOM     19  C   LEU     4      79.288  58.437  15.230  1.00  0.00           C
ATOM     20  O   LEU     4      78.660  59.507  15.241  1.00  0.00           O
ATOM     21  CB  LEU     4      79.040  56.251  13.992  1.00  0.00           C
ATOM     22  CG  LEU     4      78.445  54.825  13.832  1.00  0.00           C
ATOM     23  CD1 LEU     4      79.094  54.113  12.637  1.00  0.00           C
ATOM     24  CD2 LEU     4      78.610  53.968  15.098  1.00  0.00           C
ATOM     54  H   LEU     4      76.733  58.181  15.224  1.00  0.00           H
ATOM     25  N   TRP     5      80.577  58.469  15.232  1.00  0.00           N
ATOM     26  CA  TRP     5      81.254  59.763  15.246  1.00  0.00           C
ATOM     27  C   TRP     5      82.754  59.589  15.245  1.00  0.00           C
ATOM     28  O   TRP     5      83.281  58.478  15.233  1.00  0.00           O
ATOM     29  CB  TRP     5      80.817  60.551  16.512  1.00  0.00           C
ATOM     30  CG  TRP     5      79.341  60.958  16.536  1.00  0.00           C
ATOM     31  CD1 TRP     5      78.304  60.258  17.187  1.00  0.00           C
ATOM     32  CD2 TRP     5      78.731  62.003  15.872  1.00  0.00           C
ATOM     33  CE2 TRP     5      77.340  61.929  16.130  1.00  0.00           C
ATOM     34  CE3 TRP     5      79.257  63.017  15.031  1.00  0.00           C
ATOM     35  NE1 TRP     5      77.049  60.848  16.950  1.00  0.00           N
ATOM     36  CZ2 TRP     5      76.464  62.876  15.556  1.00  0.00           C
ATOM     37  CZ3 TRP     5      78.372  63.946  14.482  1.00  0.00           C
ATOM     38  CH2 TRP     5      76.996  63.879  14.742  1.00  0.00           C
ATOM     55  H   TRP     5      81.105  57.620  15.224  1.00  0.00           H
ATOM     39  N   PRO     6      83.505  60.638  15.257  1.00  0.00           N
ATOM     40  CA  PRO     6      84.959  60.560  15.257  1.00  0.00           C
ATOM     41  C   PRO     6      85.498  60.036  16.599  1.00  0.00           C
ATOM     42  O   PRO     6      84.724  59.735  17.507  1.00  0.00           O
ATOM     43  CB  PRO     6      85.354  62.018  14.978  1.00  0.00           C
ATOM     44  CG  PRO     6      84.201  62.829  15.588  1.00  0.00           C
ATOM     45  CD  PRO     6      82.967  61.980  15.271  1.00  0.00           C
ATOM     56  H   PRO     6      83.072  61.539  15.267  1.00  0.00           H
ATOM     46  N   THR     7      86.770  59.913  16.774  1.00  0.00           N
ATOM     47  CA  THR     7      87.263  59.414  18.054  1.00  0.00           C
ATOM     48  C   THR     7      88.771  59.342  18.062  1.00  0.00           C
ATOM     49  O   THR     7      89.436  59.668  17.100  1.00  0.00           O
ATOM     50  CB  THR     7      86.639  58.014  18.376  1.00  0.00           C
ATOM     51  CG2 THR     7      87.080  57.358  19.702  1.00  0.00           C
ATOM     52  OG1 THR     7      85.224  58.110  18.479  1.00  0.00           O
ATOM     57  H   THR     7      87.409  60.156  16.044  1.00  0.00           H
TER      58      THR     7
END

pic :
Not sure I got the code above (in question) right, can't read it , had to use comments as trace.
Spent a while banging my head with : HAAD: A Quick Algorithm for Accurate Prediction of Hydrogen Atoms in Protein Structures
and trying to figure out their alghoritm, [code (https://seq2fun.dcmb.med.umich.edu//HAAD/haad.f90) is in FORTRAN90] from the paper , but no luck, if anybody wants to help eplaining this bit:

For constructing the sp2H2 and sp2H1 H-atoms, we first decide on the orientation of the conjugated plane or the aromatic ring with respect to the neighboring heavy atoms; the normal vector of the conjugated plane is determined by taking the cross product of two vectors between the heavy atoms. For the sp2H2 H-atoms (illustrated in Fig. 1b), the normal vector of the conjugated plane is the cross product of the unit vectors B→N and B2→B; then the two H-atoms are placed at positions A1 and H, which are within the conjugated plane respected to the B→N vector with the exact bond angle from CHARMM22 force field. For sp2H1 H-atoms, two conformations are possible. The first is to place the H-atom in the peptide plane as illustrated in Fig. 1b, where A1 and B1 represent the alpha carbon atoms. The position of H in this case is decided by using the same method as the one used to determine the position of the sp2H2 H-atoms while holding the trans-conformation.

it refers to Fig.1 of linked freely available article


text:
BUSCO results apparently inconsistent
I have two assemblies obtained from two different softwares for the same run data. In one of them, I get a BUSCO score of 69,4%. In the other one, I get 68,5%. The first assembly covers a 99,7% fraction of the reference while the second one covers 100%.
The relevant thing is that the reference genome also gets 68,5% (less than the first assembly). How can I interpret this situation?
I identify these options:

The first assembly is showing a different strain.
BUSCO is imprecise
The first assembly is merging DNA from another organism
The first assembly is a misassembly, repeating part of the sequence and BUSCO interprets that as new genes
The sample really is an organism that has more genes than the reference.

The organism is a bacteria and I am using bacteria_odb10 as the lineage dataset.
This is the BUSCO result for the reference (LN554883.1 Chlamydophila abortus strain 1H):
C:68.5%[S:68.5%,D:0.0%],F:8.1%,M:23.4%,n:124
85  Complete BUSCOs (C)
85  Complete and single-copy BUSCOs (S)
0   Complete and duplicated BUSCOs (D)
10  Fragmented BUSCOs (F)
29  Missing BUSCOs (M)
124 Total BUSCO groups searched

This for the first assembly:
C:69.4%[S:69.4%,D:0.0%],F:8.9%,M:21.7%,n:124
86  Complete BUSCOs (C)
86  Complete and single-copy BUSCOs (S)
0   Complete and duplicated BUSCOs (D)
11  Fragmented BUSCOs (F)
27  Missing BUSCOs (M)
124 Total BUSCO groups searched

And this for the second assembly:
C:68.5%[S:0.0%,D:68.5%],F:8.1%,M:23.4%,n:124
85  Complete BUSCOs (C)
0   Complete and single-copy BUSCOs (S)
85  Complete and duplicated BUSCOs (D)
10  Fragmented BUSCOs (F)
29  Missing BUSCOs (M)
124 Total BUSCO groups searched

I see now that the first assembly also contains a lot of contigs with human DNA. I guess BUSCO finds a match for some genes in the human DNA?

I see now that the first assembly also contains a lot of contigs with human DNA. I guess BUSCO finds a match for some genes in the human DNA?

The core idea behind BUSCO is that there are a number of genes that are Unique Single-Copy Orthologues (USCOs) that should be present in most genomes and should only appear once in the genome. This is only a guess, however, and not a universal rule.
Sometimes genomes can have multiple copies of these USCOs, and these duplicated genes reduce the overall score (which is based on the proportion of complete single-copy genes that are present as a fraction of the entire USCO set).
Regarding human DNA, bacteria are a curiously diverse group, with genomes that have more similarities in function than in form. Given this diversity, it wouldn't surprise me if contaminating DNA could increase the BUSCO score.

source: Stack Exchange

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

sort: votes


text:
Counting hexamers in fasta sequence and identify its structure (and interruptions)
I have a lot of fasta files, each one with thousand of reads containing the hexameric motif "CCCTCT". The hexameric motif is highly continuous in most cases but interruptions may occur. I need to count the hexameric motif keeping the read ID and identify their structure (CCCTCT)n. These interruptions are real and do not come from Sequencing error
Example Sequence:(interruption in bold)

ACATTTTTTTTCCACATCTGATGTGGAAAAAAAAAAAAATGAAATAGCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCTCCCTCTCCCTCTCCCCACGGTCTCCCTCTCATGCGGAGCCAAAGCTGGACTGTACTGCT

Strategy 1: Count the CCCTCT motif
grep -on "CCCTCT" fasta.fa | cut -f 1 -d ":" | sort | uniq -c

Output: 23 CCCTCT

As you can see we have one interruption (CCT) in this sequence which the codes doesn't detect.
Strategy 2: Identify reads with the same structure and count repeats in each read
from Bio import SeqIO
from collections import defaultdict

dedup_records = defaultdict(list)
for record in SeqIO.parse("fasta.fa", "fasta"):
    # Use the sequence as the key and then have a list of id's as the value
    dedup_records[str(record.seq)].append(record.id)
with open("Output.fasta", 'w') as output:
    for seq, ids in sorted(dedup_records.items(), key=lambda t: len(t[1]), reverse=True):
        # Join the ids and write them out as the fasta
        output.write(">{}_counts{}\n".format(ids[0], len(ids)))
        output.write(seq + "\n")

Output: list ranking of reads with identical structures but does not show the structure
>ID/ccs_counts2
CCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCTCCCTCTCCCTCT.....

Can you please give me an idea of how to obtain the following output? example of a sequence of 23 hexamers with interruptions.
(CCCTCT)20 [INTERRUPTION] (CCCTCT)3

I have used expansion hunter, repeatFinder and repeatAnalysis-tool by PacBio but they were not particularly sensitive for this case. Any help would be appreciated.

I've also mulled on this sort of puzzle in the past (when dealing with microsatellites) and I feel like there should be some elegant approaches that could fit.  How strict do you want the definition of a repeat region and interruption between repeats?  For example, more than one perfect tandem repeat triggers a match, and less than one repeat-length mismatch before it resumes counts as an interruption (rather than the end of the whole repeat)?  Or even simpler, just however many perfect matches with any-length interruptions between?

Thanks for your thoughts. I have to be either flexible to detect some changes (if there is some) and strict to be able to report those sequences following the patterns. I guess I just figured it out.

I was going to write some clever code counting k-mers, but this problem can easily be solved using the handy sequentially-operating feature of python's str.replace() method.
This code doesn't work if you need to look at the motif occurring on the other strand. Just the forward strand.
#!/usr/bin/env python3

mystring = "CCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCC\
TCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCT\
CCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCTCCCTCTCCCTCT"

motif = "CCCTCT"

# This replaces the motif with 1s
newstring = mystring.replace(motif, "1")
# results in 1111111111111111111111111111111111111CCT11
# silly, right?

# This whole block just parses that string and prints everything
#  the way that OP wanted.
done = False
space=""
accumulator = 0
while not done:
    if len(newstring) == 0:
        done = True
        break
    else:
        if newstring[0] == "1":
            accumulator += 1
        else:
            if accumulator != 0:
                print("{}({}){} ".format(space, motif, accumulator), end = "")
                accumulator = 0
            print(newstring[0], end = "")
            space = " "
    newstring = newstring[1::]

# This is our check after exiting the while loop.
#  This only fires if we stepped off cleanly at the
#  end of a motif.
if accumulator > 0:
    print("{}({}){} ".format(space, motif, accumulator), end = "")
print()

Prints, (CCCTCT)37 CCT (CCCTCT)2, matching OP's spec.

That looks good, but what if I have thousand of reads so I can use fasta/bam files as input. This looks fine when I have a read, am I right?

Not sure what you mean by "this looks fine when I have a read", but this code can be put into a method and worked into your framework you already have.

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


text:
Permutations of 4 amino acid sequence oligopeptide
I have an assignment for a project to build a tetrapeptides library to screen against some target enzymes in silicon.
I am looking for a script to construct about 160000 tetrapeptide from the 20Amino acids permutations.
I am quite familiar with Bash, but am struggling to get a function to working.

No "about" 20^4 (permutation with repetition) is 160_000, no approximation.

Hi @Ruthy and welcome to Bioinformatics.SE! We try always to be a helpful community forum, but we're most eager to help those willing to help themselves. Could you please update your question to describe why you need to solve this problem? What is the context? Is this for a homework assignment? Or a new research project? What have you tried so far? This forum is not designed to teach the very basics of programming, but if you have specific questions about a bioinformatics program you're trying to write we might be able to point you in the right direction.

Noted, @Matteo. And thank you, Daniel. It's an assignment for a project. I am building a tetrapeptides library to screen virtually against some target enzymes.  I am quite familiar with Bash. But, writing a script to achieve that has just been a problem. I hope I can just get a function to work with, then...  Thank you, once again.

Please post example of code that you tried. See also [help] and [ask].

has to be necessarily in bash? can't be python for example?

The Python solution is here. IMO Bash definitely isn't a good solution for this calculation.
from itertools import permutations, combinations

one2all ={'A': ('A', 'ALA', 'alanine'),
          'R': ('R', 'ARG', 'arginine'),
          'N': ('N', 'ASN', 'asparagine'),
          'D': ('D', 'ASP', 'aspartic acid'),
          'C': ('C', 'CYS', 'cysteine'),
          'Q': ('Q', 'GLN', 'glutamine'),
          'E': ('E', 'GLU', 'glutamic acid'),
          'G': ('G', 'GLY', 'glycine'),
          'H': ('H', 'HIS', 'histidine'),
          'I': ('I', 'ILE', 'isoleucine'),
          'L': ('L', 'LEU', 'leucine'),
          'K': ('K', 'LYS', 'lysine'),
          'M': ('M', 'MET', 'methionine'),
          'F': ('F', 'PHE', 'phenylalanine'),
          'P': ('P', 'PRO', 'proline'),
          'S': ('S', 'SER', 'serine'),
          'T': ('T', 'THR', 'threonine'),
          'W': ('W', 'TRP', 'tryptophan'),
          'Y': ('Y', 'TYR', 'tyrosine'),
          'V': ('V', 'VAL', 'valine'),
          'X': ('X', 'GLX', 'glutaminx'),
          'Z': ('Z', 'GLI', 'glycine'),
          'J': ('J', 'NLE', 'norleucine'),
          'U': ('U', 'CYC', 'cysteinc')}
myprotein = [x for x in one2all.keys()]
tetrapeps = combinations(myprotein,4)
peptide = []
count = 0
for x in list(tetrapeps):
    peptide.append(''.join(x))
for tetra in peptide:
    permpeptides = list(permutations (tetra))
    for i in permpeptides:
        count +=1
        print (''.join(i))

I get 255024 permutations so more than the estimate of 160000, however the code is correct (please see Validation below).
Truncated output

....
JUXV
UVXJ
UVJX
UXVJ
UXJV
UJVX
UJXV
VZJU
VZUJ
VJZU
VJUZ
VUZJ
VUJZ
ZVJU
ZVUJ
ZJVU
ZJUV
ZUVJ
ZUJV
JVZU
JVUZ
JZVU
JZUV
JUVZ
JUZV
UVZJ
UVJZ
UZVJ
UZJV
UJVZ
UJZV
XZJU
XZUJ
XJZU
XJUZ
XUZJ
XUJZ
ZXJU
ZXUJ
ZJXU
ZJUX
ZUXJ
ZUJX
JXZU
JXUZ
JZXU
JZUX
JUXZ
JUZX
UXZJ
UXJZ
UZXJ
UZJX
UJXZ
UJZX

I cut and paste the amino dictionary from here (I wasn't going to write it by hand), AminoAcid.py by Pasteur.

Validation
The theoretical number of permutations from 4 samples of 20 amino acids is

116280

However, the amino acid chart I used comprised 24 amino acids (I just copied and pasted it). The theoretical number for 24 amino acids is therefore:

255024

This was the number observed in the output of the above code and validated the calculation.
A permutations calutation can be performed online here
I am not sure how the expected output of 160000 was derived but suspect it was a guess.

The permutation equation is as follows, the output of this equation and the code above coincide,
$$ P_{n, r}=\frac{n !}{(n-r) !}$$
n = no. amino acids
r = size of oligopeptide

Permutation by repetition
If the permutation by repetition equation  is sought (comment below). The code is:
from itertools import product
myprotein = [x for x in one2all.keys()]
print(len([p for p in product(myprotein, repeat=4)]))

for 20 amino acids this is:

160000

for 24 amino acids:

331776

I think your calculations assume permutation without replacement. If you assume replacement occurs (i.e. an amino acid can occur more than once in a tetrapeptide), then the expected number of unique tetrapeptides is 20^4 = 160,000.

Firstly thanks for your answer @acvill. However, the equation you have used is permutation with repetition not the permutation equation. The question asked for just permutations. I've added the permutation equation below the answer and the permutation by repetition code below the equation.

you are correct that I did not use the permutation equation. From what I can tell, your approach is also more efficient and readily understood. However, it seems that there is a disconnect between the language that OP used ("permutations") and the question they are actually asking, which I read as "How to generate all 160,000 tetrapeptides?" It seems odd that a peptide library built for an enzymatic screen would exclude elements that have repeated residues, though, since this question is now months old, I don't expect we'll get any clarification from OP

Thanks @acvill. Now you raise the biological perspective its makes sense. I approached this from stats theory (ad infinitum) and the 'permutation equation' means only one thing. However, from a molecular biological perspective the permutation equation also means only one thing, as you've pointed out, but its a completely different equation. The OP's lecturers should have known this difference however. Epidemiology was one field where the perspective between biologists and maths/stats resulted in some very famous errors (I will avoid details). Anyway upvoted, thanks.

Here's a bash solution using brace expansion and the single-letter codes for the 20 common amino acids:
aalist=D,E,K,R,H,N,Q,S,T,P,G,A,V,I,L,M,C,F,Y,W
eval printf '"%s\n"' {"$aalist"}{"$aalist"}{"$aalist"}{"$aalist"}

This prints 160,000 lines, one for each tetrapeptide.


text:
What to use to determine SNPs in a population? How to find significant SNP variation between generations?
I have a population in the hundreds and did an Illumina pool-seq (ended up with ~100bp reads for the whole population).
I have the reference genome for the species and now need to align all the population reads to the reference genome, figure out which are SNPs and finally determine which ones are changing from one generation to the next (about 20generations difference).
Right now I'm trying to use popoolation2 https://sourceforge.net/p/popoolation2/wiki/Main/
but it seems to be a little old as it has not been updated in ~10 years.
In Galaxy, I couldn't find any pipelines dedicated to determine population SNP variation through generation.
So two questions:

what software to use to determine significant variation in SNPs in a population (in poolSeq)?
can you suggest me a pipeline for this? if one exists.

First focus on calling SNPs accurately. See the BWA, SAMtools, BCFtools pipeline (as described here) or GATK4 (as described here).
Once you have the [BV]CFs, then just code up the analysis to see changes over generations.


text:
Best R program to use for RNA seq analysis
We have an RNA seq taken from a the same tissue in two different developmental stages. We have one gene, with the number of reads from stage A and stage B, and we want to see which genes are differentially expressed.
After performing RNA-seq, I am left with a CSV file containing a list of genes and the amount of reads in two different cell types. We have one list of genes, with the number of reads in Stage A in one column and the number of reads in Stage B in another column. We are looking to see what types of genes are differentially expressed in the two stages.
I have used DESeq2, but I'm not sure if this is the best/correct to analyze what I am looking for.
What is the best package to analyze the data for significant differential gene expression?
Thanks

DESeq2 works for that.  Unless you literally have only two samples.
There is no statistically sound way of analyzing DE genes with literally only two samples.  This is a bad experimental design, suitable for only the roughest of exploratory purposes.  Only the most obvious differences (that you probably knew already) are likely to be borne out in a legitimate experiment.


text:
How to extract the protein sequences of a genbank file using R or biopython
sorry for the question, I'm trying to extract the proteins sequences from a genbank file.
 gene            complement(516466..532086)
                 /gene="rtxA"
                 /locus_tag="VV1_RS17390"
                 /old_locus_tag="VV2_0479"
 CDS             complement(516466..532086)
                 /gene="rtxA"
                 /locus_tag="VV1_RS17390"
                 /old_locus_tag="VV2_0479"
                 /inference="COORDINATES: similar to AA
                 sequence:RefSeq:WP_011081430.1"
                 /note="Derived by automated computational analysis using
                 gene prediction method: Protein Homology."
                 /codon_start=1
                 /transl_table=11
                 /product="MARTX multifunctional-autoprocessing
                 repeats-in-toxin holotoxin RtxA"
                 /protein_id="WP_011081430.1"
                 /translation="MGKPFWRSVEYFFTGNYSADDGNNSIVAIGFGGEIHAYGGDDHV
                 TVGSIGATVYTGSGNDTVVGGSAYLRVEDTTGHLSVKGAAGYADINKSGDGNVSFAGA
                 AGGVSIDHLGNHGDVSYGGAAAYNGITRKGLSGNVTFKGAGGYNALWHETNQGNLSFA
                 GAGAGNKLDRTWFNRYQGSRGDVTFDGAGAANSISSRVETGNITFRGAGADNHLVRKG
                 KEGNHTANLANEDISSANGYHSMGKGGYSLSDLHYSVNAVRSTSETVADIDEYTDQTL
                 FKPATDSGESSGDVRFNGAGGGNVIKSNVTRGNVYFNGGGIANVILHSSQFGHTEFNG
                 GGAANVIVKSGEEGDLTFRGAGLANVLVHQSKQGKMDVYAGGAVNVLVRIGDGQYLAH"

I'm using a code from a previous answer to extract the protein and the locus_tag, but I want some extra features like:
>locus_tag|old_locus_tag|gene|product|complement|protein_id|chromosome|Accession|organism|strain
MGKPFWRSVEYFFTGNYSADDG ..........

the result will be like:
>VV1_RS17390|VV2_0479|rtxA|MARTX multifunctional-autoprocessing repeats-in-toxin holotoxin RtxA"|516466..532086|WP_011081430.1|Chromosome-II|NC_004459.3|Vibrio vulnificus|CMCP6
MGKPFWRSVEYFFTGNYSADDG ..........

this is the code that I'm using, how to add those features, or an other way using R ???
from Bio import SeqIO
file_name = 'CMCP6.gb'
# stores all the CDS entries
all_entries = []
with open(file_name, 'r') as GBFile:
    GBcds = SeqIO.InsdcIO.GenBankCdsFeatureIterator(GBFile)

    for cds in GBcds:
        if cds.seq is not None:
            cds.id = cds.name
            cds.description = ''
            all_entries.append(cds)

# write file
SeqIO.write(all_entries, '{}.fasta'.format(file_name[:-3]), 'fasta')

But the result from the code only add the locus_tag and the sequence:
>VV1_RS17380
MFIKNLSIGKKIAAAFSIIAVINIAFGIFLSTELNTVKSELLNYTEDTLPAMEKVDAVRD
KISYWRRTQFAVFAMNDENQIKQTITRNEGIRREIETELAAYGKSVWPGEEEQTYNRLMS
LWSGYLSTMDKFNDALLAGDKDAAYPILTNSLSTFESIETEVNKLVMILKGAMDSNKNQI
LSSVNGLNTTAVISNIAIFAVMVIMTLLLTRIICGPLNIVVRQANAIARGDLSQDLDRNA
IGNDELGELADATIKMQDDLRQVIDNVIAAVTQLSSAVEEMTQISEMSASGMKDQQMQVT
LVATAMTEMKAAVADVARNTEDSASQAYDANRRTQDGAKETHQMVVSIEEVADIIAKAGD
TVAELEAQSNQINVVVDVIRGIADQTNLLALNAAIEAARAGESGRGFAVVADEVRTLAGR
TQDSTGEITSIIEKLQDLAKQAKHATEDSRTSIGACVEQGNNAQELMGSIEKSIANIADM
GAQIATACGQQDSVAEELSRNIENIHMASQEVAQGSQQTAQACRELTQLSVSLQDVMSRF
KLN
>VV1_RS17385
MNFKKTLLSIAIASASLTPAFSYSAPLLLDNTVHQTSQIAGANAWLEISLGQFKSNIEQF
KSHIAPQTKICAVMKADAYGNGIRGLMPTILEQQIPCVAIASNAEAKLVRESGFEGELIR
VRSASTSEIEQALSLDIEELIGSEQQARELASLAEKYSKTIKVHLALNDGGMGRNGIDMS
TERGPKEAVAIATHPSVAVVGIMTHFPNYNAEDVRTKLKSFNQHAQWLMESAGLKREEIT
LHVANSYTALNVPEAQLDMVRPGGVLYGDLPTNPEYPSIVAFKTRVASLHSLPAGSTVGY
DSTFTTANDAVMANLTVGYSDGYPRKMGNKAQVLINGQRANVVGVASMNTTMVDVSNIKG
VLPGDEVTLFGAQKNQHISVGEMEENAEVIFPELYTIWGTSNPRFYVK
>VV1_RS17390
MGKPFWRSVEYFFTGNYSADDGNNSIVAIGFGGEIHAYGGDDHVTVGSIGATVYTGSGND
TVVGGSAYLRVEDTTGHLSVKGAAGYADINKSGDGNVSFAGAAGGVSIDHLGNHGDVSYG
GAAAYNGITRKGLSGNVTFKGAGGYNALWHETNQGNLSFAGAGAGNKLDRTWFNRYQGSR
GDVTFDGAGAANSISSRVETGNITFRGAGADNHLVRKGKVGDITLQGAGASNRIERTRQA
EDVYAQTRGNIRFEGVGGYNSLYSDVAHGDIHFSGGGAYNTITRKGSGSSFDAQGMEYAK

Thanks so much

cds.description = ''

ok tried a little bit, not sure is the most elegant/fastest way of doing it
starting from file from : https://www.ncbi.nlm.nih.gov/nuccore/NC_004460.2?report=gbwithparts&log$=seqview        {GenBank(Full)}
here my code:
from Bio import SeqIO

from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord

file_name = 'CMCP6.gb'

# stores all the CDS entries
all_entries = []

gb_record = SeqIO.read(open(file_name,"r"), "genbank")

gb_feature = gb_record.features

accession = gb_record.annotations['accessions'][0]+'.'+str(gb_record.annotations['sequence_version'])

# accession = gb_record.annotations['accessions']
# print(gb_record.annotations['accessions'])

for seq_feature in gb_feature:
    if seq_feature.type=="source":
            organism = str(seq_feature.qualifiers['organism'][0]).rsplit(" ", 1)[0]
            strain = seq_feature.qualifiers['strain'][0]
            chromosome = seq_feature.qualifiers['chromosome'][0]

for seq_feature in gb_feature:
    if seq_feature.type=="CDS":
            if 'translation' in seq_feature.qualifiers:
                if seq_feature.qualifiers['translation'][0] != '':
                    if 'gene' in seq_feature.qualifiers:
                        gene = seq_feature.qualifiers['gene'][0]

                    else:

                        gene = 'not defined'

                    if 'product' in seq_feature.qualifiers:
                        product = seq_feature.qualifiers['product'][0]

                    else:

                        product = 'not defined'

                    if 'protein_id' in seq_feature.qualifiers:
                          protein_id = seq_feature.qualifiers['protein_id'][0]

                    else:

                          protein_id = 'not defined'

                    if 'old_locus_tag' in seq_feature.qualifiers:
                          old_locus_tag = seq_feature.qualifiers['old_locus_tag'][0]

                    else:

                          old_locus_tag = 'not defined'

                    if 'locus_tag' in seq_feature.qualifiers:
                          locus_tag = seq_feature.qualifiers['locus_tag'][0]

                    else:

                          locus_tag = 'not defined'

                    if seq_feature.location.strand == 1:
                          complement = 'no complement'

                    if seq_feature.location.strand == -1:
                        complement = str(seq_feature.location).strip('[]').split('(-)')[0]

                    pippo = SeqRecord(
                                Seq(seq_feature.qualifiers['translation'][0]),

                                id = seq_feature.qualifiers['locus_tag'][0] ,

                                description = ('|'+str(old_locus_tag)+'|'+str(gene)+'|'+str(product)\
                                                +'|'+str(complement)+'|'+str(protein_id) +'|Chromosome-'+str(chromosome)\
                                                +'|'+str(accession).strip("[]").replace("'", "")+'|'+str(organism)\
                                                    +'|'+str(strain))
                                )
                    all_entries.append(pippo)

# print(all_entries)
# write file
SeqIO.write(all_entries, '{}.fasta'.format(file_name[:-3]), 'fasta')

output for your example line is :
>VV1_RS17390 |VV2_0479|rtxA|MARTX multifunctional-autoprocessing repeats-in-toxin holotoxin RtxA|516465:532086]|WP_011081430.1|Chromosome-II|NC_004460.2|Vibrio vulnificus|CMCP6
MGKPFWRSVEYFFTGNYSADDGNNSIVAIGFGGEIHAYGGDDHVTVGSIGATVYTGSGND.....

I got extra space between locus_tag and | , I dont know how to get rid of it
Its not clear to me what you mean by Accessions, tried two ways (second one is past # in code) and got NC_004460.2 instead of NC_004459.3

thanks so much it works perfectly, I really appreciate your help !!!!


text:
Using a phased reference assembly to phase sample genotypes
I've been trying without success to phase the variants of about 30 samples in a 5Mb non-recombining region. I have a phased genome assembly of the region but so far I have not been successful in using that to inform the phasing of the samples which have been sequenced to high coverage with short reads. Non-reference panel phasing using Beagle 5.2 and Shapeit4 have produced haplotypes that are composed of many small phased blocks which switch phase between them, far to many to manually correct.
I have created a "reference panel" of SNPs of known phase for the region and want to use it to guide the rest of the phasing of the non-"reference panel" SNPs but every method I have found only uses the intersection of the reference panel and the samples for analysis. Is there some technique that I can use to phase the genotypes in the samples? I believe this is a tractable problem since recombination is suppressed in this region so haplotype blocks should be large and contiguous.

I'm unsure what you mean by "but every method I have found only uses the intersection of the reference panel and the samples for analysis". Do you mean the intersection of SNPs or individuals? In either case, I don't see why that would be an issue. Also, phasing with 30 individuals without a reference panel will give very poor results at best.

Thanks for the reply. My goal was to use the reference panel as a guide to phase the genotypes between the reference panel SNPs. A sort of "reverse imputation". Even though the sample size is small my thought was that since the region in question does not recombine (it is a homozygous lethal inversion) that the nearly complete linkage would make this a workable problem.

The fact that the reference assembly is phased doesn't do anything to help actually phase the 30 other samples. There is no way of knowing what the combination of variants on single haplotypes are in the 30 samples from the one reference genome + 30-sample VCF alone. With this info, it is possible to come up with a measure of "HaplotypeA from sample 1 is more similar to HaplotypeII than HaplotypeI from my phased reference", though.
To actually phase the 30 other samples it will take long reads and/or proximity ligation data to phase the variants over a 5Mb region. Hi-C, Chicago, or DNAse Hi-C/Dovetail Omni-C reads will be necessary to correctly phase all of the variants over a region that is as large as 5Mbp.
I recommend using the tool HapCUT2 to use the reads to phase the variants already called for each sample. When mapping the reads to your genome, only map the reads to a fasta file that contains a single haploid copy of the genome (2 homologous chromosomes/haplotypes  represented by 1 scaffold), not a fasta file that contains both haplotypes from the phased assembly (2 homologous chromosomes/haplotypes represented by 2 scaffolds).


text:
Can I traverse/enumerate PDB one at a time using a Python script without downloading all of them to my local disk?
What would be the total size of all PDB files in RCSB?
Can I traverse/enumerate them one at a time using a Python script without downloading all of them to my local disk?
If YES, what would be the Python script like?

The page https://www.wwpdb.org/ftp/pdb-ftp-sites gives the addresses of interest. PDB has an FTP site downloadable via rsync and the latter has a dry run mode.
rsync --port=33444 -rlpt -v -z --dry-run --stats rsync.rcsb.org::ftp_data/structures/divided/pdb/

You get a bunch of folders names and the following:
Number of files: 182479
Number of files transferred: 0
Total file size: 38382243431 bytes
Total transferred file size: 0 bytes
Literal data: 0 bytes
Matched data: 0 bytes
File list size: 3777395
File list generation time: 52.909 seconds
File list transfer time: 0.000 seconds
Total bytes sent: 16
Total bytes received: 3777651

sent 16 bytes  received 3777651 bytes  66861.36 bytes/sec
total size is 38382243431  speedup is 10160.30

So there are 182,479 *.ent.gz in various folders. The total is 38,382,243,431 bytes, which is 38 GB, so not too much. In Python it is easy to iterate across gzip files with the gzip module, so uncompressing them is generally a thing to avoid as it will bloat up greatly. The problem arise when dealing with the density maps, which do take up a lot of space...
Download without saving
The PDBe is easier to for downloading things than the RCSB PDB and the API are much more comprehensive.
Given a list of PDB codes one could
pdbcode : str = ''

import requests
r = requests.get(f'https://www.ebi.ac.uk/pdbe/entry-files/download/pdb{pdbcode.lower()}.ent', stream=True)
assert r.status_code == 200, f'something ({r.status_code} is wrong: {r.text}'
pdb_code = r.text

There is >180,000 PDB entries, see https://www.rcsb.org/stats/growth/growth-released-structures

Thanks for pointing that out! I thought that was suspicious. I forgot the recursive flag. Now the output is 38 GB, which is more reasonable. I don't have a folder to check in PDB format, but I thought it was a lot less...

I'll only add that there are entries (currently, 2375) that have coordinates in the mmCIF format (which is the primary format) and not in the PDB format. MmCIF files are somewhat bigger, the compressed archive is 52GB.

is rsync a Python builtin ??

It's a classic Unix command

Not Python, but you can do this with BioStructures.jl in Julia:
using BioStructures
pdblist = pdbentrylist()
for pdbid in pdblist
    # Download a file and delete it after use
    downloadpdb(pdbid, format=MMCIF) do fp
        s = read(fp, MMCIF)
        # Do something with the file, e.g. print the cysteine fraction
        nres = countresidues(s, standardselector)
        if nres > 0
            frac = countresidues(s, standardselector, x -> resname(x) == "CYS") / nres
            println(pdbid, "  ", round(frac; digits=2))
        end
    end
end

If you want to get all the PDB entries you should use the mmCIF format, though it does lead to larger file sizes.

my attempt working as of the 17th January 2023 using RCSB PDB Search API:
import json

# import urllib.parse

import requests

import sys

# data = {"query":{"type":"terminal","service":"text","parameters":{"attribute":"rcsb_entry_container_identifiers.entry_id",
            # "operator":"exists","negation":False}},"return_type":"entry"}

data = {"query":{"type":"terminal","service":"text","parameters":{"attribute":"rcsb_entry_container_identifiers.entry_id",
             "operator":"exists","negation":False}},"request_options":{"return_all_hits":True},"return_type":"entry"}

data = json.dumps(data)

headers = {'content-type': 'application/json'}
url = 'https://search.rcsb.org/rcsbsearch/v2/query?json='

params = {}

risposta = requests.post(url, params=params, data = data, headers=headers)

# print(risposta.status_code,'\n', risposta.content.decode())

rispo = json.loads(risposta.content.decode())

# print(rispo["result_set"], type(rispo["result_set"]), len(rispo["result_set"]))

id_list = []
for i in rispo["result_set"] :

    id_list.append(i['identifier'])

print(id_list, type(id_list) , 'lenght : ', len(id_list), 'size bytes  : ' ,sys.getsizeof(id_list))

you get a huge list :
.............................
'2V0H', '2V0G', '2V0F', '2V0E', '2V0D', '370D', '2ULL', '2UKD', '2UGI', '2UDP', '2UCZ', '2UBP', '2UAG', '2TYS', '369D', '368D', '367D', '366D', '2TUN', '2TSS', '2TSR', '2TSY', '365D', '2TRT', '2TRS', '2TS1', '2TRY', '2TRX', '364D', '2TSC', '2TSB', '2TSA', '363D', '2TRC', '2TRA', '2TRM', '2TRH']

<class 'list'> lenght :  200069 size bytes  :  1671784

Using the requests GET endpoints :
import requests

import sys

url = 'https://search.rcsb.org/rcsbsearch/v2/query?json=%7B%22query%22%3A+%7B%22type%22%3A+%22terminal%22%2C+%22service%22%3A+%22text%22%2C+%22parameters%22%3A+%7B%22attribute%22%3A+%22rcsb_entry_container_identifiers.entry_id%22%2C+%22operator%22%3A+%22exists%22%2C+%22negation%22%3A+false%7D%7D%2C+%22request_options%22%3A+%7B%22return_all_hits%22%3A+true%7D%2C+%22return_type%22%3A+%22entry%22%7D'

risposta = requests.get(url)

print('+++++++++++++++++++++++++++++')

rispo = risposta.json()

id_list = []
for i in rispo["result_set"] :

    id_list.append(i['identifier'])

print(id_list, type(id_list) , 'lenght : ', len(id_list), 'size bytes  : ' ,sys.getsizeof(id_list))

print('--------------------------------------')

print(risposta.status_code)

print('\n\n',risposta.url)

getting as of the 19th January 2023 --> 200375 entries


text:
building complex drug-dna for AMBER software
I will appreciate if you can please clarify some of my doubts about drug-DNA complex.
I want to study the drug-DNA simulation using AMBER. I did go through all the tutorials video on youtube but everything is about protein-ligand using GROMACS. I can also use GROMACS if that gives me the results.
My DNA sequences very small (12-mer) and I want to to intercalate it with a drug molecule.
I generated the DNA sequence in Chimera. My questions are as follows.

Is it possible to insert my drug at specific sequence and then do the binding energy calculation using AMBER?
(1) How can I insert it and make a complex?
(2) If I have to create a complex in AMBER how can I do that---reference/tutorials link will be     great.

The publications I read, they downloaded the PDB file with the drug molecule and extracted the  crystallographic information where the drug is attached and then they extended the sequence from there using NAB tools in AMBER. Can somebody please explain how can we do that.

Thank you so much for your time.

I want to intercalate the drug molecule at the desired sequence I want and calculate the binding energy.

the tutorials video on youtube

Gromacs has great written documentation, so I'd suggest referring to that directly.

I can also use GROMACS if that gives me the results.

Gromacs works fine with a DNA-ligand system. There are different Amber forcefields, proteins are generally used with AMBER99SB-ILDN simply because it has the longest name and is one of the inbuilt ones. Technically, the Amber forcefield tailored towards DNA is OL15 (recommended here) and gaff2 is for a generalised ff for ligands. These aren't natively within Gromacs, but been ported to Gromacs.
Whereas pi-pi iterations are somewhat covered by the Lenard-Jones term, the orientation will be driven not only by dipole, but also quadrupole moments, which aren't covered by classical mechanics forcefields. Likewise S-pi and cation-pi interactions will not be covered. Most DNA intercalating small molecules have multiple fused rings, e.g. doxorubicin, that pi stack with the DNA. Consequently, a simulation under a classical mechanics forcefield is likely going to be gibberish.

Is it possible to insert my drug at specific sequence

See your previous question and the easiest option of dragging it in PyMOL

Can somebody please explain how can we do that.

Read the AmberTools21 documentation.


text:
Getting translation table of organisms through Entrez direct (or other CLIs)
I have a long list of NCBI taxonomic ids for organisms and I want to write a script to get their corresponding translation tables. Is there a way to do this through Entrez?
(packages in R, Python etc. would be fine too)

I came up with this idea. So, the information for the translation table of any particular organism can be found on its page on the Taxonomy Browser: https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi
If I search a taxid into the website, eg 7227, the page for Drosophila melanogaster will come up but the URL will remain unmodified. However, you can get a modified URL specific to a taxid by appending ?mode=Info&id=your_taxid to the above URL. Eg for D. melanogaster again, it has this unique URL for its Taxonomy browser page: https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?mode=Info&id=7227
Since you have a list of taxids, you can programmatically access Taxonomy browser pages by modifying the URL in the way I said. To programmatically get the translation tables of D. melanogaster with this approach, use the urllib package in Python.
import urllib.request
url = "https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?mode=Info&id=7227"
weburl = urllib.request.urlopen(url)
data = weburl.read().decode("utf-8")  # decode because page returned in binary

data_s = data.split('<')

for d in data_s:
    if "Translation table" in d:
        print(d.split('>')[1])

The output is:
Translation table 1 (Standard)
Translation table 5 (Invertebrate Mitochondrial)

This will work for one taxid. But you have a list of them. However, with this approach in mind, it's elementary to construct a Python function for getting the Translation table of all your taxids.
import urllib.request

my_taxids = ['7227','9606','10090']

def get_translation_table(taxid):
    url = "https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?mode=Info&id=" + taxid
    weburl = urllib.request.urlopen(url)
    data = weburl.read().decode("utf-8")

    data_s = data.split('<')

    tr_table = []

    for d in data_s:
        if "Translation table" in d:
            tr_table += [d.split('>')[1]]

    return tr_table

your_tables = list(map(get_translation_table,my_taxids))

print(your_tables)

The output is:
[['Translation table 1 (Standard)', 'Translation table 5 (Invertebrate Mitochondrial)'], ['Translation table 1 (Standard)', 'Translation table 2 (Vertebrate Mitochondrial)'], ['Translation table 1 (Standard)', 'Translation table 2 (Vertebrate Mitochondrial)']]

And if you want the results nicely formatted into a pandas dataframe, just add:
import pandas as pd

To the top of the script, and
t = pd.DataFrame(your_tables,columns=["Nuclear","Mitochondrial"])
print(t)

To the bottom. Output of this print (if you have too many taxids use the head method instead so not everything prints):
                          Nuclear                                     Mitochondrial
0  Translation table 1 (Standard)  Translation table 5 (Invertebrate Mitochondrial)
1  Translation table 1 (Standard)    Translation table 2 (Vertebrate Mitochondrial)
2  Translation table 1 (Standard)    Translation table 2 (Vertebrate Mitochondrial)

You can find them here at least:
https://www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprintgc.cgi
You can download the complete list of codes (using FTP) here.
However, to link a specific organism to a specific codon usage table, you'll need to do a bit of extra work. Perhaps the script here can give you a starting point?

What I was really looking for was some way to input the taxonomic id/name of an organism and have its translation table returned. NCBI does seem to have this data for specific taxa, because a button on the page you linked ('Genetic Code in Taxonomy Tree') matches organisms with translation tables.

Yeah, you'd need to code that together using an ASN.1 parser I guess?


text:
Given SNP-Chip data for a population, what tool should I use to reconstruct haplotypes?
I think Eagle or haploSep will do the job of reconstructing haplotype data, but I'm not sure what the 'standard' or best practice tool to use is.
I have a VCF/PLINK format file of ~8,000 individuals with ~150,000 SNPs each (from across a population with several sub-populations).
I want to use the resulting data to place new samples into one or more populations (e.g. ancestry analysis).
Before I go trying to make one of these tools work / look at accuracy, etc., am I doing it right in the first place?

I've answered, but you should provide some more information to get a proper answer. In particular, i) the number of individuals, ii) the number of SNPS, iii) what you plan on doing with the haplotypes once you have obtained them. All of these things matter to a different degree. Cheers.

I've tried to add more detail, but he downstream analysis isn't fixed yet... Each 'sub-population' is inbreeding, and I want to assign an individual to one or a mix of each sub-population... Perhaps I should just be using some other analysis?

That sounds reasonable, but unless you are using a method which requires phased data (e.g. ChromoPainter/fineSTRUCTURE), you don't need to infer the haplotypes. You could use something like ADMIXTURE for that (but you get less power than when using haplotype-based methods).

HaploSep uses allele frequency data - I presume you have genotype data.
In that case, I would suggest looking at either Beagle5.2, shapeit4 and Eagle2. I believe they more or less give the same accuracy (Eagle2 is slightly less accurate than the other two) and Beagle 5.2 is faster than the rest. Depending on what you plan on doing with the haplotypes, I don't think that the choice of phasing algorithm will affect downstream results a great deal.


text:
In Galaxy, using the bracken tool, how do I get a Kmer distribution?
As you can see in the image, the Kmer distribution field is empty: "No options available"
What am I missing?
Is there another tool I must use first to generate a Kmer distribution?

Try posting here instead: https://help.galaxyproject.org/


text:
Understanding exercise on file coverage with question on summary statistics
I'm doing an exercise that asks for two files:

Input 1: A target file (.bed format) contains multiple regions from
chr7:40000000-50000000 of human reference genome GRCh37 (hg19)
Input 2: Refseq exon list file (.bed format) for all human coding genes
(hg19 position)

The final goal is:
For all genes located in chr7:40000000-50000000, get the summary statistics of the target file coverage. (For each gene, get the fraction of exonic bases that was covered by the target file).

I believe what they refer in this exercise is that the target file
should be something like the whole chr sizes of hg19 hg19.chr.sizes and refseq_exon_list the list of exons from reqseq database. Both can be
downloaded from tools like table browser .Is that correct ?

I'm not sure which files I should download here to perform this task.
Once downloaded the file I believe what I need to do is to restrict
the refseq_exon_list by the requested region and then perform
coverage with something like bedtools. Something on those lines
bedtools coverage -a hg19.chr.sizes -b reqseq_exon_list.

Am i right here ? Any input is appreciated, thanks.

Based on the description in your question, I don't think this assumption is correct:

I believe what they refer in this exercise is that the target file should be something like the whole chr sizes of hg19 hg19.chr.sizes [. . .]

I understand that the question is asking you to simply produce a random selection of sub-regions from chromosome 7, all of which fall within the larger region 40000000-50000000. It can't be the list of chromosome sizes since, first of all, the question is asking you to get regions from chr17 only, and second because that just wouldn't make sense: only chr17 is relevant here, so the rest would be useless and, in any case, the length of the chromosome itself isn't helpful information.
The idea here is to demonstrate that given a set of genomic regions and a list of start and end coordinates of exons, you can calculate how many of your exons overlap the target regions. Therefore, the steps to follow would be:

Generate random intervals in the 40000000-50000000 range.
Get the bed file for all human coding genes on hg19. You should indeed be able to download this from the UCSC genome browser, among other places.
Use a tool, bedtools or build your own using awk, or anything else that takes your fancy and for each gene found in the chr7:40000000-50000000 range, calculate how many bases of its exons fall within your randomly chosen selection of regions.

Of course, the best thing to do would be to ask your teacher to clarify.

This sounds like homework - can you post what you have tried so far + a few lines of the input files + your code? Also, please add the "homework" tag if this is homework, or a task from a learning exercise.
Without giving it away think about the information in your Refseq exon list file. Is it a 3-field bed file with just coordinates, or is it the expanded file format with fields 4-12? Do you really need to download the RefSeq annotation, or can you parse the information in your bed file to answer your question?
It sounds like this is your ultimate goal:

For each gene, get the fraction of exonic bases that was covered by the target file.

What information from the output of bedtools coverage will give you this answer? Does this information exactly match the spec you provided? Are there other tools in bedtools you could use to achieve the same thing?
Side note: It isn't 100% clear what you are asking and what your inputs are. I think the policy here isn't to outright answer homework questions without seeing evidence of some attempts. Someone can provide a better answer if you can update your question.

The problem is that I don't have the input files. The question implicitly hints on what should be downloaded from public resources but I don't know what to download in the first place =(. If I knew it then surely my attempt would be here.

Download the refseq human genome annotation in gff format. You can find that in the older versions here: https://www.ncbi.nlm.nih.gov/genome/?term=human . You can use this gff file to do all of the things that you mentioned above.


text:
Perl convert GFF3 file
I have the following GFF3 file which I would like to convert to ZFF file with the below script:
$ head -n 20 traingenes.gff3
##gff-version 3
# gffread v0.12.7
# gffread -E traingenes.gtf -o traingenes.gff3
NbLab330C00 GeneMark.hmm3   mRNA    74501   76501   .   +   .   ID=1_t;geneID=1_g
NbLab330C00 GeneMark.hmm3   CDS 74501   74512   .   +   0   Parent=1_t
NbLab330C00 GeneMark.hmm3   CDS 75568   75999   .   +   0   Parent=1_t
NbLab330C00 GeneMark.hmm3   CDS 76151   76501   .   +   0   Parent=1_t
NbLab330C00 GeneMark.hmm3   mRNA    76535   77079   .   +   .   ID=2_t;geneID=2_g
NbLab330C00 GeneMark.hmm3   CDS 76535   76639   .   +   0   Parent=2_t
NbLab330C00 GeneMark.hmm3   CDS 76702   77079   .   +   0   Parent=2_t
NbLab330C00 GeneMark.hmm3   mRNA    93763   100703  .   -   .   ID=3_t;geneID=3_g
NbLab330C00 GeneMark.hmm3   CDS 93763   93837   .   -   0   Parent=3_t
NbLab330C00 GeneMark.hmm3   CDS 93915   94031   .   -   0   Parent=3_t
NbLab330C00 GeneMark.hmm3   CDS 94351   94430   .   -   2   Parent=3_t
NbLab330C00 GeneMark.hmm3   CDS 95483   95589   .   -   1   Parent=3_t
NbLab330C00 GeneMark.hmm3   CDS 95697   95746   .   -   0   Parent=3_t
NbLab330C00 GeneMark.hmm3   CDS 100464  100703  .   -   0   Parent=3_t

#!/usr/bin/env perl
use strict;
#
# source: https://biowize.wordpress.com/2012/06/01/training-the-snap-ab-initio-gene-predictor/
#
my @exons;
my $gene_count = 0;
my $current_seq = "";
while (my $line = <STDIN>) {
  if ($line =~ m/^###/) {
    flush(\@exons);
    next;
  }
  my @fields = split(/\t/, $line);
  if ($fields[2] eq "mRNA") {
    flush(\@exons);
  } elsif ($fields[2] eq "exon") {
    if ($fields[0] ne $current_seq) {
      $current_seq = $fields[0];
      printf(">%s\n", $current_seq);
    }
    push(@exons, \@fields);
  }
}
flush();

sub flush
{
  my $num_exons = scalar(@exons);
  return if($num_exons == 0);

  my $group = sprintf("%s.%d", $exons[0]->[0], $gene_count);
  $gene_count++;

  if ($num_exons == 1) {
my($start, $end) = ($exons[0]->[3], $exons[0]->[4]);
if ($exons[0]->[6] eq "-") {
  ($start, $end) = ($exons[0]->[4], $exons[0]->[3]);
}
printf("Esngl\t%lu\t%lu\t%s\n", $start, $end, $group);
  } else {
@exons = reverse(@exons) if($exons[0]->[6] eq "-");
for (my $i = 0; $i < $num_exons; $i++) {
  my $exon_type = "Exon";
  if ($i == 0) {
    $exon_type = "Einit";
  } elsif ($i == $num_exons - 1) {
    $exon_type = "Eterm";
      }

      my($start, $end) = ($exons[$i]->[3], $exons[$i]->[4]);
      if ($exons[0]->[6] eq "-") {
    ($start, $end) = ($exons[$i]->[4], $exons[$i]->[3]);
      }

      printf("%s\t%lu\t%lu\t%s\n", $exon_type, $start, $end, $group);
    }
  }
  @exons = ();
}

However, when I run it like perl gff3_to_zff.pl < traingenes.gff3 >> traingenes.ann  an empty file get generated.
What did I miss?
Thank you in advance,
Michal

Are you going for the short format or long format? Can you provide a spec for what you'd like the output to look like? ZFF is not a common file format.

I would go the long format. The specs are described here.

Please look at the responses below @user977828 and consider upvoting/accepting an answer

Your file looks like it is CDS feature based, not exon feature based, which seems to be what the script expects.
There is a tool to convert CDS based GFF to exon based that I've used successfully in the past...

Could be GFF utilities?
Actually, I think it's GenomeTools

Originally, I had a GTF but I failed with gffread --gene2exon -E traingenes.gtf -o traingenes.gff3.  What did I miss?

What do you mean 'failed'?

Sorry, I meant it still kept CDS rather than changing them to Exon. What did I miss?

The easiest is to use agat_convert_sp_gff2zff.pl from AGAT
(Edit Actually you will have to first standardize your file because you do not have gene features, neither gene_id or Parent atttributes in your mRNA features. AGAT will be confused. So, in your case, prior the gff2zff conversion you will have to run agat_convert_sp_gxf2gxf.pl -c geneID --gff file.gff -o file_clean.gff)
Otherwise you can use your script but you need to add the exon as suggested by @Dan.
To do so you can just replace the CDS by exon using a sed command, it should work.
You can also use agat_convert_sp_gxf2gxf.pl (with parameter -c geneID) it will add the gene and exon features.


text:
Use biopython to align SeqRecords stored in dict
I'd like to perform multiple alignments, where a gene from each sample was read in from fasta files. The fasta file represented one sample and had multiple genes. I have read in each sample fasta file and now have a dictionary of genes and their samples and sequences.
Here is a small example.
The main key is the gene ID for another dictionary as the value. In the dictionary 'R??' is the sample name and key,
followed by the annotation and sequence in a list.
{'PRELSG_MIT_v1:5894-6032': {'R47': ['not annotated',
                                     Seq('ATAAAAATAGTTACCATAGCTGTTGATGGATGCTTCGATCATAAAATATGTTTG...GGA')],
                             'R48': ['not annotated',
                                     Seq('ATAAAAATAGTTACCATAGCTGTTGATGGATGCTTCGATCATAAAATATGTTTG...GGA')],
                             'R49': ['not annotated',
                                     Seq('ATAAAAATAGTTACCATAGCTGTTGATGGATGCTTCGATCATAAAATATGTTTG...GGA')],
                             'R50': ['not annotated',
                                     Seq('ATAAAAATAGTTACCATAGCTGTTGATGGATGCTTCGATCATAAAATATGTTTG...GGA')],
                             'R51': ['not annotated',
                                     Seq('ATAAAAATAGTTACCATAGCTGTTGATGGATGCTTCGATCATAAAATATGTTTG...GGA')],
                             'R52': ['not annotated',
                                     Seq('ATAAAAATAGTTACCATAGCTGTTGATGGATGCTTCGATCATAAAATATGTTTG...GGA')]},
 'exon_PRELSG_MIT00300.1-E1': {'R47': ['cytochrome b',
                                       Seq('TCCAGATAATGCTATTACAGTAGATAGATATGCTACTCCTTTACATATTGTTCC...TTC')],
                               'R48': ['cytochrome b',
                                       Seq('TCCAGATAATGCTATTACAGTAGATAGATATGCTACTCCTTTACATATTGTTCC...TTC')],
                               'R49': ['cytochrome b',
                                       Seq('TCCAGATAATGCTATTACAGTAGATAGATATGCTACTCCTTTACATATTGTTCC...TTC')],
                               'R50': ['cytochrome b',
                                       Seq('TCCAGATAATGCTATTACAGTAGATAGATATGCTACTCCTTTACATATTGTTCC...TTC')],
                               'R51': ['cytochrome b',
                                       Seq('TCCAGATAATGCTATTACAGTAGATAGATATGCTACTCCTTTACATATTGTTCC...TTC')],
                               'R52': ['cytochrome b',
                                       Seq('TCCAGATAATGCTATTACAGTAGATAGATATGCTACTCCTTTACATATTGTTCC...TTC')]}}

I'd like to align each gene and write out those alignments to file without any intermediate files written to disk.
How do I do that?

do you want to align each of those Seq('ATAAAAATAGTT...... on top of each other ?

Yes, what you said moth. I want one alignment for each gene.

I don't believe this is possible. So I'm answering this to close it, since I don't have enough reputation to close it myself.

This BioPython document should get you started.
A specific example is given here.

All the examples in documentation you linked require the files to be ready to be aligned, where each file represents only one gene and multiple samples. That is not what I have here, and I don't want to write intermediate gene files. I want to do it all in memory until the final alignments are ready.

Strange requirement. No idea how to align in python without calling another application by file.

I learnt how to create a MultipleSeqAlignment object based on that dict, but I can't see how to pass that object to a muscle subprocess...
I guess I'd need a muscle python library.

I don't think muscle can work 'file free'... may be wrong.

Yeah, I've given up on this. Biopython supports this using pairwise alignments, but that's obviously not going to work here. Thanks for your help.

Good, what you're trying is very weird.


text:
Simulation of DNA sequences through substitution rates
I'm looking for a little bit of guidance.
My question is regarding the simulation of DNA sequences with a fix substitution rate.
The majority of the programs for simulating sequences use Continuous Markov Chain Models, with different instantaneous rate matrices that can be modified according to the free parameters of each model (JC69, HKY, GTR, etc..). For example, in HKY the free parameters are the transition/transversion (kappa) rate and nucleotide frequencies (pi).
In these programs (for example INDELible, Seq-Gen, Pyvolve), the substitution rate is expected to be represented by the branch length of the phylogenetic tree used to create the simulations.
INDELible paper.
"..Rate matrices are rescaled by INDELible such that the branch lengths represent the expected number of substitutions per site (or the average expected number of substitutions per site under a heterogeneous-sites model)." (P2 - Simulation of substitutions)
INDELible: a flexible simulator of biological sequence evolution
"...Each branch length is assumed to denote the mean number of nucleotide substitutions per site that will be simulated along that branch..." (P2 - Algorithm)
Seq-Gen: an application for the Monte Carlo simulation of DNA sequence evolution along phylogenetic trees
Let's say I want to simulate sequences with a mean substitution rate of 2.5 x 10^-7 per year. How would I prepare a tree that can accurately represent that substitution rate and be used in the simulation? In other words, how can the length be represented as a substitution rate accurately?
Any insight is appreciated!.
Thanks in advance!

UPDATE based on comment
Simplest answer: just rescale your substitution rate from years to units of branch length. (Stop thinking in terms of years.)
original answer
The short answer is that you have to calibrate the branch lengths of the tree to yield the substitution rates per year you want. (Leave the topology alone.)
Tree branch lengths are usually agnostic to chronological time (e.g. years). Instead it is very focused on whatever the model used to construct the tree is. This can be simple substitution counts (e.g. Jukes-Cantor), or a sophisticated model that takes into account a lot of information about relative rates of different changes, possibly learned from the dataset itself (e.g. GTR).
You yourself will have to decide what the conversion factor is between your branch lengths (presumably derived from such a substitution model or a distance metric) and units of chronological time. In some trees, it may already be calibrated to units of chronological time. But those will be rare.
So what I'd do is to make sure that the branch lengths of the tree are of a length that will, on average, yield the number of substitutions that you expect based on how long (in years) a given branch on the tree is supposed to be. For example, if you have a tree of primates, then you want root to tip distances to be on the order of tens of millions of years, probably. If you have a tree of all bacteria, you want the root-to-tip distances to be on the order of billions of years.
Note please that this is identical to re-scaling the substitution matrix to substitutions per year. What is crucial is that you know what you mean by "substitutions per year" and how you are defining that in terms of the tree.

Please correct me if I'm wrong. Considering that branch length is agnostic to time, if i have a tree with certain "value" of length (e.g 0.04),Is the substitution model that defines the "time" itself? (years or millions of years)?. If I only want to simulate DNA sequences in a defined topology, it would be better to obtain the model parameters from a real dataset? (With IQTREE for example). Thanks again for your help

@Kako Note that there is no scientific reason why chronological time has to be represented in a phylogeny or a substitution matrix. the substitution model is also usually thinking only in terms of branch length as the time scale (for phylogenetics at least). Neither branch lengths nor substitution models are usually understood in terms of years. So if you want your simulation to use years, then your value 0.04 should be years to match your substitution model. Or, you could just rescale your substitution rate to be subs/unit branch length (much simpler, what all those programs assume).


text:
list comprehension and codon translation
I'm trying to create a python program that follows the codon translation and turn mRNA sequence into amino acids (see the dictionary). The system has to return an error of the string is not a multiplier of 3, and the the function should not produce any more amino acids once a stop codon (keys with a value of X in the dictionary) is found, and should not display 'X' in the output. I do not know how to incorporate the stop codon to stop the programming entirely, as well as the multiplier of 3 aspect in the code. Can someone give me any advice? Thank you so much!
example: protfromRNA('GCCAACGCACCCACUUUUUUUUAGAAA') should return ANAPTFF
def protfromRNA(Seq):
    GeneticCode = {'AAA':'K', 'AAC':'N', 'AAG':'K', 'AAU':'N',
               'ACA':'T', 'ACC':'T', 'ACG':'T', 'ACU':'T',
               'AGA':'R', 'AGC':'S', 'AGG':'R', 'AGU':'S',
               'AUA':'I', 'AUC':'I', 'AUG':'M', 'AUU':'I',
               'CAA':'Q', 'CAC':'H', 'CAG':'Q', 'CAU':'H',
               'CCA':'P', 'CCC':'P', 'CCG':'P', 'CCU':'P',
               'CGA':'R', 'CGC':'R', 'CGG':'R', 'CGU':'R',
               'CUA':'L', 'CUC':'L', 'CUG':'L', 'CUU':'L',
               'GAA':'E', 'GAC':'D', 'GAG':'E', 'GAU':'D',
               'GCA':'A', 'GCC':'A', 'GCG':'A', 'GCU':'A',
               'GGA':'G', 'GGC':'G', 'GGG':'G', 'GGU':'G',
               'GUA':'V', 'GUC':'V', 'GUG':'V', 'GUU':'V',
               'UAA':'X', 'UAC':'Y', 'UAG':'X', 'UAU':'Y',
               'UCA':'S', 'UCC':'S', 'UCG':'S', 'UCU':'S',
               'UGA':'X', 'UGC':'C', 'UGG':'W', 'UGU':'C',
               'UUA':'L', 'UUC':'F', 'UUG':'L', 'UUU':'F'}
    return ''.join(GeneticCode[i] if GeneticCode[i]!=X else '' for i in Seq)

You need to iterate over the sequence 3 characters at a time and you need to stop when you hit the stop codon (or before for an invalid codon, e.g. of length 2). Here is a memory efficient implementation using iterators (requires python >= 3.8):
from itertools import islice

def protfromRNA(Seq):
    it = iter(Seq)
    lst = []
    while codon := "".join(islice(it, 3)):
        aa = GeneticCode.get(codon)
        if aa is None:
            raise ValueError(f"invalid codon: {codon}")
        if aa == "X":
            break
        lst.append(aa)
    return "".join(lst)

print(protfromRNA('GCCAACGCACCCACUUUUUUUUAGAAA'))
# ANAPTFF


text:
Snakemake Fastqc: "Multiple run or shell keywords in rule run_fastqc."
I am trying to check the quality of RNA-Seq data from Illumina using fastqc in snakemake in a conda environment. I get the error "Multiple run or shell keywords in rule run_fastqc".
Snakefile:
rule run_fastqc:
    input:
        "/[RNA fastq file path]/well01_S001_L001_R1_001.fastq.gz"
    output:
        "./[out path]/well01_S001_L001_R1_001.fastqc.html",
        "./[out path]/well01_S001_L001_R1_001.fastqc.zip"
    wrapper:
        "0.79.0/bio/fastqc"
    shell:
        "fastqc {input} --outdir=./[out path]/"

I've also tried:
...
    shell:
        "fastqc {input} -o ./[out path]/"
...
    shell:
        "fastqc -o ./[out path]/ {input}"
...
    shell:
        "fastqc {input}"
        "fastqc -o ./[out path]/"
...
    shell:
        "fastqc {input} > -o ./[out path]/"

Code:
snakemake --cores 1 well01_S001_L001_R1_001.fastqc.html well01_S001_L001_R1_001.fastqc.zip

Output:
SyntaxError in line 9 of /[path to Snakefile]/Snakefile:
Multiple run or shell keywords in rule run_fastqc. (Snakefile, line 9)

When I run fastqc alone (not in snakemake), it works and I get the two output files. How can I get snakemake to work?
Thank you for your time and help!

From a quick look at the manual, it looks like you are supposed to use just the wrapper, not combine the wrapper with a shell. Can you show me an example of such combined usage?
The error seems to suggest that the wrapper already has a shell directive that you shell conflicts with.

Thank you! In Output, I also had to remove the "./[out path]/" before the output file name, I guess the "--outdir" in Shell is sufficient. I got the files but still received the error: "Job Missing files after 5 seconds:
well01_S001_L001_R1_001.fastqc.html
well01_S001_L001_R1_001.fastqc.zip
This might be due to filesystem latency. If that is the case, consider to increase the wait time with --latency-wait.
Job id: 0 completed successfully, but some output files are missing. 0" and I also waited for 60s and got the same error.


text:
How to match peptide sequence from two different condition
I have to peptide output sample output i would give here ..
Sample A
dput(head(B5_ILS_Mus_PeptideSummary))
structure(list(N = c(1, 1, 1, 1, 1, 1), Unused = c(47.67, 47.67,
47.67, 47.67, 47.67, 47.67), Total = c(47.67, 47.67, 47.67, 47.67,
47.67, 47.67), `%Cov` = c(98.0000019073486, 98.0000019073486,
98.0000019073486, 98.0000019073486, 98.0000019073486, 98.0000019073486
), `%Cov(50)` = c(86.0000014305115, 86.0000014305115, 86.0000014305115,
86.0000014305115, 86.0000014305115, 86.0000014305115), `%Cov(95)` = c(86.0000014305115,
86.0000014305115, 86.0000014305115, 86.0000014305115, 86.0000014305115,
86.0000014305115), Accessions = c("sp|P20065|TYB4_MOUSE", "sp|P20065|TYB4_MOUSE",
"sp|P20065|TYB4_MOUSE", "sp|P20065|TYB4_MOUSE", "sp|P20065|TYB4_MOUSE",
"sp|P20065|TYB4_MOUSE"), Names = c("Thymosin beta-4 OS=Mus musculus OX=10090 GN=Tmsb4x PE=1 SV=1",
"Thymosin beta-4 OS=Mus musculus OX=10090 GN=Tmsb4x PE=1 SV=1",
"Thymosin beta-4 OS=Mus musculus OX=10090 GN=Tmsb4x PE=1 SV=1",
"Thymosin beta-4 OS=Mus musculus OX=10090 GN=Tmsb4x PE=1 SV=1",
"Thymosin beta-4 OS=Mus musculus OX=10090 GN=Tmsb4x PE=1 SV=1",
"Thymosin beta-4 OS=Mus musculus OX=10090 GN=Tmsb4x PE=1 SV=1"
), Used = c(NA, NA, NA, NA, NA, NA), Annotation = c(NA, NA, NA,
NA, NA, NA), Contrib = c(2, 2, 2, 2, 2, 2), Conf = c(99.0000009536743,
58.2899987697601, 97.5199997425079, 99.0000009536743, 98.9300012588501,
73.7299978733063), Sequence = c("EKNPLPSKETIEQEKQAGES", "KKTETQEKNPLPSKETIEQEKQAGES",
"KNPLPSKETIEQEKQAGES", "KTETQEKNPLPSKETIEQEKQAGES", "LKKTETQEKNPLPSKETIEQEKQAGES",
"PLPSKETIEQEKQAGES"), Modifications = c(NA_character_, NA_character_,
NA_character_, NA_character_, NA_character_, NA_character_),
    ProteinModifications = c(NA_character_, NA_character_, NA_character_,
    NA_character_, NA_character_, NA_character_), Cleavages = c(NA,
    NA, NA, NA, NA, NA), dMass = c(-0.0625222995877266, 0.180465996265411,
    -0.00865099020302296, -0.0651467964053154, -0.107698999345303,
    -0.0186963994055986), `Obs MW` = c(2241.05004882813, 2956.67944335938,
    2112.06127929688, 2828.3388671875, 3069.47534179688, 1869.91320800781
    ), `Obs m/z` = c(748.0239, 740.1771, 529.0226, 708.092, 768.3761,
    624.3117), `Theor MW` = c(2241.11254882813, 2956.4990234375,
    2112.06982421875, 2828.40405273438, 3069.5830078125, 1869.93200683594
    ), `Theor m/z` = c(748.044738769531, 740.132019042969, 529.024719238281,
    708.108276367188, 768.403015136719, 624.317932128906), `Theor z` = c(3,
    4, 4, 4, 4, 3), Sc = c(15, 8, 11, 21, 12, 10), Spectrum = c("1.1.1.372.12",
    "1.1.1.339.6", "1.1.1.358.3", "1.1.1.355.16", "1.1.1.346.14",
    "1.1.1.384.8"), `Acq Time` = c(4.2116, 3.539133, 3.91795,
    3.867617, 3.68345, 4.45125), `Intensity (Peptide)` = c(3775.35,
    110.58, 4374.93, 7285.52, 3073.51, 5776.43), PrecursorIntensityAcquisition = c(3626.33,
    37.12, 4301.05, 7124.23, 2706.13, 5533.12), `Apex Time (Peptide)` = c(4.2,
    3.59, 3.91, 3.85, 3.67, 4.44), `Elution Peak Width (Peptide)` = c(0.1,
    0.06, 0.08, 0.1, 0.08, 0.08), MS2Counts = c(3077.868, 477.5326,
    2262.827, 4694.357, 1759.267, 2618.842)), row.names = c(NA,
-6L), class = c("tbl_df", "tbl", "data.frame"))

Sample B
dput(head(B6_ILS_mus_PeptideSummary))
structure(list(N = c(1, 1, 1, 1, 1, 1), Unused = c(82.41, 82.41,
82.41, 82.41, 82.41, 82.41), Total = c(82.41, 82.41, 82.41, 82.41,
82.41, 82.41), `%Cov` = c(100, 100, 100, 100, 100, 100), `%Cov(50)` = c(87.9999995231628,
87.9999995231628, 87.9999995231628, 87.9999995231628, 87.9999995231628,
87.9999995231628), `%Cov(95)` = c(86.0000014305115, 86.0000014305115,
86.0000014305115, 86.0000014305115, 86.0000014305115, 86.0000014305115
), Accessions = c("sp|P20065|TYB4_MOUSE", "sp|P20065|TYB4_MOUSE",
"sp|P20065|TYB4_MOUSE", "sp|P20065|TYB4_MOUSE", "sp|P20065|TYB4_MOUSE",
"sp|P20065|TYB4_MOUSE"), Names = c("Thymosin beta-4 OS=Mus musculus OX=10090 GN=Tmsb4x PE=1 SV=1",
"Thymosin beta-4 OS=Mus musculus OX=10090 GN=Tmsb4x PE=1 SV=1",
"Thymosin beta-4 OS=Mus musculus OX=10090 GN=Tmsb4x PE=1 SV=1",
"Thymosin beta-4 OS=Mus musculus OX=10090 GN=Tmsb4x PE=1 SV=1",
"Thymosin beta-4 OS=Mus musculus OX=10090 GN=Tmsb4x PE=1 SV=1",
"Thymosin beta-4 OS=Mus musculus OX=10090 GN=Tmsb4x PE=1 SV=1"
), Used = c(NA, NA, NA, NA, NA, NA), Annotation = c(NA, NA, NA,
NA, NA, NA), Contrib = c(2, 2, 2, 2, 2, 2), Conf = c(99.0000009536743,
99.0000009536743, 99.0000009536743, 99.0000009536743, 94.4599986076355,
99.0000009536743), Sequence = c("DKPDMAEIEKFDKSKLK", "EKNPLPSKETIEQEKQAGES",
"ETQEKNPLPSKETIEQEKQAGES", "IEKFDKSKLKKTETQEKNPLPSKETIEQEKQAGES",
"KKTETQEKNPLPSKETIEQEKQAGES", "KNPLPSKETIEQEK"), Modifications = c(NA_character_,
NA_character_, NA_character_, NA_character_, NA_character_, NA_character_
), ProteinModifications = c(NA_character_, NA_character_, NA_character_,
NA_character_, NA_character_, NA_character_), Cleavages = c(NA,
NA, NA, NA, NA, NA), dMass = c(0.0129736997187138, -0.0542825013399124,
-0.0299271997064352, -0.107947997748852, -0.0191929005086422,
-0.0046052997931838), `Obs MW` = c(2021.06335449219, 2241.05834960938,
2599.2314453125, 4045.01391601563, 2956.47973632813, 1639.87341308594
), `Obs m/z` = c(506.2731, 748.0267, 650.8151, 810.0101, 592.3032,
547.6318), `Theor MW` = c(2021.05029296875, 2241.11254882813,
2599.26123046875, 4045.12182617188, 2956.4990234375, 1639.87805175781
), `Theor m/z` = c(506.269866943359, 748.044738769531, 650.822631835938,
810.031677246094, 592.307067871094, 547.63330078125), `Theor z` = c(4,
3, 4, 5, 5, 3), Sc = c(11, 12, 12, 26, 10, 10), Spectrum = c("1.1.1.858.2",
"1.1.1.621.5", "1.1.1.640.9", "1.1.1.643.4", "1.1.1.597.3", "1.1.1.604.4"
), `Acq Time` = c(8.74735, 4.069217, 4.43505, 4.4987, 3.615083,
3.7369), `Intensity (Peptide)` = c(1798.16, 2856.36, 1278.57,
5833.81, 758.71, 599.21), PrecursorIntensityAcquisition = c(1603.78,
443.55, 1104.72, 5265.07, 484.03, 423.02), `Apex Time (Peptide)` = c(8.72,
4.19, 4.39, 4.45, 3.55, 3.78), `Elution Peak Width (Peptide)` = c(0.14,
0.12, 0.12, 0.37, 0.1, 0.08), MS2Counts = c(456.3906, 313.663,
712.3438, 2748.186, 351.0362, 228.9387)), row.names = c(NA, -6L
), class = c("tbl_df", "tbl", "data.frame"))

Unlike genes expression where it can be intersected since i know the gene name or ensemble ID are consistent.
Here My objective is how many of the sequence are common or unique to these two condition.
In the accession column many ID are same but the differences are in the sequence column.
Is there any programmatic way of finding the elements that are common or different in both the conditions

merge(a, b, by = "Sequence") would "merge" the two tables (as the name suggests) based on the Sequence column. The caveat is that sequences have to be exactly the same to be considered "common to both samples", so even 1 aa difference would be missed.
Simplified output:
1       EKNPLPSKETIEQEKQAGES   1    47.67   47.67     98         86         86 sp|P20065|TYB4_MOUSE
2 KKTETQEKNPLPSKETIEQEKQAGES   1    47.67   47.67     98         86         86 sp|P20065|TYB4_MOUSE


text:
find conservative regions relative to a specific accession
I have 7 RefSeq accessions which I performed multiple sequence aligment on it using CLUSTALO. These are the following:
NC_026431.1
NC_007363.1
NC_004907.1
NC_007377.1
NC_007367.1
NC_002016.1
NC_026427.1

I input the fasta sequence of those 7 accessions on CLUSTALO under this link. I get the result in either the msf format or clustal format with characters. The result looks like the following:
CLUSTAL O(1.2.4) multiple sequence alignment

NC_002016.1      -------AGCGAAAGCAGGTAGATATTGAAAGATGAGTCTTCTAACCGAGGTCGAAACGT   53
NC_007377.1      -------AGCAAAAGCAGGTAGATATTGAAAGATGAGCCTTCTAACCGAGGTCGAAACGT   53
NC_007367.1      -------AGCAAAAGCAGGTAGATATTGAAAGATGAGCCTTCTAACCGAGGTCGAAACGT   53
NC_004907.1      GGGGAATTCCAAAAGCAGGTAGATATTGAAAGATGAGTCTTCTAACCGAGGTCGAAACGT   60
NC_026427.1      --------------------------------ATGAGTCTTCTAACCGAGGTCGAAACGT   28
NC_007363.1      -------AGCAAAAGCAGGTAGATATTGAAAAATGAGTCTTCTAACCGAGGTCGAAACGT   53
NC_026431.1      --------------------------------ATGAGTCTTCTAACCGAGGTCGAAACGT   28
                                                 ***** **********************

I'm aware that the * symbols means a highly conserved region and that's exatcly what I'm looking for , however I need the highly conserved regions relative to a specific accession, namely NC_007367.1. I'm not sure how to do that. I'm also not sure what the - represents.

https://www.ddbj.nig.ac.jp/faq/en/explain-three-symbols-e.html, https://abacus.bates.edu/bioinformatics1/clustalw2.html. "highly conserved" is probably not an accurate way to judge "*". All it means is that the position aligns perfectly across the input sequences, there is no statement about conservation or other evolutionary phenomena. clustal is only an alignment program.

thanks, what would be the tool to find the conservative region then ? since CLUSTAL just gives a similarity estimation.

THere are various tools- depends a bit what your goal is in finding these positions. Simple visualizations that emphasize relative stability in alignments, e.g. here: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3226196/. Gblocks is a somewhat more disciplined statistical approach, I am still unsure if it is exactly what you want: https://academic.oup.com/mbe/article/17/4/540/1127654, webserver here: http://molevol.cmima.csic.es/castresana/Gblocks_server.html

You can opt for a similarity matrix as output (here is an example). The fact that you would like to judge conservation "relative to a given sequence", you can use this matrix to rank similarities of different sequences with respect to any of these sequences.
I would advise visualizing these similarities, for example in terms of a heatmap, and seeing if anything stands out (with respect to your sequence of interest). And I would advise against selecting an arbitrary threshold (i.e. 90% or 95%), this would trigger questions based on the selection criteria of such a cutoff and this is understandable, such arbitrary cut-offs often do not mean much in biology.