# Question
Title: How to retrieve sequences from a txt file by protein names from other txt file?
My problem is that I want to retrieve the sequences line + secondary structure line that match the protein names from the `header.txt` file:

The text file(`header.txt`) contains the protein names that look like this:

```
1XXIJ
1P0EA
1H9HI

```

And the other file(`pdb.txt`) structured like this:

Line 1:the protein name

Line 2: protein length

Line 3: amino acid sequence

Line 4: secondary structure

Line 5: empty

```
4LGTD
247
M       S       E       K       L       Q       K       V       L       A       R
.       .       E       E       H       H       H       H       H       H       H

1XXIJ
334
M       R       W       Y       P       W       L       R       P       D       F
.       .       .       .       T       T       H       H       H       H       H

6U02G
116
E       V       Q       L       V       E       S       G       G       G       L
.       .       E       E       E       E       E       .       .       E       E

1P0EA
363
R       P       R       F       S       F       S       I       A       A       R
.       .       S       E       E       E       E       E       E       E       E

1H9HI
30
G       C       P       R       I       L       I       R       C       K       Q
.       E       E       .       .       .       .       B       S       S

```

So my expected output is:

```
1XXIJ
334
M       R       W       Y       P       W       L       R       P       D       F
.       .       .       .       T       T       H       H       H       H       H

1P0EA
363
R       P       R       F       S       F       S       I       A       A       R
.       .       S       E       E       E       E       E       E       E       E

1H9HI
30
G       C       P       R       I       L       I       R       C       K       Q
.       E       E       .       .       .       .       B       S       S

```

Do you guys have any suggestions for me by awk or any(Linux)code? Please let me know if you require additional information. Cheers!

# Answer
> 2 votes
If your file structure is consistent, one potential solution is to use GNU grep:

```
grep -A4 --no-group-separator -Fwf header.txt pdb.txt

```

This works with your example data and should be relatively fast on your 'real' data, due to `-Fw` ('fixed' i.e. not regex, and only match whole words) but it requires GNU grep for the `--no-group-separator` option (i.e. not BSD grep).

# Answer
> 2 votes
I'm assuming that your second file is tab-delimited and has no commas. The easiest way to do this would be to transpose each record in your second file to a single line and then grep, then transform the result again:

```
grep -Ff header.txt <(paste -d, - - - - - < pdb.txt) | tr "," "\n"

```

I just tested and it works:

<pre><code>➜ cat pdb.txt
4LGTD
247
M       S       E       K       L       Q       K       V       L       A       R
.       .       E       E       H       H       H       H       H       H       H

1XXIJ
334
M       R       W       Y       P       W       L       R       P       D       F
.       .       .       .       T       T       H       H       H       H       H

6U02G
116
E       V       Q       L       V       E       S       G       G       G       L
.       .       E       E       E       E       E       .       .       E       E

1P0EA
363
R       P       R       F       S       F       S       I       A       A       R
.       .       S       E       E       E       E       E       E       E       E

1H9HI
30
G       C       P       R       I       L       I       R       C       K       Q
.       E       E       .       .       .       .       B       S       S

➜ cat header.txt
1XXIJ
1P0EA
1H9HI

➜ grep -Ff header.txt <(paste -d, - - - - - < pdb.txt) | tr "," "\n"
1XXIJ
334
M       R       W       Y       P       W       L       R       P       D       F
.       .       .       .       T       T       H       H       H       H       H

1P0EA
363
R       P       R       F       S       F       S       I       A       A       R
.       .       S       E       E       E       E       E       E       E       E

1H9HI
30
G       C       P       R       I       L       I       R       C       K       Q
.       E       E       .       .       .       .       B       S       S
```
</code></pre>

# Answer
> 2 votes
One way using :

```
awk 'FNR==NR { a[$1]; next } FNR%5==1 { f = ($1 in a) ? 1 : 0 }f' header.txt pdb.txt

```

This adds each of the protein names in *header.txt* to an array. By adding the first field (`$1`), rather than the complete record (`$0`) to the array, we can ensure that we strip off any lagging whitespace which can sometimes cause problems depending on how the list of protein names was created. Next, we can check whether the record number (i.e. line number) of the second input file is a multiple of 5 with remainder 1. This condition returns true for the first line of each record in *pdb.txt*. We can then use a ternary operator to set a flag variable to one if the first field of the current line (i.e., the protein name) is present in the a array, *or* to zero to indicate that the current record is not one of the records in the array. The `f` by itself is just shorthand for `f != 0 { print }`, which prints the current line when the flag is non-zero.

---
Tags: proteins, protein-structure, linux, awk
---


# Question
Title: Gene symbol list for all protein coding genes in mice
How can I get a csv or a list of the gene symbol names for all the protein coding genes in mice? I have RNA sequencing data and I'm not interested in the non-coding stuff. I'm worried it could mess up my results. Can GEO or blast provide this?

Edit: I am not interested in EntrezID or Ensemble, I don't not want them I just am more interested in gene name like Fos or Penk. Etc. I work in RNA seq I don't need to map any fastq files to a genome that has been done already I just need to filter out genes that are not protien coding. If you are not familiar with what RNA seq data looks like it is simply rows of genes and columns of samples each entry giving the number of reads for that gene in that sample, single cell RNA seq is the same. Have a look at GEO if you want more examples: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE145970

# Answer
If you want *all* known mouse protein coding genes, you can get the list from Ensembl's BioMart.

First, select the mouse genome:

Then, in "Filters" limit to protein coding genes only:

And next, choose "Gene name" and, to have a more standardized and stable name and to ensure that you get allr esults (not all genes have a "name", but they will all have IDs), include the "Gene stable ID":

Finally, click on "Results":

Select "File" and "CSV" and then click on "Go" to download the list:

The end result is a file like this:

```
Gene stable ID,Gene name
ENSMUSG00000064341,mt-Nd1
ENSMUSG00000064345,mt-Nd2
ENSMUSG00000064351,mt-Co1
ENSMUSG00000064354,mt-Co2
ENSMUSG00000064356,mt-Atp8
ENSMUSG00000064357,mt-Atp6
ENSMUSG00000064358,mt-Co3
ENSMUSG00000064360,mt-Nd3
ENSMUSG00000065947,mt-Nd4l
ENSMUSG00000064363,mt-Nd4
ENSMUSG00000064367,mt-Nd5
ENSMUSG00000064368,mt-Nd6
ENSMUSG00000064370,mt-Cytb
ENSMUSG00000095019,
ENSMUSG00000094915,
ENSMUSG00000095041,
ENSMUSG00000095666,
ENSMUSG00000094741,Pramel39

```

With 21849 gene IDs, of which 21772 have a "name" as well.
> 5 votes

---
Tags: rna-seq, blast, rna, geoquery, geo
---


# Question
Title: 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/`).

# Answer
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 **mmCIF**s 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**
> 1 votes

# Answer
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 :

<pre><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

```
</code></pre>
> 1 votes

# Answer
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.
> 1 votes

# Answer
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']))

```
> 1 votes

---
Tags: python, pdb, data-download
---


# Question
Title: Plot phylogenetic tree from list of edges
I have a dataset that I wish to convert to tree or phylo format like in the `ape` package, in order to plot the phylogenetic tree. It is formatted like a list of edges and species names for the leaves. The dataset is quite large, but I have made an example dataset that represents it:

```
df <- data.frame(c(8,8,7,6,6,7,9,9,9),
                c(1,2,3,4,5,6,7,8,9),
                c(9, 9, 12, 14, 17, 19, 21, 23, 9),
                c("a","b","c","d","e", "node1", "node2", "node3", "node4"))
names(df) <- c("parent","node","branch.length", "label")

```

I wish to plot it as a rooted phylogenetic tree using the `ape` package, but I can't figure out how to convert it to for example a newick string. I have tried to look at solutions both in R and Python, but I can't figure it out. Does anyone have any idea how this can be done?

Additional information: I also have two other vectors, one contains the length of each edge, and one containing the age of the child node, equal to the sum of all edge length from that child to any leaf of the tree. It would be great if this could be included, but not neccessary for the plotting as it is only illustrative.

Thanks!

**Update**

I am using R version 4.1.2. One of my attempts at translating my data to Newick code (as i understand this could include the branch lengths), using one of \[these answers\]. I don't understand what goes wrong, but I have not worked with Newick code before so I am not very experienced with this. (https://biology.stackexchange.com/questions/37057/construct-phylogeny-form-edgelist):

```
df <- data.frame(c(8,8,7,6,6,7,9,9,9),
                c(1,2,3,4,5,6,7,8,9),
                c(9, 9, 12, 14, 17, 19, 21, 23, 9),
                c("a","b","c","d","e", "node1", "node2", "node3", "node4"))
names(df) <- c("parent","node","branch.length", "label")
tree_data <- df
str(tree_data)

tree_data$newick = paste(tree_data$label,':',signif(tree_data$branch.length,3),sep="")
nparents = 2
  while(nparents > 1){
  leaf = setdiff(tree_data$node,tree_data$parent)
  parent_leaf = unique(tree_data$parent[!is.na(match(tree_data$node,leaf))])
  index_remove = integer(0)
  for(i in 1:length(parent_leaf)){
parent_children = tree_data$node[which(tree_data$parent == parent_leaf[i])]
all_children_in_leaf = all(parent_children %in% leaf)
if(all_children_in_leaf==T){
  index_remove = c(index_remove,match(parent_children,tree_data$node))
  tree_data$newick[which(tree_data$node==parent_leaf[i])] = paste("(",paste(tree_data$newick[match(parent_children,tree_data$node)],collapse=" , "),")",tree_data$newick[which(tree_data$node==parent_leaf[i])] ,sep = " ")
}
  }
  tree_data = tree_data[-index_remove,]
  parent_next = unique(tree_data$parent[!is.na(match(tree_data$node,leaf))])
  nparents = length(parent_next)
  if (nparents==1) newick_code = paste("(",paste(tree_data$newick,collapse=","),")",colnames(tree_data)[parent_next],";",sep = " ")
}

tree <- ape::read.tree(text = newick_code)

plot.phylo(tree)

```

