text:
How to add charges to a PDB file?
I did some calculations using a PDB file as input. As a result, I got a file with the atomic partial charges (only the values in one column, no other information) and the output PDB file structure without the charge information.
My question: Is there a tool to add the charges to the output PDB file using the separate file?
My PDB file:
    TITLE     sistema05
    REMARK   4      COMPLIES WITH FORMAT V. 3.0, 1-DEC-2006
    REMARK 888
    REMARK 888 WRITTEN BY MAESTRO (A PRODUCT OF SCHRODINGER, LLC)
    HETATM    1  C1  BAS A   1      10.151  -3.071  -1.881  1.00  0.00           C
    HETATM    2  C2  BAS A   1       9.389  -3.683  -0.874  1.00  0.00           C
    HETATM    3  C3  BAS A   1       9.976  -3.936   0.382  1.00  0.00           C
    HETATM    4  C4  BAS A   1       9.194  -3.742   1.530  1.00  0.00           C
    HETATM    5  C5  BAS A   1       9.781  -3.201   2.695  1.00  0.00           C
    HETATM    6  C6  BAS A   1       9.029  -2.357   3.510  1.00  0.00           C
    HETATM    7  C7  BAS A   1       9.651  -1.313   4.236  1.00  0.00           C
    HETATM    8  C8  BAS A   1       8.948  -0.139   4.488  1.00  0.00           C
    HETATM    9  C9  BAS A   1       9.635   1.100   4.502  1.00  0.00           C
    HETATM   10  C10 BAS A   1       9.001   2.259   4.094  1.00  0.00           C
    HETATM   11  C11 BAS A   1       9.724   3.177   3.329  1.00  0.00           C
    HETATM   12  C12 BAS A   1       9.148   3.815   2.228  1.00  0.00           C
    HETATM   13  C13 BAS A   1       9.909   3.945   1.081  1.00  0.00           C
    HETATM   14  C14 BAS A   1       9.332   3.754  -0.198  1.00  0.00           C
    HETATM   15  C15 BAS A   1      10.098   3.266  -1.245  1.00  0.00           C
    HETATM   16  C16 BAS A   1       9.501   2.445  -2.232  1.00  0.00           C
    HETATM   17  C17 BAS A   1      10.234   1.431  -2.847  1.00  0.00           C
    HETATM   18  C18 BAS A   1       9.575   0.238  -3.215  1.00  0.00           C
    HETATM   19  C19 BAS A   1      10.249  -0.986  -3.121  1.00  0.00           C
    HETATM   20  C20 BAS A   1       9.529  -2.124  -2.713  1.00  0.00           C
...
...
...
    CONECT    1  212    2   20
    CONECT    1   20
    CONECT    2    1    3   21
    CONECT    2   21
    CONECT    3    2  213    4
    CONECT    3    4
    CONECT    4    3    5   23
    CONECT    4    3
    CONECT    5    4  214    6
    CONECT    5    6
    CONECT    6    5    7   25
    CONECT    6    5
    CONECT    7    6  215    8
    CONECT    7    8
    CONECT    8    7    9   27
    CONECT    8    7
    CONECT    9    8  216   10
    CONECT   10    9   11   29
    CONECT   11   10  217   12
    CONECT   12   11   13   31
    CONECT   13   12  218   14
    CONECT   14   13   15   33
    CONECT   14   15
    CONECT   15   14  219   16
    CONECT   15   14
    CONECT   16   15   17   35
    CONECT   16   17
    CONECT   17   16  220   18
    CONECT   17   16
    CONECT   18   17   19   37
    CONECT   18   19
    CONECT   19   18  221   20
    CONECT   19   18
    CONECT   20    1   19   39
    CONECT   20    1
    CONECT   21    2   22   40
    END

My CHARGE file:
   -0.02474906
    0.02291960
   -0.02789302
    0.02227110
   -0.02831805
    0.02230034
   -0.02876149
    0.02225527
   -0.02774076
    0.02216971
   -0.02074528
    0.02275309
   -0.02662750
    0.02177531
   -0.03007329
    0.02283631
   -0.02779055
    0.02312780
   -0.02312680
    0.02315866

You could use the b factor or occupancy value column to store your charge values, dont remember how many digits are available for each of the 2 proprierties. Parse a pdb file with biopython and change it with your charge files, I assume the pdb doesnt have any header or remark just atoms from lines 1 onwards that relate to the list of charges one per line in the second file ?? A small snippet of both files would be of help for anybody interested in helping you out

See https://stackoverflow.com/questions/53938256/how-to-modify-columns-contents-in-a-protein-pdb-file-by-keys-and-values-stored-i

Sorry just realized pdb has column for that https://www.cgl.ucsf.edu/chimera/docs/UsersGuide/tutorials/pdbintro.html , column 79-80

Could use biopython https://biopython.org/docs/1.79/api/Bio.PDB.Atom.html and set the charge for each atom   atom.set_charge( CHARGE) should do

Thanks @pippo1980. I edited the post, adding part of the files. I will go through your suggestions.

Mmmh your charges are a lot of digits... . You'll need a custom made pdb thats not a pdb format anymore to store all that and a parser to be able to reload that pdb, are you sure there are not other more suitable formats

OK given , input_connect.pdb , using Biopython library
TITLE     sistema05
REMARK   4      COMPLIES WITH FORMAT V. 3.0, 1-DEC-2006
REMARK 888
REMARK 888 WRITTEN BY MAESTRO (A PRODUCT OF SCHRODINGER, LLC)
HETATM    1  C1  BAS A   1      10.151  -3.071  -1.881  1.00  0.00           C
HETATM    2  C2  BAS A   1       9.389  -3.683  -0.874  1.00  0.00           C
HETATM    3  C3  BAS A   1       9.976  -3.936   0.382  1.00  0.00           C
HETATM    4  C4  BAS A   1       9.194  -3.742   1.530  1.00  0.00           C
HETATM    5  C5  BAS A   1       9.781  -3.201   2.695  1.00  0.00           C
HETATM    6  C6  BAS A   1       9.029  -2.357   3.510  1.00  0.00           C
HETATM    7  C7  BAS A   1       9.651  -1.313   4.236  1.00  0.00           C
HETATM    8  C8  BAS A   1       8.948  -0.139   4.488  1.00  0.00           C
HETATM    9  C9  BAS A   1       9.635   1.100   4.502  1.00  0.00           C
HETATM   10  C10 BAS A   1       9.001   2.259   4.094  1.00  0.00           C
HETATM   11  C11 BAS A   1       9.724   3.177   3.329  1.00  0.00           C
HETATM   12  C12 BAS A   1       9.148   3.815   2.228  1.00  0.00           C
HETATM   13  C13 BAS A   1       9.909   3.945   1.081  1.00  0.00           C
HETATM   14  C14 BAS A   1       9.332   3.754  -0.198  1.00  0.00           C
HETATM   15  C15 BAS A   1      10.098   3.266  -1.245  1.00  0.00           C
HETATM   16  C16 BAS A   1       9.501   2.445  -2.232  1.00  0.00           C
HETATM   17  C17 BAS A   1      10.234   1.431  -2.847  1.00  0.00           C
HETATM   18  C18 BAS A   1       9.575   0.238  -3.215  1.00  0.00           C
HETATM   19  C19 BAS A   1      10.249  -0.986  -3.121  1.00  0.00           C
HETATM   20  C20 BAS A   1       9.529  -2.124  -2.713  1.00  0.00           C
CONECT    1  212    2   20
CONECT    1   20
CONECT    2    1    3   21
CONECT    2   21
CONECT    3    2  213    4
CONECT    3    4
CONECT    4    3    5   23
CONECT    4    3
CONECT    5    4  214    6
CONECT    5    6
CONECT    6    5    7   25
CONECT    6    5
CONECT    7    6  215    8
CONECT    7    8
CONECT    8    7    9   27
CONECT    8    7
CONECT    9    8  216   10
CONECT   10    9   11   29
CONECT   11   10  217   12
CONECT   12   11   13   31
CONECT   13   12  218   14
CONECT   14   13   15   33
CONECT   14   15
CONECT   15   14  219   16
CONECT   15   14
CONECT   16   15   17   35
CONECT   16   17
CONECT   17   16  220   18
CONECT   17   16
CONECT   18   17   19   37
CONECT   18   19
CONECT   19   18  221   20
CONECT   19   18
CONECT   20    1   19   39
CONECT   21    2   22   40

and charges list input_charge.txt :
-0.02474906
 0.02291960
-0.02789302
 0.02227110
-0.02831805
 0.02230034
-0.02876149
 0.02225527
-0.02774076
 0.02216971
-0.02074528
 0.02275309
-0.02662750
 0.02177531
-0.03007329
 0.02283631
-0.02779055
 0.02312780
-0.02312680
 0.02315866

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

from Bio.PDB import (
                    PDBParser,
                    Atom,

                    )

                    # PDBIO)

import PDBIOcopy

# from Bio.PDB import  PDBIO

pdb = 'input'

charges = 'input_charge.txt'

charges_list = []

with open(charges , 'r') as handle :

    for i in handle.readlines() :

        charges_list.append(i.strip())

print(charges_list)

parser = PDBParser(PERMISSIVE=True , QUIET = False)

structure = parser.get_structure(pdb, pdb+'.pdb')

atoms = structure.get_atoms()

idx = 0

for atom in atoms :

    print(atom)

    atom.set_charge(charges_list[idx])

    idx += 1

atoms = structure.get_atoms()

for atom in atoms :

    print(atom , atom.get_charge())

struct_to_save = PDBIOcopy.PDBIO()

# struct_to_save = PDBIO()

struct_to_save.set_structure(structure)

struct_to_save.save('result_input.pdb')

and modifying Bio.PDB.PDBIO as :
line 22 :
_ATOM_FORMAT_STRING = (
    "%s%5i %-4s%c%3s %c%4i%c   %8.3f%8.3f%8.3f%s%6.2f      %4s%2s%12s\n"
)

line 164 :
# charge="  ", #deleted
pqr_charge     # insertes

line 213:
charge = atom.pqr_charge  #inserted

line 389:
atom.pqr_charge,  #inserted

output result_input.pdb:
HETATM    1  C1  BAS A   1      10.151  -3.071  -1.881  1.00  0.00           C -0.02474906
HETATM    2  C2  BAS A   1       9.389  -3.683  -0.874  1.00  0.00           C  0.02291960
HETATM    3  C3  BAS A   1       9.976  -3.936   0.382  1.00  0.00           C -0.02789302
HETATM    4  C4  BAS A   1       9.194  -3.742   1.530  1.00  0.00           C  0.02227110
HETATM    5  C5  BAS A   1       9.781  -3.201   2.695  1.00  0.00           C -0.02831805
HETATM    6  C6  BAS A   1       9.029  -2.357   3.510  1.00  0.00           C  0.02230034
HETATM    7  C7  BAS A   1       9.651  -1.313   4.236  1.00  0.00           C -0.02876149
HETATM    8  C8  BAS A   1       8.948  -0.139   4.488  1.00  0.00           C  0.02225527
HETATM    9  C9  BAS A   1       9.635   1.100   4.502  1.00  0.00           C -0.02774076
HETATM   10  C10 BAS A   1       9.001   2.259   4.094  1.00  0.00           C  0.02216971
HETATM   11  C11 BAS A   1       9.724   3.177   3.329  1.00  0.00           C -0.02074528
HETATM   12  C12 BAS A   1       9.148   3.815   2.228  1.00  0.00           C  0.02275309
HETATM   13  C13 BAS A   1       9.909   3.945   1.081  1.00  0.00           C -0.02662750
HETATM   14  C14 BAS A   1       9.332   3.754  -0.198  1.00  0.00           C  0.02177531
HETATM   15  C15 BAS A   1      10.098   3.266  -1.245  1.00  0.00           C -0.03007329
HETATM   16  C16 BAS A   1       9.501   2.445  -2.232  1.00  0.00           C  0.02283631
HETATM   17  C17 BAS A   1      10.234   1.431  -2.847  1.00  0.00           C -0.02779055
HETATM   18  C18 BAS A   1       9.575   0.238  -3.215  1.00  0.00           C  0.02312780
HETATM   19  C19 BAS A   1      10.249  -0.986  -3.121  1.00  0.00           C -0.02312680
HETATM   20  C20 BAS A   1       9.529  -2.124  -2.713  1.00  0.00           C  0.02315866
TER      21      BAS A   1
END

