# Question
Title: SyntaxError after end of Snakefile code?
I am trying to run Fastqc on multiple files using Snakemake but am receiving an ambiguous SyntaxError past the last line of code in the Snakefile.

Snakefile:

```
import pandas as pd

SAMPLE_TABLE=pd.read_csv('./SampleTable.txt',sep='\t')
SAMPLES=SAMPLE_TABLE["sample_name"]
SAMPLES=list(SAMPLES)
DIRECTORY=list(SAMPLE_TABLE["FolderName"])

rule run_fastqc:
    input:
        expand("/PathToInputsFolder/{dir}/{sample}_{{Rnumber}}_001.fastq.gz", zip, dir=DIRECTORY, sample=SAMPLES)
    output:
        expand("{sample}_{{Rnumber}}_001.fastqc.html", sample=SAMPLES)
        expand("{sample}_{{Rnumber}}_001.fastqc.zip", sample=SAMPLES)
    shell:
        "fastqc {input} --outdir=./qc/"

```

Code:

```
snakemake --cores 1 {.+}.fastqc.html {.+}.fastqc.zip

```

Output:

```
SyntaxError in line 16 of /PathToSnakefile/Snakefile:
invalid syntax (Snakefile, line 16)

```

How should I figure out where the error is? Thank you for your time!

# Answer
How to figure out what error messages are trying to tell you is one of the great mysteries… :-)

What version of Snakemake are you using? I get slightly different messages. With an ancient version (3.13.3) on one machine:

```
SyntaxError in line 14 of /path/to/Snakefile:
invalid syntax

```

With a more up-to-date version (6.9.1, running on Python 3.9):

```
SyntaxError in line 16 of /path/to/Snakefile:
invalid syntax. Perhaps you forgot a comma? (Snakefile, line 16)

```

Both of these are useful hints. So one piece of advice would be to update your Snakemake and/or Python if it is old.

Another useful technique when faced with a mysterious error like this is to temporarily remove parts of your *Snakefile* until the error goes away. The complicated part of your script towards the end of the file is the `output:` stanza, so you might try removing one of the `expand(…)` lines and seeing what happens. If the error goes away, then you have localised the problem and can put the line back and carefully check the syntax in that particular part of the *Snakefile*.
> 2 votes

# Answer
**You're missing a comma between the two `output` lines.** (John Marshall's more detailed syntax error message was spot-on.)

With a simplified version of your Snakefile:

```
SAMPLES=["s1", "s2", "s3"]
rule run_fastqc:
    output:
        expand("{sample}_{{Rnumber}}_001.fastqc.html", sample=SAMPLES)
        expand("{sample}_{{Rnumber}}_001.fastqc.zip", sample=SAMPLES)
    shell:
        "echo {output}"

```

You get the syntax error:

```
$ snakemake -nrp s1_R1_001.fastqc.html
SyntaxError in line 6 of /data/home/jesse/Snakefile:
invalid syntax (Snakefile, line 6)

```

With the comma:

```
SAMPLES=["s1", "s2", "s3"]
rule run_fastqc:
    output:
        expand("{sample}_{{Rnumber}}_001.fastqc.html", sample=SAMPLES), # <-- here
        expand("{sample}_{{Rnumber}}_001.fastqc.zip", sample=SAMPLES)
    shell:
        "echo {output}"

```

Now it's OK:

```
$ snakemake -nrp s1_R1_001.fastqc.html
Building DAG of jobs...
Job stats:
job           count    min threads    max threads
----------  -------  -------------  -------------
run_fastqc        1              1              1
total             1              1              1

[Fri Nov 12 08:54:54 2021]
rule run_fastqc:
    output: s1_R1_001.fastqc.html, s2_R1_001.fastqc.html, s3_R1_001.fastqc.html, s1_R1_001.fastqc.zip, s2_R1_001.fastqc.zip, s3_R1_001.fastqc.zip
    jobid: 0
    reason: Missing output files: s1_R1_001.fastqc.html
    wildcards: Rnumber=R1
    resources: tmpdir=/tmp

echo s1_R1_001.fastqc.html s2_R1_001.fastqc.html s3_R1_001.fastqc.html s1_R1_001.fastqc.zip s2_R1_001.fastqc.zip s3_R1_001.fastqc.zip
Job stats:
job           count    min threads    max threads
----------  -------  -------------  -------------
run_fastqc        1              1              1
total             1              1              1

This was a dry-run (flag -n). The order of jobs does not reflect the order of execution.

```

Each `expand()` makes a list, so at first Snakemake sees two lists back-to-back and doesn't know what to do with them. If you use a comma *and* give output names you could handle the two lists separately in your rule, but if you don't need to, you could also just us a single `expand` call, like:

```
output:
    expand("{sample}_{{Rnumber}}_001.fastqc.{ext}", sample=SAMPLES, ext=["html", "zip"])

```

**Other suggestions for troubleshooting:**

* If the above simplified Snakefile works but the real one still doesn't, edit them incrementally to be more and more similar and you'll hone in on whatever problems are left.
* I find trimming away complexity until the problem reveals itself is often easiest, so I'd keep simplifying the real case until that happens.
* `expand` is just a Python function that Snakemake provides you, so you can call it yourself outside of a rule definition and check what it returns.
* You can even comment out *all* the Snakemake-specific syntax and leave just Python code when debugging, to sanity-check the rest of the code (e.g. your inputs and outputs from `expand`).
> 1 votes

---
Tags: rna-seq, python, snakemake, quality-control, fastqc
---


# Question
Title: How to match peptide sequence from two different condition
I have to peptide output sample output i would give here ..

Sample A

```
dput(head(B5_ILS_Mus_PeptideSummary))
structure(list(N = c(1, 1, 1, 1, 1, 1), Unused = c(47.67, 47.67,
47.67, 47.67, 47.67, 47.67), Total = c(47.67, 47.67, 47.67, 47.67,
47.67, 47.67), `%Cov` = c(98.0000019073486, 98.0000019073486,
98.0000019073486, 98.0000019073486, 98.0000019073486, 98.0000019073486
), `%Cov(50)` = c(86.0000014305115, 86.0000014305115, 86.0000014305115,
86.0000014305115, 86.0000014305115, 86.0000014305115), `%Cov(95)` = c(86.0000014305115,
86.0000014305115, 86.0000014305115, 86.0000014305115, 86.0000014305115,
86.0000014305115), Accessions = c("sp|P20065|TYB4_MOUSE", "sp|P20065|TYB4_MOUSE",
"sp|P20065|TYB4_MOUSE", "sp|P20065|TYB4_MOUSE", "sp|P20065|TYB4_MOUSE",
"sp|P20065|TYB4_MOUSE"), Names = c("Thymosin beta-4 OS=Mus musculus OX=10090 GN=Tmsb4x PE=1 SV=1",
"Thymosin beta-4 OS=Mus musculus OX=10090 GN=Tmsb4x PE=1 SV=1",
"Thymosin beta-4 OS=Mus musculus OX=10090 GN=Tmsb4x PE=1 SV=1",
"Thymosin beta-4 OS=Mus musculus OX=10090 GN=Tmsb4x PE=1 SV=1",
"Thymosin beta-4 OS=Mus musculus OX=10090 GN=Tmsb4x PE=1 SV=1",
"Thymosin beta-4 OS=Mus musculus OX=10090 GN=Tmsb4x PE=1 SV=1"
), Used = c(NA, NA, NA, NA, NA, NA), Annotation = c(NA, NA, NA,
NA, NA, NA), Contrib = c(2, 2, 2, 2, 2, 2), Conf = c(99.0000009536743,
58.2899987697601, 97.5199997425079, 99.0000009536743, 98.9300012588501,
73.7299978733063), Sequence = c("EKNPLPSKETIEQEKQAGES", "KKTETQEKNPLPSKETIEQEKQAGES",
"KNPLPSKETIEQEKQAGES", "KTETQEKNPLPSKETIEQEKQAGES", "LKKTETQEKNPLPSKETIEQEKQAGES",
"PLPSKETIEQEKQAGES"), Modifications = c(NA_character_, NA_character_,
NA_character_, NA_character_, NA_character_, NA_character_),
    ProteinModifications = c(NA_character_, NA_character_, NA_character_,
    NA_character_, NA_character_, NA_character_), Cleavages = c(NA,
    NA, NA, NA, NA, NA), dMass = c(-0.0625222995877266, 0.180465996265411,
    -0.00865099020302296, -0.0651467964053154, -0.107698999345303,
    -0.0186963994055986), `Obs MW` = c(2241.05004882813, 2956.67944335938,
    2112.06127929688, 2828.3388671875, 3069.47534179688, 1869.91320800781
    ), `Obs m/z` = c(748.0239, 740.1771, 529.0226, 708.092, 768.3761,
    624.3117), `Theor MW` = c(2241.11254882813, 2956.4990234375,
    2112.06982421875, 2828.40405273438, 3069.5830078125, 1869.93200683594
    ), `Theor m/z` = c(748.044738769531, 740.132019042969, 529.024719238281,
    708.108276367188, 768.403015136719, 624.317932128906), `Theor z` = c(3,
    4, 4, 4, 4, 3), Sc = c(15, 8, 11, 21, 12, 10), Spectrum = c("1.1.1.372.12",
    "1.1.1.339.6", "1.1.1.358.3", "1.1.1.355.16", "1.1.1.346.14",
    "1.1.1.384.8"), `Acq Time` = c(4.2116, 3.539133, 3.91795,
    3.867617, 3.68345, 4.45125), `Intensity (Peptide)` = c(3775.35,
    110.58, 4374.93, 7285.52, 3073.51, 5776.43), PrecursorIntensityAcquisition = c(3626.33,
    37.12, 4301.05, 7124.23, 2706.13, 5533.12), `Apex Time (Peptide)` = c(4.2,
    3.59, 3.91, 3.85, 3.67, 4.44), `Elution Peak Width (Peptide)` = c(0.1,
    0.06, 0.08, 0.1, 0.08, 0.08), MS2Counts = c(3077.868, 477.5326,
    2262.827, 4694.357, 1759.267, 2618.842)), row.names = c(NA,
-6L), class = c("tbl_df", "tbl", "data.frame"))

```

Sample B