I don't get an error any more, but I get some unexpected output:

This is what I would expect, all nodes that are not single letters are inner nodes, not leaves:

I really appreciate any help!

# Answer
The following example uses python to generate the newick string. You can than import that in R for further analysis.

Please keep in mind that you can not use self reference for node4 (9,9,9), so I changed this to 'Root'.

```
from ete3 import Tree
import pandas as pd

df = pd.DataFrame({'parent': [8, 8, 7, 6, 6, 7, 9, 9, 'Root'],
                   'node': [1, 2, 3, 4, 5, 6, 7, 8, 9],
                   'branch.length': [9, 9, 12, 14, 17, 19, 21, 23, 9],
                   'label': ['a', 'b', 'c', 'd', 'e', 'node1', 'node2', 'node3', 'node4']})

parent_child_list = df[['parent', 'node', 'branch.length']].apply(tuple, axis=1).tolist()
tree = Tree.from_parent_child_table(parent_child_list)
newick_str = tree.write(format=1)

```

newick\_str:

```
'(((3:12,(4:14,5:17)6:19)7:21,(1:9,2:9)8:23)9:9);'

```
> 4 votes

---
Tags: r, python, phylogenetics, phylogeny
---


# Question
Title: How to tackle an Index Error in my code for finding ORFs
I have created a python code that prints out all lengths of ORFs in a DNA sequence (start codon - "ATG", stop codons - "TAG", "TAA" and "TGA". It seems to give correct answers, when compared with specialized ORF-finding programs. The thing is that besides correct answers it also gives an Index Error in the end:

```
Traceback (most recent call last):
  File "scripts.py", line 33, in <module>
    if list_start[i]>list_stop[n-1]:
IndexError: list index out of range

```

Now, I have created a dictionary from a FASTA file where *key* is a sequence name and *value* is the sequence itself and want to find ORFs for all sequences in that dictionary using a *for* loop and this code. However, the loop runs only once and I think that the Index Error to be the cause.

The code is the following:

```
#!/usr/bin/env python3
line="TCGGGCGAAGGCGGCAGCAAGTCGTCCACGCGCAGCGCGGCACCGCGGGCCTCTGCCGTGCGCTGCTTGGCCATGGCCTCCAGCGCACCGATCGGATCAAAGCCGCTGAAGCCTTCGCGCATCAGGCGGCCATAGTTGGCGCCAGTGACCGTACCAACCGCCTTGATGCGGCGCTCGGTCATCGCTGCATTGATCGAGTAGCCACCGCCGCCGCAAATGCCCAGCACGCCAATGCGTTCTTCATCCACATAGGGGAGCGTTACGAGGTAGTCGCAGACCACGCGGAAATCCTCGACGCGCAGTGTCGGGTCTTCGGTAAAACGTGGTTCGCCGCCGCTGGCACCCTGGAAGCTGGCGTCGAAGGCGATGACGACGAAACCTTCCTTGGCCAGCGCCTCGCCATACACGTTCCCCGATGTTTGCTCCTTGCAGCTGCCGATCGGATGCGCGCTGATGATGGCGGGATATTTCTTGCCTTCGTCGAAGTTCGGCGGGAAGTGGATGTCGGCTGCGATATCCCAATACACATTCTTGATCTTGACGCTTTTCATGACAGCTCCGTTCAGGGGGAGGGGGTAAGTTCGCCAGGCCGAATCGTTGGTAGCCAAGCGGCAACGACTCGAATATAGAGAGCCGATTGGAATTCCGTAAGATCGCAATCTGGACTACAGTGGTATCTTCAAATTGACAATGGCACCTACATGGATCCCTCACTGCTTCCGTCTCTCGCGTGGTTCGCCCACGTCGCACATCATCGTAGCTTCACGAAAGCGGCTGCGGAAATGGGCGTTTCTCGAGCAAACCTGTCGCAGAACGTGAAGGCGCTCGAACGCCGGTTGAACGTCAAGCTGCTGTATCGAACGACTCGCGACATGTCGCTGACCGAGGAGGGGCAGCGGCTCTACGAGGTGTGGTATCCCGCGCTGGTCGCGGTCGAGCGGACGGTCGACGCGCTGCACGAGGAGCGCGACGAGCCGTCGGGGTTGATTC"
r_frame=input('Enter a reading frame')
line=line[int(r_frame)-1:]
# making list of codons (3-nucleotide groups)
x=3
codon=[line[y-x:y] for y in range(x, len(line)+x,x)]
#defininig function  that will find all the indexes of codon in list of codons
def find_indices(list_to_check, item_to_find):
    indices = []
    for idx, value in enumerate(codon):
        if value == item_to_find:
            indices.append(idx)
    return indices
# makirng list with all the indexes of "ATG" (start codon)
list_start=find_indices(codon, "ATG")
#making list with all the indexes of "TGA", "TAA" or "TAG" (stop codons) and sorting it from smaller to bigger
list_stop=find_indices(codon, "TGA")+find_indices(codon, "TAA")+find_indices(codon, "TAG")
list_stop.sort()
i=0 # variable for indexes in list_start
n=0 #variable for indexes in list_stop
j=0 # variable that determines numer of times searching for ORF will be repeated
#first round without loop so that i and n will be >0 and we can make  if list_start[i]>list_stop[n-1] without errors
if list_start[i]<list_stop[n]:
    print(3*(list_stop[n]-list_start[i])+3)
    j=j+1
    n=n+1
    i=i+1
else:
    n=n+1
#loop that will search for suitable start codon-stop codon pairs and print the length of ORF
while j<codon.count('ATG'):
    if list_start[i]>list_stop[n-1]:
        if list_start[i]<list_stop[n]:
            print(3*(list_stop[n]-list_start[i])+3)
#multiply by 3 'cause one item in list stands for 3 nucleotides and add 3 so the whole length of stop codon will be counted
            j=j+1
            n=n+1
            i=i+1
        else:
            n=n+1
    else: i=i+1

```

For sequence given in the code the ORFs lengths should be as following:

Reading frame 1:63; 36; 36; 213; 129

Reading frame 2: 39

Reading frame 3: 186

This code might seem weird as I'm new to python and bioinformatics. I will be very thankful if anyone gives me an advice how can I get rid of the Index Error here.

# Answer
> 1 votes
Firstly, there is a bug in the code,

Original,

```
def find_indices(list_to_check, item_to_find):
    indices = []
    for idx, value in enumerate(codon):
        if value == item_to_find:
            indices.append(idx)
    return indices

```

Corrected,

```
def find_indices(list_to_check, item_to_find):
    indices = []
    for idx, value in enumerate(list_to_check):
        if value == item_to_find:
            indices.append(idx)
    return indices

```

You renamed the variable `codon` for the function, but then retained its name as `codon` within the function. It should be `list_to_check`. Also it is convenient.

If you are *certain* your code is fine, the easiest way to silence the error is via a `try` and `except` statement. The critical variable is `line` and I'm not sure where that is being declared. Without that its difficult to unpick the underlying error, which is an out of range error.

```
def find_indices(list_to_check, item_to_find):
        indices = []
        for idx, value in enumerate(list_to_check):
            if value == item_to_find:
                indices.append(idx)
        return indices

if __name__ == '__main__':
    #line= sample DNA nucleotide sequence

    try:
        r_frame= input('Enter a reading frame')
        line = line[int(r_frame)-1:]
        line=line[int(r_frame)-1:]
        x=3
        codon=[line[y-x:y] for y in range(x, len(line)+x,x)]
        list_start=find_indices(codon, "ATG")
        list_stop=find_indices(codon, "TGA")+find_indices(codon, "TAA")+find_indices(codon, "TAG")
        list_stop.sort()
        i=0 # variable for indexes in list_start
        n=0 #variable for indexes in list_stop
        j=0 # variable that determines numer of times searching for ORF will be repeated
        if list_start[i]<list_stop[n]:
            print(3*(list_stop[n]-list_start[i])+3)
            j=j+1
            n=n+1
            i=i+1
        else:
            n=n+1
        while j<codon.count('ATG'):
            if list_start[i]>list_stop[n-1]:
                if list_start[i]<list_stop[n]:
                    print(3*(list_stop[n]-list_start[i])+3)
                    j=j+1
                    n=n+1
                    i=i+1
                else:
                    n=n+1
            else:
                i = i + 1
    except AssertionError as error:
        print (error)
    except:
        pass

```

If you want to see the error you simply deleted the late two lines of code.

**However*** to resolve the second bug is quite complicated, particularly without a `line` declaration. We would simply use is `Biopython`'s `.translate()` method, we wouldn't do this from first principles. In a past question thats been done.

---
Tags: python, orf
---


title: Converting a Dante Lab VCF file

text:
# Question
Title: Converting a Dante Lab VCF file
Is there an easy way to grab the rsids (hg38) of a VCF file? I know one of the tabs can contain RSIDs, but this file doesn't contain any.

I have a software I created for mapping 23AndMe DTC testing to SNPedia, and I'm trying to integrate more DTC testing. (https://www.github.com/mentatpsi/osgenome)

The user found this https://gist.github.com/KlausTrainer/26c2996cf32677e4c107bd8aaee67794 which is a script to use variant summary data to map to the RSIDs, but I'm hesitant to perform a SQL dump to find a table that contains chromosome, postion, and rsID columns.

As a bonus question, 23AndMe reports all positive orientations. How would I go about converting strands to match this orientation? Does the VCF contain information that would allow this action autonomously?

# Answer
If your file doesn't contain the rsIDs, you can use a tool like Ensembl's VEP to annotate the file and get rsIDs as well. Of course, not all variants have an rsID so don't expect to find one for everything reported in your file.

More importantly, if these are the raw, unfiltered data from commercial direct-to-consumer sequencing products like 23&me, you need to be aware that they are incredibly untrustworthy. This is why 23&me don't include them in their official report. They only show you the quality results that have been filtered. When working with the raw data, you are very likely to have completely bogus results and will need to have a full VCF file and filter on various quality criteria.

I have looked at the raw 23&me data from 4 individuals, and *all* of them had extremely rare pathogenic variants reported that would basically have killed them long before they learned to walk. They were clearly false positives.

I stress that I am not suggesting these companies are producing crap: none of these variants was reported in the official report of 23&me, they were filtered out. The bad calls are only present when you work on the raw, unfiltered data but if that is what you are working on, you need to be extremely careful.
> 1 votes

---
Tags: vcf, gene, snp, genotyping, rsid
---


# Question
Title: How to run python request for list of url's with multiple page numbers?
Hi I am trying to get the cancer ontologies (obo\_id and label) from EBI-OLS. Earlier I have used the below code to get the `obo_id` terms and `labels`

```
descendants = requests.get(
    "https://www.ebi.ac.uk/ols/api/ontologies/efo/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FMONDO_0004992/descendants"
)

```

As there are more than 1 page so I created `total_pages`

```
total_pages = descendants.json()["page"]["totalElements"] // 1000 + 1

```

```
cancer_id = []

for i in range(total_pages):
    ## Check the number of id's present in 'totalElements'
    des = {"page": i, "size": descendants.json()["page"]["totalElements"]}
    ## Add the number of 'totalElements' and 'pages' in the new python request for descendants
    cancer_efo = requests.get(
        "https://www.ebi.ac.uk/ols/api/ontologies/efo/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FMONDO_0004992/descendants",
        params=des,
    )
    ## Get the number of obo_id's present in the cancer_efo.json
    obo = [obj["obo_id"] for obj in cancer_efo.json()["_embedded"]["terms"]]
    ## Get the number of cancer labels present in the cancer_efo.json
    label = [obj["label"] for obj in cancer_efo.json()["_embedded"]["terms"]]
    data = zip(obo, label)
    cancer_id.extend(data)

```

The above code worked pretty well with one link Now the issue is that I want to get ontologies for 37 descendent links and maximum links have more than 1 page. I want to get one file (cancer\_id) containing all obo\_id's and labels from these all 37 links.

I made a separate file called `list.txt` (Right now I am copying only 3 here)

`list.txt`

```
https://www.ebi.ac.uk/ols/api/ontologies/ncit/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FNCIT_C4861/descendants
https://www.ebi.ac.uk/ols/api/ontologies/ncit/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FNCIT_C3171/descendants
https://www.ebi.ac.uk/ols/api/ontologies/ncit/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FNCIT_C3457/descendants

```

I tried the below code but it is not working

```
import requests

with open ('list.txt', 'r') as fin:
    links = fin.readlines()

def requestIt(links):
    cancer_all = requests.get(links)
    return cancer_all.json()

data = list(map(requestIt, links))

obo = obj["obo_id"] for obj in data

```

The error is

```
obo = obj["obo_id"] for obj in data
                        ^
SyntaxError: invalid syntax

```

# Answer
> 1 votes
Okay, its done. What you are looking for is as follows:

```
import requests
from pathlib import Path

# Test code
def get_all_keys(d):
    for key, value in d.items():
        yield key
        if isinstance(value, dict):
            yield from get_all_keys(value)

def requestIt(links):
    cancer_all = requests.get(links)
    return cancer_all.json()

if __name__ == '__main__':

    path = Path('/PathToDir/test', 'link.html')
    with open (path, 'r') as fin:
        links = fin.readlines()

    data = list(map(requestIt, links))

    # Test code
    for x in get_all_keys(data[0]['_embedded']['terms'][0]):
        print(x)

    # Answer code
    for i,d in enumerate(data):
        print (f'Link {i + 1} ')
        for obo_id in d['_embedded']['terms']:
            print (obo_id['obo_id'])

```

**Output**

> **Link 1** NCIT:C8263 NCIT:C9163 NCIT:C5630

> **Link 2** NCIT:C7175 NCIT:C68700 NCIT:C68696 NCIT:C68697 NCIT:C68698 NCIT:C68699 NCIT:C9155 NCIT:C9288 NCIT:C9287 NCIT:C9019 NCIT:C9018 NCIT:C36055 NCIT:C36058 NCIT:C36057 NCIT:C36056 NCIT:C38377 NCIT:C122688 NCIT:C156718 NCIT:C156719 NCIT:C169107

> **Link 3** NCIT:C7056 NCIT:C131911 NCIT:C165799 NCIT:C68684 NCIT:C4998 NCIT:C5089 NCIT:C5095 NCIT:C8154 NCIT:C8604 NCIT:C127840 NCIT:C8489 NCIT:C8852 NCIT:C69144 NCIT:C115349 NCIT:C115367 NCIT:C165240 NCIT:C8926 NCIT:C160230 NCIT:C160231 NCIT:C182035

Where **Link 1** is ...

```
https://www.ebi.ac.uk/ols/api/ontologies/ncit/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FNCIT_C4861/descendants

```

... and so on.

If the code is listed `# Test code` that is what is used to interrogate the data structure. This doesn't contribute to your answer.

However, given the way I've done it, it's important to understand the data structure. Thus output of `# Test code` is at the base of this post. If you need another part of the data structure you simply switch `id_odo` to whatever else you are looking for (at the base of this post).

There will be an easy to use OOP method based out the `data` object. I did't want to spend time investigating the object ... the question is specifically about `obo_id` so it was easier to go for that.

**Test code output**

```
iri
lang
description
synonyms
annotation
Contributing_Source
Display_Name
Legacy Concept Name
Neoplastic_Status
Preferred_Name
Semantic_Type
UMLS_CUI
code
label
ontology_name
ontology_prefix
ontology_iri
is_obsolete
term_replaced_by
is_defining_ontology
has_children
is_root
short_form
obo_id  # <- its here
in_subset
obo_definition_citation
obo_xref
obo_synonym
is_preferred_root
_links
self
href
parents
href
ancestors
href
hierarchicalParents
href
hierarchicalAncestors
href
jstree
href
graph
href

```

# Answer
> 2 votes
@Megha The requestIt(link) function returns a list of dictionary objects. So your "data" variable is a list of dictionary objects.

Every item on this list (a dict) has a key named "\_embedded" whose value is also a dictionary. This dictionary (the \_embedded) has a key named "terms" which is a list of dictionaries. The obo\_ids you need are inside the dictionaries of the "terms" list.

The problem with your code was that you were trying to access the "obo\_id" key on the top level dictionary.

So a simple solution would be like this:

```
with open ('urls.txt', 'r') as fin:
    links = fin.readlines()
resps = [requestIt(link) for link in links] # this one holds all the data
obo_ids = []
for resp in resps:                          # for every response we got
    terms = resp['_embedded']['terms']      # choose the terms list
    for term in terms:                      # for every dict in this list
        obo_ids.append(term['obo_id'])      # extract obo_id and append it to our answer
print(obo_ids)                              # do stuff with obo_ids

```

You should probably perform some exception handling inside the requestIt() function to handle broken links etc.

# Answer
> 1 votes
@Megha, I have not looked through the code in too much detail, however.

```
obo = obj["obo_id"] for obj in data

```

This isn't formal Python code. The correction would simply be to place this in a comprehension ...

```
obo = [obj["obo_id"] for obj in data]

```

... unless you want an iterator.

I think there is a deep issue with `obj["obo_id"]`. If the syntax still throws an exception could you provide a print out of either `data` please, or else ....

```
obj for obj in data:
    print (obj)
    print (obj.keys)

```

The `data` data structure for the output of `get.requests` followed by `.json()` is kinda key here

---
Tags: python, api, pandas, ontology
---


# Question
Title: Identifying somatic mutations in cell lines
I would like to identify the somatic mutations present in a cell line and characterise the genes that are potentially affected by those mutations. For example, are there oncogenes mutated in a subpopulation in the cell line?

I have currently access to single cell WGS data and sequencing data of the bulk sample derived from the same clone. The data was generate after MDA amplification of the single cell's genome so I am currently evaluating SCCaller and Prosolo for variant calling. If I correctly understand single cell variant calling, mutations called in each cell are called relatively to the bulk background sample. So each cell separately doesn't provide much information on the cell line itself.

My questions are:

* Do i need single cell WGS to quantify the overall 'genomic health' of a cell line or can I just use deep sequencing followed by somatic variant calling against a reference genome.
* What added value can single cell WGS provide here and how can I extract that information?

Links to papers, packages or statistical models more than welcome.

Thanks!

# Answer
> 0 votes
I am not actively following this field, so take my answer with a grain of salt.

With bulk data, you can't easily tell how subclones differ from each other or characterize the mutation rates in subclones. It is also challenging to get accurate somatic mutation calls especially when you don't have the matched normal.

With single-cell data, you usually get more accurate point mutation calls in 20-80% of the genome, depending on the amplification method. You can estimate the mutation rate and call large CNVs more accurately in each cell. Sometimes you can even reconstruct the evolution of subclones. I would say single-cell data is better than bulk data all around. The only problem is that the single-cell approach is more costly as you need to sequence many cells to see the pattern.

On a side note, MDA is not that accurate for calling point mutations. This may not be a problem if the mutation rate is much higher than the error rate.

---
Tags: single-cell, genomics, variant-calling, mutations, variants
---


# Question
Title: To find genes that don't change in RNA-seq, Deseq2 has altHypothesis="lessAbs". Is there a way to make limma do the same thing for proteomics?
I love that Deseq2 has altHypothesis="lessAbs" !!!!!!! I've used it a ton for RNA-seq. However, now I'm working with mass spect data using limma. Is there a way to make limma do the same thing?

# Answer
> 2 votes
There is no direct support for this by best knowledge. Similar questions have been asked previously in the Bioconductor support forum, for example:

https://support.bioconductor.org/p/64300/#64793

They discuss options around the TOST concept.

For edgeR:

https://support.bioconductor.org/p/66283/#66287

From these posts please go back and forth to find other posts in that forum to capture what the developers of limma/edgeR responded on the matter.

Depending on your exact purpose it might even be enough to go for genes with large AveExpr and a nominal (not adj) p-value close to one.

---
Tags: r, rna-seq, differential-expression, proteomics, limma
---


title: file order after nextflow groupTuple

text:
# Question
Title: file order after nextflow groupTuple
Below is the code from nextflow site (https://www.nextflow.io/docs/edge/operator.html#grouptuple).

For the output, the first is A B C, and the second is C, A.

Is there a way to make the output like this

\[1, \[A, B, C\]\] \[2, \[A, C\]\]

That is, the stuff inside \[\] is in some order.

```
Channel
    .of(
        [1, 'A'],
        [1, 'B'],
        [2, 'C'],
        [3, 'B'],
        [1, 'C'],
        [2, 'A'],
        [3, 'D']
    )
    .groupTuple()
    .view()

[1, [A, B, C]]
[2, [C, A]]
[3, [B, D]]

```

# Answer
> 2 votes
You can change the sorting criteria for the grouped items using the `sort` keyword. To perform a *natural* sort, just use `sort: true`. For example:

```
Channel.of(
    [1, 'A'],
    [1, 'B'],
    [2, 'C'],
    [3, 'B'],
    [1, 'C'],
    [2, 'A'],
    [3, 'D']
)
.groupTuple(sort: true)
.view()

```

Results:

```
$ nextflow run main.nf
N E X T F L O W  ~  version 22.04.0
Launching `main.nf` [special_darwin] DSL2 - revision: 2e641ccfd9
[1, [A, B, C]]
[2, [A, C]]
[3, [B, D]]

```

Otherwise, a custom sorting criteria can be specified using a closure.

---
Tags: nextflow
---


# Question
Title: Replace a string in a column of a text file based on matching string listed in another file using linux
I have multiple sets of files, so I grepped the list from a column in a seperate file such as file1.txt. I want to change all other files by using this file1.txt based on a matching string.

for example I am pasting three files here

**File 1:**

```
chr - seq1 NZ_JAHWGH010000010.1 0 60562 green_a4
chr - seq3 NZ_JAHWGH010000012.1 0 466573 green_a4
chr - seq5 NZ_JAHWGH010000013.1 0 125526 green_a4
chr - seq39 NZ_CAJGBF010000006.1 0 211154 green_a4
chr - seq40 NZ_CAJGBF010000007.1 0 199119 green_a4
chr - seq41 NZ_CAJGBF010000008.1 0 190934 green_a4
chr - seq42 NZ_CAJGBF010000009.1 0 178251 green_a4

```

I extracted the required information using grep; the command is given below:

```
grep "NZ_JAHWGH" file.txt |awk '{print $3}' >file1.txt
grep "NZ_CAJGBF" file.txt |awk '{print $3}' >file2.txt

```

Now I wanted to use this file1.txt and file2.txt to change another files. Output will be like if list in file1.txt matches then all matched lists should be substituted as seq1 similarly using the second file2.txt I want to substitute all matching patterns as seq2.

For example the file is given here to be changed by using information from the previous file:

```
seq3 19301 20914 fill_color=green_a0
seq1 16726 18337 fill_color=green_a0
seq5 167934 169650 fill_color=green_a0
seq42 200621 202367 fill_color=red_a0
seq41 160164 161903 fill_color=green_a0
seq40 31356 33104 fill_color=green_a0

```

And final output after the change should be as follows:

```
seq1 19301 20914 fill_color=green_a0
seq1 16726 18337 fill_color=green_a0
seq1 167934 169650 fill_color=green_a0
seq2 200621 202367 fill_color=red_a0
seq2 160164 161903 fill_color=green_a0
seq2 31356 33104 fill_color=green_a0

```

I tried using this command:

```
for i in f2.txt; do grep -wf f1.txt $i|awk '{$1="seq1";print}'; grep -wf other.txt $i |awk '{$1="seq2";print}'; done

```

2nd attempt:

```
grep -wf f1.txt f2.txt|sed -e 's/seq[0-9]*/seq1/g'

```

Is this correct? If it is correct is there any more efficient way to do it? also it is only printing to the console. I want the output and also changes in the main file. Kindly improve my scripts.

I appreciate yout time and help.

# Answer
> 0 votes
Here's one way using . This avoids the creation of the intermediary files by matching the alphabetic characters in the fourth column. If the prefix has not been seen, it is added to an array along with a counter, which is only incremented when the prefix has not yet been seen. We then store the seq ids and the replacement strings in another array. Once we've finished reading the first file, we can then query the first column of the second file in the arguments list. If the seq id in the first column exists in the array of seq ids, we can print the replacement string and remaining columns.

With the following in a file called `script.awk`:

```
BEGIN {

   counter = 1
}

FNR==NR {
    if( match($4, /^NZ_[A-Z]+/) ) {
        prefix = substr($4, RSTART, RLENGTH)

        if( !(prefix in seen) ) {
            seen[prefix] = counter++
        }

        seq_ids[$3] = "seq" seen[prefix]
    }

    next
}

$1 in seq_ids {

    print seq_ids[$1], $2, $3, $4
}

```

Run using:

```
awk -f script.awk file1 file2

```

Results:

```
seq1 19301 20914 fill_color=green_a0
seq1 16726 18337 fill_color=green_a0
seq1 167934 169650 fill_color=green_a0
seq2 200621 202367 fill_color=red_a0
seq2 160164 161903 fill_color=green_a0
seq2 31356 33104 fill_color=green_a0

```

---
Tags: bash, linux, awk, sed, for-loop
---


# Question
Title: Given a FASTA file with a set of genes, how do I determine that set's conservation among a genetic order?
Given a set of genes in a FASTA-formatted file, for example: `exFasta.fa`, what is the best way to determine how well conserved the genes in that file are among an order of organisms, for example: Orthoptera?

BLAST is an obvious option, but I am not sure that it is the most efficient way considering there are likely thousands of organisms in that order. OrthoFinder also comes across as a valid option, but I am not sure how to implement such a thing considering the required inputs.

Any advice would be greatly appreciated!

# Answer
> 2 votes
It's possible that MOLE-BLAST could help for you:

> MOLE-BLAST is an experimental tool that helps taxonomists find closest database neighbors of submitted query sequences. It computes a multiple sequence alignment (MSA) between the query sequences along with their top BLAST database hits, and generates a phylogenetic tree. Query sequences in the tree are denoted with highlighted node labels. If the input sequences come from different genes or loci, MOLE-BLAST can cluster them and compute an MSA and a phylogenetic tree for each locus separately. This tool may eventually be incorporated into our supported suite of BLAST applications based on the feedback we receive.

I haven't actually used MOLE-BLAST, but I have used Smart BLAST, which does a somewhat similar phylogenetic search (although on a restricted diverse taxonomic set) for protein queries.

---
Tags: python, blast, blastp, orthofinder
---


title: Is my reference sequence too small?

text:
# Question
Title: Is my reference sequence too small?
I'm trying to map ONT long reads to a portion of a gene I'm looking at. The region is about 25bp long. When I search for the region in the document it pulls up the sequence in every read but when I map to reference with minimap2 it produces no contigs.

Is 25bp too small to be a reference? I can potentially extend it but it's right in an area of conservation in this region of a gene and I would rather not extend much farther.

# Answer
> 3 votes
That's not the right way to do it.

You either have reference or not. If not, you need to assemble one (that problem is called genome assembly). If yes, you map all the reads to the whole reference - the process of mapping is made for this - placing individual WGS reads to their corresponding locations and as a matter of fact, missing portions of the reference can mess it up (that was one of the reasons why there was a push for T2T for human).

Once the reads are mapped, you can use any genome browser (e.g. IGV) to zoom in to the exact location you are interested in and check the stack of reads mapping in there. Given how small the region is, manual inspection of the quality of the mapping is the best you got.

---
Tags: sequence-alignment, gene, nanopore, minimap2
---


# Question
Title: Using R `ggtree`, attach extra data next to the nodes on an illustrated tree/newick file
This is really a programming question but `ggtree` might be too obscure for this question to be posted on StackOverFlow as a general R related question.

I am having trouble extrapolating this this tutorial and the `%<+%` operator in `ggtree`.

I have a Newick file and a table that has a number for every tip and node in the Newick file. I want to produce a final graphic that shows a tree where the numbers from the table are shown next to respective tips and nodes.

What I want to do is close to the tutorial, I want to take a tree object and use a data-frame where one column is node-names and the second column holds numerical information per node. I want to attach the information from the data-frame to the tree object.

The data-frame looks something like this:

| Node / tips | Count |
| --- | --- |
| A | 12 |
| B | 14 |

So, I want to ultimately display the number 12 next to node A and number 14 next to node B while keeping the names intact. Now I have the code working up to the following point:

(I am really only just working my way through learning R's grammar of graphics and R in general, so assume I am doing something wrong.)

```
library(ggplot2)
library(ggtree)

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

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

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

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

print(p2)

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

```

The above code overlays the `Newick_label` column with the node names instead of attaching the value in the second column.

But specifically, I am at a loss as to how I can list the options available within `geom_tippoint(aes())`. I can use `?geom_tippoint` to access the documentation of `geom_tippoint` but I think my task needs to manipulate options within the `aes()` function and the the geom\_tippoint() documentation does not list all the `aes()` options available. And am stuck trying different keywords for `geom_tippoint(aes())`.

Can that be done in `ggtree` or is this specific application of `ggtree` obscure enough that there is no way to do this and produce a legible graphic? (eg issue: the graphic getting mushy overlays). I have hundred of (not uniform) genes trees so splitting data for each layer (node and tips) might be not the best way.

The files if anyone wants to try: counts and tree

Clarification on the labels for nodes and tips, I want to label both if possible. The data contains information for both.

Is annotating tips AND nodes out of the ordinary? Tips and nodes are technically different but I have quantitative data at both levels and am looking for the option to annotate at both levels.

# Answer
> 1 votes
I ended up using the following code to solve this issue:

Change the `tsv` file and the `tree` file name accordingly. Be mindful of the format the `tsv` file should be in (indicated in the code).

Also change the tree and the font scales as needed.

# Load required libraries
library(ape)
library(ggtree)

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

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

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

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

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

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

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

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

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

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

---
Tags: r, phylogenetics, ggtree
---


# Question
Title: How to programatically download SARS-CoV-2 fasta from NCBI/Genbank via API?
I would like to download a number of SARS-CoV-2 fastas from NCBI/Genbank using a web endpoint. I looked on their website but it's confusing.

I have the accession numbers, but can't figure out how to use the NCBI/Genbank Entrez API.

This is how I would do it using cov-spectrum: https://lapis.cov-spectrum.org/open/v1/sample/fasta?genbankAccession=MZ039713,MZ039682,OV522182

Is there something similar for NCBI/Genbank?

There must be a way, I just can't figure it out. The whole entrez thing is overwhelming to me, all I want is a simple API, like covSpectrum offers, for example. I scanned this but didn't help: https://www.ncbi.nlm.nih.gov/sars-cov-2/

E.g. how could I download all these SARS-CoV-2 nucleotide sequences with accession numbers: `MZ039713`,`MZ039682` and `OV522182` using a URL?

Some command line tools are advertised, but I don't want to have to use a CLI tool. Just a simple URL, does that not exist? CLI tools: https://www.ncbi.nlm.nih.gov/datasets/docs/v1/quickstarts/command-line-tools/

I also don't want to do any clicking through a web interface, just a URL so that other web apps can download it behind the scenes.

# Answer
You can get the data from the API of the European Nucleotide Archive (ENA) which is synched with NCBI/GenBank:

https://www.ebi.ac.uk/ena/browser/api/fasta/MZ039713,MZ039682,OV522182
> 2 votes

# Answer
You can try the NCBI Datasets API. We have an endpoint that would allow you to download virus sequences by accession The only thing to be aware (based on your list of accessions) is that our API needs the accession version to return results. In the case of your examples, they would need to be MZ039713.1, MZ039682.1 and OV522182.1. I hope this helps!
> 2 votes

# Answer
In the `efetch.fcgi` URL, you can set `rettype=fasta` to download the nucleotide sequence or `rettype=fasta_cds_aa` to download protein coding sequences:

```
$ curl 'https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&rettype=fasta&id=NC_045512.2'|head -n3
>NC_045512.2 Severe acute respiratory syndrome coronavirus 2 isolate Wuhan-Hu-1, complete genome
ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATCTGTTCTCTAAA
CGAACTTTAAAATCTGTGTGGCTGTCACTCGGCTGCATGCTTAGTGCACTCACGCAGTATAATTAATAAC
$ curl 'https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&rettype=fasta_cds_aa&id=NC_045512.2'|head -n3
>lcl|NC_045512.2_prot_YP_009724389.1_1 [gene=ORF1ab] [locus_tag=GU280_gp01] [db_xref=GeneID:43740578] [protein=ORF1ab polyprotein] [exception=ribosomal slippage] [protein_id=YP_009724389.1] [location=join(266..13468,13468..21555)] [gbkey=CDS]
MESLVPGFNEKTHVQLSLPVLQVRDVLVRGFGDSVEEVLSEARQHLKDGTCGLVEVEKGVLPQLEQPYVF
IKRSDARTAPHGHVMVELVAELEGIQYGRSGETLGVLVPHVGEIPVAYRKVLLRKNGNKGAGGHSYGADL

```
> 2 votes

# Answer
There are three ways to do this,

1. If you click on this link Your NCBI examples. Its done.
2. You would then just click "download"
3. Click "next", click "next"

The sequence then downloads.

Essentially all I did was in the website click on the left hand navigation bar and clicked "accession number" and type in `MZ039713 MZ039682 OV522182`.

1. The second way is to go to https://ncbi.nlm.nih.gov/nuccore
2. Type type in `MZ039713 MZ039682 OV522182`
3. Click on the "Summary" and click on "Fasta"
4. Click on "send to" and click on "File"

I'm not certain is point 3 immediately above is required.

The final approach is to use efetch to automate the download, but that requires downloading the NCBI's efetch suite, which is then flipped through on a bash loop.

```
efetch -db nuccore -id MZ039682 -format fasta

```

The link for looking to automate is: https://www.ncbi.nlm.nih.gov/labs/virus/vssi/#/virus?SeqType\_s=Nucleotide&VirusLineage\_ss=SARS-CoV-2,%20taxid:2697049&ids=MZ039713%20MZ039682%20OV522182%20

Just to highlight **MZ039713**%20**MZ039682**%20**OV522182%20**

The API is here, but I don't use them.
> 1 votes

---
Tags: fasta, phylogenetics, sars-cov-2, covid-19, genbank
---


title: Nextflow Pipeline: Unexpected Input "{"

text:
# Question
Title: Nextflow Pipeline: Unexpected Input "{"
I'm attempting to run a nextflow pipeline and have begun getting unexpected output for a few of the processes in my pipeline. I believe this to be a syntax error, but I am not sure where to correct it. Any help is finding where my error is located is appreciated.

```
#!/usr/bin/env nextflow
 nextflow.enable.dsl=2
 import groovy.io.*

  process hisat2Alignment {
    fastq1 = NS.1389.002.NEBNext_dual_i7_*_R1.fastq.gz
    fastq2 = NS.1389.002.NEBNext_dual_i7_*_R2.fastq.gz
    input:
    tuple path("fastq1"), path("fastq2")

    output:
    tuple path("*.sam"), val(fastq1.baseName.substring(0, fastq1.baseName.indexOf("_R1")))
    script:
    """
    hisat2 -p 16 -x genome -1 ${fastq1} -2 ${fastq2} -S mapped_${baseName}.sam
    """

 }

process Stringtier {
    publishDir "${params.baseDir}/assembled/$baseName", mode: "move"
    env {
     PATH = $PATH:/data/selena/stringtie-2.2.1
    }

    input:
    tuple path(mapped_bam), val(baseName)
    output:
    tuple path("*.gtf"), val(baseName)

    script:
    """
    stringtie -e -G ${params.annotation} -o ${baseName}.gtf ${mapped_bam}

    """
}

 workflow {
    samples = new File(params.configFilePath).readLines()*.split('\t');
    samples  = Channel.fromList(samples);
    hisat2Alignment(samples);
    assembled = Stringtier(samples)

}

```

# Answer
Nextflow does not have an `env` process directive. If `stringtie` is not already in your `$PATH`, you'll need to append it before running `nextflow`, for example:

```
export PATH="${PATH}:/data/selena/stringtie-2.2.1"
nextflow run de-pipeline.nf

```

However, you may prefer to avoid local tool installations altogether, and use either the `conda` or `container` process directives. The latter will however require a container runtime, like Docker or Singularity.

You may also consider including a 'sample name' or some sort of sample identifier in your config file, which I suspect is just a simple tab-separated file of filenames. This would greatly simplify your workflow. For example:

```
params.samplesheet = './samples.tsv'
params.genome = './genomes/GRCh38/GRCh38'
params.gtf = './Homo_sapiens.GRCh38.gtf'

genome_name = file( params.genome ).name
genome_path = file( params.genome ).parent

params.outdir = './results'

```

```
process HISAT2 {

    tag { sample }

    publishDir "${params.outdir}/HISAT2/${sample}", mode: 'copy'

    cpus 16

    input:
    tuple val(sample), path(reads)
    path genome_dir

    output:
    tuple val(sample), path("${sample}.sam")

    script:
    def (r1, r2) = reads

    """
    hisat2 \\
        -p ${task.cpus} \\
        -x "${genome_dir}/${genome_name}" \\
        -1 "${r1}" \\
        -2 "${r2}" \\
        -S "${sample}.sam"
    """
}

```

```
process SAMtools_sort {

    tag { sample }

    input:
    tuple val(sample), path(bam)

    output:
    tuple val(sample), path("${sample}.sorted.bam")

    """
    samtools sort \\
        -o "${sample}.sorted.bam" \\
        "${bam}"
    """
}

```

```
process StringTie {

    tag { sample }

    publishDir "${params.outdir}/StringTie/${sample}", mode: 'copy'

    input:
    tuple val(sample), path(bam)
    path gtf

    output:
    tuple val(sample), path("${sample}.gtf")

    """
    stringtie \\
        -e \\
        -G "${gtf}" \\
        -o "${sample}.gtf" \\
        "${bam}"
    """
}

```

```
workflow {

    gtf = file( params.gtf )

    Channel
        .fromPath( params.samplesheet )
        .splitCsv( sep: '\t' )
        .map { sample, fq1, fq2 ->
            def reads = tuple( file(fq1), file(fq2) )

            tuple( sample, reads )
        }
        .set { samples }

    HISAT2( samples, genome_path )

    SAMtools_sort( HISAT2.out )

    StringTie( SAMtools_sort.out, gtf )
}

```

Where the contents of the input sample sheet (e.g. `samples.tsv`) looks like:

```
sample_A    /path/to/reads/A_R1.fastq.gz    /path/to/reads/A_R2.fastq.gz
sample_B    /path/to/reads/B_R1.fastq.gz    /path/to/reads/B_R2.fastq.gz
sample_C    /path/to/reads/C_R1.fastq.gz    /path/to/reads/C_R2.fastq.gz

```

Results:

```
$ nextflow run main.nf
N E X T F L O W  ~  version 22.04.0
Launching `main.nf` [adoring_poisson] DSL2 - revision: 43e8c59949
executor >  local (9)
[60/1cba2f] process > HISAT2 (sample_A)        [100%] 3 of 3 ✔
[0a/931a07] process > SAMtools_sort (sample_A) [100%] 3 of 3 ✔
[6d/7689a2] process > StringTie (sample_A)     [100%] 3 of 3 ✔

```
> 2 votes

---
Tags: rna-seq, nextflow, errors
---


# Question
Title: sra toolkit fastq-dump comand not terminating
I ran the following command trying to download a fasta file from NIH but the command doesn't seem to terminate and I get no output file in my repository:

```
fastq-dump SRR10134980

```

I have installed SRA Toolkit following the instructions on this page.

# Answer
That's an 11GB file and fastq-dump is not a speedster. Without more info the likely answer is that it keeps downloading which might take several hours on poor connection. You can get fastq download links from sra-explorer.info -- That is more stable than this tool.
> 1 votes

---
Tags: sratoolkit
---


# Question
Title: Convert Abricate output tsv file to gff3 format
I have a tsv file generated from abricate (https://github.com/tseemann/abricate). I need to convert them to gff3 format with certain columns retained, certain columns reordered, while other columns deleted.

We are trying to use these gff3 files for downstream applications and for piping into other applications. However, we could not solve it.

Below are examples of my tsv files, what possibly we may need to do, and desired output files in gff3 format.

**Any help will be much appreciated.**

## Input tsv file:

FILE SEQUENCE START END GENE COVERAGE COVERAGE\_MAP GAPS %COVERAGE %IDENTITY DATABASE ACCESSION PRODUCT

UBird\_Cyou\_D3.fna BJCZ01000001.1 1866608 1867417 cdtB 1-810/810 =============== 0/0 100 90 vfdb CAD48850 (cdtB) cytolethal distending toxin B \[CDT (VF0185)\] \[Escherichia coli O157:H str. 493/89\] UBird\_Cyou\_D3.fna BJCZ01000001.1 1867414 1868190 cdtA 1-777/777 =============== 0/0 100 90.61 vfdb CAD48849 (cdtA) cytolethal distending toxin A \[CDT (VF0185)\] \[Escherichia coli O157:H str. 493/89\] UBird\_Cyou\_D3.fna BJCZ01000001.1 2245186 2246238 ompA 1-1041/1041 ========/====== 1/12 100 94.11 vfdb AAF37887 (ompA) outer membrane protein A \[OmpA (VF0236)\] \[Escherichia coli O18:K1:H7 str. RS218\]

## What we may need to do (there may be other ways too, I am not sure):

```
Row 1 (always starts with "#") - Need to replace with the string "##gff-version 3"
Col 1 - get rid of ".fna" and retain other data
Insert new Col - print the string from Col 11 for all rows
Col 2 - get rid of entire column
Insert new Col - print "CDS" for all rows
Col 3 - retain data
Col 4 - retain data
Insert new Col and print "." for all rows
Insert new Col and print "+" for all rows
Insert new Col and print "0" for all rows
Col 5 to Col 10 - get rid of all these columns and data
Col 11 - delete column
Col 13 - retain data except "(", ")", "[", "]"
Add new Col - Starting with "ID=" followed by the string taken from Col 1 and a underscore added (for the example data "UBird_Cyou_D3_") alongwith numerals starting from 1 and incrementing by 1. This column data needs to be appended by "product=" followed by data from the corresponding row of the modified Col 13. The separator between ID string and product string should be ";". After completion this column should be like "ID=UBird_Cyou_D3_1;product=cdtB cytolethal distending toxin B CDT VF0185 Escherichia coli O157:H str. 493/89"

```

## Desired final output (\*.gff3) considering the example data:

##gff-version 3

UBird\_Cyou\_D3 vfdb CDS 187 756 . + 0 ID=UBird\_Cyou\_D3\_1;product=cdtB cytolethal distending toxin B CDT VF0185 Escherichia coli O157:H str. 493/89

# Answer
Here's one way using , that I think fulfills the requirements. It adds each of the column names (on the first line) to an array to make accessing each of the fields a bit easier. This approach isn't strictly necessary, but it does make for a more readable solution in my opinion. With the following in a file called `tsv2gff3.awk`:

```
BEGIN {
    FS=OFS="\t"
    c=1
}

NR==1 && sub(/^#/, "") {
    for(i=1; i<=NF; i++) {
        a[$i]=i
    }

    print "##gff-version 3"
    next
}

{
    sub(/\.fna$/, "", $a["FILE"])
    gsub(/[\(\)\[\]]/, "", $a["PRODUCT"])

    id = "ID=" $a["FILE"] "_" c++
    product = "product=" $a["PRODUCT"]

    print \
        $a["FILE"],
        $a["DATABASE"],
        "CDS",
        $a["START"],
        $a["END"],
        ".",
        "+",
        "0",
        id ";" product
}

```

And the following in `file.tsv`:

```
#FILE   SEQUENCE    START   END GENE    COVERAGE    COVERAGE_MAP    GAPS    %COVERAGE   %IDENTITY   DATABASE    ACCESSION   PRODUCT
UBird_Cyou_D3.fna   BJCZ01000001.1  1866608 1867417 cdtB    1-810/810   =============== 0/0 100 90  vfdb    CAD48850    (cdtB) cytolethal distending toxin B [CDT (VF0185)] [Escherichia coli O157:H str. 493/89]
UBird_Cyou_D3.fna   BJCZ01000001.1  1867414 1868190 cdtA    1-777/777   =============== 0/0 100 90.61   vfdb    CAD48849    (cdtA) cytolethal distending toxin A [CDT (VF0185)] [Escherichia coli O157:H str. 493/89]
UBird_Cyou_D3.fna   BJCZ01000001.1  2245186 2246238 ompA    1-1041/1041 ========/====== 1/12    100 94.11   vfdb    AAF37887    (ompA) outer membrane protein A [OmpA (VF0236)] [Escherichia coli O18:K1:H7 str. RS218]

```

Run using `awk -f tsv2gff3.awk file.tsv`, results:

```
##gff-version 3
UBird_Cyou_D3   vfdb    CDS 1866608 1867417 .   +   0   ID=UBird_Cyou_D3_1;product=cdtB cytolethal distending toxin B CDT VF0185 Escherichia coli O157:H str. 493/89
UBird_Cyou_D3   vfdb    CDS 1867414 1868190 .   +   0   ID=UBird_Cyou_D3_2;product=cdtA cytolethal distending toxin A CDT VF0185 Escherichia coli O157:H str. 493/89
UBird_Cyou_D3   vfdb    CDS 2245186 2246238 .   +   0   ID=UBird_Cyou_D3_3;product=ompA outer membrane protein A OmpA VF0236 Escherichia coli O18:K1:H7 str. RS218

```
> 2 votes

---
Tags: gff3, tsv
---


# Question
Title: Why are asymmetric distributions an issue for GSEA?
I'm reading Subramanian et. al's (2005) original GSEA paper. In the paper, the authors mention the following:

> We noticed that the use of weighted steps could cause the distribution of observed ES scores to be asymmetric in cases where many more genes are correlated with one of the two phenotypes. We therefore estimate the significance levels by considering separately the positively and negatively scoring gene sets

I understand what they're saying, but I simply don't understand why it's a problem. Basically, they're saying that in some cases, most genes in the list we're testing could appear mostly at the top or bottom of the list. But I don't understand why having asymmetric distributions is a problem. What is the issue here?

EDIT: Here are some example GSEA plots I found via a web search with such gene sets. Where is the problem?

One thing I don't get though is why this problem would be introduced through "weighting" of the steps. I feel like it would also be a problem without weighting. The way it is phrased in the paper seems to suggest that it's the weighting that causes an issue.

# Answer
> 1 votes
The problem is that the enrichment score is modified based on whether the sets are positively or negatively associated, which means that an unbalanced set via either associative statistic or gene count (as is almost always the case) will lead to an incorrect total score.

Put in other words, the calculated score from the positive genes alters the calculated score from the negative genes.

Given this, it's a better idea to either separately calculate scores for the positive set and the negative set (this is the easiest approach, and works with existing methods), or (equivalently) recalculate enrichment score modifiers for positive and negative genes separately (i.e. adjust the graph) so that the zero point for the enrichment score crosses the zero point for association.

A common argument against this (when I have mentioned it to others) is that GSEA has *always* been done on the total gene set, so it makes sense to do it on the total gene set because it has always been done that way. I'm not a huge fan of the status-quo argument.

I alluded to this issue in this previous answer.

With regards to your figures:

**a**. Zero-point crossover is at an enrichment score of ~0.42, so the negative set will be underemphasised. It's not contributing much anyway (except for the one most negative gene at the end), so wouldn't affect the outcome much.

**b**. Zero-point crossover is at an enrichment score of ~0.45. Similar situation as in **a**.

**c**. Zero-point crossover is at an enrichment score of ~0.32. There aren't really any positive genes with substantial association scores, so no issue.

**d**. Plot looks to be identical to **c**, so the same comments apply.

The scores are most affected when they have a lot of area on both sides of the Y axis. The top right plot of this sub-figure is an example of that, showing a sharp-ish tail on the negative side, so the calculated enrichment score for the negative peak is less extreme than it would be if only the negative genes were considered \[source\]:

> The way it is phrased in the paper seems to suggest that it's the weighting that causes an issue.

Yes. If the positive genes are weighted separately from the negative genes, then the bias can be removed. I did briefly attempt a non-parametric weighting for GSEA (based on rank), but got too confused by the maths involved. Someone answered my question about it, but I haven't had to do any GSEA since then, so haven't actually verified if their solution works for me:

https://math.stackexchange.com/questions/3567714/probability-that-the-nth-highest-of-a-set-of-k-numbers-selected-from-m-numbers-h

---
Tags: gsea
---


title: complete the res type in sdf

text:
# Question
Title: complete the res type in sdf
I use the diffdock (link to diffdock on github) to dock a peptide with protein and this tool give me a ligand in sdf file format. I traslate this ligand sdf format file to a pdb format file and find the pdbfile lack res information. I wonder if tools/code can automaticly complete this file. I need the UNK become the residue types like pdb.This my error pdb file（this structure is a peptide）

```
HETATM    1  C   UNK     0       9.797  22.501  -6.484  0.00  0.00           C
HETATM    2  C   UNK     0      10.366  21.086  -6.446  0.00  0.00           C
HETATM    3  O   UNK     0       9.550  20.218  -7.135  0.00  0.00           O
HETATM    4  N   UNK     0      11.710  21.084  -6.978  0.00  0.00           N
HETATM    5  C   UNK     0      11.793  20.201  -8.123  0.00  0.00           C
HETATM    6  C   UNK     0      11.151  20.760  -9.346  0.00  0.00           C
HETATM    7  C   UNK     0      11.586  20.031 -10.581  0.00  0.00           C
HETATM    8  C   UNK     0      13.000  20.308 -10.986  0.00  0.00           C
HETATM    9  C   UNK     0      13.397  19.716 -12.309  0.00  0.00           C
HETATM   10  N   UNK     0      12.636  18.523 -12.633  0.00  0.00           N
HETATM   11  C   UNK     0      11.861  18.392 -13.830  0.00  0.00           C
HETATM   12  O   UNK     0      11.837  19.364 -14.611  0.00  0.00           O
HETATM   13  C   UNK     0      11.148  17.131 -14.065  0.00  0.00           C
HETATM   14  C   UNK     0      11.800  16.246 -15.085  0.00  0.00           C

```

# Answer
An SDF file lacks residue information, including atom name. Although for the latter, the `molFileAlias` atomic property is frequently used.

## Adding residue info

There are several ways to do this, such as in RDKit (blog post):

```
from rdkit import Chem
from rdkit.Chem import AllChem

from typing import List

with Chem.SDMolSupplier('.sdf') as mol_store:
     mols: List[Chem.Mol] = list(mol_store)

mol: Chem.Mol
for mol in mols:
    atom: Chem.Atom
    for i, atom in enumerate(mol.GetAtoms()):
        info = AllChem.AtomPDBResidueInfo(atomName=f'ATM{i}',
                                          residueName='LIG',
                                          residueNumber=1,
                                          chainId='X')
        atom.SetMonomerInfo(info)

mol: Chem.Mol
for mol in mols:
    name: str = mol.GetProp('_Name')
    Chem.MolToPDBFile(mol, f'{name}.pdb')

```

Notes:

* Atom names do have somewhat of a standard format, C1, C2, C3. I am being lazy in the above. If your ligand is present in the RCSB PDB, it will have a chemical component info on the PDB: example and it is often imperative to obey the atom names and atom index from there.
* PyMOL or pymol2 module can also add PDB information, therein the residue index is not an integer but a string as it combines w/ insertion code.
* In some tools, you lose the bond order, but if you have a SDF and not a PDB, then you are safe. But if you have a reference it's easy to fix.

```
from rdkit import Chem
from rdkit.Chem import AllChem
template: Chem.Mol = Chem.MolFromSmiles('')
AllChem.AssignBondOrdersFromTemplate(template, mol)

```
> 1 votes

---
Tags: pdb, peptide-shaker
---


# Question
Title: Any existing workflow on Terra that can run multinomial logistic regression?
I need to run multinomial logistic regression on a smaller set of variants e.g. using Trinculo. May I know if there's any existing workflow on Terra that can do this? Thanks.

# Answer
What I think is occurring is Terra is the Broad Institute's cloud service. Within that there is a Jupiter notebook and you then just use scikit-learn's (Python) LogisticRegression() using multinomial option. Then is simply ....

```
from sklearn.linear_model import LogisticRegression
import pandas as pd
model = LogisticRegression(multi_class='multinomial', solver='lbfgs')
df = pd.DataFrame ({'colA':[1,2,3], 'colB':[1,2,5] ... etc ...)
X = df[['colA', 'colB', 'colC']]
y = df[['colD']]
model.fit(X, y)
preds = model.predict(X)
print (preds)

```

This sort of operation. Does that make sense?
> 1 votes

---
Tags: logistic-regression
---


# Question
Title: Nextflow Error: failed to read header from "-"
I am trying to run my nextflow pipeline, and have gotten this error:
`samtools sort: failed to read header from "-"`

I'm not sure why this error is occurring.

```
#!/usr/bin/env nextflow
nextflow.enable.dsl=2
import groovy.io.*

params.t = 8

process joinFastq {
input:
path(ds_directory)

output:
tuple path("joined_${ds_directory.baseName}.fastq"), val("$ds_directory.baseName")

script:
"""
gunzip -c ${ds_directory}/*/fastq_pass/*.fastq.gz >
 joined_${ds_directory.baseName}.fastq
"""

}
 process minimap2 {
cpus 8
memory "20 GB"
input:
tuple path(fastq), val(baseName)

output:
tuple path("*.bam"), val(baseName)

  script:
   """
  minimap2 -ax splice -uf -k14 -t 8 ${params.genome} ${fastq} > mapped_${baseName}.sam
    samtools view -h -b mapped_${baseName}.sam > ${baseName}.bam
   samtools sort -o ${baseName}.bam > ${baseName}_sorted.bam
   """
   }

```

I have tried applying the `-h` flag to `samtools sort` in the script

```
   samtools sort -h -o ${baseName}.bam > ${baseName}_sorted.bam

```

and the same error is produced. I'm not sure why the header flag didn't solve this issue.

# Answer
This is not a Nextflow error. The problem is that without an input file, samtools `sort` tries to read from stdin. Using a recent `samtools`, you can however coordinate sort the SAM and write a sorted BAM using:

```
samtools sort -o "${baseName}.bam" "mapped_${baseName}.sam"

```

You may have been intending to pipe the output to `samtools sort`, which would avoid writing large SAM files and is usually preferable. A minimal example might look like:

```
params.reads = './*.fastq.gz'
params.reference = './gencode.v43.transcripts.fa'

```

```
process minimap2 {

    tag { sample }

    cpus 8
    memory 20.GB

    input:
    tuple val(sample), path(fastq)
    path ref_fasta

    output:
    tuple val(sample), path("${sample}.bam")

    script:
    def cpus = task.cpus > 1 ? task.cpus - 1 : task.cpus

    """
    minimap2 \\
        -ax splice \\
        -uf \\
        -k14 \\
        -t ${cpus} \\
        "${ref_fasta}" \\
        "${fastq}" |
    samtools sort \\
        -o "${sample}.bam"
    """
}

```

```
workflow {

    reads = Channel.fromFilePairs( params.reads, size: 1 )
    reference = file( params.reference )

    minimap2( reads, reference )
}

```

Also, be sure to avoid using absolute paths in your workflow scripts. For your workflow to be portable, processes should only access files in the working directory. The example above avoids this, by using a value channel, to stage the reference FASTA file into the working directory.

In case anyone is wondering, the ternary is to ensure that we do not oversubscribe the node(s) that these jobs land on. Otherwise, each job could use (up to) 900% (800% + 100% = 900%) CPU.
> 3 votes

---
Tags: rna-seq, samtools, nextflow
---


# Question
Title: Download all pubmed central article ids that have a keyword appearing in their titles/abstracts and also falling between two dates
I am using the following code to download all pubmed central article ids that have a keyword appearing in their titles/abstracts and also falling between two dates. The code seems to work when I select database = "pubmed" but when I switch to database = "pmc" I get empty list. For example, using the start and end dates as mentioned in the code I get 4 PubMed ids: \['36902449', '36880658', '36583446', '35239250'\] and when I check manually, the first one among them (which is '36902449') is also present in PMC database, of course, the id is different (10003631) but same article. The code is expected to fetch that id. I'd appreciate it if someone could point out what I'm missing here.

```
import os
import datetime
from Bio import Entrez

# Set the email address to identify yourself to NCBI
Entrez.email = "example@gmail.com"

# Set the date range for the search
start_date = datetime.datetime(2023, 3, 1)
end_date = datetime.datetime(2023, 3, 22)

# Format the dates as strings in the correct format
start_date_str = start_date.strftime("%Y/%m/%d")
end_date_str = end_date.strftime("%Y/%m/%d")

# Build the search query with the date range
keyword = "dosage compensation"
# database = "pubmed"
database = "pmc"

# Use esearch to fetch the list of matching article IDs
if database == "pubmed":
    search_term = f'("{start_date_str}"[Date - Publication] : "{end_date_str}"[Date - Publication]) AND {keyword}[Title/Abstract]'
    handle = Entrez.esearch(db=database,
                            term=search_term,
                            retmax=100,
                            api_key = os.environ['ENTREZ_API_KEY'])
elif database == "pmc":
    search_term = f'("{start_date_str}"[Publication Date] : "{end_date_str}"[Publication Date]) AND {keyword}[Title/Abstract] AND open_access[filter]'
    handle = Entrez.esearch(db=database,
                        term=search_term,
                        retmax=1000,
                        api_key = os.environ['ENTREZ_API_KEY'])

record = Entrez.read(handle)
id_list = record['IdList']
handle.close()

# Print the number of matching articles
print(f"Found {len(id_list)} articles published between {start_date_str} and {end_date_str} in {database}")
# print(id_list)

```

# Answer
> 1 votes
There seems to be a problem with using the `[Title/Abstract]` field. What you have looks correct and should produce the expected result, since it should capture articles where the search term appears in either the title or in the abstract, or in both. But it does not. Strangely, you get the expected result using `{keyword}[Abstract]`.

A workaround could be to `OR` the fields with:

```
f'("{keyword}"[Title] OR "{keyword}"[Abstract])'

```

Perhaps something like:

```
search_term = ' AND '.join([
    f'("{start_date_str}"[Publication Date] : "{end_date_str}"[Publication Date])',
    f'("{keyword}"[Title] OR "{keyword}"[Abstract])',
    f'(open_access[filter])',
])

```

---
Tags: python, biopython
---


# Question
Title: RNAseq: Why would using a more complete scaffold reduce DEGs?
I have a set of RNAseq data for a couple of relatively under-studied species.

They're fungi, of the order Hypocreales. The assemblies we're working with are de novo - reference seqs don't actually exist yet for most of our species. What the collaborator provided was updated genome assemblies in 60-90 fragments as opposed to the original 250-300, but obviously that's still not full resolution. As for the gene families, not a whole lot is known... in fact our species' genes are still not well described and we're mostly working off of predicted orthology.

During the analysis, a collaborator created new DNA scaffolds for these species which we'd now like to use.

However when I redid our RNAseq using the existing reads with the new scaffolds, one of our species had the number of differentially-expressed genes drop by anywhere from 20-50% per condition.

As the reads themselves haven't changed, we weren't expecting significant changes in differential gene expression, and indeed for our other species changes have been minimal.

Does anybody know what might cause this kind of change?

Thank you!

# Answer
Fungi, ain't a kingdom most of us will know much about. The only thing I can guess at is that in bacteria there has resurgence in methods for *de novo* assembly specifically aimed at reducing the total number of contigs in an assembly. The 'new wave' is summarised in this response here, Improving prokaryotic assembly with other contig/scaffold-level data?. This is potentially very exciting, except the OP reported they didn't work ;-(

A reasonable guess is your colleague has deployed these approaches to reduce the number of contigs. HOWEVER, what is appropriate for bacteria, may not be appropriate for fungi, alternatively it may be an improvement. Both bateria and fungi are small genomes, but beyond that there could notable differences.

If fungi are dominated by multigene families these might be being removed in the new wave *de novo* assemblies (they are designed for bacteria). The collapse of multigene families absolutely occurs in other eukaryotes with small genomes, but they have to be really big on multigene families (\>10% of genes in their genome). If that was true you would probably know about it.

Alternative, this could simply be poor gene identification in the first assembly resulting in an artificial jump in the number of genes.

**Solution?** Either way you should be able to figure which genes are different between the assemblies and precisely what type of genes are differing, e.g. are these large multi-gene families. Personally, I would pipe the predicted genes from both assemblies through `eggNOG`. This might be useful, but not really what it's designed for, but it will return the KEGG nos. More explicit methods based around comparing the RNAseq reads against the assembly via e.g. `bedtools` or these guys here recommend  `htseq`. This might be a better way to rephrase the question.

Thus I recommend:

1. Finding out which *de novo* methods were used between first and second fungi genome assemblies. Do any of those methods match this post, Improving prokaryotic assembly with other contig/scaffold-level data?
2. Consider a more explicit genomics study to examine RNAseq reads against the genome, or comparison between predicted no. genes between genomes.
3. Having said all that ... once you've done your DEG you identify the KEGG no. and then simply compare the lists directly to find out which genes are different between the assemblies. This is an easy database operation to perform, e.g. via a `merge` function (separate question).

I'd check point 1 and then go for point 3 - because it's something I easily know how to do. Thats just my opinion, a purist would likely explore this via point 2.

**Conclusion** What I am certain about is that for a species where very little might be known it is better to perform an investigation rather than make an assumption that might prove later to be wrong.
> 1 votes

---
Tags: rna-seq
---


# Question
Title: Problem with samtools command with a final hyphen
```
bwa mem ../01.data/GCF_000146045.2_R64_genomic.fna ../01.data/${sra}_1.fastq.gz ../01.data/${sra}_2.fastq.gz | \
    samtools view -S -b -t ../01.data/GCF_000146045.2_R64_genomic.fna.fai - | \
    samtools sort - >$sra.sorted.bam
samtools index $sra.sorted.bam

```

I am mapping the downloaded yeast sequences in sra format to the reference genome in a bash shell script. And I am using samtools to change to bam file. Then I try to sort the resulted bam file.

I could follow the steps but I could not understand the use of hyphen after the .fai (index file) and the use of this hyphen before the pipe. Moreover, what is the use of the hyphen after the samtools sort?

# Answer
It allows you to chain commands together to avoid creating intermediate files and simplify your workflow:

> Samtools is designed to work on a stream. It regards an input file \`-' as the standard input (stdin) and an output file \`-' as the standard output (stdout). Several commands can thus be combined with Unix pipes. Samtools always output warning and error messages to the standard error output (stderr).

This is a common idiom in the \*nix world, and is actually part of the POSIX specification:

> **Guideline 13:**
>
> For utilities that use operands to represent files to be opened for either reading or writing, the '-' operand should be used to mean only standard input (or standard output when it is clear from context that an output file is being specified) or a file named -.
> 4 votes

# Answer
The hyphen "-" character is commonly used as a placeholder in Unix commands to represent the standard input (stdin) or standard output (stdout) streams. From the POSIX specification:

> **Guideline 13:**
>
> For utilities that use operands to represent files to be opened for either reading or writing, the '-' operand should be used to mean only standard input (or standard output when it is clear from context that an output file is being specified) or a file named -.

In your pipeline, the hyphen after the index file name indicates that the output of the `bwa mem` command should be piped to `samtools view`, which will then read it from the standard input stream. The hyphen after the `samtools sort` command, is not strictly necessary since an input SAM/BAM/CRAM is entirely optional. Note that recent versions of `samtools` let you coordinate sort the SAM and write a sorted BAM, so you can do away with the call to `samtools view` entirely. For example:

```
set -euo pipefail

bwa mem \
    ../01.data/GCF_000146045.2_R64_genomic.fna \
    "../01.data/${sra}_1.fastq.gz" \
    "../01.data/${sra}_2.fastq.gz" |
samtools sort \
    -o "${sra}.sorted.bam"
samtools index \
    "${sra}.sorted.bam"

```

If you wanted to, you could replace the hyphens with the `/dev/stdin` special file, which represents the standard input streams. Sometimes this can make more complex commands a bit more readable, but note that not all Unix systems support the `/dev/stdin` and `/dev/stdout` special files. Most modern Unix systems do however.

The above uses `set -eo pipefail` to ensure pipeline errors are not masked. This exits immediately if any command in the pipeline exits with a non-zero status.
> 4 votes

---
Tags: samtools
---


# Question
Title: Histogram plot using R/python with two variables[BC and Control] and their corresponding gene names
How do we plot histogram in R\Python with two variables\[BC and Control\] and their corresponding gene names examples is below:

Thank you,

# Answer
Real data is better, particularly if genes names are wanted (complicated simulation) .... anyway there you go. Simulated values under a negative binomial which are pretty realistic to the data supplied (that'll be negative binomial, just not ranked).

```
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

np.random.seed(2)
s1 = list(np.random.negative_binomial(1, 0.1, 500))
s2 = list(np.random.negative_binomial(1, 0.5, 500))
mydict = {'A':s1, 'B':s2}
df = pd.DataFrame(mydict)
fig, ax = plt.subplots()

heights, bins = np.histogram(df['A'])
heights2, bins2 = np.histogram(df['B'], bins=bins)

width = (bins[1] - bins2[0])/3

ax.bar(bins[:-1], heights, width=width, facecolor='coral')
ax.bar(bins2[:-1]+width, heights2, width=width, facecolor='blue')

```

---

Here's horizontal bars (minor differences) with slightly different negative binomial values used ...

```
# same imports as before

np.random.seed(2)
s1 = list(np.random.negative_binomial(1, 0.1, 500))
s2 = list(np.random.negative_binomial(1, 0.3, 500))
mydict = {'A':s1, 'B':s2}
df = pd.DataFrame(mydict)
fig, ax = plt.subplots()

heights, bins = np.histogram(df['A'])
heights2, bins2 = np.histogram(df['B'], bins=bins)

width = (bins[1] - bins2[0])/3

ax.barh(bins[:-1], heights, width, facecolor='coral')
ax.barh(bins2[:-1]+width, heights2, width, facecolor='blue')
ax.invert_yaxis()
plt.show()

```

---

To adapt it you just need to provide an appropriate set of values that will then be fitted to a dataframe.
> 1 votes

---
Tags: r, python, proteins, ggplot2
---


title: Variants list to aminoacids

text:
# Question
Title: Variants list to aminoacids
I have the sequencing data of a Monkeypox virus sample. I have the list of variants extracted from Nextclade and I want to find the corresponding aminoacids substitution. For example I have the variant C1092T (nucleotide annotation) and I want to automatically convert it into OPG001:D162N (aminoacid annotation).

# Answer
> 1 votes
Neither `C1092T` nor `OPG001:D162N` seem to follow HGVS conventions. You could try to convert them to HGVS nomenclature standards then use modules such as hgvs to get amino acid changes from nucleotide changes. You will need reference sequence identifiers though, so it might be easier to go from chromosomal level changes to AA changes, as cDNA changes depend on the transcript (which may or may not be an issue with viral transcripts).

---
Tags: ngs, snp, format-conversion
---