That unfortunately is not a pdb file format anymore:
PDB Format uses for RECORD ATOM CHARGE columns 79-80 , see Protein Data Bank Contents Guide: Atomic Coordinate Entry Format Description
At the end of all of this I realized that using python and just adding the charge values at the end of the pdb file read line by line would be a lot easier & faster but anyway I had the chance to explore PDBIO
If you want to add the charges in the bfactor or occupancy columns just use :atom.set_bfactor(charges_list[idx]) or atom.set_occupancy(charges_list[idx]) be sure to modify PDBIO to make space for your values
There are other examples on SO about same problem too :
How to modify columns contents in a protein PDB file by keys and values stored in a text file using only bash tools
Just Python code :
pdb = 'input_connect.pdb'

charges = 'input_charge.txt'

charges_list = []

with open(charges , 'r') as handle :

    for i in handle.readlines() :

        charges_list.append(i.strip())

print(charges_list)

with open('input_connect.pdb' , 'r') as handle :

    pdb_file = handle.readlines()

print(pdb_file)

idx = 0

with open('result_python.pdb' , 'w') as handle :

    for i in pdb_file :

        i = i.strip()

        if i[0:6] == 'HETATM' :

            i = i +'  '+  charges_list[idx] +'\n'

            idx += 1

            handle.writelines(i)

        else :

            i = i + '\n'

            handle.writelines(i)

Am sure there is a bash or perl one liner for the code above.
... How to merge two files consistently line by line

The existing answer works I'm sure, however there are existing tools that do this or something quite similar.

OpenBabel is able to add charges using the option --partialcharge, though the docs are not impressive.
Autodock vina has a utility script that adds charges and extends the pdb to pdbqt, which is admittedly not exactly what you want.

From the comments: @pippo1980 suggests this page which shows how to use the openbabel Python API to set charges
Openbabel models (note that it includes fromfile meaning you can supply your own):
obabel -L charges
eem    Assign Electronegativity Equilization Method (EEM) atomic partial charges. Bultinck B3LYP/6-31G*/MPA
eem2015ba    Assign Electronegativity Equilization Method (EEM) atomic partial charges. Cheminf B3LYP/6-311G/AIM
eem2015bm    Assign Electronegativity Equilization Method (EEM) atomic partial charges. Cheminf B3LYP/6-311G/MPA
eem2015bn    Assign Electronegativity Equilization Method (EEM) atomic partial charges. Cheminf B3LYP/6-311G/NPA
eem2015ha    Assign Electronegativity Equilization Method (EEM) atomic partial charges. Cheminf HF/6-311G/AIM
eem2015hm    Assign Electronegativity Equilization Method (EEM) atomic partial charges. Cheminf HF/6-311G/MPA
eem2015hn    Assign Electronegativity Equilization Method (EEM) atomic partial charges. Cheminf HF/6-311G/NPA
eqeq    Assign EQEq (charge equilibration) partial charges.
fromfile    Assign charges from file containing {'atom-name', charge} pairs
gasteiger    Assign Gasteiger-Marsili sigma partial charges
mmff94       Assign MMFF94 partial charges
none    Clear all partial charges
qeq    Assign QEq (charge equilibration) partial charges (Rappe and Goddard, 1991)
qtpie    Assign QTPIE (charge transfer, polarization and equilibration) partial charges (Chen and Martinez, 2007)

The issue here is that the charges were already calculated by other software (using a semi-empirical approach that is much better than force-field).

@Camps the indicated fromfile approach using openbabel should accommodate inputs of the sort that you indicate in your question. Maybe I am misunderstanding that you have separately computed some estimates that need to be added to a PDB file.

nice use of obabel here too https://leela-dodda-yale.medium.com/atomic-charges-from-openbabel-8943394f2d62


text:
What do you use the CLANS program for?
I am reading this paper called Birth and death in the Terminal Complement pathway, doi : https://doi.org/10.1101/2022.04.21.489004 where they have performed CLANS analysis.
I understand that CLANS performs blast and visually represents the results as clusters.
What does these clusters mean and is there meaning to the distance between them?
And in the perspective of the paper, what does the CLANS analysis tell the researchers?
Thank you, I am still a student, sorry if this is a simple question. And I have attached a screenshot from the supplementary materials of the article

CLANS (Clustering Locus-specific ANnotationS) is a bioinformatics tool that uses BLAST to compare a set of nucleotide or protein sequences to a reference database and visualizes the results as clusters, where each cluster represents a group of similar sequences, the distance between the clusters represents the degree of similarity between the sequences in each cluster. In the context of the paper you are reading, the CLANS analysis likely helped the researchers to identify functional domains within the Terminal Complement pathway, and the clusters formed by the CLANS analysis could represent different functional domains or sub-systems within the pathway, and the distance between the clusters could indicate the degree of similarity or dissimilarity between these domains.

source: Stack Exchange

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

sort: votes


text:
What are common ways to calculate gene length for TPM calculations?
I am aware of this similar question. But the accepted answer there answers how to calculate TPM given a mapping from gene name to gene length.
My question is, given an annotation file of genes (for example from Ensembl here), how do you calculate gene length for the purpose of calculating TPM on a gene count matrix coming from RNA seq data. That is, I understand the actualy calculation of TPM given gene lengths, but I do not understand what the gene lengths are.
From what I understand, to calculate gene lengths from annotation file for RNA seq data, you only consider the exons of each gene, but even so, many exons overlap. What is the standard method of dealing with this?
For example, I have read that using the length of the union of all exons of the gene is acceptable, but from what I understand TPM is meant to represent the "expected number of reads coming from gene X when randomly selecting 1e6 transcripts". So taking the union of all exons will be an estimate of this but it might under-estimate this.
EDIT: From what I can tell this question is not a duplicate of RNASeq: Normalization, stabilization, gene length and rlog
That question does discuss how to perform normalization that takes gene length into account, but it does not ask "what is gene length? How does one determine the length of each gene?" Possible answers could be like what I already mentioned: length of the union of all exons in the gene, or another possible answer is one from the linked question "length of the longest transcript isoform" but I feel that it is not really discussed in that answer/question.

I think the most common way is to use something like Kallisto which figures out what proportion of reads are likely assigned to each transcript variant, and calculated TPM  based on that.

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


text:
Using GATK4.2 to create .idx using IndexFeatureFile for multiple vcf files
I have used GATK script : gatk IndexFeatureFile --input sampleid_hc.g.vcf to create .idx files for individual vcf files. However I am trying to create a loop and submit as a slurm job to create .idx for 200 vcf files. The script I used for this is
for file in *.g.vcf; do ./gatk IndexFeatureFile -I ${file} -O ${file}.idx; done

When I do this, I keep getting an error saying bad substitution. Please assist.

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

If you're running an older version of GATK you may need to specify the input using -F (re https://gatk.broadinstitute.org/hc/en-us/articles/360037428111-IndexFeatureFile) e.g. for file in *.g.vcf; do ./gatk IndexFeatureFile -F "${file}" -O "${file}.idx"; done. Also, is this the exact command you're running? Not for f in *.g.vcf; do ... $file ... or for f in *.g.vcf.gz? Might be worth adding in an echo before the command to make sure you are 'finding' the files you want to find (e.g. for file in *.g.vcf; do echo "$file"; ./gatk IndexFeatureFile -I "${file}" -O "${file}.idx"; done)

What is the output of ls *.g.vcf in the directory you run the command from? Can your file names contain spaces? Does it work if you quote your variables correctly with for file in *.g.vcf; do ./gatk IndexFeatureFile -I "$file" -O "$file".idx; done? What is the exact error you get?

Thank you. I will try these suggestions


text:
Contamination on genome assembly
I had a question for the community.
I have a genome of a new species that has been sequenced via 150pb Illumina paired-end.
To verify the quality of the assembly I used the BUSCO pipeline which allows me to search for 99% of the BUSCO genes (gene not duplicated) present in a taxon in the genome of my species.
Globally I only miss 18% of the BUSCO genes, which means that the assembly is not too bad.
Now when I plot the G+C% and the coverage of all the scaffolds containing BUSCO genes (so a priori non-contaminating) I should find a fairly homogeneous cloud.
Instead I end up with two different coverage profiles (noted red and blue below):

I then thought that the low coverage sequencing profile around 4x might be a eukaryotic contaminant like a acarians, however when I blast the BUSCOs present on these low coverage scaffolds I get hits that are very close to the genome of my species...
Does anyone have any idea what these are?
here is an output of the contigsTable.csv where the two first row are scaffold with low coverage in the figure, and two last are scaffold with high coverage in the figure

FRC output

Coverage distribution

Here is an Icarus summary of the assembly :
Assembly
# contigs (>= 0 bp) 1038631
# contigs (>= 1000 bp) 120132
# contigs (>= 5000 bp) 22245
# contigs (>= 10000 bp) 9986
# contigs (>= 25000 bp) 1322
# contigs (>= 50000 bp) 97
Total length (>= 0 bp) 655530043
Total length (>= 1000 bp) 451497861
Total length (>= 5000 bp) 260894112
Total length (>= 10000 bp) 174183054
Total length (>= 25000 bp) 45006012
Total length (>= 50000 bp) 5956145
# contigs 214542
Largest contig 120620
Total length 519138286
GC (%) 32.94
N50 5062
N75 1624
L50 21982
L75 71030
# N's per 100 kbp 22.28

Here is the BUSCO summary as well :
2544 Complete BUSCOs (C)
2312 Complete and single-copy BUSCOs (S)
232 Complete and duplicated BUSCOs (D)
1060 Fragmented BUSCOs (F)
811 Missing BUSCOs (M)
4415 Total BUSCO groups searched
C:57.7%[S:52.4%,D:5.3%],F:24.0%,M:18.3%,n:4415

EDIT
I ran the fastq sequences with another aligner (MEGAHIT) to see if it was not because of an assembly issue and I got the same profils coverage see here :