```
dput(head(B6_ILS_mus_PeptideSummary))
structure(list(N = c(1, 1, 1, 1, 1, 1), Unused = c(82.41, 82.41,
82.41, 82.41, 82.41, 82.41), Total = c(82.41, 82.41, 82.41, 82.41,
82.41, 82.41), `%Cov` = c(100, 100, 100, 100, 100, 100), `%Cov(50)` = c(87.9999995231628,
87.9999995231628, 87.9999995231628, 87.9999995231628, 87.9999995231628,
87.9999995231628), `%Cov(95)` = c(86.0000014305115, 86.0000014305115,
86.0000014305115, 86.0000014305115, 86.0000014305115, 86.0000014305115
), Accessions = c("sp|P20065|TYB4_MOUSE", "sp|P20065|TYB4_MOUSE",
"sp|P20065|TYB4_MOUSE", "sp|P20065|TYB4_MOUSE", "sp|P20065|TYB4_MOUSE",
"sp|P20065|TYB4_MOUSE"), Names = c("Thymosin beta-4 OS=Mus musculus OX=10090 GN=Tmsb4x PE=1 SV=1",
"Thymosin beta-4 OS=Mus musculus OX=10090 GN=Tmsb4x PE=1 SV=1",
"Thymosin beta-4 OS=Mus musculus OX=10090 GN=Tmsb4x PE=1 SV=1",
"Thymosin beta-4 OS=Mus musculus OX=10090 GN=Tmsb4x PE=1 SV=1",
"Thymosin beta-4 OS=Mus musculus OX=10090 GN=Tmsb4x PE=1 SV=1",
"Thymosin beta-4 OS=Mus musculus OX=10090 GN=Tmsb4x PE=1 SV=1"
), Used = c(NA, NA, NA, NA, NA, NA), Annotation = c(NA, NA, NA,
NA, NA, NA), Contrib = c(2, 2, 2, 2, 2, 2), Conf = c(99.0000009536743,
99.0000009536743, 99.0000009536743, 99.0000009536743, 94.4599986076355,
99.0000009536743), Sequence = c("DKPDMAEIEKFDKSKLK", "EKNPLPSKETIEQEKQAGES",
"ETQEKNPLPSKETIEQEKQAGES", "IEKFDKSKLKKTETQEKNPLPSKETIEQEKQAGES",
"KKTETQEKNPLPSKETIEQEKQAGES", "KNPLPSKETIEQEK"), Modifications = c(NA_character_,
NA_character_, NA_character_, NA_character_, NA_character_, NA_character_
), ProteinModifications = c(NA_character_, NA_character_, NA_character_,
NA_character_, NA_character_, NA_character_), Cleavages = c(NA,
NA, NA, NA, NA, NA), dMass = c(0.0129736997187138, -0.0542825013399124,
-0.0299271997064352, -0.107947997748852, -0.0191929005086422,
-0.0046052997931838), `Obs MW` = c(2021.06335449219, 2241.05834960938,
2599.2314453125, 4045.01391601563, 2956.47973632813, 1639.87341308594
), `Obs m/z` = c(506.2731, 748.0267, 650.8151, 810.0101, 592.3032,
547.6318), `Theor MW` = c(2021.05029296875, 2241.11254882813,
2599.26123046875, 4045.12182617188, 2956.4990234375, 1639.87805175781
), `Theor m/z` = c(506.269866943359, 748.044738769531, 650.822631835938,
810.031677246094, 592.307067871094, 547.63330078125), `Theor z` = c(4,
3, 4, 5, 5, 3), Sc = c(11, 12, 12, 26, 10, 10), Spectrum = c("1.1.1.858.2",
"1.1.1.621.5", "1.1.1.640.9", "1.1.1.643.4", "1.1.1.597.3", "1.1.1.604.4"
), `Acq Time` = c(8.74735, 4.069217, 4.43505, 4.4987, 3.615083,
3.7369), `Intensity (Peptide)` = c(1798.16, 2856.36, 1278.57,
5833.81, 758.71, 599.21), PrecursorIntensityAcquisition = c(1603.78,
443.55, 1104.72, 5265.07, 484.03, 423.02), `Apex Time (Peptide)` = c(8.72,
4.19, 4.39, 4.45, 3.55, 3.78), `Elution Peak Width (Peptide)` = c(0.14,
0.12, 0.12, 0.37, 0.1, 0.08), MS2Counts = c(456.3906, 313.663,
712.3438, 2748.186, 351.0362, 228.9387)), row.names = c(NA, -6L
), class = c("tbl_df", "tbl", "data.frame"))

```

Unlike genes expression where it can be intersected since i know the gene name or ensemble ID are consistent.

Here My objective is how many of the sequence are common or unique to these two condition.

In the **accession** column many ID are same but the differences are in the sequence column.

Is there any programmatic way of finding the elements that are common or different in both the conditions

# Answer
> 1 votes
`merge(a, b, by = "Sequence")` would "merge" the two tables (as the name suggests) based on the `Sequence` column. The caveat is that sequences have to be exactly the same to be considered "common to both samples", so even 1 aa difference would be missed.

Simplified output:

```
1       EKNPLPSKETIEQEKQAGES   1    47.67   47.67     98         86         86 sp|P20065|TYB4_MOUSE
2 KKTETQEKNPLPSKETIEQEKQAGES   1    47.67   47.67     98         86         86 sp|P20065|TYB4_MOUSE

```

---
Tags: proteomics
---


# Question
Title: How do I increase the minimum dot size in Seurat's DotPlot function?
I want to use the DotPlot function from Seurat v3 to visualise the expression of some genes across clusters. However when the expression of a gene is zero or very low, the dot size is so small that it is not clearly visible when printed on paper. Please is there a possibility to increase the minimum dot size in the DotPlot function to make the dot sizes more visible when printed?

For example, I would like to have a minimum dot size set to be like

from this publication

Rather than:

from this other publication

Thank you in advance for your helpful hint.

# Answer
Did you try to use `DotPlot(..., scale.by = "size")`?

In contrast to the default `scale.by= "radius"`, this will link the *area* (`==2*pi*r^2`), not the radius, of the circles to the fraction of cells expressing the feature. This corresponds much better to our perception of size and will make differences in low values easier to see.

Reading `?Seurat::DotPlot` the `scale.min` parameter looked promising but looking at the code it seems to censor the data as well.

Since Seurat's plotting functionality is based on `ggplot2` you can also adjust the color scale by simply adding `scale_fill_viridis()` etc. to the returned plot. This might also work for `size`. Try something like:

```
DotPlot(...) + scale_size(range = c(5, 10)) # will likely warn about supplying the same scale twice

```

Apart from this, Seurat's plotting system is not very hackable and I find it much easier to extract the relevant data and plot them myself with `ggplot2`.

Yet another comment: Your plot with the strong differences looks much more convincing to me wrt. to the marker property of these genes than the cited plot.
> 2 votes

# Answer
You can define your breaks with:

```
p + scale_size(breaks = c(0, 25, 50, 75))

```

You can adjust the size of the points with:

```
p + scale_size(range = c(1,8))     # Note: c(1,6) is the default

```

For more info, see `??scale_size`
> 0 votes

---
Tags: seurat, ggplot2
---


# Question
Title: find conservative regions relative to a specific accession
I have 7 `RefSeq` accessions which I performed `multiple sequence aligment` on it using `CLUSTALO`. These are the following:

```
NC_026431.1
NC_007363.1
NC_004907.1
NC_007377.1
NC_007367.1
NC_002016.1
NC_026427.1

```

I input the fasta sequence of those 7 accessions on `CLUSTALO` under this link. I get the result in either the `msf` format or `clustal` format with characters. The result looks like the following:

```
CLUSTAL O(1.2.4) multiple sequence alignment

NC_002016.1      -------AGCGAAAGCAGGTAGATATTGAAAGATGAGTCTTCTAACCGAGGTCGAAACGT   53
NC_007377.1      -------AGCAAAAGCAGGTAGATATTGAAAGATGAGCCTTCTAACCGAGGTCGAAACGT   53
NC_007367.1      -------AGCAAAAGCAGGTAGATATTGAAAGATGAGCCTTCTAACCGAGGTCGAAACGT   53
NC_004907.1      GGGGAATTCCAAAAGCAGGTAGATATTGAAAGATGAGTCTTCTAACCGAGGTCGAAACGT   60
NC_026427.1      --------------------------------ATGAGTCTTCTAACCGAGGTCGAAACGT   28
NC_007363.1      -------AGCAAAAGCAGGTAGATATTGAAAAATGAGTCTTCTAACCGAGGTCGAAACGT   53
NC_026431.1      --------------------------------ATGAGTCTTCTAACCGAGGTCGAAACGT   28
                                                 ***** **********************

```

I'm aware that the `*` symbols means a highly conserved region and that's exatcly what I'm looking for , however I need the highly conserved regions relative to a specific accession, namely `NC_007367.1`. I'm not sure how to do that. I'm also not sure what the `-` represents.

# Answer
> 2 votes
You can opt for a similarity matrix as output (here is an example). The fact that you would like to judge conservation "relative to a given sequence", you can use this matrix to rank similarities of different sequences with respect to any of these sequences.

I would advise visualizing these similarities, for example in terms of a heatmap, and seeing if anything stands out (with respect to your sequence of interest). And I would advise against selecting an arbitrary threshold (i.e. 90% or 95%), this would trigger questions based on the selection criteria of such a cutoff and this is understandable, such arbitrary cut-offs often do not mean much in biology.

---
Tags: multiple-sequence-alignment
---


# Question
Title: How to get UCSC Genome Browser to show the actual sequence?
Is it possible to make UCSC Genome Browser to show the sequence (nucleotide chain) of the current region?

I fiddled with the various viewing options and couldn't figure it out.

# Answer
If you zoom in at one specific region, it shows the nucleotide chain on the third line (below the chromosome number). It can't be there by default as the size of default region is more than 100,000bp. You can use the scale on the top of the page to see the size of the region.
> 5 votes

# Answer
In case it’s of any benefit to you, I’ll often check IGV when looking at tracks with that kind of zoomage. Not that it zooms any better than the UCSC browser (although the bases are colour coded!) but you can see your aligned reads individually if you want, along with directional information etc etc. Loads of stuff.
> 0 votes

# Answer
From the top bar: View -\> DNA
> 0 votes

---
Tags: genome-browser
---


# Question
Title: How can I get a list of neighboring Hydrogen atoms of an alpha-carbon?
Suppose, I want to find the neighboring Hydrogen atoms of a specific C-alpha atom.

The following was my effort. However, it produces a bunch of empty lists.

```
from typing import List
from Bio.PDB import Chain, Atom, NeighborSearch, PDBParser

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

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

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

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

```

How can I correct this?

# Answer
You posted similar question here: https://stackoverflow.com/questions/69676943/how-can-i-get-a-list-of-neighboring-hydrogen-atoms-of-an-alpha-carbon/69815509#69815509

ok got some code :

```
from Bio.PDB import Chain, Atom, NeighborSearch, PDBParser, Selection

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

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

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

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

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

hydrogen_list = get_hydrogen_atoms(c_alpha_list[0], po)

print('hydrogen_list : \n',hydrogen_list)

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

```

had to use a different PDB file because I needed the H atoms,

pH8\_5\_1bq9\_short.pdb :

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

```

output :

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

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

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

```

the selection list contains both ATOM and HETATM

have a look at : https://www.biostars.org/p/269579/ How to use biopython neighborsearch?

and Bio.PDB.NeighborSearch module :https://biopython.org/docs/1.75/api/Bio.PDB.NeighborSearch.html ;

Bio.PDB.Atom module :https://biopython.org/docs/1.75/api/Bio.PDB.Atom.html
> 2 votes

---
Tags: python, protein-structure, biopython
---


# Question
Title: Error while parsing gene bank file using Biopython
I am trying to extract the protein sequence of specific genes from gene bank like format file obtained from antismash part of which looks like this

