text:
Problems with BED format for ClinCNV
I have some problems with ClinCNV(https://github.com/imgag/ClinCNV):
I made a BED file from a BAM, then i made a GC count with ngs-bits (link: https://github.com/imgag/ngs-bits )
The main idea of redacting a BED file is making GC content, but i think i need to remove all other columns. Here an example:
My BED:
chr1    10041   10086   NB501410:34:H23TGAFXY:1:11102:23474:15441/1     0       +       0.4889

In description says:
chr1    12171   12245   0.4595

I tried to run script with my BED file and there is an issue:
 Rscript ClinCNV/clinCNV.R --normal myfolder/myfile.cov --out myfolder/ --bed myfolder/myfileGCannotated.bed --folderWithScript ClinCNV/

and there is an output:
"We run script located in folder /root/myfolder . Please, specify ABSOLUTE paths, relative paths do not work for every machine. If everything crashes, please, check the correctness of this path first."

[1] "START cluster allocation."

[1] "Cluster allocated."

[1] "END cluster allocation."

[1] "We are started with reading the coverage files and bed files 2023-01-02 13:49:43"

Error in fread(fileName, header = F, stringsAsFactors = F, sep = "\t") :
  File 'myfolder/myfileGCannotated.bed' does not exist or is non-readable. getwd()=='/root/ClinCNV'
Calls: ReadFileFast -> as.data.frame -> fread
Execution halted

Now I started the script with absolute paths and it works but I have a new problem:
Error in setnames(localDf <- as.data.frame(fread(fileName, skip = 1, header = F,  :
  Can't assign 4 names to a 7 column data.table
Calls: ReadFileFast -> setnames
Execution halted

Did you try using absolute paths instead of relative ones as the error message says? Something like Rscript /home/yourUser/ClinCNV/clinCNV.R --normal /home/yourUser/myfolder/myfile.cov --out /home/yourUser/myfolder/ --bed /home/yourUser/myfolder/myfileGCannotated.bed --folderWithScript /home/yourUser/ClinCNV/?

Also, are you running this as root? Why?

I tried with absolute paths and it works, but i have a new problem, check the question pls

Can't assign 4 names to a 7 column data.table looks like you need to specify a four column bed file

Please don't do that. If the original issue has been solved, then please ask a new question about any further issues.  I will post an answer with the absolute paths.

The error message is telling you:

"We run script located in folder /root/myfolder . Please, specify ABSOLUTE paths, relative paths do not work for every machine. If everything crashes, please, check the correctness of this path first."

So it sounds like you simply need to use absolute instead of relative paths. Something like this:
Rscript /home/yourUser/ClinCNV/clinCNV.R \
    --normal /home/yourUser/myfolder/myfile.cov \
    --out /home/yourUser/myfolder/ \
    --bed /home/yourUser/myfolder/myfileGCannotated.bed \
    --folderWithScript /home/yourUser/ClinCNV/

Also, don't run things as root unless absolutely necessary. if anything goes wrong or if any of your scripts or tools has a bug, you could end up breaking your machine if you run as root.


text:
How to download ligands for PDB structure
I have some PDB IDs and for every structure, I need all its ligands, so I want to automate the process.
The ligands are different for every chain:
$ cat ./7rdx_N_1N7.smi
C1[C@@]2(C)[C@@H](C[C@@H](C1)O)[C][C]([C]1[C@H]3[C@]([C@H](C[C]21)O)([C@H](CC3)[C]([C])[C][C])C)O   1N7
$ cat ./7rdx_U_1N7.smi
C1[C@@]2(C)[C@@H](C[C@@H](C1)O)C[C@H]([C@H]1[C@H]3[C@]([C@H](C[C@H]21)O)([C@H](CC3)[C@H](C)CCC(=O)NC[C][C][N+]([C])([C])[C])C)O 1N7

So, I cannot just download the ligands by its name (e.g. https://files.rcsb.org/ligands/download/ZN_ideal.sdf).
I've tried to find the structures in .pdb file for a given structure, but I see only their names:

Also, as you could see above, I'll convert the .sdf / .mol2 files to SMARTS which seems to be a simpler representation - maybe this information'll make the task easier.

I want to make a script for the process in Python, and the only way to access the RSCB PDB with a library I've found was using biopython but I found there an option to download just a .pdb file (Bio.PDB.PDBList.retrieve_pdb_file), not ligands.
I've also tried to make my own script, but still I don't know how to access ligands files automatically for every structure. The .pdb files I can access simply by ftp (e.g. https://ftp.wwpdb.org/pub/pdb/data/biounit/PDB/divided/rd/).

How to save each ligand from a PDB file separately with Bio.PDB?  ---> https://stackoverflow.com/questions/61390035/how-to-save-each-ligand-from-a-pdb-file-separately-with-bio-pdb

have a look at ModelServer API too https://models.rcsb.org/#/ and for Pyhton https://github.com/OpenAPITools/openapi-generator

there is a (new for me) https://data.rcsb.org/graphql/index.html GraphQL API : The GraphQL interface offers more flexible data retrieval, essentially making it possible to grab any piece of data from any level of the hierarchy in a single query. To use it programmatically POST your GraphQL queries under the data.rcsb.org/graphql endpoint. No idea how it works but should be the one you are looking for to download all Ligands ??

dont understand what does:  " /v1/{id}/ligand
Coordinates of the first group satisfying the given criteria."  means for the ModelServerAPI https://models.rcsb.org/#/

Thank you for your response. It seems they meant first ligand for the given structure, but not unique one. I tired it on 7RDZ and got ASP which is not an unique ligand for the structure.
Still trying to get all unique ligands.

Played a little bit on modelserver with no luck, graphQL should be easier because the test page has autocompletion  but wont work with OpenAPI generator for Python bindings. Ligands used to belong to a Chain in the old PDBs standard now they are non polymer entity or chemical compound so buffer components get mixed up with 'proper' ligands. Need to figure out how to get these entity out from API just sending entry ID and then recover them one by one (for all the Asymmetric units !?) ; not sure it'll be a,less hassle then parsing the full PDBx/mmCIF file

ADDENDUM https://www.apollographql.com/blog/graphql/python/complete-api-guide/ : Using GraphQL with Python – A Complete Guide

Thanks, I'll try to use a query I've made earlier in GraphiQL. It seems to provide the ligands in .smi format right away.
Still, I'm curious if it's possible to get the files from the webpage (.sdf, .mol2) with any GrapgQL query.

Found this on rcsb pdb : Ligand of Interest (LOI)
A small molecule is termed a Ligand of Interest (LOI) when it is the subject of the author’s research rather than a byproduct of the experimental conditions. For example, in entry 1OBV, the small molecule palmitic acid. (PLM) is designated as an LOI as it is bound inside the cavity of the the fatty-acid-binding protein EgFABP1 and is integral to the study. Now need to figure out where LOI is written

Pdb 101 has a typo example pdb is 1O8B not 1OBV

More, I am learning mmcif :  residues/ligands do not have a residue number in mmCIF

my memoruy is failing please find this post I tried to answer about 1 year ago https://bioinformatics.stackexchange.com/questions/18076/help-me-in-understanding-the-pdb-file

Ok, cannot try the API approach, no resources during vacation. Got the download and parse/split approach, added conversion to SMILES by Openbabel. I believe it is exactly the same question of : https://stackoverflow.com/questions/61390035/how-to-save-each-ligand-from-a-pdb-file-separately-with-bio-pdb from where I stole the juicy bit of the code:
## stealing from https://stackoverflow.com/questions/61390035/how-to-save-each-ligand-from-a-pdb-file-separately-with-bio-pdb

from Bio import __version__

print('BIOPYTHON Version : ' , __version__)

from Bio.PDB import MMCIFParser, PDBIO, Select, PDBList

import subprocess

filename = '7rdx.cif'

pdblist = PDBList()

pdblist.retrieve_pdb_file(filename.upper().split('.')[0] , pdir = '.')

pars = MMCIFParser(QUIET = True)

struct = pars.get_structure('7rdx',  filename)

print(struct, type(struct))

def is_het(residue):
    res = residue.id[0]

    if res != " " and res != "W":

        return True #1

    else:

        return False #0

class ResidueSelect(Select):
    def __init__(self, chain, residue):
        self.chain = chain
        self.residue = residue

    def accept_chain(self, chain):

        if chain.id == self.chain.id:

            return True #1

        else:

            return False #0

    def accept_residue(self, residue):
        """ Recognition of heteroatoms - Remove water molecules """

        if residue == self.residue and is_het(residue):

            return True #1

        else:

            return False #0

for model in struct :

    dizio = model.child_dict

    print('\n'+ str(dizio))

for model in struct:

    i = 1

    print(model, model.id)

    for chain in model:

        print(chain, chain.id)

        for res in chain:

            # print('\n'+'res : ', res, 'res.resname : ', res.resname)

            if res.id[0] != " " :

                print('\n'+'res : ', res, res.id, res.resname)

                io = PDBIO()

                io.set_structure(struct)

                io.save(f"lig_{chain.id}_{res.resname}_{i}.pdb", ResidueSelect(chain, res))

                try:
                    subprocess.run(["openbabel.obabel", "-ipdb", f"lig_{chain.id}_{res.resname}_{i}.pdb", "-osmi" ,"-O",

                                    f"lig_{chain.id}_{res.resname}_{i}.smi"])

                except Exception as excp:

                    print('\n ', excp, '\n')

                i += 1

Problem is Biopython parses mmCIFs files; that are now the official standard for RCSB PDB Submission of PDBx/mmCIF format files for crystallographic depositions to the PDB will be mandatory from July 1st 2019 onward. PDB format files will no longer be accepted for deposition of structures solved by MX techniques. A nice review about the long history of the 'new' format :[ PDBx/mmCIF Ecosystem: Foundational Semantic Tools for Structural Biology.
Unfortunately the Bio.PDB.MMCIFParser module return a a structure that is:
<Structure id=7rdx> <class 'Bio.PDB.Structure.Structure'> so all the ligands are related to PDB chains classification and not PDBx entity.
I could be wrong here, I am a novice in both Biopython and Structural Bioinformatics but this is what I've noticed, please feel free to elaborate on this and add a solution using the RCSB APIes ... before curiosity kills the cat
PS:
NEEDs to be tested with sample containing water

adding a different way to solve the question , still without using RCSB PDB database, but using Biopython.PDB.MMCIF2Dict module to start to get a feeling on how info are stored in .cif format files. Ligands are selected amongs the ones in the structure for their Mw.
Ligands are then saved using the same approach of first answer or a new one that takes advantage of Biopython.PDB.Dice module to extract and save portion of PDB. Again Openbabel is used to translate PDBs to SMILES.
Code :
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Jan  5 11:43:07 2023

@author: bob

https://bioinformatics.stackexchange.com/questions/20301/how-to-download-ligands-for-pdb-structure?noredirect=1#comment29972_20301

NEEDs to be tested with sample containing water

"""

## stealing from https://stackoverflow.com/questions/61390035/how-to-save-each-ligand-from-a-pdb-file-separately-with-bio-pdb

from Bio import __version__

print('BIOPYTHON Version : ' , __version__)

from Bio.PDB import (
                    MMCIFParser,
                    PDBIO,
                    Select,
                    PDBList,
                    MMCIF2Dict,
                    Dice
                    )

import subprocess

filename = '7rdx.cif'

# filename = '7RDX.cif'  # this will  work , but biopython saves in lower case

# pdblist = PDBList()

# pdblist.retrieve_pdb_file(filename.split('.')[0] , pdir = '.')

pars = MMCIFParser(QUIET = True)

struct = pars.get_structure('7rdx',  filename)

# print(struct, type(struct))

mmcifdict = MMCIF2Dict.MMCIF2Dict(filename)

# print(mmcifdict)

# for i in mmcifdict:

#     print('\n', i ,'   ---> ', mmcifdict[i])

print('##############################################################')

# print(mmcifdict['_entity'])

print(mmcifdict['_pdbx_nonpoly_scheme.entity_id'])
print(mmcifdict['_pdbx_nonpoly_scheme.asym_id'])
print(mmcifdict['_pdbx_nonpoly_scheme.entity_id'])
print(mmcifdict['_pdbx_nonpoly_scheme.mon_id'])

ligands_map_mmcif = list(zip(mmcifdict['_pdbx_nonpoly_scheme.entity_id'],

                             mmcifdict['_pdbx_nonpoly_scheme.asym_id'],

                             mmcifdict['_pdbx_nonpoly_scheme.entity_id'],

                             mmcifdict['_pdbx_nonpoly_scheme.mon_id']))

print('#########################################')

print(ligands_map_mmcif)

print(mmcifdict['_pdbx_nonpoly_scheme.entity_id'])
print(mmcifdict['_pdbx_nonpoly_scheme.pdb_seq_num'])
print(mmcifdict['_pdbx_nonpoly_scheme.pdb_strand_id'])
print(mmcifdict['_pdbx_nonpoly_scheme.pdb_mon_id'])

ligands_map_pdb = list(zip(mmcifdict['_pdbx_nonpoly_scheme.entity_id'],

                           mmcifdict['_pdbx_nonpoly_scheme.pdb_seq_num'],

                           mmcifdict['_pdbx_nonpoly_scheme.pdb_strand_id'],

                           mmcifdict['_pdbx_nonpoly_scheme.pdb_mon_id']))

print('#########################################')

print(ligands_map_pdb)

print('WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWw')

print(mmcifdict['_entity.id'])
print(mmcifdict['_entity.formula_weight'])

entity_Mw = dict(zip(mmcifdict['_entity.id'],

                     mmcifdict['_entity.formula_weight']))

print(entity_Mw)

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

ligands_map_pdb_red = [i for i in ligands_map_pdb if float(entity_Mw[i[0]]) >= 200 ]

print(ligands_map_pdb_red)

#########################################################   first approach

# class ResidueSelect(Select):
#     def __init__(self, chain, residue):
#         self.chain = chain
#         self.residue = residue

#     def accept_chain(self, chain):

#       if chain.id == self.chain.id:

#           return True #1

#       else:

#           return False #0

#     def accept_residue(self, residue):
#         """ Recognition of heteroatoms - Remove water molecules """

#         if residue == self.residue :

#             return True #1

#         else:

#             return False #0

# for model in struct:

#     print(model)

#     for chain in model:

#         if chain.id == 'A':

#             for res in chain:

#                 print(res, res.id, res.get_resname(), type(res))

# i = 1

# for res in ligands_map_pdb_red :

#     print('\n'+'res : ', res[3], res[0], res[1], type(res[1]), res[2])

#     print(struct[0][res[2]] , type(struct[0][res[2]]) )

#     print(struct[0][res[2]][('H_'+res[3], int(res[1]), ' ')] , type(struct[0][res[2]][('H_'+res[3], int(res[1]), ' ')]) )

#     io = PDBIO()

#     io.set_structure(struct)

#     io.save(f"lig_{res[3]}_{res[2]}_{str(i).zfill(5)}.pdb", ResidueSelect(struct[0][res[2]],

#                                                     struct[0][res[2]][('H_'+res[3], int(res[1]), ' ')]))

#     try:
#         subprocess.run(["openbabel.obabel", "-ipdb", f"lig_{res[3]}_{res[2]}_{str(i).zfill(5)}.pdb", "-osmi" ,"-O",

#         f"lig_{res[3]}_{res[2]}_{str(i).zfill(5)}.smi"])

#     except Exception as excp:

#         print('\n ', excp, '\n')

#     i += 1

########################################################   second approach

import warnings

from Bio import BiopythonWarning

def accept_residue_mod(self, residue):
        """Verify if a residue sequence is between the start and end sequence."""
        # residue - between start and end
        hetatm_flag, resseq, icode = residue.get_id()
        # if hetatm_flag != " ":
        #     # skip HETATMS
        #     return 0
        if icode != " ":
            warnings.warn(
                f"WARNING: Icode {icode} at position {resseq}", BiopythonWarning
            )
        if self.start <= resseq <= self.end:
            return 1
        return 0

Dice.ChainSelector.accept_residue = accept_residue_mod

i = 1

for res in ligands_map_pdb_red :

    print(struct, res[2], int(res[1]), int(res[1]))

    Dice.extract(struct, res[2], int(res[1]), int(res[1]), f"lig_{res[3]}_{res[2]}_{str(i).zfill(5)}.pdb")

    try:
        subprocess.run(["openbabel.obabel", "-ipdb", f"lig_{res[3]}_{res[2]}_{str(i).zfill(5)}.pdb", "-osmi" ,"-O",

        f"lig_{res[3]}_{res[2]}_{str(i).zfill(5)}.smi"])

    except Exception as excp:

        print('\n ', excp, '\n')

    i += 1

```

OK, seems ligands and PDBs are not really a trending topic here.... (or my posts are really poor).
Nevertheless the question is literally eating my brain so I went into a second direction: "the API" already skimmed here Help me in understanding the PDB file.
So inspecting and copying directly from the PDB entry 7RDX webpage the download value: https://models.rcsb.org/v1/7rdx/ligand?auth_seq_id=1005&label_asym_id=M&encoding=sdf&filename=7rdx_M_1N7.sdfI got this bit of code:
import requests

def main():
    risposta = requests.get(f"https://models.rcsb.org/v1/7rdx/ligand?auth_seq_id=1005&label_asym_id=M&encoding=sdf",
                      allow_redirects=True)
    open(filename, 'wb').write(risposta.content)
    print(risposta.status_code)

if __name__ == '__main__':
   filename = '7rdx_M_1N7.sdf'
   main()

That is equivalent to feeding this:
{
  "atom_site": [
    {
      "auth_seq_id": 1005,
      "label_asym_id": "M"
    }
  ]
}

to the 'https://models.rcsb.org/#/General/ligand-post' test webpage Request body.
now need to figure out how to get the auth_seq_id and the label_asym_id  parameters for all the
'non_polymer_entity_ids' of my targed PDB 7RDX from the https://data.rcsb.org/rest/v1/ DATA API.
Went till here:
import requests
import json

def main(entry_id: str):
    risposta = requests.get(f"https://data.rcsb.org/rest/v1/core/entry/{entry_id}/")
    parsed = json.loads(risposta.content.decode())
    print(risposta.status_code)
    print((parsed))
    print(parsed["rcsb_entry_container_identifiers"]["non_polymer_entity_ids"])
    print('\n', type(parsed))

    parsed_lig = parsed["rcsb_entry_container_identifiers"]["non_polymer_entity_ids"]
    print('\n\n',parsed_lig, type(parsed_lig))

if __name__ == '__main__':
   main('7rdx')

output :
 ['7', '8', '9', '10', '11'] <class 'list'>

The list values are the ids of the 'non_polymer_entity_ids' of my targed PDB 7RDX.
Right now I am trying to recover the the auth_seq_id and the label_asym_id  for every instance of the 'non_polymer_entity_ids' inside my mmCIF/PDBx using the DATA API.
Not sure this is the fastest and right way using the RCSB PDB API to accomplish the OP question (download all ligands of the model structure (not the ligands models)) from the PDB database. Please feel free to add more or point me towards the best solution.
EDITED:
It took a while, I am sure this is not the right way to make good use of all the databases, there should be a way to navigate the API schemas better and faster, and to get and manage the intermediate values needed to get the two keys for the download from ModelServer "Post /v1/{id}/ligand Coordinates of the first group satisfying the given criteria", anyway:
import requests
import json

parsed_lig_dict = {}

def main(entry_id: str):
    risposta = requests.get(f"https://data.rcsb.org/rest/v1/core/entry/{entry_id}/")
    parsed = json.loads(risposta.content.decode())
    # print(risposta.status_code)
    # print((parsed))
    # print(parsed["rcsb_entry_container_identifiers"]["non_polymer_entity_ids"])
    # print('\n', type(parsed))

    parsed_lig = parsed["rcsb_entry_container_identifiers"]["non_polymer_entity_ids"]
    # print('\n\n',parsed_lig, type(parsed_lig))

    # for i in parsed:
    #     print('\n', i  , parsed[i])

    for lig in parsed_lig:
        risposta = requests.get(f"https://data.rcsb.org/rest/v1/core/nonpolymer_entity/{entry_id}/{lig}")
        parsed = json.loads(risposta.content.decode())
        # print('\n',risposta.status_code, risposta.content)

        # print(parsed["rcsb_nonpolymer_entity_container_identifiers"]["asym_ids"],
              # parsed["pdbx_entity_nonpoly"]["comp_id"])

        parsed_lig_dict[lig] = [parsed["pdbx_entity_nonpoly"]["comp_id"] ,
                                dict.fromkeys(parsed["rcsb_nonpolymer_entity_container_identifiers"]["asym_ids"])]

    for lig in parsed_lig_dict:
        for chain in parsed_lig_dict[lig][1]:
            risposta = requests.get(f"https://data.rcsb.org/rest/v1/core/nonpolymer_entity_instance/{entry_id}/{chain}")
            # print(risposta.status_code)

            parsed = json.loads(risposta.content.decode())
            # print('\n',risposta.status_code, risposta.content)
            # print(chain , '--> ',parsed["rcsb_nonpolymer_entity_instance_container_identifiers"]["auth_seq_id"])

            parsed_lig_dict[lig][1][chain] = parsed["rcsb_nonpolymer_entity_instance_container_identifiers"]["auth_seq_id"]

    # print(parsed_lig_dict)

    # for i in parsed_lig_dict:
    #     print(parsed_lig_dict[i])

    # for i in parsed_lig_dict:
    #     for y in parsed_lig_dict[i][1]:
    #         print(i,parsed_lig_dict[i][0], y , parsed_lig_dict[i][1][y])

    cnt = 1
    for i in parsed_lig_dict:
        for y in parsed_lig_dict[i][1]:
            seq_id = parsed_lig_dict[i][1][y]
            chain = y
            comp_id = parsed_lig_dict[i][0]
            print(seq_id , y)

            risposta = requests.get(f"https://models.rcsb.org/v1/7rdx/ligand?auth_seq_id={seq_id}&label_asym_id={chain}&encoding=sdf" ,
                              allow_redirects=True)
            print(risposta.status_code)

            filename = f"lig_{comp_id}_{chain}_{seq_id}_{cnt}.sdf"
            open(filename, 'wb').write(risposta.content)
            cnt += 1

if __name__ == '__main__':
   main('7rdx')

This approach is very slow on my machine, the other two answers that parse the downloaded full mmCIF file via Biopython Structure Object are faster even adding the translation .sdf files to SMILES via Openbabel. I think this is because I translate the DATA API JSON formatted schema into a dictionary to extract the data I need, there should be another way to make use of that data. This way doesn't need installed Biopython though.

Here can use GraphiQL API provided by RCSB.
Here is an example
import requests
from rdkit import Chem
pdbid = '7RDX'
ls_smi = []

query = f'''{{
  entry(entry_id: "{pdbid}") {{
    nonpolymer_entities {{
      nonpolymer_comp {{
        rcsb_id
        pdbx_chem_comp_descriptor {{
          type
          comp_id
          program
          descriptor
        }}
      }}
    }}
  }}
}}'''
url = "https://data.rcsb.org/graphql"
r = requests.post(url, json={'query': query})
for i in r.json()['data']['entry']['nonpolymer_entities']:
    if i['nonpolymer_comp']['pdbx_chem_comp_descriptor'][0]['type'] == 'SMILES':
        ls_smi.append(Chem.CanonSmiles(i['nonpolymer_comp']['pdbx_chem_comp_descriptor'][0]['descriptor']))

Thanks @AustinApple and welcome to the site. A good answer

source: Stack Exchange

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

sort: votes


text:
Why are conversion output SMARTS strings are not the same for different file types of the same structure?
I've just used Open Babel to convert .sdf and .mol2 files to SMARTS, and I'm a bit surprised the output can be different for the same structure. It's confusing me that one structure, one piece of information can be different depending on the file format in which it was saved.
I converted a ligand 1N7 in two available formats (.sdf and .mol2)  of the 7RDX structure from RSCB PDB for chain N, and also for chain U.
After the conversion, the SMARTS strings for the ligand in both chains are not the same for both input files types (i.e. the difference between the same ligand, from the same chain, just in two different formats available to download from RSCB PDB):
example no. 1:
[additional @H added in places, some square brackets added]
$ cat ./7rdx_N_1N7_from_mol2.smi

C1[C@@]2(C)[C@@H](C[C@@H](C1)O)[C][C]([C]1[C@H]3[C@]([C@H](C[C]21)O)([C@H](CC3)[C]([C])[C][C])C)O   1N7

$ cat ./7rdx_N_1N7_from_sdf.smi

C1[C@@]2(C)[C@@H](C[C@@H](C1)O)C[C@H]([C@H]1[C@H]3[C@]([C@H](C[C@H]21)O)([C@H](CC3)[C@H](C)CC)C)O   1N7

example no. 2:
[some square brackets added]
$ cat ./7rdx_U_1N7_from_mol2.smi

C1[C@@]2(C)[C@@H](C[C@@H](C1)O)C[C@H]([C@H]1[C@H]3[C@]([C@H](C[C@H]21)O)([C@H](CC3)[C@H](C)CCC(=O)NC[C][C][N+]([C])([C])[C])C)O 1N7

$ cat ./7rdx_U_1N7_from_sdf.smi

C1[C@@]2(C)[C@@H](C[C@@H](C1)O)C[C@H]([C@H]1[C@H]3[C@]([C@H](C[C@H]21)O)([C@H](CC3)[C@H](C)CCC(=O)NCCC[N+](C)(C)C)C)O   1N7

In both examples, the structures look quite similar to each other, but the difference is a little worrisome, because there is no chance to check every output by hand.
Is it a normal behaviour during the conversion?
Is it related to a different notation in the input files?
Which input format is better?

Mmmhhh be sure that each of the 3 modeled ligands correspond to the full one, could be that the density is poor and pdb shows only what is really ((??) there and not what they soaked the crystal with, actually 1N7 is CHAPSO I believe a detergent

Prior to grid freezing, 3-([3-cholamidopropyl]dimethylammonio)-2-hydroxy-1-propanesulfonate (CHAPSO, Anatrace) was added to the sample (8 mM final), resulting in a final complex concentration of 8 μM. The final buffer condition for the cryo-EM sample was 20 mM HEPES pH 8.0, 150 mM K-Acetate,10 mM MgCl2, 2 mM DTT, 1 mM ADP, 1 mM AlF3, 8 mM CHAPSO.

to me your 2 strings seem different, no N in first molecule.
Here pictures, I used a PDB version of the original mmCIF file
but dindt use the RCSB PDB available one, imported .cif into
Biopython and parsed the resulting Struct object [so there could be
mistakes].
so chain N --> A and chain U --> E

as per in my comment:

Mmmhhh be sure that each of the 3 modeled ligands correspond to the full one, could be that the density is poor and pdb shows only what is really ((??) there and not what they soaked the crystal with, actually 1N7 is CHAPSO I believe a detergent

... I could be wrong though, this mmCIF format is driving me crazy

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


text:
If fastp output is not a good measure of FASTQ correctness, what is?
In the beginning of my pipeline, I just fed the paired reads (2 files) into fastp, with the default options, and assumed it would do a good job preparing the reads for the next step:  alignment
But I now checked the files in FastQC and it doesn't show much difference from before the fastp processing: Some extra yellow and red signals appear. I think the end result is even worse.
fastqc doesn't analyse paired reads files together. So it might be wrong regarding the paired quality, taken as a whole.
Anything else I should use together with FastQC, to verify the output quality?
I'm aligning a bunch of poolseqs and then comparing them to find consistency in the change of snp frequencies
As requested here is one of the examples that made me ask this:
Before fastp. only one problem noted by fastqc

After fastp. overall quality seems to have decreased

The raw data is a pool seq from illumina.

Please edit your question and give more details about your data (genomic/RNA-seq, whole transcriptome/3' end) and which fastqc modules are showing problems.

Hi fullmooninu a screenshot would be a good start, particularly if it sheds light on the QC. Would it be sufficient ... I'm sure @Cloudberry will comment.

Please give more context around your problem. For example, what are you trying to do with your pipeline and fastp?

updated the question with your requests.

What makes you think that quality of reads went down after running fastp? You have two warnings after running fastp. One is GC content, and other is length distribution. Both of them can be easily explained by the removal of adapter and low quality reads. You should also read fastp manual and see if running the default is a good idea.

Two extra warnings is worse than no warnings, right? I saw it consistently over the tests I did: more warnings. Dont remember if I saw any reds. The question is, if fastqc output is not a good measure of fastp correctness, what is?

The sequence length distribution warning is triggered if the reads have variable lengths, which you'd expect after adapter removal and/or quality trimming. I'm not familiar with fastp but apparently it tries to guess adapter sequences from the reads by default. If you know the adapter sequence, you'd probably get more accurate results by providing it to fastp. You didn't mention what type of sequencing library the data is from. If this is targeted sequencing, that could explain a mild bias in base composition.

edited. it's a pool-seq. i'm interested in validating what the cleaning process. fasqc doesnt seem to perform well, in this scenario.

Tools like fastqc are useful in helping you determine the quality of your generated FASTQ files, but they are not the determining factor in whether your data is "good", or your experiment is a success. How does your data look when run through your analysis pipeline? That is the critical question.

On a regular basis, using FastQC is quite enough as to assess the quality of FastQ-formatted data. It gives you plenty of details as can be seen in the screenshots you shared.
As others have commented  or suggested, there is  nothing wrong with having more 'yellow' warnings after processing the FASTQ files using fastp. At lest not as long as what you observe in the data looks fine and the warning has some explanation.
The sequence length distribution warning is very often present post-processing since at the begging all your reads are of the same length but this is not the case any more after you trim the adapters and the low quality bases at the end of the read for example. You have to go to the plots specific to that section and make sure all looks fine. It should be the case.
The GC sequence content warning. It depends a lot what organism you are sequencing, which regions of the genomes, etc. Importantly, what is the expected GC content of those regions. So you can compare to these observed numbers. Once again the warning can be easily triggered after the processing, since for a start Illumina sequencing has known GC bias
Finally, I think actually your reads look better  after the processing with fastp. You are showing us the per base sequence quality and before, there are many reads with bases on the yellow and even the red sections of the quality scale. But this is not the case any more after the processing !

Thanks. I still suspect a tool like fasqc, pair ended aware, will give better results.

You're wrong about fastqc being PE aware as well as fastp not being PE aware. Are you even reading the manuals, @fullmooninu?

I am always amazed by users new to the field and then calling out on of the most extensively used tools in all of bioinformatics. There is nothing wrong with your data in the first place. Quality does not go down after trimming as seen in ghe per-base boxplots, the length warning is expected due after any trim and GC content is just deviation from a theoretical distribution. There is no basis to call results worse or regarding fastqc as not accurate.

@RamRS i didn't find out that in the manual. can you point it out for me, please?

@fullmooninu Literally in the second usage example here: https://github.com/OpenGene/fastp#simple-usage

And it is common knowledge that fastqc analyzes FASTQ files independently. It has no options for PE, you can check using fastqc --help 2>&1 | grep -ic pair: This will give you 0.

@RamRS I appreciate you trying to be helpful. But I'm literally complaining about fastqc not being PE aware, unlike fastp. And am looking for a replacement for fasqc. please delete your comments to decrease entropy.

@fullmooninu in that case, you might want to edit your question since you state "fastp doesn't analyse paired reads files together" which is presumably why Ram linked you to the fastp manual which shows that it can actually process paired end reads together.

@fullmooninu I recognise a lot of time has been invested here and it's a great question. Short-read QC is not my area, however as general advice sometimes different viewpoints may not be helpful in the short term. However, in the longterm a different perspective may help you. Ram RS is a very experienced bioinformatician, but simply Ram RS's  approach is direct and this may cause unintentional offence if there's a misunderstanding in the question,  cultures or possibly issues such pressure to complete a project (which we all understand).

@fullmooninu you said "I still suspect a tool like fasqc, pair ended aware, will give better results".  Now you say you're complaining about it not being PE aware. Am I missing something? And yes my tone is quite blunt so I apologize for that.

@terdon corrected, thank you. To the others. Repeating: i'm looking for something to validate fastp's work. Probably like fasqc but that suppports PE. Having more warnings is not a good sign.


text:
Is there a package that can filter the `raw_feature_bc_matrix` of 10X's `cellranger multi` pipeline output?
There is no filtered_feature_bc_matrix in 10X Genomics' cellranger multi pipeline output when multiplexing analysis was performed.
There are indeed sample_filtered_feature_bc_matrix for each sample after demultiplexing. But there are still some No tag-assigned cells that are not included in any sample, and I have to filter those cells before analysis.
Is there a package that can filter the raw_feature_bc_matrix of cellranger multi pipeline output?
Is this package scCB2(https://github.com/zijianni/scCB2) suitable for this purpose?
I tried scCB2 with the following parameters, I am not sure whether these parameters will produce similar output with cellranger’s filter method.
library(scCB2)
library(SummarizedExperiment)

for (i in data_file) {
    print(i)
    tmp =  QuickCB2(dir = i, FDR_threshold = 0.01, lower = 100, Ncores = 12)
    DropletUtils::write10xCounts(paste0("./data/filtered/", basename(i)), tmp, version = "3")
}

I got 83836 cells after filtering, but the tag_calls_per_cell.csv file of the multiplexing_analysis only contains 19427 cells.
It seems I have to use a more stringent filter.
Is there a python package that can do a similar thing?
Thanks

I can work around this by only selecting cells in the assignment_confidence_table.csv. (https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/output/cellplex?src=social&lss=facebook&cnm=soc-fb-ra_g-program&cid=7011P000000y06x/#:~:text=Description-,assignment_confidence_table.csv,-A%20table%20that)


text:
Should be RT-qPCR values standardized before PCA analysis?
I come from this question of gene exression analysis: Should PCA be standardized for gene expression?
My experiment is based in 151 x 51 (individuals/samples x genes), in which, patients are subjected to 3 possible groups (~ balanced groups), so we have 2 measurements of gene expression: one in the baseline, and one after the intervention. The gene expression is actually summarised in an scalar (number) called fold-change. Fold-change is calculated through the Ct (cycle threshold), which is the number of the cycle in which you are able to detect a certain amount of the target gen. This value or cycle is normalised using Ct values of 3 reference genes which are supposed to remain stable, and each patient is measured using oneself baseline value, in the same run, to remove run-to-run variation. As you can see I am working with RT-qPCR and relative quantification values.
As you can see, "normalization" or "standardization" is performed, not in the strict meaning of mathematical foundation, but to normalize quantitation of targets in the experiments. I was wondering if, taking into account that normalization is done, when I try to perform PCA to see how these genes behave, if they form cluster, if the intervention received has influence of the expression, patterns...
Reading about, I found quite divergent opinions about scaling or standardizing the results. The conflict, as far as I can see in my own experiment, lies in all my variables are genes, measured with the same high-throughput analyzer, different batches, but normalized with reference genes, and their own baseline values.
On one hand I found quite homogeneous working with scaled values, but on the other hand, I am working with values that seem to be working on the same scale, and have been through several filters to normalize results. Not sure if executing more mathematical operations could modify in the wrong way. What is your opinion about it? (all is about using the covariance matrix or correlation matrix based PCA, as far as I understand)
Furthermore, the dataset contain missing values, so I am working with the missMDA package to impute them. In case I scaled the data, what should I do first: scale or impute the missing values?
The steps to calculate fold-change:

Update 1: About the technology used, is a OpenArray plate (analyzed with QuantStudio™ 12K Flex Real-Tim) which contains subarrays with customised target genes.It is a 3,072-reaction high-throughput real-time PCR. The reaction are individually performed (as far as I know, getting Ct for each sample and assay). About the 2nd point, it is 2x time series study (baseline and post-intervention)- But if you see the equations to calculate the fold-change, I am going to obtain a single value of gene expression combining the 2 Ct or measurements for each patient (as I told: baseline and post-intervention).
About the "genes aim" you can see that my interest is purely human genes, after an intervention (inflammation-related genes).
About the pippeting issue. As you know is never discarded, but the internal controls or reference genes (as I am adressing them) goal is to normalize results and remove potential quantification errors. However, as you said the RNA concentration (A260) and purity were calculated spectrophotometrically from peripheral blood samples origin
As far as I know, the terms correlation and covariance related to PCA, are meant to discrimanate the standardization process (in correlation matrix) and not standardized matrix (scale function in R). Not sure if I understood this is what you imply?
My guess is nowadays everything is RNAseq and microarrays, whose mechanisms, and therefore, statistical preprocesses are different from RT-qPCR. But there are papers of PCA working wiht data coming from RT-qPCR. Shouldn't it work though?
Update2
Covariance and correlation: "Covariance analysis of PCA: PCA is dimensionality reduction based on maximising multi-dimensional variance so that would be unusual, but correlation matrix of PCA is understandable. Normally covariance analysis is performed directly on standardised data."
Let me insist on this topic a little bit more.I am not sure if I understood correctly what you stated, but the information I've reviewed seems to say otherwise.
As far as I was able to find, https://aedin.github.io/PCAworkshop/articles/b_PCA.html : "PCA was computed as a singular value decomposition (SVD) of a column centered, scaled matrix. This was PCA of the correlation matrix. If the matrix is centered but not scaled, it is PCA of the covariance matrix." I am not sure if I understood correctly what you stated, but the information I've reviewed seems to say otherwise
"The correlation matrix is the standardized version of the covariance matrix. Analysing the correlation matrix is a useful default method because it takes the standardized form of the matrix; therefore, if variables have been measured using different scales this will not affect the analysis. Often you will want to analyse variables that use different measurement scales. Analysing the correlation matrix ensures that differences in measurement scales are accounted for. In addition, even variables measured using the same scale can have very different variances and this too creates problems for principal component analysis. Using the correlation matrix eliminates this problem also. There are statistical reasons for preferring to analyse the covariance matrix (the reason being that correlation coefficients are insensitive to variations in the dispersion of data whereas covariance is and so produces better-defined factor structures (Tinsley & Tinsley, 1987)) and generally the results will differ from analysis on the correlation matrix. However, the covariance matrix should be analysed only when your variables are commensurable."
Conclusion: shouldn't it be the covariance matrix the one without standardized values, because variables share same scale?
Missing values: this is an issue, at least with several functions in R (prcomp(), princomp())which don't work with NAs in the matrix.
About this statement: "The key question, which is standardisation and I would say it is essential for RT-qPCR derived data at a minimum to at least investigate consistency of signal against a control. That minimum is what the second equation describes (why it works is complicated)."
I don't really understand what you mean standardize at a minimum. The endogenous control genes we are using are individualized for every sample. Let say we are investigating the fold-change in the individual i and the variable k. We are using endogenous controls in the i individual to "standardize" values. The

This is as good an answer I can give based the information. The issue here is that the experiment is derived via RT-qPCR and thats not the same as RNAseq type approach that the normal basis for questions on this site. Generally, we sort of assume RNAseq for gene expression where all the gene expression within a sample are automatically standardised against each other. Thus PCA can be performed directly on that data type. I will summarise this at the end because it's the key point I'm trying to make.
In your experiment there are 51 separate reactions per patient (or is 51 the total time series e.g. for one gene??) and it's not clear whether this constitutes a time series (but there must be some time-series aspect to the study: I assume the 'baseline' is at least one point in a time series). That is a lot of opportunity for variation within one patient sample, even if it is the same time point. Maybe the tech doing the RT-qPCR makes a pipetting error for 1 in 10 samples - the results could collapse. The other issue is whether the RT-qPCR is internally standardised, i.e. are you assessing two amplicons per PCR and the same control amplicon is used for each different gene under investigation? This traps the 1 in 10 hypothetical error rate because that is the first point of standardisation. The question is whether a further standardisation is required beyond that, i.e. per patient - thats tricky but basically its difficult to avoid because of variation between different samples from the same patient (again you must be doing some sort of time series). So for example the RNA extraction could have different amounts of blood etc ... However, it needs some thought because here you are standardisation a standardised result against a backdrop of an unstated hypothesis (comments).
Summary The single point I want to make is that large scale RT-qPCR is quite different from the standard RNA-seq approach on this site from a statistical perspective, i.e. it appear to be 51 reactions per sample whereas normally the approach is 1 reaction per sample which generates all the expression data across thousands of genes, thus PCA is a fantastic analysis for this data-type. This difference creates a layer of statistical complexity and I what I would do is carefully consider my hypothesis and formulate an analytical plan in accordance with the aims.
The question is whether a simple standardisation could be performed and then shift straight into PCA? Yes. Is that sufficient? I dunno. I think you might be avoiding time-series analysis by standardisation and then subtraction of the baseline which is then placed into PCA. At a minimum I would perform two PCA analyses one for baseline and one for post-treatment.
Generally, a complex standardisation (like I'm inferring) is not cool if the next stage of the analysis is PCA.
Final thoughts Although PCA is considered the primary analysis here. I think it takes serious thinking given all the available information to formulate an analytical plan from my experience in this type of data, but as a rough guide (given I dunno the experimental aims) a complex standardisation and then I personally would consider covariance thereafter for a time series rather than PCA as general approach here. However, two PCAs before and after  using a simple standardisation and compare presence/absence, that could work and very clear results would be accepted.

Note Just to note that "correlation matrix based PCA" is basically PCA, it's just trying to sort out the groups. However "covariance" is a different analysis to PCA and statistically and the aims are very different. Also just to mention, but I assume you know this, "correlation" is different from "correlation matrix based PCA".
Missing values
Missing values for PCA isn't an issue - you just miss them out. Missing values in a time series can have lots of different solutions, the key is justification.

Update
Covariance analysis of PCA: PCA is dimensionality reduction based on maximising multi-dimensional variance so that would be unusual, but correlation matrix of PCA  is understandable. Normally covariance analysis is performed directly on standardised data.
Overall, I am now clearer one part of the calculation: this is essentially an ad hoc statistic to assess fold-change of the raw Ct via a single /minimal PCAs following one or more standardisations.
What you are measuring is the optimised multi-dimensional variance of the fold-change of raw Ct after at least one or else numerous standardisations. So this is not the actual fold-change of Ct, but variance of the fold-change derived via Eigen values.
If you are clear about that and what this represents and the reader is clear about this then it would be up to the reviewers.
The key question, which is standardisation and I would say it is essential for RT-qPCR derived data at a minimum to at least investigate consistency of signal against a control. That minimum is what the second equation describes (why it works is complicated).


text:
How to concatenate 2 `mdata` of `muon` package?
I used muon(https://github.com/scverse/muon) for single cell Multimodal omics analysis.
How to concatenate 2 mdata from different samples? How can I combine more than 2 samples for analysis together?
import anndata as ad
import muon as mu
import os

mdata = mu.read_10x_h5(os.path.join(samples[0], "filtered_feature_bc_matrix.h5"))
mdata.var_names_make_unique()
mdata2 = mu.read_10x_h5(os.path.join(samples[1], "filtered_feature_bc_matrix.h5"))
mdata2.var_names_make_unique()

mdata3 = mdata.concatenate(mdata2, batch_key='sample_id')
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[25], line 1
----> 1 mdata3 = mdata.concatenate(mdata2, batch_key='sample_id')
      2 mdata3

AttributeError: 'MuData' object has no attribute 'concatenate'

mdata3 = ad.AnnData.concatenate(mdata, mdata2)
mdata3
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[26], line 1
----> 1 mdata3 = ad.AnnData.concatenate(mdata, mdata2)
      2 mdata3

File ~/miniconda3/lib/python3.8/site-packages/anndata/_core/anndata.py:1770, in AnnData.concatenate(self, join, batch_key, batch_categories, uns_merge, index_unique, fill_value, *adatas)
   1767     adatas = adatas[0]  # backwards compatibility
   1768 all_adatas = (self,) + tuple(adatas)
-> 1770 out = concat(
   1771     all_adatas,
   1772     axis=0,
   1773     join=join,
   1774     label=batch_key,
   1775     keys=batch_categories,
   1776     uns_merge=uns_merge,
   1777     fill_value=fill_value,
   1778     index_unique=index_unique,
   1779     pairwise=False,
   1780 )
   1782 # Backwards compat (some of this could be more efficient)
   1783 # obs used to always be an outer join
   1784 out.obs = concat(
   1785     [AnnData(sparse.csr_matrix(a.shape), obs=a.obs) for a in all_adatas],
   1786     axis=0,
   (...)
   1790     index_unique=index_unique,
   1791 ).obs

File ~/miniconda3/lib/python3.8/site-packages/anndata/_core/merge.py:873, in concat(adatas, axis, join, merge, uns_merge, label, keys, index_unique, fill_value, pairwise)
    870 X = concat_Xs(adatas, reindexers, axis=axis, fill_value=fill_value)
    872 if join == "inner":
--> 873     layers = inner_concat_aligned_mapping(
    874         [a.layers for a in adatas], axis=axis, reindexers=reindexers
    875     )
    876     concat_mapping = inner_concat_aligned_mapping(
    877         [getattr(a, f"{dim}m") for a in adatas], index=concat_indices
    878     )
    879     if pairwise:

File ~/miniconda3/lib/python3.8/site-packages/anndata/_core/merge.py:459, in inner_concat_aligned_mapping(mappings, reindexers, index, axis)
    456 def inner_concat_aligned_mapping(mappings, reindexers=None, index=None, axis=0):
    457     result = {}
--> 459     for k in intersect_keys(mappings):
    460         els = [m[k] for m in mappings]
    461         if reindexers is None:

File ~/miniconda3/lib/python3.8/site-packages/anndata/_core/merge.py:68, in intersect_keys(ds)
     67 def intersect_keys(ds: Collection) -> OrderedSet:
---> 68     return reduce(and_, map(OrderedSet, ds))

File ~/miniconda3/lib/python3.8/site-packages/anndata/_core/merge.py:32, in OrderedSet.__init__(self, vals)
     31 def __init__(self, vals=()):
---> 32     self.dict = OrderedDict(zip(vals, repeat(None)))

TypeError: 'NoneType' object is not iterable

Thanks


text:
How to apply a mathematical function to a matrix in R
I have a matrix that looks something like this:

chr
pos
het

chrI
27
2

chrI
55
0

chrI
27
0

chrI
55
1

It's essentially a VCF file that I have converted into a TSV. I want to be able to take all the values in the het column and apply this function:
het_df['het'] = het_df['het'].transform(lambda x: -1*(x-1))
That line is from a python script and it works really well in transforming that column. However, as I am making plots from this data I want to be able to do this in R. When I went to convert it to an R script I started with this function:
convert_df <- lapply(het_df[,3], FUN=function(x) formula=(-1*(x-1)))
As I saw it in this issue: https://stackoverflow.com/questions/7009333/how-to-use-lapply-with-a-formula
However, I keep getting error messages such as this one:
Error in x - 1 : non-numeric argument to binary operator

Which is odd because when I run het_df[,3] it gives me the correct output:
   [1] "2" "0" "1" "1" "1" "1" "1"

I think I'm just a little lost and I don't know what to do to run this function...everything worked fine in my python script (I was able to generate plots with matplotlib) but I like manipulating plots more in R than python so I wanted to be able to convert the script. If anyone is able to see what I'm clearly doing wrong please let me know! Thank you!
Here is the python code I've been using in case this helps (I'm transforming the last column so I always get a 1 or a 0 for heterozygosity/homosygosity and then using a rolling window to obtain the variability in the data series and plotting them):
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns

# read in file
het_df = pd.read_csv("het_file.tsv", sep="\t", header=None).set_axis(['chr', 'pos', 'het'], axis=1, copy=False)

# getting rid of homozygous alt allele
# always get 1 for heterozygosity
het_df['het'] = het_df['het'].transform(lambda x: -1*(x-1))

windows = [25, 50, 100, 200] # 100 bp, get rid of a lot of noise
steps = [5]

for window in windows:
    for step in steps:
        het_ratio = {}
        for chr_het_df in het_df.groupby('chr'):
            # rolling over the window
            # center = have the count in the center of the window
            # 50bp window count in the front or back of window or count in the center
            # ratio of heterozygosity
            het_ratio[chr_het_df[0]] = chr_het_df[1].set_index('pos').rolling(window=window, step=step, center=True).mean()
        fig, ax = plt.subplots(len(het_ratio), 1, sharex=True)
        for index in range(len(het_ratio)):
            chromosome = list(het_ratio.keys())[index]
            sns.scatterplot(het_ratio[chromosome], ax=ax[index], s=4)
            ax[index].set_ylim(bottom=0, top=1)
            ax[index].legend([], [], frameon=False)
            ax[index].set(xlabel="Chromosome position", ylabel=chromosome)
            ax[index].set_yticks([0, 0.5, 1])
            ax[index].set_yticklabels([0, 0.5, 1], size = 6)
        fig.subplots_adjust(top=0.98, bottom=0.05, left=0.05, right=0.98)
        fig.set_size_inches(20, 15)
        plt.savefig("het_plot" + str(window) + "_step" + str(step) + ".png", dpi=150)
```

If the numbers are in quotes it means that it has read them as text. A numeric vector will be:
[1] 2 0 1 1 1 1 1
To replace that column with the new value you can use this:
het_df[,3] <- sapply(as.numeric(het_df[,3]), FUN=function(x){-1*(x-1)})

I used het_df[, 3] <-  to replace the content of the column.
I changed lapply to sapply, to get a vector and not a list.
I forced the input to be numeric via as.numeric.
The syntax of the functions was wrong, to create a function in R you need function(x){#code}. The linked question is for a function that has an argument called formula, which is not the case here, because you want to define your own (anonymous) function.

Thank you so much that's exactly what I needed! That makes a lot more sense too - thank you for explaining what was going on.
If you wouldn't mind can I ask you one more question? How would you work a rolling window into this data? From my above python code I was able to complete it using a for loop and the pandas package. Everywhere I am looking to try and replicate it in R is telling me a lot of different things so if you know a good place to look or if you have any advice moving forward that would be great!

That's a different question entirely; please try to limit posts to one question only, and ask a new question if you have follow-ups.

Doing mathematical operations on a matrix will only work if all the matrix elements are numerical. In your case, you have mixed character and numerical variables, and it looks like a data frame or tibble would be more appropriate (calling it het_df is a strong indicator that this is actually what you want to do, as well as the pandas structures in the python code). You can load a TSV file into a tibble in R using the read_tsv function from tidyverse:
library(tidyverse);
het_df <- read_tsv("input_file_name.tsv");

After this, numerical mathematical operations can be carried out column-by-column, because the columns are treated as separate structures internally:
het_df$het <- -1 * (het_df$het - 1); # approach 1 [base R]
het_df <- mutate(het_df, het = -1 * (het - 1); # approach 2 [tidyverse]


text:
cromwell fails to download url input file
I have an input file that is a URL. Cromwell (v84) fails to pull it when I'm not using "Local" as a backend.  The only difference between a working case and failed case is the name of the backend in the cromwell config.
For example,
this works
default = "Local"
  providers
  {
    Local

And this fails
default = "MyLocal"
  providers
  {
    MyLocal

Command to reproduce error
export _JAVA_OPTIONS="--add-opens=java.base/sun.security.util=ALL-UNNAMED"

java -Dconfig.file=cromwell_docker.conf \
-Dbackend.providers.MyLocal.config.dockerRoot=$(pwd)/cromwell-executions \
-Dbackend.providers.MyLocal.config.root=$(pwd)/cromwell-executions \
-jar ~/cromwell/cromwell-84.jar run fq_count.wdl -i fq_count.json

Error:
java.lang.IllegalArgumentException: Could not build the path "https://portal.nersc.gov/cfs/m342/jaws/test_data/sample.fastq". It may refer to a filesystem not supported by this instance of Cromwell. Supported filesystems are: MacOSXFileSystem. Failures:

WDL file
version 1.0
workflow fq_count {
    input {
        File fastq_file
    }
    call count_reads {
      input: fastq = fastq_file
    }
}

task count_reads {
  input {
    File fastq
  }

  command <<<
    ls ~{fastq}
  >>>

  runtime {
    docker: "ubuntu:latest"
  }
}

json
{
  "fq_count.fastq_file": "https://portal.nersc.gov/cfs/m342/jaws/test_data/sample.fastq"
}

Cromwell conf
include required(classpath("application"))

webservice
{
  port = 50110
}

workflow-options
{
  workflow-log-dir: "cromwell-workflow-logs"
  workflow-log-temporary: false
}

call-caching
{
  enabled = true
  invalidate-bad-cache-result = true
}

backend
{
  default = "MyLocal"

  providers
  {

    MyLocal
    {
      actor-factory = "cromwell.backend.impl.sfs.config.ConfigBackendLifecycleActorFactory"

      config
      {

        concurrent-job-limit = 5
        run-in-background = true

        # The list of possible runtime custom attributes.
        runtime-attributes = """
          String? docker
          String? docker_user
        """

        # Submit string when there is no "docker" runtime attribute.
        submit = "/usr/bin/env bash ${script}"

        # Submit string when there is a "docker" runtime attribute.
        submit-docker = "docker run --volume=${cwd}:${docker_cwd} ${docker} /usr/bin/env bash ${script}"

        root = "/cromwell-executions"
        dockerRoot = "/cromwell-executions"

        filesystems
        {
          local
          {
            localization: [ "soft-link", "hard-link", "copy" ]
            http { enabled: true }

            caching {
                duplication-strategy: [ "soft-link", "hard-link", "copy" ]
                hashing-strategy: "xxh64"
            }
          }
        }

        default-runtime-attributes
        {
          failOnStderr: false
          continueOnReturnCode: 0
        }
      }
    }
  }
}


text:
how to download all metadata for a single TCGA dataset and link to data file UUIDs?
I've been able to download specific data types from the TCGA portal (https://portal.gdc.cancer.gov/) using the Manifest files and the GDC Tool.
But it's not clear to me how to download the metadata (e.g. cancer type, stage, etc.) with the links to the UUIDs of a specific data type, e.g. the UUID of a methylation data file to the case UUID:
https://portal.gdc.cancer.gov/files/14dc6f1f-21d0-4b05-9ad5-1d52f5020b12
Any ideas?

To download metadata associated with a specific data file, such as the cancer type and stage, you can use the GDC API (Application Programming Interface). The API allows you to query the GDC database and retrieve information in a structured format, such as JSON or XML.
To download the metadata for a specific data file, you will need the file's unique identifier, also known as its UUID. In your example, the UUID of the methylation data file is "14dc6f1f-21d0-4b05-9ad5-1d52f5020b12".
To use the GDC API, you can construct a query using the GDC API endpoint and the file's UUID. For example, to retrieve the metadata for the methylation data file with the UUID "14dc6f1f-21d0-4b05-9ad5-1d52f5020b12", you could use the following query:
https://api.gdc.cancer.gov/files/14dc6f1f-21d0-4b05-9ad5-1d52f5020b12?expand=metadata
This query uses the "files" endpoint of the GDC API and specifies the file's UUID. The "expand=metadata" parameter tells the API to return the metadata associated with the file. The response will be in json format:
{"data": {"data_format": "TXT", "access": "open", "file_name": "5b6353ee-af39-4bb4-a801-1454dfbc338b.methylation_array.sesame.level3betas.txt", "submitter_id": "d67719e6-4149-4492-a265-1f5d5689b904", "data_category": "DNA Methylation", "acl": ["open"], "type": "methylation_beta_value", "platform": "Illumina Human Methylation 450", "file_size": 13168876, "created_datetime": "2022-02-15T17:17:23.357435-06:00", "md5sum": "ab0c90db97a3e12634ef3f54f3ec835d", "updated_datetime": "2022-02-16T19:13:44.087132-06:00", "file_id": "14dc6f1f-21d0-4b05-9ad5-1d52f5020b12", "data_type": "Methylation Beta Value", "state": "released", "experimental_strategy": "Methylation Array", "version": "1", "data_release": "32.0 - 36.0"}, "warnings": {"expand": "unrecognized values: [metadata]"}}


text:
How many percents explained variance by the first 50 principal components of PCA analysis is suitable for downstream `KNN` and `umap`?
In single-cell RNA sequencing analysis, I am wondering how many percents explained variance by the first 50 principal components of PCA analysis is suitable for downstream KNN and umap?
If there are lots of samples (>40) and cells (>100,000) to analyze together, there are a lot of highly_variable_genes (>4000) can be found
import scanpy as sc
sc.pp.highly_variable_genes(adata, flavor='cell_ranger', n_top_genes=None)

The explained variance by the first 50 principal components of PCA analysis sc.pp.pca(adata, n_comps=50, use_highly_variable=True, svd_solver='arpack') is very low (<12%), increase the number of principal components does not help a lot.
I assume that the low explained variance of PCs will make KNN and umap difficult to separate the cells well. Am I correct?
sc.pp.neighbors(adata, n_neighbors=15)
sc.tl.umap(adata, min_dist=0.3, spread=2.0)

I can reduce the number of highly_variable_genes to get more explained variance by the first 50 principal components, but KNN and umap still can not separate the cells well.
I checked the first 5 PCs when I reduce the number of highly_variable_genes, and it doesn't look very high:
>>> np.cumsum(adata.uns['pca']['variance_ratio'])

array([0.08419679, 0.12635055, 0.13610176, 0.14490107, 0.15326007])

Any suggestions would be greatly appreciated.

Edit To answer the question, you're looking for 80% variance .... from the comments the actual command is np.cumsum which changes things a lot.
Previous response ...
Firstly it's still good work. It wouldn't be sensible to proceed with KNN if 12% explained 50 PCs, now there are 15% in 5 PCs.
You wrote,

0.08419679, 0.12635055, 0.13610176, 0.14490107, 0.15326007

PC1 ~ 0.085, PC2 ~ 0.068
I think that you've cross checked the analysis and that is a good sign of good analysis.
The only easy possibility is to completely tighten down on all "variable genes" parameter. The total variance for the analysis think 50% total variance rather than an idealistic 80%.
The emergency measures ...

Simply use say 10 samples
Reducing the total number of genes, e.g. removing more variable genes looks good, anything else
Check the input internally to make sure there is no formatting bug
Check via Seraut

The thing I would cross check is reducing the sample size a lot ... is there a stray sample(s) here that is throwing everything out? That would be thing I'd really look into. What happens to the total variance when say 20 samples, or 10 samples are used?
Just to finalise, are there any natural population divisions in the sample that could be used? If this results in a huge jump in variance that would be an answer. The only thing that is unusual is that there are a lot of samples for this type of study (from my perspective).

The array I showed is the accumulated sum of the first 5 PCs. Sorry for the confusion. Thanks.

Okay I missed it cumsum @DanLi


text:
Sanger sequencing annotation error
I am a student in a Cancer lab. Working with sanger is new to me. While analyzing a report we found an insertion that has not been reported in any databases so far, we were working on checking if the insertion is pathogenic or not but when we reanalyzed the sanger annotation through variant reporter of the insertion it was giving it one position ahead and we analyzed that with the VEP as well that is were the confusion started.
The sequencer showed annotation of an insertion not between the bases but a step/position next to it:
c.5532-5533insGAA is the insertion position came after position 32 and before 33 but the sequencer annotates it like c.5533-5534insGAA.
We are working on human genome hg38.
Should I consider it a machine error or should I take that position for the result? Or is there a nomenclature error I am making here?
And
the annotation is in this form [c.33-34insAUU]+[=] and I don't know what the +[=] represents here.

How do you know it should be c.32-33insAUU and not c.33-34insAUU? Could this be an issue of 0-based vs 1-based position counting?

I haven't heard of the 0-based and 1-based system. I read the linked post and it would make sense if that is the case. The insertion in that case would be 0-based, I assume. But it would be hard to say if that really is the case since I don't about this system.
And it would possibly affect the result analysis, right? bcz the 1st annotation in my data tells its frameshift insertion while the 2nd annotation shows protein termination even when there is no stop codon there, according to Ensembl VEP results.
(the annotation in the question are examples)

Please [edit] your question and add more detail. You said you were performing Sanger sequencing, but then mention VEP. VEP requires specific variants, so can only be used after sequencing. Please explain what you are really doing so we can understand. Basically, who is telling you c.32-33insAUU and who is telling you c.33-34insAUU? What input data is given in each case?

Please EDIT your question and explain what you have. Sanger sequencing doesn't give variants, it gives sequences. So presumably you are comparing to some reference genome and inferring the HGVS notation, right? Tell us. Then you are using VEP? OK that needs to be in the question. If possible, you should also give us the actual, full HGVS notation so we can look it up and see what the reference genome has in that position. I am guessing you are working with human? Which genome? hg19? hg38? Did VEP use the same one as you did originally?

Actually, I dont have full information on that right now. I'll have to look it up and be more specific in the details of the questions and let me edit it. Thank you

Try looking up the variant on a tool like https://varsome.com (disclaimer: I work for the company that builds that tool, but it's free) where you can see the underlying genome sequence, and can easily switch between hg19 and hg38. VarSome will also tell you if there are multiple, equivalent indels (various variants that all cause the same genotype, think of a deletion of A from AAA: no matter which A you delete, you end up with AA).

Thank you for your guidance, I'll redo my post once I gather more specific information.

Just to give my thoughts:

The missing information is the actual sequence

Is this a conflict in the change in position of the insert?
Could it simply be a consequence of down stream numbering, either a logical error (coding stuff), or lab artefact.

My bet is its point 2, but really we don't know except for the OP because we don't know the sequence.
Sources of error, Point 2

Sanger sequencing was traditionally known for "frame-shift" errors if there were mono-mers, e.g. AAAA -> AAAAA, it can mess up over the number of bases in a mono-mer.
PCR amplification can also cause slippage against monomers, same thing as the above point - but it's real.
1 vs 0 ... this is because in a Pythonic list or any C/Perl-based array numbering starts at 0 (not 1) - hence @terdon's comments.

If upstream mono-mers are present that might be the source of the confusion here and given the amount of sequence data thats very likely.
What I'd do:

Just Blast the new sequence and see what emerges: if it is genuinely novel - they'll be zero 100% identity hits.
VCFs ... I think this is VarSome's stuff (forget).

If there are Blast hits it will identify whether there is any upstream lab artefacts.
Conclusion
If this was any other situation in human genetics, or even mammalian genetics I'd be excited. However, in this case - I think its a routine result. Firstly, the chance of finding a new "driver" mutation are slim. Secondly, I've looked at this chromosome in a cancer patient (only once) and it was nuts (and thats an under statement) - I've a fairly good memory for numerical patterns, and I'd never seen a pathogenic species behaving like that. I think this insert is just a consequence of the entire mutational chaos that are cancer cells.
Just to make it clear: I don't do cancer but if you want an authoritative assessment and think you've something major I'd contact (or your PI) the Broad Institute via Geraldine A. Van der Auwera who also has an account on this site via @Geraldine_VdAuwera. Sanger sequencing is not what they do, but they'd know at the drop of a hat what the significance and if they don't respond - thats your answer ;-). Alternatively, VarSome would know, albeit they're commercial enterprise offering freeware.

@terdon
I have the reads generated by sanger. If there is a tool I can use to get the insertion coordinates to confirm? The sanger generated  results are doubtful because according to them the protein terminates due to insertion between a codon that would indicate severe condion which is not suggested in the patient's history.

The insertion is GAA so is it possible that the sequence shifted the results one base ahead? I have a picture of graph if you want to have alook

One thing you can try, @sanaamir, (and, again, I am pushing my own company's product here; I'm sorry, but those are the tools I am most familiar with) simply pasting the read into VarSome. That will align it against hg19 or hg38 (your choice), and identify any variants found in that read. That will help you see the full notation of the variant.

Thank you I will do that. But do i put the files from sanger or just fasta format of nucleotides with the insertion? I did have a look at the website today. @terdon

Please come into the chat room linked to the question (https://chat.stackexchange.com/rooms/141809/discussion-between-sana-amir-and-terdon) so we don't spam M_.


text:
How to change the position of groups colors in DimPlot of Seurat
I have a Seurat object and plotted the Dimplot for UMAP visualization for 2 variables, as shown in the image below.

Now, the problem is that I want the group by variables such as Non-responder and Responder and anti-CLTA4, anti-CLTA4+PD1, anti-PD1 on the top of the UMAP plot and not on the right side. My desired output would look like the following.

The timely help is highly appreciated, many thanks.

If you plan to export the graph anyway, the easiest way could be to edit it with a vector graphics program. The other option is to save the ggplot object returned by Dimplot() and modify it to suit your needs.

Could you please provide an example of Seurat DimPlot() to ggplot ? @Cloudberry

DimPlot returns a ggplot object, so ggplot functions can be applied to it. You can change the legend position like this:
DimPlot(data) + theme(legend.position = "top")
# or
DimPlot(data) + theme(legend.position = c(.1, .9))

Thank you so much, problem solved. I was trying to solve it from DimPlot() method.


text:
Genomic inflation factor distribution (lambda)
Considering that the p-value of a chi statistic(df=1) is uniform, the expected genomic inflation factor is λ=1. I show its calculation in R with the next chunk of code:
chi = qchisq(df$pval,df=1,lower.tail=FALSE)
factor = median(chi)/qchisq(p=0.5,df=1)

My question comes here. Imagine that I compute λ on several (e.g. >30 independent GWAS) studies. My null hypothesis is that all of them lack true-polygenic signals (i.e. uniform p-value distribution). What test can I do to reject H0? What is the distribution that I can assume for the inflation factors?

Why is p fixed at 0.5? Anyway it is unlikely to matter here: the answer could only be worked out from simulation studies in my opinion and a 0.05% threshold assigned against the resulting distribution.

You are right! I can simply do the simulation. On the other matter, I am not a statistician but I guess is related to the fact that 0.5 is the expected (p) value.

Given the distribution is bell-shaped (looks normal) then p=0.5 makes sense, that can be justified. The key is getting the right simulation to test your statistic to model what your statistic is testing.  If that test works well, it doesn't really matter how you got there. Obviously simulations need replications and that sort of thing.

The distribution as I suspected seems very much like a normal distribution. However if someone knows to explain the reason (limit central theorem?) I will give select their answer.
    # Do simulation
    factor_calc = function(pvals){
      chi = qchisq(pvals,df=1,lower.tail=FALSE)
      return(median(chi)/qchisq(p=0.5,df=1))}
    pvals_studies = sapply(1:10000, function(i){runif(10000)})
    factor_distribution = apply(pvals_studies, 2, function(pvals){factor_calc(pvals)})

    # Plot resulting distribution
    hist(factor_distribution)


text:
Merge the dataframes with python iterator
I have 3 experiment files (File1.tsv, File2.tsv, File3.tsv) from three experiments and 1 library file (library_file.tsv). I want to merge the files by creating a python iterator because the number of experiment files will increase with every additional experiment and the library file will be always the same. All the experiment files are in the same directory.
File 1
id  gene    Exp_1
ID000000001 EGFR    0.127
ID000000002 ADAD1   0.421
ID000000003 ADAD1   0.301

File2
id  gene    Exp2
ID000000001 EGFR    0.303
ID000000002 ADAD1   0.236
ID000000003 ADAD1   0.071

File3
id  gene    Exp_3
ID000000001 EGFR    0.847
ID000000002 ADAD1   0.992
ID000000003 ADAD1   0.912

library_file
id  value1  value2
ID000000001 Single  Guides_Single
ID000000002 Single  Guides_Single
ID000000003 Single  Guides_Single

Expected output
id  gene    Exp_1   Exp2    Exp_3   value1  value2
ID000000001 EGFR    0.127   0.303   0.847   Single  Guides_Single
ID000000002 ADAD1   0.421   0.236   0.992   Single  Guides_Single
ID000000003 ADAD1   0.301   0.071   0.912   Single  Guides_Single

Its a simple join operation.
import pandas as pd

df1 = pd.DataFrame({'id':['ID000000001','ID000000002','ID000000003'], 'gene1':['EGFR', 'ADAD1','ADAD1'],'score1':[0.127, 0.421, 0.301]}).set_index('id')
df2 = pd.DataFrame({'id':['ID000000001','ID000000002','ID000000003'], 'gene2':['EGFR', 'ADAD1','ADAD1'], 'score2':[0.303, 0.236, 0.071]}).set_index('id')
df3 = pd.DataFrame({'id':['ID000000001','ID000000002','ID000000003'], 'gene3':['EGFR', 'ADAD1','ADAD1'], 'score3':[0.847,0.992,0.912]}).set_index('id')
df4 = pd.DataFrame({'id':['ID000000001','ID000000002','ID000000003'], 'v1':['Single', 'Single', 'Single'], 'v2':['Guide', 'Guide', 'Guide']}).set_index('id')

df = df1.join([df2,df3,df4])
df.drop(['gene2', 'gene3'],axis=1,inplace=True)
print (df)

Output
             gene1  score1  score2  score3      v1     v2
id
ID000000001   EGFR   0.127   0.303   0.847  Single  Guide
ID000000002  ADAD1   0.421   0.236   0.992  Single  Guide
ID000000003  ADAD1   0.301   0.071   0.912  Single  Guide

'Iterator' just pickle via
path='/home/desktop/df.pkl'
df.to_pickle(path)
# Reload
df = pd.read_pickle(path)
# Add another data set
df = df.join(newdf)
df.drop(['gene4'],axis=1,inplace=True)

... then pickle it again

@Riya there are quite a few aspects to this question so maybe worth posting a new question. For example, pickle is one approach there are 3 others depending on what you want (readable, fast, long-term storage and compression). There is also the read csv and argparse issues. So you might consider refining the code towards a versatile application.


text:
Drawing synteny using CIRCOS from progressive MAUVE alignment
I am using Mauve app to align two whole genomes. I've aligned these genomes (gbk files) with the progressive Mauve aligner, and I got sp.sslist, sp.xmfa, sp.guide_tree,sp.backbone file, sp.bbcols etc file. I'm trying to generate a circos plot from MAUVE output. However, The output file format is different from the format required for CIRCOS. Here is the
output of sp.backbone:
seq0_leftend    seq0_rightend   seq1_leftend    seq1_rightend
-14446737   -14446957   89204   89424
-14446680   -14446734   89445   89499

but the CIRCOS input format is different,  as shown here:
chr - hs1 1 0 249250621 chr1
chr - hs2 2 0 243199373 chr2
chr - hs3 3 0 198022430 chr3

I want to draw synteny in a similar way as this example figure, which containg contigs of two chromosomes. Kindly help; I am new to python and all. A python script would be a great help for me. Thanks

What research have you done? Various tutorials exist, e.g. here. Your Mauve backbone file does not look like what that tutorial shows as expected input. I am not a Mauve person but I am pretty sure you want the alignment file output. https://darlinglab.org/mauve/user-guide/files.html


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.


text:
Improving prokaryotic assembly with other contig/scaffold-level data?
I have what at first sight appears to be a high-quality MAG (~10 pieces, high completion%) that I built from a hybrid assembly (Illumina + Nanopore data) from a cyanobacterium.
Workflow:
Quality control (BBDuk) > Hybrid assembly (Unicycler) > co-binning multiple different samples (vamb) > Quality control (CheckM, anvi'o)
*The Nanopore data is pretty old, and therefore I'd thought Unicycler might be a safer bet rather than long-read-first assemblers, due to the higher error rate
The same strain also has 1-2 very fragmented genomes (100s of pieces) available in NCBI at scaffold level. The reads themselves are not available, so I can't just add them to the assembly.
Is there a way to use the available NCBI genome(s) to try to check/polish mine? (I've checked out Improve a reference genome with sequencing data but it's already a few years old so I imagine different tools might be applicable)
Thank you for your time!

In bacteria de novo assemblies are the norm because they shift their genomes about a lot (which @Laura will know). A reference genome would be easy - even Spades would do this and that is likely a long way from the optimal approach. First question, which genome? Secondly, this is an impressive data set - long read, short read why not simply use advanced gap closing algorithms and the latest de novo assemblers?

Can you please explain how you built this genome in the first place? That will help to avoid answers that duplicate things you've already tried.

@M__ Thank you! Edited the question with more details, including workflow. I've used de novo assembly (never used reference-guided assembly before). But I can't help but wonder whether I can close some of the gaps still remaining by considering the NCBI genomes as well. Perhaps not - I imagine the gaps are characterized by repetitive sequences that would break any assembly...but maybe? One of the NCBI genomes is the same strain from the same culture collection that I used, so I thought there's a chance that rearrangements would not have happened. I'm not sure if this is a good idea.

@Laura I will put an answer together this week. Mine is old technology in this area for this technique. However, I personally don't think any reference genome approach will be an improvement, I think it will result in alternative gaps. I'll highlight an alternative modern de novo approach as well.

What is your input data? Is it isolate or metagenome? On one hand, you mentioned MAG (metagenome-assembled genome) and vamb. On the other hand, you are running Unicycler – I thought Unicycler wouldn't work well with metagenome samples? Also, how "old" is your Nanopore data? Flye/Canu can handle >10% error rate. If you have FAST5 files for the Nanopore data, redo base calling with the latest basecaller. This might dramatically improve the base accuracy. Your assembly seems to have much higher quality than the NCBI one. I would ignore the NCBI assembly.

@user172818 Isolate or metagenome - both? Low-diversity metagenome. For my field of work (cyanobacteria), the vast majority of isolates are not axenic. (in this case, the metagenome had 3 MAGs). I know Unicycler isn't technically recommended for those, but I wasn't sure what else to use that could handle older long-read data (~2018) as well as incorporate Illumina reads (which I don't think Flye/Canu do). No fast5 files anymore, only fastq. I'm probably missing something obvious, because I still have a lot of gaps in my knowledge. This is why I ask ; )

many metagenome assemblers can handle this kind of thing (old NP + ILMN), for example, hybridspades, OPERA-MS, probably others

If your goal is to compare to the NCBI version I'd suggest using QUAST to compare the two assemblies, supplying the NCBI assembly as the reference genome (-r option I think). This will not guarantee anything, but you would expect the two assemblies to be fairly syntenic and to have good overlap of sequences. The tool includes visualization, but you could also use an approach like circoletto or dot plots for something a little more intuitive/simpler.
I'd expect your genome to be quite a bit higher quality even as a metagenome-derived assembly due to nanopore contiguity, I'm not sure that using the NCBI reference to improve yours is a great idea. I'd instead suggest following @M__'s suggestion to scaffold with your nanopore reads (you will have to map them onto your full metagenome first, though approach may differ based on what the scaffolder wants). I'm not up on the current methods here but here is a recent tool. If you have 10 pieces I'd expect to be able to close at least some of that using nanopore (even old).
If you do choose to use the reference to improve, the tool that people seem to be using now is RagTag.

I've tried a couple of things that were suggested: in the end, LongStitch, which you've linked to, proved the best option, and it looks like it got the main genome to 3-4 pieces with some potential plasmids on the side (will have to check). Thank you!

@Laura great to know, glad that this seems to have been helpful!

Summary
What I am saying in this post is there approximates to 3 "generations" of in silico gap-closers/polishes:

GapFiller - amongst other "early"  generations (2012);
medaka/racon/Nanopolish/Pilon (2014 - 2018 inclusively)
Nextpolish/PEPPER/Apollo/Homopolish/NeuralPolish (2020 onwards).

The trendy approach is combining 2. and 3. I can't guarantee it will work for cyanobacteria gap-closing, but the shear volume of algorithms dedicated to this purpose, infers they must be delivering results.
Background
Whats the issue with a reference genome? The key thing about many  bacteria is their ability to "dump" genes and then regain them via plasmid associated recombination. This can occur extensively within a single species. This makes reference associated assembly essentially impossible in many bacteria - particularly against the official NCBI reference genome. Hence, the official NCBI reference gene, is not like a human genome reference - which is a gold standard. With due to respect to those who have dedicated time to presenting the official bacterial reference genomes are not very useful.
This BioStars post here demonstrates that even Mycobacteria tuberculosis would not assemble via a reference genome and TB in context is one of the most conserved genomes in bacteria, the SNP behaviour is famously limited.
If attempting a reference genome it is worth considering the following:

Any amount of genetic divergence can result genetic island phenomena and non-contiguous genomes that will make assembly difficult, therefore a complete genome as close to the target genome as possible is a good reference.

For example, simply Blasting sections of you genome to look for a consistent reference that is close to yours. A custom database of complete genomes for your species my be useful, but just sending bits of your target onto the NCBI website will be useful. The idea is to get lots of bits of the genome matching the closest strain.

QUAST was going to be my suggestion too.
Mapping via e.g. Mummer4 against the reference genome as the "scaffold".

Key question:

The same strain also has 1-2 very fragmented genomes (100s of pieces) available in NCBI at scaffold level. The reads themselves are not available, so I can't just add them to the assembly.

I'm not sure what "not available" means.
Key suggestion: Polishing
The other approach are recent de novo assemblers and in silico "polishing". What I'm recommending is to first try "polishing". This doesn't exclude subsequent reference genome assemblies - not at all - but it's a way to improve the de novo assemblies and if there are still issues then consider a reference genome.
Firstly I would read someone who is really up to date, e.g.

Comparative evaluation of Nanopore polishing tools for microbial genome assembly and polishing strategies for downstream analysis

They are medaka fans and are using this as part of a combination approach. medaka was part of the established stable of "polishers"

GapFiller - I don't think this is useful.
racon is specifically for long reads
Nanopolish here

The authors prefer medaka then goes on to describe

Nextpolish, PEPPER, Apollo, Homopolish, and NeuralPolish

Much more recent from 2020 ... thats a lot of polishing tools. Note, PEPPER and Apollo are often used for algorithms names - make sure its 2020 onwards, with PEPPER being 2021. They suggest medaka x PEPPER.
I really think this combination polishing approach with an established (medaka) and one the very latest polishing algorithms, such as PEPPER (2021) is really worth trying. If it doesn't work for the specific purpose it will certainly improve the contig quality. However, given the shear volume of work dedicated to this area then there has to be good reason.
By the same token that

Money doesnt give happiness, but sure gives a better standard of misery.

Then ...

Polishing might not give gap closure, but sure gives a better standard of quality to the same old bunch of miserable contigs

de novo Assembly
Maybe Miniasm? I'd look at "polishing" the de novo assembly you have.
RagTag as a SPAdes alternative is interesting @MaximilianPress

I'm not convinced long-read polishing is as useful for hybrid data, but I've given medaka and pilon (short-read polishing) a try. It didn't work in this particular case (at best moved some indels around, at worst it made more), but that's more because I've only tested things on my absolutely worst sample. Either way, it didn't work in this case, but now I know about these tools, so I wanted to say thank you!

@Laura one day the algorithms will catch up with the problem, they always do, but clearly not today. So, I'm really sorry it's not worked. Keep in mind the 2020+ tools. From a distance its seems an exciting future, where all de novo assemblies would be resolved via combination gap-filling and polishing. I was too optimistic ... for now. That means complex mapping - and a lot of work, via QUAST, Mummer and Daniel Standage's recommendations (his suggests looked cool).


text:
What is changed in this PDB file?
I used a tool called reduce to add Hydrogen bonds information to the PDB file.
I used the command:
$ ./reduce -NOFLIP 2GB1.pdb > 2GB1_with_H.pdb

The output file is supplied as follows:

https://textdoc.co/ZiuA0tY2elnoh5VK

I haven't understood the output.
What special change was made by the tool to the 2GB1.pdb file?

Did you read the paper and maybe contact the authors (http://kinemage.biochem.duke.edu/about/) if the paper doesn't make it clear?

Are you sure its pdt format and not pdbqt ? (Didnt check , but never found it so far) :https://userguide.mdanalysis.org/1.1.1/formats/reference/pdbqt.html

So the 'right' answer here would be for you to use a diff tool of some kind (WinDiff, or just the diff command in a Unix-like environment). That would show you the exact added lines. (Although, trying this myself does not produce very clear results...)
However, I know from running reduce on PDB files that it's added a few things:

The "USER MOD" lines, there's one on the first line and a large block just before the ATOM records
The hydrogen atoms, of course, needed for correct H-bond calculations

As some commenters have pointed out, this should really be in the paper or other documentation for the tool.


text:
Run nextflow pipeline multiple times depending on length of a list
This question was also asked on Stack Overflow
I have a Nextflow pipeline that runs a BUSCO search on input data, returns the output and then creates a report based on the output. One of my parameters is the lineage dataset. I would like to have my pipeline run multiple times (preferably at the same time) based on how long my list of lineage datasets is. For example, this is my nextflow.config file. The goal is that my pipeline would run twice in this case. If there are n lineages specified I would like it to run n times.
params {
    inputFilePath           = "$PWD/input_data/test_data.fasta"
    lineage_datasets        = [
        "fungi_odb10",
        "microsporidia_odb10"
        ]
    augustus_species        = "fumigatus"
    mode                    = "geno"

    threads                 = 16
    help                    = false
    outdir {
        main                = "$PWD/results"
    intermediate        = "$PWD/busco_outputs"
    }
}

The first part of my busco process looks like this:
  publishDir params.outdir.intermediate, mode: 'copy'

  input:
    path inputFile
    string lineage_dataset

  output:
    path 'hap1/*hap1.txt'

The idea is that each of the summary.txt files will be collected inside of my outdir.intermediate folder which is where my createReport process will target. For context, this is my createReport process:
process createReport {

  publishDir params.outdir.main, mode: 'copy'

  input:
    path 'busco_outputs/*'

  output:
    path 'report.html'

  script:
    """
    source "$launchDir/venv/bin/activate"
    report.py > report.html
    """
}

I have tried the following in my workflow process but it does not seem to work:
workflow {
  runBuscoChannel = Channel.fromPath(params.inputFilePath).cross(Channel.of(params.lineage_datasets))
  | mix
  | collect
  | set { outputs }

  createReport(outputs)
}

When I try to run the pipeline it does not run properly
Does anyone have any ideas of how I could change my file to allow for this? Thanks very much!

I think the answer to your question is the each input repeaters. You can also check Nextflow FAQ for the question on how to iterate over a process or a file.
Not in front of a PC now but you could try something like each lineage_dataset.


text:
Extract lines from a file based on a range of coordinates listed on a second file
I would be highly appreciated to have a script able to perform the following operation:
Only extract scaffold_1 lines from file1.vcf in which POS number (column 2) falls within column 2 and column 3 numbers of scaffold_1 lines from listfile.txt (inclusive). Similarly, only extract scaffold_2 lines from file1.vcf in which POS number falls within column 2 and column 3 numbers of scaffold_2 lines from listfile.txt. In reality the vcf I have continues until scaffold_12 and has 17,000 lines.
file1.vcf
#CHROM      POS ID  REF ALT QUAL    FILTER  INFO    FORMAT          COUNTS
scaffold_1  160 .   C   A   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_1  210 .   AC  A   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_1  340 .   A   C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_1  470 .   T   C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_1  580 .   A   AC  40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_1  630 .   CA  C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_2  40  .   C   A   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_2  100 .   AC  A   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_2  370 .   A   C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_2  600 .   T   C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_2  800 .   A   AC  40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_2      850     .       C       CTA     40      PASS    .       GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_2  1120    .   CA  C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40

listfile.txt
scaffold_1  50  100
scaffold_1  200 350
scaffold_2  350 450
scaffold_2  750 850

The desired output file is the following:
desired_file.vcf
#CHROM      POS ID  REF ALT QUAL    FILTER  INFO    FORMAT          COUNTS
scaffold_1  210 .   AC  A   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_1  340 .   A   C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_2  370 .   A   C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_2  800 .   A   AC  40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_2      850     .       C       CTA     40      PASS    .       GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40

To do this job I found the below awk code from another similar post but they don't differentiate scaffold_1 from scaffold_2, and I don't how to add this extra variable properly. Thank you.
awk 'NR==FNR{rng++;start[rng]=$2;end[rng]=$3;next} {for (i=1;i<=rng;i++) if (($2>=start[i])&&($2<=end[i])) {print; next}}' listfile.txt file1.vcf > output.vcf

<listfile.txt awk '{ print "$2 >= " $2 " && $2 <= " $3 }' | awk -f - file1.vcf > output.vcf

I tried to modify the last awk code like this, but it didn't work:
<listfile.txt awk '{ print "$1 == " $1 " && $2 >= " $2 " && $2 <= " $3 }' | awk -f - file1.vcf > output.vcf

Using GNU awk for arrays of arrays:
$ cat tst.awk
NR==FNR {
    ranges[$1][$2 FS $3]
    next
}
FNR == 1 {
    print
    next
}
$1 in ranges {
    for ( range in ranges[$1] ) {
        split(range,r)
        if ( (r[1] <= $2) && ($2 <= r[2]) ) {
            print
            next
        }
    }
}

$ awk -f tst.awk listfile.txt file1.vcf
#CHROM      POS ID  REF ALT QUAL    FILTER  INFO    FORMAT          COUNTS
scaffold_1  210 .   AC  A   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_1  340 .   A   C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_2  370 .   A   C   40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_2  800 .   A   AC  40  PASS    .   GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40
scaffold_2      850     .       C       CTA     40      PASS    .       GT:DP:AD:RO:QR:AO:QA    0:1:1,2:1:40:2:40

Using any awk it'd be:
$ cat tst.awk
NR==FNR {
    ranges[$1] = ranges[$1] $2 FS $3 RS
    next
}
FNR == 1 {
    print
    next
}
$1 in ranges {
    n = split(ranges[$1],list,RS)
    for ( i=1; i<n; i++ ) {
        range = list[i]
        split(range,r)
        if ( (r[1] <= $2) && ($2 <= r[2]) ) {
            print
            next
        }
    }
}

Hi Ed, your awk code works very well doing this job. Very clever solution. Thanks again for your timely help. That's right, the question I posted in the usual Unix page was moved to this site by the administrator.

@rseg whilst I am sure this is a great solution (ed definitely knows awk well!), I would strongly advise to use dedicated tools like bcftools. It is much easier to use, less prone to errors as they will have accounted for any of the peculiarities of vcf files, and is more robust in a pipeline. There's no need to use awk for such a simple example.

If you have a VCF file with an intact header, you can use bcftools for this. All you need to do first is to bgzip and index your VCF:
$ bgzip file.vcf
$ bcftools index file.vcf.gz

Then, assuming your list of regions are actually 1-based tab-delimited, use bcftools view and the --regions-file option to restrict the search to the regions listed. Use also the --no-header option to suppress the header in the output:
bcftools view -HR regions.txt file.vcf.gz

From the common options documentation:

-R, --regions-file FILE
Regions can be specified either on command line or in a VCF, BED, or tab-delimited file (the default). The
columns of the tab-delimited file can contain either positions
(two-column format) or intervals (three-column format): CHROM, POS,
and, optionally, END, where positions are 1-based and inclusive. The
columns of the tab-delimited BED file are also CHROM, POS and END
(trailing columns are ignored), but coordinates are 0-based,
half-open. To indicate that a file be treated as BED rather than the
1-based tab-delimited file, the file must have the ".bed" or ".bed.gz"
suffix (case-insensitive). Uncompressed files are stored in memory,
while bgzip-compressed and tabix-indexed region files are streamed.
Note that sequence names must match exactly, "chr20" is not the same
as "20". Also note that chromosome ordering in FILE will be respected,
the VCF will be processed in the order in which chromosomes first
appear in FILE. However, within chromosomes, the VCF will always be
processed in ascending genomic coordinate order no matter what order
they appear in FILE. Note that overlapping regions in FILE can result
in duplicated out of order positions in the output. This option
requires indexed VCF/BCF files. Note that -R cannot be used in
combination with -r.

Hi Steve, I converted my list file to tab delimited and following your suggestions bcftools worked great doing this job, and it's very fast too. It's handy to have a dedicated tool to parse vcf files. Thanks for your help.