This could be organism-specific. We don't have a lot of info so far, so I would check a few more things:
Run something like FRC_align. Check if there's a clear signal between regions flagged as suspicious by it and your coverage graph.
Is it a eukaryote? Plant? Check where mitchondria and chloroplasts are on the plot. They will have different GC/coverage signals than the rest of your assembly, which is fine, expected even.
Worried about contamination? Run kraken2 on your raw data and see what you get there.
Run RepeatMasker and check if repetive regions contribute to the effect you see.
Finally, unless your specimen is low heterozygosity, haploid, or inbred, you probably have some effect from ploidy. There are numerous ways to check that. A self-self dot plot with MUMmer and looking coverage histograms are a good start. The idea is that you should see two peaks (let's assume diploid). One for the split haplotigs, and one for collapsed.
I checked the GC distribution of 45 genomes in your order. It's kind of all over the place, but there are some bimodal ones:

What does GC vs Length look like? Here's the one from the 45 genomes mentioned above:

Thank you for the tips I will try them and let you know. It is a genome from a male Hymenoptera (so haploid genome).

I ran FRC_align but to be honest I do not really understand what I have to look into the FRC_align output (they are many). Is there one that gives you suspicious candidats?

Go by OUTPUT_HEADER_Features.txt. It will have features for abnormal coverage as well as abnormal PE alignments - given your input size. Compression is where they align too close together, stretch is the opposite. I also edited my post with some more info on genome GC densities from these genomes.

Thank you very much for that. I added the outputs I got, I do not have OUTPUT_HEADER_Features.txt... But I'm wondering is the issue does not come from the fact that there are 2 coverage profils rather than GC content ? I also added content of the file contigTable.csv to show you 4 scaffolds, maybe you will se what is going on there ?

The OUTPUT_HEADER part is just a placeholder. In your image, the file is Features.txt. I suggest if there are any more concerns of that nature, we should move this to chat, as comments are discouraged from becoming extended discussions.

Yes sure: here is a room where we can talk  https://chat.stackoverflow.com/rooms/224463/chippycentra

That is indeed puzzling. My best guess would be contamination, even if they carry some of the BUSCOs. I run all de-novo genomes through blobtools, which creates a similar plot to the one you made manually, but also adds taxonomical annotations. It can look like this:

Also, what is the proportion of the genome that is low-coverage? Does it have a significant span? I almost always see a bunch of tiny contigs with low coverage in de novo assemblies, but they usually represent a very small fraction of the genome.
-- edit --
One more thought, are you sure you don't have a tetraploid species? The two coverage peaks kind of correspond to 8 and 32x right? One is 1/4 of the other. We developed some tools to investigate this kind of problems in data.
-- edit 2 --
Just a few more comments. Coverage differences can have various reasons, but all the sequencing biases would cause a skew of the distribution, hence the apparent bimodality of your coverage (more apparent from the two 2d plots) suggests two distinct sources. In general, could be due to various ploidy levels or contamination either by bacteria, endosymbionts or even other species getting into DNA isolation. Now that polyploidy is ruled out, it boils down to contamination, and to figure out what contamination it is, I recommend the tool I liked above.

Hello, I have 553 265 scaffolds with coverage >0 in total in which 283 812 have covdepth between [0.1-15] (so more than half of scaffolds), in fact when I plot the coverage I get the plot I added on the post.

That sounds like a lot. Are the scaffolds long? What is the sum of the 283 812 scaffolds?

I also added one more thought to the answer. Is it possible you deal with a polyploid?

I don't have much experience but my common sense tells me I would focus on getting larger scaffolds with higher coverage before trying to make fireworks with the content of these scaffolds. If you can't do that for all the genome maybe you should focus on a region and analyze in-depth that region with better data. (addressed to the OP)

In my experience, sequencing data must make sense. If you have extra coverage peak that affects half of your genome, it's a rather big deal and I do think OP is on right track to investigate where it comes from.

I added summary about assembly metrics on the post. And I'm sure it is an haploïd genome since it comes from an Hymenoptera male. I tought it would be contamination but as I said, even if the genome is quite fragmented it still found 82% BUSCO genes. And when I blast the BUSCO in both high and low (two pics) coverage scaffolds, I found an Hymenoptera hit...

Just wanted to show my skepticism on the suitability of methods for assembly evaluation based on evaluating the genes obtained or in the GC percent. Especially with non-robust data. I am not sure how these methods work but seems to me that there are many requirements that could be missed leading to misleading data. But that's an interesting thread and I'm interested to learn.  Please go ahead...

I did the assembly with another assembler (MEGAHIT) and it seems that it does not come from here since I see the same coverage patterns (see figure in edit post)...

Ha, having a hymenopteran male is really useful for ploidy doubts :-). Did you use a single male for the library prep? I thinking if you could have contamination of different lab strains? Is the BUSCO evenly distributed between high and low coverage scaffolds?


text:
Measure the purifying selection for certain taxa along a phylogeny
I am posting this message because I need clarity about the analysis I want to perform.
In my analysis, I have a homologous gene in 13 species, and I would like to evaluate the selection pressure of this gene (its dN/dS) in 3 species of the tree and test if this value is significantly different from the neutral model (dN/dS = 1).
So here is the protocol I am considering, and I would like you to tell me if I am right or not:
To calculate this selection pressure, the ratio of non-synonymous to synonymous substitution rates (dN/dS) must be calculated.
For this I use CODEML implemented in ete3.
So I need :
A phylogenetic tree of the gene (performed with the AA alignment of the gene via Iqtree)

A codon alignment (performed with the DNA alignment of the gene with pal2nal)

Mark branches in the tree.

Here is the exemple where I want to get the dN/dS value of SP5,6 and 7.:
         /-SP1
        /-|
       |   \-SP2
       |
     /-|      /-SP3
    |  |   /-|
    |  |  |   \-SP4
    |   \-|
    |     |   /-SP5
  /-|      \- |  /-SP6
 |  |         \-|
 |  |           |
 |  |            \-SP7
 |  |
 |  |   /-SP8
-|   \-|
 |      \-SP9
 |
 |     /-SP10
 |  /-|
 | |   \-SP11
 |-|
   |   /-SP12
    \-|
       \-SP13

And I want to calculate the dN/dS of SP5,SP6 and SP7 in order to confirme the hypotesis that the gene is under strong purifying selection on this 3 species.
I mark the SP5, SP6 and SP7, run a free-branch model and look at the Codeml output (the dN/dS value should be the same for the 3 marked species right?) Then this value is kind of the mean dN/dS value between three leafs?
Suppose I get an omega = 0.15
If I want to test that this value is significantly different from 1, I need to mark the same sheets and run the neutral model, then compare the two models (neutral vs. free) using chi2, and if the p-value is <0.5, then the dN/dS value (0.15) is significantly different from 1.
I'm right  ?

Hi @chippycentra ... thats a funny username mixing.crossing cultural colloquialisms from very different regions, Ireland and the Northern GB! Anyway your question is fine but you need to provide some code. I get what you are doing, this is Python. Otherwise the response will be equally as vague. I do not think codeml is implemented in ETE3 ... I may have missed recent developments, but I'd be surprised if ETE3's developers have gone down that path.

Hello, ahah I did not know that it was a random username ^.^. Well I used this protocol http://etetoolkit.org/docs/latest/tutorial/tutorial_adaptation.html#test-on-branches) ete3 implemented codeml (they still use codeml but with a python code).  I edited the post in order to be more precise.

There is no way that is a random username ... anyway you are right about CODEML. I want to read how they have done this first because it would imply CODEML (UCL, UK) have teamed up with the ETE3 guys (Spain).  Its possible CODEML code was just given to ETE3 and has changed its user license. I can give you the answer in terms of CODEML outside ETE3, but inside it ... this is new to me.

Well, since ete3 simply write a script for codeml and then run this script on it, you can give my an answer with codeml and i will understand :-)

Okay, not till the weekend, its just that Biopython was the standard wrapper for this sort of stuff so I'm astonished by the ETE3 thing.

it is juste that ete3 allows to deal easily with trees I guess

In my opinion .. it is very good for handling and traversing trees and great in data pipelines. As a CODEML wrapper I never knew it existed.

@Michael Just as an aside, do you know a software that allows, from an alignment of n sequences, to calculate the dN, dS, dN/dS between each pairwise sequence? But with a robust method such as likelihood (and not a counting method).

Please see this post https://bioinformatics.stackexchange.com/questions/14801/multiple-alignment-cost-application/14813#14813 ... I'm not a fan of what I think you are proposing, but I might be wrong in what you are asking

The following is not the clearest explanation of dN/dS humankind has seen, but I hope you get the gist and core idea of likelihood is absolutely correct. I am wondering if you are thinking of a different test to the way I would do this calculation. This ain't a text book clarity answer.
Anyway what you are proposing is a likelihood ratio test (LRT) for  selection based on the parameterisation of omega. Your idea isn't quite correct ... its not dN/dS value you are looking at to assess significance, its the likelihood. This point is critical and is the most important difference in our thinking on this test.
LRT works because its a log likelihood and subtracting two log likelihood is the same as dividing them (I forget the theory - but you subtract the likelihood to work out the p-value via Chi squared). There's some theory about why the Chi-squared distribution works but I can't remember it.
You need to work out the degrees of freedom which is the parameter differences between 'free branch model' and the neutral model and this value is LOW. If the neutral model is the "global model", then its 2 degrees of freedom - IN YOUR TREE. This is because the SP you are looking at form a single clade (monophyly). From memory you stick a #1 on the branch leading into the clade, ie. defining the node and thats your 'model' on the treefile. Here you have two dN/dS one for the whole tree excluding your clade and one for your taxa. This gives you a likelihood.
Note you could stick a #1 on each branch within the clade, i.e. the treefile. That means something different - like every branch was subject to ... in this case purifying selection ... not just the common ancestor.
You now work out the likelihood of the 'global' model... i.e. remove the #1 and make sure it gives you only 1 dN/dS value.
Subtract the likelihoods and multiple by two - if its LRT (super important).
You now need to work out the degrees of freedom. If its the same model but one has an additional node with a different dN/dS ... the df = 2. If you've dumped a #1 on every branch in the clade ... I still think df = 2 - you need to check.
If the difference in likelihood results in p < 0.05 you are correct the difference is significant ... whatever that difference might be.
If the neutral model is some different calculation of omega within the population genetics process (see here) then it will be the parameter differences between your model and the 'neutral' model.
To reiterate
I like simulations so much I type in reiterations ... So again the LRT usual will usually subtract one likelihood from another and THEM multiply the probability by two ... and compare the result on the Chi squared table for the difference in degrees of freedom. And again what you are looking at is not the value of dN/dS but the difference in LIKELIHOOD. That is very important. You need to be careful with the degrees of freedom on an LRT because they can introduce an additional subtraction constant, test dependings (its ages since I've done one so I can't remember).
If the test is significant and one dN/dS is x for the clade and y for the rest of the tree, then that differences is significant and if y is >1 (bit unlikely) and x < 1 then sure thats purifying selection. Usually clades evolve via positive selection rather than purifying selection.
From memory that is my understanding ... but if you have a neutral model of omega and a purifying model of omega ... yeah sort of, but usually thats for a 'global' model, i.e. no differences within the tree. In your case your interested in a specific clade. The idea is still the same its measured through likelihood.
Its so long since I've done this, but I do caution not to confuse two different tests.

I do understand Yang et al did an MCMC for dN/dS ... what I'm thinking of here is maximum likelihood. You'd need to watch for that carefully and thats what you are referring to when you say Chi-squared, because that is LRT 100%. LRT is maximum likelihood.

