# Question
Title: When using cutadapt, can I specify an R2 adapter as optional when I have a required R1 adapter?
I'm using cutadapt 3.5 to trim adapters and perform some filtering on paired-end data. Both R1 and R2 sequences have 3' adapters that *might* be found depending on the sequence length, but R1 also has an anchored 5' adapter that I want to require for the pair to be kept in the output. For example, if that required R1 5' adapter is AAAAAAAA, the R1 3' adapter is TTTTTTTT, and the R2 3' adapter is GGGGGGGG, and I have this input:

R1:

```
@read1 should be kept
AAAAAAAACGTCCTGGGATTTGTAATAATATTTTAGTTCTGAGCGACAAGTAAGGGATAATTTTTTTT
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
@read2 should be filtered
TACCAACTACATTTAGCTTCAGGCTAGTGATGCCCGCCGTCGGCACACTGGACACATGGCTTTTTTTT
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
@read3 should be kept
AAAAAAAACAATTTTTACTCTAGAAATGTCTGTGCTCATCACCCCGACACCGAATAGCTATGACTCAC
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF

```

and R2:

```
@read1
AGGCTCACCGCCACTGTTGTACCTTCCTATCGCCACTCTAAGATCTATGTAACCTCTCCCGGGGGGGG
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
@read2
TGAATAATACTGCGGTTTAATCGATATATTCGGGATTATGCAAGACACCCTACGTACTTAGGGGGGGG
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
@read3
TACTGGAGGGTACAAGCCCGGACTATCCACGGCGTCAGCTGGCTTAGCATTGAAGGTCGGCGGTGTGT
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF

```

I'm aiming to get this R1 out:

```
@read1 should be kept
CGTCCTGGGATTTGTAATAATATTTTAGTTCTGAGCGACAAGTAAGGGATAA
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
@read3 should be kept
CAATTTTTACTCTAGAAATGTCTGTGCTCATCACCCCGACACCGAATAGCTATGACTCAC
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF

```

And this R2 out:

```
@read1
AGGCTCACCGCCACTGTTGTACCTTCCTATCGCCACTCTAAGATCTATGTAACCTCTCCC
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
@read3
TACTGGAGGGTACAAGCCCGGACTATCCACGGCGTCAGCTGGCTTAGCATTGAAGGTCGGCGGTGTGT
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF

```

This approach doesn't work as it discards read3:

```
cutadapt -a '^AAAAAAAA...TTTTTTTT' -A 'GGGGGGGG' example.R{1,2}.fastq --discard-untrimmed -o out.R1.fastq -p out.R2.fastq

```

I expect that's because the required/optional behavior only applies to the linked adapter for R1, and R2 is considered untrimmed for read3, which then filters the whole pair because of the default `--pair-filter=any` behavior (if *either* read is untrimmed, the pair is discarded). I found I can get the filtering I want with `--pair-filter first`, but then the other filtering criteria like quality and length will only apply to R1, too.

My question is, can I get the trimmed/untrimmed behavior I'm looking for (read1 and read3 trimmed, read2 untrimmed) while also applying cutadapt's filtering to both reads in each pair? In other words, can I specify to cutadapt that R1's anchored adapter is required and the other two optional, without resorting to changing the pair filter setting?

# Answer
> 0 votes
This is not possible in a single cutadapt call, but can be done fairly seamlessly with a pipe between two calls and with cutadapt's interleaving feature. (Thanks to the author for help with this.)

---

Giving the R2 3' adapter to one cutadapt command and interleaving the read pairs to piped output, a second cutadapt command can take that interleaved output and use both of R1's adapters as a linked adapter:

```
cutadapt -A GGGGGGGG --interleaved example.R{1,2}.fastq | \
cutadapt -a "^AAAAAAAA...TTTTTTTT" --interleaved --discard-untrimmed -o out.R1.fastq -p out.R2.fastq -

```

Now I can add other filters to the second call.

An additional read pair that should be filtered:

R1:

```
@read4 would be kept but R2 is too short
AAAAAAAACAATTTTTACTCTAGAAATGTCTGTGCTCATCACCCCGACACCGAATAGCTATGACTCAC
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF

```

R2:

```
@read4
CTACTAGACATATACCCTGGGGGGGGTTCCCCGTGAAGACGGTTACGATGTGGACAAGCACTCCAGCG
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF

```

I can now give a minimum length that applies to both in each pair:

```
cutadapt -A GGGGGGGG --interleaved example.R{1,2}.fastq | \
cutadapt -a "^AAAAAAAA...TTTTTTTT" --interleaved --discard-untrimmed --minimum-length 30 -o out.R1.fastq -p out.R2.fastq -

```

With that I get read1 and read3 in the output, as intended.

---
Tags: ngs, fastq, trimming, cutadapt
---


# Question
Title: How to identify Father-Daughter Relationship from Phylogenetic analysis?
I have performed a phylogenetic analysis using CLUSTALW and obtained a phylogenetic tree. I know how to read a phylogenetic tree and how to relate the sequence similarity between two sequences by looking at the tree, but what I am concerned about is how one can decipher the father-daughter relationship by looking at the tree? Also, if I want to isolate the sequence of the father corresponding to a daughter sequence how can I do that? I have attached the corresponding image. In the image, I know that sequences "QXN18196" and "QXN18436" are closely related. But will it be safe to say that sequence "QXN18520" is the father sequence for both the sequences? Or if I want to establish a father-daughter relationship for the whole tree, what will be the relation between the sequences?

# Answer
> 1 votes
Unless you have wildly differing mutation rates on different branches, you can't say that QXN18520 has the same sequence as the most recent common ancestor of QXN18196 and QXN18436, because QXN18520 was collected today, not back whenever QXN18196 and QXN18436 shared a common ancestor. All that you can say from this phylogeny is that there probably existed a common ancestor of QXN18196 and QXN18436, marked here by 5e-9\[86.8\], and that this common ancestor and QXN18520 shared a common ancestor even further back, marked here by 5e-9\[0\].

(General observation: Note also that for a number of reasons, phylogenies often don't reflect reality -- they are simply a maximum likelihood estimate of the actual ancestral pedigree, often under some very unrealistic assumptions, and/or according to an incomplete or incorrect view of reality -- e.g. the samples that you have collected may have a lot of allele matches due to identity by state (IBS) rather than identity by descent (IBD). In other words, let phylogenies inform you, but don't assume they are the absolute truth.)

---
Tags: phylogenetics
---


# Question
Title: How to de novo hybrid assemble with Pacbio CCS and Illumina PE reads
I would like to perform *de novo* genome assembly on a diploid microalgal strain.

I have two datasets:

* PacBio CCS/HiFi reads, low coverage.
* Illumina PE 2x150 (standard shotgun)

Does anybody have any software recommendations for how best to use both in a hybrid *de novo* assembly?

Most options out there are set up for Pacbio CLR reads, not Pacbio CCS/HiFi, for using in hybrid assembly. Perhaps I should assemble short read contigs first then combine with wengan, or go back to the PacBio subreads instead of CCS to use in a hybrid method?

# Answer
> 1 votes
If you have reasonable coverage with HiFi (30X per haplotype), then there usually isn't any reason to also use Illumina. My understanding is that Illumina can make the assembly worse if you try to polish with it, because it can erroneously induce haplotype switching errors.

I would use HiCanu or HiFiAsm, depending on whether you expect significant heterozygosity, and also assuming that you have non-trivial HiFi coverage.

If your HiFi coverage is very low (\<15X per haplotype, maybe), you might consider going back to subreads, assembling those, and polishing with Pilon and your Illumina reads, or one of the workflows that you have found. But otherwise I don't think the Illumina will be very useful. If you had Hi-C reads you could use them directly in HiFiAsm to assist phasing and scaffolding, but I'm assuming that your Illumina is shotgun.

I know it's a PB promo material, but you see that PB promotes HiFi as not needing polishing, as HiFi has accuracy \>99%, so shotgun Illumina is not going to do much. There is some argument for still polishing with your original long reads.

I can recommend any of the linked papers for ideas or further explanation.

---
Tags: assembly, software-recommendation, illumina, pacbio
---


title: Converter between PDB or mmCIF and MMTF

text:
# Question
Title: Converter between PDB or mmCIF and MMTF
I'd like to test MMTF, a new format for storing biomolecular structures which is promoted by RCSB as a more compact alternative to mmCIF and PDB.

From MMTF FAQ:

> * How do I convert a PDBx/mmCIF file to an MMTF file?
>
> The BioJava library contains methods to read and write PDBx/mmCIF files and MMTF files.

Can I do such a conversion, ideally from command line, but without writing my own Java program?

# Answer
I wrote a very quick and dirty script to handle conversion between file types using BioJava.

https://github.com/eedlund/Utils/tree/master/BioUtils

Download the jar file here

To run: `java -jar BioUtils.jar $FILE $TYPE`

where \\$FILE is a PDB or mmCIF file you'd like to convert and \\$TYPE is the format of the output file \[PDB, CIF, MMTF\].
> 2 votes

# Answer
You looked at this https://github.com/rcsb/mmtf-python

"The python implementation of the MMTF API, decoder and encoder."
> 2 votes

# Answer
You can do this with BioStructures.jl in Julia. All the 6 transformations between PDB/mmCIF/MMTF are possible.

For example, PDB to MMTF:

```
using BioStructures
struc = read(in_filepath, PDB)
writemmtf(out_filepath, struc)

```

mmCIF to MMTF:

```
using BioStructures
struc = read(in_filepath, MMCIF)
writemmtf(out_filepath, struc)

```
> 2 votes

# Answer
You can use PyMol to change the PDB to mmtf or other format. First open the PDB file and then choose file-\>export molecule-\> choose which format you want to save as
> 2 votes

# Answer
This python script uses Biopython and does the job well even on large structures. For example:

`python cif2pdb.py 4ckh-assembly-1.cif 4ckh-assembly-1.pdb`

This generates the structure below.
> 1 votes

---
Tags: file-formats, pdb, 3d-structure
---


# Question
Title: Do I have to download all nt.xx files to perform a search in blastn 2.12.0
I downloaded and extract blast nt.24 database, but when I run :

blastn -task blastn-short -query consulta.fa -db **nt.24** -out misalida.txt

it says: mdb\_env\_open: No such file or directory

and if a try with nt instead of nt.24 it says:

BLAST Database error: Could not find volume or alias file (nt.00) referenced in alias file (/root/my-documents/blastdb/nt/nt So my question is :

Do I have to download all nt.xx files?

And : Will it show the organisms in the output result?

# Answer
In general, yes, you do need to download and extract all (nt.00 - nt.46 as of now) of the nt portions in order to do a comprehensive BLAST against all of nt. The database has been split up for convenience (to make smaller individual files), so all are necessary to be comprehensive.

Once you download all the pieces, you’ll want to use `nt` as the database parameter to blastn.

`blastn -task blastn-short -query consulta.fa -db nt -out misalida.txt`

As for the second question - this depends on how you run the BLAST. Each accession is associated with a taxid, which will map to its taxonomy. This can be mapped back to the organism, but you may need to do some interpretation in order to convert the BLAST hits to the proper organism - this depends highly on what kind of sequences you are aligning.
> 1 votes

# Answer
Thanks, I downloaded nt.00 (by http firefox because wget via ftp always returned a corrupted file when decompressing it) and it works fine locally. The command is:

**blastn -task blastn-short -query consulta.fa -db nt.00 -out misalida.txt**

I also added a first line \>Seq1 to the query file consulta.fa (if not didn't work):**\>Seq1 GTGCC...**

And the results in misalida.txt first print the taxid with the organism name and chromosome. At the bottom it printed the positions details. What I'm doing is a genome editor, so when you enter some nucleotids(like 14 to 20 or more) it will search with blastn possible next nucleotids paths from many organisms with blastdbcmd.
> 0 votes

---
Tags: blastn
---


# Question
Title: How to Find List of Genomic Coordinates from List of Genes using UCSC Website?
Suppose I have a list of genes and I want to find their genomic coordinates. To be specific and to keep things consistent with the rest of my project, I'm using hg38 on the website genome.ucsc.edu. So far, I've been typing one gene symbol at a time into the gene browser and clicking the correct match. This sends me to the gene's page, which shows me the start and end genome coordinates for the gene's transcript and coding region. That's what I'm looking for.

This is a tedious process when I have hundreds of genes to go through, so I was wondering if I can input a list and get out a list instead. I know about the page https://genome.ucsc.edu/cgi-bin/hgTables. As an example, suppose I click "paste list" and just type the following two lines:

```
CASZ1
MTHFR

```

Here's a screenshot of the output (click to see full sized image):

I think this page is generally intended to go the opposite direction (genomic coordinates -\> genes)?

Edit: Ok, here's an example of what I want. I'd like to just enter the gene symbols and get the transcript and coding regions shown.

```
CASZ1, Transcript chr1:10,636,604-10,796,646, Coding Region chr1:10,638,942-10,693,889
MTHFR, Transcript chr1:11,785,723-11,805,964, Coding Region chr1:11,790,680-11,803,116
ST7L, Transcript chr1:112,525,883-112,618,928, Coding Region chr1:112,526,013-112,598,043

```

I obtained these by searching through genome.ucsc.edu one by one manually.

# Answer
> 1 votes
The current release of Gencode for `hg38` is v38. This is one way you could get a list of Gencode gene annotations:

```
$ wget -qO- ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_38/gencode.v38.annotation.gff3.gz \
    | gunzip --stdout - \
    | awk '$3 == "gene"' - \
    | convert2bed -i gff - \
    > genes.bed

```

To get `convert2bed`, you could install a toolkit called BEDOPS, which is available via Github or package managers. Refer to the installation instructions for more detail.

Once you have `genes.bed`, you don't need to make it again, unless Gencode publishes a new set of annotations.

You could create a text file called `genes_of_interest.txt` that contains your HGNC names of interest, e.g.:

```
$ cat > genes_of_interest.txt
CASZ1
MTHFR
ST7L

```

Press `Ctrl` \+ `D` at the end to close the file.

To search `genes.bed` for entries where there are matches to your genes of interest, you can use `grep`:

```
$ grep -Fwf genes_of_interest.txt genes.bed
chr1    10636603    10796650    ENSG00000130940.14  .   -   HAVANA  gene    .   ID=ENSG00000130940.14;gene_id=ENSG00000130940.14;gene_type=protein_coding;gene_name=CASZ1;level=1;havana_gene=OTTHUMG00000002035.5
chr1    11785722    11806920    ENSG00000177000.11  .   -   HAVANA  gene    .   ID=ENSG00000177000.11;gene_id=ENSG00000177000.11;gene_type=protein_coding;gene_name=MTHFR;level=2;havana_gene=OTTHUMG00000002277.5
chr1    112523517   112620825   ENSG00000007341.18  .   -   HAVANA  gene    .   ID=ENSG00000007341.18;gene_id=ENSG00000007341.18;gene_type=protein_coding;gene_name=ST7L;level=2;havana_gene=OTTHUMG00000011753.7

```

The above is just sample output. Your result may differ.

---
Tags: genome
---


# Question
Title: Build RepeatMasker default libraries beforehand (bioconda installation)
I'm using RepeatMasker on human sequences. When running RepeatMasker this way:

```
RepeatMasker -e rmblast -species human -dir . myFile.fa

```

RepeatMasker builds some libraries:

```
Building general libraries in: /path/to/share/RepeatMasker/Libraries/CONS-Dfam_3.3/general
Building species libraries in: /path/to/share/RepeatMasker/Libraries/CONS-Dfam_3.3/human

```

(Paths edited) On subsequent runs, this does not happen. I want to check if these libraries exist beforehand (which I can easily do with `ls /path/to/share/RepeatMasker/Libraries/CONS-Dfam_3.3/`), and if they don't, build them.

I could do a run with a decoy fasta to trigger that step, but it's a bit hacky. Is there a better way to force Repeat Masker to build those libraries?

(I mention that I'm using bioconda in the title because I'm working on a remote server with restrictions, so that should be taken into account, no because there is already a solution for other installation methods).

# Answer
> 1 votes
Something like this:

```
db_dir=$CONDA_PREFIX/share/RepeatMasker/Libraries/CONS-Dfam_3.3/

if [ -d "${db_dir}" ] && [ -n "$(ls -A ${db_dir})" ]; then
    echo "CONS-Dfam folder is not empty"
    # ... do something
else
    echo > sequence.fasta
    RepeatMasker -e rmblast -species human -dir . sequence.fasta
fi

```

If the folder exists and is not empty, otherwise build it using an empty file.

---
Tags: repeatmasker
---


# Question
Title: How to convert amino acid sequence to SMILES format
I have an amino acids string "MKGK" and I would like to convert it into SMILES format.

I tried two approaches.

First with PepSMI which gives me:

```
N[C@@]([H])(CCSC)C(=O)N[C@@]([H])(CCCCN)C(=O)NCC(=O)N[C@@]([H])(CCCCN)C(=O)O

```

Second, using ChemAxon Marvin's **molconvert** command

```
/Applications/MarvinSuite/bin/molconvert --peptide MKGK smiles

```

I get:

```
[H]N[C@@H](CCSC)C(=O)N[C@@H](CCCCN)C(=O)NCC(=O)N[C@@H](CCCCN)C(O)=O

```

Both of them give different result. Are they chemically the same? If not, which result should I trust?

# Answer
SMILES are meant for machines, but reading up on the format is very handy and one totally can read and write them —I do it often and could not recommend it enough.

You can have SMILES written differently and they give the same result —hence why you should never compare SMILES, but Inchi or better still a parsed molecule —e.g. in RDKit.

Looking at them, I can see that

* they are the same except for one explicit proton on the N-terminus.
* Neither are canonical SMILES —which is fine.
* Both are written N to C like biologist write peptides —chemists write peptides C to N due to solid phase synthesis.
* they are L-amino acids: the pattern `N[C@@H](*)C(=O)O` is the give away. Technically, the implicit hydrogen in the Cα \[C@@H\] has `@@`, which means counter-clockwise —Voet & Voet taught everyone to read C(=O)RN counter-clockwise with the N-C (the first bond) sticking out.

The first has rather oddly declared chiral protons, the second has an unnecessary N-terminus proton, which is not adding a charge: `[H]N*` is the same as `[NH2]*` as it is assumed to be `[H][NH]*`. So the correct version **at pH 7** should be:

`[NH3+][C@@H](CCSC)C(=O)N[C@@H](CCCCN)C(=O)NCC(=O)N[C@@H](CCCCN)C([O-])=O`
> 3 votes

---
Tags: computational-biochemistry, conversion, cheminformatics
---


# Question
Title: In Seurat, how can I export cluster IDs to csv files
I am doing scRNAseq analysis with Seurat.

I clustered the cells using the FindClusters() function.

What I want to do is to export information about which cells belong to which clusters to a CSV file.

In a Seurat object, we can show the cluster IDs by using Idents(・), but I have no idea how to export this to CSV files.

I would be grateful if you could show this by using the PMBC data (https://satijalab.org/seurat/articles/pbmc3k\_tutorial.html) as an example.

# Answer
> 2 votes
The clusters information, together with additional information on each cells (%mito, sample name, group, orig.ident,...) are stored inside `seuratObject@meta.data`. `write.csv(seuratObject@meta.data,"./seurat_metadata.csv")` will dump everything to a file.

The method @osmoc suggested, is using built in seurat functions to interact with `seuratObject@meta.data`. `seuratObject$graph_name_resolution` points to `seuratObject@meta.data$graph_name_resolution`.

I personally prefer to use the metadata table, as I have more control on how I interact with the data.

# Answer
> 2 votes
I've recently done this with my own Seurat data; here's the process:

```
# load tidyverse packages
library(tidyverse)
# create cluster ID table
tibble(cellID = colnames(pbmc), clusterID = Idents(pbmc)) %>%
  # write out to a file, with today's date
  write_csv(file = sprintf("cluster_mappings_%s.csv", Sys.Date()))

```

# Answer
> 0 votes
The clusters are saved in the field `object$seurat_clusters` where `object` represents your Seurat object.

You can also get them running the command `Idents(object)`.

**NOTE**: if you run the command `FindClusters` on different SNN graphs, or different resolutions, you need to specify the clusters by `object$graph_name_resolution`. Every time you run `FindClusters`, the idents of the cells are updated to the latest results, same for `seurat_clusters`.

For instance, let's say that you have run the command `FindNeighbors` on the `RNA` assay (default value).
Then you run `FindClusters` with `graph = "RNA_snn"` and `resolution=1`.
In that case, the clusters are in `object$RNA_snn_1`.

Once you have the clusters, you can save them as a CSV:

```
clusters <-  <here you use one of the described methods>
# data.frame
clust.df <- data.frame(cellID = colnames(object), clusters = clusters)
write.csv(clust.df, <path.of.the.csv.file>)

```

---
Tags: seurat, single-cell, clustering
---


# Question
Title: Download proteomes from NCBI based only on binomial names
I'm trying to download around 2000 proteomes from NCBI, but all I have is the genus and species e.g: `lophodermium_seditiosum trichoderma_harzianum`

Does anyone know a good way of doing this in script form?
My code at the moment is:

```
from Bio import Entrez
from Bio import SeqIO
import os

file_name = "proteome1.faa"
net_handle = Entrez.efetch(db="protein", id="trichoderma_harzianum", rettype="gp", retmode="text")
out_handle = open(file_name, "w")
out_handle.write(net_handle.read())
out_handle.close()
net_handle.close()
print("Saved")

for seq_record in SeqIO.parse("proteome1.faa", "genbank"):
    print(seq_record.id)
    print(repr(seq_record.seq))
    print(len(seq_record))

```

But I think it wants a number as the id, I don't know if it's possible to convert species name to a numerical ID within the script, or if that's the best way of going about this.

# Answer
The general procedure is explained here:

First step is to search the Taxonomy database with the organism name. Accepted common names usually work at all taxonomic levels. Use the scientific name or formal name if no results are obtained with the common name.

In code:

```
from Bio import Entrez

Entrez.email = 'user@example.org' # Put your email here

def get_tax_id(species):
    species = species.replace('_', '+').strip()
    search = Entrez.esearch(term=species, db='taxonomy', retmode='xml')
    record = Entrez.read(search)
    return record['IdList'][0]

if __name__ == '__main__':
    organisms = ['lophodermium_seditiosum', 'trichoderma_harzianum']
    taxids = [get_tax_id(organism) for organism in organisms]
    print(taxids)

```

This gives us the `taxid`s for our species, respectively `128023` and `5544`. Now we can use these to get the proteomes. For *Lophodermium seditiosum* we find two results, but for *Trichoderma harzianum* there are 65900 results, so you might want to tweak the query to filter the results. For example, if you only want the sequences appearing in PDB you add `AND pdb[filter]` so that the query becomes `txid5544[Organism:exp] AND pdb[filter]`.

```
def get_proteomes(taxid):
    query = 'txid{}[Organism:exp]'.format(taxid)
    handle = Entrez.esearch(db='protein', term=query)
    result = Entrez.read(handle)
    if int(result['Count']) <= 2000:
        ids = ','.join(result['IdList'])
        proteomes = Entrez.efetch(db='protein', id=ids, rettype='fasta').read()
        return proteomes
    else:
        raise RuntimeError('Too many results!')

if __name__ == '__main__':
    print(get_proteomes('128023'))

```

Returns:

```
>ADI44294.1 actin, partial [Lophodermium seditiosum]
CKAGFAGDDAPRAVFPSIVGRPRHHGIMIGMGQKDS

>ADI44293.1 actin, partial [Lophodermium seditiosum]
CKAGFAGDDAPRAVFPSIVGRPRHHGIMIGMGQKDS

```
> 3 votes

# Answer
You can use the NCBI Datasets to download the protein sequences of all annotated protein-coding genes for an organism starting with the organism name. There is a web-interface that allows exploring by taxonomy or specific organism search as well as a command line utility and an API. For example, the command line utility can be used as follows:

```
$ datasets download genome taxon 'trichoderma harzianum' \
    --annotated --exclude-genomic-cds --exclude-gff3 \
    --exclude-rna --exclude-seq --reference
Collecting 1 genome accessions [================================================] 100% 1/1
Downloading: ncbi_dataset.zip    3.92MB done

```

This downloads a zip archive containing the `protein.faa` file for the *Trichoderma harzianum* reference assembly.

Note, unlike the Entrez method described above, this will not download all proteins for a given organism from the NCBI Protein portal. Thus, *Lophodermium seditiosum* does not return any results as there is no available genome for this assembly in the NCBI Assembly database.

On the other hand, the NCBI Datasets method will be significantly faster and return cleaner results. I would argue that if you are interested in obtaining the 'proteomes' of a given set of organisms, this is a better choice as it includes only those cases where a complete proteome is available based on the annotation of the genome and does not include cases where only a few proteins from an organism are available.
> 1 votes

---
Tags: proteins, python, phylogenetics
---


# Question
Title: How can I identify acceptor and donor atoms using BioPython?
I know how to traverse through structures, models, residues, and atoms in a protein chain using BioPython (which is very easy).

How can I identify donor and acceptor atoms in a chain?

# Answer
> 3 votes
A hydrogen bond donor heavy atom is one that is protonated, while the acceptor has a lone pair. If you talking about amino acids, then it's simple as they are fixed.

| atom name | description | role | note |
| --- | --- | --- | --- |
| `N` | backbone nitrogen | donor | SS |
| `CA` | backbone Cα | apolar |  |
| `C` | backbone carboxyl C | sterically hindered and electrophilic | nucleophilic attack |
| `O` | backbone carboxyl O | acceptor | SS |
| `CB` etc. | sidechain carbons | apolar |  |
| `OXT` | C-terminal oxygen | acceptor | may be protonated especially if a membrane embedded C-terminus |
| Serine's `OG` / Threonine's `OG1` | sidechain hydroxyl | donor | SS: may form ST turns. Fairly uncommon: may be a donor |
| Asp `OD1`/`OD2`/Asn `OD1` | sidechain carboxyl/amide O | acceptor | SS: may form Asx turns |
| Glu `OG1`/`OG2`/Gln `OG1` | sidechain carboxyl/amide O | acceptor |  |
| Lysine `NZ` | sidechain amine | donor |  |
| Asn `ND2` / Gln `NG2` | terminal admide N | donor | the direction of an Asn/Gln may be accidentally filled if manually assigned as N/O have the same density in X-ray crystals |

Water can be both a donor and an acceptor. Two charged residues form a salt bridge more so than a hydrogen bond, because the interaction is driven by the Coulomb electrostatic forces of the two heavy atoms.

---
Tags: python, proteins, small-molecules
---


# Question
Title: In Pygraphviz, How can I assign edge colors according to edge weights?
I am now drawing a weighted directed graph using Pygraphviz.

The adjacency matrix A of the weighted directed graph is

```
A = pd.DataFrame({'X': [0.0,0.1,-0.8], 'Y': [-0.7,0.0,-0.1], 'Z':
    [0.4,0,-0.1]} ,index=["X", "Y", "Z"])

```

where A\_{i,j} indicates the edge weight from the node i to the node j.

I was able to draw the graph without edge color.

```
import networkx as nx
from networkx.drawing.nx_agraph import to_agraph
from IPython.display import SVG, display

G_nx = nx.from_pandas_adjacency(A,create_using=nx.DiGraph())
G_nx.graph['edge'] = {'arrowsize': '1.0', 'splines': 'curved'}
G_nx.graph['graph'] = {'scale': '3'}

Agraph_eg = to_agraph(G_nx)
Agraph_eg .node_attr["height"] = 0.3
Agraph_eg .node_attr["width"] = 0.3
Agraph_eg .node_attr["shape"] = "circle"
Agraph_eg .node_attr["fixedsize"] = "true"
Agraph_eg .node_attr["fontsize"] = 8
Agraph_eg.layout(prog="neato")
Agraph_eg.draw('graph_eg1.png')

```

But I have no idea how to assign edge colors and widths according to the following two rules.

#1 If weights are more than 0, edge color should be red. If weight are less than 0, edge color should be blue.

#2 Edge widths should be proportional to the absolute values of weights.

Do you have any idea how to do this??

# Answer
> 2 votes
This might help you

```
for edge in G_nx.edges(data=True):
    color = "black"
    weight = edge[2]["weight"]
    if weight > 0:
        color = "red"
    elif weight < 0:
        color = "blue"
    edge[2]["color"] = color
    edge[2]["penwidth"] = abs(weight)

```

This gives you the following graph:

Depending on your network size this solution might take some time. "penwidth" is the graphviz information for the edge width.

---
Tags: python, networks, graphs
---


# Question
Title: What statistics can I use on qPCR data composed of several runs if all normality tests fail on the dataset?
I have \>100 qPCR experiments that I'd like to analyze together, each containing the same set of genes (10 genes of interest and 2 reference genes). I have four different samples (untreated & 3 treatments) from the same cell culture. Our setup allows two of them to be run at once, i.e. the samples from the same cell culture are split to two qPCR runs.

1. I would like to do statistics with the dCq values after normalizing to the reference genes, but the dataset doesn't pass any normality tests. What methods are suitable for finding valid differences between treatments (or any other features like age)?
2. If the same amount of cDNA was used and reference genes were added to every run, is it still necessary to normalize between runs or between cell cultures?

# Answer
1. A wilcoxon test is what you're looking for.
2. Yes, since there will inevitably be some small batch effects to compensate for.
> 1 votes

---
Tags: statistics, normalization, qpcr
---


title: Using python to call Emboss wordcount

text:
# Question
Title: Using python to call Emboss wordcount
Hi guys I would like to run wordcount from Emboss in a bunch of fasta files. I tried this:

```
#!/usr/bin/env python
import sys
import subprocess

filenames = ['Downloads/Data/GCF_000016725.1_chr.fna']
for filename in filenames:
    name = filename.split('/')[2]
    print(filename)
    subprocess.call (['wordcount', 'Input sequence(s): ' + filename, 'Word size [4]: ' + sys.argv[1]])

```

But I got error messages:

```
Count and extract unique words in molecular sequence(s)
Error: Failed to open filename 'Input sequence(s): Downloads/Data/GCF_000016725.1_chr.fna'
Error: Unable to read sequence 'Input sequence(s): Downloads/Data/GCF_000016725.1_chr.fna'
Died: wordcount terminated: Bad value for '-sequence' and no prompt

```

But if a use shell with Downloads/Data/GCF\_000016725.1\_chr.fna it works fine!

Any tip would be helpful! I am using python, but I would have a method to compare my results!

Thank you by your time!

Paulo

# Answer
> 3 votes
You have to use command line options, not interactive prompts, like so:

```
subprocess.call(['wordcount', '-sequence', filename, '-wordsize', sys.argv[1]], '-outfile', filename + '.wordcount.txt')

```

Or, using the f-strings or formatted string literals:

```
subprocess.call(f'wordcount -sequence {filename} -wordsize {sys.argv[1]} -outfile {filename}.wordcount.txt')

```

Note that I had to add `-outfile` option to your system call, otherwise `wordcount` prompts for the output file. See `wordcount` manual.

---
Tags: python, emboss
---


# Question
Title: How to show electrostatic interactions in Pymol on the wanted residue
I have the following Pymol visual, The green molecule is receptor and red/yellow is the ligand. In the ligand I highlighted Lysine (K) residue as yellow.

What I want to do is to highlight and draw the electro static contact of the ligand especially with Lysine residue. Like this in Pymol:

How can I achieve that?

# Answer
Depends whether you work with named selections which is preferable and helpful if you want to use this Pymol session for future images. In this case, use the distance command for your two selections, one being the N of the K in your ligand and the other the O of the E of your protein.

Full syntax and examples here: https://pymolwiki.org/index.php/Distance

If you need a quick solution: in the menu Wizard, click on Measurement. Then click on the two atoms where you want the dotted line to appear in between.
> 4 votes

---
Tags: protein-structure, docking, pymol, structural-biology
---


# Question
Title: Constructing a phylogenetic tree for classical MHC genes from different organisms
I am interested in constructing a phylogenetic tree containing the classical HLA α1 and α2 domains and potential classical MHC α1 and α2 domains of a different organism. I am using the amino acid sequences for the tree.

However, I am unsure of what bioinformatic tools and software to use and what parameters to fill in.

Lastly any advice on how to read a phylogenetic tree and how to relate the sequence similarity between two sequences by looking at the tree would be very helpful.

My aim is to investigate the classical MHC genes of a non human organism.

Any help and advice would be greatly appreciated!

# Answer
This is a very unusual situation as a main effort is finding the outgroup where these two domains aren't duplicated.

## General situation

In animals a gene tree nearly always follows a species tree, with the exception of gene duplication events, which frequently stem from genome duplication events. Making a tree with lots of genes gets hard to follow and actually becomes pointlessly problematic when you get a branch that does not quite follow the species phylogram —say, you have rabbit, mouse and human genes and the inferred tree, due to poor signal, comes out with mouse+human together as opposed to rabbit+mouse. As a result a nice tree needs to present the species before the genome duplication, often sea lamprey —I wrote a blog post recently were I discuss nice species to pick for a multiple sequence alignment.

In this simple case, one can get the protein sequence of interest and blastp it against human, sea lamprey and fruit fly, see how divergent are the paralogues in human and if the split is more recent that the latter two, if so, the latter will be an outgroup and after getting a few representative species in between, one can make a MSA with Muscle, which contains also a tool to make a neighbour joining tree.

## Domain duplication

In your case you have something more complicated in that you have a series of paralogues, each gene encoding a protein (e.g. HLA-A) containing two globular domains, which form the MHC-1 fold, and you want make a tree of this... which means you want find out when the two domains duplicated (i.e. outgroup)? If the domains work together (correct?) and have the same fold (as you alluded to, correct?), this would mean that before the duplication event of the domains, the ancestral protein may have worked as a homodimer.

I am not familiar with the literature, so do not know if this has already been looked at —but I suspect it has.

If not, clicking on the sequence of one of the two domains in the Uniprot entry of one of the paralogues and doing a few iterations of psi-blast against human and a few other distant species (such as zebra fish, sea urchins, which are basal to vertebrates, and fruit fly, which I do not believe have a MHC etc.) most of these will have two hits per protein... but one might not. This is the outgroup!
> 1 votes

---
Tags: phylogenetics, phylogeny
---


# Question
Title: I have two fasta files READ1 and READ 2. How to identify which one is the reverse and the forward?
As you can see in the picture I have two files from the paired end reads. I would assume READ1 would be forward reads and READ2 would be reverse reads, but my lab colleague told me that that might not be the case. How do I confirm?

# Answer
> 3 votes
Unless this is RNA-seq prepared with a directional library, or something else that you're sure is directional, it is irrelevant.

DNA is complementary, so a DNA molecule can be turned into a sequence-able fragment with either strand given any pair of Illumina adapters and indices.

As far as Illumina lingo goes, R1 is the first read that is read off the sequencer, or what people call the "forward read". R2 is the second pass of the same clusters, and some people call it the "reverse read". I recommend just calling them R1 and R2 reads to avoid confusion. Sequencing by synthesis produces the cleanest reads in the first pass, so R1 reads usually have better Q values. You can check this easily by comparing the quality scores with fastqc.

You can learn more about Illumina sequencing chemistry here, or check out some other answers on the bioinformatics SE or biostars.

# Answer
> 0 votes
I think @conchoecia answer is correct but maybe not clear enough for someone making himself that question for the first time (like was my case recently).

I think the simple way to explain this is that in F\_NL\_READ1.fq you have random reads, in mixed directions, and in F\_NL\_READ2.fq you have the paired read in the opposed direction for every read in F\_NL\_READ1.fq, but not necessarily and almost certainly, not every read in every file will go in the same sense (forward or reverse).

At least for paired-end reads from Illumina.

# Answer
> -1 votes
If this is RNASeq, then the question might have an answer (or it might be nonsense), but you have to know what kind of library prep was used.

Or just align and inspect the results to figure it out.

---
Tags: fasta
---


title: BUSCO results apparently inconsistent

text:
# Question
Title: BUSCO results apparently inconsistent
I have two assemblies obtained from two different softwares for the same run data. In one of them, I get a BUSCO score of 69,4%. In the other one, I get 68,5%. The first assembly covers a 99,7% fraction of the reference while the second one covers 100%.

The relevant thing is that the reference genome also gets 68,5% (less than the first assembly). How can I interpret this situation?

I identify these options:

* The first assembly is showing a different strain.
* BUSCO is imprecise
* The first assembly is merging DNA from another organism
* The first assembly is a misassembly, repeating part of the sequence and BUSCO interprets that as new genes
* The sample really is an organism that has more genes than the reference.

The organism is a bacteria and I am using bacteria\_odb10 as the lineage dataset.

This is the BUSCO result for the reference (LN554883.1 Chlamydophila abortus strain 1H):

```
C:68.5%[S:68.5%,D:0.0%],F:8.1%,M:23.4%,n:124
85  Complete BUSCOs (C)
85  Complete and single-copy BUSCOs (S)
0   Complete and duplicated BUSCOs (D)
10  Fragmented BUSCOs (F)
29  Missing BUSCOs (M)
124 Total BUSCO groups searched

```

This for the first assembly:

```
C:69.4%[S:69.4%,D:0.0%],F:8.9%,M:21.7%,n:124
86  Complete BUSCOs (C)
86  Complete and single-copy BUSCOs (S)
0   Complete and duplicated BUSCOs (D)
11  Fragmented BUSCOs (F)
27  Missing BUSCOs (M)
124 Total BUSCO groups searched

```

And this for the second assembly:

```
C:68.5%[S:0.0%,D:68.5%],F:8.1%,M:23.4%,n:124
85  Complete BUSCOs (C)
0   Complete and single-copy BUSCOs (S)
85  Complete and duplicated BUSCOs (D)
10  Fragmented BUSCOs (F)
29  Missing BUSCOs (M)
124 Total BUSCO groups searched

```

I see now that the first assembly also contains a lot of contigs with human DNA. I guess BUSCO finds a match for some genes in the human DNA?

# Answer
The core idea behind BUSCO is that there are a number of genes that are Unique Single-Copy Orthologues (USCOs) that should be present in most genomes and should only appear *once* in the genome. This is only a guess, however, and not a universal rule.

Sometimes genomes can have multiple copies of these USCOs, and these duplicated genes reduce the overall score (which is based on the proportion of complete *single-copy* genes that are present as a fraction of the entire USCO set).

Regarding human DNA, bacteria are a curiously diverse group, with genomes that have more similarities in function than in form. Given this diversity, it wouldn't surprise me if contaminating DNA could increase the BUSCO score.
> 1 votes

---
Tags: assembly, busco
---


# Question
Title: Best R program to use for RNA seq analysis
We have an RNA seq taken from a the same tissue in two different developmental stages. We have one gene, with the number of reads from stage A and stage B, and we want to see which genes are differentially expressed.

After performing RNA-seq, I am left with a CSV file containing a list of genes and the amount of reads in two different cell types. We have one list of genes, with the number of reads in Stage A in one column and the number of reads in Stage B in another column. We are looking to see what types of genes are differentially expressed in the two stages.

I have used DESeq2, but I'm not sure if this is the best/correct to analyze what I am looking for.

What is the best package to analyze the data for significant differential gene expression?

Thanks

# Answer
> 1 votes
DESeq2 works for that. Unless you literally have only two samples.

There is **no** statistically sound way of analyzing DE genes with literally only two samples. This is a bad experimental design, suitable for only the roughest of exploratory purposes. Only the most obvious differences (that you probably knew already) are likely to be borne out in a legitimate experiment.

---
Tags: rna-seq, deseq2
---


# Question
Title: How to extract the protein sequences of a genbank file using R or biopython
sorry for the question, I'm trying to extract the proteins sequences from a genbank file.

```
 gene            complement(516466..532086)
                 /gene="rtxA"
                 /locus_tag="VV1_RS17390"
                 /old_locus_tag="VV2_0479"
 CDS             complement(516466..532086)
                 /gene="rtxA"
                 /locus_tag="VV1_RS17390"
                 /old_locus_tag="VV2_0479"
                 /inference="COORDINATES: similar to AA
                 sequence:RefSeq:WP_011081430.1"
                 /note="Derived by automated computational analysis using
                 gene prediction method: Protein Homology."
                 /codon_start=1
                 /transl_table=11
                 /product="MARTX multifunctional-autoprocessing
                 repeats-in-toxin holotoxin RtxA"
                 /protein_id="WP_011081430.1"
                 /translation="MGKPFWRSVEYFFTGNYSADDGNNSIVAIGFGGEIHAYGGDDHV
                 TVGSIGATVYTGSGNDTVVGGSAYLRVEDTTGHLSVKGAAGYADINKSGDGNVSFAGA
                 AGGVSIDHLGNHGDVSYGGAAAYNGITRKGLSGNVTFKGAGGYNALWHETNQGNLSFA
                 GAGAGNKLDRTWFNRYQGSRGDVTFDGAGAANSISSRVETGNITFRGAGADNHLVRKG
                 KEGNHTANLANEDISSANGYHSMGKGGYSLSDLHYSVNAVRSTSETVADIDEYTDQTL
                 FKPATDSGESSGDVRFNGAGGGNVIKSNVTRGNVYFNGGGIANVILHSSQFGHTEFNG
                 GGAANVIVKSGEEGDLTFRGAGLANVLVHQSKQGKMDVYAGGAVNVLVRIGDGQYLAH"

```

I'm using a code from a previous answer to extract the protein and the locus\_tag, but I want some extra features like:

```
>locus_tag|old_locus_tag|gene|product|complement|protein_id|chromosome|Accession|organism|strain
MGKPFWRSVEYFFTGNYSADDG ..........

```

the result will be like:

```
>VV1_RS17390|VV2_0479|rtxA|MARTX multifunctional-autoprocessing repeats-in-toxin holotoxin RtxA"|516466..532086|WP_011081430.1|Chromosome-II|NC_004459.3|Vibrio vulnificus|CMCP6
MGKPFWRSVEYFFTGNYSADDG ..........

```

this is the code that I'm using, how to add those features, or an other way using R ???

```
from Bio import SeqIO
file_name = 'CMCP6.gb'
# stores all the CDS entries
all_entries = []
with open(file_name, 'r') as GBFile:
    GBcds = SeqIO.InsdcIO.GenBankCdsFeatureIterator(GBFile)

    for cds in GBcds:
        if cds.seq is not None:
            cds.id = cds.name
            cds.description = ''
            all_entries.append(cds)

# write file
SeqIO.write(all_entries, '{}.fasta'.format(file_name[:-3]), 'fasta')

```

But the result from the code only add the locus\_tag and the sequence:

```
>VV1_RS17380
MFIKNLSIGKKIAAAFSIIAVINIAFGIFLSTELNTVKSELLNYTEDTLPAMEKVDAVRD
KISYWRRTQFAVFAMNDENQIKQTITRNEGIRREIETELAAYGKSVWPGEEEQTYNRLMS
LWSGYLSTMDKFNDALLAGDKDAAYPILTNSLSTFESIETEVNKLVMILKGAMDSNKNQI
LSSVNGLNTTAVISNIAIFAVMVIMTLLLTRIICGPLNIVVRQANAIARGDLSQDLDRNA
IGNDELGELADATIKMQDDLRQVIDNVIAAVTQLSSAVEEMTQISEMSASGMKDQQMQVT
LVATAMTEMKAAVADVARNTEDSASQAYDANRRTQDGAKETHQMVVSIEEVADIIAKAGD
TVAELEAQSNQINVVVDVIRGIADQTNLLALNAAIEAARAGESGRGFAVVADEVRTLAGR
TQDSTGEITSIIEKLQDLAKQAKHATEDSRTSIGACVEQGNNAQELMGSIEKSIANIADM
GAQIATACGQQDSVAEELSRNIENIHMASQEVAQGSQQTAQACRELTQLSVSLQDVMSRF
KLN
>VV1_RS17385
MNFKKTLLSIAIASASLTPAFSYSAPLLLDNTVHQTSQIAGANAWLEISLGQFKSNIEQF
KSHIAPQTKICAVMKADAYGNGIRGLMPTILEQQIPCVAIASNAEAKLVRESGFEGELIR
VRSASTSEIEQALSLDIEELIGSEQQARELASLAEKYSKTIKVHLALNDGGMGRNGIDMS
TERGPKEAVAIATHPSVAVVGIMTHFPNYNAEDVRTKLKSFNQHAQWLMESAGLKREEIT
LHVANSYTALNVPEAQLDMVRPGGVLYGDLPTNPEYPSIVAFKTRVASLHSLPAGSTVGY
DSTFTTANDAVMANLTVGYSDGYPRKMGNKAQVLINGQRANVVGVASMNTTMVDVSNIKG
VLPGDEVTLFGAQKNQHISVGEMEENAEVIFPELYTIWGTSNPRFYVK
>VV1_RS17390
MGKPFWRSVEYFFTGNYSADDGNNSIVAIGFGGEIHAYGGDDHVTVGSIGATVYTGSGND
TVVGGSAYLRVEDTTGHLSVKGAAGYADINKSGDGNVSFAGAAGGVSIDHLGNHGDVSYG
GAAAYNGITRKGLSGNVTFKGAGGYNALWHETNQGNLSFAGAGAGNKLDRTWFNRYQGSR
GDVTFDGAGAANSISSRVETGNITFRGAGADNHLVRKGKVGDITLQGAGASNRIERTRQA
EDVYAQTRGNIRFEGVGGYNSLYSDVAHGDIHFSGGGAYNTITRKGSGSSFDAQGMEYAK

```

Thanks so much

# Answer
ok tried a little bit, not sure is the most elegant/fastest way of doing it

starting from file from : https://www.ncbi.nlm.nih.gov/nuccore/NC\_004460.2?report=gbwithparts&log$=seqview {GenBank(Full)}

here my code:

```
from Bio import SeqIO

from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord

file_name = 'CMCP6.gb'

# stores all the CDS entries
all_entries = []

gb_record = SeqIO.read(open(file_name,"r"), "genbank")

gb_feature = gb_record.features

accession = gb_record.annotations['accessions'][0]+'.'+str(gb_record.annotations['sequence_version'])

# accession = gb_record.annotations['accessions']
# print(gb_record.annotations['accessions'])

for seq_feature in gb_feature:
    if seq_feature.type=="source":
            organism = str(seq_feature.qualifiers['organism'][0]).rsplit(" ", 1)[0]
            strain = seq_feature.qualifiers['strain'][0]
            chromosome = seq_feature.qualifiers['chromosome'][0]

for seq_feature in gb_feature:
    if seq_feature.type=="CDS":
            if 'translation' in seq_feature.qualifiers:
                if seq_feature.qualifiers['translation'][0] != '':
                    if 'gene' in seq_feature.qualifiers:
                        gene = seq_feature.qualifiers['gene'][0]

                    else:

                        gene = 'not defined'

                    if 'product' in seq_feature.qualifiers:
                        product = seq_feature.qualifiers['product'][0]

                    else:

                        product = 'not defined'

                    if 'protein_id' in seq_feature.qualifiers:
                          protein_id = seq_feature.qualifiers['protein_id'][0]

                    else:

                          protein_id = 'not defined'

                    if 'old_locus_tag' in seq_feature.qualifiers:
                          old_locus_tag = seq_feature.qualifiers['old_locus_tag'][0]

                    else:

                          old_locus_tag = 'not defined'

                    if 'locus_tag' in seq_feature.qualifiers:
                          locus_tag = seq_feature.qualifiers['locus_tag'][0]

                    else:

                          locus_tag = 'not defined'

                    if seq_feature.location.strand == 1:
                          complement = 'no complement'

                    if seq_feature.location.strand == -1:
                        complement = str(seq_feature.location).strip('[]').split('(-)')[0]

                    pippo = SeqRecord(
                                Seq(seq_feature.qualifiers['translation'][0]),

                                id = seq_feature.qualifiers['locus_tag'][0] ,

                                description = ('|'+str(old_locus_tag)+'|'+str(gene)+'|'+str(product)\
                                                +'|'+str(complement)+'|'+str(protein_id) +'|Chromosome-'+str(chromosome)\
                                                +'|'+str(accession).strip("[]").replace("'", "")+'|'+str(organism)\
                                                    +'|'+str(strain))
                                )
                    all_entries.append(pippo)

# print(all_entries)
# write file
SeqIO.write(all_entries, '{}.fasta'.format(file_name[:-3]), 'fasta')

```

output for your example line is :

```
>VV1_RS17390 |VV2_0479|rtxA|MARTX multifunctional-autoprocessing repeats-in-toxin holotoxin RtxA|516465:532086]|WP_011081430.1|Chromosome-II|NC_004460.2|Vibrio vulnificus|CMCP6
MGKPFWRSVEYFFTGNYSADDGNNSIVAIGFGGEIHAYGGDDHVTVGSIGATVYTGSGND.....

```

I got extra space between locus\_tag and | , I dont know how to get rid of it

Its not clear to me what you mean by Accessions, tried two ways (second one is past # in code) and got `NC_004460.2` instead of `NC_004459.3`
> 3 votes

---
Tags: r, python, biopython
---


# Question
Title: Detection of CNV(InDel) of intermediate size
My impression is that small InDel (a couple of bp) is identified through cigar string in BAM and typical CNV (at least thousands of bp) is detected through read depth.

What about InDel or CNV with size these them, say hundreds of bp? These CNV is too long to be covered by read and too short to be detected by statistically comparing the difference of read depth.

Do we have consensus that this kind of CNV is hard to detect? Or is there a common strategy of detection for them? Do we generally think them having impact on protein function (especially in oncology background) ?

# Answer
> 1 votes
> *My impression is that small InDel (a couple of bp) is identified through cigar string in BAM...*

Well, yes, by gaps in the alignment, easily detected by examining the CIGAR string.

> *...and typical CNV (at least thousands of bp) is detected through read depth.*

Also by deviations from expected read pair distance or orientation.

> *What about InDel or CNV with size these them, say hundreds of bp? These CNV is too long to be covered by read and too short to be detected by statistically comparing the difference of read depth. Do we have consensus that this kind of CNV is hard to detect?*

Yes, what you describe has certainly been recognized and acknowledged in the literature. It was one of the primary motivations for a recently published paper I contributed to on *de novo* (germline) variant discovery. Here's an excerpt from that paper:

> *In a reference-mapping context, calling indels with confidence requires accurate mapping of each read spanning the indel, with all gaps arranged consistently. This is possible only for short indels and tends to be prone to error and misalignment. Thus prediction of indels with length 10 bp has proved to be very challenging and accompanied by high false-positive and false-negative rates. Furthermore, the prediction of SVs via read mapping is only possible through indirect signatures such as alterations in read depth or read-pair signatures. These signatures can be quite noisy and result in high rate of false-negative and false-positive prediction.*

As for your question:

> *is there a common strategy of detection for them?*

Numerous mapping-free k-mer based methods for variant detection are coming up in the literature, and all seem to perform pretty well. The trick is how to define which k-mers are interesting—then you grab those k-mers (or the reads containing them) and assemble them into variant-spanning contigs. A few examples:

* In the Kevlar paper, we defined "interesting" k-mers as those with high abundance in an individual of interest but zero (or very low) abundance in the individual's parents.
* The HAWK paper, which describes a k-mer based GWAS-like method, defines "significant" k-mers as those significantly over- or under-represented in one class of samples versus a second class.
* The NovoBreak paper compares paired tumor and normal WGS data to find k-mers that span mutatino breakpoints in cancer genomes.

> *Do we generally think them having impact on protein function (especially in oncology background)?*

There is definitely *potential* that indels can have a substantial impact on protein function. But since their detection has been so difficult for so long, a lot remains unknown about their *actual* collective impact on function.

# Answer
> 1 votes
There is a trade off between read length and the limit of the indel that can be accurately mapped. Intermediate-size indels specifically are better detected with longer read lengths.

---
Tags: variant-calling, human-genome, cnv, indel
---


# Question
Title: Given SNP-Chip data for a population, what tool should I use to reconstruct haplotypes?
I think Eagle or haploSep will do the job of reconstructing haplotype data, but I'm not sure what the 'standard' or best practice tool to use is.

I have a VCF/PLINK format file of ~8,000 individuals with ~150,000 SNPs each (from across a population with several sub-populations).

I want to use the resulting data to place new samples into one or more populations (e.g. ancestry analysis).

Before I go trying to make one of these tools work / look at accuracy, etc., am I doing it right in the first place?

# Answer
HaploSep uses allele frequency data - I presume you have genotype data.

In that case, I would suggest looking at either Beagle5.2, shapeit4 and Eagle2. I believe they more or less give the same accuracy (Eagle2 is slightly less accurate than the other two) and Beagle 5.2 is faster than the rest. Depending on what you plan on doing with the haplotypes, I don't think that the choice of phasing algorithm will affect downstream results a great deal.
> 1 votes

---
Tags: snp, phylogenetics, haplotypes, phasing
---


# Question
Title: Using a phased reference assembly to phase sample genotypes
I've been trying without success to phase the variants of about 30 samples in a 5Mb non-recombining region. I have a phased genome assembly of the region but so far I have not been successful in using that to inform the phasing of the samples which have been sequenced to high coverage with short reads. Non-reference panel phasing using Beagle 5.2 and Shapeit4 have produced haplotypes that are composed of many small phased blocks which switch phase between them, far to many to manually correct.

I have created a "reference panel" of SNPs of known phase for the region and want to use it to guide the rest of the phasing of the non-"reference panel" SNPs but every method I have found only uses the intersection of the reference panel and the samples for analysis. Is there some technique that I can use to phase the genotypes in the samples? I believe this is a tractable problem since recombination is suppressed in this region so haplotype blocks should be large and contiguous.

# Answer
> 0 votes
The fact that the reference assembly is phased doesn't do anything to help actually phase the 30 other samples. There is no way of knowing what the combination of variants on single haplotypes are in the 30 samples from the one reference genome + 30-sample VCF alone. With this info, it is possible to come up with a measure of "HaplotypeA from sample 1 is more similar to HaplotypeII than HaplotypeI from my phased reference", though.

To actually phase the 30 other samples it will take long reads and/or proximity ligation data to phase the variants over a 5Mb region. Hi-C, Chicago, or DNAse Hi-C/Dovetail Omni-C reads will be necessary to correctly phase all of the variants over a region that is as large as 5Mbp.

I recommend using the tool HapCUT2 to use the reads to phase the variants already called for each sample. When mapping the reads to your genome, only map the reads to a fasta file that contains a single haploid copy of the genome (2 homologous chromosomes/haplotypes represented by 1 scaffold), not a fasta file that contains both haplotypes from the phased assembly (2 homologous chromosomes/haplotypes represented by 2 scaffolds).

---
Tags: haplotypes, phasing
---


# Question
Title: Counting hexamers in fasta sequence and identify its structure (and interruptions)
I have a lot of fasta files, each one with thousand of reads containing the hexameric motif "CCCTCT". The hexameric motif is highly continuous in most cases but interruptions may occur. I need to count the hexameric motif keeping the read ID and identify their structure (CCCTCT)n. These interruptions are real and do not come from Sequencing error

Example Sequence:(interruption in bold)

> ACATTTTTTTTCCACATCTGATGTGGAAAAAAAAAAAAATGAAATAGCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCT**CCT**CCCTCTCCCTCTCCCCACGGTCTCCCTCTCATGCGGAGCCAAAGCTGGACTGTACTGCT

Strategy 1: Count the CCCTCT motif

```
grep -on "CCCTCT" fasta.fa | cut -f 1 -d ":" | sort | uniq -c

```

> Output: `23 CCCTCT`

As you can see we have one interruption (`CCT`) in this sequence which the codes doesn't detect.

Strategy 2: Identify reads with the same structure and count repeats in each read

```
from Bio import SeqIO
from collections import defaultdict

dedup_records = defaultdict(list)
for record in SeqIO.parse("fasta.fa", "fasta"):
    # Use the sequence as the key and then have a list of id's as the value
    dedup_records[str(record.seq)].append(record.id)
with open("Output.fasta", 'w') as output:
    for seq, ids in sorted(dedup_records.items(), key=lambda t: len(t[1]), reverse=True):
        # Join the ids and write them out as the fasta
        output.write(">{}_counts{}\n".format(ids[0], len(ids)))
        output.write(seq + "\n")

```

Output: list ranking of reads with identical structures but does not show the structure

```
>ID/ccs_counts2
CCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCTCCCTCTCCCTCT.....

```

Can you please give me an idea of how to obtain the following output? example of a sequence of 23 hexamers with interruptions.

```
(CCCTCT)20 [INTERRUPTION] (CCCTCT)3

```

I have used expansion hunter, repeatFinder and repeatAnalysis-tool by PacBio but they were not particularly sensitive for this case. Any help would be appreciated.

# Answer
> 1 votes
I was going to write some clever code counting *k-mer*s, but this problem can easily be solved using the handy sequentially-operating feature of python's str.replace() method.

This code doesn't work if you need to look at the motif occurring on the other strand. Just the forward strand.

```
#!/usr/bin/env python3

mystring = "CCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCC\
TCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCT\
CCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCCTCTCCTCCCTCTCCCTCT"

motif = "CCCTCT"

# This replaces the motif with 1s
newstring = mystring.replace(motif, "1")
# results in 1111111111111111111111111111111111111CCT11
# silly, right?

# This whole block just parses that string and prints everything
#  the way that OP wanted.
done = False
space=""
accumulator = 0
while not done:
    if len(newstring) == 0:
        done = True
        break
    else:
        if newstring[0] == "1":
            accumulator += 1
        else:
            if accumulator != 0:
                print("{}({}){} ".format(space, motif, accumulator), end = "")
                accumulator = 0
            print(newstring[0], end = "")
            space = " "
    newstring = newstring[1::]

# This is our check after exiting the while loop.
#  This only fires if we stepped off cleanly at the
#  end of a motif.
if accumulator > 0:
    print("{}({}){} ".format(space, motif, accumulator), end = "")
print()

```

Prints, `(CCCTCT)37 CCT (CCCTCT)2`, matching OP's spec.

---
Tags: fasta, python, phylogenetics, structural-variation
---


# Question
Title: Allele Count and Allele Frequency in VCF files
I'm working in bioinformatics, but my computational skills far outstrip my knowledge of biology or genomics. So forgive the noobish question.

According to the VCF specification, the `INFO` column for each coordinate `[CHROM, POS]` can contain values for

**AC :** allele count in genotypes, for each ALT allele, in the same order as listed
**AF :** allele frequency for each ALT allele in the same order as listed: use this when estimated from primary data, not called genotypes
**AN :** total number of alleles in called genotypes

For example, I'm working with the the VCF file you can download from Kaviar.

I think I understand how to calculate those numbers from raw sequences from a pileup. If for a given coordinate `[CHROM, POS]`, $N\_0$ is the number of hits for the reference allele, $\textrm{AC}\_i \equiv N\_i$ is the number of hits for the $i^{th}$ variant (maybe only one), and AN is the total number of alleles counted: $\textrm{AN} \equiv \sum\_0^n N\_i$, then the allele frequency for the $i^{th}$ variant is

$$\textrm{AF}\_i = \frac{\textrm{AC}\_i}{\textrm{AN}} = \frac{N\_i}{\sum\_0^n N\_i}$$

So the phrase "total number of alleles in the dataset" would mean "total number of hits to coordinate `[CHROM, POS]` for *any* allele. And by "hit" I mean an aligned sequence fragment that includes that coordinate.

But I may not be taking into account the caveat for AF: "use this when estimated from primary data, not called genotypes."

Am I correctly understanding AF, AC, and AN in VCF files?

# Answer
> 2 votes
No, that is not quite right. specifically, this part:

> And by "hit" I mean an aligned sequence fragment that includes that coordinate.

To get data on how many reads contributed to calling a specific allele at a given position, one would use the DP field, for example, which isn't included in the Kaviar VCF.

The AC, AN, and AF fields in a VCF file are meant to be used to show the allele's frequency in the context of all of the samples used when making that VCF file.

**Investigation process**

This is one line from the Kaviar vcf file:

```
1       10230   rs200279319     AC      AA,A    .       .       AF=0.0002654,0.0048525;AC=7,128;AN=26378;END=10231

```

If we look at the README attached with the Kaviar vcf file, it says,

> "Kaviar is a compilation of SNVs, indels, and comlex variants observed in humans, designed to facilitate testing for the novelty and frequency of observed variants."

How many humans though? If we go to the Kaviar website it says:

> Kaviar contains 162 million SNV sites (including 25M not in dbSNP) and incorporates data from 35 projects encompassing 77,781 individuals (13.2K whole genome, 64.6K exome).

If we look at the Kaviar VCF header we see that the fields are defined:

```
##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency">
##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele Count">
##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in data sources">

```

If we look at these fields again, `AF=0.0002654,0.0048525;AC=7,128;AN=26378`, the AN=26378 tells us that for this locus, the Kaviar dataset contains information from 26378/2=13189 humans. That's the 13.2K whole genomes the Kaviar website was talking about. This site I'm looking at must be a non-coding region, because there is no exome data (protein-coding regions).

For this same site, AC=7,128 means that out of 13189 humans (26378 haplotypes), we saw the variants 7 and 128 times. Out of 26378 haplotypes, this is a frequency of AF=0.0002654,0.0048525.

---
Tags: vcf, allele-frequency
---


# Question
Title: How to compare sequences of genes obtained after whole-genome sequencing to reference genome?
My idea is to align whole-genome sequencing data (as fastq files after, 30× coverage, gDNA) to the reference mouse genome (NCBI), extract the immunoglobulin locus and compare it to the reference. I think the question is too broad but maybe a hint how to start ? (I'm new to the field.)

Task is: to check for difference between reference genome and a sample for 100-200 genes.

I checked these: question (but it's closed), answer (gave a hint about whole reference mapping not only the locus of interest), question (not answered), question (useful I guess but not sure).

Additional info: Instrument Platform: ILLUMINA
Instrument Model: HiSeq X Ten
Library Layout: PAIRED
Library Strategy: WGS
Library Source: GENOMIC

Typical sequence length distribution from FastQC: 150nt.

I have Win10 with WSL 1 installed.

# Answer
> 2 votes
Surprisingly there aren't actually many helpful one-liners to map Illumina WGS reads on the internet - everyone just sort of assumes that people can figure it out from the manual. Many people end up doing it in a not-so-ideal way, in that they stop at the `.sam` files that `bwa mem` makes that take up a lot of disk space.

I like this incantation:

```
# index the reference
bwa index reference_genome.fasta

# map your reads to the reference (bwa mem)
# compress that output to a bam file (samtools view -hb)
# sort the bam file (samtools sort)

bwa mem -t <number_of_threads> reference_genome.fasta your_sample.R1.fastq.gz your_sample.R2.fastq.gz | \
samtools view -hb - | \
samtools sort - > your_reads_to_reference.sorted.bam

# index the bamfile. makes a .bai file

samtools index your_reads_to_reference.sorted.bam

```

With this sorted `.bam`, the `.bai`, and the reference genome you can now use IGV to directly look at the reads mapped to your genome.

Also with the bam, you can generate a VCF file using `freebayes` or another tool like `samtools mpileup`. This will give you a file of how your sample is different from the reference genome. You can use `bedtools intersect` to get the variants only from your gene/regions of choice.

# Answer
> 1 votes
This is very easy. First map all reads to the reference sequence for the organism. If they are short reads use BWA mem. If long reads use minimap2 or ngmlr. Once mapped view the locus of interest in igv (GUI). If you need to extract the reads at the locus pass the coordinates to samtools view (cli).

---
Tags: sequence-alignment, fastq, reference-genome
---


# Question
Title: Understanding exercise on file coverage with question on summary statistics
I'm doing an exercise that asks for two files:

* Input 1: `A target file` (.bed format) contains multiple regions from `chr7:40000000-50000000` of human reference genome GRCh37 (hg19)
* Input 2: `Refseq exon list file` (.bed format) for all human coding genes (hg19 position)

The final goal is:

For all genes located in `chr7:40000000-50000000`, get the summary statistics of the target file coverage. (For each gene, get the fraction of exonic bases that was covered by the target file).

* I believe what they refer in this exercise is that the target file should be something like the whole chr sizes of hg19 `hg19.chr.sizes` and `refseq_exon_list` the list of exons from reqseq database. Both can be downloaded from tools like table browser .Is that correct ?
* I'm not sure which files I should download here to perform this task. Once downloaded the file I believe what I need to do is to restrict the `refseq_exon_list` by the requested region and then perform coverage with something like `bedtools`. Something on those lines `bedtools coverage -a hg19.chr.sizes -b reqseq_exon_list`.

Am i right here ? Any input is appreciated, thanks.

# Answer
> 1 votes
Based on the description in your question, I don't think this assumption is correct:

> I believe what they refer in this exercise is that the target file should be something like the whole chr sizes of hg19 hg19.chr.sizes \[. . .\]

I understand that the question is asking you to simply produce a random selection of sub-regions from chromosome 7, all of which fall within the larger region 40000000-50000000. It can't be the list of chromosome sizes since, first of all, the question is asking you to get regions from chr17 only, and second because that just wouldn't make sense: only chr17 is relevant here, so the rest would be useless and, in any case, the length of the chromosome itself isn't helpful information.

The idea here is to demonstrate that given a set of genomic regions and a list of start and end coordinates of exons, you can calculate how many of your exons overlap the target regions. Therefore, the steps to follow would be:

1. Generate random intervals in the 40000000-50000000 range.
2. Get the bed file for all human coding genes on hg19. You should indeed be able to download this from the UCSC genome browser, among other places.
3. Use a tool, `bedtools` or build your own using `awk`, or anything else that takes your fancy and for each gene found in the chr7:40000000-50000000 range, calculate how many bases of its exons fall within your randomly chosen selection of regions.

Of course, the best thing to do would be to ask your teacher to clarify.

# Answer
> 0 votes
This sounds like homework - can you post what you have tried so far + a few lines of the input files + your code? Also, please add the "homework" tag if this is homework, or a task from a learning exercise.

Without giving it away think about the information in your `Refseq exon list file`. Is it a 3-field bed file with just coordinates, or is it the expanded file format with fields 4-12? Do you really need to download the RefSeq annotation, or can you parse the information in your bed file to answer your question?

It sounds like this is your ultimate goal:

> For each gene, get the fraction of exonic bases that was covered by the target file.

What information from the output of `bedtools coverage` will give you this answer? Does this information exactly match the spec you provided? Are there other tools in `bedtools` you could use to achieve the same thing?

Side note: It isn't 100% clear what you are asking and what your inputs are. I think the policy here isn't to outright answer homework questions without seeing evidence of some attempts. Someone can provide a better answer if you can update your question.

---
Tags: bedtools, data-download, coverage, refseq
---


# Question
Title: What to use to determine SNPs in a population? How to find significant SNP variation between generations?
I have a population in the hundreds and did an Illumina pool-seq (ended up with ~100bp reads for the whole population).

I have the reference genome for the species and now need to align all the population reads to the reference genome, figure out which are SNPs and finally determine which ones are changing from one generation to the next (about 20generations difference).

Right now I'm trying to use popoolation2 https://sourceforge.net/p/popoolation2/wiki/Main/ but it seems to be a little old as it has not been updated in ~10 years.

In Galaxy, I couldn't find any pipelines dedicated to determine population SNP variation through generation.

So two questions:

* what software to use to determine significant variation in SNPs in a population (in poolSeq)?
* can you suggest me a pipeline for this? if one exists.

# Answer
> 1 votes
First focus on calling SNPs accurately. See the BWA, SAMtools, BCFtools pipeline (as described here) or GATK4 (as described here).

Once you have the \[BV\]CFs, then just code up the analysis to see changes over generations.

---
Tags: snp, phylogenetics
---


# Question
Title: building complex drug-dna for AMBER software
I will appreciate if you can please clarify some of my doubts about drug-DNA complex.

I want to study the drug-DNA simulation using AMBER. I did go through all the tutorials video on youtube but everything is about protein-ligand using GROMACS. I can also use GROMACS if that gives me the results.

My DNA sequences very small (12-mer) and I want to to intercalate it with a drug molecule. I generated the DNA sequence in Chimera. My questions are as follows.

1. Is it possible to insert my drug at specific sequence and then do the binding energy calculation using AMBER? (1) How can I insert it and make a complex? (2) If I have to create a complex in AMBER how can I do that---reference/tutorials link will be great.
2. The publications I read, they downloaded the PDB file with the drug molecule and extracted the crystallographic information where the drug is attached and then they extended the sequence from there using NAB tools in AMBER. Can somebody please explain how can we do that.

Thank you so much for your time.

1. I want to intercalate the drug molecule at the desired sequence I want and calculate the binding energy.

# Answer
> the tutorials video on youtube

Gromacs has great written documentation, so I'd suggest referring to that directly.

> I can also use GROMACS if that gives me the results.

Gromacs works fine with a DNA-ligand system. There are different Amber forcefields, proteins are generally used with `AMBER99SB-ILDN` simply because it has the longest name and is one of the inbuilt ones. Technically, the Amber forcefield tailored towards DNA is `OL15` (recommended here) and `gaff2` is for a generalised ff for ligands. These aren't natively within Gromacs, but been ported to Gromacs.

Whereas pi-pi iterations are somewhat covered by the Lenard-Jones term, the orientation will be driven not only by dipole, but also quadrupole moments, which aren't covered by classical mechanics forcefields. Likewise S-pi and cation-pi interactions will not be covered. Most DNA intercalating small molecules have multiple fused rings, e.g. doxorubicin, that pi stack with the DNA. **Consequently, a simulation under a classical mechanics forcefield is likely going to be gibberish.**

> Is it possible to insert my drug at specific sequence

See your previous question and the easiest option of dragging it in PyMOL

> Can somebody please explain how can we do that.

Read the AmberTools21 documentation.
> 1 votes

---
Tags: phylogenetics, docking, cancer, simulated-data
---


title: Perl convert GFF3 file

text:
# Question
Title: Perl convert GFF3 file
I have the following GFF3 file which I would like to convert to ZFF file with the below script:

```
$ head -n 20 traingenes.gff3
##gff-version 3
# gffread v0.12.7
# gffread -E traingenes.gtf -o traingenes.gff3
NbLab330C00 GeneMark.hmm3   mRNA    74501   76501   .   +   .   ID=1_t;geneID=1_g
NbLab330C00 GeneMark.hmm3   CDS 74501   74512   .   +   0   Parent=1_t
NbLab330C00 GeneMark.hmm3   CDS 75568   75999   .   +   0   Parent=1_t
NbLab330C00 GeneMark.hmm3   CDS 76151   76501   .   +   0   Parent=1_t
NbLab330C00 GeneMark.hmm3   mRNA    76535   77079   .   +   .   ID=2_t;geneID=2_g
NbLab330C00 GeneMark.hmm3   CDS 76535   76639   .   +   0   Parent=2_t
NbLab330C00 GeneMark.hmm3   CDS 76702   77079   .   +   0   Parent=2_t
NbLab330C00 GeneMark.hmm3   mRNA    93763   100703  .   -   .   ID=3_t;geneID=3_g
NbLab330C00 GeneMark.hmm3   CDS 93763   93837   .   -   0   Parent=3_t
NbLab330C00 GeneMark.hmm3   CDS 93915   94031   .   -   0   Parent=3_t
NbLab330C00 GeneMark.hmm3   CDS 94351   94430   .   -   2   Parent=3_t
NbLab330C00 GeneMark.hmm3   CDS 95483   95589   .   -   1   Parent=3_t
NbLab330C00 GeneMark.hmm3   CDS 95697   95746   .   -   0   Parent=3_t
NbLab330C00 GeneMark.hmm3   CDS 100464  100703  .   -   0   Parent=3_t

```

```
#!/usr/bin/env perl
use strict;
#
# source: https://biowize.wordpress.com/2012/06/01/training-the-snap-ab-initio-gene-predictor/
#
my @exons;
my $gene_count = 0;
my $current_seq = "";
while (my $line = <STDIN>) {
  if ($line =~ m/^###/) {
    flush(\@exons);
    next;
  }
  my @fields = split(/\t/, $line);
  if ($fields[2] eq "mRNA") {
    flush(\@exons);
  } elsif ($fields[2] eq "exon") {
    if ($fields[0] ne $current_seq) {
      $current_seq = $fields[0];
      printf(">%s\n", $current_seq);
    }
    push(@exons, \@fields);
  }
}
flush();

sub flush
{
  my $num_exons = scalar(@exons);
  return if($num_exons == 0);

  my $group = sprintf("%s.%d", $exons[0]->[0], $gene_count);
  $gene_count++;

  if ($num_exons == 1) {
my($start, $end) = ($exons[0]->[3], $exons[0]->[4]);
if ($exons[0]->[6] eq "-") {
  ($start, $end) = ($exons[0]->[4], $exons[0]->[3]);
}
printf("Esngl\t%lu\t%lu\t%s\n", $start, $end, $group);
  } else {
@exons = reverse(@exons) if($exons[0]->[6] eq "-");
for (my $i = 0; $i < $num_exons; $i++) {
  my $exon_type = "Exon";
  if ($i == 0) {
    $exon_type = "Einit";
  } elsif ($i == $num_exons - 1) {
    $exon_type = "Eterm";
      }

      my($start, $end) = ($exons[$i]->[3], $exons[$i]->[4]);
      if ($exons[0]->[6] eq "-") {
    ($start, $end) = ($exons[$i]->[4], $exons[$i]->[3]);
      }

      printf("%s\t%lu\t%lu\t%s\n", $exon_type, $start, $end, $group);
    }
  }
  @exons = ();
}

```

However, when I run it like `perl gff3_to_zff.pl < traingenes.gff3 >> traingenes.ann`  an empty file get generated.

What did I miss?

Thank you in advance,

Michal

# Answer
Your file looks like it is `CDS` feature based, not `exon` feature based, which *seems* to be what the script expects.

There is a tool to convert CDS based GFF to exon based that I've used successfully in the past...
> 1 votes

# Answer
The easiest is to use `agat_convert_sp_gff2zff.pl` from AGAT
(**Edit** Actually you will have to first standardize your file because you do not have gene features, neither gene\_id or Parent atttributes in your mRNA features. AGAT will be confused. So, in your case, prior the *gff2zff* conversion you will have to run `agat_convert_sp_gxf2gxf.pl -c geneID --gff file.gff -o file_clean.gff`)

Otherwise you can use your script but you need to add the exon as suggested by @Dan.
To do so you can just replace the CDS by exon using a sed command, it should work.
You can also use `agat_convert_sp_gxf2gxf.pl` (with parameter `-c geneID`) it will add the gene and exon features.
> 1 votes

---
Tags: perl, gff3
---


# Question
Title: How to import a large amount of .bed, .gff, .vcf, .paf, .sam files into an SQL database?
Are there best practices to load different bioinformatics file formats such as VCF, BED, GFF, and SAM to SQL databases? I am wondering how people out there do that efficiently.

All of these three formats are tab-separated files, so basically the following should work. I feel weird about it since most people I know don't use MySQL to work with these files.

```
LOAD DATA LOCAL INFILE 'bed.bed'
INTO TABLE bed-file
FIELDS TERMINATED BY '\t'
LINES TERMINATED BY '\n'
IGNORE 1 ROWS (list of the columns)
SET creation_date  = STR_TO_DATE(@creation_date, '%m/%d/%y');

```

# Answer
The clodius tool puts tab-delimited data from a BED file into an rtree in a sqlite3 database:

https://github.com/higlass/clodius/blob/c98bb16ade93402fcea3b749d705c52ea165b609/clodius/cli/aggregate.py#L457-L777

Data are stored in this database at multiple "zoom levels", which allows for fast queries at wide genomic extents:

https://github.com/higlass/clodius/blob/develop/clodius/tiles/beddb.py

Aggregation (input) and querying (output) seem efficient in practice, but efficiency may depend on aspects of your specific dataset, i.e., whether you mean efficiency in execution time of `INSERT` or `SELECT` statements, or in space efficiency in the overall database size.

Other data can be converted to BED via `vcf2bed`, `sam2bed`, `gff2bed`, etc. and then aggregated. This is not the only way to do things, of course, but hopefully this can suggest directions to explore.
> 2 votes

# Answer
You might also take a look at this genomics-focused sqlite extension: https://github.com/mlin/GenomicSQLite

As described in the repository README:

> This SQLite3 loadable extension adds features to the ubiquitous embedded RDBMS supporting applications in genome bioinformatics:
>
> * genomic range indexing for overlap queries & joins
> * in-SQL utility functions, e.g. reverse-complement DNA, parse "chr1:2,345-6,789"
> * automatic streaming storage compression (also available standalone)
> * reading directly from HTTP(S) URLs (also available standalone)
> * pre-tuned settings for "big data"
> 1 votes

---
Tags: vcf, database, bed, sam, gff
---


# Question
Title: list comprehension and codon translation
I'm trying to create a python program that follows the codon translation and turn mRNA sequence into amino acids (see the dictionary). The system has to return an error of the string is not a multiplier of 3, and the the function should not produce any more amino acids once a stop codon (keys with a value of X in the dictionary) is found, and should not display 'X' in the output. I do not know how to incorporate the stop codon to stop the programming entirely, as well as the multiplier of 3 aspect in the code. Can someone give me any advice? Thank you so much! example: protfromRNA('GCCAACGCACCCACUUUUUUUUAGAAA') should return ANAPTFF

```
def protfromRNA(Seq):
    GeneticCode = {'AAA':'K', 'AAC':'N', 'AAG':'K', 'AAU':'N',
               'ACA':'T', 'ACC':'T', 'ACG':'T', 'ACU':'T',
               'AGA':'R', 'AGC':'S', 'AGG':'R', 'AGU':'S',
               'AUA':'I', 'AUC':'I', 'AUG':'M', 'AUU':'I',
               'CAA':'Q', 'CAC':'H', 'CAG':'Q', 'CAU':'H',
               'CCA':'P', 'CCC':'P', 'CCG':'P', 'CCU':'P',
               'CGA':'R', 'CGC':'R', 'CGG':'R', 'CGU':'R',
               'CUA':'L', 'CUC':'L', 'CUG':'L', 'CUU':'L',
               'GAA':'E', 'GAC':'D', 'GAG':'E', 'GAU':'D',
               'GCA':'A', 'GCC':'A', 'GCG':'A', 'GCU':'A',
               'GGA':'G', 'GGC':'G', 'GGG':'G', 'GGU':'G',
               'GUA':'V', 'GUC':'V', 'GUG':'V', 'GUU':'V',
               'UAA':'X', 'UAC':'Y', 'UAG':'X', 'UAU':'Y',
               'UCA':'S', 'UCC':'S', 'UCG':'S', 'UCU':'S',
               'UGA':'X', 'UGC':'C', 'UGG':'W', 'UGU':'C',
               'UUA':'L', 'UUC':'F', 'UUG':'L', 'UUU':'F'}
    return ''.join(GeneticCode[i] if GeneticCode[i]!=X else '' for i in Seq)

```

# Answer
> 1 votes
You need to iterate over the sequence 3 characters at a time and you need to stop when you hit the stop codon (or before for an invalid codon, e.g. of length 2). Here is a memory efficient implementation using iterators (requires python \>= 3.8):

```
from itertools import islice

def protfromRNA(Seq):
    it = iter(Seq)
    lst = []
    while codon := "".join(islice(it, 3)):
        aa = GeneticCode.get(codon)
        if aa is None:
            raise ValueError(f"invalid codon: {codon}")
        if aa == "X":
            break
        lst.append(aa)
    return "".join(lst)

print(protfromRNA('GCCAACGCACCCACUUUUUUUUAGAAA'))
# ANAPTFF

```

---
Tags: python, codon
---


# Question
Title: Simulation of DNA sequences through substitution rates
I'm looking for a little bit of guidance.

My question is regarding the simulation of DNA sequences with a fix substitution rate.

The majority of the programs for simulating sequences use Continuous Markov Chain Models, with different instantaneous rate matrices that can be modified according to the free parameters of each model (JC69, HKY, GTR, etc..). For example, in HKY the free parameters are the transition/transversion (kappa) rate and nucleotide frequencies (pi).

In these programs (for example INDELible, Seq-Gen, Pyvolve), the substitution rate is expected to be represented by the branch length of the phylogenetic tree used to create the simulations.

INDELible paper.

*"..Rate matrices are rescaled by INDELible such that the branch lengths represent the expected number of substitutions per site (or the average expected number of substitutions per site under a heterogeneous-sites model)."* (P2 - Simulation of substitutions)

INDELible: a flexible simulator of biological sequence evolution

*"...Each branch length is assumed to denote the mean number of nucleotide substitutions per site that will be simulated along that branch..."* (P2 - Algorithm)

Seq-Gen: an application for the Monte Carlo simulation of DNA sequence evolution along phylogenetic trees

Let's say I want to simulate sequences with a mean substitution rate of 2.5 x 10^-7 per year. How would I prepare a tree that can accurately represent that substitution rate and be used in the simulation? In other words, how can the length be represented as a substitution rate accurately?

Any insight is appreciated!. Thanks in advance!

# Answer
## UPDATE based on comment

Simplest answer: just rescale your substitution rate from years to units of branch length. (Stop thinking in terms of years.)

## original answer

The short answer is that you have to calibrate the branch lengths of the tree to yield the substitution rates per year you want. (Leave the topology alone.)

Tree branch lengths are usually agnostic to chronological time (e.g. years). Instead it is very focused on whatever the model used to construct the tree is. This can be simple substitution counts (e.g. Jukes-Cantor), or a sophisticated model that takes into account a lot of information about relative rates of different changes, possibly learned from the dataset itself (e.g. GTR).

You yourself will have to decide what the conversion factor is between your branch lengths (presumably derived from such a substitution model or a distance metric) and units of chronological time. In some trees, it may already be calibrated to units of chronological time. But those will be rare.

So what I'd do is to make sure that the branch lengths of the tree are of a length that will, on average, yield the number of substitutions that you expect based on how long (in years) a given branch on the tree is supposed to be. For example, if you have a tree of primates, then you want root to tip distances to be on the order of tens of millions of years, probably. If you have a tree of all bacteria, you want the root-to-tip distances to be on the order of billions of years.

Note please that this is identical to re-scaling the substitution matrix to substitutions per year. What is crucial is that you know what you mean by "substitutions per year" and how you are defining that in terms of the tree.
> 2 votes

---
Tags: phylogeny, phylogenetics
---


# Question
Title: Use biopython to align SeqRecords stored in dict
I'd like to perform multiple alignments, where a gene from each sample was read in from fasta files. The fasta file represented one sample and had multiple genes. I have read in each sample fasta file and now have a dictionary of genes and their samples and sequences. Here is a small example.

The main key is the gene ID for another dictionary as the value. In the dictionary 'R??' is the sample name and key, followed by the annotation and sequence in a list.

```
{'PRELSG_MIT_v1:5894-6032': {'R47': ['not annotated',
                                     Seq('ATAAAAATAGTTACCATAGCTGTTGATGGATGCTTCGATCATAAAATATGTTTG...GGA')],
                             'R48': ['not annotated',
                                     Seq('ATAAAAATAGTTACCATAGCTGTTGATGGATGCTTCGATCATAAAATATGTTTG...GGA')],
                             'R49': ['not annotated',
                                     Seq('ATAAAAATAGTTACCATAGCTGTTGATGGATGCTTCGATCATAAAATATGTTTG...GGA')],
                             'R50': ['not annotated',
                                     Seq('ATAAAAATAGTTACCATAGCTGTTGATGGATGCTTCGATCATAAAATATGTTTG...GGA')],
                             'R51': ['not annotated',
                                     Seq('ATAAAAATAGTTACCATAGCTGTTGATGGATGCTTCGATCATAAAATATGTTTG...GGA')],
                             'R52': ['not annotated',
                                     Seq('ATAAAAATAGTTACCATAGCTGTTGATGGATGCTTCGATCATAAAATATGTTTG...GGA')]},
 'exon_PRELSG_MIT00300.1-E1': {'R47': ['cytochrome b',
                                       Seq('TCCAGATAATGCTATTACAGTAGATAGATATGCTACTCCTTTACATATTGTTCC...TTC')],
                               'R48': ['cytochrome b',
                                       Seq('TCCAGATAATGCTATTACAGTAGATAGATATGCTACTCCTTTACATATTGTTCC...TTC')],
                               'R49': ['cytochrome b',
                                       Seq('TCCAGATAATGCTATTACAGTAGATAGATATGCTACTCCTTTACATATTGTTCC...TTC')],
                               'R50': ['cytochrome b',
                                       Seq('TCCAGATAATGCTATTACAGTAGATAGATATGCTACTCCTTTACATATTGTTCC...TTC')],
                               'R51': ['cytochrome b',
                                       Seq('TCCAGATAATGCTATTACAGTAGATAGATATGCTACTCCTTTACATATTGTTCC...TTC')],
                               'R52': ['cytochrome b',
                                       Seq('TCCAGATAATGCTATTACAGTAGATAGATATGCTACTCCTTTACATATTGTTCC...TTC')]}}

```

I'd like to align each gene and write out those alignments to file without any intermediate files written to disk. How do I do that?

# Answer
> 0 votes
I don't believe this is possible. So I'm answering this to close it, since I don't have enough reputation to close it myself.

# Answer
> 0 votes
This BioPython document should get you started.

A specific example is given here.

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


# Question
Title: Snakemake Fastqc: "Multiple run or shell keywords in rule run_fastqc."
I am trying to check the quality of RNA-Seq data from Illumina using fastqc in snakemake in a conda environment. I get the error "Multiple run or shell keywords in rule run\_fastqc".

Snakefile:

```
rule run_fastqc:
    input:
        "/[RNA fastq file path]/well01_S001_L001_R1_001.fastq.gz"
    output:
        "./[out path]/well01_S001_L001_R1_001.fastqc.html",
        "./[out path]/well01_S001_L001_R1_001.fastqc.zip"
    wrapper:
        "0.79.0/bio/fastqc"
    shell:
        "fastqc {input} --outdir=./[out path]/"

```

I've also tried:

```
...
    shell:
        "fastqc {input} -o ./[out path]/"
...
    shell:
        "fastqc -o ./[out path]/ {input}"
...
    shell:
        "fastqc {input}"
        "fastqc -o ./[out path]/"
...
    shell:
        "fastqc {input} > -o ./[out path]/"

```

Code:

```
snakemake --cores 1 well01_S001_L001_R1_001.fastqc.html well01_S001_L001_R1_001.fastqc.zip

```

Output:

```
SyntaxError in line 9 of /[path to Snakefile]/Snakefile:
Multiple run or shell keywords in rule run_fastqc. (Snakefile, line 9)

```

When I run fastqc alone (not in snakemake), it works and I get the two output files. How can I get snakemake to work?

Thank you for your time and help!

# Answer
From a quick look at the manual, it looks like you are supposed to use just the `wrapper`, not combine the `wrapper` with a `shell`. Can you show me an example of such combined usage?

The error seems to suggest that the wrapper already has a `shell` directive that you `shell` conflicts with.
> 1 votes

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


# Question
Title: Ensembl Variant Effect Predictor (VEP) issue during execution
I got this error during execution, pls guide me.

vep --cache -i xxxx.vcf -o xxxxx_new.vcf

Possible precedence issue with control flow operator at /software/xxxx_pkgs/envs/ensemble-vep/lib/site_perl/5.26.2/Bio/DB/IndexedBase.pm line 805.
DBI connect('host=ensembldb.ensembl.org;port=3306','anonymous',...) failed: Can't connect to MySQL server on 'ensembldb.ensembl.org' (110) at /software/xxxx_pkgs/envs/ensemble-vep/share/ensembl-vep-88.9-1/Bio/EnsEMBL/Registry.pm line 1771.

-------------------- EXCEPTION --------------------
MSG: Cannot connect to the Ensembl MySQL server at ensembldb.ensembl.org:3306; check your settings & DBI error message: Can't connect to MySQL server on 'ensembldb.ensembl.org' (110)
STACK Bio::EnsEMBL::Registry::load_registry_from_db /software/xxxx_pkgs/envs/ensemble-vep/share/ensembl-vep-88.9-1/Bio/EnsEMBL/Registry.pm:1771
STACK Bio::EnsEMBL::VEP::BaseVEP::registry /software/xxxx_pkgs/envs/ensemble-vep/share/ensembl-vep-88.9-1/modules/Bio/EnsEMBL/VEP/BaseVEP.pm:163
STACK Bio::EnsEMBL::VEP::BaseRunner::setup_db_connection /software/xxxx_pkgs/envs/ensemble-vep/share/ensembl-vep-88.9-1/modules/Bio/EnsEMBL/VEP/BaseRunner.pm:80
STACK Bio::EnsEMBL::VEP::Runner::init /software/xxxx_pkgs/envs/ensemble-vep/share/ensembl-vep-88.9-1/modules/Bio/EnsEMBL/VEP/Runner.pm:72
STACK Bio::EnsEMBL::VEP::Runner::run /software/xxxx_pkgs/envs/ensemble-vep/share/ensembl-vep-88.9-1/modules/Bio/EnsEMBL/VEP/Runner.pm:111
STACK toplevel /software//envs/ensemble-vep/bin/vep:209
Ensembl API version = 88
```

# Answer
> 2 votes
You might prefer to run VEP in offline mode:

> Using `--cache` (without `--offline`) uses the local cache on disk to fetch most annotations, but allows database connections for some features (see cache limitations)

If you need HGVS names in your output, you'll also need to supply a FASTA file using `--fasta`:

> --offline
>
> Enable offline mode. No database connections will be made, and a cache file or GFF/GTF file is required for annotation. Add --refseq to use the refseq cache (if installed). Not used by default.

> --fasta \[file|dir\]
>
> Specify a FASTA file or a directory containing FASTA files to use to look up reference sequence. The first time you run VEP with this parameter an index will be built which can take a few minutes. This is required if fetching HGVS annotations (--hgvs) or checking reference sequences (--check\_ref) in offline mode (--offline), and optional with some performance increase in cache mode (--cache). See documentation for more details. Not used by default.

---
Tags: vcf, sequence-annotation, ensembl, errors
---


title: R help in Jupyter notebook?

text:
# Question
Title: R help in Jupyter notebook?
I'm using this Jupyter notebook to work on R... How do I do the equivalent of `?someFun`?

When I run that in the notebook, I get:

```
Warning message in seq_len(head.end.idx):
“first element used of 'length.out' argument”
ERROR while rich displaying an object: Error in seq_len(head.end.idx): argument must be coercible to non-negative integer

Traceback:
1. FUN(X[[i]], ...)
...

```

### Update

When I execute `docker run jupyter/r-notebook` I notice that it **does** pop up the widget (after throwing the error). Not it isn't a popup window, but a panel at the bottom of the workbook. My notebook is 'Trusted', which I think is significant, but I don't know how to turn it off.

When I execute `docker run -e JUPYTER_ENABLE_LAB=yes jupyter/r-notebook` I get the error and the widget **does not** pop up. I'm not sure how to set or inspect the JS trust in that environment. I suspect it's related. I don't see any blocked popups in the browser.

Strange that even when it works (the widget opens at the bottom of the browser) I still get an error in the output.

# Answer
If you launch like `docker run -p 8888:8888 jupyter/r-notebook` you will get the help page like a widget and a warning in your cell output. Check your browser, maybe is blocking pop-up elements.
> 3 votes

# Answer
I met the same problem. R after v4.1.1 seem to have problems with `repr_html`. I fixed the problem by following the quotes.

> This tiny commit solves a long standing problem while trying to read help from inside jupyterlab.
>
> In order not to wait for the PR to be merged, you may just install from @renkun-ken 's fork: `devtools::install_github("renkun-ken/repr", ref = "fix-html-filter")`

#### Reference

https://github.com/IRkernel/repr/pull/144#event-5589475824

https://github.com/IRkernel/repr/issues/143
> 1 votes

# Answer
I think error is not related to `?`. This is what I see when I try `?seq`:
> 0 votes

---
Tags: r, jupyter
---