text:
How can I get the coordinates of Hydrogen atoms attached to alpha-carbon?
I am trying to extract the hydrogen atoms attached to the Ca carbon in N no. of proteins. I would be needing the coordinates to calculate its centroid. How should I do it in python ?

OK, tried to figure out a way using Biopython library, I believe there are tens of faster ways, having 1lkk_pdb_test.pdb as input :
ATOM      1  N   LEU A 122      49.672   7.391  11.364  1.00 20.16           N
ANISOU    1  N   LEU A 122     2561   2767   2050    670   -700     42       N
ATOM      2  CA  LEU A 122      50.030   7.561  12.734  1.00 20.25           C
ANISOU    2  CA  LEU A 122     2020   3335   2058   -226   -571    629       C
ATOM      3  C   LEU A 122      48.886   7.480  13.726  1.00 16.62           C
ANISOU    3  C   LEU A 122     1863   2277   1944   -333   -710    706       C
ATOM      4  O   LEU A 122      47.816   8.070  13.517  1.00 15.23           O
ANISOU    4  O   LEU A 122     2132   1556   1887   -122   -845    184       O
ATOM      5  CB  LEU A 122      51.116   8.488  13.074  1.00 21.98           C
ANISOU    5  CB  LEU A 122     1769   3688   2587   -148   -358    491       C
ATOM      6  CG  LEU A 122      51.381   8.978  14.477  1.00 22.50           C
ANISOU    6  CG  LEU A 122     2191   3459   2586   -263   -523    565       C
ATOM      7  CD1 LEU A 122      52.229   8.017  15.285  1.00 23.90           C
ANISOU    7  CD1 LEU A 122     2795   3452   2500   -216   -104    138       C
ATOM      8  CD2 LEU A 122      51.959  10.377  14.502  1.00 26.42           C
ANISOU    8  CD2 LEU A 122     3125   3322   3222    375   -323    398       C
ATOM      9  H1  LEU A 122      48.995   6.786  11.300  1.00 30.24           H
ANISOU    9  H1  LEU A 122     3829   3829   3829      0      0      0       H
ATOM     10  H2  LEU A 122      49.399   8.188  11.020  1.00 30.24           H
ANISOU   10  H2  LEU A 122     3829   3829   3829      0      0      0       H
ATOM     11  H3  LEU A 122      50.392   7.090  10.897  1.00 30.24           H
ANISOU   11  H3  LEU A 122     3829   3829   3829      0      0      0       H
ATOM     12  HA  LEU A 122      50.481   6.686  12.910  1.00 24.30           H
ANISOU   12  HA  LEU A 122     3077   3077   3077      0      0      0       H
ATOM     13  HB2 LEU A 122      51.955   8.066  12.760  1.00 26.37           H
ANISOU   13  HB2 LEU A 122     3339   3339   3339      0      0      0       H
ATOM     14  HB3 LEU A 122      50.982   9.295  12.515  1.00 26.37           H
ANISOU   14  HB3 LEU A 122     3339   3339   3339      0      0      0       H
ATOM     15  HG  LEU A 122      50.493   9.028  14.933  1.00 27.00           H
ANISOU   15  HG  LEU A 122     3419   3419   3419      0      0      0       H
ATOM     16 HD11 LEU A 122      51.771   7.153  15.351  1.00 35.85           H
ANISOU   16 HD11 LEU A 122     4540   4540   4540      0      0      0       H
ATOM     17 HD12 LEU A 122      53.095   7.895  14.844  1.00 35.85           H
ANISOU   17 HD12 LEU A 122     4540   4540   4540      0      0      0       H
ATOM     18 HD13 LEU A 122      52.368   8.382  16.184  1.00 35.85           H
ANISOU   18 HD13 LEU A 122     4540   4540   4540      0      0      0       H
ATOM     19 HD21 LEU A 122      52.806  10.390  14.010  1.00 39.63           H
ANISOU   19 HD21 LEU A 122     5019   5019   5019      0      0      0       H
ATOM     20 HD22 LEU A 122      51.327  10.999  14.083  1.00 39.63           H
ANISOU   20 HD22 LEU A 122     5019   5019   5019      0      0      0       H
ATOM     21 HD23 LEU A 122      52.116  10.649  15.431  1.00 39.63           H
ANISOU   21 HD23 LEU A 122     5019   5019   5019      0      0      0       H
ATOM     22  N   GLU A 123      49.124   6.750  14.812  1.00 14.09           N
ANISOU   22  N   GLU A 123     1166   1937   2055   -326   -424    468       N
ATOM     23  CA  GLU A 123      48.215   6.514  15.890  1.00 12.65           C
ANISOU   23  CA  GLU A 123     1143   1584   1903   -273   -331    -29       C
ATOM     24  C   GLU A 123      46.790   6.330  15.399  1.00 10.88           C
ANISOU   24  C   GLU A 123     1138   1238   1605   -155   -376    -53       C
ATOM     25  O   GLU A 123      45.875   7.059  15.809  1.00 10.93           O
ANISOU   25  O   GLU A 123     1300   1048   1652     37   -173   -144       O
ATOM     26  CB  GLU A 123      48.256   7.537  17.028  1.00 14.54           C
ANISOU   26  CB  GLU A 123     1615   1778   1930   -299   -450   -280       C
ATOM     27  CG  GLU A 123      47.696   6.998  18.341  1.00 19.07           C
ANISOU   27  CG  GLU A 123     2450   2721   1808   -148   -132   -221       C
ATOM     28  CD  GLU A 123      47.620   7.948  19.490  1.00 21.48           C
ANISOU   28  CD  GLU A 123     3199   2805   1858    113    106      6       C
ATOM     29  OE1 GLU A 123      47.461   9.187  19.266  1.00 25.80           O
ANISOU   29  OE1 GLU A 123     4100   2736   2608   1002   -257     -4       O
ATOM     30  OE2 GLU A 123      47.491   7.497  20.672  1.00 22.45           O
ANISOU   30  OE2 GLU A 123     3024   3373   1819    545    690    367       O
ATOM     31  H   GLU A 123      49.919   6.375  14.873  1.00 16.91           H
ANISOU   31  H   GLU A 123     2141   2141   2141      0      0      0       H
ATOM     32  HA  GLU A 123      48.484   5.640  16.294  1.00 15.18           H
ANISOU   32  HA  GLU A 123     1922   1922   1922      0      0      0       H
ATOM     33  HB2 GLU A 123      49.194   7.821  17.171  1.00 17.45           H
ANISOU   33  HB2 GLU A 123     2210   2210   2210      0      0      0       H
ATOM     34  HB3 GLU A 123      47.737   8.336  16.761  1.00 17.45           H
ANISOU   34  HB3 GLU A 123     2210   2210   2210      0      0      0       H
ATOM     35  HG2 GLU A 123      46.784   6.654  18.167  1.00 22.88           H
ANISOU   35  HG2 GLU A 123     2897   2897   2897      0      0      0       H
ATOM     36  HG3 GLU A 123      48.252   6.227  18.618  1.00 22.88           H
ANISOU   36  HG3 GLU A 123     2897   2897   2897      0      0      0       H
ATOM     37  N   PRO A 124      46.576   5.340  14.522  1.00  9.86           N
ANISOU   37  N   PRO A 124      893   1140   1575     13   -285      8       N
ATOM     38  CA  PRO A 124      45.182   5.138  14.041  1.00  8.89           C
ANISOU   38  CA  PRO A 124      875   1151   1229     -6   -227     93       C
ATOM     39  C   PRO A 124      44.244   4.845  15.208  1.00  8.67           C
ANISOU   39  C   PRO A 124      858   1071   1245    -81   -189    -53       C
ATOM     40  O   PRO A 124      44.593   4.143  16.165  1.00  8.98           O
ANISOU   40  O   PRO A 124     1088   1065   1134    -47   -410    127       O
ATOM     41  CB  PRO A 124      45.309   3.946  13.121  1.00  9.41           C
ANISOU   41  CB  PRO A 124     1048   1085   1312     99     -1    -37       C
ATOM     42  CG  PRO A 124      46.555   3.221  13.556  1.00  9.55           C
ANISOU   42  CG  PRO A 124     1054   1161   1282     62   -243      6       C
ATOM     43  CD  PRO A 124      47.489   4.326  13.988  1.00 10.06           C
ANISOU   43  CD  PRO A 124      911   1179   1591    122   -319   -125       C
ATOM     44  HA  PRO A 124      44.868   5.939  13.533  1.00 10.67           H
ANISOU   44  HA  PRO A 124     1351   1351   1351      0      0      0       H
ATOM     45  HB2 PRO A 124      44.517   3.357  13.204  1.00 11.29           H
ANISOU   45  HB2 PRO A 124     1429   1429   1429      0      0      0       H
ATOM     46  HB3 PRO A 124      45.388   4.241  12.180  1.00 11.29           H
ANISOU   46  HB3 PRO A 124     1429   1429   1429      0      0      0       H
ATOM     47  HG2 PRO A 124      46.363   2.605  14.307  1.00 11.46           H
ANISOU   47  HG2 PRO A 124     1451   1451   1451      0      0      0       H
ATOM     48  HG3 PRO A 124      46.943   2.702  12.807  1.00 11.46           H
ANISOU   48  HG3 PRO A 124     1451   1451   1451      0      0      0       H
ATOM     49  HD2 PRO A 124      48.119   4.010  14.683  1.00 12.07           H
ANISOU   49  HD2 PRO A 124     1528   1528   1528      0      0      0       H
ATOM     50  HD3 PRO A 124      48.005   4.678  13.220  1.00 12.07           H
ANISOU   50  HD3 PRO A 124     1528   1528   1528      0      0      0       H
END

my code is:
import Bio

from Bio.PDB import PDBParser

print('imported Biopython module vesion : ' , Bio.__version__)

parser = PDBParser (QUIET = True)

structure = parser.get_structure('test',  '1lkk_pdb_test.pdb')

atoms = structure.get_atoms()

coordz = []

for i in [i for i in atoms]:

    if i.get_id() == 'HA' :

        print(i.get_parent(), i.get_parent().id, i, i.coord)

        coordz.append(i.coord)

print('coordz list : ')

for i in coordz :
    print(i , type(i))

output:
imported Biopython module vesion :  1.80
imported Biopython module vesion :  1.80
<Residue LEU het=  resseq=122 icode= > (' ', 122, ' ') <Atom HA> [50.481  6.686 12.91 ]
<Residue GLU het=  resseq=123 icode= > (' ', 123, ' ') <Atom HA> [48.484  5.64  16.294]
<Residue PRO het=  resseq=124 icode= > (' ', 124, ' ') <Atom HA> [44.868  5.939 13.533]
coordz list :
[50.481  6.686 12.91 ] <class 'numpy.ndarray'>
[48.484  5.64  16.294] <class 'numpy.ndarray'>
[44.868  5.939 13.533] <class 'numpy.ndarray'>

as per the reason HA represent Cα hydrogen, please refer to _

Westbrook, John D et al. “The chemical component dictionary: complete descriptions of constituent molecules
in experimentally determined 3D macromolecules in the Protein Data Bank.”
Bioinformatics (Oxford, England) vol. 31,8 (2015): 1274-8. doi:10.1093/bioinformatics/btu789

httpswww.ncbi.nlm.nih.govpmcarticlesPMC4393513pdfbtu789.pdf
Chemical Component Dictionary
as explamples pages for Leucine links to :
https://www.ebi.ac.uk/pdbe-srv/pdbechem/chemicalCompound/show/LEU
and http://ligand-expo.rcsb.org/reports/L/LEU/LEU_D3L3.gif
lot of trouble to find a diagram:

Just loop over all your proteins in a similar way and store the coords as you need/wish, remember coordinates are retrieved as numpy arrays.