```
LOCUS       scaffold_10            47160 bp    DNA     linear   UNK 01-JAN-1980
DEFINITION  scaffold_10.
ACCESSION   scaffold_10
VERSION     scaffold_10
KEYWORDS    .
SOURCE
  ORGANISM
            .
COMMENT     ##antiSMASH-Data-START##
            Version      :: 6.0.1-25e43a7(changed)
            Run date     :: 2021-11-05 17:39:45
            NOTE: This is a single cluster extracted from a larger record!
            Orig. start  :: 2049841
            Orig. end    :: 2097001
            ##antiSMASH-Data-END##
FEATURES             Location/Qualifiers
     protocluster    1..47160
                     /aStool="rule-based-clusters"
                     /contig_edge="False"
                     /core_location="join{[2069841:2070128](+),
                     [2070175:2070449](+), [2070507:2071973](+),
                     [2072023:2074178](+), [2074232:2074784](+),
                     [2074838:2077001](+)}"
                     /cutoff="20000"
                     /detection_rule="cds(PKS_AT and (PKS_KS or ene_KS or mod_KS
                     or hyb_KS or itr_KS or tra_KS))"
                     /neighbourhood="20000"
                     /product="T1PKS"
                     /protocluster_number="1"
                     /tool="antismash"
     proto_core      join(20001..20287,20335..20608,20667..22132,22183..24337,
                     24392..24943,24998..27160)
                     /aStool="rule-based-clusters"
                     /tool="antismash"
                     /cutoff="20000"
                     /detection_rule="cds(PKS_AT and (PKS_KS or ene_KS or mod_KS
                     or hyb_KS or itr_KS or tra_KS))"
                     /neighbourhood="20000"
                     /product="T1PKS"
                     /protocluster_number="1"
     cand_cluster    1..47160
                     /SMILES="CC(=O)C(=O)O"
                     /candidate_cluster_number="1"
                     /contig_edge="False"
                     /detection_rules="cds(PKS_AT and (PKS_KS or ene_KS or
                     mod_KS or hyb_KS or itr_KS or tra_KS))"
                     /kind="single"
                     /product="T1PKS"
                     /protoclusters="1"
                     /tool="antismash"
     region          1..47160
                     /candidate_cluster_numbers="1"
                     /contig_edge="False"
                     /product="T1PKS"
                     /region_number="1"
                     /rules="cds(PKS_AT and (PKS_KS or ene_KS or mod_KS or
                     hyb_KS or itr_KS or tra_KS))"
                     /tool="antismash"
     CDS             1899..2720
                     /ID="CDS_13110"
                     /NRPS_PKS="Domain: PKS_KR (26-149). E-value: 7.6e-08.
                     Score: 24.2. Matches aSDomain:
                     nrpspksdomains_jgi.p_FusspF11_1_766668_PKS_KR.1"
                     /NRPS_PKS="type: other"
                     /gene="jgi.p_FusspF11_1_766668"
                     /gene_functions="biosynthetic-additional (smcogs)
                     SMCOG1001:short-chain dehydrogenase/reductase SDR (Score:
                     161.6; E-value: 3.9e-49)"
                     /gene_kind="biosynthetic-additional"
                     /phase="0"
                     /source="prediction"
                     /transl_table=1
                     /translation="MTPPPKSALESDPAVGAQRRFSVTGNAVVTGGAGVLGLHACDALL
                     EHGLEGLMILDVNPAQSQGQITSLQNKFPRAKIMALKVDVTDENVVNAAMEETARVLGS
                     IDTLICFVGVVGCVETLDMPVPQWRKILDINTTGSFICAQAAARQMVKRGRGGSIVFVA
                     SISAHRVNYPQPQAAYNVSKSALLMLKSCLAAEWARYGIRTNSISPGYMDTILNEGDGI
                     AEHRKIWAEHNPSGRMGAPSELTGTVVLLASSAGSYINGADIVVDGGGIVL"
     aSDomain        1977..2345
                     /aSDomain="PKS_KR"
                     /aSTool="nrps_pks_domains"
                     /database="nrpspksdomains.hmm"
                     /detection="hmmscan"
                     /domain_id="nrpspksdomains_jgi.p_FusspF11_1_766668_PKS_KR.1
                     "
                     /evalue="7.60E-08"
                     /label="jgi.p_FusspF11_1_766668_PKS_KR.1"
                     /locus_tag="jgi.p_FusspF11_1_766668"
                     /protein_end="149"
                     /protein_start="26"
                     /score="24.2"
                     /specificity="KR activity: inactive"
                     /specificity="KR stereochemistry: C1"
                     /tool="antismash"
                     /translation="AVVTGGAGVLGLHACDALLEHGLEGLMILDVNPAQSQGQITSLQN
                     KFPRAKIMALKVDVTDENVVNAAMEETARVLGSIDTLICFVGVVGCVETLDMPVPQWRK
                     ILDINTTGSFICAQAAARQ"
     aSModule        1977..2345
                     /domains="nrpspksdomains_jgi.p_FusspF11_1_766668_PKS_KR.1"
                     /incomplete
                     /locus_tags="jgi.p_FusspF11_1_766668"
                     /tool="antismash"
                     /type="pks"
     CDS_motif       1980..2042
                     /aSTool="nrps_pks_domains"
                     /database="abmotifs"
                     /detection="hmmscan"
                     /domain_id="nrpspksmotif_jgi.p_FusspF11_1_766668_0001"
                     /evalue="1.90E-05"
                     /label="PKSI-KR_m1"
                     /locus_tag="jgi.p_FusspF11_1_766668"
                     /protein_end="48"
                     /protein_start="27"
                     /score="17.3"
                     /tool="antismash"
                     /translation="VVTGGAGVLGLHACDALLEHG"
     CDS             join(3252..3263,3321..3443,3493..4428)
                     /ID="CDS_13111"
                     /NRPS_PKS="Domain: PKS_ER (28-277). E-value: 6.4e-13.
                     Score: 40.2. Matches aSDomain:
                     nrpspksdomains_jgi.p_FusspF11_1_794272_PKS_ER.1"
                     /NRPS_PKS="type: other"
                     /gene="jgi.p_FusspF11_1_794272"
                     /gene_functions="biosynthetic-additional (smcogs)
                     SMCOG1040:alcohol dehydrogenase (Score: 275.7; E-value:
                     7.6e-84)"
                     /gene_kind="biosynthetic-additional"
                     /phase="0"
                     /source="prediction"
                     /transl_table=1
                     /translation="MSQTNLSCVLYGPGKARFENRPVPSLKDAHDVIIRISYVGVCGSD
                     VHFWTDGGFARKVSEDQPLVMGHEASGIVRSIGPDVTLLKPGDRVAIEPGFSCRRCKQC
                     KDGRYNLCPKMKFAADPPLTQGTLSRFFSIPEDFAYKIPDSLSLEEAVLVEPLAVAVHG
                     IRLAGLEVGQRVLVQGSGTIGLLTAAVAKAYGAKQVYITDVNLDKIKFAKKYLECSAFI
                     PDLGSTPEENAARFKTETGLDDGVDAVIECTGVEASTQTGLLALSAGGVLVQVGLGKPV
                     QAIPIHAMSEKEIVLKTSFRYGPGDYEIALELLESGKVSVRPLISSITPFEKATEAWEK
                     TRKGKGIKNLIRGVQD"

```

Now from this I want to extract the protein sequence of gens having PKS domain (PKS\_ER/PKS\_AT/PKS\_KS ...). So far I have tried the following biopython code but it is giving me error

```
from Bio import SeqIO
from Bio.SeqFeature import SeqFeature

gbk = "myfile_predicted_SMclusters.gbk"
fa = "myfile_PKS_protein.fa"
input_handle = open(gbk, "r")
output_handle = open(fa, "w")

for record in SeqIO.parse(input_handle, "genbank"):
    features = [feature for feature in record.features if feature.type == "CDS"]
    for feature in features:
        if feature.qualifiers["NRPS_PKS"][1] == "PKS":
        assert len(seq_feature.qualifiers['translation'])==1
        output_handle.write(">%s %s\n%s\n" % (
            record.description,
            feature.location,
            seq_feature.qualifiers['translation'][0]))

output_handle.close()
input_handle.close()

```

The error is as follows:

```
if feature.qualifiers["NRPS_PKS"][1] == "PKS":
KeyError: 'NRPS_PKS'

```

The link to the sample inputfile is below: sample my\_file.gbk

I'm using python 3.8.8

Please help!!!

# Answer
> 1 votes
I see 2 problems.

First, you need to check if the qualifier your are searching for exists in the given CDS feature.

Second, the qualifier holds more then just the string `PKS` for example `Domain: PKS_AT (537-846). E-value: 1.5e-98`.

You could modify your code like this:

```
feature_type = "CDS"
feature_qualifier = "NRPS_PKS"
target = "PKS"

def has_pks_value(feature, qualifier, part):
    return any(part in target for target in feature.qualifiers[qualifier])

for record in SeqIO.parse(input_handle, "genbank"):
    for feature in record.features:
        if feature.type == feature_type and feature_qualifier in feature.qualifiers.keys() and has_pks_value(feature, feature_qualifier, target) and len(feature.qualifiers['translation'])==1:
            output_handle.write(f">{record.description} {feature.location.start} {feature.location.end} {feature.location.strand} \n {SeqFeature(feature.location).extract(record.seq)}")

```

I've created some variables which for your string. This might prevent you from making some typos.

Also instead of a list comprehension, I just run over all features and checked with a single `if` statement if the asked qualifier exists and contains the given value.

I also changed the string formatting to f-string. But I think you should create the Fasta Files using the `SeqIO` module of Biopython.

As result I got 37 sequences

---

Update - I added a missing method to check for the correct value of the qualifier

---
Tags: python
---


# Question
Title: How to map list items which are in groups in the dictionary values and return keys corrosponding it?
```
dict1 = {'chr1:1000-2000': ['A', 'B', 'C'], 'chr2:999-2000': ['A', 'X', 'Y'], 'chr3:1200-3000': ['B', 'A', 'C']}

list1 = [('A', 'B'), ('A', 'C'), ('A', 'X'), ('A', 'Y'), ('B', 'C'), ('B', 'X'), ('B', 'Y'), ('C', 'X'), ('C', 'Y'), ('X', 'Y')]

```

I want to map the items in list1 eg('A','B') both in dict1 values and return keys in wich this group occur

# Answer
I am assuming you only want keys where *both* tuple items are present. If so, this should work:

```
no_match = []
for item in list1:
    item1, item2 = item
    wanted_keys = [
        key for key in dict1.keys() if
        item1 in dict1[key] and
        item2 in dict1[key]
    ]
    if len(wanted_keys) > 0:
        print("%s:  %s" % (item, wanted_keys))
    else:
        no_match.append(item)
print("No matches found for items:", no_match)

```

The output looks like this:

```
('A', 'B'):  ['chr1:1000-2000', 'chr3:1200-3000']
('A', 'C'):  ['chr1:1000-2000', 'chr3:1200-3000']
('A', 'X'):  ['chr2:999-2000']
('A', 'Y'):  ['chr2:999-2000']
('B', 'C'):  ['chr1:1000-2000', 'chr3:1200-3000']
('X', 'Y'):  ['chr2:999-2000']
No matches found for items: [('B', 'X'), ('B', 'Y'), ('C', 'X'), ('C', 'Y')]

```

There is probably a better, more elegant way of doing this, I am not a python person, but this is the best I could come up with.
> 1 votes

---
Tags: python, sequence-analysis
---


# Question
Title: Does anyone know what tools are used to construct this beautiful genome map?
Does anyone know what tools are used to construct this beautiful genome map? The research paper didn't mention how to construct this map. Link to paper: https://www.sciencedirect.com/science/article/pii/S2052297515000773?via%3Dihub

