text:
curl: (22) The requested URL returned error: 400 Bad Request when using entrez direct
I've just dowloaded entrez direct on my Mac and I was trying to perform some simple queries like:
esearch -db pubmed -query "tn3 transposition immunity" | efetch -format medline

but I keep getting this error
curl: (22) The requested URL returned error: 400 Bad Request
 ERROR:  curl command failed ( Mar 19 Ott 2021 13:02:26 CEST ) with: 22
https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi -d db=pubmed&id=Unable%2Cto%2Clocate%2Cxtract%2Cexecutable.%2CPlease%2Cexecute%2Cthe%2Cfollowing%2Cnquire%2Cdwn%2Cftp.ncbi.nlm.nih.gov%2Centrez%2Centrezdirect%2Cxtract.Darwin.gz%2Cgunzip%2Cf%2Cxtract.Darwin.gz%2Cchmod%2Cx%2Cxtract.Darwin&rettype=medline&retmode=text&tool=edirect&edirect=16.1&edirect_os=Darwin&email=mariocatalano%40MBP-di-Mario.lan"

I'm using bash scripts on my Mac's terminal
How should I fix this?

update your curl or esearch version, use the miniconda/brew versions. If you already have conda, you can create an environment and install it on it.

Looks like you haven't installed it correctly. The error is saying it cannot locate/extract the executable.

@zorbax -- it is recommended to use the official version of Edirect instead of the  ones from miniconda, brew, etc as they often tend to be quite behind the latest, official version.

It's not a curl problem.  (The 22 error is just passing along the server-side HTTP 400 error, which is because it dumped an error message into a URL.)  I think terdon is right and there's something wrong with your edirect installation.  You have esearch and efetch on your PATH but apparently not xtract and/or xtract.Darwin.

More detail:
Look toward the end of the error message, at the text being passed to curl's -d argument:
db=pubmed&id=Unable%2Cto%2Clocate%2Cxtract%2Cexecutable.%2CPlease%2Cexecute%2Cthe%2Cfollowing%2Cnquire%2Cdwn%2Cftp.ncbi.nlm.nih.gov%2Centrez%2Centrezdirect%2Cxtract.Darwin.gz%2Cgunzip%2Cf%2Cxtract.Darwin.gz%2Cchmod%2Cx%2Cxtract.Darwin&rettype=medline&retmode=text&tool=edirect&edirect=16.1&edirect_os=Darwin&email=mariocatalano%40MBP-di-Mario.lan"

It's putting this in the id= part:
Unable%2Cto%2Clocate%2Cxtract%2Cexecutable.%2CPlease%2Cexecute%2Cthe%2Cfollowing%2Cnquire%2Cdwn%2Cftp.ncbi.nlm.nih.gov%2Centrez%2Centrezdirect%2Cxtract.Darwin.gz%2Cgunzip%2Cf%2Cxtract.Darwin.gz%2Cchmod%2Cx%2Cxtract.Darwin

which is maybe something like (trying to guess which characters got turned into %2C):

Unable to locate xtract executable. Please execute the following: nquire -dwn ftp.ncbi.nlm.nih.gov entrez/entrezdirect xtract.Darwin.gz; gunzip -f xtract.Darwin.gz; chmod +x xtract.Darwin

But you should already have xtract and xtract.Darwin from the same edirect install that gave you esearch and efetch.
For what it's worth the official Entrez Direct install on Mac OS looks a little goofy.  When I try it myself it ends with:
Entrez Direct has been successfully downloaded and installed.

-e In order to complete the configuration process, please execute the following:

  echo "export PATH=\\\${PATH}:" >> $HOME/.bashrc
  echo "export PATH=\\\${PATH}:" >> $HOME/.zshrc

or manually edit the PATH variable assignments in your .bashrc and .zshrc files.

Would you like to do that automatically now? [y/N]

Those commands don't look right (what's getting added to the path? an empty string?) and if I say y, I get this at the end of .bashrc and .zshrc which definitely isn't right:
export PATH=${PATH}:

It seems like an empty string at the end of PATH behaves like . (current directory), so I can go into ~/edirect and run your command and it does actually work.  So this isn't exactly your problem... but still, that seems wrong.

Note that macOS uses login shells by default for historical reasons, which means that if the OP is using bash (which is no longer the default on macOS), then ~/.bashrc is irrelevant and they should use ~/.profile instead. Either way though, that command is quite wrong, as you point out.

This was an issue with the Edirect installation script. The install-edirect.sh script has been updated to fix this issue. Please follow the instructions here to install.

Thanks!  The install still seems to misbehave on Mac OS, though.  When I say ymy ~/.bashrc is now untouched and my ~/.zshrc contains literallyexport PATH=${PATH}:/Users/jesse echo export PATH=${PATH}:/Users/jesse

@Jesse thank you for trying it out. This is yet another issue that needed a quick fix. The updated script is already on FTP so if you try again it should work. Sorry for the inconvenience.

No trouble, thanks for working on this.  Still behaving strangely for me but when I look at the script I think I see one reason.  It references DIR when setting up PATH, but that's the original working directory / the directory containing the install script, not the destination directory, right?  Should that instead be a $PWD or something once inside the ~/edirect location, maybe? We might be drifting from OP's question but I'm happy to help debug this via email or a github gist or something.

@Jesse, turns out this is yet another issue that had to be fixed. The new installation script was derived from the old setup.sh file, which was run within the edirect folder. In that context the DIR variable was correct. The new code now sets DIR after calling cd edirect, so its value is now also correct. Sorry for the inconvenience, and thanks again for your reports.

Yep, works as expected now on Mac OS.  Thanks!


text:
Length of Contigs in Transcriptome and Whole Genome Assembly
Why are there shorter contigs from transcriptome assembly than from a whole genome assembly? I know the difference between transcriptome and genome, but don't really understand what contigs are in the context of sequencing in bioinformatics, and why are the contigs shorter in transcriptome assembly than in whole genome assembly. Anyone mind explaining? Much appreciated.

The longest possible transcriptome contig reflects the longest possible transcript. In the human genome, that's possibly Titin, which is ~35K aa (that's ~105Kb) and will be longer pre-splicing.
The longest possible genome contig reflects the longest chromosome. In the human genome that's chr1 at ~250Mb.
In practise your contigs/scaffolds will be a lot shorter for various reasons (e.g. hard to sequence the telomeres and repeats) but it illustrates the difference since all the chromosomes consist of contiguous DNA but many separate RNA transcripts (and not all DNA is transcribed).

source: Stack Exchange

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

sort: votes


text:
How to normalise indel of different transcript direction?
I am trying to reproduce an analysis result, and there is a requirement regarding the format of the result: the '3' rule ': all mutations should be expressed in a position close to the end of the gene transcription direction (3' end).
The result I obtained is(annotate from VEP)
Chromosome  Start_Position  End_Position    Strand  Variant_Classification  Variant_Type    Reference_Allele    Tumor_Seq_Allele1   Tumor_Seq_Allele2   HGVSc   HGVSp   HGVSp_Short Transcript_ID
10  89692835    89692840    +   In_Frame_Del    DEL GATCTT  GATCTT  -   c.321_326del    p.Leu108_Asp109del  p.L108_D109del  NM_000314.8

the standard answer
chr10   89692837    89692842    TCTTGA      PTEN    NM_00314    c.321_326del    p.Leu108_Asp109del

I speculate that ：
the HGVS standard followed by VEP includes the 3'rule . So ultimately, there is the same p.Leu108_ Asp109del。 But the output of vep did not adjust Start_ Position and End_ Position？（Perhaps it thinks that after obtaining the correct amino acid mutation , it doesn't matter anymore）
But for me, I hope the ultimate Start_ Position and End_ The Position also complies with the regulations, which seems a bit troublesome because the different  transcripts may have different directions and always rely on the IGV to correct it , which is also tiring.
Thanks for any suggestions

I have added an IGV image to better illustrate the problem. Reference genome regions 89692830 to 89692845 for the example in the figure.
bwa(or other software，like bowtie2) believes that the position of del should be
ref:GTGAAGATCTTGACCA
alt:GTGAA------GACCA

But 3'rule believes that the position of del should be
ref:GTGAAGATCTTGACCA
alt:GTGAAGA------CCA

To verify whether the results match 3'rule, it is necessary to know the direction of the transcripts, and each gene always has multiple transcripts, which is where I am confused

Thanks for the update @tanqiang

Can you clarify which representation you need? Is it the genomic level, the coding level or the protein level? Also, what is the input here? Do you start with variants or with VEP output? How many do you have? Do you only want to get one result, or does this need to be automated? Why are other transcripts relevant if you already have NM_000314.8?

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


text:
Eigen phred score interpretation
I need help. I am not able to find any reference for Eigen_phred_coding prediction tools score from dbNSFP.
I know that it may be scaled in a similar manner to CADD_phred - ranking raw values from highest to smallest; get rank of that value in the queue; phred scale that value with the same formula as sequencing quality score.
But I couldn't find the phred scaling of Eigen values directly, not to mention its suggested cutoff (same as CADD_phred > 15 or different?).
How are the Eigen scores transformed into Eigen_phred?
Is it in same manner as Cadd_phred?
If so, can same cutoff be used in both cases?

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

@terdon Thank you for that reference. I don't know why, but I always thought that CADD_phred is phred scaled value of ratio rank of value in sorted queue over the total number of used SNV. At least in read me I found out that everything phred is  ratio of rank of the variat over total number of scores in database. Still I am not sure with the cutoff of 15. Anyway I am going to try that.


text:
How can I get a list of neighboring Hydrogen atoms of an alpha-carbon?
Suppose, I want to find the neighboring Hydrogen atoms of a specific C-alpha atom.
The following was my effort. However, it produces a bunch of empty lists.
from typing import List
from Bio.PDB import Chain, Atom, NeighborSearch, PDBParser

def get_hydrogen_atoms(c_alpha_atom):
    ns = NeighborSearch([c_alpha_atom])
    return ns.search_all(190.0, "A")

p = PDBParser()
structure = p.get_structure('16pk', '16pk.pdb')

c_alpha_list = [atom for atom in structure.get_atoms()]

for c_alpha in c_alpha_list:
    if c_alpha.name == "CA":
        hydrogen_list = get_hydrogen_atoms(c_alpha)
        print(hydrogen_list)

How can I correct this?

hydrogen_list = get_hydrogen_atoms(c_alpha) should be hydrogen_list.append(get_hydrogen_atoms(c_alpha)).
Also your structure needs to be protonated —most crystal structures are not. radius is in Å, not pm. Cα atoms are no hydrogen bond donors/acceptors. Why no search for nearby N atomnames to the atomname O, which are donor - acceptor pairs...

You posted similar question here: https://stackoverflow.com/questions/69676943/how-can-i-get-a-list-of-neighboring-hydrogen-atoms-of-an-alpha-carbon/69815509#69815509
ok got some code :
from Bio.PDB import Chain, Atom, NeighborSearch, PDBParser, Selection

def get_hydrogen_atoms(c_alpha_atom: Atom.Atom, po):
    ns = NeighborSearch(po)
    return [h for h in ns.search( c_alpha_atom.get_coord(), 10.0, 'A') if h.element == 'H']

p = PDBParser()
# structure = p.get_structure('16pk', '16pk.pdb')
structure = p.get_structure('16pk', 'pH8_5_1bq9_short.pdb')

c_alpha_list = [atom for atom in structure.get_atoms() if atom.name == "CA"]

print('c_alpha_list[0]  : ', c_alpha_list[0], c_alpha_list[0].id , c_alpha_list[0].get_serial_number(),'\n')

po = Selection.unfold_entities(structure, 'A')

hydrogen_list = get_hydrogen_atoms(c_alpha_list[0], po)

print('hydrogen_list : \n',hydrogen_list)

print('\nSelection (po) : \n', po)