Thank you very much for your time and all the explanation it wad really usefull !!! I also enjoy it since you know PAML implemented in biopython, I'm currently doing a dN/dS pairwise analysis (runmode=-2), but it doesn't seem to work, do you have an idea ? I've created a post here with more informations  : https://www.biostars.org/p/474752/

Hi @chippycentra much appreciated. I'm afraid I've never done Biostars - I don't even have an account - ... but I will have a look and post my thoughts here. PAML does work but can be a bit fiddly.


text:
gmx mdrun error: "Error in user input: Invalid input values In option s"
I'm new to bioinformatics and gromacs I tried to run the code and got this error
Program:     gmx mdrun, version 2020
Source file: src/gromacs/options/options.cpp (line 179)
Function:    void gmx::internal::OptionSectionImpl::finish()

Error in user input:
Invalid input values
In option s
Required option was not provided, and the default file 'topol' does not
exist or is not accessible.
The following extensions were tried to complete the file name:
  .tpr

I already have the file topol. I searched and found that I need to run mdrun-h to ensure that the command work but when I check
mdrun -h

I got
mdrun: command not found

try 'which mdrun', it might not simply be installed

thanks for replying, I tried it and got nothing

@ATpoint has a good point! The problem is your environment is wrongly configured

If the command line says "command not found", it means that the binary is not executable or not found in any of the paths mentioned in the $PATH environmental variable. You can check for example this tutorial for details on how it works and how to resolve problems related to this variable.
Which probably means, that you need to go a step back, to the installation instructions of gromacs and try to install the mdrun binary. However, if the mdrun was not installed successfully, it's very likely that more stuff of the package will be missing, so it might be worth reinstalling it all.
P.S. It's not obvious to me if getting mdrun it will resolve the first posted problem (I have never worked with gromacs). If you would share the link to the research you have done, that would help others to help you.

This was a great answer Kamil! We have a few unanswered GROMACS questions here, are you able to help us out at all with those? I'm trying to help take care of the unanswered queue.

I think there was a misunderstanding here. Two errors are in user5520049's post.

A .tpr is not present in the working
directory so gromacs complains, if your .tpr is topol then you have to give this file
to gmx mdrun:  gmx mdrun -s topol.
mdrun is not a command, you should have called instead: gmx mdrun -h. That's why you had command not found and not because gmx is not in your $PATH. You couldn't have run gmx mdrun in the first place otherwise.

Side tip: run gmx mdrun -h -hidden to see hidden options.

thanks for helping . the path is also was a problem but your answer is important for me . i tried to execute those commands gmx grompp -f run_$count.mdp -p topol.top -c em.gro -o run_$count.tpr -nb gpu;
gmx mdrun -deffnm run_$count -nb gpu and got the same error although topol.tpr is exists in the directory