I would import the PDB into the molecular dynamics package GROMACS and use the built-in utilities to calculate the centroid (https://manual.gromacs.org/current/onlinehelp/gmx-rms.html).  This is much simpler than trying to write your own utility for this.  Other pre-written calculations to perform on your PDB files:
https://manual.gromacs.org/current/user-guide/cmdline.html


text:
Mismatch between mutant gene and reference gene outside the site of mutation
I am tasked with designing primers for a particular mutant (Target Gene Locus: At1g28490 i.e. SYP61) of Arabidopsis, obtained from Gabi-Kat:

The NCBI reference sequence for the genome goes:
GGCTTTATAAGATTCTCCAGTAGCCCAAAAATAAACTTTCAGCGAATCTCAAAGCCAACCGCAGGTGATGCGTTTTGCTATGTAAGAGAGACGAATCTCACGCCGGAATCTCTTTTACTTCACCGTCGTTTGTTCTCCGTCGCGTCACCTTTCGATTAATTCCTTCACAGGGTAGCAAAATTTTGTTGATTTCGGG$\color{red}{\textrm{[}}$ATGTCTTCAGCTCAAGATCCATTCTACATTGTTAAAGAAGAGATTCAGGATTCTATTGATAAGTTGCAATCTACATTCCACAAATGGGAGCGTATCTCTCCCGACATGGGTGATCAAGCACATGTGGCGAAGGAGCTTGTTGCTACTTGTGGAAGCATTGAGTGGCAGGTA͟G͟A͟T͟G͟A͟G͟C͟T͟G͟G͟A͟A͟A͟A͟A͟G͟C͟G͟A͟T͟T͟A͟C͟T͟G͟T͟T͟G͟C͟A͟G͟C͟A͟A͟A͟A͟G͟A͟T͟C͟C͟T͟T͟C͟ATGGTATGGCATTGATGAAGCTGAGCTTGAAAAACGGAGGAGATGGACTAGTAATGCTAGGACACAGGTGCGGAATGTGAAGTCTGGTGTTTTAGCTGGAAAGGTTAGTTCTGGAGCAGGTCATGCAAGTGAAGTGCGTAGAGAACTAATGAGAATGCCAAACTCGGGTGAAGCAAGTAGATATGATCAGTATGGAGGGAGAGACGATGATGGTTTTGTACAGTCAGAATCAGATAGGCAAATGCTGTTGATAAAGCAACAAGACGAAGAATTGGATGAGCTTAGTAAAAGTGTACAGAGAATTGGAGGAGTGGGGCTCACTATACACGACGAACTCGTTGCACAGGAGAGAATAATAGATGAATTGGATACGGAGATGGACAGCACAAAGAACCGTCTTGAATTTGTTCAGAAAAAAGTGGGAATGGTAATGAAGAAAGCTGGAGCAAAAGGTCAAATGATGATGATATGTTTCTTGCTCGTCTTGTTCATCATCCTATTCGTTCTTGTCTTCTTGACCTAA$\color{red}{\textrm{]}}$AGAACAGAGCCGCGGGATCTTATAGATGGGAGTTTAATGTCTTACATTATCATGTTTGTGGTGAAACAAATGGAATATTTACATAAGGTTGTCTGTTTTTTCTTCTCCTATCCAATCCCACATATAATTTTGAACTTTGGCTGTGTGATCTCTAGAGAACCTGAAATTAAATACTTTCATAGCTTGATCCCATGGCTAAGTTCCT
Gabi-Kat provides the following sequence from the mutant as a confirmatory match with the genome locus:
AGCGGATCTCTCCCGACATGGCTGATCAAGCACATGTGGCGAAGGAGCTTGTTGCTACTTGTGGAAGCATTGAGTGGCAGGTC͟C͟T͟C͟C͟T͟A͟T͟T͟C͟A͟C͟C͟C͟A͟C͟G͟G͟T͟T͟C͟T͟A͟T͟T͟A͟A͟T͟C͟C͟C͟C͟T͟C͟T͟T͟T͟A͟T͟A͟T͟A͟T͟A͟A͟

I have bold-faced the common region to both sequences. What explains the complete mismatch of the trailing part (u͟n͟d͟e͟r͟l͟i͟n͟e͟d͟) of the Gabi-Kat sequence with the ref genome?
It is not the case that the mutation happened there; we are in the coding region (demarcated by red square brackets in the NCBI sequence) while Gabi-Kat's mutation is in an untranslated region.
Update 1:
Running Ensembl Blast on the GabiKat seq returns a 92.9% identicality with SYP61. Their sequence for the gene (taken from TAIR10.1) goes:
AGCGGATCTCTCCCGACATGGCTGATCAAGCACATGTGGCGAAGGAGCTTGTTGCTACTTGTGGAAGCATTGAGTGGCAGGTA͟C͟T͟C͟G͟T͟A͟T͟T͟T͟A͟G͟T͟C͟A͟C͟G͟G͟T͟T͟C͟T͟A͟T͟T͟A͟A͟T͟C͟T͟C͟G͟T͟C͟T͟T͟T͟A͟T͟A͟T͟A͟T͟A͟A͟
Note the difference in the trailing part from NCBI ref; this almost coincides with the Gabi-Kat seq!
So, is the NCBI ref-genome and the data for the particular TAIR genome inaccurate? Maybe, TAIR updated their whole genome set but forgot to update data for individual genes?
Update 2:
My proposed explanation for the anomaly makes no sense. TAIR 7 (2004) has the same sequence. It looks like I am missing out on something more fundamental.

Conclusion The evidence suggests the region you are looking at Gabi-kat appears either 27045 downstream or 192250 upstream of syntaxin SNARE (SYP61), but is definitely not within the gene. It is chance homology (e.g. a pseudo-gene) that appears to occur within syntaxin (SYP61). The intergenic region could be much closer to the SYP61 gene if the same reference genome was used for both analyses Genbank LR782542.1 versus Genbank NC_003070.9 (different reference genomes in my analysis). Hence, the University is probably correct.
What I suggest is the genome of A. thaliana is downloaded and Blast analysis is performed locally. This way differences in the reference genome are removed.
The analysis is as follows:
The target protein in question is involved in vesicle trafficking (here) called SNARE, comprising coiled-coil domains and are highly conserved, for example show 100% conservation between Arabidopsis thaliana, Arabidopsis suecica and   Camelina sativa

MSSAQDPFYIVKEEIQDSIDKLQSTFHKWERISPDMGDQAHVAKELVATCGSIEWQVDELEKAITVAAKDP

If Gabi-kat was a protein its sequence would be ...

RISPDMADQAHVAKELVATCGSIEWQVLLFTHGSINPLFIY

Here's a blast of the whole regions of the Gabi-kat. The region in green is the SNARE proteins ...

So the issue is what is the terminal 3' protein and Blast that, i.e. the final bit that matches no known plant and the locus you flagged as being anomalous ..

This is basically zero result via Blast, the protein with low homology are lower metazoans (animals) and a few bacteria ... a long way from plants. In other words this isn't a known protein of natural origin at least on NCBI.
The locus is literally like no protein on earth.
However ... Blasting the nucleotide ...

Reference Arabidopsis thaliana = syntaxin (SNARE; SYP61)

Gabi-Kat Arabidopsis thaliana top two Blast hits are ...

Chromosome 1 at 10039082 - 10039208, or 9819787 to 981991 or
Both value    1e-42
Both dentity is 93%

So Gabi-Kat is a real part of the Arabidopsis thaliana, Gabi-kat is a non-coding region. Thats the major conclusion.
The question where is the location of SYP61? And that is

Chromosome 1 at 10012037 to 10020311

Described here from NCBI's full report

Thus the mutant is 27045bp downstream or 192250 bp upstream IF the same reference genome is used for the blasting (Genome LR782542.1) and the Genbank report NC_003070.9. They could easily be much closer to each other is the same reference genome was used for both sets of analysis.

To answer the comments

I [the OP] do not follow your [my] blast analysis
My [the OP] blasts gives ... ['hits within the gene']
etc ....

It appears to me there is a psuedo-gene phenomenon occurring which is making you think a mutant is located within a gene. However, if you disagree thats perfectly okay and the original version of this answer I had assumed the Gaba-kat team had got this wrong (I always from the basis of supporting the OP).
NCBI was the database repository in the remote analysis, but there's an issue with multiple references ... so to reiterate I suggest you perform this analysis as a local blast to remove multiple reference phenomena.
To conclude I personally don't think there is a basis to believe the University responsible for Gaba-kat has got this wrong https://www.gabi-kat.de. However, if you feel otherwise then I would recommend contacting them with your analysis and any supportive analysis for example from this thread or elsewhere. Please note, I don't have a connection with this university.

I do not follow your BLAST of the Gabi-Kat Nucleotide Seq. My blast gives a 92.9% alignment for the entire stretch (127 bp long) with SYP61 (1:10016601 to 1:10016721). E = 4.6e-45

Indeed, SYP61 is located from 1:10015856 to 10018047 in NC_003070.9 (GCA_000001735.2; TAIR v. 10.1). Which is inclusive of 1:10016601 to 1:10016721. You have got the specifics wrong.

Running a nucleotide blast on the underlined part of the NCBI ref-seq (that is the part, which diverged from Gabi-Kat seq; 45 bp) gives 100% match to SYP61 (1:10016871 to 1:10016915). E = 3.5e-18.

To conclude, the 45 bp stretch is non-contiguous to the 127 bp stretch in NC_003070.9 (GCA_000001735.2; TAIR v. 10.1).

source: Stack Exchange

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

sort: votes


text:
Availability of information on genes in Gnomad VCF data
I was wondering does the gnomad genome data that is downloaded in the vcf format on variants contain information on what is the nearest gene and is the genomic location available exome/intron?
if yes then what are the columns that contain this data? Is it the INFO column?
I have a list of SNPs  with chromosome position and the reference and alternative allele information and for each, I would like to know what is he nearest gene and the location:exome/intron.
example data

We would need to see all the columns of the gnomad VCF, see e.g. here the mtDNA VCF (the smallest file), the information you are looking for should be described in the header. As stated on the other question, if you want to know the nearest gene, the intersect operation via bedtools is pretty trivial to combine this VCF with a human genome GFF3 file (of the appropriate genome build).

This question has gone through some iterations and commentary, so I've decided to just answer it, though it's more of a "homework" question. What you need to do is to look at the VCF file, the main thing you want is already in there.
Go to the VCF of interest. I chose the mtDNA because it's smaller and more convenient. You can zless into the bgz file to see what it looks like, but if you want uncompressed you can run this (head -1000 just to be a tractable example):
gzip -dc gnomad.genomes.v3.1.sites.chrM.vcf.bgz | grep -v intergenic | head -1000 > mt.vcf

I grep -v intergenic because one of the first things that pops up is that most of the first variants in the file are specifically annotated as intergenic_variant, so that answers part of your question. For those variants you can just run a bedtools closest against the right GFF3 version, if you think the nearest gene will be informative.
For variants in genes, what you want is already in the VCF. See the entry for this variant:
chrM    1492    .       A       C       .       npg     filters=npg;variant_colla
psed=A1492C;vep=C|non_coding_transcript_exon_variant|MODIFIER|MT-RNR1|ENSG0000021
1459|Transcript|ENST00000389680|Mt_rRNA|1/1||ENST00000389680.2:n.845A>C||845|||||
1||1|SNV||HGNC|HGNC:7470|YES||||||||1||||||||||||;base_qual_hist=0|0|0|0|0|0|0|0|
0|0;position_hist=0|0|0|0|0|0|0|0|0|0;strand_bias_hist=0|0|0|0|0|0|0|0|0|0;weak_e
vidence_hist=0|0|0|0|0|0|0|0|0|0;contamination_hist=0|0|0|0|0|0|0|0|0|0;heteropla
smy_below_min_het_threshold_hist=2|0|0|0|0|0|0|0|0|0;excluded_AC=2;AN=56432;AC_ho
m=0;AC_het=0;hl_hist=0|0|0|0|0|0|0|0|0|0;dp_mean=2976.53;mq_mean=.;tlod_mean=.;AF
_hom=0.00000;AF_het=0.00000;max_hl=0.00000;hap_AN=2680|1537|868|603|34|282|91|147
83|701|934|3144|2732|663|2977|4724|5672|126|1|1298|366|7|393|3079|6037|1234|819|5
46|12|89;hap_AC_het=0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0;hap
_AC_hom=0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0;hap_AF_hom=0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00;hap_AF_het=0.00000e+00|0.00000e+00|0.00 ...
# this is of course a very long line!!! it goes on but all the interesting stuff already happened.

You can see that this variant lies within the RNR1 rRNA gene of the mitochondrial genome (MT-RNR1 in the file). Granted it mangles it a little bit by describing a variant in a non-protein-coding gene as a non_coding_transcript_exon_variant, but that's the kind of thing that you can work around.
Indeed, if you look at the header of the file, you can see that quite extensive annotations are included (where available):
##INFO=<ID=vep,Number=.,Type=String,Description="Consequence annotations from Ensembl VEP; note that the SINGLE_EXON flag and END_TRUNC filters have been removed from the LOFTEE annotations to avoid misinterpretation in context of the mitochondrial genome. Format: Allele|Consequence|IMPACT|SYMBOL|Gene|Feature_type|Feature|BIOTYPE|EXON|INTRON|HGVSc|HGVSp|cDNA_position|CDS_position|Protein_position|Amino_acids|Codons|ALLELE_NUM|DISTANCE|STRAND|VARIANT_CLASS|MINIMISED|SYMBOL_SOURCE|HGNC_ID|CANONICAL|TSL|APPRIS|CCDS|ENSP|SWISSPROT|TREMBL|UNIPARC|GENE_PHENO|SIFT|PolyPhen|DOMAINS|HGVS_OFFSET|MOTIF_NAME|MOTIF_POS|HIGH_INF_POS|MOTIF_SCORE_CHANGE|LoF|LoF_filter|LoF_flags|LoF_info">

In conclusion, look at your data!

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


text:
How to run a nextflow process for each file generated by another process separately - Input tuple does not match input set cardinality error
I have a process called FINDTAIL that generates different number of files depending on the input data. Its either 2 files or four files i.e.

1. read_1{,.pr,.sl}.fasta and read_2{,.pr,.sl}.fasta

or

2. read_1.pr.fasta, read_1.sl.fasta, read_2.pr.fasta and read_2.sl.fasta

/*
 * the `FINDTAIL` process takes "${pair_id}_strand_check.txt" file
 * and executes commands according to the reads strandedness.
 * The output generated is FASTA files.
 */

process FINDTAIL {

    tag { pair_id }

    debug true

    input:
    path strand_check_file
    tuple val(pair_id), path(reads)

    output:
    tuple val(pair_id), path("tail_${pair_id}_{1,2}{,.pr,.sl}.fasta")

    script:
    """
#!/usr/bin/env python3

import subprocess
import pandas as pd
import os

result = pd.read_csv("${strand_check_file}", sep="\\r\\n", header=None, engine='python')

failed = float(result.iloc[1,0].replace('Fraction of reads failed to determine: ', ''))
fwd = float(result.iloc[2,0].replace('Fraction of reads explained by "1++,1--,2+-,2-+": ', ''))
rev = float(result.iloc[3,0].replace('Fraction of reads explained by "1+-,1-+,2++,2--": ', ''))
fwd_percent = fwd/(fwd+rev)
rev_percent = rev/(fwd+rev)

if float(result.iloc[1,0].replace('Fraction of reads failed to determine: ', '')) > 0.50:
    cmd1 = "findtail_v1.01 --input_file " + "${reads[0]}" + " --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > " + "tail_${pair_id}_1.sl.fasta"
    cmd2 = "findtail_v1.01 --input_file " + "${reads[1]}" + " --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > " + "tail_${pair_id}_2.sl.fasta"
    print(cmd1)
    subprocess.call(cmd1, shell=True)
    print(cmd2)
    subprocess.call(cmd2, shell=True)

if fwd_percent > 0.9:
    #Over 90% of reads explained by "1++,1--,2+-,2-+
    #Data is likely FR/fr-secondstrand
    cmd1 = "findtail_v1.01 --input_file " + "${reads[0]}" + " --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type pr --d > " + "tail_${pair_id}_1.pr.fasta"
    cmd2 = "findtail_v1.01 --input_file " + "${reads[1]}" + " --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type pr --d > " + "tail_${pair_id}_2.pr.fasta"
    print(cmd1)
    subprocess.call(cmd1, shell=True)
    print(cmd2)
    subprocess.call(cmd2, shell=True)

elif rev_percent > 0.9:
    # Over 90% of reads explained by "1+-,1-+,2++,2--
    # Data is likely RF/fr-firststrand
    cmd1 = "findtail_v1.01 --input_file " + "${reads[0]}" + " --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > " + "tail_${pair_id}_1.sl.fasta"
    cmd2 = "findtail_v1.01 --input_file " + "${reads[1]}" + " --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > " + "tail_${pair_id}_2.sl.fasta"
    print(cmd1)
    subprocess.call(cmd1, shell=True)
    print(cmd2)
    subprocess.call(cmd2, shell=True)

elif max(fwd_percent, rev_percent) < 0.6:
    #Under 60% of reads explained by one direction
    #Data is likely unstranded
    cmd1 = "findtail_v1.01 --input_file " + "${reads[0]}" + " --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > " + "tail_${pair_id}_1.sl.fasta"
    cmd2 = "findtail_v1.01 --input_file " + "${reads[0]}" + " --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type pr --d > " + "tail_${pair_id}_1.pr.fasta"
    cmd3 = "findtail_v1.01 --input_file " + "${reads[1]}" + " --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > " + "tail_${pair_id}_2.sl.fasta"
    cmd4 = "findtail_v1.01 --input_file " + "${reads[1]}" + " --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type pr --d > " + "tail_${pair_id}_2.pr.fasta"
    print(cmd1)
    subprocess.call(cmd1, shell=True)
    print(cmd2)
    subprocess.call(cmd2, shell=True)
    print(cmd3)
    subprocess.call(cmd3, shell=True)
    print(cmd4)
    subprocess.call(cmd4, shell=True)

else:
    #if strand couldnt be detected
    # we assume it is first-stranded
    # since most illumina reads are
    cmd1 = "findtail_v1.01 --input_file " + "${reads[0]}" + " --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > " + "tail_${pair_id}_1.sl.fasta"
    cmd2 = "findtail_v1.01 --input_file " + "${reads[1]}" + " --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > " + "tail_${pair_id}_2.sl.fasta"
    print(cmd1)
    subprocess.call(cmd1, shell=True)
    print(cmd2)
    subprocess.call(cmd2, shell=True)
    """
}

I want to pass these fasta file outputs of the FINDTAIL process into an ALIGN process that runs on each file separately.
The ALIGN process:
/*
 * define the `ALIGN` process that outputs unmapped files in FASTA format
 * It will take the FASTA files creates by the `FINDTAIL` process
 * and try to align them INDIVIDUALLY.
 */

fasta_file_ch = find_tail_ch.flatten()

process ALIGN {

    tag { pair_id }

    input:
    path index_files, stageAs: "bt2_index/*"
    tuple val(pair_id)
    file(fasta_file) from fasta_file_ch

    output:
    tuple val(pair_id), path("unmapped_tail_${pair_id}_{1,2}{,.pr,.sl}.fasta")

    script:
    def idx = index_files[0].getBaseName(2)

    """
    bowtie2 -x "${idx}" -f "${fasta_file}" --un "unmapped_tail_${pair_id}/${it.baseName}.fasta"
    """
}

I have tried using the flatten() function to format my find_tail_ch but it still does not work.
Best wishes~

If the output channel produces a tuple where the first element is a simple value and the second element is a collection of files, like in this case, you can simply use the transpose operator to produce a transposition where the first element in each resulting tuple is the pair_id and the second is a FASTA file. You'll then need to update your process definition to receive this tuple, for example:
process ALIGN {

    tag { pair_id }

    input:
    path index_files, stageAs: "bt2_index/*"
    tuple val(pair_id), path(fasta_file)

    output:
    tuple val(pair_id), path("${fasta_file.baseName}.unmapped.fasta")

    script:
    def idx = index_files[0].getBaseName(2)

    """
    bowtie2 \\
        -x "${idx}" \\
        -f "${fasta_file}" \\
        --un "${fasta_file.baseName}.unmapped.fasta"
    """
}

Your workflow block might look like:
workflow {

    bt2_index = ...

    ...

    align_input_ch = find_tail_ch.transpose()

    ALIGN( bt2_index, align_input_ch )
}

Thank you very very much Steve!

If you could clarify one thing: since we have used the transpose function for align_input_ch, if I were to pass down the contents i.e. two separate fasta files, of the ALIGN process to another process i.e. file by file separately, will it automatically pass separately?

@pubsurfted Yes, if I understand correctly. But if you're ever in doubt, you can use the view() operator to see what's going on. In the above example, you could try: ALIGN.out.view(). I think you should be able to see either two or four process with the sample pair_id in the first element of each tuple. Be sure to add -ansi-log false to your command line options.

Thanks again! I'm so grateful for your help.


text:
Using R `ggtree`, attach extra data next to the nodes on an illustrated tree/newick file
This is really a programming question but ggtree might be too obscure for this question to be posted on StackOverFlow as a general R related question.
I am having trouble extrapolating this this tutorial and the %<+% operator in ggtree.
I have a Newick file and a table that has a number for every tip and node in the Newick file. I want to produce a final graphic that shows a tree where the numbers from the table are shown next to respective tips and nodes.
What I want to do is close to the tutorial, I want to take a tree object and use a data-frame where one column is node-names and the second column holds numerical information per node. I want to attach the information from the data-frame to the tree object.
The data-frame looks something like this:

Node / tips
Count

A
12

B
14

So, I want to ultimately display the number 12 next to node A and number 14 next to node B while keeping the names intact. Now I have the code working up to the following point:
(I am really only just working my way through learning R's grammar of graphics and R in general, so assume I am doing something wrong.)
library(ggplot2)
library(ggtree)

x <- read.tree("tree.txt")
info <- read.csv("counts.csv",sep = '\t',header=TRUE)

p <- ggtree(x) %<+% info + xlim(-.1, 5)

p2 <- p + geom_tiplab(offset = 0.2) + geom_tippoint(aes(caption = Increase)) + theme_tree2()+geom_tiplab(align=TRUE, linesize=.5)+ xlim(0,25)

ggsave("test_2.pdf", width = 60, height = 80, units = "cm", limitsize = FALSE)

print(p2)

?geom_tippoint::aes() # is there a way to list all the option inside geom_tippoint::aes()

The above code overlays the Newick_label column with the node names instead of attaching the value in the second column.
But specifically, I am at a loss as to how I can list the options available within geom_tippoint(aes()). I can use ?geom_tippoint to access the documentation of geom_tippoint but I think my task needs to manipulate options within the aes() function and the the geom_tippoint() documentation does not list all the aes() options available. And am stuck trying different keywords for geom_tippoint(aes()).
Can that be done in ggtree or is this specific application of ggtree obscure enough that there is no way to do this and produce a legible graphic? (eg issue: the graphic getting mushy overlays). I have hundred of (not uniform) genes trees so splitting data for each layer (node and tips) might be not the best way.
The files if anyone wants to try: counts and tree
Clarification on the labels for nodes and tips, I want to label both if possible. The data contains information for both.
Is annotating tips AND nodes out of the ordinary? Tips and nodes are technically different but I have quantitative data at both levels and am looking for the option to annotate at both levels.

The answer is geom_tiplab()  if I get the chance this week I will post a response on a hypothetical data set, but can't do this today. Oh wait, maybe not, there is still the tip versus node issue - those are different solutions - might be easier just to copy part of the tree file to see what 'A' and 'B' are. At present the table is called node/tips: a tip is a subset of the nodes.

I will look into geom_tiplab() for now, thanks!

There's nothing 'wonky' with any procedure, it's just I don't quite know what you are looking for.  Main point, geom_tiplab() only works on tips and the tutorial (and the code) you provided is just for tips. It's different code for nodes and the labels are inserted at different places. Naming internal nodes is not straight forward and I'm not sure how you are doing that. It could be you assume what I would call the tip, taxa or terminal node are what you are calling a node and my taxa label or tip label are what you call a tip. If that is the case thats easy.

I ended up using the following code to solve this issue:
Change the tsv file and the tree file name accordingly. Be mindful of the format the tsv file should be in (indicated in the code).
Also change the tree and the font scales as needed.
# Load required libraries
library(ape)
library(ggtree)

# Read the TSV file with species and node names and numbers
# Make sure the TSV file has headers: "name", and "number"
data <- read.table("counts.tsv", header = TRUE, sep = "\t", fill = TRUE, colClasses = "character")

# Read the Newick file
newick_file <- "species_tree.newick"
tree <- read.tree(newick_file)

# Add labels and numbers to the tree
add_labels <- function(tree, data) {
  missing_numbers <- c()

  for (i in 1:nrow(data)) {
    label_name <- data$name[i]
label_number <- data$number[i]
    if (label_name %in% tree$tip.label) {
  tip_index <- match(label_name, tree$tip.label)
  tree$tip.label[tip_index] <- paste(label_name, label_number, sep = ":")
} else if (label_name %in% tree$node.label) {
      node_index <- match(label_name, tree$node.label)
  tree$node.label[node_index] <- paste(label_name, label_number, sep = ":")
    } else {
      missing_numbers <- c(missing_numbers, label_name)
    }
  }

  return(list(tree = tree, missing_numbers = missing_numbers))
}

# Modify the tree object to include labels and numbers
result <- add_labels(tree, data)
tree <- result$tree
missing_numbers <- result$missing_numbers

# Write missing labels to a log file
if (length(missing_numbers) > 0) {
  writeLines(c("Missing numbers for the following labels:", missing_numbers), "missing_labels.log")
}

# Convert factors to characters
tree$tip.label <- as.character(tree$tip.label)
tree$node.label <- as.character(tree$node.label)

# Visualize the tree
p <- ggtree(tree, layout = "rectangular") +
  geom_tiplab(size = 3) +
  geom_nodelab(size = 3, na.rm = TRUE, nudge_x = 0.008) +
  theme_tree2() +
  theme(legend.position = "none")

# Save the tree as a PDF
ggsave("species_tree_with_numbers.pdf", p, width = 20, height = 20, dpi = 300)
```


text:
What is the correct method of identifying the interacting residues for specific molecular docking?
I'm trying to do a specific molecular docking using Autodock Vina in PyRx. So, one way to do that is to mark all of the interacting residues and make sure the grid box encompasses all of that interacting residues.
So, my protein already has a co-crystalized ligand attached to it. However, I'm not sure what method to use.
After searching online, I found out that I should identify the interacting residues of that ligand. Then, I have opened this protein-ligand complex in 2 applications and 1 web server.
Here are the information of the protein-ligand complex:

The protein PDB ID is: 5TL8 (Naegleria fowleri CYP51-posaconazole complex)
The ligand is Posaconazole
The cofactor and other small molecules: HEME, 1,2-ETHANEDIOL, calcium ion

I tried opening the complex in Discovery Studio
Ligplus
and Protein Plus webserver (this server has more interacting residues). All of them have different residues.
Which one should I follow? Should we included the HEME or do we just look at the Posaconazole itself?

have a look at Rupp paper too as alwayis in Sciences: "trust, but verify" https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4531100/pdf/nihms662686.pdf

@pippo1980 hello, thank you very much.

I think I've found the answer to my own question.
So, I tried using this PRANK web server for predicting the binding site of my protein. The results is suprisingly very similar to the one from the Protein Plus web server. So, I think it's better to use the one from Protein Plus.
You can see the results on the link below.
References:

https://prankweb.cz/viewer?id=5tl8&database=v3-conservation-hmm&created=2022-01-27T07:50:59&lastChange=2022-01-27T07:50:59%20%20&structureName=structure.pdb&predictionName=5TL8&predictedStructure=undefined
https://proteins.plus/5tl8

Not an answer either. But a series of caveats.
First, HEM parameterisation. For molecular thermodynamics tools you will need to generate a topology file for this bad-boi. The iron has 2 donor bonds and 4 single bonds, one of which with the coordinating cysteine and the other with an oxygen, which will attach the substrate, in this case a steroid —let's ignore the radical state of the activated bound oxygen, but do keep in mind that the partial charges and atom type auto-assigned will be wrong. A cheat would be constraining said molecular oxygen and using a porphyrin ring with no iron: not reality, but would make the experiment work, is as physically questionable as any iron-bound topology you'll get and the oxygen is the key player here.
Second, cytochromes act on greasy substrates, there are few if any polar interactions generally, so constraining to polymer residues may not be that useful. Specificity is controlled by the gating and entranceway (think of trucks reversing into a loading bay), which is annoying the gating loops are floppy loops that are generally not crystallised, so tethered MD from a bound to unbound state (undocking) is sometimes done.
There's a lot of literature on P450 and P411 cytochromes, but often it's not as clear cut and simple as say a hydrolase.
However, if your substrate you want placed is similar to one present you may be interested in a tool of mine Fragmenstein, which can be used to convert a compound to a similar one (stripping the haem and adding it later: I had a Rosetta params file for HEM somewhere but not sure where it is nowadays)

Thank you for your comments, and that's a pretty cool tool. I will look more into it.

not an answer, just more inputs, but https://proteins.plus/5tl8#poseview
gives me this:

references:
ProteinsPlus: interactive analysis of protein–ligand binding interfaces
Molecular complexes at a glance: automated generation of two-dimensional complex diagrams
Drawing the PDB: Protein−Ligand Complexes in Two Dimensions
from this last paper :

Hydrogen bonds are implemented mainly following the measures published by Desiraju and Steiner in 2001.14 The optimal distance between two atoms connected by a hydrogen bond is set to 1.9 Å with a tolerance of 0.5 Å. Additional to this measure, the acceptor−hydrogen−donor angle must not fall below 120°. Hydrogen atoms that are bound to a noncarbon atom are treated as hydrogen donor candidates. Potential acceptor atoms are either nitrogen, oxygen, or sulfur atoms provided that they are uncharged or negatively charged and that their surface is accessible.

Metal interactions are calculated between metal atoms embedded in the receptor and metal acceptor atoms, which are identical to the hydrogen bond acceptor atoms. Their geometry is based on the calculated coordination geometry of the metal.15 Each coordination point that is not occupied by a receptor atom is checked for close ligand atoms, and the maximal distance deviation is set to 0.8 Å. If no geometry can be calculated, a sphere with a radius of 2 Å is placed around the metal. In this case, all atoms lying on the sphere, again with a tolerance of 0.8 Å, are regarded as potential interaction partners.

In contrast to the formerly mentioned interactions, hydrophobic contacts are estimated based on the distance between two hydrophobic atoms only. They are visualized not by a dashed interaction line but by drawing the label of the contacting residue and a spline segment denoting the hydrophobic part of the ligand. Because many atoms typically form a hydrophobic subpocket, this representation reflects the interaction geometry better. A prerequisite for a hydrophobic contact is that at least three hydrophobic ligand atoms lie in the range of the currently examined receptor residue. Hydrophobic atoms are in this context carbon atoms with accessible surface and halogens. The maximum distance is set to the sum of the van der Waals radii of atoms in question and a tolerance of 0.8 Å.

EDITED :
Adding more  on a side note I came accross this paper, I am completely naive about the topic but the title :

“Observing Noncovalent Interactions in Experimental Electron Density for Macromolecular Systems: A Novel Perspective for Protein-Ligand Interaction Research.” Journal of chemical information and modeling vol. 62,7 (2022): 1734-1743. doi:10.1021/acs.jcim.1c01406

is kind of intriguing, I am not an expert and cannot vauch for its strictness but to me sounds like a good why to validate an interaction like seeing it even if in abstract terms:

... Non-covalent interactions (NCI) can be observed by pinpointing the ED saddle point, i.e. (3,−1) critical points and further quantified by measuring ED deviation from a homogeneous electron distribution using density and its first derivative (s = [1/(2(3π2)1/3)]|∇ρ|/ρ4/3)....

Please feel free to add to this, its not my field.
https://www.biorxiv.org/content/10.1101/2022.01.24.468575v2.full

no prob. I would be interested in the answer though, I mean just looking at interatomic distances and angles which are the most trusted definition of interactions, and more interesting how the cystallization-complex resolution influence the detection by alghoritm of such interactions (i.e.  how LigPlot+ or other software takes into account the resolution of the protein-ligand crystal under investigation).

And still I have problem to figure out how/if crystal resolution relate to the common concept of optical resolution

thank you for your elaborated comments


text:
Merge data frame on either of the 2 columns in R or python
I have two data frames. One data frame with one column and second data frame with two columns. I need to merge first data frame with either column of the second data frame and returns the values. After merging the values,the order of the first column should be same as the input. Please find the sample input below
first_column = c("ENSG00000165588","ENSG00000213551")
df1 = data.frame(first_column)
first_column = c("ENSG00000142192", "ENSG00000140575", "ENSG00000165588", "ENSG00000165588", "ENSG00000213551", "ENSG00000213551","ENSG00000197153")
second_column = c("ENSG00000165588", "ENSG00000165588", "ENSG00000186908", "ENSG00000135446", "ENSG00000273983", "ENSG00000274267","ENSG00000213551")
df2 = data.frame(first_column,second_column)

Please find the output below.
first_column = c("ENSG00000165588","ENSG00000165588","ENSG00000165588","ENSG00000165588","ENSG00000213551","ENSG00000213551","ENSG00000213551")
second_column = c("ENSG00000142192","ENSG00000140575","ENSG00000186908","ENSG00000135446","ENSG00000273983","ENSG00000274267","ENSG00000197153")
output = data.frame(first_column,second_column)

Just to confirm are you happy for a Python pandas solution? I recognise this is the implication of the title, but the code is in R

Yes, I am happy with python

This does not look like something I'd call a data frame merge. Are you able to give a bit more context around this problem (i.e. what you want to do with this data, rather than the intermediate problem you're trying to solve)? This helps us to wrap our heads around the problem and give better answers.

History This question was shifted from 'deleted' to 'closed' by me, which apparently was a mistake - apologies therein and then moved from 'closed' to 'opened'. I've now fully addressed the question. There was a misunderstanding: I originally provided a column merge solution when the OP wanted an index merge. Both solutions have now been provided. All code is vectorised and will be fast over a large data set.

I think this question is more about subsetting than merging. You can try this solution in R:
first_column = c("ENSG00000165588","ENSG00000213551")
df1 = data.frame(first_column)
first_column = c("ENSG00000142192", "ENSG00000140575", "ENSG00000165588", "ENSG00000165588",
                 "ENSG00000213551", "ENSG00000213551","ENSG00000197153")
second_column = c("ENSG00000165588", "ENSG00000165588", "ENSG00000186908", "ENSG00000135446",
                  "ENSG00000273983", "ENSG00000274267","ENSG00000213551")
df2 = data.frame(first_column,second_column)

output = rbind(df2[df2$first_column %in% df1$first_column,],
               df2[df2$second_column %in% df1$first_column,])

output <- transform(output, second_column = ifelse(output$second_column%in%df1$first_column,
                                                   first_column, second_column),
                            first_column = ifelse(output$second_column%in%df1$first_column,
                                                  second_column, first_column))
output <- output[ order(row.names(output)), ]

Edit The OP in hindsight is requesting an index merge. The first solution, presented here as Column merge solution was rejected by the OP.
From the comments below this post the order requested by the OP was:
# 0 ENSG00000165588 ENSG00000142192
# 1 ENSG00000165588 ENSG00000140575
# 2 ENSG00000165588 ENSG00000186908
# 3 ENSG00000165588 ENSG00000135446
# 4 ENSG00000213551 ENSG00000273983
# 5 ENSG00000213551 ENSG00000274267
# 6 ENSG00000213551 ENSG00000197153,

Column merge solution

The first operation is a concatenation and inner merge *.
The second operation are two outer merges, the second requiring a 'left-hand' merge *.

The code is here:
import pandas as pd

df1 = pd.DataFrame({'A':["ENSG00000165588","ENSG00000213551"]})
df2i = pd.DataFrame({'A':["ENSG00000142192", "ENSG00000140575", "ENSG00000165588",\
                          "ENSG00000165588", "ENSG00000213551", "ENSG00000213551", "ENSG00000197153"]})
df2ii = pd.DataFrame({'A':["ENSG00000165588", "ENSG00000165588", "ENSG00000186908", \
                            "ENSG00000135446", "ENSG00000273983", "ENSG00000274267","ENSG00000213551"]})

# Column 'A'
df3in = pd.concat([df2i, df2ii])
df4inner = pd.merge(df3in,df1, on='A', how='inner')

# Column 'A+'
df3out = pd.merge(df2i, df2ii, on='A', how='outer')
df4outer = pd.merge(df3in,df1, on='A', how='outer', indicator = True).query('_merge=="left_only"')
df5 = df4outer.drop('_merge', axis=1).reset_index(drop=True).rename(columns={'A':'A+'})
df6 = df4inner.join(df5)

print (df6)

Output
                 A               A+
0  ENSG00000165588  ENSG00000142192
1  ENSG00000165588  ENSG00000140575
2  ENSG00000165588  ENSG00000197153
3  ENSG00000165588  ENSG00000186908
4  ENSG00000213551  ENSG00000135446
5  ENSG00000213551  ENSG00000273983
6  ENSG00000213551  ENSG00000274267

This approach was rejected (comments) because of the position of ENSG00000197153 at row 2, Column A+ in relation to ENSG00000165588.

Index merge solution
import pandas as pd

def indexmerge(dfa, dfb, df1):
    df3out = dfa.reset_index().merge(dfb, how="left").set_index('index')
    df4outer = pd.merge(df3out,df1, on='A', how='outer', indicator = True).query('_merge=="left_only"')
    return df4outer

if __name__ == '__main__':

    df1 = pd.DataFrame({'A':["ENSG00000165588","ENSG00000213551"]})
    df2i = pd.DataFrame({'A':["ENSG00000142192", "ENSG00000140575", "ENSG00000165588",\
                              "ENSG00000165588", "ENSG00000213551", "ENSG00000213551", "ENSG00000197153"]})
    df2ii = pd.DataFrame({'A':["ENSG00000165588", "ENSG00000165588", "ENSG00000186908", \
                                "ENSG00000135446", "ENSG00000273983", "ENSG00000274267","ENSG00000213551"]})
    df3in = pd.concat([df2i, df2ii])
    df4inner = pd.merge(df3in,df1, on='A', how='inner')

    df4outer = indexmerge(df2i,df2ii,df1)
    df4outer2 = indexmerge(df2ii,df2i,df1)
    df4outer2.rename(columns={'A':'_A'},inplace=True)
    df6 = pd.concat([df4outer,df4outer2],axis=1)
    df6['A+'] = df6['A'].astype(str).replace('nan','') + df6['_A'].astype(str).replace('nan','')
    df7 = df6.sort_index().reset_index()
    df8 = df4inner.join(df7['A+'])
    print (df8)

Output
                 A               A+
0  ENSG00000165588  ENSG00000142192
1  ENSG00000165588  ENSG00000140575
2  ENSG00000165588  ENSG00000186908
3  ENSG00000165588  ENSG00000135446
4  ENSG00000213551  ENSG00000273983
5  ENSG00000213551  ENSG00000274267
6  ENSG00000213551  ENSG00000197153

This solution is identical to the order requested in the comments and is now also presented at the beginning of this post. The solution used here pivots around line 1 of indexmerge() and thus is an index merge operation. The remaining code is common pandas manipulation supporting this operation.
Note the new position of ENSG00000197153 row 6 column A+. The index merge retains ENSG00000197153 at the base of the column and therefore in association, i.e. same row, as ENSG00000213551.

* technical pandas terminology

Can you please check the output. The output must be in the dataframe. The order may be different but column A and column B entries should match. "ENSG00000213551" match in one case is wrong. If possible check my sample output above.

0 ENSG00000165588 ENSG00000142192
2 ENSG00000165588 ENSG00000140575
3 ENSG00000165588 ENSG00000186908
4 ENSG00000165588 ENSG00000135446
5 ENSG00000213551 ENSG00000273983
6 ENSG00000213551 ENSG00000274267
7 ENSG00000213551 ENSG00000197153, this is the output which i need.  For eg. ENSG00000213551 should interact with ENSG00000273983, ENSG00000274267, ENSG00000197153, but your program output shows ENSG00000213551 interact with ENSG00000135446, ENSG00000274267,  ENSG00000197153.

I have added the sample output and my question is correct.


text:
What should GC coverage bias plot of exome data look like?
I want to see if there is any GC coverage bias in some paired-end Illumina human exome data. Update This is for a variant calling project. I used the CollectMultipleMetrics from GATK4 toolkit with -PROGRAM CollectGcBiasMetrics. An example of the command I used is:
gatk CollectMultipleMetrics -I input.bam -O prefix_dir_name -R input.ref_genome -PROGRAM CollectGcBiasMetrics -PROGRAM CollectAlignmentSummaryMetrics -PROGRAM CollectInsertSizeMetrics -PROGRAM QualityScoreDistribution -PROGRAM MeanQualityByCycle -PROGRAM CollectBaseDistributionByCycle

Update This is my result for the QC coverage bias with one sample highlighted:

There seems to be a large bias. I would expect the peak to be around 50% for human. And those weird peaks around 90% are also puzzling to me. On the other hand, I am working with exome data so it might be logical to have a bias towards the higher values of GC% because you have reads on GC-rich exons?
Update I have another plot with the bins produced by the above mentioned tool (.gc_bias.pdf)

Update first part of useful file (gc_bias.detail_metrics) which shows the reads per bin, as suggested by @Maximilian Press. This I will plot.
## htsjdk.samtools.metrics.StringHeader
# CollectMultipleMetrics --INPUT WESP_run/results/aligned/RM8398_MF_S8_MD.sorted.bam --OUTPUT WESP_run/results/aligned/post_alignment_QC/multiple_metrics/RM8398_MF_S8 --PROGRAM CollectAlignmentSummaryMetrics --P
ROGRAM CollectBaseDistributionByCycle --PROGRAM CollectInsertSizeMetrics --PROGRAM MeanQualityByCycle --PROGRAM QualityScoreDistribution --PROGRAM CollectGcBiasMetrics --REFERENCE_SEQUENCE /mnt/scratch_dir/deute
koe/projects/HP/WESPipe/workflow/resources/reference_genomes/hg38/alt_plus/GCA_000001405.15_GRCh38_full_plus_hs38d1_analysis_set.fna --ASSUME_SORTED true --STOP_AFTER 0 --METRIC_ACCUMULATION_LEVEL ALL_READS --IN
CLUDE_UNPAIRED false --VERBOSITY INFO --QUIET false --VALIDATION_STRINGENCY STRICT --COMPRESSION_LEVEL 2 --MAX_RECORDS_IN_RAM 500000 --CREATE_INDEX false --CREATE_MD5_FILE false --GA4GH_CLIENT_SECRETS client_sec
rets.json --help false --version false --showHidden false --USE_JDK_DEFLATER false --USE_JDK_INFLATER false
## htsjdk.samtools.metrics.StringHeader
# Started on: Sat Dec 03 14:25:15 CET 2022

## METRICS CLASS        picard.analysis.GcBiasDetailMetrics
ACCUMULATION_LEVEL      READS_USED      GC      WINDOWS READ_STARTS     MEAN_BASE_QUALITY       NORMALIZED_COVERAGE     ERROR_BAR_WIDTH SAMPLE  LIBRARY READ_GROUP
All Reads       ALL     0       139942  1326    19      0.413544        0.011357
All Reads       ALL     1       102412  61      17      0.025996        0.003328
All Reads       ALL     2       118957  74      16      0.02715 0.003156
All Reads       ALL     3       152181  929     20      0.266429        0.008741
All Reads       ALL     4       167732  248     16      0.06453 0.004098
All Reads       ALL     5       173319  232     16      0.058421        0.003836
All Reads       ALL     6       196365  304     15      0.067567        0.003875
All Reads       ALL     7       211217  343     17      0.070875        0.003827
All Reads       ALL     8       241125  441     18      0.079822        0.003801
All Reads       ALL     9       270272  1068    18      0.172463        0.005277

I have read a paper that uses the CollectGcBiasMetrics on exome data and they magically produce a plot without bias. I say magically because they don't explain how they use CollectGcBiasMetrics.
Update I added the target coverage plot.

Update Am I missing something? Should I, for instance, use my target regions in instead of the reference genome when using CollectGcBiasMetrics? I will try it, but I am not sure this is the correct thing to do... It is hard to find any workflow for exome data with gatk and particularly CollectGcBiasMetrics, which specifies to use the reference genome.
From the comment @Maximilian Press

Do you have on-instrument diagnostics? Did you run FastQC on the raw reads? What do the other picard metrics look like? Did you manually inspect your alignments (strongly recommended)? Do you actually have coverage of the exome?

I will see if I can get some instrument diagnostics. I did look at the alignments, and as far as I can tell I don't see anything weird? Just that it is very deeply sequenced. I do see some intronic regions, with some reads, but only a few reads compared to the many on the targets, which I think is a consequence of the (too) deep sequencing. The fastqc of the raw reads didn't look troubling, and it did not give me any warnings. The picard metrics I thought did not look out of the ordinary. There is good target coverage. Yes, it is sequenced quite deep. Also, the lab tech did tell me they over-cycled because they wanted to be sure to have all the regions/targets.
Update So this bias might be due to too deep sequencing? And another important question, can I still reliably use this data? If so, I guess I can perform some corrections like suggested here
Any help/suggestions would be very welcome. Thank you in advance.

A bit of upvoting would be cool, its very good QC indeed and clearly presented. @supertech's answer makes sense to me. The GC% bias in living organisms could only represent bacteria IMO, I've seen this in soil bacteria so this has to be a lab artefact.

Thank you @M__. About the bacterial IMO (which I think means a contamination?). There is a separate run (another day a protocol slightly different) that shows the same bias. So either the lab tech is very sloppy, or the protocol was not correct (too many PCR cycles or something).

Thank you for incorporating the comments (and your answers to the comments) into your question. This is how comments should be used, and makes my job as a moderator much easier.

Am I missing something? Should I, for instance, use my target regions
instead of the reference genome?

Yes you should do this. Why do you think it's wrong doing that? It should be the first thing to check. And that GC bias might be effected by multiple factors. Most prominent one would be PCR over amplification. Most PCR polymerases give horrible GC bias when library is over cycled. There are many articles published on this.
I also suggest that you look at the number of reads in each bin. Peaks around 90% GC should not puzzle you so much. There is chance that you have only few reads in these bins out of millions of reads.

Just to formalise my comments.
@Supertech's is the right answer.
Background The one thing that is unusual is that the bias peaks at around 75%. The only time I've seen strong GC% bias is for particular groups of bacteria. 75% is about the maximum bias that can be achieved in a protein coding gene and these "extreme bacteria" hit that limit.
I better try and explain that a bit better. The primary selection pressure on a gene is the amino acid composition and the GC% is a second issue, i.e. selection pressure, for some organisms (e.g. certain bacteria, insects ..). Thus amino acids are not mutating, but the triplet codon degeneracy permits a bias in the GC% without alternating a given amino acids. In other words this GC% bias is occurring only for "neutral" mutations (clearly there must be a selection pressure, so they are not really neutral).
Summary and personal thoughts In summary the GC% limit is around 75% for a protein coding gene and beyond that it requires the amino acids to be changed. The peak in @Dandelion's data on that limit is weird because GC% artefacts from Taq errors are not subject to selection pressure: there should not be a peak.
"Contamination" I would avoid as an explanation because in my experience these groups of bacteria are unlikely to invade a lab environment. Sample mix up is a possible explanation. However, I would be cautious about assigning blame because in the chain of command within the wet-lab it could have occurred anywhere.
Official response To reiterate the formal and official answer would be @supertech's response and the 75% GC peak remains a mystery.

Please note The second plot has a distribution around GC content around 40% ... There appears a discrepancy between plots. Neither output appears compatible with the GC% of humans, but 40% GC with a reasonable proportion of that data on the 50% mark is better the 75% peak plot. Personally I would investigate the empirical basis for the differences.


text:
Manually calculating log2 fold change values from DESeq2 normalized counts
This question was also asked on Biostars and Bioconductor Support
My aim is to perform gene clustering based on RNA seq datasets (raw expression values) from NCBI database. After clustering I realized that it is hard for me to get the biological understanding of the clusters. So, I thought it will be better if I start with data such as expression change in different conditions with respect to control (log fold change) for my later interpretation.
I need to calculate log2 fold change values for lot of different experimental conditions when compared to their corresponding controls. Just to mention, I am not going to use these for differential expression analysis but for some other downstream analysis like clustering and stuff. Traditionally in my field, counts are normalized by TPM method and then fold change values are calculated by log2(TPM_exp+1)-log2(TPM_control+1) [using 1 or 0.5 as pseudo counts for log transformation]. In my case, I realized that TPM is not a good way to normalize the data as I have few samples with lot of reads mapping to only one or two genes [RNA composition bias]. DESeq2 median of ratios nomalization seems to take care of that issue. So, I prefer using DESeq2 normalization. But I cannot use DESeq2 for getting log2 fold change values because I don't have replicates for some of the experimental conditions and DESeq2 needs replicates to estimate log2 fold change values. So, I want to manually calculate log2 fold change values from DESeq2 normalized counts.
So, I am using log2(DESeq2norm_exp+0.5)-log2(DESeq2norm_control+0.5) for calculating log2 fold change values. I am not sure whether it is a good idea or the choice of pseudo-count here is very critical. Any comments or help is really appreciated.
[This question was originally posted on StackOverflow, and someone commented suggesting other communities that would be more appropriate, in particular, here and Bioconductor Support]

Can you please provide more information/detail about your experiment? In particular, how many experimental conditions there are, whether there are any interactions (e.g. condition x genetic variant), and the number of replicates per condition. Depending on the experimental design, it might be possible to use other conditions as replicates to help with normalisation and/or L2FC shrinkage.

A semantic comment: if you're calculating log2 fold change (i.e. the log of the fold change between two conditions), you're doing a differential expression analysis. It's possible that what you actually want are log-normalised expression values, but there's not enough detail to determine that.

Thanks, that's a little bit more information, but what's most important in this case is your experimental design. DESeq2 doesn't need replicates for every condition, but it needs enough biological variation to gat a reasonable idea of variability.


text:
Interpreting short indel calls in 1000 Genomes Project VCFs
Consider the following short indel polymorphism rs59679400 on chr7.
$ rsidx search dbSNP/dbSNP_GRCh38.{vcf.gz,rsidx} rs59679400 | cut -f 1-7
[rsidx] running version 0.3
NC_000007.14    12510038        rs59679400      GAAAA   G,GA,GAA,GAAA,GAAAAA    .       .

In the Phase III (2015) 1000 Genomes VCFs (GRCh37), I see the following record for the variant.
$ rsidx search 1000GenomesPhaseIII/ALL.chr7.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.{vcf.gz,rsidx} rs59679400 | cut -f 1-7
[rsidx] running version 0.3
7       12549664        rs59679400      GA      G       100     PASS

That's easy enough to interpret, although it looks like only two of the several possible alleles at that locus are represented.
When I look at the VCFs from the latest NYGC update to the 1000 Genomes data, I see the following. (I'll note that it's really inconvenient that RS numbers are no longer included for any variants.)
$ tabix 1000GenomesNYGC/1kGP_high_coverage_Illumina.chr7.filtered.SNV_INDEL_SV_phased_panel.vcf.gz chr7:12510038-12510038 | cut -f 1-7
chr7    12474780        HGSV_106040     G       <DUP>   .       .
chr7    12510038        7:12510038:GA:G GA      G       .       .
chr7    12510038        7:12510039:GAA:G        GAA     G       .       .

Ignoring the first record, which is a larger structural variant, we see two records at the position in question. My initial question is why they've used two records to encode variation for a single indel variant? As shown in the dbSNP record in the first snippet, VCF allows one to encode an arbitrary number of alternate alleles in a single record.
My second—and more important—question is how one is to interpret these records and resolve the sequence for any particular individual genotyped at the locus? Selecting an arbitrary sample HG00103, I get the following allele calls for the two records? What am I to make of this?
>>> from pysam import VariantFile
>>> vcf = VariantFile("1000GenomesNYGC/1kGP_high_coverage_Illumina.chr7.filtered.SNV_INDEL_SV_phased_panel.vcf.gz")
>>> records = list(vcf.fetch("chr7", 12510037, 12510038))
>>> len(records)
3
>>> records[1].samples["HG00103"].alleles
('G', 'G')
>>> records[2].samples["HG00103"].alleles
('GAA', 'GAA')

this is Marta from the NYGC. Thanks for reaching out with your questions.
The full SNV/INDEL callset that we released (listed as "SNV/INDEL VCFs" in the key resource table in our paper) includes multiallelic representation of all multiallelic sites (i.e. multiple ALT alleles are listed together in column 5 separated by commas) and the record you showed as an example above is listed there in the following way:
chr7    12510038    .   GAAA    GAA,GA,G
The VCF file you were looking at is the phased panel that we generated and released alongside the full unphased callset. The phasing software tools that we used in the paper (SHAPEIT2, SHAPEIT4, Eagle2) as well as most other phasing tools do not support multiallelic sites and require that such sites are split into separate lines and represented as biallelic for the purposes of phasing. As you noted, this is not ideal and can lead to difficulties in interpretation (as others explained above) but it's the best we can do with the existing approaches. We mention this in the Methods section of the paper as well as in the README that's included with the phased panel on the IGSR ftp (all linked in the paper). Please let us know if you have any further questions related to the high-coverage 1kGP resources. Thanks!

Thank you Marta!

Older GATK was using the so-called variant model where overlapping variants were combined during calling. More recent GATK is using the allelic model where one non-reference allele is evaluated at a time.
With the allelic model, a reference allele doesn't mean the sample possesses a reference allele; it means the sample doesn't have the alternate allele. For this particular example, records[1] means a homozygous -A deletion. records[2] means the sample doesn't have the -AA deletion. The two records are consistent with each other.
It is still possible to have conflictive records (for example, homozygous -A deletion and homozygous -AA deletion at the same location). At least this used to be a problem with GATK. Not sure if/how much this has been improved recently.

Often multi-allelic sites would be split into several lines (one per alternate allele) to make filtering easier for specific types of alternate alleles (e.g. filter out MNV vs SNV).
VCF-level, one thing you can do to identify the actual genotype over the locus easier is to join the split alleles into a single multi-allelic record. You can easily do this using bcftools:
bcftools view --samples HG00103 --regions chr7:12510038-12510038 1kGP_high_coverage_Illumina.chr7.filtered.SNV_INDEL_SV_phased_panel.vcf.gz | bcftools norm -m + | bcftools view -H

From here is clear that starting on position chr7:12510038, instead of "GAA", you have homozygous "GA" (1|1), therefore a hmozygous "A" deletion, which could be either the first or second "A" at positions 12510039 and 12510040.


text:
How to preserve colour assignments in matplotlib when a category is empty?
Some of my samples have empty clusters. I want to keep the same color for each cluster in different plots.
How do I get sc.pl.umap to keep the same colors for each cluster, even when a cluster is empty?
from matplotlib.pyplot import rc_context
with rc_context({'figure.figsize': (3, 3)}):
    fig, (ax1, ax2) = plt.subplots(2, 2, figsize=(6,6))
    ax1[0] = sc.pl.umap(adata[adata.obs['sample'] == adata.obs['sample'].cat.categories[0]],
                 color = 'leiden_r1', size=20,
                 legend_loc = 'on data', ax=ax1[0],
                 palette = mycolormap_26,
                 frameon=True, legend_fontsize = 8,title = adata.obs['sample'].cat.categories[0], show=False)
    ax1[0].set_xlabel('')
    ax1[0].set(xlim=xlim, ylim=ylim)
    ax1[1] = sc.pl.umap(adata[adata.obs['sample'] == adata.obs['sample'].cat.categories[1]],
                 color = 'leiden_r1', size=20,
                 legend_loc = 'on data', ax=ax1[1],
                 palette = mycolormap_26,
                 frameon=True, legend_fontsize = 8,title = adata.obs['sample'].cat.categories[1], show=False)
    ax1[1].set_xlabel('')
    ax1[1].set_ylabel('')
    ax1[1].set(xlim=xlim, ylim=ylim)
    ax2[0] = sc.pl.umap(adata[adata.obs['sample'] == adata.obs['sample'].cat.categories[2]],
                 color = 'leiden_r1', size=20,
                 legend_loc = 'on data', ax=ax2[0],
                 palette = mycolormap_26,
                 frameon=True, legend_fontsize = 8,title = adata.obs['sample'].cat.categories[2], show=False)
    ax2[0].set(xlim=xlim, ylim=ylim)
    #ax2[0].set_xlabel('')
    ax2[1] = sc.pl.umap(adata[adata.obs['sample'] == adata.obs['sample'].cat.categories[3]],
                 color = 'leiden_r1', size=20,
                 legend_loc = 'on data', ax=ax2[1],
                 palette = mycolormap_26,
                 frameon=True, legend_fontsize = 8,title = adata.obs['sample'].cat.categories[3], show=False)
    ax2[1].set(xlim=xlim, ylim=ylim)
    #ax2[1].set_xlabel('')
    ax2[1].set_ylabel('')

Thanks

Good question. I assume there's a import scanpy as sc . What about inserting a dumbie value, say single 0,0 to prevent the skip? It appears UMAP1 cluster 1 is absent causing 2 and 3 ... thereon to change colour (colour) schemes for the proceeding colour (colour)

I can customize the color palette mycolormap_26 for each specific sample, which is tedious.

This could be scripted, which would remove the tedium, so it would only need to be done once. What you've hit is a shortcoming in ScanPy and could be reported (answer below).

My recommendation is to report the issue directly on the ScanPy Github page here: https://github.com/scverse/scanpy/issues
This is a shortcoming with the testing of the algorithm, without any question. Replicates with experimental modification are standard procedure in this line of analysis. For example, one in vitro cell culture does not receive the drug therapy, i.e. a control, and the rest do, which results in one group with a complete absence of a transcript profile. The authors simply have not described this possibility in any of their documentation and therefore may not have originally considered it.
What I suspect has happened is the authorship only tested extensive and complete data - like really big data sets. "BigData" is what ScanPy is targeting because they look to have adapted to large data in Wolf et al Genome Biology 2018 paper, notably by specifying how quick their algorithm performs with increasing data set size.
What I suspect has happened is that smaller data sets, with experimental replication and modification, were not tested so these issues simply get missed. It's easy to overlook, the authors can't think of everything.
The authors of ScanPy are responding to issues arising unevenly. So it will be hit and miss whether they respond. Similar issues should have been reported to them before because its a common wet-lab experimental design. However, I've not checked their full intray.
They should provide the bit of code that will prevent the tedium described above, but then again they may not. If they do, please do write back. The whole ScanPy project is interesting.


text:
Is there a way I can combine geological data with population genetics data?
So I have ran SAMOVA and made calculations regarding haplotype diversity and nucleotide diversity of several populations. I have also made raster data regarding geographical data like slope, ruggedness etc. Is there a program/package that allows me to combine these to types of data into a map, so I can study if there is a correlation between these indexes and the diversity?

Update I looked into modern use of Mantel's test and it will not accept categorical variables. I now remember the null hypothesis of Mantel's ... the association is random between the matrices, i.e. no genetic structure. Therefore p < 0.05 excluding degrees of freedom (so in practice its smaller).
The options are:

Fst implementing that is easy but interpreting is complicated because it's about what constitutes genetic structure

this can be done in Arlequin https://www.youtube.com/watch?v=U2fj3pfPUlE as a tutorial
Fstat is the more famous package, but notably old

General linear modelling;
Machine learning for example random forest, SVM etc ...

I would use points 1 and point 3, 2 is old school now. I would use point 1 to datamine the ML.
The key to Fst is that 0.2 comprises population structure for complicated reasons. Implementing ML is not trivial but is the trendiest analysis going. It is not simply implementation its checking and verifying the associations and avoiding overtraining, which would be the major criticism.
The ML would permits the genetic structure to be predicted across broad geographic regions and those predictions are tested and used to reformulate the model. However, you'd need significant expertise in ML (and coding).

The technical word for this is allopatry and if that doesn't mean much, please do read on the subject because that is central to the theory behind the calculations you are performing. There is a HUGE bank of theory on the subject and is the central alternative to Darwinian selection, yep it's that big. So this is not a passing calculation but a fundamental calculation at the heart of the model of speciation of the organisms within your study. What you want to show is your species does not demonstrate allopatry but niche selection. Ultimately this may not be precisely decipherable (depends on the geographic distances).
Very briefly the analysis to perform now is the Mantel's test. This should be within the Arlequin package. Mantel's test will assess if there is a 'cline' in your data, i.e. an association with gene frequency and geographic distance - in an allele frequency package (genetic distance of point mutation otherwise).
Geological data versus geographic distance ... Mantel's test will cover both calculations, however I think that a simple numerical code will be needed for categorical data. I'd need to check that carefully, because I get that wrong the calculation with be wrong. Basically you need to tell the matrix correlation test that this categorical data being tested.
The next level of analysis in allele frequency calculations involves plotting the pairwise (not global) Fst against against geographic distance and geological type - if the association in Mantel's test is not linear. If the association is linear the Mantel's will be positive, i.e. null hypothesis will be accept, p > 0.05.
I'll check the location of the calculations in various packages later (not necessarily today). I need to check the categorical variable for Mantel's - I've never done that calculation so am not sure.
What you are looking for is Mantel's to reject geographic distance and accept geological data type. That would be a decisive and very publishable result. Fst will be the fall-back to shift between the hypotheses.
If you have loads of SNPs there are phylogenetic tree solutions based on modelling point mutation. The next layer of biophylo-geographic analysis involves packages within the Beast framework. I'm sure that would work here because the number of mutations would be low, also this is really just geographic distance rather than niche. If you've loads of SNPs this would be a good analysis to perform however, but for intraspecifies mtDNA particularly involving cox1 I personally wouldn't think it would be informative. It is perfectly doable though.

A positive Mantel's test is a rare result, but generally allopatry is very common so it is a VERY strong null hypothesis for the data. What I am trying to say is you can't just reject allopatry without testing it. ... well you can but not if you want to publish in e.g. Molecular Ecology.
Good luck with the study and again I'll update the post later. It looks really good and is more than interesting.

I forgot to mention in my original post that I did run a mantel's test, and in fact it's results are what prompted me to check geographical data in the first place. For two of my 4 species there was no correlation between genetic and geographical distance (p value>0.05) , so my thought process is that diversity comes from other geological factors.

I did the test in R using the packages ape vegan and geodist, due to not wanting to complicate things even more with arlequin. I am pretty sure the p value needs to be above 0.005 for it to be significant but I might be wrong.

@Nickmofoe I've answered the question, as far as I practically can.


text:
Searching motifs in sequence and their frequencies
This is a two part question.
I am searching for a motif, and in that search I wanted to also find the total number of sequences in my FASTA file, but the code I wrote is not yielding that please see the attempt below. The second part is when I have found the matching motifs, as well as the frequencies of the matches how do i go about to exclude the reference sequence matches and only get matches and frequncies from sequence2-sequence5.
> Ref
MPPRRSIVEVKVLDVQKRRVPNKHYVYIIRVTWSSGATEAIYRRYSKFFDLQMQMLDKFP
MEGGQKDPKQRIIPFLPGKILFRRSHIRDVAVKRLIPIDEYCKALIQLPPYISQCDEVLQ
FFETRPEDLNPPKEEHIGKKKSGNDPTSVDPMVLEQYVVVADYQKQESSEISLSVGQVVD
IIEKNESGWWFVSTAEEQGWVPATCLEGQDGVQDEFSLQPEEEEKYTVIYPYTARDQDEM
NLERGAVVEVVQKNLEGWWKIRYQGKEGWAPASYLKKNSGEPLPPKLGPSSPAHSGALDL
DGVSRHQNAMGREKELLNNQRDGRFEGRLVPDGDVKQRSPKMRQRPPPRRDMTIPRGLNL

>Sequence1
MAEVRKFTKRLSKPGTAAELRQSVSEAVRGSVVLEKAKLVEPLDYENVITQRKTQIYSDP
LRDLLMFPMEDISISVIGRQRRTVQSTVPEDAEKRAQSLFVKECIKTYSTDWHVVNYKYE
DFSGDFRMLPCKSLRPEKIPNHVFEIDEDCEKDEDSSSLCSQKGGVIKQGWLHKANVNST
ITVTMKVFKRRYFYLTQLPDGSYILNSYKDEKNSKESKGCIYLDACIDVVQCPKMRRHAF
ELKMLDKYSHYLAAETEQEMEEWLIMLKKIIQINTDSLVQEKKDTVEAIQEEETSSQGKA
ENIMASLERSMHPELMKYGRETEQLNKLSRGDGRQNLFSFDSEVQRLDFSGIEPDVKPFE
EKCNKRFMVNCHDLTFNILGHIGDNAKGPPTNVEPFFINLALFDVKNNCKISADFHVDLN
PPSVREMLWGTSTQLSNDGNAKGFSPESLIHGIAESQLCYIKQGIFSVTNPHPEIFLVVR

>Sequence2
GDDSEWLKLPVDQKCEHKLWKARLSGYEEALKIFQKIKDEKSPEWSKYLGLIKKFVTDS
NAVVQLKGLEAALVYVENAHVAGKTTGEVVSGVVSKVFNQPKAKAKELGIEICLMYVEIE
KGESVQEELLKGLDNKNPKIIVACIETLRKALSEFGSKIISLKPIIKVLPKLFESRDKAV
RDEAKLFAIEIYRWNRDAVKHTLQNINSVQLKELEEEWVKLPTGAPKPSRFLRSQQELEA
KLEQQQSAGGDAEGGGDDGDEVPQVDAYELLDAVEILSKLPKDFYDKIEAKKWQERKEAL
EAVEVLVKNPKLEAGDYADLVKALKKVVGKDTNVMLVALAAKCLTGLAVGLRKKFGQYAG
HVVPTILEKFKEKKPQVVQALQEAIDAIFLTTTLQNISEDVLAVMDNKNPTIKQQTSLFI
ARSFRHCTSSTLPKSLLKPFCAALLKHINDSAPEVRDAAFEALGTALKVVGEKSVNPFLA

> Sequence3
MPPRRSIVEVKVLDVQKRRVPNKHYVYIIRVTWSSGATEAIYRRYSKFFDLQMQMLDKFP
MEGGQKDPKQRIIPFLPGKILFRRSHIRDVAVKRLIPIDEYCKALIQLPPYISQCDEVLQ
FFETRPEDLNPPKEEHIGKKKSGNDPTSVDPMVLEQYVVVADYQKQESSEISLSVGQVVD
IIEKNESGWWFVSTAEEQGWVPATCLEGQDGVQDEFSLQPEEEEKYTVIYPYTARDQDEM
NLERGAVVEVVQKNLEGWWKIRYQGKEGWAPASYLKKNSGEPLPPKLGPSSPAHSGALDL
DGVSRHQNAMGREKELLNNQRDGRFEGRLVPDGDVKQRSPKMRQRPPPRRDMTIPRGLNL

>Sequence5
GDDSEWLKLPVDQKCEHKLWKARLSGYEEALKIFQKIKDEKSPEWSKYLGLIKKFVTDS
NAVVQLKGLEAALVYVENAHVAGKTTGEVVSGVVSKVFNQPKAKAKELGIEICLMYVEIE
KGESVQEELLKGLDNKNPKIIVACIETLRKALSEFGSKIISLKPIIKVLPKLFESRDKAV
RDEAKLFAIEIYRWNRDAVKHTLQNINSVQLKELEEEWVKLPTGAPKPSRFLRSQQELEA
KLEQQQSAGGDAEGGGDDGDEVPQVDAYELLDAVEILSKLPKDFYDKIEAKKWQERKEAL
EAVEVLVKNPKLEAGDYADLVKALKKVVGKDTNVMLVALAAKCLTGLAVGLRKKFGQYAG
HVVPTILEKFKEKKPQVVQALQEAIDAIFLTTTLQNISEDVLAVMDNKNPTIKQQTSLFI
ARSFRHCTSSTLPKSLLKPFCAALLKHINDSAPEVRDAAFEALGTALKVVGEKSVNPFLA

The code to search the peptides is below
import re
import collections
m_list =[]
from Bio import SeqIO
input_file = 'Sequences2.fasta'
num_sequences = 0
fasta_sequences = SeqIO.parse(open(input_file),'fasta')
for fasta in fasta_sequences:
    num_sequence=+1
    name, sequence = fasta.id, str(fasta.seq)
    matches=re.finditer(r"P[A-Z]{2}P..",sequence)
    for m in matches:
        match_id = f'{m.start()}_{m.end()}_{m.group()}'
        m_list.append(match_id)
        print(name, m.start(), m.end(), m.group())
print("total sequences",num_sequences)
print("frequency",collections.Counter(m_list))

It yields the following output below, which I am happy with but now it doesn't give me the total number of sequences:
Ref 73 79 PFLPGK
Ref 281 287 PLPPKL
Ref 288 294 PSSPAH
Sequence3 73 79 PFLPGK
Sequence3 281 287 PLPPKL
Sequence3 288 294 PSSPAH
total sequences 0
frequency Counter({'73_79_PFLPGK': 2, '281_287_PLPPKL': 2, '288_294_PSSPAH': 2})

The total sequences should be 5. I used the code below to find out. I do not know where I went wrong since it is not giving me that on the more compact code above.
sequence = open('Sequences2.fasta')
x = 0
for line in sequence:
    if line.startswith('>'):
        x += 1
sequence.close()
print(x)

The second part of my question is how do I modify my code such that, the matches and frequencies of the reference sequence(Ref) are excluded. The idea would be to get an output such as the one below. I do not even know whether that is possible (A disclaimer this second part is a question on a tutorial that I don't know how to do. I am still learning)
Sequence3 73 79 PFLPGK
Sequence3 281 287 PLPPKL
Sequence3 288 294 PSSPAH
total sequences 5  #this would be the total number of sequences in file including the refrence
frequency Counter({'73_79_PFLPGK': 1, '281_287_PLPPKL': 1, '288_294_PSSPAH': 1}) #frequency of matches, excluding the reference

The first part appears to be a simple bug if there is no copy/paste error in the code:
Original code,
num_sequence=0
fasta_sequences = SeqIO.parse(open(input_file),'fasta')
for fasta in fasta_sequences:
    num_sequence=+1 # < its here
print (num_sequence)

The code I'd recommend is
num_sequence=0
fasta_sequences = SeqIO.parse(input_file,'fasta')
for fasta in fasta_sequences:
    num_sequence += 1 # corrected
print (num_sequence)

I will look at the second part later.

for the second part re-arrange the for m in matches block to :
for m in matches:

        match_id = f'{m.start()}_{m.end()}_{m.group()}'

        print(name, m.start(), m.end(), m.group())

        if num_sequences > 1:       #### added if clause so that appends only after first sequence to m_list

                m_list.append(match_id)


text:
Display row and column names in R
I have plotted a heat-map with 235 rows and 570 columns which is below: how do we make clear names of row and columns on this graph using R?

Please show the command you used to create that heatmap. Answers will differ depending on what library / function you're using.

You don't, simple as that, it is too many in both columns and rows to show them all. You could shrink size, but then it becomes unreadable. You can highlight individual genes (rows) that are important to support the message you want to send with that experiment, with packages such as ComplexHeatmap, see the Mark annotation section.
Alternatively, you can use something like hierarchical clustering to form groups, then make functional annotations (such as REACTOME terms) and then annotate eaqch cluster with selected terms, again to support the message you want to send. Go through the ComplexHeatmap manual, lots of examples with code there to get inspiration from.
As for the columns, you could use block annotations to color-highlight columns bhelonging to the same experimental group, or split the heatmap by these groups to put some visual separation/guidance.

Thank you @ATpoint

[note: This is a generic answer; the post was not specific enough for me to be able to include any more useful detail]
A heatmap of 235 rows and 570 columns is within the realm of what should be explorable / visible with a bit of tweaking. I have needed to do this with large heatmap plots - for example, single cell data with 500 genes and 2000 cells, where I care about the genes but not the cell IDs. My usual approach is to make the dimensions of the output image larger, and also reduce the text size.
The method of doing that depends a lot on how the plot was generated initially. For example, if using a base R package:
pdf("output_file_name.pdf", width=57, height=23.5); # about 10 squares per inch
## <plot generation function>
dev.off()

If using ggplot:
library(tidyverse);
input_data %>%
  ggplot() +
  aes(...) +
  # <plot generation function>
ggsave("output_file_name.pdf", width=57, height=23.5);

Changing the function used to display the heatmap can sometimes help as well (e.g. pheatmap instead of heatmap).
Within the functions, text size can often be reduced to help with this as well (e.g. cex = 0.2). There is a limit to how much this can be done, as most PDF viewing programs have a zoom limit - I think Adobe Reader's limit is 6400% (i.e. 64 times magnification) - as well as maximum page dimension limits.

Thank you @gringer♦

This was useful, gringer.


text:
ERR: error while loading shared libraries: libjulia.so.1: cannot open shared object file [Julia] - Ubuntu 22.04
I'm trying to build a Dockerfile. It has a tool called Atria that utilizes Julia. Since Ubuntu 22.04 does not have Julia package, I had to resort to installing it. The following is the Dockerfile:
FROM ubuntu:22.04

# Install tzdata
RUN apt-get update &&\
    DEBIAN_FRONTEND=noninteractive TZ=Etc/UTC apt-get -y install tzdata

# Install Python, Bowtie2 and Java
RUN apt-get install -y python3.10 python3-pip \
    openjdk-8-jdk  \
    bowtie2 \
    wget

RUN  apt-get install --yes --no-install-recommends \
     zlib1g-dev \
     libbz2-dev \
     liblzma-dev

# Install RSeQC
RUN pip3 install RSeQC

# Install biopython=1.80
RUN pip3 install biopython

#Install Julia
RUN wget https://julialang-s3.julialang.org/bin/linux/x64/1.8/julia-1.8.1-linux-x86_64.tar.gz && \
    tar zxvf julia-1.8.1-linux-x86_64.tar.gz && \
    mv julia-1.8.1/bin/julia /usr/local/bin/julia && \
    chmod +x /usr/local/bin/julia

RUN /sbin/ldconfig -v

# Install Atria
RUN wget https://github.com/cihga39871/Atria/releases/download/v3.1.2/atria-3.1.2-linux.tar.gz && \
    tar -zxf atria-3.1.2-linux.tar.gz && \
    mv atria-3.1.2/bin/atria /usr/local/bin/atria && \
    chmod +x /usr/local/bin/atria

#Atria dependencies
RUN  apt-get install pigz pbzip2

# Install findtail
 RUN wget https://storage.googleapis.com/google-code-archive-downloads/v2/code.google.com/findtail/findtail_v1.01 && \
     mv findtail_v1.01 /usr/local/bin/findtail_v1.01 && \
     chmod +x /usr/local/bin/findtail_v1.01
# Cleanup
RUN apt clean

What I have tried so far is :

Julia v1.6.7 and v1.8.4

RUN /sbin/ldconfig -v as suggested here

How to fix it? Should I try to install Julia using PPA?
edit- The error:
root@b0f4b17a7c18:/# atria
/usr/local/bin/julia: error while loading shared libraries: libjulia.so.1: cannot open shared object file: No such file or directory

It is recognizing the atria command tho.

You just need to make sure you include the required library files when you install Atria. Since you're not building any software (i.e. using make), the usual way is to simply move the directory under /opt and add the executable(s) to your PATH using the ENV keyword. Note that Julia comes bundled with Atria, so all you need is:
FROM ubuntu:22.04

RUN apt-get update && apt-get install -y \
    bowtie2 \
    openjdk-8-jdk \
    pbzip2 \
    pigz \
    python3.10 \
    python3-pip \
    python3-biopython \
    wget \
    && apt-get clean \
    && rm -rf /var/lib/apt/lists/*

# Install RSeQC
RUN pip3 install --no-cache-dir RSeQC

# Install Atria
RUN wget -q https://github.com/cihga39871/Atria/releases/download/v3.1.2/atria-3.1.2-linux.tar.gz \
    && tar xf atria-3.1.2-linux.tar.gz \
    && mv atria-3.1.2 /opt/atria

# Install findtail
RUN wget -q https://storage.googleapis.com/google-code-archive-downloads/v2/code.google.com/findtail/findtail_v1.01 \
    && chmod +x findtail_v1.01 \
    && mv findtail_v1.01 /usr/local/bin/findtail

# Add the Atria executable to the PATH
ENV PATH="/opt/atria/bin:${PATH}"

Hello @Steve, thank you for that. I'm still a novice (obviously lol) trying to wrap my head around this stuff and I so appreciate your time and effort in helping me out.

I don't understand what modifications to make, as the question and answer are rather sizeable files, can you show a more minimal example?

Just to note RUN apt clean  ... is much better performed using && apt clean, like in @Steve's approach, recommended for complicated reasons.
Suggestion: what about using the official Julia Docker as the base rather than create it from scratch? There isn't a Ubuntu Julia but there is a Debian base.
docker run -it --rm Julia

or
docker pull julia

or via a Dockerfile,
FROM Julia:buster

or
FROM Julia:latest

Doesn't work when you want to layer Julia on top of an existing image that contains a lot of stuff you don't wont to re-layer on the Julia base image


text:
CONECT doesnt include all the connections the 1vfl macro molecule have in a pdb file
I have a 1vfl macro molecule with 3000 atoms in it. When i open it in DC accelrys visualizer. It has much bonds. but when i open the pdb file it just have 8 connections. You can download the pdb here(I'm linking because the pdb is much) : https://www.rcsb.org/structure/1VFL
I'm using this pdb for my program so i need to know all the connections between atoms.
Any help would be good; I'm a complete beginner at using PDB.

The CONECT entries specify the bonding of the ligands (HETATM), not your polymeric residues (ATOM). By residue I mean a collection of atoms that forms a unit, not only your standard amino acid definition.
You can export in PyMOL all connections if you wanted, but the connectivity of the polymeric residues is standardised and does not need to be specified as they ought to obey the definitions (the 'chemical component' dataset), for example https://www.rcsb.org/ligand/ILE.
If I were to have a ligand, LIG, it would be HETATM entered and with CONECT entries as most likely it is not the RCSB's https://www.rcsb.org/ligand/LIG.
Crosslinks between certain residues ought to be entered via the LINK record, not CONECT so they will not feature.
Some programs connect atoms by proximity regardless of stated connectivity —spiderwebs in PyMOL.
An odd corner case are non-standard amino acids, such as phosphorylations (e.g. https://www.rcsb.org/ligand/SEP) which are ATOM not HETATM, which are not meant to have CONECT entries, but some were deposited with them.
Water (HOH) (from neutron diffraction structures) does not need connections.


text:
Best method to compare differentially expressed genes between 2 single cell clusters for GSEA
I want to get a ranked gene list for GSEA analysis from single cell data, which needs a group differential expression statistic for ranking.
Which is the best method to compare the differently expressed genes between 2 clusters scanpy.tl.rank_genes_groups (https://scanpy.readthedocs.io/en/stable/generated/scanpy.tl.rank_genes_groups.html) or MAST.zlm (https://rdrr.io/bioc/MAST/man/zlm.html) or FindMarkers function of Seurat (https://satijalab.org/seurat/reference/findmarkers)?
Thanks

Just to add I personally think this is a cool question, others may disagree. There is quite a bit of discussion on DEGs per se ... ScanPy hits that R to Python threshold, which is no trivial issue IMO.

I wish to focus my contribution singly on ScanPy vs. Seurat FindMarkers.
ScanPy's claim is it is essentially a speeded up version of Seurat FindMarkers with better performance (discussed below) written in Python.
What they are doing are essentially datamining the expression signal using multivariate statistics (PCA) focused through tSNE. There is no controversy or dispute with this approach at all (please see Caveat I shouldn't use the term "datamining" in context - bad habit).
What ScanPy is doing is using Graph theory Louvian groups to extract the differential signal, presumably from the tSNE focused PCA. Cool.
If you've already got to grips with ScanPy then leveraging it as a data mining approach - to me - looks sensible. The statistics will not be hugely powerful, that doesn't mean they are bad, it approximately means they are solid.
What is ScanPy doing?
The speeding up is via parallelisation, which is very easy in Python. The implication is Seurat is not parallelised, which would amaze me. The authors talk about some form of 'chunking' (they say 'partially loading the object') via HDF5. This is officially what they say

ANNDATA is similar to R’s EXPRESSIONSET [26], but supports sparse data and allows HDF5-based backing of ANNDATA objects on disk, a format independent of platform, framework, and language. This allows operating on an ANNDATA object without fully loading it into memory

I am confused by this, because HDF5 is a rapid I/O format - its not a datachunking format which his what the authors imply. Yes it facilitates datachunking but simply via rapid I/O. HDF5 is of course independent Python, but why discuss that at length (we all know that) when the critical information is how the partially loaded object is achieved. There's a lot of strategies to achieve the last bit - none of which are discussed.
The overall levels of high performance Python  (outside ScanPy) - because the reality is Python's core code is slow - are three steps:

Parallelisation and code optimisation.
Data chunking - using algorithm strategies to partial load data e.g. Hadoop strategies.
Compiling via JIT, numba and static typing via  Cython.

The authors have done point 1, however again HDF5 doesn't facilitate point 2 above, all it does is enable faster read/write. HDF5 is not actually the fastest I/O either both Pacquet (sp?) and feather are faster but normally for pandas dataframes. HDF5 is actually just as fast as pickle under many tests (!).
So the paper is a bit weird in that it is sparsely described - but this does often happen admittedly. The partial loading stuff via HDF5 is weird, because they are not describing how that is done.
The other bit that is weird is the ML and ANNDATA class. They state,

ANNDATA is similar to R’s EXPRESSIONSET [26]

I rarely use R so to me that isn't clear. The other weird thing is the authors talk about ML, they use an ML library, but give no further details nor why it is used. What I would assume is either: that they use it to assist with pseudotemporal ordering analysis OR they use it for unsupervised learning (PCA, tSNE). The latter is more likely. Louvian groups are from the igraph package they use
So in my personal view, despite the paper being a bit (to very) obscure I see no reason why ScanPy cannot replace Seurat purely in terms of the datamining being performed. Python is just as good at PCA and tSNE etc ...  The use of Louvain groups is cool to trap that information. The key benefit of ScanPy is:

your already familiar with it.
faster - it definitely could be accelerated further.
good for Python

The drawbacks I can see are:

RNA analysis is dominated by R and you'll likely pick up a R focused reviewer (which could be fun).
As per your previous post there appear  shortcomings with ScanPy - but like the development team can't think of everything first time round.
ScanPy2 is not released.
its not generating a critical (test) probability to distinguish between expression levels just like Seurat - there are real advantages to this however.

Caveat  I keep saying "data-mining" and I shouldn't. If the reviewers are perfectly happy with the above analysis alone, that's cool. This is because what I repeatedly call "data-mining" are other investigators final results, particularly in RNA expression.  If your aim is to publish a big RNA-seq data set and identify the core features of e.g. a cellular immune response - good approach. I don't generate data, mine is purely analytics and algorithms, so my perspective (and reviewing standards) are different.

Just to state my interest in DEGs (or DGE). I use unsupervised learning, PCA, tSNE to generate a hypothesis and get a feel for the data. Then however, I do customised full-scale ML (supervised learning) to interrogate the data. Finally, I use DEG to assess which way around the ML result is pointing (under- or over-represented). I will not let the DEG override the ML (supervised learning) analysis - which I see as much stronger analysis. My context is different however and I'm consolidating the result is via extensive custom databases, i.e. very different.

In the ScanPy paper Wolf et al 2018, when the authors say "HDF5" I wonder if they mean "HDFS"? This would be the Hadoop Distributed File System and would actually make sense. I can't imagine it - it's a HUGE error - but the rationale with partial object loading would be completely understandable. It does demonstrate how sparsely this paper is written some algorithms in Genome Biology papers are so detailed it's incredible. This is the opposite spectrum.

Additional weird point The authors calling a class ANNdata for machine learning is also weird. ANN stands for Artificial Neural Network and is a basic type of deep learning. The authors definitely do not use deep learning in ScanPy ... their analysis appear to be singly unsupervised learning. It's a weirdly written paper. I guess the authors must be carrying kudos in algorithm development - that does happen. The reviewers can wave through a paper if the authorship is carries kudos and editor doesn't pick it up.