had to use a different PDB file because I needed the H atoms,
pH8_5_1bq9_short.pdb :
HETATM    1  N   FME A   1      23.447  -5.765   5.119  1.00 19.07           N
HETATM    2  CN  FME A   1      23.308  -6.776   4.291  1.00 20.46           C
HETATM    3  O1  FME A   1      22.504  -6.713   3.357  1.00 24.75           O
HETATM    4  CA  FME A   1      22.488  -4.745   5.410  1.00 18.15           C
HETATM    5  CB  FME A   1      22.433  -4.342   6.827  1.00 21.15           C
HETATM    6  CG  FME A   1      22.198  -5.618   7.751  1.00 22.84           C
HETATM    7  SD  FME A   1      20.445  -6.032   7.772  1.00 23.40           S
HETATM    8  CE  FME A   1      20.476  -7.587   8.656  1.00 28.76           C
HETATM    9  C   FME A   1      22.699  -3.530   4.487  1.00 17.12           C
HETATM   10  O   FME A   1      23.804  -3.248   4.077  1.00 20.25           O
HETATM    0 3HE  FME A   1      19.574  -7.938   8.726  1.00 28.76           H   new
HETATM    0 2HG  FME A   1      22.713  -6.369   7.418  1.00 22.84           H   new
HETATM    0 2HE  FME A   1      21.032  -8.221   8.177  1.00 28.76           H   new
HETATM    0 2HB  FME A   1      23.260  -3.901   7.077  1.00 21.15           H   new
HETATM    0 1HG  FME A   1      22.508  -5.438   8.652  1.00 22.84           H   new
HETATM    0 1HE  FME A   1      20.838  -7.447   9.545  1.00 28.76           H   new
HETATM    0 1HB  FME A   1      21.717  -3.701   6.961  1.00 21.15           H   new
HETATM    0  HN  FME A   1      24.199  -5.711   5.532  1.00 19.07           H   new
HETATM    0  HCN FME A   1      23.841  -7.577   4.415  1.00 20.46           H   new
HETATM    0  HA  FME A   1      21.619  -5.138   5.231  1.00 18.15           H   new
ATOM     11  N   ALA A   2      21.633  -2.939   4.043  1.00 15.50           N
ATOM     12  CA  ALA A   2      21.644  -1.884   3.062  1.00 15.08           C
ATOM     13  C   ALA A   2      20.899  -0.661   3.519  1.00 13.72           C
ATOM     14  O   ALA A   2      20.278  -0.702   4.584  1.00 13.49           O
ATOM     15  CB  ALA A   2      21.083  -2.339   1.710  1.00 17.86           C
ATOM      0  H   ALA A   2      20.842  -3.144   4.312  1.00 15.50           H   new
ATOM      0  HA  ALA A   2      22.579  -1.652   2.951  1.00 15.08           H   new
ATOM      0  HB1 ALA A   2      21.110  -1.599   1.083  1.00 17.86           H   new
ATOM      0  HB2 ALA A   2      21.618  -3.072   1.368  1.00 17.86           H   new
ATOM      0  HB3 ALA A   2      20.166  -2.634   1.823  1.00 17.86           H   new
ATOM     16  N   LYS A   3      21.005   0.412   2.774  1.00 12.48           N
ATOM     17  CA  LYS A   3      20.341   1.658   3.029  1.00 12.61           C
ATOM     18  C   LYS A   3      19.415   1.991   1.879  1.00 11.53           C
ATOM     19  O   LYS A   3      19.784   1.753   0.701  1.00 11.66           O
ATOM     20  CB  LYS A   3      21.338   2.806   3.273  1.00 13.38           C
ATOM     21  CG  LYS A   3      22.062   2.721   4.593  1.00 16.70           C
ATOM     22  CD  LYS A   3      23.062   3.847   4.814  1.00 18.05           C
ATOM     23  CE  LYS A   3      23.806   3.576   6.136  1.00 20.42           C
ATOM     24  NZ  LYS A   3      24.878   4.503   6.403  1.00 22.91           N
ATOM      0  H   LYS A   3      21.494   0.433   2.067  1.00 12.48           H   new
ATOM      0  HA  LYS A   3      19.821   1.557   3.842  1.00 12.61           H   new
ATOM      0  HB2 LYS A   3      21.992   2.812   2.556  1.00 13.38           H   new
ATOM      0  HB3 LYS A   3      20.862   3.650   3.230  1.00 13.38           H   new
ATOM      0  HG2 LYS A   3      21.411   2.732   5.312  1.00 16.70           H   new
ATOM      0  HG3 LYS A   3      22.527   1.871   4.644  1.00 16.70           H   new
ATOM      0  HD2 LYS A   3      23.689   3.892   4.076  1.00 18.05           H   new
ATOM      0  HD3 LYS A   3      22.606   4.702   4.851  1.00 18.05           H   new
ATOM      0  HE2 LYS A   3      23.170   3.611   6.868  1.00 20.42           H   new
ATOM      0  HE3 LYS A   3      24.166   2.676   6.116  1.00 20.42           H   new
ATOM      0  HZ1 LYS A   3      25.266   4.294   7.176  1.00 22.91           H   new
ATOM      0  HZ2 LYS A   3      25.479   4.458   5.748  1.00 22.91           H   new
ATOM      0  HZ3 LYS A   3      24.551   5.330   6.449  1.00 22.91           H   new
ATOM     25  N   TRP A   4      18.274   2.598   2.165  1.00 10.75           N
ATOM     26  CA  TRP A   4      17.283   3.019   1.181  1.00 10.69           C
ATOM     27  C   TRP A   4      16.875   4.472   1.502  1.00 11.37           C
ATOM     28  O   TRP A   4      16.812   4.872   2.651  1.00 17.19           O
ATOM     29  CB  TRP A   4      16.012   2.143   1.352  1.00 11.20           C
ATOM     30  CG  TRP A   4      16.164   0.734   0.997  1.00 11.80           C
ATOM     31  CD1 TRP A   4      16.913  -0.218   1.669  1.00 12.50           C
ATOM     32  CD2 TRP A   4      15.535   0.036  -0.079  1.00 11.53           C
ATOM     33  NE1 TRP A   4      16.806  -1.435   1.050  1.00 12.76           N
ATOM     34  CE2 TRP A   4      15.959  -1.302  -0.013  1.00 12.61           C
ATOM     35  CE3 TRP A   4      14.642   0.421  -1.079  1.00 11.88           C
ATOM     36  CZ2 TRP A   4      15.531  -2.258  -0.929  1.00 13.18           C
ATOM     37  CZ3 TRP A   4      14.217  -0.527  -1.989  1.00 12.54           C
ATOM     38  CH2 TRP A   4      14.668  -1.858  -1.894  1.00 14.26           C

output :
c_alpha_list[0]  :  <Atom CA> CA 4

hydrogen_list :
 [<Atom HB1>, <Atom HB3>, <Atom H>, <Atom 2HE>, <Atom 3HE>, <Atom 1HE>, <Atom HB3>, <Atom HA>, <Atom HG2>, <Atom HB2>, <Atom HCN>, <Atom HA>, <Atom HN>, <Atom 1HB>, <Atom 2HB>, <Atom 2HG>, <Atom 1HG>, <Atom H>, <Atom HB2>, <Atom HA>, <Atom HD2>, <Atom HG3>, <Atom HD3>, <Atom HZ2>, <Atom HE3>, <Atom HE2>, <Atom HZ1>]

Selection (po) :
 [<Atom N>, <Atom CN>, <Atom O1>, <Atom CA>, <Atom CB>, <Atom CG>, <Atom SD>, <Atom CE>, <Atom C>, <Atom O>, <Atom 3HE>, <Atom 2HG>, <Atom 2HE>, <Atom 2HB>, <Atom 1HG>, <Atom 1HE>, <Atom 1HB>, <Atom HN>, <Atom HCN>, <Atom HA>, <Atom N>, <Atom CA>, <Atom C>, <Atom O>, <Atom CB>, <Atom H>, <Atom HA>, <Atom HB1>, <Atom HB2>, <Atom HB3>, <Atom N>, <Atom CA>, <Atom C>, <Atom O>, <Atom CB>, <Atom CG>, <Atom CD>, <Atom CE>, <Atom NZ>, <Atom H>, <Atom HA>, <Atom HB2>, <Atom HB3>, <Atom HG2>, <Atom HG3>, <Atom HD2>, <Atom HD3>, <Atom HE2>, <Atom HE3>, <Atom HZ1>, <Atom HZ2>, <Atom HZ3>, <Atom N>, <Atom CA>, <Atom C>, <Atom O>, <Atom CB>, <Atom CG>, <Atom CD1>, <Atom CD2>, <Atom NE1>, <Atom CE2>, <Atom CE3>, <Atom CZ2>, <Atom CZ3>, <Atom CH2>]

the selection list contains both ATOM and HETATM
have a look at : https://www.biostars.org/p/269579/ How to use biopython neighborsearch?
and Bio.PDB.NeighborSearch module :https://biopython.org/docs/1.75/api/Bio.PDB.NeighborSearch.html ;
Bio.PDB.Atom module :https://biopython.org/docs/1.75/api/Bio.PDB.Atom.html


text:
Is there an alternative to bulked segregant analysis for insects?
This strategy seems to be most commonly used for plants. When crossing animals isn't possible, how can I do a similar study to identify a particular locus responsible for a polymorphic trait (specifically, the pigmentation of a beetle)?