# Answer
Circos: http://circos.ca/ (or a copy of it written for another language, depending on what the paper's authors used — probably best to check with them).
> 6 votes

# Answer
I tried to used Circos for my MSc Thesis but didn't have enough time to actually learn how to use so I managed to find a less fancy but also cool version of it called ClicO.FS which is an interactive web-based service of Circos (as the title of the article below says).
> 2 votes

# Answer
Possibly BRIG (blast ring image generator) was used

http://brig.sourceforge.net/

I also catalog other tools here if interested https://cmdcolin.github.io/awesome-genome-visualization/?latest=true&selected=%23BRIG
> 2 votes

---
Tags: sequence-annotation
---


# Question
Title: Synteny test and tools to use for mammalian organisms?
I am trying to do a synteny test for MHC genes in human and another mammalian organism

For this reason, I need to download all paralogs of rabbits from humans and then construct position the coordinates of paralog to test the synteny. However, I am unsure how to download the paralogs from Ensembl and how to perform this test?

My main aim is to research synteny between MHC genes of humans, rabbits and mice.

Any suggestions on how to do this or alternative bioinformatic methods on how to perform a synteny test would be highly appreciated!

# Answer
> 1 votes
This is not an easy ask. Firstly, very few species (I think only human and mouse) have alternative haplotypes annotated for MHC. Most species only have the primary assembly, only one possible haplotype of this locus. This means that the haplotype in the primary assembly is not necessarily the same haplotype as you'll find in another species. Since the MHC haplotypes can contain completely different genes, you will necessarily be able to find orthologues and synteny.

Secondly, we in Ensembl do not calculate homologues of the genes on the alternative loci, so we only have homologues of the genes that just happen to be on the primary assembly.

There is already pre-calculated synteny in Ensembl. Here's the human MHC locus compared to Rabbit. This includes all the homologous genes across this locus. However, this is just one arbitrary MHC vs another arbitrary MHC.

---
Tags: ensembl, human-genome
---


title: Peptide data visualization

text:
# Question
Title: Peptide data visualization
I'm trying to replicate the figure given in this paper using this table which they have given.

So what i understand from the figure is they have calculated the zscore between two condition and plotted the same.

Now How is the length of each bar which is the peptide here is variable is it the zscore or the actual peptide length?

Small subset of the data

```
dput(head(d2))
structure(list(Peptide = structure(1:6, .Label = c("AAHLPAEFTPAVH",
"AAHLPAEFTPAVHASLD", "AAWGKVG", "AEFTPAVHA", "AEFTPAVHASLDKFLASVSTVLTSKYR",
"AFSDGLAH", "AFSDGLAHLD", "AFSDGLAHLDNLK", "AGEYG", "AHAGEYG",
"AHAGEYGAEALE", "AHAGEYGAEALERMF", "AHAGEYGAEALERMFLSFPTTK",
"AHGKKVLG", "AHLDNLK", "AHLPAEFTPAVH", "AHLPAEFTPAVHASLD", "AHLPAEFTPAVHASLDKF",
"AHVDDMPNA", "AHVDDMPNALS", "ALSDLH", "ALWGKVNVDEVG", "ALWGKVNVDEVGGEA",
"ALWGKVNVDEVGGEALG", "ALWGKVNVDEVGGEALGRLL", "ANALAHKYH", "ASLDKF",
"ASLDKFLA", "ASLDKFLASVS", "ASLDKFLASVSTVLT", "ASLDKFLASVSTVLTS",
"ASLDKFLASVSTVLTSKYR", "ATLS", "AVAHVDD", "AVAHVDDMPN", "AVAHVDDMPNA",
"AVAHVDDMPNAL", "AVAHVDDMPNALS", "AVAHVDDMPNALSA", "AYQKVVAGVANALAHK",
"DALT", "DALTN", "DDMPNA", "DDMPNALS", "DEVG", "DEVGGEA", "DEVGGEALG",
"DEVGGEALGRLL", "DGLA", "DGLAH", "DGLAHLD", "DGLAHLDNLK", "DKLHVD",
"DKLHVDPEN", "DKLHVDPENF", "DKLHVDPENFR", "DKLHVDPENFRL", "DKLHVDPENFRLL",
"DKLHVDPENFRLLG", "DKTNVK", "DLH", "DLS", "DLSTPD", "DLSTPDA",
"DMPNALS", "DPENF", "DPENFR", "DPVNF", "DPVNFK", "EFTPAVH", "EFTPAVHASLD",
"EFTPAVHASLDKFLASVSTVLTSKYR", "EFTPPVQ", "EFTPPVQA", "EFTPPVQAA",
"ESFG", "ESFGDL", "ESFGDLSTPDAVM", "ESFGDLSTPDAVMG", "EVGG",
"EVGGEALG", "FDLS", "FESF", "FESFG", "FESFGDLS", "FESFGDLSTPDAVMGNPKVK",
"FGDL", "FGDLS", "FPHF", "FPHFD", "FPHFDL", "FPHFDLS", "FPT",
"FPTT", "FPTTKT", "FPTTKTY", "FPTTKTYFPHF", "FSDGLA", "FSDGLAH",
"FSDGLAHLD", "FSDGLAHLDNLK", "FTPAVH", "FTPAVHASLD", "FTPAVHASLDKFLASVSTVLT",
"FTPPVQ", "FTPPVQA", "FTPPVQAA", "GAFSDGLA", "GAFSDGLAH", "GAFSDGLAHLDNLK",
"GEALG", "GEYG", "GGEAL", "GKKVL", "GKKVLG", "GKVNVD", "GKVNVDE",
"GKVNVDEVG", "GKVNVDEVGGEA", "GKVNVDEVGGEAL", "GKVNVDEVGGEALG",
"GKVNVDEVGGEALGRL", "GKVNVDEVGGEALGRLL", "GLAHLDNLK", "GNPKV",
"GNPKVK", "GNVL", "GRLL", "GTFAT", "GVANALAHKYH", "HAGEYG", "HAGEYGAEALE",
"HAGEYGAEALERMF", "HFDLS", "HHFG", "HKYH", "HLD", "HLDNLK", "HLPAE",
"HLPAEF", "HLPAEFTPAVH", "HLPAEFTPAVHASLD", "HLTPEEK", "HVDD",
"HVDDMP", "HVDDMPN", "HVDDMPNA", "HVDDMPNAL", "HVDDMPNALS", "HVDDMPNALSA",
"HVDDMPNALSALS", "HVDP", "HVDPEN", "HVDPENF", "HVDPENFR", "HVDPENFRL",
"HVDPENFRLLG", "KAAWGKVG", "KEFTPPVQ", "KEFTPPVQA", "KEFTPPVQAA",
"KEFTPPVQAAYQ", "KLHVDPENFRLLG", "KTYFP", "KTYFPH", "KTYFPHF",
"KTYFPHFD", "KVNVDEVGG", "KVNVDEVGGEA", "KVNVDEVGGEALG", "LDNL",
"LDNLK", "LGAFSDGLA", "LHVDP", "LHVDPENF", "LHVDPENFRLLG", "LKGTFATL",
"LLVVYPWTQRF", "LPAE", "LPAEF", "LPAEFTPAVH", "LPAEFTPAVHASLD",
"LPAEFTPAVHASLDK", "LPAEFTPAVHASLDKF", "LSALS", "LSFPTTK", "LSFPTTKT",
"LSFPTTKTY", "LSFPTTKTYFPHF", "LSFPTTKTYFPHFDL", "LSFPTTKTYFPHFDLS",
"LSFPTTKTYFPHFDLSHG", "LSTPDA", "LTPEEK", "LTPEEKSAVT", "LTPEEKSAVTA",
"LTPEEKSAVTAL", "LTPEEKSAVTALW", "LTPEEKSAVTALWG", "LTPEEKSAVTALWGKVNVD",
"LTPEEKSAVTALWGKVNVDEVG", "LTPEEKSAVTALWGKVNVDEVGGEA", "LTPEEKSAVTALWGKVNVDEVGGEALG",
"LVVYPWTQRF", "LVVYPWTQRFF", "LWGKVNVDEVGGEALG", "NAVAHVDD",
"NAVAHVDDMPNA", "NAVAHVDDMPNAL", "NAVAHVDDMPNALS", "NAVAHVDDMPNALSA",
"NAVAHVDDMPNALSALS", "NAVAHVDDMPNALSALSDL", "NAVAHVDDMPNALSALSDLHAHKLR",
"NVDEVG", "NVDEVGG", "NVDEVGGEA", "NVDEVGGEALG", "PADKTNVK",
"PAEF", "PAEFTPAVH", "PAEFTPAVHASLD", "PAVH", "PDAV", "PDAVM",
"PDAVMG", "PEE", "PEEK", "PEEKSAVT", "PEEKSAVTAL", "PENF", "PENFR",
"PENFRLLG", "PHF", "PHFD", "PHFDL", "PHFDLS", "PKVK", "PNAL",
"PNALS", "PPVQ", "PPVQA", "PPVQAA", "PVNF", "PVNFK", "PVQ", "PVQA",
"PWTQ", "PWTQRF", "PWTQRFF", "QVKGHGKKVADALTN", "RVDPVNF", "SDGLA",
"SDGLAHLDN", "SDGLAHLDNLK", "SDLH", "SFGD", "SFGDL", "SFGDLS",
"SFGDLSTPDAVMG", "SLDKF", "SPADK", "SPADKTNVK", "SVSTVLT", "TKTYFPHF",
"TNVK", "TPAVH", "TPAVHASLD", "TPAVHASLDKFLASVSTVLT", "TPDAVM",
"TPDAVMG", "TPDAVMGNPKVK", "TPEEKSAVT", "TPEEKSAVTAL", "TPEEKSAVTALW",
"TPPVQ", "TSKY", "TVLT", "TYFP", "TYFPH", "TYFPHF", "TYFPHFD",
"TYFPHFDLS", "VAHVDD", "VAHVDDMPNA", "VAHVDDMPNALS", "VAHVDDMPNALSA",
"VDDMP", "VDDMPN", "VDDMPNA", "VDDMPNAL", "VDDMPNALS", "VDEVG",
"VDEVGG", "VDEVGGEA", "VDEVGGEALG", "VDPENF", "VDPENFR", "VDPENFRL",
"VDPENFRLL", "VDPENFRLLG", "VDPVN", "VDPVNF", "VDPVNFK", "VDPVNFKL",
"VDPVNFKLL", "VDPVNFKLLS", "VHL", "VHLT", "VHLTP", "VHLTPE",
"VHLTPEE", "VHLTPEEK", "VHLTPEEKSAVT", "VHLTPEEKSAVTAL", "VHLTPEEKSAVTALWGKVNVD",
"VHLTPEEKSAVTALWGKVNVDE", "VHLTPEEKSAVTALWGKVNVDEVG", "VHLTPEEKSAVTALWGKVNVDEVGGEA",
"VHLTPEEKSAVTALWGKVNVDEVGGEALG", "VHLTPEEKSAVTALWGKVNVDEVGGEALGRLL",
"VKAAWGKVGAHAGEY", "VLSP", "VLSPA", "VLSPAD", "VLSPADK", "VLSPADKT",
"VLSPADKTN", "VLSPADKTNVK", "VLSPADKTNVKAA", "VLSPADKTNVKAAWGKVG",
"VNVDE", "VNVDEVG", "VNVDEVGG", "VNVDEVGGEA", "VNVDEVGGEAL",
"VNVDEVGGEALG", "VNVDEVGGEALGRL", "VVAGVANALA", "VVYPWTQRF",
"VVYPWTQRFF", "VVYPWTQRFFE", "VVYPWTQRFFESFGDLS", "VYPWTQRF",
"VYPWTQRFF", "VYPWTQRFFE", "WGKVG", "WGKVNVD", "WGKVNVDE", "WGKVNVDEVG",
"WGKVNVDEVGGEA", "WGKVNVDEVGGEAL", "WGKVNVDEVGGEALG", "WGKVNVDEVGGEALGRL",
"WGKVNVDEVGGEALGRLL", "YFPHFDLS", "YPWT", "YPWTQ", "YPWTQR",
"YPWTQRF", "YPWTQRFF", "YPWTQRFFE"), class = "factor"), C2 = c(-1.13005927266738,
-1.15451123584162, -0.771975424114207, -0.891503793435265, 0.686928914375861,
-1.15469327326147), C4 = c(0.770515719059141, 0.559148858503814,
1.12965461223673, 1.08129555353272, 0.46033682419062, 0.58089396147208
), C6 = c(0.359543553608241, 0.595362377337809, -0.357679188122518,
-0.18979176009746, -1.14726573856648, 0.573799311789392)), row.names = c(NA,
6L), class = "data.frame")

```