I can assure you that gmx is in your $PATH, you could not launch any gmx commands otherwise (gmx grompp or mdrun etc.).
The problem is that you use -deffnm (which you didn't mention in your orginal post btw) with gmx mdrun which is the option to setup default input and output names. So with -deffnm un_count gromacs is looking for a tpr named: un_count.tpr and not topol.tpr. Either:

use -deffnm un_count and name your tpr un_count.tpr
do not use -deffnm and not use -s and have default topol.tpr in the working directory.
do not use -deffnm & use -s yourcustomtpr.tpr

Here is a link to good gromacs tutorials: http://www.mdtutorials.com/gmx/ , maybe you could take a look to be familiarized with this MD engine ;)

This was a great answer Jeremy! We have a few unanswered GROMACS questions here, are you able to help us out at all with those? I'm trying to help take care of the unanswered queue.


text:
Is it possible to integrate multiple gene expression datasets and use it for WGCNA?
I have 8 RNA-seq datasets and am interested in looking at genes co-expressed with a specific gene.
Among 8 RNA-seq datasets, 6 have less than 20 samples. Rather than working on each dataset individually for WGCNA, do you think I can merge all datasets and use them for WGCNA?

Yes, if the expression datasets are from the same platform, you can merge all the datasets for common genes among them. However, if they are from different platforms then some batch correction needs to be done before analysis. You can plot a PCA to check how much variability there is among datasets from different platforms and then decide whether batch correction is needed.

thanks a lot nitesh

I would try to process all data sets in a similar way, then combine them and batch-correct, e.g., using ComBat or WGCNA's own empiricalBayesLM.

thanks a lot, peter. Can I use edgeR removebatcheffect?

It should work as well.

Hi @Peter could you please help me with this https://bioinformatics.stackexchange.com/questions/19922/how-to-find-closely-related-genes-for-a-specific-gene-from-wgcna-modules


text:
Finding SNP/Indel Markers of Interspecific Variation

This is a revised version of my original question, edited to include the portions already answered by @Maximilian Press.
Goal:
Find SNP and indel alleles which are fixed between but lack polymorphism within two closely related species (I'm looking for the BLUE region of the Venn diagram)
Desired Outputs:

All SNPs/Indels of Species A which are exclusive to and lack variation within Species A

All SNPs/Indels of Species B which are exclusive to and lack variation within Species B

Scratch Materials:

Species A reference genome

Species B reference genome

Illumina short reads for Species A

Illumina short reads for Species B

My current plan is to...

Map Species B to Species A:

minimap2 -ax asm5 referenceA.fasta referenceB.fasta > BtoA.sam

Find all SNP/Indel alleles of Species B relative to Species A (Represents the entire circle in the above Venn diagram):

bcftools mpileup -f referenceA.fasta -s BtoA.sam | bcftools call -o B_compared_to_A.vcf

Find all SNPs/Indel of Species B reads relative to Species B reference (Represents the red region of the Venn diagram):

bcftools mpileup -f referenceB.fasta -b bam_listB | bcftools call -o B_compared_to_B.vcf

Perform a liftover on "B_compared_to_B.vcf" to correct for coordinate differences between Species A and Species B references (I haven't completely figured out how to do this - I want to make sure I have the general pipeline nailed down before I get into the weeds)
Write a program which extracts all loci found in "B_compared_to_A.vcf" but NOT in "B_compared_to_B_liftedOver.vcf" (This should yield the BLUE region of the Venn Diagram
Repeat reciprocally to get alleles exclusive to Species A

Does the logic for this new pipeline check out?

That seems like a valid approach to me.

I might be missing something, but does this approach not require that the coordinate systems of referenceA.fasta and referenceB.fasta be the same? I would expect a liftover step to map variant coordinates between the two. For example, the BAM headers in bam_listA will not match the fasta index of referenceB.fasta. Maybe these tools don't care about that, I don't know that off the top of my head. but I worry about coordinate frameshifting.

At the risk of  sounding dumb, is there a reason that you don't use reciprocal whole-genome alignment via e.g. minimap2 to identify the interspecific differences, which you can then check with reads in each species separately? i don't have a specific workflow fleshed out, but i would hope that the reference assemblies were good enough to get the (vast) majority of things that would be picked up by reciprocal shotgun read mapping.

@MaximilianPress could you link me to an article describing coordinate frameshifting? I am relatively new to bioinformatics and haven't encountered this concept before. I will also give minimap2 a look, I haven't heard of this software before.

not really an article-level concept (and others would use different term probably), it's just that the two genomes probably don't have the exact same coordinates (e.g. chr1 is 10.5 Mbp in one genome, and 10.6 Mbp in the other genome). If this is true, then when you call a variant in one genome, those coordinates will usually differ from the coordinates for the homologous variant in the other genome. Stated differently: it is very unlikely that two species' genomes can be aligned without gaps. Wherever there is a gap, variant coordinates will shift.

@MaximilianPress I apologize for dredging up a nearly month-old conversation, but I am still struggling with this problem. When you mentioned "reciprocal whole-genome mapping" with minimap2, what did you intend for me to do with the reads once I have mapped the Species A reference genome to the Species B reference genome? I understand you didn't have a pipeline worked out, but I'm trying to figure out the logic of the approach.

I will post as a short answer, easier to format etc.

Ok, with invaluable aid from @Maximillian Press, I've figured out how to find alleles which are variable between two species but fixed within one species.
Here is how you find SNP/Indel loci which are distinct between species A (hereafter sA) and species B (hereafter sB) but are fixed WITHIN sA:
# Maximillian Press's Contribution: Map sA reference to sB reference using minimap2, then convert to bam.

minimap2 -a sB_ref.fasta sA_ref.fasta > AtoB.sam
samtools view -bh AtoB.sam > AtoB.bam

# Sort bam by coordinates, then use bcftools to generate a pileup and call all alleles of sA relative to sB

java -jar $EBROOTPICARD/picard.jar SortSam I=AtoB.bam O=sorted_AtoB.bam CREATE_INDEX=true SORT_ORDER=coordinate TMP_DIR=$TMPDIR
bcftools mpileup -O b -o AtoB.bcf -f sB_ref.fasta sorted_AtoB.bam
bcftools call -m -v -o AtoB.vcf AtoB.bcf

 # Call all alleles of sA relative to sA.
# This represents the SNP/Indel loci which are variable WITHIN sA and thus ought to be discarded from the final output.
# "bam_listA.txt" contains names of all the coordinate-sorted bam files for species A
bcftools mpileup -O b -o AtoA.bcf -f sA_ref.fasta -b bam_listA.txt
bcftools call -m -v -o AtoA.vcf AtoA.bcf

I then use the "sed" function to replace chromosome names with their analogs (happy to give more details if anyone else needs to do this) and also removed the VCF headers prior to making a C++ program which extracts all the rows found in AtoB.vcf but not in AtoA.vcf. My final output was thus all SNP/Indel loci which vary between sA and sB but NOT within sA.
Happy to share the C++ program if anyone else needs to do this.

Problem statement
You want to find a set of SNPs in populations of two species, for each of which you have a distinct reference genome. You then want to know which of those SNPs are shared between the two genomes (there is variation in each species), and which ones are private to only one of the species (in the other genome there is no variant at that location).
To accomplish this task, you need to first find SNP variation in each genome, and then figure out which SNPs are the same between the two genomes. In other words, line up two VCF files for different genomes and figure out which SNPs are the same between them. In comparing genome assemblies, this is sometimes called a liftover operation.
Genomes can have different coordinate systems.
Let's first think about what exactly you are looking at with 2 different reference genomes. For example, here are outputs of samtools faidx for 2 different S. cerevisiae strains, S288C and Kyokai-7:
# S288C
ref|NC_001133|  230203  95  60  61
CHR2.19980913   813142  234173  60  61
CHR3.19980521   315339  1060907 60  61
CHR4.19990210   1531929 1381540 60  61
CHR5.19970727   576870  2939039 60  61
CHR6.19970727   270148  3525562 60  61
CHR7.19970703   1090936 3800252 60  61
CHR8.19970727   562638  4909411 60  61
ref|NC_001141|  439885  5481523 60  61
CHR10.19970727  745440  5928778 60  61
CHR11.19970727  666448  6686681 60  61
CHR12.19970730  1078172 7364277 60  61
ref|NC_001145|  924430  8460517 60  61
ref|NC_001146|  784328  9400452 60  61
CHR15.19970811  1091283 10197892    60  61
ref|NC_001148|  948061  11307461    60  61
ref|NC_001224|  85779   12271444    60  61

# Kyokai first 10 lines
gi|347729946|dbj|BABQ01000001.1|    5491    137 80  81
gi|347729945|dbj|BABQ01000002.1|    4704    5834    80  81
gi|347729944|dbj|BABQ01000003.1|    203046  10734   80  81
gi|347729943|dbj|BABQ01000004.1|    1036    216456  80  81
gi|347729942|dbj|BABQ01000005.1|    126714  217642  80  81
gi|347729941|dbj|BABQ01000006.1|    121900  346077  80  81
gi|347729940|dbj|BABQ01000007.1|    20385   469638  80  81
gi|347729939|dbj|BABQ01000008.1|    243458  490415  80  81
gi|347729938|dbj|BABQ01000009.1|    189876  737054  80  81
gi|347729937|dbj|BABQ01000010.1|    81973   929441  80  81
# ...
# goes on for many thousands of lines

What you will notice is that there are significant differences between the 2 assemblies. The first is almost all in chromosomes, the second is largely mapped to various chromosomes but exists in thousands of tiny contigs. There is no way to just map reads to them, call variants, and assume that you (or bcftools) can figure out what is going on.
This is possibly an extreme example, but even in the case of two genomes assembled to chromosome pseudomolecules, there is basically no chance that those two chromosomes are on the same coordinate system or have similar metadata. That means again that bcftools probably won't be able to figure it out, because a SNP at position 112929 in genome A will be homologous to position 110201 in genome B.
Liftover: a way to solve the mapping problem between 2 genomes.
You can create a mapping between the two genomes. This is 100% a solved problem in bioinformatics, and it is more or less the same as doing liftovers between different assembly versions of the same genome. I think that you can more or less just use the UCSC LiftOver tool to convert SNP coordinates, unless I am missing something. That how-to guide that I linked is going to give you much better instructions than I can give you, here is a more minimal version that might be simpler.
The only extra suggestion that I have is that the guide uses BLAT to align the two assemblies, and minimap2 is a much faster tool. You may need to fight the input formatting a little bit, but if you run into runtime issues you might make that change. This is why I suggested it in the comment, though I wasn't aware that the liftover docs were so complete.
Short pipeline:

Align reads (individuals of species A to reference A, individuals of species B to reference B).
call variants-->get VCF files using preferred workflow for each of A and B independently.
compute liftover mapping (chain file, usually) from reference A to reference B.
use computed liftover mapping to actually lift VCF over from one genome to another.
align A and B references, call SNPs between them (possibly can be done from (3) output). Ensure that you are using as reference the liftover target genome so that the VCF is informative.
For each SNP private to each species (output of 5), query the BAMs to ensure that no variation exists.

Thank you for your detailed response! I still am not sure at which point in the pipeline I perform variant calling. Do I use the liftover output in lieu of the Species A and Species B reference genomes in the SNP-calling steps to get the SNPs of each species? Here's an example of what I'm thinking of:
bcftools mpileup -f liftover  -b bam_listA | bcftools call -o ACompatibleCoordinates.vcf bcftools mpileup -f liftover -b bam_listB | bcftools call -o BCompatibleCoordinates.vcf
I would then search the output VCF files for alleles which are unique to each species

I have added a short suggested workflow that indicates what I saw as the best option. I am not sure that I totally understand your suggestion, it looks like you are passing a liftover to bcftools but I don't totally understand what that is supposed to look like? You could pass a FASTA which is the output of a liftover to bcftools I guess, or just use it as a reference in alignment in the first place. Those would likely also work, though I don't have a great sense for which way is best.

Ok, thank you! I see a problem with the proposed workflow. In my understanding, to be called as an SNP/Indel, the locus in question must have different alleles in the reference and sam/bam. I am trying to find all SNPs/indels which are exclusive to Species A and all SNPs/indels which are exclusive to Species B (i.e. completely lack variation within species). If I call SNPs/Indels of A relative to A, won't the resulting VCF file not include alleles which are 100% fixed in A? I am about to post an answer to my own question which solves this issue, though I would appreciate critique

@annabelperry in this case, can you not align the two genomes to call such variants, and then confirm post hoc that they do not have variation in the BAM files? e.g. minimap2 -a genomea.fa genomeb.fa | samtools view -bh | bcftools etc., and then scan those against whichever VCF is the reference? that should be quite simple.

Ah, thank you @Maximilian Press! I will try that and upvote your answer if it works. Else, I will be back here shortly haha


text:
Can HTSlib extract bam reads occurring in a specific region without iterating through the whole file?
I am using htslib/sam.h to write a C++ program. As part of this program, I must extract reads occurring on specific scaffolds in specific regions from a bam file. Essentially, I want to perform the following samtools command in C++:
samtools view name.bam -o name.sam scaffold:pos-pos

Since I must extract reads from thousands of regions, I do not want to iterate through the whole bam file each time I extract reads from a region. I am wondering if HTSlib can extract reads occurring in a specified scaffold and range without iterating through the whole bam. I found the function bam_itr_querys() in the htslib/sam.h documentation on GitHub, but don't know how to format the inputs for this function, or what the "idx" input is. I also saw somewhere that I would need to pair this function with bam_itr_next(), but don't know the format to do so.
EDIT: My bam files are sorted and indexed, if that matters.

Have you looked at how samtools view works with a bed file of coordinates as input?

@swbarnes2 I have not, but my issue is that I cannot use 'samtools view'. I should provide some context (I will edit the original question to include this info): in my first draft of my C++ program, I simply used the system() command to run 'samtools view' where needed. I am now wrapping the original program in an R package, and simply running the "system" command does not work. So, I need to run something with the same functionality as "samtools view" from within C++ itself.

@annabelperry I can strongly recommend looking around at the internals of samtools if you need to work with htslib. Samtools is very well put together and I've found scrolling through its code to be a great way to gain insight into how to structure my own programs using htslib.

#include <stdio.h>
#include "htslib/sam.h"
int main(int argc, char *argv[])
{
  samFile *fp = sam_open(argv[1], "r");
  hts_idx_t *idx = sam_index_load(fp, argv[1]);
  bam_hdr_t *h = sam_hdr_read(fp);
  bam1_t *b = bam_init1();
  hts_itr_t *itr = sam_itr_querys(idx, h, argc < 3? "." : argv[2]);
  while (sam_itr_next(fp, itr, b) >= 0) {
    if (b->core.tid < 0) continue;
    printf("%s\t%d\t%d\n", h->target_name[b->core.tid], b->core.pos,
      b->core.pos + bam_cigar2rlen(b->core.n_cigar, bam_get_cigar(b)));
  }
  bam_destroy1(b);
  hts_itr_destroy(itr);
  bam_hdr_destroy(h);
  sam_close(fp);
}

The usage would then be ./a.out alignments.bam chr1:123-456, assuming your BAM file is sorted an indexed. Importantly this is making use of the index via sam_itr_querys() so it's not iterating over the entire file.

So I'm looking for a solution which doesn't require iterating through the whole bam file. Unless I am missing something, this solution requires that you iterate through the whole bam file and add "if-else" statements within the "for" loop to extract out your reads of interest. I am trying to avoid this, as it takes so much time. I'm looking for a solution which searches the binary bam file for reads in specific scaffolds/regions WITHOUT iterating through every read

@annabelperry That isn't reading through the entire file. sam_itr_querys is accepting a string like chr1:123-456 and then making an iterator over just those using the index (idx). This is all that samtools is doing internally. I've added a bit onto the end of this message to hopefully clarify that.

Thank you, @DevonRyan, that clears it up!

Have you considered using a sorted and indexed bam file? e.g.
samtools sort in.bam > out.sorted.bam
samtools index out.sorted.bam

This makes arbitrary region lookups trivial. I am more familiar with pysam (HTSlib for python), which has the AlignmentFile.fetch() functionality.
Surely there is an analog in HTSlib API, though I don't know what to search to get it. It is in samtools as demonstrated here.
This answer seems to put the components together that you are interested in, among them sam_index_load(), sam_itr_querys().

My bam files are already sorted and indexed - I will check out this other answer


text:
Unable to open .bam file in C++ using SeqAn due to 'seqan::UnknownExtensionError'
I am trying to open .bam files in C++ to extract reads occurring at specific scaffolds and loci. I essentially want to call "samtools view sample.bam -o sample.sam scaffold:pos-pos" from C++. I have tried system("samtools view sample.bam -o sample.sam scaffold:pos-pos") and can explain the issues if you want to know, but for the sake of conciseness I will leave the details out of this post. So, I am using SeqAn. When I run my script, I get the following error:
terminate called after throwing an instance of 'seqan::UnknownExtensionError'
  what():  Unknown file extension of copyPolly_SRR6511930.bam: iostream error

I am running my program on my institution's supercomputer clusters (Texas A&M Terra cluster), to which I do not have sudo access. I have loaded the following modules:
GCC/9.3.0
SeqAn/2.4.0-GCCcore-9.3.0
zlib/1.2.11-GCCcore-9.3.0

Here is my minimal script:
#include <seqan/bam_io.h>
#include <zlib.h>

using namespace seqan;

int main(){
    CharString bamFileName = "copyPolly_SRR6511930.bam";
// Open input file.
    BamFileIn bamFileIn(toCString(bamFileName));
    if (!open(bamFileIn, toCString(bamFileName))){
        std::cerr << "ERROR: Could not open " << bamFileName << std::endl;
        return 1;
    }

    unsigned numUnmappedReads = 0;

    try{
        // Read header.
        BamHeader header;
        readHeader(header, bamFileIn);

        // Read records.
        BamAlignmentRecord record;
        while (!atEnd(bamFileIn)){
            readRecord(record, bamFileIn);
            if (hasFlagUnmapped(record))
                numUnmappedReads += 1;
        }
    }
    catch (Exception const & e)
    {
        std::cout << "ERROR: " << e.what() << std::endl;
        return 1;
    }
    return 0;
}

I have tried the following solutions:
-Adding #define SEQAN_HAS_ZLIB (1) to my script's header
-Adding -D SEQAN_HAS_ZLIB to my compilation such that it reads g++ script.cpp -D SEQAN_HAS_ZLIB
-Adding -lz to my compilation such that it reads g++ script.cpp -lz
I am not dead-set on using SeqAn for this project, so if anyone knows of another C++ library which  extracts specific reads from .bam files, please share.
Update: I've never used HTSlib and it seems to have a fairly steep learning curve, so if anyone has an answer for the SeqAn error, that would still be very much appreciated.

The helpdesk for my institution's supercomputers showed me the following flags:
g++ seqan_test.cpp -DSEQAN_HAS_ZLIB -lz -lpthread

This fixed the error. I'm now running into different errors and am suspicious that I am getting them because I am using SeqAn2 rather than SeqAn3. I'm going to try HTSlib.

Htslib is the official library to access alignment files. It supports more formats. The typical way to read through a SAM/BAM/CRAM:
#include "htslib/sam.h"
int main(int argc, char *argv[])
{
  samFile *fp = sam_open(argv[1], "r");
  bam_hdr_t *h = sam_hdr_read(fp);
  bam1_t *b = bam_init1();
  while (sam_read1(fp, h, b) >= 0) {
    if (b->core.tid < 0) continue;
    int end_pos = bam_endpos(b);
    printf("%s\t%d\t%d\n", h->target_name[b->core.tid], b->core.pos, end_pos);
  }
  bam_destroy1(b);
  bam_hdr_destroy(h);
  sam_close(fp);
}

PS: the above started as a comment for "relevant but minor or transient information to" the original question. I didn't consider htslib a direct answer to the question, but because this had been changed to an answer, I added some details anyway. By the way, seqan2 has largely been discontinued. Use seqan3.

Thank you for this! I tested this script and it worked. I will try to edit it to suit my needs. Also, I am using SeqAn2 because SeqAn3 is not available on my institution's supercomputers. If HTSlib doesn't pan out or if no one can figure this error out, I will email them to request that it be installed.

I'm not quite sure why assigning a string to a variable would lead to a complaint about file extensions. From the SeqAn BAM tutorial, it looks like you need to use the BamFileIn object for reading BAM/SAM files:
#include <seqan/bam_io.h>
using namespace seqan;

int main(){
    CharString bamFileName = "copyPolly_SRR6511930.bam";
    // Open input file, BamFileIn can read SAM and BAM files.
    BamFileIn bamFileIn;
    if (!open(bamFileIn, toCString(bamFileName))) {
        std::cerr << "ERROR: Could not open " << bamFileName << std::endl;
        return 1;
    }
    return 0;
}

Edit: after seeing the updated code... I've got nothing, sorry. It all seems like reasonable code that shouldn't be generating errors.

Sorry about that, I did not include all of my code. I have updated my post with the rest of the code which led to my error.

I think SEQAN_HAS_ZLIB flag not works correctly as it supposed to be.
So you can check the seqan library to see if SEQAN_HAS_ZLIB works correctly, and seek the reason of extension error.
Seqan library is a mature library which using template to define general-purpose classes and structs.
1.The reason for extension error.
Step into the seqan library:
When you call BamfileIn: BamFileIn bamFileIn(toCString(bamFileName));
An object of BamFileIn struct is created, and bamFileName file is opened by the struct construction function.
And when this object is created, some extra information is wrapped with the object to check the type of file.
This information include suffix categories. In bam_file.h header file.
template <typename TDirection, typename TSpec>
struct FileFormat<FormattedFile<Bam, TDirection, TSpec> >
{
#if SEQAN_HAS_ZLIB
    typedef TagSelector<
                TagList<Bam,
                TagList<Sam
                > >
            > Type;
#else
    typedef Sam Type;
#endif
};

As you can see, if SEQAN_HAS_ZLIB is setted, "Bam" file will be added to TagSelector which will be wrapped into the BamFileIn struct and redefined by typedef as format.
Then _open() function is called to open Bam file, and extension is checked in guess_format.h
template <typename TFilename, typename TTagList>
inline bool guessFormatFromFilename(TFilename const &fname, TagSelector<TTagList> &format)
{
    typedef typename TTagList::Type TFormat;

    if (value(format) == -1 || isEqual(format, TFormat()))
    {
        // if tagId is set to -1 (auto-detect) or the current format (TFormatTag) then test for TFormatTag format
        if (guessFormatFromFilename(fname, TFormat()))
        {
            assign(format, TFormat());
            return true;
        }
    }
    return guessFormatFromFilename(fname, static_cast<typename TagSelector<TTagList>::Base &>(format));
}

template <typename TFilename, typename TFormat_>
inline bool guessFormatFromFilename(TFilename const & fileName, Tag<TFormat_> /*formatTag*/)
{
    typedef typename Value<TFilename>::Type                                     TValue;
    typedef ModifiedString<TFilename const, ModView<FunctorLowcase<TValue> > >    TLowcase;
    typedef Tag<TFormat_>                                                       TFormat;

    TLowcase lowcaseFileName(fileName);
    for (unsigned i = 0; i < sizeof(FileExtensions<TFormat>::VALUE) / sizeof(char*); ++i)
        if (endsWith(lowcaseFileName, lowerString(FileExtensions<TFormat>::VALUE[i])))
            return true;

    return false;
}

The above code is two instance of  guessFormatFromFilename function.
It's pretty complex here, but if SEQAN_HAS_ZLIB is set correctly, TagSelector will include bam and sam, and the first guessFormatFromFilename above will return true in if condition, and the error will not happend. But if SEQAN_HAS_ZLIB is not working, second formt guessFormatFromFilename function will be called and only "sam" suffix is exist in TFormat, and it will be exhibited by FileExtensions<TFormat>::VALUE[i]. You can using std::cerr to print the suffix on your screen to see if there is only "sam" suffix is exist in TFormat.
Althrough you add -DSEAN_HAS_ZLIB, I think you could check the code in library to see if it really works fine.
If SEQAN_HAS_ZLIB is working, you could check other variables or functiones related with open bam process.These function are in bam_file.h, formatted_file.h, fundamental_tag.h and so on.


text:
Diploid Consensus Sequences
My goal is to generate a diploid consensus genome. By "diploid consensus", I mean that I want to merge all supporting reads for a genomic region into a single haplome such that the resulting sam file contains TWO sets of reads, one from each haplome.
I have Illumina short-read files from the NCBI's SRA. These short-reads were obtained from a diploid species.
I am worried that, by simply using Picard's MarkDuplicates with the "REMOVE_DUPLICATES" option set to "TRUE", I will retain only the homolog with the highest quality and output a haploid version of my genome.
Here is my current workflow:

## Map Illumina reads to the reference genome. Include read IDs
bwa mem -M -R '@RG\tID:BioSample.LibraryName\tSM:BioSample\tLB:LibraryName\tPU:SingleFlowcell.SingleLane\tPL:ILLUMINA\t' reference.fasta SRR[number]_1.fastq SRR[number]_2.fastq > SRR[number].sam

## Sort Illumina reads by coordinate
java -jar $EBROOTPICARD/picard.jar SortSam I=SRR[number].sam O=sorted_SRR[number].sam CREATE_INDEX=true SORT_ORDER=coordinate TMP_DIR=$TMPDIR

## Remove duplicate reads from sam file
## *Note: How do I remove sampling error reads?
java -jar $EBROOTPICARD/picard.jar MarkDuplicates REMOVE_DUPLICATES=true REMOVE_SEQUENCING_DUPLICATES=true I= 'sorted_SRR[number].sam' O= nodup_SRR[number].sam METRICS_FILE= metrics_SRR[number].txt TMP_DIR=$TMPDIR

## Convert duplicate-free sam file to bam file
samtools view -S -b nodup_SRR[number].sam > nodup_SRR[number].bam

## Re-align reads around indels to reduce slippage.
gatk LeftAlignIndels -R reference.fasta -I nodup_SRR[number].bam -O aligned_SRR[number].bam

Does anyone know of a program which will whittle my sam files down to ONLY two homologs per region?

What you are asking for doesn't make sense. How on earth is any software supposed to determine if a discrepancy is a polymorphism or an error if you remove all the supporting reads?

@swbarnes2 It's not that I want to simply remove the supporting reads - I should clarify my intent. My end goal is to use the consensus genome as input in a homebrew microsat-calling software. I don't want this software to call sequencing errors as microsats, so I want to merge all reads of the same sequence into a single consensus read based on base quality scores, generating two sets of sequences which are the pipeline's best guess at the individual's haplomes. I will edit my original post to make this clear.

@swbarnes2 I want to do what is described in this article, but I don't have PacBio long-reads to compliment my Illumina short-reads: https://doi.org/10.1093/bioinformatics/bty279

FWIW, I do not think that removing lines from a bam file is what you want to do at all.  I think you need to be calling variants with phasing data, and using that vcf to make a modified fasta.  A quick google turned this up  https://www.biostars.org/p/225262/

This is a little hard to answer as is, could you give an example of the kind of data you want to get out? Do you want to phase the haplotypes, like in the linked paper? I think swbarnes2 has the right of it, that calling variants is the way to go. Also- duplicates are usually defined as reads with the same coordinates; unless you have very high coverage this is unlikely to remove more than a very small fraction of informative reads (they are more likely PCR or optical duplicates, i.e. artifacts).

@MaximilianPress Ah I have failed to clarify my end goal. I am designing a software for microsatellite detection from short-reads (no such software currently exists which meets my needs). I only want to call variants in "real" reads rather than duplicates. However, I have found an alternative to diploid consensus creation - I am instead going to run my software on short-reads (including duplicates) and use variant frequencies to filter out false positives.

How to separate reads, one for each reference haplotype.

Map all the reads against the reference genome
Use GATK to call SNPs
Use GATK tool [FastaAlternateReferenceMaker][1] to build the alternative reference
You now have the two "haploid references" (quotes needed because you have to take into account that they will not necessarily reflect the exact haplotypic structure, it could be a mix of the two haplotypes)
Align the reads against the two references, and attribute each read to the reference to which they map better. My colleagues developed for a different purpose (detection of Allelic Imbalance) a pipeline that includes the step for determining to which of two references a read belongs. The tool is called [BayesASE][2]

In general you may study all the work that has been done in the field of allelic imbalance and allele specific expression, where this problem has been encountered and solved.
[1]: https://gatk.broadinstitute.org/hc/en-us/articles/360037594571-FastaAlternateReferenceMaker
[2]: https://github.com/McIntyre-Lab/BayesASE


text:
BWA: Detecting Variation between Reference Genome and Short-Read Sequences
I need to identify all loci in the short-read sequence at which the number of microsatellite repeats (i.e. number of copies of "AA," "GTC," etc.) differ from the reference genome, as well as the locations at which the number of repeats is the same in both the reference genome and the short-read sequence.
I used the Burrows-Wheeler Aligner (BWA mem) to map a high-coverage short-read sequence (obtained from the NCBI's short-read sequence archive) to a reference genome. The output is in .sam format. I have also used a separate program to identify the loci in the reference genome at which microsatellites occur.
I would like to identify all the loci in the short-read sequence at which the microsatellite length and loci differ from the reference genome. Does anyone know any tools or packages I could use to read a .sam/.bam file of a short-read sequence mapped to a reference genome and identify specific loci at which the short-read sequence differs from the reference genome? I am using RStudio and have access to my university's supercomputer clusters.
For info on microsatellites, see here: https://en.wikipedia.org/wiki/Microsatellite#:~:text=A%20microsatellite%20is%20a%20tract,locations%20within%20an%20organism's%20genome.
My end goal is to determine the loci at which microsatellite alleles tend to differ within a species, as well as the loci at which microsatellite alleles tend to be identical between all members of the same species.

Can you please provide more information about your specific problem? This will greatly help in working through a solution. There is presumably some end-goal that this STR detection is used for. No algorithm will detect all variants, most will produce some amount of false negative and false positive results.

I haven't done any whole-genome STR analysis from NGS data myself, but are aware of others that have used lobSTR for this. There's also a recent paper [here] that compares a few different STR analysis packages (i.e. RepeatSeq, LobSTR, HipSTR, GangSTR). Here's the concluding paragraph:

In conclusion, all these tools are built to genotype STRs but
have different strengths and weaknesses. Based on our analysis there is
no clear overall winner. RepeatSeq and HipSTR are the best when
considering genotyping error rate even with low coverage. On the other
hand, GangSTR has an advantage because it is the only tool among them
that can call alleles longer than the read length but shows higher error
rate unless looking at only the enclosed class of reads, which in turn
would lose the GangSTR's advantage of picking up long genotypes. In
addition, GangSTR is the newest tool and so comes with reference files for
different reference builds that are periodically updated according to
the tool's webpage. The correct choice of a tool and the
subsequent filtering depends on the aim of the analysis, and might be
influenced by available hardware resources and time limit for running
tools.

In other words, without more information about the specific problem you're trying to solve, it's not clear what the best tool to use is (this is a common issue with Bioinformatics problems).

ah I see - thank you so much for linking the paper! I will give it a read. My specific problem is that I need to identify all loci in the short-read sequence at which the number of microsatellite repeats (i.e. number of copies of "AA," "GTC," etc.) differ from the reference genome, as well as the locations at which the number of repeats is the same in both the reference genome and the short-read sequence

Honestly, if anyone knows of a software that works on R and detects indels using BWA-generated files, that would work for my purposes

Yep its a cool answer, I like it too. Your comment above is a separate question. My personal thoughts are you'd hit a RAM bottleneck pretty quickly via R and you should think about a more complex data pipeline.


text:
Getting all variants in a region from ensembl
I want to get all the variants (SNPs) in a given region of the human genome. Is there a way on how to do this with the ensembl-API? It should be a representation similar to this:
There I added dbSNP to the region details. Having this data in a JSON or XML-Format would be really helpful! Here's the link to the image.
Edit: I know that the export-function on ensembl is a thing. But then I only get the rs-Numbers. To get all the info I need I would have to run (potentially thousands) all the rs-Numbers through some additional API. Getting that info in bulk right away would be a huge time-saver.

REST API overlap region endpoint. For example.

Thats exaclty what I need! Thanks Emily!


text:
Calculating average coverage for .bam files (sequence data)
(Full discolosure that this is my first time working with sequence data, and with the bash scripting.)
I need to calculate the average coverage for any .bam file.
After some searching I wrote the following script:
# Script to calculate the average coverage of a genome sample

# compute the total length of the sample
tot=$(samtools view -H $1 | grep -P '^@SQ' | cut -f 3 -d ':' | awk '{sum+=$1} END {print sum}')
echo $tot

#compute the coverage at each point
samtools depth -ao cov_$1 $1 #your .bam file has depth calculated for each point

#calculate the average coverage
sum=$(awk '{sum+=$3} END {print sum}' cov_$1)
echo $sum
echo

avg=$(echo "$sum/$tot" | bc -l)
echo ''The average coverage of the sample $1 is $avg x.''

The logic behind the calculation of tot is as follows:
samtools view -H  gets the header of the bam file
grep -P gets all lines in this header with SQ
cut -f 3 -d ':' removes fields that are lists of more than 3 letters while -d uses delimiter ':' instead of tab and this leaves only the numbers i.e. lengths of the chromosomes
awk '{sum+=$1} END {tot=sum}' calculates the length/total bases in the chromosome
However, my output seems incorrect. Any advice on where I might be making an error (in code or thought process), or a better method would be appreciated.
One thing that worries me is that my tot variable is the total number of bases, so if the samtools depth function does not calculate coverage for each base then my sum and tot will not be in the same units and my calculated coverage will be wrong...

For a quick estimate you’re making it more complicated than necessary.

The theoretical average coverage is $\frac{n \cdot \hat l}{N}$ where $n$ is the number of reads, $\hat l$ is the average read length and $N$ is the genome size.

samtools idxstats gives you the chromosome lengths and number of mapped reads in one convenient list.

Putting this together, we get
samtools idxstats data.bam \
| awk -vreadlen=120 '
    {
        len += $2
        nreads += $3
    }
    END {
        print nreads * readlen / len
    }
'

(In this case I’ve used 120 as the read length.)
However, note that this over-estimates empirical sequencing depth, mainly due to the fact that a given read may not be fully mapped; thus the the average read length is an upper bound on the average mapped length. If you have an estimate for the average percentage of mapped bases per read, you can multiple that by the above coverage estimate for a more accurate number.

If you need a more accurate estimate, your approach is fundamentally sound (and the result is probably correct!). For low-coverage data you can speed up your analysis by dropping the -a argument to samtools depth: you don’t need to record zero-coverage bases, they don’t contribute to the total. You also don’t need to store the coverage into an intermediate file, thus reducing IO cost:
sum=$(samtools depth "$1" | awk '{sum+=$3} END {print sum}')

Some more remarks on your code:

The -P flag for grep is not supported everywhere and thus limits portability, so you want to avoid using it if possible. It turns out that you don’t need it here: grep '^@SQ' works just as well.
While the philosophy of UNIX is based on solving problems by combining individual tools that each solve one part, it is sometimes convenient to merge individual steps:

You can often combine cut and awk, in this case by setting the field separator:
awk -vFS=: '{sum+=$3} END {print sum}'

You can often combine grep and awk:
awk -vFS=: '/^@SQ/ {sum+=$3} END {print sum}'

Always double-quote your shell variables, otherwise your script will break if the values contain unexpected spaces or other special characters. This is especially important for filenames, because spaces and other special characters are valid characters in paths!
Your code is clean and readable. Well done.

With the current versions of samtools, you can also use:
samtools coverage sample.bam

The documents are here,

samtools coverage – produces a histogram or table of coverage per
chromosome

A further example from the site,
samtools coverage -r chr1:1M-12M input.bam

Output
#rname  startpos  endpos    numreads  covbases  coverage  meandepth  meanbaseq  meanmapq
chr1    1000000   12000000  528695    1069995   9.72723   3.50281    34.4       55.8

The histogram option is
samtools coverage -A -w 32 -r chr1:1M-12M input.bam

... and output available on the website link.


text:
How to extract a list of hydrogen bonds from a PDB?
I am analyzing the 7VF2 structure, from the RCSB-database. I noticed that some residues belonging to one of its chains interact with the other chains that are present within the assembly. My question is if there is a method to extract a list of interacting residues that are present in a PDB.
EDIT:

Just as an example, we have Lys 1616 belonging to ZC3H13 chain that interacts with Glu 621 of the protein virilizer homolog (the first chain)
Mol* 3D visualizer can be used for major information in the 3D View in the link page I provided

do you need alpha chain hydrogen interactions or side chains hydrogen interactions ? I believe most of the alpha helics and beta sheets emerge from apha carbon hydrogen bonds so thet'll be a huge number of interaction to be shown, but could be wrong

https://files.rcsb.org/pub/pdb/validation_reports/vf/7vf2/7vf2_full_validation.pdf : Atom-1 Atom-2 Interatomic distance (Å) Clash overlap (Å) : 2:B:1616:LYS:HA 2:B:1616:LYS:HD2 1.71 0.42 there is something wrong with k1616 as per validation report states ... see atoms nomenclature in Lys http://ligand-expo.rcsb.org/reports/L/LYS/LYS_D3L3.gif

DSSP gives you hydrogen bonds, not sure how it works for complexes though.

You can do that with GROMACS. You probably need to check the manual for this, but roughly these are the steps:
gmx pdb2gmx -f input.pdb -o output.gro -water spc
gmx solvate -cp output.gro -cs spc216.gro -o output_solv.gro -p output.top
gmx pdb2gmx -f output_solv.gro -o output_processed.gro -p output.top -ignh
gmx make_ndx -f output_processed.gro -o index.ndx
gmx hbond -f output_processed.gro -n index.ndx -s output.tpr -num hbond.xvg -hbn hbond.ndx

First, you convert the PDB file into a GROMACS topology file using the pdb2gmx command. Then you can solvate the system with water molecules and add hydrogen atoms. Then you cCreate an index file containing the groups of atoms you want to analyze. Then you can use the hbond command to calculate the hydrogen bonds.
I have not used GROMACS in a while. Hope that provides a starting point. With GROMACS you can even run a simulation and monitor the h-bonds over time.
You actually might need to run an energy minimziation first, as the PDB might not in an optimal state for the h-bonds detection.

Great comment! I will keep this question open so that more possible answers can be added. Meanwhile, thank you for the suggestion! I will try and let you know!


text:
How to reduce the occupied RAM when you are dealing with a very sparse matrix in a single-cell Experiment in R?
I'm dealing with a very large and sparse dataset and the first issues I met occurred when I tried to use quickCluster that reported me this error:
                'cannot allocate vector of size 156.6 Mb'

So, given that I cannot wait to change the RAM of my computer and I can't afford to use a cluster, I want to rely on some other strategies like some package that would allow me to handle sparse matrices. I'm thinking about sparseM but given that I don't know well this package I'd like to know how to shrink the ram allocation for these kind of matrices. Any suggestion will be very appreciated!
Link to the dataset

This is with the scran package, yes? Can you update your tags to mention this? To my knowledge scran doesn't support sparse matrix input, which will be a bit of a show-stopper unless you can get access to a system with more RAM (maybe you can pay a few bucks for an AWS instance?).

Yes it is. So you are telling me to utilize AWS?This could be a problem since it is a thesis and I do not want to pay.

You could go to GCP which is cheaper than AWS (Google Cloud). How many cores do you have ?

4 cores I presume

Ah, looks like I can't even procrastinate on StackExchange anymore without seeing work-related stuff. Oh well.
Anyway, the other answers and comments are way off. scran has supported sparse matrices for years, ever since we switched over to the SingleCellExperiment class as our basic data structure. quickCluster does no coercion to dense format unless you tell it so explicitly, e.g., with use.ranks=TRUE (in which case you're asking for ranks, so there's little choice but to collapse to a dense matrix).
You don't provide an MWE or your session information, but this is how it rolls for me:
# Using the raw counts in the linked dataset. Despite being
# called a CSV, it's actually space delimited... typical.
library(scater)
mat <- readSparseCounts("GBM_raw_gene_counts.csv", sep=" ")

# Making an SCE just for fun. Not strictly necessary for
# this example, but you'll find it useful later.
sce <- SingleCellExperiment(list(counts=mat))

library(scran)
system.time(clust <- quickCluster(sce))
##   user  system elapsed
##  3.170   0.174   3.411

This is running on my laptop - 16 GB RAM but I'm definitely not using all of it. I only go full throttle when I'm working on some real data, e.g., the 300k HCA bone marrow dataset. Check out the book for more details.
Session info below, I don't know quite enough SO-fu to know collapse it.
R version 4.0.0 Patched (2020-04-27 r78316)
Platform: x86_64-apple-darwin17.7.0 (64-bit)
Running under: macOS High Sierra 10.13.6

Matrix products: default
BLAS:   /Users/luna/Software/R/R-4-0-branch/lib/libRblas.dylib
LAPACK: /Users/luna/Software/R/R-4-0-branch/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets
[8] methods   base

other attached packages:
 [1] scran_1.16.0                scater_1.16.1
 [3] ggplot2_3.3.2               SingleCellExperiment_1.10.1
 [5] SummarizedExperiment_1.18.1 DelayedArray_0.14.0
 [7] matrixStats_0.56.0          Biobase_2.48.0
 [9] GenomicRanges_1.40.0        GenomeInfoDb_1.24.2
[11] IRanges_2.22.2              S4Vectors_0.26.1
[13] BiocGenerics_0.34.0

loaded via a namespace (and not attached):
 [1] beeswarm_0.2.3            statmod_1.4.34
 [3] tidyselect_1.1.0          locfit_1.5-9.4
 [5] purrr_0.3.4               BiocSingular_1.4.0
 [7] lattice_0.20-41           colorspace_1.4-1
 [9] vctrs_0.3.1               generics_0.0.2
[11] viridisLite_0.3.0         rlang_0.4.6
[13] pillar_1.4.4              glue_1.4.1
[15] withr_2.2.0               BiocParallel_1.22.0
[17] dqrng_0.2.1               GenomeInfoDbData_1.2.3
[19] lifecycle_0.2.0           zlibbioc_1.34.0
[21] munsell_0.5.0             gtable_0.3.0
[23] rsvd_1.0.3                vipor_0.4.5
[25] irlba_2.3.3               BiocNeighbors_1.6.0
[27] Rcpp_1.0.4.6              edgeR_3.30.3
[29] scales_1.1.1              limma_3.44.3
[31] XVector_0.28.0            gridExtra_2.3
[33] dplyr_1.0.0               grid_4.0.0
[35] tools_4.0.0               bitops_1.0-6
[37] magrittr_1.5              RCurl_1.98-1.2
[39] tibble_3.0.1              crayon_1.3.4
[41] pkgconfig_2.0.3           ellipsis_0.3.1
[43] Matrix_1.2-18             DelayedMatrixStats_1.10.0
[45] ggbeeswarm_0.6.0          viridis_0.5.1
[47] R6_2.4.1                  igraph_1.2.5
[49] compiler_4.0.0

Not a direct solution but some workarounds:

As far as I know, Seurat can work with sparse matrices.

The particular function of scran that you are using eats up quite some memory. I believe it is needed for the "normalization" step (that is how I used it anyway). Whereas the scaling normalization performed by this function is superior to crude "log normalization", you can give a try with the latter, which is far less computationally-intensive (does not do clustering). Seurat, once again, can help with this.

You can downsample your data to the extent that it fits your into your RAM.

You can give a try with Python. More and more single cell packages are written in Python, to some extent because of the problem you have experienced. For example Scanpy output is comparable to that of Seurat, although I am not sure if you can use scaling normalization with Scanpy.

I solved the issue. Anyway I will try what you suggests also because it is instructive. Thank you very much!

Please let the community know how have you solved the issue, this is important for people who will encounter the same problem. You can just "answer" your own question and "accept" the answer.

Essentially you've hit a RAM bottleneck and the calculation will slow to zero, or in this instance refuse to go forward. The way to do this normally is to parallelize the calculation across the cores of your machine. This will likely remove the RAM bottlenack, dont ask me the computer archetectural reasons for why it works - but its works.
However, my knowledge of R is minimal. I wouldn't know how to parallelize an R calculation. It is certainly doable in Perl and Python, but the calculation needs to be written to ensure parallelisation.
The other way in is to reconfigure your calculation to remove sparse matrices OR find someone doing NGS where they've confgured their machine around heavy RAM.
Looking at your calculation I don't quite understand why you need to use a specific package, it looks like unsupervised machine learning, PCA - tSNE that sort of thing and you don't need a given package to do that you just need to vectorise the inputs. If you worked out the statistical components of Scran then there are a few extremely strong R statisticians/bioinformaticians on site who would have zero problem in replicating this within a few lines of code. Its not hard in Python's Sci-kit learn either. At a guess they perform PCA and resolve it via tSNE and this give nice clear clusters.
GCP is free for 3 months so it will cost you zero in context to a single calculation.

Even PCA gives problem if it ingests this sce-object. Anyway I solved changing the computer with one more performative. Very good to know that GCP is free for three months. Thank you!

I prefer @haci 's work around, thanks are not a cool way to say thanks.but anyway

good suggestions @Michael. It's not a big dataset, 3500+ cells, so you can do PCA or tSNE on the dimension reduced data. I don't see the OP's calculation though, or what is the larger question..

.. me neither but they are doing Plasmodium (malaria, P. falciparum?) so there will be 24-hour periodicity in the RNA expression depending on what in vitro assay they are doing. That would make it a bit complicated. They also appear to do martial arts - its a tough field malaria ;-)

Not a direct solution to your problem because it does not use SCE. But if you (and others who ended on this post) are facing the RAM bottleneck when analysing single-cell data, then Scarf might be useful.
Here is the Github link:
https://github.com/parashardhapola/scarf
Documentation:
https://scarf.readthedocs.io/
Paper:
https://www.nature.com/articles/s41467-022-32097-3


text:
Question about Using MDAPackmol (packmol in python)
I would like to use the python wrapper for packmol (MDAPackmol) for my research. In order to do so, I have tried to work with the example code and pdb files provided on GitHub so I can gain an understanding of the how to use the wrapper. However, when I run the code, I receive the following error message: ValueError: invalid literal for int() with base 10: '' .
Here is the complete trace back.

The error originates at start.resname being a zero-length string. It seems to be a problem with the resname having a number and a letter and thus is a problem with Universe (the object returned by load_packmol_output). Thus, there is a problem with the residues. Any advice or insight would be greatly appreciated. Thanks!

H @bioinf99 thanks for your post. The error is simple the parser is expecting an integer and not getting one. Tracing it on the information given, is not going to be easy. The location doesn't help much. It ain't packmol thats the problem. Maybe someone here has used the wrapper ... It could lots of things ...I ain't optimistic here, but lets see. My best guess is its the 'argparse' command and you're not entering what it wants to see. In that case you'd need to look at the opening code. If that doesn't mean much to you I'd contact the developers. Just my views however.

Thanks for your reply. I agree with your assessment. I have looked at the source code but am unsure about making changes directly to it. Unfortunately, I have emailed the developer, but he has not answered.

You'd need to send it through an IDE with breakpoints (or at least that's how I'd do it, to look at all the variables and their parsing to the package) and thats lots of work. Someone on here might know, or the developer might be on site. Consider Snakemake if this doesn't work out?? Snakemake is well supported (I don't use it BTW)..

Thanks for the ideas. Definitely will look into it if there isn't an easier solution.

One thing worth checking is if test files provided crash similarly. Also, about the email radio silence, a package with 7 stars could be one that is not really intended for the general public. An alternative is raising issues on GitHub —I for one answer those quicker as watchers get notified too and is public.

Thank you for your comment. The test file does not appear to produce any errors. Thank you for your insight regarding raising the issues on GitHub. My lab really wants this to work.

Sorry, I just re-read your Q and I think you've inadvertently misdirected everyone. The first arrow points at the function/class call not at a specific argument. The error is from start.resname being a zero-length string —not sure why resname is a letter + number. It is a property of the Universe object returned by load_packmol_output. So the packmol call seems to have worked, but the residues have some oddity. So you might want to alter your Q...

Thanks. So don't worry about the instructions? That would make sense, because it is what was shown in the documentation. It's actually a problem with MDAnalysis? Thanks again.

Sorry for another question, but I am new to working through these types of errors. I was wondering how you knew it started at start.resname. Thanks!


text:
gene-gene correlation from two different Tissues
I am stuck on how to do correlation for two independent data sets with common row and column names.
A and B are datasets that contain as many rows as genes and as many columns as samples.
The rows in A and B represent a common set of genes but measured in two different tissues.
The columns represent measurements in the same 5 samples in both A and B.
I want to do a correlation between the set of genes in A and B. This is to see if the same genes in both tissues are correlated or not.
Since the matrix would be big in my actual data, I only want to retain a correlation coefficient higher than 0.5.
Here I simulate the data set.
set.seed(1)
A <- data.frame(rnorm(100),
                rnorm(100),
                rnorm(100),
                rnorm(100),
                rnorm(100))
row.names(A) <- paste0("G_", 1:100)
colnames(A) <- paste0("M_", 1:5)

set.seed(42)
B <- data.frame(rnorm(100),
                rnorm(100),
                rnorm(100),
                rnorm(100),
                rnorm(100))
row.names(B) <- paste0("G_", 1:100)
colnames(B) <- paste0("I_", 1:5)

Thank you!

You can use mapply(). ta and tb being transposed data frames of your A and B data frames respectively:
> mapply(cor, ta, tb)[mapply(cor, ta, tb) > 0.5]

      G_3       G_5       G_9      G_10      G_11      G_15      G_20      G_23      G_25      G_26      G_33      G_40      G_43      G_48      G_57      G_60
0.5346591 0.8066507 0.8379777 0.6752681 0.7221359 0.5285787 0.7333045 0.5627962 0.6533379 0.7256878 0.5996492 0.6486557 0.5108215 0.7386332 0.6596823 0.6919915
     G_63      G_72      G_76      G_80      G_81      G_90      G_97      G_98      G_99
0.5589583 0.8391917 0.7608801 0.8003665 0.6364557 0.5030968 0.7298439 0.5693024 0.5709411

However, it returns NAs in my studio!
tb = t(B); ta = t(A)mapply(cor, ta, tb)[mapply(cor, ta, tb) > 0.5]

It must be something to do with mapply() treating matrices and data frames differently, not sure. t() outputs matrices not data frames so you can do this: tb = as.data.frame(t(B)).


text:
Kraken2 or metaphlan2 report to phyloseq
I wish to convert a kraken2 or metaphlan2 report into the R package to analyze using Phyloseq.
I have performed obtained an importation and Phyloseq analysis, but the result was not correct.
library("phyloseq")
library("biomformat")
otu_Rat11 <- read_biom("Rat11.biom")

otumat_Rat11 = as(biom_data(otu_Rat11), "matrix")
OTU_Rat11 <- otu_table(otumat_Rat11, taxa_are_rows=TRUE)
str(OTU_Rat11)
#Formal class 'otu_table' [package "phyloseq"] with 2
#slots ..@ .Data : num [1:1018, 1] 942 433 3007 116765 15103 ... .. ..-
#attr(*, "dimnames")=List of 2 .. .. ..:ℎ[1:1018]"186826""33958""1578"
#"33959".........:chr[1:1018] "186826" "33958" "1578" "33959"......... :
# chr "sa1" ..@ taxa_are_rows: logi TRUE

taxmat_Rat11 = as.matrix(observation_metadata(otu_Rat11), rownames.force=TRUE)
TAX_Rat11 = tax_table(taxmat_Rat11)

physeq_Rat11 = phyloseq(OTU_Rat11, TAX_Rat11)
str(physeq_Rat11)
#Formal class 'phyloseq' [package "phyloseq"] with 5 slots ..@ otu_table:Formal
#class 'otu_table' [package "phyloseq"] with 2 slots .. .. ..@ .Data : num
#[1:1018, 1] 942 433 3007 116765 15103 ... .. .. .. ..- attr(*, "dimnames")=List
# of 2 .. .. .. .. ..:ℎ[1:1018]"186826""33958""1578""33959".............: chr
#[1:7] "taxonomy1" "taxonomy2" "taxonomy3" "taxonomy4" ... ..@ sam_data : NULL ..
#@ phy_tree : NULL ..@ refseq : NULL –

Now I try to obtain @ sam_data
sampledata = sample_data(otu_Rat11)
#Error in access(object, "sam_data", errorIfNULL) : sam_data slot is empty.

What do with sample data and how to merge several phyloseq-format datasets from different metagenomes?

Is it possible to have an example file for Rat11.biom? It's easier to debug with a reproducible result.
What do you mean by "the result was not correct"?
What's your code for merging several phyloseq-format datasets from different metagenomes?

I'm working on this procedure with some difficult. Could you please tell me how do you convert kraken2 report files into the biom format?

I had received the same error. My metadata file was in a matrix format, I converted it into a data.frame and it worked. Pls try and see if it works.