I'll have a go at answering this. I don't understand beetles within Insecta: I am well aware they very distinct, obviously in phenotype and ecology, but also genetically. I think they have unique compositional bias even amongst holometabolous insects.
Personally I think if you've a recent genetic trait, e.g. insecticide resistance, that evolved within the last 30 years or less, you're dealing with pretty homogeneous backgrounds between traits. In this example trait, bulk segregant analysis might work.
Beetles are really ancient. I suspect their families are particularly ancient, because they are just so successful and I'd guess Scarabaeidae are ancient too, i.e. the family scarab beetles. This family could easily be 50 million years old, or older. Each species therein could a few million years old. For insects thats normal and these are beetles.
The concern here is over million+ years of evolution (I don't know the study size), they will comprise loads of neutral signal that would act as noise, i.e. there will be excessive third codon wobble which will mask SNPs/genes causing colour (color) variation.
My approach, given I know nothing about colour (color) variation in insects other than "white-eyed" phenotype in Drosophila melanogaster,

assess genetic complexity of colour variation in other insect Orders;
get hold of as many proteins involved in colour in insects as possible;
particularly get hold of colour variation proteins in beetles;
sequence the genome of a representative scarab beetle (which I think your doing)
perform psi-blast using the hypotheses from 1., 2. and particularly 3. to assess which proteins in the scarab beetle might be responsible.
perform deep sequencing on those targets, because they could be multigene families *.

The next step is phylogenetics to assess convergent and/or parallel evolution **. That in my view is a separate question because its quite complicated and you'd need to write your own pipeline to automate it (you could do it by hand however).
What I would be cautious about is that beetles are not any other insect so inter-Order homology needs to be treated with caution.

*, Just to point out, there are lots sequencing strategies its not my thing, e.g. low coverage draft genomes. PCR deep sequencing for a given gene target, without any draft genomes and isn't conserved e.g. rRNA types genes, will need small PCR fragments because of the neutral mutations which will disrupt primer specificity. If you're good at sequencing a draft genome could be cool (just guessing).
**, A phylogenetics approach would dovetail with the COI and rRNA gene sequencing (28S?) already underway, because these genes should approximate to the species phylogeny. What is sought are residues  within the colour genes (or perhaps presence/absence genes) that will cluster not according to their evolutionary descent but according to their colours (colors).


text:
Off-target % for whole-exome sequencing panel
My samples have been sequenced with the Twist Exome 2.0 sequencing panel, and I want to assess how efficient the sequencing has been. By efficient, I mean examining metrics such as the uniformity (fold-80 score), coverage, on/off-target % etc. as documented in their webpage. I checked the % off-target as follows:

Count the # of reads in the bam file (samtools view -c "$dir/$base.bam)

Count the # of reads that intersect with the bed region (bedtools intersect -a "$dir/$base.bam" -b $bed -u | samtools view -c -)

Calculate the % off-target rate as (bc -l <<< "(1 - ($bed_reads / $total_reads)) * 100")

I have two questions:

Is there a cleaner/simpler way to calculate the off-target rate? I know I can calculate the on-target rate with bedtools coverage, but I am not sure of how to calculate the off-target rate with a single command line.
I obtained an off-target rate of approximately 20-25% for my samples. Is this within the "good" off-target % range?


text:
Read PDB file, extract dihedral angles, modify dihedral angles, reconstruct Cartesian coordinates, and write PDB file
As the title summarizes, I am trying to:

Read a PDB file (for example, 1enh.pdb).
Extract the backbone dihedral angles (phi, psi).
Modify the dihedral angles (phi, psi) arbitrarily (for example, replace them with np.random.uniform(-np.pi, np.pi, size = (56, 2)) in 1enh, since there are 56 pairs of dihedrals).
Reconstruct the Cartesian coordinates of the modified backbone.
Write a PDB file based on the modified backbone.

From this question it seems the tools for what I want are available in BioPython. There is also a set_angle() method. But I only have a very limited knowledge of Python and BioPython. I am stuck at Steps 3 and 4, and dubious on Step 5.
This is what I have tried so far:
#!/usr/bin/python

import numpy as np
from Bio.PDB import PDBParser, PICIO, PDBIO

## Step 1

parser = PDBParser()
structure = parser.get_structure("", "1enh.pdb")

## Step 2

structure.atom_to_internal_coordinates()
chain = list(structure.get_chains())[0]
ic_chain = chain.internal_coord
dihedrals = ic_chain.dihedra

## Step 3 - ??

# Modifiy (phi,psi). For example, replace (phi,psi) with
newdihedrals = np.random.uniform(-np.pi, np.pi, size = (56, 2))

## Step 4 - ??

modified = PICIO.read_PIC('internal')
modified.internal_to_atom_coordinates()

## Step 5

io = PDBIO()
io.set_structure(modified)
io.save('1enh_modified.pdb', preserve_atom_numbering = True)

Many thanks in advance.
EDIT: I do not care about possible clashes if the side-chain is included in the modified PDB. In any case, writing a modified PDB file that only includes the backbone is also fine.

POINT 5 READS : Write a PDB file based on the modified backbone. --> Just C-alpha or side-chains too ?? How to avoid clashes in case of the latter ?

Thanks for the comment. C-alpha would suffice. If side-chains are included, I do not really care about clashes, as I could determine them later and reject the modified structure.

Any Python package that checks PDB geometry  just for clashes ??

In PyMOL I remember it is possible to detect clashes, I do not know other functionality beyond it. In any case, my main interest is on tweaking the (phi,psi) angles and reconstruct the protein backbone in a PDB file. Side-chains could be present or not. Thanks.

ok, here my input pdb, example_short.pdb  :
ATOM      1  N   MET A   1      14.067  31.664  -1.639  1.00 80.87           N
ATOM      2  HN1 MET A   1      14.250  31.259  -0.699  1.00  0.00           H
ATOM      3  HN2 MET A   1      14.891  32.220  -1.945  1.00  0.00           H
ATOM      4  HN3 MET A   1      13.229  32.279  -1.596  1.00  0.00           H
ATOM      5  CA  MET A   1      13.824  30.564  -2.612  1.00 79.17           C
ATOM      6  HA  MET A   1      14.702  29.921  -2.671  1.00  0.00           H
ATOM      7  C   MET A   1      13.556  31.158  -3.989  1.00 78.52           C
ATOM      8  O   MET A   1      13.301  32.365  -4.085  1.00 79.40           O
ATOM      9  CB  MET A   1      12.634  29.724  -2.148  1.00 83.38           C
ATOM     10  HB1 MET A   1      12.600  29.755  -1.059  1.00  0.00           H
ATOM     11  HB2 MET A   1      11.725  30.170  -2.552  1.00  0.00           H
ATOM     12  CG  MET A   1      12.680  28.264  -2.580  1.00 85.71           C
ATOM     13  HG1 MET A   1      12.066  28.145  -3.473  1.00  0.00           H
ATOM     14  HG2 MET A   1      12.274  27.650  -1.776  1.00  0.00           H
ATOM     15  SD  MET A   1      14.354  27.697  -2.951  1.00100.24           S
ATOM     16  CE  MET A   1      15.137  27.675  -1.338  1.00 97.10           C
ATOM     17  HE1 MET A   1      15.118  28.678  -0.913  1.00  0.00           H
ATOM     18  HE2 MET A   1      14.597  26.991  -0.683  1.00  0.00           H
ATOM     19  HE3 MET A   1      16.170  27.342  -1.441  1.00  0.00           H
ATOM     20  N   GLU A   2      13.604  30.306  -5.018  1.00 75.32           N
ATOM     21  H   GLU A   2      13.628  29.291  -4.794  1.00  0.00           H
ATOM     22  CA  GLU A   2      13.626  30.684  -6.444  1.00 72.66           C
ATOM     23  HA  GLU A   2      13.508  29.732  -6.961  1.00  0.00           H
ATOM     24  C   GLU A   2      14.903  31.384  -6.843  1.00 65.44           C
ATOM     25  O   GLU A   2      15.489  31.127  -7.892  1.00 67.46           O
ATOM     26  CB  GLU A   2      12.516  31.672  -6.810  1.00 72.71           C
ATOM     27  HB1 GLU A   2      12.094  32.068  -5.886  1.00  0.00           H
ATOM     28  HB2 GLU A   2      12.960  32.486  -7.383  1.00  0.00           H
ATOM     29  CG  GLU A   2      11.397  31.081  -7.622  1.00 77.72           C
ATOM     30  HG1 GLU A   2      10.901  31.877  -8.178  1.00  0.00           H
ATOM     31  HG2 GLU A   2      11.809  30.354  -8.321  1.00  0.00           H
ATOM     32  CD  GLU A   2      10.385  30.393  -6.731  1.00 79.89           C
ATOM     33  OE1 GLU A   2       9.171  30.574  -6.987  1.00 82.69           O
ATOM     34  OE2 GLU A   2      10.824  29.702  -5.779  1.00 72.73           O
ATOM     35  N   ASN A   3      15.194  32.399  -6.043  1.00 60.12           N
ATOM     36  H   ASN A   3      14.712  32.444  -5.123  1.00  0.00           H
ATOM     37  CA  ASN A   3      16.140  33.452  -6.363  1.00 56.21           C
ATOM     38  HA  ASN A   3      16.072  33.715  -7.419  1.00  0.00           H
ATOM     39  C   ASN A   3      17.548  32.968  -6.047  1.00 50.40           C
ATOM     40  O   ASN A   3      18.513  33.666  -6.331  1.00 50.95           O
ATOM     41  CB  ASN A   3      15.821  34.667  -5.493  1.00 57.10           C
ATOM     42  HB1 ASN A   3      15.913  34.381  -4.445  1.00  0.00           H
ATOM     43  HB2 ASN A   3      16.538  35.457  -5.717  1.00  0.00           H
ATOM     44  CG  ASN A   3      14.425  35.196  -5.732  1.00 62.07           C
ATOM     45  OD1 ASN A   3      14.220  36.001  -6.641  1.00 69.12           O
ATOM     46  ND2 ASN A   3      13.464  34.748  -4.923  1.00 57.21           N
ATOM     47 1HD2 ASN A   3      12.486  35.079  -5.045  1.00  0.00           H
ATOM     48 2HD2 ASN A   3      13.693  34.068  -4.170  1.00  0.00           H
ATOM     49  N   PHE A   4      17.643  31.764  -5.492  1.00 47.14           N
ATOM     50  H   PHE A   4      16.781  31.193  -5.378  1.00  0.00           H
ATOM     51  CA  PHE A   4      18.914  31.219  -5.038  1.00 46.77           C
ATOM     52  HA  PHE A   4      19.696  31.915  -5.340  1.00  0.00           H
ATOM     53  C   PHE A   4      19.239  29.888  -5.695  1.00 48.73           C
ATOM     54  O   PHE A   4      18.490  28.919  -5.594  1.00 49.50           O
ATOM     55  CB  PHE A   4      18.922  31.095  -3.507  1.00 42.32           C
ATOM     56  HB1 PHE A   4      18.091  30.456  -3.207  1.00  0.00           H
ATOM     57  HB2 PHE A   4      19.862  30.635  -3.201  1.00  0.00           H
ATOM     58  CG  PHE A   4      18.787  32.415  -2.806  1.00 39.42           C
ATOM     59  CD1 PHE A   4      17.545  32.878  -2.408  1.00 40.76           C
ATOM     60  HD1 PHE A   4      16.673  32.232  -2.509  1.00  0.00           H
ATOM     61  CD2 PHE A   4      19.881  33.258  -2.648  1.00 35.13           C
ATOM     62  HD2 PHE A   4      20.857  32.941  -3.015  1.00  0.00           H
ATOM     63  CE1 PHE A   4      17.397  34.138  -1.889  1.00 36.49           C
ATOM     64  HE1 PHE A   4      16.400  34.520  -1.669  1.00  0.00           H
ATOM     65  CE2 PHE A   4      19.757  34.474  -2.043  1.00 32.88           C
ATOM     66  HE2 PHE A   4      20.640  35.090  -1.873  1.00  0.00           H
ATOM     67  CZ  PHE A   4      18.507  34.929  -1.643  1.00 43.95           C
ATOM     68  HZ  PHE A   4      18.401  35.892  -1.144  1.00  0.00           H
END

here my code :
from Bio import PDB

from Bio.PDB import PICIO, PDBIO

from typing import TypedDict, Dict, Tuple

parser = PDB.PDBParser(PERMISSIVE=1, QUIET=1)

structure: PDB.Structure.Structure = parser.get_structure("prova", "example_short.pdb")

structure.atom_to_internal_coordinates()
chain: PDB.Chain.Chain = list(structure.get_chains())[0]
ic_chain: PDB.internal_coords.IC_Chain = chain.internal_coord
d:  Dict[Tuple[PDB.internal_coords.AtomKey,
                       PDB.internal_coords.AtomKey,
                       PDB.internal_coords.AtomKey,
                       PDB.internal_coords.AtomKey],
                 PDB.internal_coords.Dihedron] = ic_chain.dihedra

cnt = 1
for key in d:

    if key[0].akl[3] == 'N':
        if key[1].akl[3] == 'CA':
            if key[2].akl[3] == 'C':
                if key[3].akl[3] == 'N':

                    print ('\n',cnt,' :   ',  [x.akl[3] for x in key], d[key].angle)

                    d[key].angle += 45

                    cnt += 1

cnt = 1
for key in d:

    if key[0].akl[3] == 'N':
        if key[1].akl[3] == 'CA':
            if key[2].akl[3] == 'C':
                if key[3].akl[3] == 'N':

                    print ('\n',cnt,' :   ',  [x.akl[3] for x in key], d[key].angle)

                    cnt += 1

structure.internal_to_atom_coordinates(verbose = True)

io = PDBIO()

io.set_structure(structure)

io.save('atom_coord.pdb',  preserve_atom_numbering=True)

here an image of the original vs new quadripeptid (cartoon repr, just backbone) :
original is rainbowed one (blue to red), magenta is the new one
saved as atom_coord.pdb

I believe I changed φ (Phi) adding 45° since

Now, a φ angle is N-CA-C-N

It is worth mentioning that changing
 `d[key].angle += 45` to `d[key].angle += 6645`

works too. Apparently, need confirmation from someone more knowledgeable,
`structure.internal_to_atom_coordinates`

can cope with more than 360° angle values.

A nice and simple experiment is making a spaghetti (spaghetto?) into a helix with φ=-57.8 and ψ=-47.0

For the record, I wanted to award the bounty to this answer, but it was auto-assigned :(

I do not know if there is a function in Biopython.PDB for dihedral space transformations, but due to the lack of answers I would say no. Biopython.PDB is nice and clean, but I would use something with more molecular mechanics power... or should I say potential (:drum::laughing:).
Personally, I use PyRosetta, but it is not open source (code available to academics), so others may opt for OpenForceField.
PyRosetta requires a licence (free for academia), which is different from the Rosetta proper licence.
The code below uses a helper module I wrote, so I'll use the installer therein (which is shorthand for the official pip wheel, which provides different options etc., while the official shorthand is actually pyrosettacolabsetup)
pip3 install pyrosetta-help
install_pyrosetta -u  -p

Ligands?
First lets make sure there's nothing interesting ligand wise.
NB. RCSB has a new shiny API which may be better... Not tried it: currently using European.
import requests

def retrieve_nonpolymer_data(entry_id: str) -> dict:
    """cf. https://www.ebi.ac.uk/pdbe/api/doc/"""
    response: requests.Response = requests.get(f' https://www.ebi.ac.uk/pdbe/api/pdb/entry/molecules/{entry_id.lower()}')
    response.raise_for_status()
    return response.json()[entry_id.lower()]

info = retrieve_nonpolymer_data('1enh')
chem_resns = [c for entity in info for c in entity.get('chem_comp_ids', [])]
print(chem_resns)  # ['HOH']

If is were something that is not water, I made a list here, which is based off Peter Curran's repo for the HotSpots algorithm, which allows the elimination of worthless ligands like glycol.
I mention this because the route differs if there is a ligand.
Init and read
import pyrosetta
import pyrosetta_help as ph
from types import ModuleType

logger = ph.configure_logger()
pyrosetta.init(extra_options=ph.make_option_string(no_optH=False,
                                                ex1=None,
                                                ex2=None,
                                                #mute='all',
                                                ignore_unrecognized_res=True,
                                                load_PDB_components=False,
                                                ignore_waters=True)
                               )

# shorthand for sanity
prc: ModuleType = pyrosetta.rosetta.core

# this route strips ligands
pose: Pyrosetta.Pose = pyrosetta.toolbox.pose_from_rcsb('1enh')

Note ignore_unrecognized_res, load_PDB_components and ignore_waters.
Angles
This will not happen to an non-Rosetta user, but I get tripped out by the fact that ψ & φ data is not stored in residue, but in the pose.
for residue in pose.residues:  #: prc.conformation.Residue
    # residue.get_phi() # ---> Error!
    pass

For the residue index in pose, Rosetta was written in Fortran so is one-indexed.
import numpy as np
import numpy.typing as npt
phi:npt.NDArray[float] = np.fromiter(map(pose.phi, range(1, pose.total_residue()+1)), dtype=float)
psi:npt.NDArray[float] = np.fromiter(map(pose.psi, range(1, pose.total_residue()+1)), dtype=float)

Change φ/ψ angles
Time for some fun in dihedral space, say:
r = 2
phi = 180
pose.set_phi(r, phi)

For modelling, I wrote some code for setting spans with the archetypical angles for different SS: here. I do not recall whence the data came (maybe Wikipedia?), but it may be handy anyway.
Save
To save a PDB file:
pose.dump_pdb('.pdb')

There is also pose.dump_scored_pdb, but that is only for friends of J. Willard Gibbs.
Minimisation
It is said in the comments that energy minimisation is not required. But here is a snippet that minimises the sidechains only
scorefxn = pyrosetta.get_fa_scorefxn()  # ref2015 will do.
before: float = scorefxn(pose)
cycles = 5
movemap = pyrosetta.MoveMap()
movemap.set_bb(False)  # backbone
movemap.set_chi(True)  # sidechain
relax = pyrosetta.rosetta.protocols.relax.FastRelax(scorefxn, cycles)
relax.set_movemap(movemap)
relax.apply(pose)
after: float = scorefxn(pose)
print(after - before)  # Willard says "Remember free energy is a potential so it should be negative"

If the BB tweaks are minimal, there's also the option of scoring/minimising against the electron density etc.

Thanks for the answer. In BioPython there is functionality for setting angles in the internal_to_atom_coordinates() method. To test your solution, could you edit to a single script that does Steps 1-5 and sets all the phi/psi angles? I am not familiar with PyRosetta. Something that worries me: are you doing some behind-the-scenes optimization? I just want to make the transformation of going from dihedral to Cartesian space. A test of correctness: get back the modified dihedrals from the reconstructed atom coordinates.

Uhm... It does strike me as odd that you cannot alter the dihedrals in that case —let me check. Bio.PDB is rather basic so I would assume it might be convoluted hence the lack of Bio.PDB replies.
PyRosetta does add protons during reading and strips waters if told to do so. If told to autoload the chemical components (ligands) and the structure has an anomaly in its ligand (i.e. it is not what it ought to be) it will not be distorted on load unless minimised.

I answered about doing it with PyRosetta as I assumed that as BioPython PDB module is basic and the fact it was being asked made me think it was not possible.
Turns out I was wrong. So I am making a separate answer —as I do not actually believe doing it with Bio.PDB is the best solution.
The provided code from above with type hints is:
import numpy as np
from Bio import PDB

from typing import TypedDict, Dict, Tuple

parser = PDB.PDBParser()
structure: PDB.Structure.Structure = parser.get_structure("", "/Users/matteo/Downloads/1enh.pdb")
structure.atom_to_internal_coordinates()
chain: PDB.Chain.Chain = list(structure.get_chains())[0]
ic_chain: PDB.internal_coords.IC_Chain = chain.internal_coord
dihedrals:  Dict[Tuple[PDB.internal_coords.AtomKey,
                       PDB.internal_coords.AtomKey,
                       PDB.internal_coords.AtomKey,
                       PDB.internal_coords.AtomKey],
                 PDB.internal_coords.Dihedron] = ic_chain.dihedra

I got the types simply by using type and the attributes of note with dir. The latter helps me explore the Bio.PDB.internal_coords.Dihedron telling me of the angle attribute:
d: PDB.internal_coords.Dihedron = list(dihedrals.values())[0]
d.angle: np.float64

A change here appears it the main PDB.Structure.Structure instance, i.e. it is not a disconnected copy.
The class PDB.internal_coords.AtomKey appears to be hard to instantiate by humans and may not work as expected —help is gibberish so it is a case for inspect.getsource. So iterating across the dihedrals for the Ramachandran dihedrals is the solution. Now, a φ angle is N-CA-C-N and a ψ is C-N-CA-C (in PDB atom names the carboxylic carbon is C ). When the __str__ method of one of these PDB.internal_coords.AtomKey instances is called some values are given, which are available in the tuple attribute akl
atom_key: PDB.internal_coords.AtomKey = list(ic_chain.dihedra.keys())[0][0]
print(atom_key)  #3_R_N
print(atom_key.akl) #('3', None, 'R', 'N', None, None)

The first entry is the residue number and is a string because of insertion codes.
The fourth entry is the atom name. So iterating across the keys for the correct sequence is the trick —I know there's bounty asking for a complete solution, but it's simple, so anyone care to fill out the blanks and get the meaningless points, please do!

Thanks. Indeed I am actually looking for a complete solution that takes care also of the 'meaningless points' and that has a test of correctness.

I am at a loss here, one question so biopython doesnt give you just the backbone/Ramachandran dihedrals but you need to sort them out ?


text:
Lower case vs. upper case nucleotids in sequence vs. dots at the end
What is the difference between lower case and upper case nucleotides in a sequence? My other question is what are the dots at the end of the sequence? Some examples are shown below:

GGgG,GGgG,GGgG,GGgG,GGGG,GGGG…
GGGG,GGGG,GGGG,GGGG,GGGG,GGGG…
GGGG,gggg,GGGG,GGGG,UUUU
ggGg,ggGg,ggGg,GGgG,GGgG,GGgG…

These sequences are real experimental samples. I obtained them from QuadBase. Here is the link below:
https://quadbase.igib.res.in/
Hopefully you will provide some literature within the thread in order to understand this concept in bioinformatics.
BTW I'm a beginner in this field and hopefully this
Bioinformatics Stack Exchange site will help me find out the answers.

Please [edit] your question and tell us exactly how you got these sequences. You have linked to a page with various tools, which one gave you this output? What was your input? The meaning of upper vs lower case can vary enormously depending on context. for example, in genomic sequences, lowercase letters often mean low complexity regions, but we need to know more about these specific sequences and how you got them to be able to help you. Have you looked at the quadbase paper? http://nar.oxfordjournals.org/content/early/2016/05/16/nar.gkw425.full

GGGG,gggg,GGGG,GGGG,UUUU

Quadbase appears to be about DNA secondary structures, albeit the above example is RNA. The website states,

G-Quadruplexes or G4s are secondary DNA structures that are now established to have functionally crucial roles in various cellular processes.

Two guesses, firstly its all just a bit weird, secondly its a representation of hypoxanthine. Its more likely to be the second guess than the first.

The nomenclature is a bit weird: gggg with UUUU would infer non-canonical binding/bounding in WUSS (Washington University Secondary Structure) format for a psuedoknot  (type of secondary structure). If that wasn't the interpretation it would be kinda weird to share nomenclature with a major secondary structure format (WUSS) for secondary structure notation.

The 'dots', i.e. '.' appear to represent continuation, although in dot-bracket format (secondary structure format) this represents a nucleotide that doesn't hybridise.

Its worth reading the paper for hypoxanthine representation (normally its 'I'), because 'I' will behave like a guanine residue and would then make sense if it was being represented as 'g'. Thus, the following means no hypoxanthine,

GGGG,GGGG,GGGG,GGGG,GGGG,GGGG…

My bet is its point 2, otherwise it doesn't make sense.


text:
How to rename a file based on a directory name?
I have a lot of FASTA files that I recently scaffolded with RagTag. The output of RagTag includes a directory for each sample and the corresponding FASTA file written as "ragtag.scaffold.fasta"
I would like to go through each file and rename it like so:
Original Path with File name:
directoryName_ID/ragtag.scaffold.fasta

Changed File name:
directoryName_ID/ID_ragtag.scaffold.fasta

As you can see, the directory name includes the ID number I would like to use to rename the fasta file with (along with an underscore to space it out). Is there a way to loop through all of the directories, using that last number after the "_" in the directory name and to use that to rename my fasta file?
Here is what I have tried so far:
for filename in ./directory*/ragtag.scaffold.fasta; do
    mv filename ${./directory%*_*}_ragtag.scaffold.fasta;
done

It doesn't work at all but I know I need to use mv to rename the filename. However, when it comes to referencing the directory name the file is located in I get lost. I also have a .txt file that contains a list of IDs if that might be easier to reference. I'm still kinda new to the command line so if anyone could help that would be greatly appreciated!

What operating system are you running on?

Just execute the following immediately above the directories of interest, (post edited to include OPs "constant term")
for d in *; do
  mv "$d"/ragtag.scaffold.fasta
 "$d/${d/000mergedCluster_/}_ragtag.scaffold.fasta";
done

The key thing I assumed is that the directory name comprises a "constant term" and is suffixed by a "variable term".
Just to note the $d here is directory not file name - because that is constant in each directory.

Comments...
$d/${d/000mergedCluster_/}

The key part of the code is above.
What its doing is saying from the variable $d delete the term 000mergedCluster_ (or whatever the "constant term" is)

Thank you this works great! Only thing is that it includes the "constant term" in the new renamed file...is there a way to change it so only the "variable term" (whats after the _ in the original file name) is added to the file name?

@rimo What is your "constant term" in the directory name?

The whole name is something like "000mergedCluster_1234" so the constant term is "000mergedCluster_" and the variable term is "1234"

@rimo Can you try it now? If you copy and pasted the edited code that should work

She's perfect! Sorry that was silly of me but I completely understand what this command is doing now. Thank you so much for your help!! This works great!

The asterisk wildcard by itself will match all files and directories (except dot files/folders). There's no guarantee that $d is actually a directory. Better to use a more specific pattern, for example: for d in 000mergedCluster_*.

@Steve agreed. That should have been updated, but at the time the directory name was not known, however for d in *_* would have been better.

Assuming the output is only a list of directories and each directory includes a file called "ragtag.scaffold.fasta", you could use brace expansion to generate the SOURCE and DEST strings:
for i in directoryName_* ; do
    echo mv "${i}"/{,"${i#directoryName_}"_}ragtag.scaffold.fasta
done

But the above will only work using Bash and Zsh, and maybe some other shells. Without Bash/Zsh, just use:
fasta="ragtag.scaffold.fasta"
for i in directoryName_* ; do
    echo mv "${i}/${fasta}" "${i}/${i#directoryName_}_${fasta}"
done

Remove the echo to perform the rename.

Perl one-liner:
ls */filename | perl -lane 'print "mv $1/$3 $1/$2_$3" if /(\S+_([^\/_\s]+))\/(\S+)/'

If you want to use Perl:
#!/usr/bin/perl
use 5.030;
use warnings;

use File::Find::Rule qw( );

for my $file (File::Find::Rule->name("ragtag.scaffold.fasta")->in(".")){
  my $newfile = $file;
  $newfile =~ s{(.+_(\d+)/)(.+.fasta)}{$1$2_$3};

  if (-e $newfile) {
    warn "can't rename $file to $newfile: $newfile exists\n";
  }elsif (rename $file, $newfile) {
    # do nothing
  }else{
    warn "rename $file to $newfile failed: $!\n";
  }
}


text:
How to identify matching and unique rows across different columns in R data-frame?
I am working on a large data-frame that has 4 columns and each column has variable rows. I am trying to find unique and identical pathways between the 4 columns (each column represents a particular day of treatment with a drug). Here below is a small example.
Pathways_1Day <- c("blood","kidney","testis","No","bone","liver","intestine","lungs","ABC","pancreas","Yes","spleen")

Pathways_2Day <- c("blood","kidney","testis","eyes","bone","cells","intestine","cervix","ABC","pancreas","None")

Pathways_3Day <- c("blood","kidney","vessels","lymph","t-cells","liver","intestine","lungs","ABC","epidermis","None")

df<-data.frame(Pathways_1Day,Pathways_2Day,Pathways_3Day)

I want to get a summary of the no of pathways that are common between different timepoints (1, 2 and 3 days).
Important: The no of pathways is not the same for each day.
I have tried this:
All_pathwayNames <- df%>%group_by_all%>%count

But the desired output is not what I am trying to get.
There can be different ways to address that. It will be great if I can get matching rows infront of each other across all columns.
Regards

Are pathways on different day with the same genes? Then It is just a matter to compare the names of each day via setdiff(Pathways_3Day, Pathways_1Day) for example. If not I would recommend it to work with lists instead of data.frame for this kind of data or use GeneSetCollection from GSEABase to read it and work with this.

When I create the dataframe using the commands you have provided, I get an error, and the data frame can't be created:
Error in data.frame(Pathways_1Day, Pathways_2Day, Pathways_3Day) :
  arguments imply differing number of rows: 12, 11

What I think you were trying to do is something more like this:
library(tidyverse)
tibble() %>%
  bind_rows(tibble(day=1, pathway = Pathways_1Day)) %>%
  bind_rows(tibble(day=2, pathway = Pathways_2Day)) %>%
  bind_rows(tibble(day=3, pathway = Pathways_3Day)) -> combined.tbl

Producing a long-format data structure that looks like this:
> combined.tbl
# A tibble: 34 x 2
     day pathway
   <dbl> <chr>
 1     1 blood
 2     1 kidney
 3     1 testis
 4     1 No
 5     1 bone
 6     1 liver
 7     1 intestine
 8     1 lungs
 9     1 ABC
10     1 pancreas
# ... with 24 more rows

When laid out in this way, you can do something like grouping by pathway, and then do a pivot_wider to show the sets:
> combined.tbl %>%
     group_by(pathway) %>%
     arrange(pathway) %>%
     pivot_wider(names_from=day, values_from=1)
# A tibble: 20 x 4
# Groups:   pathway [20]
   pathway     `1`   `2`   `3`
   <chr>     <dbl> <dbl> <dbl>
 1 ABC           1     2     3
 2 No            1    NA    NA
 3 None         NA     2     3
 4 Yes           1    NA    NA
 5 blood         1     2     3
 6 bone          1     2    NA
 7 cells        NA     2    NA
 8 cervix       NA     2    NA
 9 epidermis    NA    NA     3
10 eyes         NA     2    NA
11 intestine     1     2     3
12 kidney        1     2     3
13 liver         1    NA     3
14 lungs         1    NA     3
15 lymph        NA    NA     3
16 pancreas      1     2    NA
17 spleen        1    NA    NA
18 t-cells      NA    NA     3
19 testis        1     2    NA
20 vessels      NA    NA     3

And insert a filter to select only pathways that appear in all three days (i.e. n() == 3):
> combined.tbl %>%
     group_by(pathway) %>%
     filter(n() == 3) %>%
     arrange(pathway) %>%
     pivot_wider(names_from=day, values_from=1)
# A tibble: 4 x 4
# Groups:   pathway [4]
  pathway     `1`   `2`   `3`
  <chr>     <dbl> <dbl> <dbl>
1 ABC           1     2     3
2 blood         1     2     3
3 intestine     1     2     3
4 kidney        1     2     3

Or at least two days (n() >= 2):
> combined.tbl %>%
     group_by(pathway) %>%
     filter(n() >= 2) %>%
     arrange(pathway) %>%
     pivot_wider(names_from=day, values_from=1)
# A tibble: 10 x 4
# Groups:   pathway [10]
   pathway     `1`   `2`   `3`
   <chr>     <dbl> <dbl> <dbl>
 1 ABC           1     2     3
 2 None         NA     2     3
 3 blood         1     2     3
 4 bone          1     2    NA
 5 intestine     1     2     3
 6 kidney        1     2     3
 7 liver         1    NA     3
 8 lungs         1    NA     3
 9 pancreas      1     2    NA
10 testis        1     2    NA

You're misusing data frames. Split the df into 4 vectors and use a Venn diagram generator to denerate a Venn diagram with the 4 vectors.
Adapting the first example from here
#untested
venn.diagram(
  x = list(df[,1], df[,2], df[,3], df[,4]),
  category.names = colnames(df),
  filename = 'my_venn_diagram.png',
  output=TRUE
)

Good call @RamRS


text:
somatic SNV tool for ONT samples
I'm looking for a tool to call somatic single nucleotide variations on reads from an Oxford Nanopore Technologies (ONT) run.
I have both tumor and normal samples and have already aligned them to a reference genome.
Popular tools like mutect2 have been developed for short reads, and according to this paper it is not recommended to use tools specifically designed for short reads with long reads.
I know of tools that work with ONT reads like clair3 which is used in the epi2melabs workflow for SNPs wf-human-snp,  but clair3 takes only one sample at a time, so seems for germline calling only.
Another approach would be to simply call the SNVs in both tumor and normal samples independently with clair3, and then subtract normal SNVs from tumor SNVs, but according to what I have read in a few papers and online resources like this from gatk, this is a "naive approach", and a somatic caller should use both normal and tumor samples together.

Do you know of any tool that takes both tumor and normal long read samples to detect somatic variants (similar to mutect2 but for long reads)

Have you used tools for short reads like mutect2 with long reads?

Have you attempted the "naive approach" of simply subtracting SNVs?

Thanks,
Walter.

Note that tumor/normal is only one way of doing somatic variant calling. It is unquestionably the best way, but in the clinical space, it isn't as common as you'd expect and in many cases you only have the somatic sample without the normal control. The main reason you don't want to use standard variant callers is that you need to be able to call variants with very low allelic balance (proportion of reads supporting the variant) because of the heterogeneity of somatic samples which will be a mix of tumor and normal cells, as well as the heterogeneity of the tumor cells themselves.

ClairS is now released for long-read somatic variant calling. https://github.com/HKU-BAL/Clairs

Your answer could be improved with additional supporting information. Please [edit] to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers in the help center.


text:
Retrieval protein definition from accession ID
I have a list of protein accessions:
XP_025107904.1
PVD22034.1
AJA37836.1

How can I obtain the definition of the protein, and the related organism?
Output:
XP_025107904.1  receptor-type tyrosine-protein phosphatase kappa-like isoform X2  Pomacea canaliculata

Thanks for the help,
Marco

I solved with a Python script (through biopython):
from Bio import Entrez
from Bio import SeqIO
handle = Entrez.efetch(db="protein", id= input_file, rettype="gb", retmode="text")
for seq_record in SeqIO.parse(handle, "genbank"):
print(seq_record.description + "\t" +  seq_record.id)

great! Could you please post your script as an answer and accept it so the question can be marked as answered?

@Marco could you please consider upvoting/accepting the answer. It does help with site stats and its only fair to terdon's answer

You can do this using the Entrez API. For example, using the three identifiers from your example, you can do:
curl "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=protein&id=XP_025107904.1,PVD22034.1,AJA37836.1" > file.xml

This will bring an XML file with the following contents:
$ cat file.xml
<?xml version="1.0" encoding="UTF-8" ?>
<!DOCTYPE eSummaryResult PUBLIC "-//NLM//DTD esummary v1 20041029//EN" "https://eutils.ncbi.nlm.nih.gov/eutils/dtd/20041029/esummary-v1.dtd">
<eSummaryResult>
<DocSum>
    <Id>1397721498</Id>
    <Item Name="Caption" Type="String">XP_025107904</Item>
    <Item Name="Title" Type="String">receptor-type tyrosine-protein phosphatase kappa-like isoform X2 [Pomacea canaliculata]</Item>
    <Item Name="Extra" Type="String">gi|1397721498|ref|XP_025107904.1|[1397721498]</Item>
    <Item Name="Gi" Type="Integer">1397721498</Item>
    <Item Name="CreateDate" Type="String">2018/06/07</Item>
    <Item Name="UpdateDate" Type="String">2018/06/07</Item>
    <Item Name="Flags" Type="Integer">512</Item>
    <Item Name="TaxId" Type="Integer">400727</Item>
    <Item Name="Length" Type="Integer">1271</Item>
    <Item Name="Status" Type="String">live</Item>
    <Item Name="ReplacedBy" Type="String"></Item>
    <Item Name="Comment" Type="String"><![CDATA[  ]]></Item>
    <Item Name="AccessionVersion" Type="String">XP_025107904.1</Item>
</DocSum>

<DocSum>
    <Id>1383748625</Id>
    <Item Name="Caption" Type="String">PVD22034</Item>
    <Item Name="Title" Type="String">hypothetical protein C0Q70_17837 [Pomacea canaliculata]</Item>
    <Item Name="Extra" Type="String">gi|1383748625|gb|PVD22034.1||gnl|WGS:PZQS|C0Q70_17837[1383748625]</Item>
    <Item Name="Gi" Type="Integer">1383748625</Item>
    <Item Name="CreateDate" Type="String">2018/04/26</Item>
    <Item Name="UpdateDate" Type="String">2018/04/26</Item>
    <Item Name="Flags" Type="Integer">0</Item>
    <Item Name="TaxId" Type="Integer">400727</Item>
    <Item Name="Length" Type="Integer">602</Item>
    <Item Name="Status" Type="String">live</Item>
    <Item Name="ReplacedBy" Type="String"></Item>
    <Item Name="Comment" Type="String"><![CDATA[  ]]></Item>
    <Item Name="AccessionVersion" Type="String">PVD22034.1</Item>
</DocSum>

<DocSum>
    <Id>735659541</Id>
    <Item Name="Caption" Type="String">AJA37836</Item>
    <Item Name="Title" Type="String">fucolectin-related molecule [Littorina littorea]</Item>
    <Item Name="Extra" Type="String">gi|735659541|gb|AJA37836.1|[735659541]</Item>
    <Item Name="Gi" Type="Integer">735659541</Item>
    <Item Name="CreateDate" Type="String">2014/12/25</Item>
    <Item Name="UpdateDate" Type="String">2014/12/25</Item>
    <Item Name="Flags" Type="Integer">0</Item>
    <Item Name="TaxId" Type="Integer">31216</Item>
    <Item Name="Length" Type="Integer">722</Item>
    <Item Name="Status" Type="String">live</Item>
    <Item Name="ReplacedBy" Type="String"></Item>
    <Item Name="Comment" Type="String"><![CDATA[  ]]></Item>
    <Item Name="AccessionVersion" Type="String">AJA37836.1</Item>
</DocSum>

</eSummaryResult>

You can now parse this in your favorite language/parser and get the information you need. I... am not too good with XML so my attempts with xmllint failed, but here's a naive awk approach:
$ awk -F'[][><]' '{ if(/Name="Title"/){desc=$3; spc=$4}else if(/Name="AccessionVersion"/){name=$3} if($0=="</DocSum>"){print name,desc,spc; name=""; desc=""; spc=""}}' file.xml
XP_025107904.1 receptor-type tyrosine-protein phosphatase kappa-like isoform X2  Pomacea canaliculata
PVD22034.1 hypothetical protein C0Q70_17837  Pomacea canaliculata
AJA37836.1 fucolectin-related molecule  Littorina littorea

I stress that this is a very bad idea and you should use a proper XML parser, but it might work...


text:
heatmap of specific sequence motifs in aligned fasta files
I have a collection of fasta files, each containing three aligned sequences. I am interested in understanding the distribution of a specific sequence motif PGP, RGP and KGP in all the alignments. I was wondering if there is possible way in R to create a heatmap highlighting the location of these motifs.
I have tried the code in this post to create a heatmap but it is not working as expected. One issue is that ggmsa does not import the fasta file. The fasta file has to be copied into the installation directory for it to work. I have close to 50 aligned files, so not very efficient. The other major problem is it creates a heatmap compared to a reference sequence and not the motifs. Can't figure how to change this.
Any help is much appreciated. Two fasta files containing aligned sequence are as follows,
>chain A CO1A1_HUMAN/1-1464 Collagen alpha-1(I) chain
MFSFVDLRLLLLLAATALLTHGQEEGQVEGQDEDIPPITCVQNGLRYHDRDVWKPEPCRICVCDNGKVLCDDVICDETKNCPGAEVPEGECCPVCPDGSESPTDQETTGVEGPKGDTGPRGPRGPAGPPGRDGIPGQPGLPGPPGPPGPPGPPGLGGNFAPQLSYGYDEKSTGGISVPGPMGPSGPRGLPGPPGAPGPQGFQGPPGEPGEPGASGPMGPRGPPGPPGKNGDDGEAGKPGRPGERGPPGPQGARGLPGTAGLPGMKGHRGFSGLDGAKGDAGPAGPKGEPGSPGENGAPGQMGPRGLPGERGRPGAPGPAGARGNDGATGAAGPPGPTGPAGPPGFPGAVGAKGEAGPQGPRGSEGPQGVRGEPGPPGPAGAAGPAGNPGADGQPGAKGANGAPGIAGAPGFPGARGPSGPQGPGGPPGPKGNSGEPGAPGSKGDTGAKGEPGPVGVQGPPGPAGEEGKRGARGEPGPTGLPGPPGERGGPGSRGFPGADGVAGPKGPAGERGSPGPAGPKGSPGEAGRPGEAGLPGAKGLTGSPGSPGPDGKTGPPGPAGQDGRPGPPGPPGARGQAGVMGFPGPKGAAGEPGKAGERGVPGPPGAVGPAGKDGEAGAQGPPGPAGPAGERGEQGPAGSPGFQGLPGPAGPPGEAGKPGEQGVPGDLGAPGPSGARGERGFPGERGVQGPPGPAGPRGANGAPGNDGAKGDAGAPGAPGSQGAPGLQGMPGERGAAGLPGPKGDRGDAGPKGADGSPGKDGVRGLTGPIGPPGPAGAPGDKGESGPSGPAGPTGARGAPGDRGEPGPPGPAGFAGPPGADGQPGAKGEPGDAGAKGDAGPPGPAGPAGPPGPIGNVGAPGAKGARGSAGPPGATGFPGAAGRVGPPGPSGNAGPPGPPGPAGKEGGKGPRGETGPAGRPGEVGPPGPPGPAGEKGSPGADGPAGAPGTPGPQGIAGQRGVVGLPGQRGERGFPGLPGPSGEPGKQGPSGASGERGPPGPMGPPGLAGPPGESGREGAPGAEGSPGRDGSPGAKGDRGETGPAGPPGAPGAPGAPGPVGPAGKSGDRGETGPAGPAGPVGPVGARGPAGPQGPRGDKGETGEQGDRGIKGHRGFSGLQGPPGPPGSPGEQGPSGASGPAGPRGPPGSAGAPGKDGLNGLPGPIGPPGPRGRTGDAGPVGPPGPPGPPGPPGPPSAGFDFSFLPQPPQEKAHDGGRYYRADDANVVRDRDLEVDTTLKSLSQQIENIRSPEGSRKNPARTCRDLKMCHSDWKSGEYWIDPNQGCNLDAIKVFCNMETGETCVYPTQPSVAQKNWYISKNPKDKRHVWFGESMTDGFQFEYGGQGSDPADVAIQLTFLRLMSTEASQNITYHCKNSVAYMDQQTGNLKKALLLQGSNEIEIRAEGNSRFTYSVTVDGCTSHTGAWGKTVIEYKTTKTSRLPIIDVAPLDVGAPDQEFGFDVGPVCFL
>chain B CO1A1_HUMAN/1-1464 Collagen alpha-1(I) chain
MFSFVDLRLLLLLAATALLTHGQEEGQVEGQDEDIPPITCVQNGLRYHDRDVWKPEPCRICVCDNGKVLCDDVICDETKNCPGAEVPEGECCPVCPDGSESPTDQETTGVEGPKGDTGPRGPRGPAGPPGRDGIPGQPGLPGPPGPPGPPGPPGLGGNFAPQLSYGYDEKSTGGISVPGPMGPSGPRGLPGPPGAPGPQGFQGPPGEPGEPGASGPMGPRGPPGPPGKNGDDGEAGKPGRPGERGPPGPQGARGLPGTAGLPGMKGHRGFSGLDGAKGDAGPAGPKGEPGSPGENGAPGQMGPRGLPGERGRPGAPGPAGARGNDGATGAAGPPGPTGPAGPPGFPGAVGAKGEAGPQGPRGSEGPQGVRGEPGPPGPAGAAGPAGNPGADGQPGAKGANGAPGIAGAPGFPGARGPSGPQGPGGPPGPKGNSGEPGAPGSKGDTGAKGEPGPVGVQGPPGPAGEEGKRGARGEPGPTGLPGPPGERGGPGSRGFPGADGVAGPKGPAGERGSPGPAGPKGSPGEAGRPGEAGLPGAKGLTGSPGSPGPDGKTGPPGPAGQDGRPGPPGPPGARGQAGVMGFPGPKGAAGEPGKAGERGVPGPPGAVGPAGKDGEAGAQGPPGPAGPAGERGEQGPAGSPGFQGLPGPAGPPGEAGKPGEQGVPGDLGAPGPSGARGERGFPGERGVQGPPGPAGPRGANGAPGNDGAKGDAGAPGAPGSQGAPGLQGMPGERGAAGLPGPKGDRGDAGPKGADGSPGKDGVRGLTGPIGPPGPAGAPGDKGESGPSGPAGPTGARGAPGDRGEPGPPGPAGFAGPPGADGQPGAKGEPGDAGAKGDAGPPGPAGPAGPPGPIGNVGAPGAKGARGSAGPPGATGFPGAAGRVGPPGPSGNAGPPGPPGPAGKEGGKGPRGETGPAGRPGEVGPPGPPGPAGEKGSPGADGPAGAPGTPGPQGIAGQRGVVGLPGQRGERGFPGLPGPSGEPGKQGPSGASGERGPPGPMGPPGLAGPPGESGREGAPGAEGSPGRDGSPGAKGDRGETGPAGPPGAPGAPGAPGPVGPAGKSGDRGETGPAGPAGPVGPVGARGPAGPQGPRGDKGETGEQGDRGIKGHRGFSGLQGPPGPPGSPGEQGPSGASGPAGPRGPPGSAGAPGKDGLNGLPGPIGPPGPRGRTGDAGPVGPPGPPGPPGPPGPPSAGFDFSFLPQPPQEKAHDGGRYYRADDANVVRDRDLEVDTTLKSLSQQIENIRSPEGSRKNPARTCRDLKMCHSDWKSGEYWIDPNQGCNLDAIKVFCNMETGETCVYPTQPSVAQKNWYISKNPKDKRHVWFGESMTDGFQFEYGGQGSDPADVAIQLTFLRLMSTEASQNITYHCKNSVAYMDQQTGNLKKALLLQGSNEIEIRAEGNSRFTYSVTVDGCTSHTGAWGKTVIEYKTTKTSRLPIIDVAPLDVGAPDQEFGFDVGPVCFL
>chain C CO1A2_HUMAN/1-1366 Collagen alpha-2(I) chain
----------------------------------------------------------------------------------------MLSFVDTRTLLLLAVTLCLATCQSLQEETVRKGPAGDRGPRGERGPPGPPGRDGEDGPTGPPGPPGPPGPPGLGGNFAAQYDGKGVGLGPGPMGLMGPRGPPGAAGAPGPQGFQGPAGEPGEPGQTGPAGARGPAGPPGKAGEDGHPGKPGRPGERGVVGPQGARGFPGTPGLPGFKGIRGHNGLDGLKGQPGAPGVKGEPGAPGENGTPGQTGARGLPGERGRVGAPGPAGARGSDGSVGPVGPAGPIGSAGPPGFPGAPGPKGEIGAVGNAGPAGPAGPRGEVGLPGLSGPVGPPGNPGANGLTGAKGAAGLPGVAGAPGLPGPRGIPGPVGAAGATGARGLVGEPGPAGSKGESGNKGEPGSAGPQGPPGPSGEEGKRGPNGEAGSAGPPGPPGLRGSPGSRGLPGADGRAGVMGPPGSRGASGPAGVRGPNGDAGRPGEPGLMGPRGLPGSPGNIGPAGKEGPVGLPGIDGRPGPIGPAGARGEPGNIGFPGPKGPTGDPGKNGDKGHAGLAGARGAPGPDGNNGAQGPPGPQGVQGGKGEQGPPGPPGFQGLPGPSGPAGEVGKPGERGLHGEFGLPGPAGPRGERGPPGESGAAGPTGPIGSRGPSGPPGPDGNKGEPGVVGAVGTAGPSGPSGLPGERGAAGIPGGKGEKGEPGLRGEIGNPGRDGARGAPGAVGAPGPAGATGDRGEAGAAGPAGPAGPRGSPGERGEVGPAGPNGFAGPAGAAGQPGAKGERGAKGPKGENGVVGPTGPVGAAGPAGPNGPPGPAGSRGDGGPPGMTGFPGAAGRTGPPGPSGISGPPGPPGPAGKEGLRGPRGDQGPVGRTGEVGAVGPPGFAGEKGPSGEAGTAGPPGTPGPQGLLGAPGILGLPGSRGERGLPGVAGAVGEPGPLGIAGPPGARGPPGAVGSPGVNGAPGEAGRDGNPGNDGPPGRDGQPGHKGERGYPGNIGPVGAAGAPGPHGPVGPAGKHGNRGETGPSGPVGPAGAVGPRGPSGPQGIRGDKGEPGEKGPRGLPGLKGHNGLQGLPGIAGHHGDQGAPGSVGPAGPRGPAGPSGPAGKDGRTGHPGTVGPAGIRGPQGHQGPAGPPGPPGPPGPPGVSGGGYDFGYDGDFYRADQPRSAPSLRPKDYEVDATLKSLNNQIETLLTPEGSRKNPARTCRDLRLSHPEWSSGYYWIDPNQGCTMDAIKVYCDFSTGETCIRAQPENIPAKNWYRSSKDKKHVWLGETINAGSQFEYNVEGVTSKEMATQLAFMRLLANYASQNITYHCKNSIAYMDEETGNLKKAVILQGSNDVELVAEGNSRFTYTVLVDGCSKKTNEWGKTIIEYKTNKPSRLPFLDIAPLDIGGADQEFFVDIGPVCFK

and
>chain A CO4A1_HUMAN/1-1669 Collagen alpha-1(IV) chain
---------MGPRLSVWLLLLPAALLLHEEHSRAAA--KGGCAGSGC-GKCDCHGVKGQKGERGLPGLQGVIGFPGMQGPEGPQGPPGQKGDTGEPGLPGTKGTRGPPGASGYPGNPGLPGIPGQDGPPGPPGIPGCNGTKGERGPLGPPGLPGFAGNPGPPGLPGMKGDPGEILGHVPGMLLKGERGFPGIPGTPGPPGLPGLQGPVGPPGFTGPPGPPGPPGPPGEKGQMGLSFQGPKGDKGDQGVSGPPGVPGQA-------QVQEKGDFATKGEKGQKGEPGFQGMPGVGEKGEPGKPGPRGKPGKDGDKGEKGSPGFPGEPGYPGLIGRQGPQGEKGEAGPPGPPGIVIGTGPLGEKGERGYPGTPGPRGEPGPKGFPGLPGQPGPPGLPVPGQAGAPGFPGERGEKGDRGFPGTS-LP-GPSGRDGLPGPPGSPGPPGQPGYTNGIVECQPGPPGDQGPPGIPGQPGFIGEIGEKGQKGESCLICDIDGYRGPPGPQGPPGEIGFPGQPGAKGDRGLPGRDGVAGVPGPQGTPGLIGQPGAKGEPGEFYFDLRLKGDKGDPGFPGQPGMPGRAGSPGRDGHPGLPGPKGSPGSVGLKGERGPPGGVGFPGSRGDTGPPGPPGY---GPAGPIGDKGQAGFPGGPGSPGLPGPKGEPGKIVP--------------------LPGPPGAEGLPGSPGFPGPQGDRGFPGTPGRPGLPGEKGAVGQPGI-GFPGPPGPKGVDGLPGDMGPPGTPGRPGFNGLPGNPGVQGQKGEP---GVGLPGLKGLPGLPGIPGTPGEKGSIGVPGVPGEHGAIGPPGLQGIRGEPGPPGLPGSVGSPGVPGI-GPPGARGPPGGQGPPGLSGPPGIKGEKGFPGFPGLD-MPGPKGDKGAQGLPGITGQSGLPGLPGQQGAPGIPGFPGSKGEMGVMGTPGQPGSPGPVGAPGLPGEKGDHGFPGSSGPRGDPGLKGDKGDVGLPGKPGSMDKVDMGSMKGQKGDQGEKGQIGPIGEKGSRGDPGTPGVPGKDGQAGQPGQP-GPKGDPGISGTPGAPGLPGPKGSVGGMGLPGTPGEKGVPGIPGPQGSPGLPGDKGAKGEKGQAGPPGIGIPGLRGEKGDQGIAGFPGSPGEKGEKGSIGIPGMPGSPGLKGSPGSVGYPGSPGLPGEKGDKGLPGLDGIPGVKGEAGLPGTPGPTGPAGQKGEPGSDGIPGSAGEKGEPGLPGRGFPGFPGAKGDKGSKGEVGFP-GLAGSPGIPGSKGEQGFMGPPGPQGQPGLPGSPGHA-TEGPKGDRGPQGQPGLPGLPGPMGPPGLPGIDGVKGDKGNPGWPGAPGVPGPKGDPGFQGMPGIGGSPGITGSKGDMGPPGVPGFQGPKGLPGLQGIKGDQGDQGVPGAKGLPGPPGPPGPYDIIKGEPGLPGPEGPPGLKGLQGLPGPKGQQGVTGLVGIPGPPGIPGFDGAPGQKGEMGPAGPTGPRGFPGPPGPDGLPGSMGPPGTPSVDHGFLVTRHSQTIDDPQCPSGTKILYHGYSLLYVQGNERAHGQDLGTAGSCLRKFSTMPFLFCNINNVCNFASRNDYSYWLSTPEPMPMSMAPITGENIRPFISRCAVCEAPAMVMAVHSQTIQIPPCPSGWSSLWIGYSFVMHTSAGAEGSGQALASPGSCLEEFRSAPFIECHG-RGTCNYYANAYSFWLATIERSEMFKKPTPSTLKAGELRTHVSRCQVCMRRT
>chain B CO4A1_HUMAN/1-1669 Collagen alpha-1(IV) chain
---------MGPRLSVWLLLLPAALLLHEEHSRAAA--KGGCAGSGC-GKCDCHGVKGQKGERGLPGLQGVIGFPGMQGPEGPQGPPGQKGDTGEPGLPGTKGTRGPPGASGYPGNPGLPGIPGQDGPPGPPGIPGCNGTKGERGPLGPPGLPGFAGNPGPPGLPGMKGDPGEILGHVPGMLLKGERGFPGIPGTPGPPGLPGLQGPVGPPGFTGPPGPPGPPGPPGEKGQMGLSFQGPKGDKGDQGVSGPPGVPGQA-------QVQEKGDFATKGEKGQKGEPGFQGMPGVGEKGEPGKPGPRGKPGKDGDKGEKGSPGFPGEPGYPGLIGRQGPQGEKGEAGPPGPPGIVIGTGPLGEKGERGYPGTPGPRGEPGPKGFPGLPGQPGPPGLPVPGQAGAPGFPGERGEKGDRGFPGTS-LP-GPSGRDGLPGPPGSPGPPGQPGYTNGIVECQPGPPGDQGPPGIPGQPGFIGEIGEKGQKGESCLICDIDGYRGPPGPQGPPGEIGFPGQPGAKGDRGLPGRDGVAGVPGPQGTPGLIGQPGAKGEPGEFYFDLRLKGDKGDPGFPGQPGMPGRAGSPGRDGHPGLPGPKGSPGSVGLKGERGPPGGVGFPGSRGDTGPPGPPGY---GPAGPIGDKGQAGFPGGPGSPGLPGPKGEPGKIVP--------------------LPGPPGAEGLPGSPGFPGPQGDRGFPGTPGRPGLPGEKGAVGQPGI-GFPGPPGPKGVDGLPGDMGPPGTPGRPGFNGLPGNPGVQGQKGEP---GVGLPGLKGLPGLPGIPGTPGEKGSIGVPGVPGEHGAIGPPGLQGIRGEPGPPGLPGSVGSPGVPGI-GPPGARGPPGGQGPPGLSGPPGIKGEKGFPGFPGLD-MPGPKGDKGAQGLPGITGQSGLPGLPGQQGAPGIPGFPGSKGEMGVMGTPGQPGSPGPVGAPGLPGEKGDHGFPGSSGPRGDPGLKGDKGDVGLPGKPGSMDKVDMGSMKGQKGDQGEKGQIGPIGEKGSRGDPGTPGVPGKDGQAGQPGQP-GPKGDPGISGTPGAPGLPGPKGSVGGMGLPGTPGEKGVPGIPGPQGSPGLPGDKGAKGEKGQAGPPGIGIPGLRGEKGDQGIAGFPGSPGEKGEKGSIGIPGMPGSPGLKGSPGSVGYPGSPGLPGEKGDKGLPGLDGIPGVKGEAGLPGTPGPTGPAGQKGEPGSDGIPGSAGEKGEPGLPGRGFPGFPGAKGDKGSKGEVGFP-GLAGSPGIPGSKGEQGFMGPPGPQGQPGLPGSPGHA-TEGPKGDRGPQGQPGLPGLPGPMGPPGLPGIDGVKGDKGNPGWPGAPGVPGPKGDPGFQGMPGIGGSPGITGSKGDMGPPGVPGFQGPKGLPGLQGIKGDQGDQGVPGAKGLPGPPGPPGPYDIIKGEPGLPGPEGPPGLKGLQGLPGPKGQQGVTGLVGIPGPPGIPGFDGAPGQKGEMGPAGPTGPRGFPGPPGPDGLPGSMGPPGTPSVDHGFLVTRHSQTIDDPQCPSGTKILYHGYSLLYVQGNERAHGQDLGTAGSCLRKFSTMPFLFCNINNVCNFASRNDYSYWLSTPEPMPMSMAPITGENIRPFISRCAVCEAPAMVMAVHSQTIQIPPCPSGWSSLWIGYSFVMHTSAGAEGSGQALASPGSCLEEFRSAPFIECHG-RGTCNYYANAYSFWLATIERSEMFKKPTPSTLKAGELRTHVSRCQVCMRRT
>chain C CO4A2_HUMAN/1-1712 Collagen alpha-2(IV) chain
MGRDQRAVAGPALRRWLLLGTVTVGFLAQSVLAGVKKFDVPCGGRDCSGGCQCYPEKGGRGQPGPVGPQGYNGPPGLQGFPGLQGRKGDKGERGAPGVTGPKGDVGARGVSGFPGADGIPGHPGQGGPRGRPGYDGCNGTQGDSGPQGPPGSEGFTGPPGPQGPKGQKGEP-YALPKEERDRYRGEPGEPGLVGFQGPPGRPGHVGQMGPVGAPGRPGPPGPPGPKGQQGNRGLGFYGVKGEKGDVGQPGPNGIPSDTLHPIIAPTGVTFHPDQYKGEKGSEGEPGIRGISLKGEEGIMGFPGLRGYPGLSGEKGSPGQKGSRGLDGYQGPDGPRGPKGEAGDPGPPGLP--AYSPHPSLAKGARGDPGFPGAQGEPGSQGEPGDPGLPGPPGLSIGDGDQRRGLPGEMGPKGFIGDPGIPALYGGPPGPDGKRGPPGPPGLPGPPGPDGFL-FGLKGAKGRAGFPGLPGSPGARGPKGWKGDAGECRCTEGDEAIKGLPGLPGPKGFAGINGEPGRKGDRGDPGQHGLPGFPGLKGVPGNIGAPGPKGAKGDS-RTITTKGERGQPGVPGVPGMKGDDGSPGRDGLDGFPGLPGPPGD-GIKGPPGDPGYPGIPGTKGTPGEMGPPGLGLPGLKGQRGFPGDAGLPGPPGFLGPPGPAGTPGQIDCDTDVKRAVGGDRQEAIQPGCIGGPKGLPGLPGPPGPTGAKGLRGIPGFAGADGGPGPRGLPGDAGREGFPGPPGFIGPRGSKGAVGLPGPDGSPGPIGLPGPDGPPGERGLPGEVLGAQPGPRGDAGVPGQPGLKGLPGDRGPPGFRGSQGMPGMPGLKGQPGLPGPSGQPGLYGPPGLHGFPGAPGQEGPLGLPGIPGREGLPGDRGDPGDTGAPGPVGMKGLSGDRGDAGFTGEQGHPGSPGFKGIDGMPGTPGLKGDRGSPGMDGFQGMPGLKGRPGFPGSKGEAGFFGIPGLKGLAGEPGFKGSRGDPGPPGPP-PVILPGMKDIKGEKGDEGPMGLKGYLGAKGIQGMPGIPGLSGIPGLPGRPGHIKGVKGDIGVPGIPGLPGFPGVAGPPGITGFPGFIGSRGDKGAPGRAGLYGEIGATGDFGDIGDT-INLPGRPGLKGERGTTGIPGLKGFFGEKGTEGDIGFPGITGVTGVQGPPGLKGQTGFPGLTGPPGSQGELGRIGLPGGKGDDGWPGAPGLPGFPGLRGIRGLHGLPGTKGFPGSPGSDIHGDPGFPGPPGERGDPGEANTLPGPVGVPGQKGDQGAPGERGPPGSPGLQGFPGITPPSNISGAPGDKGAPGIFGLKGYRGPPGPPGSAALPGSKGDTGNPGAPGTPGTKGWAGDSGPQGRPGVFGLPGEKGPRGEQGFMGNTGPTGAVGDRGPKGPKGDPGFPGAPGTVGAPGIAGIPQKIAVQPGTVGPQGRRGPPGAPGEMGPQGPPGEPGFRGAPGKAGPQGRGGVSAVPGFRGDEGPIGHQGPIGQEGAPGRPGSPGLPGMPGR-SVSIGYLLVKHSQTDQEPMCPVGMNKLWSGYSLLYFEGQEKAHNQDLGLAGSCLARFSTMPFLYCNPGDVCYYASRNDKSYWLSTTAPLP--MMPVAEDEIKPYISRCSVCEAPAIAIAVHSQDVSIPHCPAGWRSLWIGYSFLMHTAAGDEGGGQSLVSPGSCLEDFRATPFIECNGGRGTCHYYANKYSFWLTTIPEQSFQGSPSADTLKAGLIRTHISRCQVCMKNL

When I want an MSA visualization in R that I can't find pre-made, my approach is to put the text into a data frame of position/character combos and use generic plotting tools like pheatmap or ggplot2's geom_tile.  Looks like that'd get tricky here because of the gaps and multi-AA motifs though.  Maybe it'd work with a mapping of gapped/ungapped locations of all AA trios instead of individual characters?  And for simple FASTA reading maybe check out this minimalistic package from a colleague of mine.  It just reads the text into data frames.

Here's a solution based on Jesse's comment. The steps are:

read in MSA with Biostrings::readAAStringSet()
right-fill sequences so they are all the same length (assumes left alignment)
search for motif substrings in each amino acid sequence with gregexpr(). If a motif in the motif set is found, record position of center of motif.
bind the results from each sequence to a single table

find_motifs <- function(msa, motifs) {

  require(tibble)
  require(Biostrings)

  # read in fasta as Biostrings object
  seqs <- Biostrings::readAAStringSet(filepath = msa, format = "fasta")
  maxseqlen <- max(width(seqs))

  # initiate results tibble
  results <- tibble::tibble(position = seq(1:maxseqlen))

  # readAAMultipleAlignment requires equal length of sequences
  # instead, use readAAStringSet and fill "-" on right
  # assumes left-aligned MSA
  for (i in 1:length(seqs)) {

    seqlen <- nchar(seqs[i])

    if (seqlen < maxseqlen) {
      fill <- maxseqlen - seqlen
      seqs[[i]] <- paste0(as.character(seqs[[i]]),
                          paste(rep("-", fill), collapse = ""))
    }

    # initiate vector of zeros to store hit positions
    positions <- rep(0, maxseqlen)

    lapply(motifs, function(motif) {

      # factor to add to gregexpr output to center hit
      # for motifs with an even number of characters, use left center
      shift <- ceiling(nchar(motif)/2)

      hits <- unlist(gregexpr(motif, seqs[[i]])) + shift

      lapply(hits, function(hit) {
        positions[hit] <<- 1
      }) # end apply over hits

    }) # end apply over motifs

    results <- dplyr::bind_cols(results, positions)

  } # end loop over sequences

  names(results) <- c("position", paste0("seq", seq(1:length(seqs))))
  results

} # end function

Running this function on the provided fastas, we get a tables where positions in sequences matching motifs are recorded as 1.
setwd("/path/to/files")
ngp <- c("PGP","RGP","KGP")
msa1 <- find_motifs(msa = "msa1.fasta",
                    motifs = ngp)
msa2 <- find_motifs(msa = "msa2.fasta",
                    motifs = ngp)

dplyr::slice(msa1, 120:130)

# A tibble: 11 × 4
   position  seq1  seq2  seq3
      <int> <dbl> <dbl> <dbl>
 1      120     0     0     0
 2      121     0     0     0
 3      122     1     1     1
 4      123     0     0     0
 5      124     0     0     0
 6      125     1     1     0
 7      126     0     0     0
 8      127     0     0     0
 9      128     0     0     1
10      129     0     0     0
11      130     0     0     0

And here's a function to make the plots.
heatmap_separate <- function(data) {
  require(ggplot2)
  require(tidyr)
  numseqs <- ncol(data) - 1
  numrows <- nrow(data)
  toplot <- data |> tidyr::pivot_longer(cols = -position,
                                        names_to = "sequence",
                                        values_to = "hits")
    ggplot(toplot) +
    geom_tile(mapping = aes(x = position,
                            y = sequence,
                            fill = as.factor(hits),
                            color = as.factor(hits))) +
    scale_x_continuous(limits = c(1,numrows),
                       n.breaks = 10,
                       expand = c(0,0)) +
    scale_y_discrete(expand = c(0,0),
                     limits = rev) +
    scale_fill_manual(values = c("#f2f0f7","#54278f")) +
    scale_color_manual(values = c("#f2f0f7","#54278f")) +
    geom_hline(yintercept = 0.5 + 0:numseqs,
               color = "black",
               linewidth = 1.5) +
    theme_minimal() +
    theme(legend.position = "none",
          axis.title.y = element_blank(),
          panel.grid = element_blank(),
          axis.ticks.x = element_line(color = "black"),
          panel.background = element_blank())
}

heatmap_separate(msa1)
heatmap_separate(msa2)

Plot for fasta #1

Plot for fasta #2

The lines representing hits are all plotted as the same color (#54278f), but they appear to be different colors as an artifact of plotting lots of narrow tiles. This can be rectified by viewing the plots as svg files, though I can't display that file type here.


text:
Command line based interaction with the UCSC genome browser, automated data upload and track generation
Is there a way to submit data to ucsc genome browser for result/track display, on the fly, as a local bedgraph file with example contents like below:
browser position chr19:49302001-49304701
browser hide all
browser pack refGene encodeRegions
browser full altGraph
#   300 base wide bar graph, autoScale is on by default == graphing
#   limits will dynamically change to always show full range of data
#   in viewing window, priority = 20 positions this as the second graph
#   Note, zero-relative, half-open coordinate system in use for bedGraph format
track type=bedGraph name="BedGraph Format" description="BedGraph format" visibility=full color=200,100,0 altColor=0,100,200 priority=20
chr19 49302000 49302300 -1.0
chr19 49302300 49302600 -0.75
chr19 49302600 49302900 -0.50
chr19 49302900 49303200 -0.25
chr19 49303200 49303500 0.0
chr19 49303500 49303800 0.25
chr19 49303800 49304100 0.50
chr19 49304100 49304400 0.75
chr19 49304400 49304700 1.00

Normally, you can input the data in the data box like below:

and submit - and display results after another [GO] button that's pushed manually on the following web page to produce the output below:

I am hoping there is a command line option for this where it takes the data input above as a local file name and generates the output that normally one can generate with a few manual clicks (including pasting of the data text into the window or upload of the data file) on the ucsc web page.
Thanks in advance


text:
BCL files half the size of fastq files (10x chromium single cell)
I was just curious about this. The BCL files (~100GB) we have for a single cell experiment (10x chromium), are less than half the size of the produced fastq files (240GB, gzipped). This is excluding the uncharacterized reads fastqs.
Is this usual? All other types of sequencing data demultiplexes do not have this increase in file size. Is this because there are 4 fastqs per sample (indexes and reads).
Thanks in advance,

I do not have direct experience with BCL, so this is just a guess. FASTQ is text-based, while BCL is in binary format and probably optimized to be efficiently stored. The conversion BCL2FASTQ might require some filtering. The level of compression can also be a factor.

yes, this would be my first guess too.


text:
What can be the bias of aligning paired-reads in a single-end mode?
I don't know if I'm in the right place but I have a technical problem to fix. I would have to align paired-end reads from Illumina sequencing to compare a normal genome with a tumor one.
When I align with bwa-mem (default parameters) with paired-end mode I get a bigWig file where the coverage of the two samples is quite similar. When I align with bwa-mem (default parameters) with single-end mode I get coverage holes (definitely does not align any read) in some specific regions. But this only happens in tumor genomes.
Is it possible that something is really going on in that region? Or is it just a bias? And if it's a bias, why in the normal condition, in that specific region, do they align?
I think I understand that by treating paired-end reads as single-end I lose the directionality of mate2 if I align mate 1. But I still can't understand why I have this alignment defect only in the tumor genome and not in the normal genome (I observe this phenomenon in all the samples analyzed, not just in one).
Could someone give me some explanation or recommend any checks to do?
spoiler - these holes fall exactly in highly repetitive regions

In principle one should simply not align paired-end reads as if they were single end reads, not only because you are not taking advantage of the forward/reverse information contained in the reads, but also you are not taking into account the expected distance between the two pairs, both of which result in better mapping rates/accuracy.
Having said that, and not recalling from the top of my head right now the relevant BWA-mem algorithmic implementation details, it is generally difficult to map (even) paired-reads to repetitive regions. Still there will be some rate of success. In tumor samples however, we could reasonably imagine that chromosomal rearrangements, INDELs, CNVs, etc. will make it even harder to confidently assign these reads to a location in the reference, since these extra variation will increase e.g. multiple mapping locations and mapping qualities for some fraction of the reads, thus resulting in regions with low or no coverage compared to non-tumor samples.
Repetitive regions might exacerbate the problem and get affected even further compared to other regions.

I tried to give me this explanation: with single-end alignment of paired-end reads, losing the directionality of the reads and also the expected fragment length, the R2 no longer maps if it maps the R1, but can map as a single reads anywhere, or don't map at all.
I thought that therefore I don't see alignment (in tumors) in repetitive regions because of the high mutation rate I have in the tumor sample. While in normal I don't have this situation.
However, aligning with paired-end reads, I see coverage in these regions again despite the high mutation rate because R1 and R2 help each other map in the best position, and so even if I have mutations I still see the alignment.
Could it be appropriate to make a variant call despite being on repetitive regions?

As it’s currently written, your answer is unclear. Please [edit] to add additional details that will help others understand how this addresses the question asked. You can find more information on how to write good answers in the help center.

Check on the topic of 'mapability' (e.g. 1000 genomes or UCSC tracks, snpable, ..). Often it is practice to mask out regions of poor mapability in order to decrease the rate of false positives. Repetitive regions fall often within poorly mappable regions, your read length becomes somewhat relevant here


text:
BAF and LRR calculation and usage
I'm  working with a vcf file, in particular I'm trying to identify a gene duplication in a sequenced genome.
I studied bcftools cnv command, but for this tool is needed a vcf file containing a BAF and LRR values, that there are not present in my vcf file.
And so the first question is, can I calculate them starting from read depth and allele depth? Once that I have calculate BAF and LRR, can I modify the original vcf file inserting the calculated values? To have a vcf file containing BAF and LRR values I have to ask to the geneticist to use particular pipeline in NGS instrument?

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


text:
Creating a PDB from a figure in paper
I want to create a PDB file of a drug that is given in a paper.
Can somebody share a tutorial or similar resources for a beginner-level user?

Thank you.

Could you add some context? What do you want the PDB file for? Do you have the structure of this molecule in another file format? (I didn't read the paper).

yep have same question. Wouldn' t a cif file better serve the use if needed for fitting a ligand into electron density ? Or you just nred to depict molecule in 3D ?

You can use Image2Smiles tool. It produces smi format, then you can convert smi to pdb using openbabel.
If you just have a few figures, I suggest using PubChem Sketcher.

or you could taket it from :
https://www.rcsb.org/structure/7WMV pdb file
HETATM 4865  C13 1YI A 701     106.585 114.218 112.698  1.00 74.13           C
HETATM 4866  C15 1YI A 701     106.050 112.633 114.547  1.00102.13           C
HETATM 4867  C17 1YI A 701     104.024 111.140 110.766  1.00 89.07           C
HETATM 4868  C20 1YI A 701     103.626 108.779 110.527  1.00 83.00           C
HETATM 4869  C21 1YI A 701     104.626 108.507 111.437  1.00 84.86           C
HETATM 4870  C22 1YI A 701     102.234 110.466 109.205  1.00 85.32           C
HETATM 4871  C24 1YI A 701     100.911 110.073 109.846  1.00 66.83           C
HETATM 4872  C26 1YI A 701      92.545 108.407 114.139  1.00 95.68           C
HETATM 4873  C28 1YI A 701     100.141 109.053 109.319  1.00 63.43           C
HETATM 4874  C11 1YI A 701     105.694 113.343 111.806  1.00 75.91           C
HETATM 4875  C12 1YI A 701     105.840 111.929 112.354  1.00 91.76           C
HETATM 4876  C14 1YI A 701     105.990 114.087 114.099  1.00 85.43           C
HETATM 4877  C16 1YI A 701     105.029 110.853 111.663  1.00 99.03           C
HETATM 4878  C18 1YI A 701     105.319 109.548 112.006  1.00100.47           C
HETATM 4879  C19 1YI A 701     103.333 110.096 110.188  1.00 81.96           C
HETATM 4880  C23 1YI A 701     104.801 110.649 116.200  1.00 88.38           C
HETATM 4881  C25 1YI A 701     102.892 107.590 109.921  1.00 99.79           C
HETATM 4882  C27 1YI A 701     100.517 110.745 110.979  1.00 88.80           C
HETATM 4883  C29 1YI A 701      98.556 109.383 111.080  1.00 77.66           C
HETATM 4884  C30 1YI A 701      97.245 109.006 111.756  1.00 83.68           C
HETATM 4885  C31 1YI A 701      97.108 107.503 111.935  1.00 85.27           C
HETATM 4886  C32 1YI A 701      99.328 110.410 111.602  1.00 90.75           C
HETATM 4887  C33 1YI A 701      98.954 108.713 109.937  1.00 76.18           C
HETATM 4888  C34 1YI A 701      96.107 107.127 113.039  1.00 98.03           C
HETATM 4889  C35 1YI A 701      94.780 107.913 113.042  1.00 92.89           C
HETATM 4890  C36 1YI A 701      91.696 107.883 115.298  1.00113.05           C
HETATM 4891  C37 1YI A 701      92.889 109.868 114.425  1.00 87.90           C
HETATM 4892  C38 1YI A 701      91.711 108.277 112.861  1.00 77.85           C
HETATM 4893  N8  1YI A 701      93.804 107.673 114.094  1.00 93.06           N
HETATM 4894  N9  1YI A 701      90.968 106.633 115.149  1.00106.00           N
HETATM 4895  O2  1YI A 701     105.306 111.885 113.648  1.00 99.43           O
HETATM 4896  O3  1YI A 701     106.050 113.444 110.457  1.00 80.66           O
HETATM 4897  O4  1YI A 701     106.560 115.552 112.284  1.00 86.11           O
HETATM 4898  O5  1YI A 701     106.634 114.878 115.050  1.00 99.48           O
HETATM 4899  O6  1YI A 701      94.512 108.703 112.208  1.00 84.41           O
HETATM 4900  O7  1YI A 701      91.625 108.487 116.316  1.00101.67           O
HETATM 4901  S1  1YI A 701     105.267 112.404 116.160  1.00121.62           S

picture :

obtained with : https://litemol.org/Viewer/
OOOPSSSSSSSSSS
realized that ligand in pdb its not the full ligand of your picture,
need to get it from Ligand page on RCSB_PDB : https://www.rcsb.org/ligand/1YI
download links highlighted in pics :

download /1YI_ideal.sdf  convert via openbabel :
http://www.cheminfo.org/Chemistry/Cheminformatics/FormatConverter/index.html
result
WARNING: damped steplength 0.0184262 to 0.00725371
WARNING: damped steplength 0.00508674 to 0.00100138
COMPND    1YI
AUTHOR    GENERATED BY OPEN BABEL 2.3.2
HETATM    1  C   LIG     1       3.537  -2.101   2.235  1.00  0.00           C
HETATM    2  C   LIG     1       2.812  -4.346   1.266  1.00  0.00           C
HETATM    3  C   LIG     1       2.262  -2.902  -2.712  1.00  0.00           C
HETATM    4  C   LIG     1       0.447  -1.602  -3.663  1.00  0.00           C
HETATM    5  C   LIG     1       0.405  -0.869  -2.474  1.00  0.00           C
HETATM    6  C   LIG     1       1.522  -3.533  -5.065  1.00  0.00           C
HETATM    7  C   LIG     1       0.305  -3.796  -5.950  1.00  0.00           C
HETATM    8  C   LIG     1      -5.920  -1.833 -12.368  1.00  0.00           C
HETATM    9  C   LIG     1      -0.581  -4.848  -5.684  1.00  0.00           C
HETATM   10  C   LIG     1       2.908  -1.583   0.938  1.00  0.00           C
HETATM   11  C   LIG     1       3.139  -2.484  -0.313  1.00  0.00           C
HETATM   12  C   LIG     1       3.714  -3.615   2.273  1.00  0.00           C
HETATM   13  C   LIG     1       2.202  -2.172  -1.496  1.00  0.00           C
HETATM   14  C   LIG     1       1.247  -1.156  -1.414  1.00  0.00           C
HETATM   15  C   LIG     1       1.366  -2.666  -3.802  1.00  0.00           C
HETATM   16  C   LIG     1       1.091  -5.111   3.310  1.00  0.00           C
HETATM   17  C   LIG     1      -0.531  -1.162  -4.707  1.00  0.00           C
HETATM   18  C   LIG     1      -0.009  -2.939  -7.020  1.00  0.00           C
HETATM   19  C   LIG     1      -2.131  -4.074  -7.418  1.00  0.00           C
HETATM   20  C   LIG     1      -3.479  -4.118  -8.115  1.00  0.00           C
HETATM   21  C   LIG     1      -3.492  -4.837  -9.473  1.00  0.00           C
HETATM   22  C   LIG     1      -1.224  -3.051  -7.715  1.00  0.00           C
HETATM   23  C   LIG     1      -1.778  -4.991  -6.415  1.00  0.00           C
HETATM   24  C   LIG     1      -4.883  -4.859 -10.129  1.00  0.00           C
HETATM   25  C   LIG     1      -5.503  -3.491 -10.484  1.00  0.00           C
HETATM   26  C   LIG     1      -7.501  -1.753 -12.406  1.00  0.00           C
HETATM   27  C   LIG     1      -5.427  -1.719 -13.834  1.00  0.00           C
HETATM   28  C   LIG     1      -5.340  -0.600 -11.644  1.00  0.00           C
HETATM   29  N   LIG     1      -5.408  -3.127 -11.821  1.00  0.00           N
HETATM   30  N   LIG     1      -8.129  -2.821 -13.038  1.00  0.00           N
HETATM   31  O   LIG     1       3.072  -3.897  -0.073  1.00  0.00           O
HETATM   32  O   LIG     1       3.421  -0.244   0.755  1.00  0.00           O
HETATM   33  O   LIG     1       2.698  -1.683   3.334  1.00  0.00           O
HETATM   34  O   LIG     1       5.094  -3.922   1.972  1.00  0.00           O
HETATM   35  O   LIG     1      -6.061  -2.821  -9.650  1.00  0.00           O
HETATM   36  O   LIG     1      -8.120  -0.753 -12.076  1.00  0.00           O
HETATM   37  S   LIG     1       1.026  -4.388   1.646  1.00  0.00           S
HETATM   38  C   LIG     1      -9.774  -4.156 -16.922  1.00  0.00           C
HETATM   39  C   LIG     1      -9.839  -1.879 -16.189  1.00  0.00           C
HETATM   40  N   LIG     1      -9.383  -3.231 -15.850  1.00  0.00           N
HETATM   41  C   LIG     1      -9.964  -3.665 -14.556  1.00  0.00           C
HETATM   42  C   LIG     1      -9.568  -2.827 -13.310  1.00  0.00           C
CONECT    1   10   12   33
CONECT    2   31   12   37
CONECT    3   15   13
CONECT    4   17   15    5
CONECT    5    4   14
CONECT    6    7   15
CONECT    7    6    9   18
CONECT    8   28   29   27   26
CONECT    8
CONECT    9    7   23
CONECT   10   11    1   32
CONECT   11   13   10   31
CONECT   12    1    2   34
CONECT   13    3   14   11
CONECT   14    5   13
CONECT   15    6    4    3
CONECT   16   37
CONECT   17    4
CONECT   18    7   22
CONECT   19   23   22   20
CONECT   20   19   21
CONECT   21   20   24
CONECT   22   18   19
CONECT   23    9   19
CONECT   24   21   25
CONECT   25   35   24   29
CONECT   26    8   30   36
CONECT   27    8
CONECT   28    8
CONECT   29   25    8
CONECT   30   26   42
CONECT   31   11    2
CONECT   32   10
CONECT   33    1
CONECT   34   12
CONECT   35   25
CONECT   36   26
CONECT   37    2   16
CONECT   38   40
CONECT   39   40
CONECT   40   38   39   41
CONECT   41   40   42
CONECT   42   41   30
MASTER        0    0    0    0    0    0    0    0   42    0   42    0
END

picture here :

more about the compound at:
PDB Chemical Component 1YI
UPDATE
Last link doesnt work anymore :

Ligand Expo
Ligand Expo is down from January 1-31, 2023. Requests to this service will return a 410 http error code (“Gone”). Please help us better understand which Ligand Expo features are utilized and should be integrated with RCSB.org by taking this survey.

Users can access CCD data in mmCIF format at https://files.wwpdb.org/pub/pdb/refdata/chem_comp/. For other data download during server downtime, please contact chenghua.shao@rcsb.org.

Can use : https://www.ebi.ac.uk/pdbe-srv/pdbechem/chemicalCompound/show/1YI


text:
finding matching motifs in a sequence and their position
I am trying to find some matching motifs on a sequence, as well as the position that the motif is located in and then output that into a fasta file. The code below shows that the motif [L*L*L] is present in the sequence, when I run it returns as "YES" but I do not know where it is positioned.
peptide1="MKFSNEVVHKSMNITEDCSALTGALLKYSTDKSNMNFETLYRDAAVESPQHEVSNESGSTLKEHDYFGLSEVSSSNSSSGKQPEKCCREELNLNESATTLQLGPPAAVKPSGHADGADAHDEGAGPENPAKRPAHHMQQESLADGRKAAAEMGSFKIQRKNILEEFRAMKAQAHMTKSPKPVHTMQHNMHASFSGAQMAFGGAKNNGVKRVFSEAVGGNHIAASGVGVGVREGNDDVSRCEEMNGTEQLDLKVHLPKGMGMARMAPVSGGQNGSAWRNLSFDNMQGPLNPFFRKSLVSKMPVPDGGDSSANASNDCANRKGMVASPSVQPPPAQNQTVGWPPVKNFNKMNTPAPPASTPARACPSVQRKGASTSSSGNLVKIYMDGVPFGRKVDLKTNDSYDKLYSMLEDMFQQYISGQYCGGRSSSSGESHWVASSRKLNFLEGSEYVLIYEDHEGDSMLVGDVPWELFVNAVKRLRIMKGSEQVNLAPKNADPTKVQVAVG"

if re.search(r"L*L*L", peptide1):
    print("YES")
else:
    print("NO")

The code that I wrote to find the position is below, but when I run it says invalid syntax. Could you please assist as I am still new in Bioinformatics, and I have no clue whether in the right track or not.
for position in range(len(peptide1)):
    if peptide1[position:].startswith(re.search(r"L*L*L")):
        print(position+1)

Lastly once the position of these motifs has been identified, how do I highlight them such that they are visible in the output fasta file.
I am searching for LxLxL where x can be any amino acid. I used the * to indicate to python that it should use any amino acid, but I think I am incorrect. Could you please assist on how I should actually be writing that motif?

What is the actual pattern you are looking for? L*L*L* will match anything at all since L* means "0 or more L", so it will match L but it will also match X or anything else.

@terdon, I am searching for LxLxL where x can be any amino acid. I used the * to indicate to python that it should use any amino acid , but I think I am incorrect. Could you please assist on how I should actually be writing that motif.

@terdon If for example I have a motif LLxxxLxPP, could you please advice on how I am supposed to write it inside re.search or re.finditer? this is where I get confused

re.search method do not return a simply boolean but a useful object contains many information:
#! /usr/bin/env python3

import re

peptide1="MKFSNEVVHKSMNITEDCSALTGALLKYSTDKSNMNFETLYRDAAVESPQHEVSNESGSTLKEHDYFGLSEVSSSNSSSGKQPEKCCREELNLNESATTLQLGPPAAVKPSGHADGADAHDEGAGPENPAKRPAHHMQQESLADGRKAAAEMGSFKIQRKNILEEFRAMKAQAHMTKSPKPVHTMQHNMHASFSGAQMAFGGAKNNGVKRVFSEAVGGNHIAASGVGVGVREGNDDVSRCEEMNGTEQLDLKVHLPKGMGMARMAPVSGGQNGSAWRNLSFDNMQGPLNPFFRKSLVSKMPVPDGGDSSANASNDCANRKGMVASPSVQPPPAQNQTVGWPPVKNFNKMNTPAPPASTPARACPSVQRKGASTSSSGNLVKIYMDGVPFGRKVDLKTNDSYDKLYSMLEDMFQQYISGQYCGGRSSSSGESHWVASSRKLNFLEGSEYVLIYEDHEGDSMLVGDVPWELFVNAVKRLRIMKGSEQVNLAPKNADPTKVQVAVG"

re_search_result = re.search(r'L*L*L', peptide1)
if re_search_result:
    print("YES")
    motif_start_position=re_search_result.start()
    motif_end_position=re_search_result.end()
    motif_length = re_search_result.end() - re_search_result.start()
    print("Motif start at position: %d" % motif_start_position)
    print("Motif end at position  : %d" % motif_end_position)
    print("SO length is ......... : %d" % motif_length)
    print("The motif is : <%s>" % peptide1[motif_start_position:motif_end_position])
else:
    print("NO")

Result is:
YES
Motif start at position: 20
Motif end at position  : 488
SO length is ......... : 468
The motif is : <LTGALLKYSTDKSNMNFETLYRDAAVESPQHEVSNESGSTLKEHDYFGLSEVSSSNSSSGKQPEKCCREELNLNESATTLQLGPPAAVKPSGHADGADAHDEGAGPENPAKRPAHHMQQESLADGRKAAAEMGSFKIQRKNILEEFRAMKAQAHMTKSPKPVHTMQHNMHASFSGAQMAFGGAKNNGVKRVFSEAVGGNHIAASGVGVGVREGNDDVSRCEEMNGTEQLDLKVHLPKGMGMARMAPVSGGQNGSAWRNLSFDNMQGPLNPFFRKSLVSKMPVPDGGDSSANASNDCANRKGMVASPSVQPPPAQNQTVGWPPVKNFNKMNTPAPPASTPARACPSVQRKGASTSSSGNLVKIYMDGVPFGRKVDLKTNDSYDKLYSMLEDMFQQYISGQYCGGRSSSSGESHWVASSRKLNFLEGSEYVLIYEDHEGDSMLVGDVPWELFVNAVKRLRIMKGSEQVNL>

The motif L*L*L is a quite strange: L zero or many times followed by L zero or many times followed by followed by a L. It's an equivalent of L*L, equivalent of L+.
Do you search for L.*L.*L, L[^L]*L[^L]*L or simply LLL ?
Notes:

When you post code, use block delimited by two ``` (no by two simple `)
In Python, code indentation is very important, always use spaces (4)

If you want to loop over all motifs (like strings), use findall method:
#! /usr/bin/env python3

import re

peptide1="MKFSNEVVHKSMNITEDCSALTGALLKYSTDKSNMNFETLYRDAAVESPQHEVSNESGSTLKEHDYFGLSEVSSSNSSSGKQPEKCCREELNLNESATTLQLGPPAAVKPSGHADGADAHDEGAGPENPAKRPAHHMQQESLADGRKAAAEMGSFKIQRKNILEEFRAMKAQAHMTKSPKPVHTMQHNMHASFSGAQMAFGGAKNNGVKRVFSEAVGGNHIAASGVGVGVREGNDDVSRCEEMNGTEQLDLKVHLPKGMGMARMAPVSGGQNGSAWRNLSFDNMQGPLNPFFRKSLVSKMPVPDGGDSSANASNDCANRKGMVASPSVQPPPAQNQTVGWPPVKNFNKMNTPAPPASTPARACPSVQRKGASTSSSGNLVKIYMDGVPFGRKVDLKTNDSYDKLYSMLEDMFQQYISGQYCGGRSSSSGESHWVASSRKLNFLEGSEYVLIYEDHEGDSMLVGDVPWELFVNAVKRLRIMKGSEQVNLAPKNADPTKVQVAVG"

re_findall_result = re.findall(r'L[^L]*L[^L]*L', peptide1)
for element in re_findall_result:
    print("Element: <%s>" % element)

Result will be:
Element: <LTGALL>
Element: <LYRDAAVESPQHEVSNESGSTLKEHDYFGL>
Element: <LNLNESATTL>
Element: <LGPPAAVKPSGHADGADAHDEGAGPENPAKRPAHHMQQESLADGRKAAAEMGSFKIQRKNIL>
Element: <LDLKVHL>
Element: <LSFDNMQGPLNPFFRKSL>
Element: <LVKIYMDGVPFGRKVDLKTNDSYDKL>
Element: <LEDMFQQYISGQYCGGRSSSSGESHWVASSRKLNFL>
Element: <LIYEDHEGDSMLVGDVPWEL>

If you want all find objects like with search method, use finditer method:
#! /usr/bin/env python3

import re

peptide1="MKFSNEVVHKSMNITEDCSALTGALLKYSTDKSNMNFETLYRDAAVESPQHEVSNESGSTLKEHDYFGLSEVSSSNSSSGKQPEKCCREELNLNESATTLQLGPPAAVKPSGHADGADAHDEGAGPENPAKRPAHHMQQESLADGRKAAAEMGSFKIQRKNILEEFRAMKAQAHMTKSPKPVHTMQHNMHASFSGAQMAFGGAKNNGVKRVFSEAVGGNHIAASGVGVGVREGNDDVSRCEEMNGTEQLDLKVHLPKGMGMARMAPVSGGQNGSAWRNLSFDNMQGPLNPFFRKSLVSKMPVPDGGDSSANASNDCANRKGMVASPSVQPPPAQNQTVGWPPVKNFNKMNTPAPPASTPARACPSVQRKGASTSSSGNLVKIYMDGVPFGRKVDLKTNDSYDKLYSMLEDMFQQYISGQYCGGRSSSSGESHWVASSRKLNFLEGSEYVLIYEDHEGDSMLVGDVPWELFVNAVKRLRIMKGSEQVNLAPKNADPTKVQVAVG"

re_findall_result = re.finditer(r'L[^L]*L[^L]*L', peptide1)
for element in re_findall_result:
    print("peptide1[%s:%s]=<%s>" % (element.start(), element.end(), peptide1[element.start():element.end()]))

Result will be:
peptide1[20:26]=<LTGALL>
peptide1[39:69]=<LYRDAAVESPQHEVSNESGSTLKEHDYFGL>
peptide1[90:100]=<LNLNESATTL>
peptide1[101:163]=<LGPPAAVKPSGHADGADAHDEGAGPENPAKRPAHHMQQESLADGRKAAAEMGSFKIQRKNIL>
peptide1[248:255]=<LDLKVHL>
peptide1[278:296]=<LSFDNMQGPLNPFFRKSL>
peptide1[378:404]=<LVKIYMDGVPFGRKVDLKTNDSYDKL>
peptide1[407:443]=<LEDMFQQYISGQYCGGRSSSSGESHWVASSRKLNFL>
peptide1[449:469]=<LIYEDHEGDSMLVGDVPWEL>

thank you so much.  I have another question, if for example the motif is written as LLxxLxPP, where x can be any amino acid, how must I write it? I think that is one of the things that confuses me

So if I have LxxLxPP where x can be any amino acid, how must I write it inside the re.finditer ? This is where I get confused

@thole in regular expressions, * means "0 or more of the previous character". To match any character, you want .. Also, [ABC] means "find either A or B or C". So the actual pattern you want is r'L.L.L' but of course that will also match LLLLL.

@terdon thank you, that makes a lot of sense. So if I had a motif LxxLxP then it would be r'L..L.P'

Yes, exactly right @thole.

if you are interested by Regular Expression (RegEx) in Python, see https://www.w3schools.com/python/python_regex.asp