Any suggestion or help would be really appreciated

# Answer
The length of the bars and their placement on the x-axis represents the length of the peptide and mapping to the source protein, respectively. The coloring represents the average difference in peptide MS intensity between the control and two experiments.
> 2 votes

---
Tags: visualization, proteomics
---


# Question
Title: FindMarkers from Seurat returns p values as 0 for highly significant genes
I have been working on FindMarkers function for identifying significant genes in the cluster. But some Significant genes have very low p values, so they are returned as 0 in the output.Any value less than the .machine$double.xmin limit on R(that is 1e-305). So this output cannot be used for further plots like Volcano plots for visualization. Any solution to get p values for all the genes?

```
cd4.markers <- FindMarkers(sc.combined, ident.1 = "CD4+ T cell", min.pct = 0.25)

```

Part of the cd4.markers output

# Answer
> 1 votes
I prefer to apply a threshold when showing Volcano plots, displaying any points with extreme / impossible p-values (e.g. \< 10e-20) with a different symbol at the top of the graph.

However, a better approach is to avoid using p-values as quantitative / rankable results in plots; they're not meant to be used in that way. You can show results as a MA plot instead, plotting log2 fold change vs average expression:

It's a bit trickier to get expression values out of Seurat because they're not currently calculated in the FindMarkers results tables, so you'll need to manually subset the cells and calculate mean expression on a per-marker basis.

---
Tags: r, scrnaseq, seurat
---


# Question
Title: STRINGdb: reverse map STRING identifiers to gene names?
Is there a functionality within STRINGdb package or in another environment to reverse map the STRING identifiers to Gene names?

For example, the entire STRING network can be downloaded using the commands:

```
library(STRINGdb)
string_db <- STRINGdb$new(version="11", species=9606,
                          score_threshold=0, input_directory="" )
full.graph <- string_db$load_all()
netFull
# IGRAPH c07d6e9 UN-- 19354 11759454 --
# + attr: name (v/c), combined_score (e/n)
# + edges from c07d6e9 (vertex names):
#  [1] 9606.ENSP00000000233--9606.ENSP00000272298 9606.ENSP00000000233--9606.ENSP00000253401
#  [3] 9606.ENSP00000000233--9606.ENSP00000401445 9606.ENSP00000000233--9606.ENSP00000418915
#  [5] 9606.ENSP00000000233--9606.ENSP00000327801 9606.ENSP00000000233--9606.ENSP00000466298
#  [7] 9606.ENSP00000000233--9606.ENSP00000232564 9606.ENSP00000000233--9606.ENSP00000393379
#  [9] 9606.ENSP00000000233--9606.ENSP00000371253 9606.ENSP00000000233--9606.ENSP00000373713
# [11] 9606.ENSP00000000233--9606.ENSP00000248929 9606.ENSP00000000233--9606.ENSP00000356737
# [13] 9606.ENSP00000000233--9606.ENSP00000338766 9606.ENSP00000000233--9606.ENSP00000386135
# [15] 9606.ENSP00000000233--9606.ENSP00000319192 9606.ENSP00000000233--9606.ENSP00000432568
# + ... omitted several edges

```

`V(netFull)$name` will provide all the STRING identifiers within this network `netFull`. How to get the corresponding gene names?

For example, `string_db$map()` maps the gene identifiers of the input dataframe to STRING identifiers. It returns the input dataframe with the "STRING\_id" as an additional column. Now, I would like to explore in the reverse direction going from STRING identifiers to gene names.

# Answer
> 2 votes
I don't know if it's too late, but I met same question and got an answer. You can use `add_proteins_description()`. It works in same way with `map()`. If you put a data frame having `STRING_id` column, it adds `preferred_name`, `protein_size` and `annotation`. You may find your answer on `preferred_name` column.

```
require(STRINGdb)
string_db <- STRINGdb$new()
string_db$add_proteins_description(data.frame(STRING_id = [vector containing string ids]))

```

# Answer
> 2 votes
This answer does not rely on STRINGdb package.

STRING provides supplementary identifier mapping files here: https://string-db.org/mapping\_files/.

It looks like you are interested in human STRING name one - https://string-db.org/mapping\_files/STRING\_display\_names/human.name\_2\_string.tsv.gz. You could read it as a separate lookup.

---
Tags: r, bioconductor, networks, stringdb
---


# Question
Title: How do you replace multiple sequence names with new names using python?
I have into a dictionary some sequences. I would like to replace the sequences name for new names contained into a dataframe.

```
#my sequences into a dictionary
sequence = {"seq01": "agatcatggctaactatgcacgtgca",
           "seq02": "atatcatggctaactatgcacgtgca",
           "seq03": "tgatcatggctaactatgcacgtgca"}

#my data frame containing old names and new names
data = {"old_name": ["seq01", "seq02","seq03"], "new_name": ["sample1", "sample2", "sample3"]}
data = pd.DataFrame(data)
print(data)

```

I tried to replace the old names with a for loop, but it was imposible to do it.

Thank you

# Answer
> 4 votes
I am pretty sure there would be more elegant, in-place solutions but here is a quick solution to your problem. Should be OK if your `sequence` dictionary is not way too large.

```
seq_names = zip(data['old_name'], data['new_name']) # brings together old and new names in tuples
renamed_sequence = {}
for i in seq_names:
    renamed_sequence[i[1]] = sequence[i[0]]
renamed_sequence

```

EDIT: The section below is not quite accurate, please see the comments below.

And regarding your comment:

> I tried to replace the old names with a for loop, but it was imposible to do it.

You are not supposed to loop over and at the same time modify an element at the same time anyway (there might be side effects that you would not think of), hence a brand new dictionary in my code.

---
Tags: fasta, python
---


# Question
Title: how can I get expression of an inserted foreign genes?
hi we have transgenic mice with human gene inserted. if we profom rnaseq for the mice, how can we get the expression value of the human gene inserted to the mice? and use this expression to compare with other mice gene expression?

the mice ref genome doesn't have this piece of human gene sequence. anybody knows how?

# Answer
You will need to create (and probably index) a "custom" genome/transcriptome, mouse + whatever human genes inserted, depending on the alignment method of choice. And how to do so would depend on the your tool of choice.
> 1 votes

---
Tags: rna-seq, transcriptome, reference-genome
---


# Question
Title: Print first colum is the condition is true using awk
I have a binary file which looks like this,

```
Orthogroup  A   B   C   D   E   F
OG0000000   1   1   1   0   0   0
OG0000001   1   1   1   0   0   0
OG0000002   1   1   1   0   0   0
OG0000003   1   1   1   0   0   0
OG0000004   0   0   0   1   1   0
OG0000005   0   0   0   1   1   0
OG0000006   0   0   0   1   1   0
OG0000007   0   0   0   1   1   0
OG0000008   1   0   0   0   0   0
OG0000009   1   0   0   0   0   0
OG0000010   1   0   0   0   0   0

```

I want to use awk to print the first column if column D and E is 1 and rest is 0.

```
OG0000004
OG0000005
OG0000006
OG0000007

```

I tried to do this with individual files of each column using the awk command

```
awk 'NR==FNR{arr[$0];next} $0 in arr' file1 file2

```

but doing this I had to generate multiple files and I thought using a binary file will be a more clean way.

Please Help,

# Answer
> 2 votes
I would do something like this:

```
$ cat test.tsv
groupA  1       1       1       0       0       0
groupB  1       1       0       0       0       0
groupC  1       0       1       0       0       0
groupD  1       1       1       0       0       0
groupE  0       0       0       1       1       1
groupF  0       0       0       1       1       0

$ awk '($5==1 && $6==1 && $2+$3+$4+$7==0){print $1}'
groupF

```

# Answer
> 2 votes
Straightforwardly

```
awk '$2==0 && $3==0 && $4==0 && $5==1 && $6==1 && $7==0 {print $1}' bin.file

```

A little spiffier

```
awk '$5 && $6 && !($2 || $3 || $4 || $7) {print $1}' bin.file

```

---
Tags: linux, awk
---


title: Help me in understanding the PDB file

text:
# Question
Title: Help me in understanding the PDB file
Can anybody please explain the below line to me.

"We took the structure and coordinates of nogalamycin from the X-ray structure determined in PDB code = 1D17"

What do I need to download from there....a pdb file or what?

Thank you

# Answer
> 1 votes
well wasnt able to figure out how to make openapigenerator work with rcsb pdb ModelServer https://models.rcsb.org/, I ended up with this code copying the webserver manual input results for ligands:

```
import requests

def main():

  risposta = requests.get('https://models.rcsb.org/v1/1d17/ligand?encoding=cif&copy_all_categories=false&download=true&filename=out',

                          allow_redirects=True)

  open('out.cif', 'wb').write(risposta.content)

  print(risposta.status_code)

if __name__ == '__main__':

   main()

```

unfortunately the server gives you the wrong ligand (try it out)

it says :' /v1/{id}/ligand Coordinates of the first group satisfying the given criteria.'

tried with:

```
risposta = requests.get('https://models.rcsb.org/v1/1d17/ligand?label_atom_id=NGM&encoding=cif&copy_all_categories=false&download=true&filename=out',

                       allow_redirects=True)

```

but no luck again, think problem is in downloaded file :

```
_model_server_params.name     atom_site
_model_server_params.value    '{"label_atom_id":"NGM"}'  ###### here

```

but not an expert of APis and openAPI at all

EDITED:

eventually I got it right using in the code above:

```
risposta = requests.get('https://models.rcsb.org/v1/1d17/ligand?label_comp_id=NGM&encoding=cif&copy_all_categories=false&download=true&filename=out',

                      allow_redirects=True)

```

so I get a out.cif file containing NGM

using another database and FTP :

```
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Nov 19 19:29:30 2021

@author: Pietro

https://bioinformatics.stackexchange.com/questions/18076/help-me-in-understanding-the-pdb-file#comment25289_18076

Ligand extract:

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

"""

from ftplib import FTP

import gzip

import os

from Bio.PDB import PDBParser, PDBIO, Select

def is_het(residue):
    res = residue.id[0]
    return res != " " and res != "W"

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

    def accept_chain(self, chain):
        return chain.id == self.chain.id

    def accept_residue(self, residue):
        """ Recognition of heteroatoms - Remove water molecules """
        return residue == self.residue and is_het(residue)

def extract_ligands(path):
    """ Extraction of the heteroatoms of .pdb files """

    for pfb_file in os.listdir(path):
        i = 1
        if pfb_file.endswith('.pdb') and not pfb_file.startswith("lig_"):
            pdb_code = pfb_file[:-4]
            pdb = PDBParser().get_structure(pdb_code, path + pfb_file)
            io = PDBIO()
            io.set_structure(pdb)
            for model in pdb:
                for chain in model:
                    for residue in chain:
                        if not is_het(residue):
                            continue
                        print(f"saving {chain} {residue}")
                        io.save(f"lig_{pdb_code}_{i}.pdb", ResidueSelect(chain, residue))
                        i += 1

def getFile(ftp, filename):
    try:
        # ftp.retrbinary("RETR", filename ,open(code+'.gz', 'wb').write)
        ftp.retrbinary("RETR " + filename ,open(code+'.gz', 'wb').write)
    except:
        print("Error somewhere !!!!")

code = '1D17'
filename = 'pdb'+code.lower()+'.ent.gz'                 ## pdb1D17.ent.gz

print(filename)

ftp = FTP('ftp.wwpdb.org')   # connect to host, default port

# ftp.login()               # user anonymous, passwd anonymous@
print(ftp.login())

ftp.cwd('/pub/pdb/data/structures/all/pdb/')

print(ftp.pwd()  )

getFile(ftp, filename)

zipped = gzip.open(code+'.gz', 'rb')

unzipped = zipped.read()

# print(unzipped)

saved = open(code+'.pdb', 'wb')

saved.write(unzipped)
saved.close()

path = './'

extract_ligands(path)

```

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

and using FTP from ftplib on http://ftp.wwpdb.org/pub/pdb/data/structures/all/pdb/ : PDB - FTP Archive over HTTP

NOGALAMYCIN is NGM third ligand downloaded (shuold be the only one according to https://www.rcsb.org/structure/1d17)

In any case I am failing at ubderstandung \*,-) how https://models.rcsb.org/

Parameters are related to https://pdb101.rcsb.org/learn/guide-to-understanding-pdb-data/beginner%E2%80%99s-guide-to-pdb-structures-and-the-pdbx-mmcif-format , the mmcif format

but as I told before I can't read inside openapi.json OpenAPI Specification

---
Tags: protein-structure, phylogenetics, pdb, drugs
---


# Question
Title: How to extract part of protein complex separated with HEADER with two molecules from one PDB file
I have a PDB file which contain two molecules (receptor and ligand). Each molecule will have its own header. All in **ONE** PDB file.

The header of receptor section looks like this (line 1-6 of the PDB file):

```
HEADER rec.pdb
REMARK original generated coordinate pdb file
ATOM      1  N   GLY A   1     -51.221 -13.970  37.091  1.00  0.00      RA0  N
ATOM      2  H   GLY A   1     -50.383 -13.584  37.482  1.00  0.00      RA0  H
ATOM      3  CA  GLY A   1     -50.902 -15.071  36.197  1.00  0.00      RA0  C
ATOM      4  C   GLY A   1     -49.525 -15.659  36.443  1.00  0.00      RA0  C

```

and this ligand section (line 11435 to 11440) of the PDB file

```
HEADER lig.000.00.pdb
ATOM      1  N   MET A   1      27.318 -26.957  12.663  1.00  0.00      LA0  N
ATOM      2  H   MET A   1      27.313 -27.570  11.870  1.00  0.00      LA0  H
ATOM      3  CA  MET A   1      28.374 -27.102  13.668  1.00  0.00      LA0  C
ATOM      4  CB  MET A   1      28.531 -28.564  14.090  1.00  0.00      LA0  C
ATOM      5  CG  MET A   1      27.224 -29.154  14.628  1.00  0.00      LA0  C

```

Note that the receptor and ligand section also contain the string `RA0` and `LA0` at the 11th column of PDB file.

What I want to do is to rename the chain of receptor as chain A and ligand as chain B. In order to do that I intended to extract both part into two different objects first, then rename the chain and finally put them together again.

I tried this code with Bio3D. But doesn't work:

```
library(bio3d)
pdb_infile <- "myfile.pdb"

pdb <- read.pdb(pdb_infile)
receptor_segment.sele <- atom.select(pdb, segid = "RA0",  verbose = TRUE)
receptor_pdb <- trim.pdb(pdb, receptor_segment.sele)
ligand_segment.sele <- atom.select(pdb, segid = "LA0",  verbose = TRUE)
ligand_pdb <- trim.pdb(pdb, ligand_segment.sele) # showed no entry

```

What's the way to do it? I'm open to solution with R or Python.

# Answer
> 2 votes
Since this is not a correct PDB file, it can be interpreted differently by different PDB reading libraries. The most reliable way to change it would be to read it as plain text and change the chain ID to `B` in the lines that have segment ID `LA0`.

Note: Strings RA0 and LA0 start at column 73. Columns in the PDB format correspond to characters. This is not a whitespace-delimited file.

Chain ID is in column 22.
Segment ID is in columns 73-76.

As a side note, segment ID was only briefly part of the official PDB spec (around 1996), nevertheless it's supported by many programs.

Of course, if you have only one file that you want to modify, the fastest way is to do it manually in a text editor.

---
Tags: r, python, protein-structure, pdb
---


# Question
Title: Why does BWA MEM orientation contradict my library prep method
I have some RNA-seq data from a stranded paired end library prep, with dUTP and UDG preparation, so the orientation should be RF (confirmed with sequencing provider). I assembled the reads with Trinity, specifying RF for the strandedness, and I am mapping the reads to the transcripts with bwa-mem (I am aware bowtie2 or STAR might be better for this, but bwa is what I presently have installed, so just using it to check).

In the bwa output, I see lines `# candidate unique pairs for (FF, FR, RF, RR): (0, 5372, 2, 0)`, which seems to contradict the expected orientation. Is there a reason that this is expected, or does it indicate a problem?

# Answer
> 2 votes
Not familiar with the protocol, but I guess RF means different things here. If I am right, with your RF protocol, the first read comes from the reverse strand of the transcript and the second read from the forward strand, but the pair orientation is still FR:

```
------------------->  transcript
  R2 ---->   <---- R1

```

The RF pair orientation should be

```
<----   ---->

```

You can confirm this by checking the mapping strand of read1.

---
Tags: rna-seq, sequence-alignment, bwa
---


title: HG19 Position meaning

text:
# Question
Title: HG19 Position meaning
I am looking at genomic data (HG19). I have many SNPs, their chromosomes and positions. I want to look at certain SNP (suppose its chromosome is 1 and position 77,226,919), and extract all SNPs in distance 1mb from that SNP.

Are the HG19 position in terms of base-pairs? Meaning that I need to take all SNPs in position $77,226,919 \pm100,000$? Or is it something else?

# Answer
Yes, they are base pair positions
> 2 votes

# Answer
HG19 refers to the *build* used as the reference genome. To extend on what @Emily\_Ensembl said, base pair positions may change between different builds.

Please look at the following link for a more detailed explanation: https://gatk.broadinstitute.org/hc/en-us/articles/360035890951-Human-genome-reference-builds-GRCh38-or-hg38-b37-hg19
> 1 votes

---
Tags: snp, gwas, hg19
---


# Question
Title: What does 'Human Alternative sequence Gene' mean in Ensembl?
I have downloaded a dataset containing RNA-Seq RPKM (reads per kilobase per million) values for 52376 genes (in humans). To find RPKM values for the gene KIR2DL2 in the dataset, I have to search for it using its Ensembl ID.

After searching for this gene using the Ensembl browser, I obtain many results.

For all the results in the first page, it says *Human Alternative sequence Gene* after the gene name. When I have searched for other genes using their Ensembl ID I get only one result:

What does *Human Alternative sequence Gene* mean in Ensembl? I am not sure how to tell which is the correct Ensembl ID for the KIR2DL2 gene.

Thanks.

# Answer
> 4 votes
These are genes on the patches and haplotypes. Patches are repairs in the genome, where they have discovered the sequence or assembly is incorrect, so there is an alternative sequence that is created on top of the reference assembly (the reason they don't change the reference is that then there's a new reference and this causes a massive headache for everyone). Haplotypes are alternative versions of genomic loci that exist in some individuals, such as the MHC locus, which are different and common enough to require a proper assembly change, more than just genetic variants.

The genes that are found on the alt sequences are listed as *Human Alternative sequence Gene*. They are both real genes and correct for the gene name.

If the gene is on a patch (the location will include the word PATCH in its name) you should probably use the patch gene because the primary assembly gene may be dodgy.

It gets more complicated if you're working with haplotypes, as all the versions are correct. Most people take the easy option of using just the primary assembly, but if you're variant calling, it may be worth trying to determine which haplotype the individual has and calling against that. This is particularly true for highly variable regions like MHC, where there are genes only found on the alternative haplotypes and not found on the primary assembly.

There's more information about alternative sequences from Ensembl here. And here from the GRC.

---
Tags: sequence-annotation, ensembl, human-genome
---


title: replacing periods with hyphen

text:
# Question
Title: replacing periods with hyphen
I am a software engineer and not a bioinformatician. I am looking at some code where they are looking at aligned proteins that are in the A2M alignment format (https://compbio.soe.ucsc.edu/a2m-desc.html).

According to the page, the

```
Uppercase characters and "-" represent alignment columns,
and there must be exactly the same number of alignment columns in each sequence

```

Now, the code I am looking at is doing something like:

```
sequence = sequence.replace(".","-")

```

So replacing the dots with hyphen and I am not sure what the purpose of that might be.

For anyone interested, the code is here (See line 182): https://github.com/debbiemarkslab/DeepSequence/blob/master/DeepSequence/helper.py

# Answer
In A2M format, upper case letters represent matches, lower case letters represent inserts, dashes represent deletions, and dots (or spaces outside the identifier lines) represent gaps aligned to inserts. So, both "-" and "." are essentially gaps, but assumed to have different origin. This information is supplementary, and most MSA programs ignore dots and spaces or treat them incorrectly.

The code you're looking at, in lines 187–197, removes sequences that contain characters not from `self.alphabet` (defined in lines 74–86) by checking `if letter not in alphabet_set and letter != "-"` (line 192). Therefore, just for the sake of parsimony, it looks natural to transform "." into "-" beforehand.
> 2 votes

---
Tags: sequence-alignment, data-preprocessing, sequence-annotation
---


# Question
Title: At what point is a gene different between species that we call it a different gene?
I am new to genetics and I know humans share many genes with mice for instance but that there are slight differences in conserved nucleotide sequences. Is there a community consensus around at what point a gene in one species is a different gene? How are regulatory elements like enhancers factored into this definition?

# Answer
From an evolutionary angle, you are describing an ortholog: "Orthologs are genes in different species that evolved from a common ancestral gene by speciation, and, in general, orthologs retain the same function during the course of evolution." quoting here. A lot of effort has been made to construct ortholog databases (e.g. see quest for orthologs). They are distinguished from paralogs (approximating duplicated genes) and lateral gene transfer events, there is however, no magic cutoffs for defining orthology.

Your 2nd point about regulatory elements, well they could fit in the same framework, any sequence can be clustered. What is a gene anyway? Is a noncoding RNA a gene? How do we go from a collection of transcripts to this unit we term gene? There is some ambiguity...
> 3 votes

---
Tags: phylogenetics, genomics
---


# Question
Title: Generating PDB file for custom DNA sequence
I asked this question in biology.stack exchange and was referred to this community. I am going to request migration for this question. https://biology.stackexchange.com/q/105024/67415

I want a PDB file of DNA sequence I am interested in (10 base pairs) for docking. I will appreciate if any of you can please clarify some of my doubts.

Do you think, I can generate PDB file by using this link http://www.scfbio-iitd.res.in/software/drugdesign/bdna.jsp

How can I get the pdb file of the DNA sequence I want from protein data bank, I did search it but did not find the one I want. Is it a way to extract the PDB file?

What are other methods I can use to generate PDB file of B-DNA sequence?

Thank you for your time.

Please delete if it is against the group for posting again and redirect me to where I can request migration of this question.

# Answer
> 1 votes
I normally use PyMOL. In PyMOL you can generate a model of a sequence of DNA with the command fnab — a sister command to the better known fab, which makes a protein sequence of a given secondary structure. I've not encountered that site, but it says it uses experimental data, which is a better approach.

## Side note 1

OP is asking for DNA, but if anyone stumbles across this, do note that RNA is a totally different ballgame and requires more advanced solutions (example).

## Side note 2

I should say docking a known sequence is fine, but not docking to guess the sequence is not.

For docking, if you have a zinc finger, helix-loop-helix, leucine zipper, MADS box or similar domain which binds DNA in a consistent manner then it is straightforward. However, predicting the bound sequence is far from trivial or accurate. For example, in MEF2C the N-terminal amine there are three reader-head dictating specificity and one is a floppy terminus, which most forcefield-based methods would struggle to get right.

## Testing

So let's say we want to test PyMOL's fnab in the GUI. We could pick a protein, MEF2A example from above:

```
fetch 1EGW, type=pdb1
print cmd.get_fastastr('chain E') # to get the sequence
fnab AAAGCTATTATTAGCTT
align obj01 and chain A, 1EGW and chain E
show sticks, chain E+F
hide cartoon
remove solvent

```

The RMSD is 1.6 Å which is not great, but pretty decent considering that it is with a generated sequence that does not take into account backbone distortion from the bound protein.

# Answer
> 1 votes
There is a tool for making DNA segments of any curvature. It exports as a PDB, and the people behind the tool are quite committed and helpful.

https://www.lifeexplorer.info/

Here is a curved DNA segment I generated using their tool; they can make much more dramatic bends.

I regularly use the site you cited (scfbio...) to generate idealized straight DNA segments for illustration purposes. It is not sequence-dependent. For an actual solved structure you can extract specific chains from existing structures containing DNA (e.g. RNA polymerase II structures by the Cramer group). These structures will have their own sequence, not your specific one, and will be distorted by protein-DNA interactions. Other than that, there likely won't be a PDB structure of your specific sequence. And even if there was, the conditions likely would differ from the ones your structure was acquired at. I may be entirely off, but it sounds like you should be fitting individual bases, sugars and phosphates, not continuous poly-d-nucleotides.

Shown below is the pymol visualization of a 100bp dsDNA input for B-DNA from the resource you cited It does not take into account the DNA flexibility (persistence length etc)

---
Tags: phylogenetics, pdb, docking
---


title: Trinity Multi-species compatibility

text:
# Question
Title: Trinity Multi-species compatibility
I am using Trinity from last couple-up months or so. Initially, using same species in different condition was worked fine for me. Now i am trying to do a DE analysis using multiple species. Actually, i have 3 fish performs similar behaviour and another 3 performs similar strategy. Here i considered 3 similar species as replicate for one type (and totally i have two type) . **Though these are 6 different species but belong to same family**.

I have already made a combine assembly according to trinity protocol (--samples\_file). After assembly i have checked the separate alignment rate for each species and it was above 85%. And lastly i checked heat-map and it gave me a result that somehow match with my expectation.

But i am a bit confused about the possible generation of "assembly chimera". In most of the cases i found people do DE analysis with same species. Also didn't find any reference article that did multiple species DE analysis using Trinity pipeline. So it will be a great help for me if i get some suggestion that either i should go with existing result or i need to find another way.

# Answer
Given that you're using multiple species / samples, the results from this are going to be quite weird.

For doing such a comparison, I think the best approach is to do a combined assembly (as it looks like you've done) to generate a single chimeric transcriptome (which will not be a good representation of any of the component species), then do a differential expression analysis using that combined assembly as the reference.

Ideally you would have multiple replicates for each different species, and then add the species and conditions into the statistical model, so that any species-specific variation can be identified and/or filtered out.

Six separate assemblies that are later merged will be less useful for a differential expression, because it will bias towards differences in the assemblies. This is a similar rationale to that used for generating UMAPs from multiple single-cell datasets where a combined UMAP is generated for the full dataset with identical parameters, so that sample-specific differences in the UMAP parameters are less of an issue.

I think it will be okay to go ahead with the combined assembly, bearing in mind that the lack of replicates for the differential expression study is going to make it more difficult to identify biological significance. Using split assemblies will not make that better; if anything, it'll make it worse.
> 2 votes

---
Tags: trinity
---


title: Visualizing top expressed genes

text:
# Question
Title: Visualizing top expressed genes
I have a spreadsheet of CPM values

I want to visualize the top 20 genes expressed across my samples

This is how my data look like

```
> head(b[1:5,1:5])
                 30566-001 - CPM 30566-002 - CPM 30566-003 - CPM 30566-004 - CPM 30566-005 - CPM
hsa-miR-548ap-3p               0               0               0               0               0
hsa-miR-548t-3p                0               0               0               0               0
hsa-miR-548e-3p                0               0               0               0               0
hsa-miR-548az-3p               0               0               0               0               0
hsa-miR-548f-3p                0               0               0               0               0
> dim(b)
[1] 2632   42
>

```

How I can do that?

# Answer
Gene selection is entirely on you and on the question you ask. Most expressed genes would be selected by `rowMeans` or `rowMedians` or even `rowSums` but this is not very informative as digital expression levels (=counts) is influenced by many technical factors and there is (to my knowledge) no direct link to biological relevance.

Alternatively, you can select by genes that vary a lot between experimental groups. Naively that could be the rowwise variance (`rowVars`) or standard deviation of the logCPMs (use logscale, otherwise variance is a function of the expression level) or the output of a differential analysis. Ideally plug the raw counts into DESeq2 or similar frameworks. If you only have CPMs then you probably want to feed the logCPMs into limma-trend pipeline, see for a discussion on ways to get going with normalized counts for differential expression in this Bioconductor thread.

Visualization is then again on you. Either make a heatmap (see for example ComplexHeatmap package) or aggregate the genes into a PCA (see for example PCAtools package at Bioconductor), all mentioned packages contain extensive example code and background information which is fully sufficient to run such an analysis.
> 2 votes

---
Tags: r, differential-expression, deseq2
---


title: DEseq2 contrast in two time points

text:
# Question
Title: DEseq2 contrast in two time points
I need help with some clarifications please

I have control versus treatment in two time points like

```
> y
                             X condition time
1       CTRL_24_hrs_replicate1   control   24
2       CTRL_24_hrs_replicate2   control   24
3       CTRL_24_hrs_replicate3   control   24
4  treatment_24_hrs_replicate1        t   24
5  treatment_24_hrs_replicate2        t   24
6  treatment_24_hrs_replicate3        t   24
7       CTRL_48_hrs_replicate1   control   48
8       CTRL_48_hrs_replicate2   control   48
9       CTRL_48_hrs_replicate3   control   48
10 treatment_48_hrs_replicate1        t   48
11 treatment_48_hrs_replicate2        t   48
12 treatment_48_hrs_replicate3        t   48
>

```

I want to see the difference of treatment versus control from 24 hours to 48 hours, so this is not just the difference of treatment versus control or difference between 24 versus 48 hours rather the treatment versus control from 24 hours to 48 hours

This is what I have done so far

```
>  dds <- DESeqDataSetFromMatrix(countData=count,
+                                 colData=design,
+                                 design=~condition, tidy = FALSE)
Warning message:
In DESeqDataSet(se, design = design, ignoreRank) :
  some variables in design formula are characters, converting to factors
>   dds$group <- factor(paste0(dds$time, dds$condition))
>   design(dds) <- ~ group
>   dds <- DESeq(dds)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
>   resultsNames(dds)
"Intercept"                    "group_24IT_vs_24control"         "group_48control_vs_24control" "group_48IT_vs_24control"
>

```

How to find the difference at time point 24 vs 48 in condition ?

I searched and I found a similar thing but I don't know how to adopt that to my data

```
results(dds, contrast=list(c("time_4_vs_0", "conditionLD.time4")), test="Wald")

```

I tried this and I got error

```
> results(dds, contrast=list(c("time_48_vs_24", "conditionIT.time24")), test="Wald")
Error in checkContrast(contrast, resNames) :
  all elements of the contrast as a list of length 2 should be elements of 'resultsNames(object)'
>

```

Can you help me in adopting this to my data please?

# Answer
> 1 votes
You can't figure out how to adapt this to your data?

http://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#contrasts

Ignore resultsNames. You don't need them for simple comparisons like you want to do.

# Answer
> 1 votes
you can do:

```
data$condition_new <- as.factor(paste0(data$condition, " - ", data$time))
dds <- DESeqDataSetFromMatrix(countData=count,
                             colData=design,
                             design=~condition_new
                             , tidy = FALSE)
enter code here
results(dds, contrast=c("condition_new", "control - 24", "condition - 4"))

```

I hope it helps.

---
Tags: deseq2
---


# Question
Title: Where can I get info about A. thaliana proteins' locations?
I have a list of *A. thaliana* proteins and I'd like to know their location in the cell. I've found SUBA but they seem not to share their dataset and only allow queries through their website. It's inconvenient for me, because I need to locate 1.6k of proteins and process the results further. Is there a dataset or an API I could use instead?

# Answer
> 1 votes
Given the fact that it is a simple query, in my opinion, the best approach would be an advanced search on AmigoGO. There is no need to bother with complex methods like programmatic access APIs, as you basically have only two criteria to meet: organism and cellular localization. After you specify those, you can download the file using the appropriate button, or mess a bit with the columns that you are interested in. Once you have the file, you can then manipulate it further to your specific need with any tool of your preference, for example, to retain only proteins of your interest. Hope this helps.

---
Tags: proteins, database
---


# Question
Title: Changing the order of colors in pheatmap
I plotted a heatmap

I want cancer to be red and healthy to be blue

```
   pheatmap(heat1, annotation_col=df,color=colorRampPalette(c("navy", "white", "red"))(50))

```

You are seeing cancer is in blue in annotation, how I can change that

# Answer
> 5 votes
It is possible with `annotation_colors` argument :

First create a list with your conditions and then add the argument :

```
annot_colors=list(HealthStatus=c(Cancer="#F0978D",Healthy="#63D7DE"))

pheatmap(heat1,
         annotation_col=df,
         color=colorRampPalette(c("navy", "white", "red"))(50),
         annotation_colors=annot_colors)

```

---
Tags: r, heatmap
---


# Question
Title: Updating the GFF3 + Fasta to GeneBank code
I'm trying to convert gff3 and fasta into a gbk file for usage in Mauve. I've found a solution but the code is outdated:

```
"""Convert a GFF and associated FASTA file into GenBank format.

Usage:
gff_to_genbank.py <GFF annotation file> <FASTA sequence file>
"""
import sys
import os

from Bio import SeqIO
from Bio.Alphabet import generic_dna
from BCBio import GFF

def main(gff_file, fasta_file):
    out_file = "%s.gb" % os.path.splitext(gff_file)[0]
    fasta_input = SeqIO.to_dict(SeqIO.parse(fasta_file, "fasta", generic_dna))
    gff_iter = GFF.parse(gff_file, fasta_input)
    SeqIO.write(gff_iter, out_file, "genbank")

if __name__ == "__main__":
    main(*sys.argv[1:])

```

But the Bio.Alphabet module is removed now. I tried to read the documentation on these changes and how to update the code (https://biopython.org/wiki/Alphabet) but I'm still confused about how to update the code. If I remove the Bio.Alphabet and the generic\_dna, the code will return ValueError

```
ValueError: missing molecule_type in annotations

```

I know that I have to type the molecule\_type = "DNA" but I'm confused as to where should I put it because the code directly inputs the fasta\_file into the SeqIO.parse() function.

Anyone knows how to update code related to SeqIO.parse() with the removal of Bio.Alphabet? Or anyone knows how to convert GFF3 + Fasta into gbk files?

Note: I can't use the seqret because my files is just too big for it, same problem as this question.

Thanks in advance!

# Answer
You can set the molecule type on the existing annotations attribute of the record objects just before they get to `SeqIO.write`.

There might be an even leaner way to do it, but if open the output file (so that you can do multiple writes to one file handle) and then iterate over each record from the GFF iterator, you can do it like this:

```
with open(out_file, "wt") as f_out:
    for record in GFF.parse(gff_file, fasta_input):
        record.annotations["molecule_type"] = "DNA"
        SeqIO.write(record, f_out, "genbank")

```

(The documentation is vague on that point but if you check you'll see that `record.annotations` is a dictionary that you can modify to set the molecule type.)
> 4 votes

---
Tags: python, parsing, seqio
---


# Question
Title: Publicly available genotype / phenotype dataset?
I'm looking for publicly available genotype data (eg. 1000 Genomes or HapMap) that also has associated phenotypes (any traits). I'm wanting to use it as a test run for a GWAS pipeline. Does anyone know if this exists anywhere?

# Answer
See openSNP here:

https://opensnp.org/

It is a web site for people to upload their own genetic datasets, providing freely-accessible public data with annotated phenotypes.

> openSNP allows customers of direct-to-customer genetic tests to publish their test results, find others with similar genetic variations, learn more about their results by getting the latest primary literature on their variations, and help scientists find new associations.

The large collaborative research genetic data projects typically have very strict access constraints due to the realisation that genetic data is personally-identifiable data, and the expectation of anonymity. That means there's a low likelihood that a large research project will make their phenotype + genotype data publicly available.
> 1 votes

---
Tags: gwas, plink
---


title: Remove gaps from alignment?

text:
# Question
Title: Remove gaps from alignment?
I have an MSA (protein sequence) and I have used various programmes (Clustal, Aliview, MEGA11 etc) to align. However, in all programmes I get many gaps which is not ideal as I am trying to construct a tree from this alignment.

My question is are there any tools to delete gaps in alignments?

Your help would be greatly appreciated.

# Answer
You may also find useful the BMGE tool (stands for Block Mapping and Gathering with Entropy). BMGE selects regions in a multiple sequence alignment that are suitable for phylogenetic inference by computing their entropy-like score. This allows to distinguish among biologically expected and unexpected variability for each aligned character. BMGE also implements different types of trimming (for instance, deletes gaps) to minimize phylogeny reconstruction artifacts.

Available from here and also installable from Bioconda.
> 3 votes

# Answer
If you want an ungapped alignment, you could *do* an ungapped alignment (the older read aligners were ungapped, e.g. bowtie). However, I suggest that you don't want an ungapped alignment, you just want a stricter one, here's a few options...

* Remove/trim poorly aligning sequences. You could try GUIDANCE, which has an explicit model or GBLOCKS for a more simple filtering criteria
* Tune your alignment parameters for ClustalOmega/MAFFT etc. to make it stricter
* Finally, consider that your sequences must be homologous otherwise you'll trim away sequences until you have a random alignment that appears good
> 2 votes

# Answer
Gblocks is a tool for removing gaps/poorly conserved positions from alignments. Available here. Here is a webserver that runs it for you.

It's a bit old, possibly there are others at this point.
> 1 votes

# Answer
This new tool, CIAlign, talks about some automated removal of gaps in multiple sequence alignments

https://twitter.com/thekatybrown1/status/1306215132670889984

https://github.com/KatyBrown/CIAlign
> 1 votes

---
Tags: sequence-alignment, multiple-sequence-alignment
---


title: bowtie No alignments

text:
# Question
Title: bowtie No alignments
*This question was also asked on Biostars*

I am working on miRNA alignment using bowtie. But it seems to fail to map. Here are the commands:

```
bowtie-build -f miRNA.fa miRNA
bowtie --norc -k 1 -v 0 /public/home/wschen_gibh/cpa-seq/ref/miRNA/miRNA -1 SRR16643593_1.NEW.trim.fq -2 SRR16643593_2.NEW.trim.fq -S test.sam

```

and here is the result:

and the output sam file:

# Answer
You should be using Bowtie 2 instead of bowtie:

> Bowtie 1 was released in 2009 and was geared toward aligning the relatively short sequencing reads (up to 25-50 nucleotides) prevalent at the time. Since then, technology has improved both sequencing throughput (more nucleotides produced per sequencer per day) and read length (more nucleotides per read).

http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml#how-is-bowtie-2-different-from-bowtie-1

In particular, Bowtie 2 supports local alignment, which is *essential* when mapping microRNA. All microRNA sequences are short, which means that there will be adapter read-through in the sequenced reads. This adapter sequence will not map to the genome, and that will cause global aligners (like Bowtie) to fail.

An alternative approach would be to trim off palindromic adapter sequences prior to aligning (which it is possible you've tried to do given the file name). Because of adapter read through, standard adapter trimming will not be sufficient. What needs to be trimmed off is the *reverse-complement* of the adapter sequences. I have no experience using *fastp*, so don't know if it is capable of doing that.

However, even after trimming there will still be some additional short sequences left on some reads, so a local alignment would still be a good idea.

Doing a local alignment with Bowtie 2 is much easier and faster than trying to troubleshoot adapter read-through issues. The command line you would run is similar, but not identical:

```
bowtie2-build -f miRNA.fa miRNA
bowtie2 --local -x miRNA -1 SRR16643593_1.NEW.trim.fq -2 SRR16643593_2.NEW.trim.fq | samtools sort > test.bam

```
> 2 votes

---
Tags: sequence-alignment
---


# Question
Title: How should I downsample and normalize R1s and R2s and incorporate this into Lexogen's QSPA tool
I am using Lexogen's Quantitative Sequencing Pool Analysis tool (here: https://github.com/Lexogen-Tools/quantseqpool\_analysis) to analyze R1 and R2 files.

I have been able to successfully run this tool. However, I need to downsample and then normalize these results so I need to do this all again.

To downsample, I put the R1s (the ones generated after idemultiplexing) through the SeqTK\_sample tool on UseGalaxy. Afterwards, I ran the QSPA tool with the new R1s and old R2s. I skipped the idemultiplexing part by # out the code. This is generating errors (screenshots are attached). I don't think these errors are from me # out the code but rather from the new R1s with the old R2s. I didn't receive any errors the first time I ran this tool.

I will be normalizing my data with DESeq2 using the summary\_unique.tsv file.

If anyone has any tips/suggestions it would be much appreciated. I am sure there is a better way to do this. If you require anything from me just ask!

Cheers,

# Answer
The issue was that we downsampled with useGalaxy's SeqTK tool prior to running UMI-tools and cutadapt (the 2nd and 3rd process, respectively, of Quantitative Sequencing Pool Analysis).

I tried downsampling with the R1 data, available after cutadapt (located in the "trimmed" folder), and everything ran fine. I believe the issues were caused by the barcode sequences being removed by the SeqTK tool prior to running UMI-tools.

You can (somewhat) customize the QSPA script by #'ing out parts of the script that you don't want to run. In this case, I #'d out Star Aligner, UMI-tools (deduplication), and Feature Counts so that the script would stop after cutadapt. After I downsampled with SeqTK, I ran the script using only the parts that I had previously #'d out.
> 0 votes

---
Tags: sequence-analysis
---


# Question
Title: Python script to simultaneously generate multiple pdbqt files for AutoDockTools?
I have a folder of pdb ligand structures that I'd like to test in some docking experiments with a certain protein using AutoDock Vina. I am not familiar with Python so I've always been using the interactive procedure on AutoDockTools to prepare the pdbqt ligand files instead of the prepare\_ligand4.py script. Could you please help me understand how to iterate this python script in a way that it can be applied to a set of pdb ligand files in order to create the corresponding set of pdbqt ligand files? Thank you!

# Answer
> 1 votes
Answer from @matteo-ferla, converted from comment:

If command line interface is more suitable for you, you may consider OpenBabel, which can predict Gasteiger charges, protonate at pH 7, accept a smi file of multiple SMILES or a multi-entry sdf file and can output in PDQBT format, but is not quite the same as the python script allows merging and various settings. But to use the script iteratively (a horrendous thing to do as it loads up the python kernel each time) check out how to loop files via bash

---
Tags: python, pdb, docking, structural-biology
---


# Question
Title: Best way to dock subset of protein PDB and definition of a ligand
I have extracted residue 45 to 88 from a PDB file (**2YRQ**). With that I have a new smaller PDB file. Let's call it 45\_88.pdb And I want to dock it to some other full protein.

My question is what's the best way to dock it. I'm considering **Autodock Vina**. But I'm not sure if it's appropriate since it is for ligand-protein docking.

What's the definition of a ligand? Can a protein with only 42 residues be considered a ligand?

# Answer
Answer from @pippo1980, converted from comment:

https://pubs.acs.org/doi/10.1021/acs.jctc.9b01208 Comprehensive Evaluation of Fourteen Docking Programs on Protein–Peptide Complexes
> 1 votes

---
Tags: pdb, docking, 3d-structure, structural-biology
---


# Question
Title: BioPython - Retrieve sequence records from pubmed database
I've been assigned a task to fetch sequence records from "**pubmed**" database:

> fetch 30 record of type "fasta" from (pubmed database, with term hemoglobin AND alpha) as a sequence records

It doesn't make any sense to me!

I know the PubMed database is for publications and citations, but how can I retrieve those records from it? is that even possible?

I tried to search for the IDs with a given term in it and then retrieve the records from another database (like "**nuccore**"), but I don't know if this is right.

Any ideas?

# Answer
Answer from @devon-ryan, converted from comment:

They meant the nucleotide database.
> 1 votes

---
Tags: python, homework
---