# Question
Title: less than 80% of your list has mapped to your chosen identifier type
I am trying to do GO analysis. However, when I run DAVID, I am getting this error:

> You are either not sure which identifier type your list contains, or less than 80% of your list has mapped to your chosen identifier type. Please use the Gene Conversion Tool to determine the identifier type. In addition, when I am trying to add a column with gene names next to esemble ids in R in the file with differential expression genes, many genes say NA under gene name column.

I have used STAR and Salmon for the alignment,and DESEQ2 for differential expression analysis.
I have used premade indexes for both tools.
Thank you in advance!

Some of the ENSEMBL codes I used:

> ENSG00000249200 ENSG00000239920 ENSG00000242295 ENSG00000235763 ENSG00000226933 ENSG00000283632 ENSG00000235838 ENSG00000232429 ENSG00000281076 ENSG00000271042 ENSG00000229119 ENSG00000253384 ENSG00000285645 ENSG00000240627 ENSG00000213590 ENSG00000270739 ENSG00000240201 ENSG00000203531 ENSG00000271626 ENSG00000250316 ENSG00000218713 ENSG00000231548 ENSG00000231049 ENSG00000286185 ENSG00000214535 ENSG00000212961 ENSG00000219747 ENSG00000259100 ENSG00000231767 ENSG00000243175 ENSG00000236804 ENSG00000240634 ENSG00000270773 ENSG00000230267

# Answer
> 2 votes
DAVID last received an update in May 2016 (https://david.ncifcrf.gov/content.jsp?file=update.html), so any additions to the Gene Ontology database since then won't be present. I checked the examples you posted in Biomart and found that only 2 have GO annotations, so you're likely not missing an immense number of genes. Having said that, you'd be best served to use a more frequently updated tool than DAVID. Around 2015 there was an outcry because it hadn't been updated in 6 years and the database was so out of date that people were getting questionable results. Have a look at one of the many Bioconductor packages for convenient free solutions that are better maintained.

---
Tags: rna-seq, go-enrichment
---


# Question
Title: How to sort multiple FASTA files based on their content?
I have around 10,000 FASTA files of Influenza A virus.

These files contains sequences of each of the 8 segments of the viral genome and I want to separate each genome segment into a different folder and subdivide each directory into different files based on their geographic region of isolation.

In each FASTA file for each segment,the first line has the segment number, for example.

> KM368312.1 Influenza A virus (A/swine/Shandong/01/2009(H1N1)) segment 3 polymerase PA (PA) and PA-X protein (PA-X) genes, complete cds

To be clear I want

* file1.txt has segment 1
* file2.txt has segment 2 etc ...

I want to ensure all segment 1 sequences are grouped into one folder, and each file is grouped according to its geographic origin. The geographic groupings are mirrored for all 8 segements and each placed into their own directory.

# Answer
The general idea is:

1. iterate over all files
    1. read the first line of each file, extract the `segment` bit
    2. move the file based on that

In Python, this can be done in many ways. Here’s one way:

```
import os
import re

def find_segment(id):
    return re.search(r'segment (\d)', id).group(1)

for segment in range(8):
    os.mkdir(f'segment-{segment + 1}')

for entry in os.scandir():
    if not entry.is_file():
        continue

    with open(entry.path, 'r') as file:
        segment = find_segment(next(file))
    os.rename(entry.path, f'segment-{segment}/{entry.path}')

```

Note the absence of error handling — this code assumes that every operation succeeds; empty files, or files without “`segment X`” in their FASTA ID, will trigger an error, as will the usual IO failures (missing permissions, etc.).

---

But most people would probably use a shell scripting language instead of Python to solve this. For example, here’s something equivalent in Bash (untested):

```
mkdir segment-{1..8}

for file in *.fa *.fasta; do
    segment=$(sed 's/.*segment \([[:digit:]]\).*/\1/; q' "$file")
    if [[ -n $segment ]]; then
    mv "$file" "segment-$segment"
    done
done

```
> 4 votes

# Answer
One "out there" approach you could take is to create a fasta file containing your regions of interest (as separate sequences), and then map your fasta files to that reference fasta file, generating BAM/SAM output. After that, position-sort the BAM file, and you've now got a content-sorted list of sequences.

The process to do this would be something like this:

```
minimap2 -a regions.fasta input1.fasta input2.fasta ... | samtools sort > mapped.bam
samtools fastq mapped.bam > sorted.fasta

```
> 1 votes

---
Tags: python, covid-19, virus
---


# Question
Title: How do I download the mitochondrial haplogroup datasets for human genetics online?
I seem to have landed on the mitomap.org site, but I don't know what to make of it or what do with it / how to get the genomes onto my computer. It sounds like the genomes are stored in GenBank, but that mitomap simply lists which of the GenBank genomes are human mitochondrial haplogroup genomes.

First of all, I can't find a simple list of links or list of names of the haplogroups the mitomap has... Can you show me where I can get or download that? They say:

> This brings our total number of FL sequences to 51,673, and the number of CR sequences to 74,660. Our SNVs now total 19,227.

I don't know what FL or CR (Control Region sequences?) means, but I am expecting to see lists of genomes or something. What am I to do with this information? Where do I find this data? Is this it? If I click on one of the results, I see this. It has path `nuccore/MT742594.1`, can I find this in some NCBI FTP server somewhere? I don't see anything of that folder structure in the GenBank FTP server...

Basically my question is:

1. Where is the complete list of human mitochondrial haplogroup genomes?
2. Where can I download all of those human mitochondrial haplogroup genomes so I can get them onto my computer?

# Answer
> 1 votes
Hmmm, I am not experienced enough in biology or bioinformatics to truly understand the question, but I am guessing you want to download the FASTA file or genebank file containing the sequence of the accession number you listed above. There are 2 ways to do this -

1. Using your browser, if you want to download a small sample of files. On the right hand side of your screen you can see the "send to" button on this link, click on it and choose the what's appropriate for you and then click on create file.

`https://www.ncbi.nlm.nih.gov/nuccore/MT742594.1`

2. If you familiar with Biopython and python programming you can download multiple files at once in whatever format you want, there is an easy script to do this and has been taken from Biopython's documentation on biopython.org

import os

from Bio import SeqIO

from Bio import Entrez

Entrez.email = "A.N.Other@example.com" # Always tell NCBI who you are

filename = "MT742594.1.gb"

if not os.path.isfile(filename):

```
 # Downloading...

 net_handle = Entrez.efetch(

     db="nucleotide", id="MG762674", rettype="gb", retmode="text"

 )
 out_handle = open(filename, "w")

 out_handle.write(net_handle.read())

 out_handle.close()

 net_handle.close()

 print("Saved")

```

print("Parsing...")

record = SeqIO.read(filename, "gb")

print(record)

You can make a list of the accession numbers of the sequences you want to download and pass the list in the id, if you are curious as to how the code works visit biopython.org .

---
Tags: genome, data-download, public-dataset
---


# Question
Title: DESEQ2:Error in rownames, what is the problem?
I am using DESEQ2 for RNAseq analysis but i do not understand why I get this error.

```
> dds <- DESeqDataSetFromMatrix(countData = data,
+                               colData = labeldata,
+                               design =~ subtype, tidy = TRUE)

**Error** in `.rowNamesDF<-`(x, value = value) :
  duplicate 'row.names' are not allowed

```

I really do not know what is the issue. Could you please help?

Thanks

# Answer
The hint is here:

`duplicate 'row.names' are not allowed`

`rownames(data)` should return a vector of only unique values, since it doesn't make sense to have the same gene/transcript in two different rows of your count matrix. How did you generate it?

What does `which(duplicated(rownames(data)))` return? Maybe that'll give you a hint on what went wrong. In general, if you use `salmon` or `kallisto` for counting, you should use the excellent `tximport` package to get the counts into R, then create a DESeq dataset with `DESeqDataSetFromTximport()`.

EDITS:

The replies to my answer prompted me to re-read your post and I saw that you had set `tidy=TRUE`, which means the duplicates are actually in your first column of your `data`.
> 3 votes

---
Tags: r, rna-seq
---


# Question
Title: Adding multiple dataframes to compare against genotypes - DATE/LOCATION - geotemporal incidences of SARS2 genotypes
I'm trying to reuse some code I have previously used to add extra parameters. In the initial code, it was used to give a 1/0 counter if a specific mutation was present (first code below)

```
import re
import os
import pandas as pd
import logging
import sys
import csv

root = "C:"

os.chdir(root)

lineages = os.listdir('Lineages')
mutations = os.listdir('Mutations')

for lineage in lineages:

    df = pd.read_csv('Lineages/' + lineage).dropna()
    df.columns = ["lineage"]
    df['lineage'] = [x.split('|')[0].split('_')[-1]
                     for x in df['lineage']]

    for mutation in mutations:

        mutation_df = pd.read_csv('Mutations/' + mutation).dropna()
        mutation_header = mutation.replace(".txt", "")
        mutation_df.columns = [mutation_header]
        mutation_df[mutation_header] = [x.split('_')[-1].split('|')[0]
                          for x in mutation_df[mutation_header]]

        mutation_set = set(mutation_df[mutation_header])

        df[mutation_header] = [1 if x in mutation_set else 0
                               for x in df['lineage']]

    df.to_csv('Results9/' + lineage.replace(".txt", "") + '.csv',
                 header = True,
                 index = False)

```

However, I now wish to add extra parameters such as date of sample and location of sample. The first code produces this type of counter system:

```
lineage,A222V,D614G,E484Q,E780Q,G476S,L18F,N439K,S477,S477N,T478I,V483A
417941,0,1,0,0,0,0,0,0,0,0,0
417942,0,1,0,0,0,0,0,0,0,0,0
417944,0,1,0,0,0,0,0,0,0,0,0
417946,0,1,0,0,0,0,0,0,0,0,0
417947,0,1,0,0,0,0,0,0,0,0,0
417948,0,1,0,0,0,0,0,0,0,0,0
417950,0,0,0,0,0,0,0,0,0,0,0
417955,0,1,0,0,0,0,0,0,0,0,0
431011,0,1,0,0,0,0,0,0,0,0,0

```

This code is fine and working, but if I try and add the second parameter (date/location) with the following code:

```
import re
import os
import pandas as pd
import logging
import sys
import csv

root = "C:/Users/jones/Downloads/SARSCOV2/"

os.chdir(root)

lineages = os.listdir('Lineages')
mutations = os.listdir('Mutations')
dates = os.listdir('Dates')
locations = os.listdir('Locations')

for date in dates:

    date_df = pd.read_csv('Dates/' + date).dropna()
    date_df.columns = ["Date"]
    date_df['Date'] = [x.split('|')[-1].split('/')[-1]
                       for x in date_df['Date']]

    for location in locations:

        location_df = pd.read_csv('Locations/' + location).dropna()
        location_df.columns = ["Location"]
        location_df['Location'] = [x.split('/')[1].split('_')[0]
                                   for x in location_df['Location']]

        for lineage in lineages:

            df = pd.read_csv('Lineages/' + lineage).dropna()
            df.columns = ["lineage"]
            df['lineage'] = [x.split('_')[-1].split('|')[0]
                             for x in df['lineage']]

            for mutation in mutations:

                mutation_df = pd.read_csv('Mutations/' + mutation).dropna()
                mutation_header = mutation.replace(".txt", "")
                mutation_df.columns = [mutation_header]
                mutation_df[mutation_header] = [x.split('_')[-1].split('|')[0]
                                                for x in mutation_df[mutation_header]]

                mutation_set = set(mutation_df[mutation_header])

                df[mutation_header] = [1 if x in mutation_set else 0
                                       for x in df['lineage']]

            df.to_csv('Results_3/' + lineage.replace(".txt", "") + '.csv',
                 header = True,
                 index = False)

```

I retrieve this data, where the counting is maintaibned, but the new parameters (DATE/LOCATION) are not added

```
lineage,A222V,D614G,E484Q,E780Q,G476S,L18F,N439K,S477,S477N,T478I,V483A
417941,0,1,0,0,0,0,0,0,0,0,0
417942,0,1,0,0,0,0,0,0,0,0,0
417944,0,1,0,0,0,0,0,0,0,0,0
417946,0,1,0,0,0,0,0,0,0,0,0
417947,0,1,0,0,0,0,0,0,0,0,0
417948,0,1,0,0,0,0,0,0,0,0,0
417950,0,0,0,0,0,0,0,0,0,0,0
417955,0,1,0,0,0,0,0,0,0,0,0
431011,0,1,0,0,0,0,0,0,0,0,0

```

In short, my end goal is to:

* Add Lineage ID|Location of Sample|DD/MM/YYY|Mutations| -- as the headers
* still maintain the mutation counter
* Use the mutation counter as a guide to determine where and when the genotypes (mutation combinations) arose.

The files where I am retrieving the parameters looks like this also:

```

hCoV-19/Singapore/4/2020|EPI_ISL_410535|2020-02-03

hCoV-19/USA/WA13-UW9/2020|EPI_ISL_413601|2020-03-02

hCoV-19/USA/WA-UW142/2020|EPI_ISL_416680|2020-03-11

hCoV-19/USA/WA-UW143/2020|EPI_ISL_416681|2020-03-11

hCoV-19/USA/WA-UW144/2020|EPI_ISL_416682|2020-03-09

hCoV-19/USA/WA-UW145/2020|EPI_ISL_416683|2020-03-15

hCoV-19/USA/WA-UW146/2020|EPI_ISL_416684|2020-03-14

hCoV-19/USA/WA-UW147/2020|EPI_ISL_416685|2020-03-15

hCoV-19/USA/WA-UW148/2020|EPI_ISL_416686|2020-03-14

```

A very laborious, and poorly explained question; I apologise in advance!

# Answer
> 2 votes
Just to clarify, when I asked you to share some of your files, I meant to share them in a way that I can reproduce your issue from start to finish. You read from at least four different files in a quadruple for-loop and I'm pretty sure your goal could be achieved much more straightforwardly than your current implementation.

Without a reproducible example, I can offer a best-guess on what's happening: you create `date_df` and `location_df`, but after you instantiate them, you never do anything with the objects. Your final write happens from `df`, so you'll need to add those columns, e.g.:

```
df["Date"] = date_df["Date"][0]
df["Location"] = location_df["Location"][0]

```

The above is just an example, you'll need to make sure the *correct* value gets assigned to your `df` object instead of simply the first one.

---
Tags: covid-19, sars-cov-2, python, pandas, dataframe
---


# Question
Title: I have really skewed RNA-seq data, what's the best way to normalise it? Preferably in R!
My data frame compares the RNA-seq reads from many genes in different tissues. The reads look as follows.

I tried using log to make it better but still looks pretty skewed.

Are there any better transformations/functions for this in R?

This is part of an assignment I have, I have used in the past edge software and agree it will do a better job but in this work I'm not meant to use any of the pipelines.

# Answer
> 6 votes
Welcome to RNA-seq. This is not at all unusual, as RNA-seq raw data are not normally distributed. People typically approximate it via the Negative Binomial distribution. I suggest you check existing packages as Devon suggests in his comment for normalization and differential expression (or anything else). Typical choices are DESeq2, edgeR or limma, but there are others as well. This guide might help you get started. Please also use google to explore existing threads both here at SE and over at biostars.org. Most basic questions have already been discussed there and are informative to get started.

# Answer
> 2 votes
I like using DESeq2. There's a great document written by the developers about how to process gene count tables and look at differential expression:

https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html

The document also mentions other tools that have similar approaches:

> Other Bioconductor packages with similar aims are edgeR, limma, DSS, EBSeq, and baySeq.

# Answer
> 1 votes
Another approach, that does not replace the methods already mentions, but could be useful if you don't have raw counts to start with, or if you want to screen hundreds of genes, and you want to focus on the most interesting.

1. either `log transform` or `zscale` your expression values
2. do a boxplot to see how your gene is expressed in the different tissues
3. run a non-parametric test to compare the distributions. I usually use a `kruskal.test` with multiple comparisons and a Mann-Whitney U test (`wilcox.test`) if I only have two groups.

This should give you an idea if the expression of a gene is different between tissue types.

# Answer
> 1 votes
Thanks for all your answers. Although the pipeline packages would have been useful I was trying to normalise the data before analysing it later and so I used the bestnormalize package which looked like it did a pretty good job.

---
Tags: r, rna-seq, homework
---


# Question
Title: Need help with Pairwise Alignment module in iterating over the alignment
```
from Bio import SeqIO
from Bio import Align

ref_seq_1 = SeqIO.read('C:/Users/King_Arthur/Desktop/ref_seq/segment 1/ref_seq_8.fasta','fasta')

seq1 = SeqIO.read('C:/Users/King_Arthur/Desktop/file/segment 1/Myfile_1 (1).fasta','fasta')

aligner = Align.PairwiseAligner()
aligner.mode = 'global'
aligner.match_score = 1
aligner.mismatch_score = -2
aligner.gap_score = -2

alignments = aligner.score(ref_seq_1.seq , seq1.seq)

print(alignments)
for alignment in sorted(alignments):
    print(alignment)

```

So this is my code and as you can see in the last section i am trying to iterate over my alignment but I am getting this error

```
TypeError: 'float' object is not iterable

```

I have tried various things like using `str()` but it gives some strange values and I also tried to read the source code by using the `inspect` module but I can't figure out the problem.

Any help would be really appreciated.

My final objective is to find out how many matches, mismatches and gaps are present in the final alignment using biopython.

if there is any other better way to do it in python please feel free to suggest.

# Answer
`aligner.score(ref_seq_1.seq , seq1.seq)` doesn't produce an alignment. It just gives you the score of the alignment as float value.

Instead you need to run:

`alignments = aligner.align(ref_seq_1.seq , seq1.seq)`

Now `alignments` is an `Bio.Align.PairwiseAlignments` you can iterate over `alignments` or get the score with `alignments.score` which gives you the same as aligner.score(ref\_seq\_1.seq , seq1.seq)\`

If you need more information, you could check the documentation.
> 2 votes

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


title: Nonbonded interactions on the GPU

text:
# Question
Title: Nonbonded interactions on the GPU
i'm using this command

```
gmx mdrun -deffnm nvt0 -nb gpu

```

but got

```
Nonbonded interactions on the GPU were requested with -nb gpu, but the GROMACS
binary has been built without GPU support. Either run without selecting GPU
options or recompile GROMACS with GPU support enabled

```

I've Nvidia rtx2060 with cuda 11.1 on ubnutu 18.4.5

# Answer
As the error tells you: `the GROMACS binary has been built without GPU support`. Did you install gromacs with `apt-get`? Usually these kind of programs are shipped with a `basic` built, with a lot of features disabled. If you did built gromacs from source, then cuda was not found during configuration. I would try building gromacs from source and:

1. be sure cuda is in your path and check that `nvcc` and the drivers have a matching version
2. run cmake and check the output, to see if the video card is recognised
3. run cmake with `-DGMX_GPU=ON`

Note that it seems that gromacs doesn't support multiple gpu accelerations( e.g. it is either `openCL` or `cuda`, you can't have both)
> 4 votes

---
Tags: python, gromacs
---


# Question
Title: gmx mdrun error: "Error in user input: Invalid input values In option s"
I'm new to bioinformatics and gromacs I tried to run the code and got this error

```
Program:     gmx mdrun, version 2020
Source file: src/gromacs/options/options.cpp (line 179)
Function:    void gmx::internal::OptionSectionImpl::finish()

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

```

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

```
mdrun -h

```

I got

```
mdrun: command not found

```

# Answer
If the command line says "command not found", it means that the binary is not executable or not found in any of the paths mentioned in the `$PATH` environmental variable. You can check for example this tutorial for details on how it works and how to resolve problems related to this variable.

Which probably means, that you need to go a step back, to the installation instructions of gromacs and try to install the `mdrun` binary. However, if the `mdrun` was not installed successfully, it's very likely that more stuff of the package will be missing, so it might be worth reinstalling it all.

P.S. It's not obvious to me if getting `mdrun` it will resolve the first posted problem (I have never worked with gromacs). If you would share the link to the research you have done, that would help others to help you.
> 2 votes

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

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

Side tip: run `gmx mdrun -h -hidden` to see hidden options.
> 1 votes

---
Tags: gromacs
---


# Question
Title: What are the coordinates of the extended major histocompatibility complex (xMHC) in GRCh38?
Is there an accepted boundary of the *extended* major histocompatibility complex (xMHC) in human? Specifically, I am interested in the boundary coordinates for Genome Reference Consortium Human Build 38 (GRCh38 / hg38).

I'd like to use these coordinates to query MyGene.info to retrieve a list of genes in the xMHC. For example, King et al 2019 used the following definition:

> we excluded drugs with nonhuman or xMHC gene targets. We considered xMHC to include HIST1H2AA and KIFC1 and all genes between them (Chromosome 6 25.7 Mb-33.4 Mb).

But I was curious whether there was an authoritative definition of the xMHC region in humans that would specify down the the base pair rather than just megabase pair.

# Answer
Here is one definition of the classical MHC region: https://www.ncbi.nlm.nih.gov/grc/human/regions/MHC?asm=GRCh38.p13

which defines it as the 5mb region chr6:28510120-33480577 in GRCh38 coordinates.

The "extended MHC" defined here \[1\] is 7.2 mb and encompasses an additional "extended Class I" region upstream of the region given above. The xMHC is defined in \[1\] as the regions between *HIST1H2AA* and *RPL12P1* genes, which would give you chr6:25726063-33400644 in Grch38. This is very similar to your HIST1H2AA - KIFC1 definition (chr6:25726063-33410226). Of course if you're trying to reproduce a result, probably better to use the original paper's definition.

\[1\] https://www.nature.com/articles/nrg1489
> 4 votes

---
Tags: gene, genomics, human-genome, human
---


# Question
Title: Interpreting contig alignments to a reference genome
I have applied two de novo genome assembly tools to data from the run SRR12707453, corresponding to a phage (I downloaded the data and have no relation with the authors of the study).

Using rnaSPAdes I get many long contigs (\>50.000 bp). When analyzing the contigs using Quast against a reference genome (MN813694.1) most long contigs are discarded and the seventh-longest one is the only one selected. If I blast the other contigs I find correspondence with other phages, which makes me conclude that many phages were present in the same probe.

This seventh contig has 47686 bp length, while the reference has 47713 bp length. But they are not equivalently arranged. Apparently, the genome was circularized and then it was broken around the position 32400 in the reference. Also, the sequence does not match between positions 45000 and 45951 (in the reference), and a little misassembly seems to occur around position 46300.

See here some screen captures for the Icarus visualization obtained with Quast:

In contrast, using Quast with the contigs obtained from the other assembly tool, I get some differences. There is also a contig with almost the same size as the reference genome (47573 bp) that has exactly the same structure than the contig obtained with rnaSPAdes: circularized, broken around position 32400, does not match between position 45000 and 45951, and small misassemblie aronund position 46300. But there are also some contigs in which there is no break around position 32400.

May I interpret that these contigs indicate that the phage genome is sometimes circularized and sometimes broken at position ~32400 and that's why some contigs don't find a sequence break at that position?

Am I interpreting anything wrong? Any other conclusion can be obtained from the data?

You can (temporaly) download here the Quast results.

# Answer
> 1 votes
some additional information to your request.

https://ftp.ncbi.nlm.nih.gov/genomes/Viruses/FamilyPhylogeneticTree/Graph/SRR12707453.MN813694.1.gfa

it contain assembled by reference graph using ncbi tool saute(skesa) which will be public later.

install bandage app and load this graph https://rrwick.github.io/Bandage/

picture looks intriguing,I would say , run contain many errors located in one position. Reason? I don’t really know.

---
Tags: sequence-analysis, assembly, reference-genome
---


title: Find CRISPR PAM-sites with python

text:
# Question
Title: Find CRISPR PAM-sites with python
I am a college student and just started bioinformatics. I am trying to write a script (not for school) that finds all potential NGG and TTTN protospacer adjacent motif sequences in a genome string. I wrote the following script:

```
TargetFile = open('file.fasta')
readTF = TargetFile.read()
useTF = "".join(readTF[66:].split())

#no_of_basepairs = len(useTF)

sites = {"ngg_strand" : [],
         "ngg_antistrand" : [],     #location counted from reverse direction
         "tttn_strand" : [],
         "tttn_antistrand" : []}    #location counted from reverse direction

for i in range(0, len(useTF)-2):
    gg_mer = useTF[i:i + 2]
    ttt_mer = useTF[i:i + 3]

    if gg_mer == "GG":
        sites["ngg_strand"].append(i)
    elif gg_mer == "CC":
        #look for reverse complement, record location on antistrand from reverse direction
        sites["ngg_antistrand"].append(len(useTF)-i-2)
    elif ttt_mer == "TTT":
        sites["tttn_strand"].append(i+1)
    elif ttt_mer == "AAA":
        # look for reverse complement, record location on antistrand from reverse direction
        sites["tttn_antistrand"].append(len(useTF)-i-2)

print(len(sites["ngg_strand"]))
print(useTF.count("GG"))

```

At the bottom I used the .count() method to compare results. I noticed that the method returns a smaller number than my code. I believe this is due to the .count() method not counting overlap, whereas my code does. I.e. in the string "NNNGGGGNNN" my code would count 3 occurences of the "GG" substring, whereas the .count() method returns 2 occurences.

My Questions:

1. If I want to find all possible PAM sites, which way of counting is correct? I believe my way would be correct, since as far as I know the RNA molecule doesn't advance in "triplet-steps". However, I am very much new to this and might be missing something.
2. If the reverse complement is found, is the location usually recorded from the same direction as on the strand, or in the reverse direction (as I implemented in the code)?
3. Like I said, I am new to this so if you see any other glaring problems I would appreciate any and all help. :-)

Thank you!

# Answer
The easiest solution is `finditer`

```
use re

# Compile you search term where [A|T|G|C] means A or G or C or T
p = re.compile("[A|T|G|C]GG")

iterator = p.finditer(useTF)
for match in iterator:
    print(match.span())

```

\** Example Output\**

```
(0, 2)
(22, 24)
(29, 31)

```

Easy does it. Its worth looking at the documentation for `re` here, https://docs.python.org/3/howto/regex.html for alternative approaches using `match` or `search` below. `finditer` is definitely the solution in context.

---

I support your coding approach because it provides lots of flexibility. For example, you might want to target GG or TTT around a given locus and you can easily modify the code to accomodate this. Someone elses code would be less flexible.

In summary, I kinda like your substrings and do like the count method. I'm not a fan of your `==` in this instance and your C-style matching should be substituted for reg-ex which in Python is implemented via `re`.

The basic issues are,

* Biopython would be used via AlignIO or SeqIO. This will allow any given format to be analysed. Yours approach is okay but prone to bugs (below). \[probably a minor issue\]
* Your sequence is in some sort of weird format, like extended phylip format and it is not fasta (that wouldn't work). It needs verifying because removing the first 65 characters might remove bits of sequence \[probably a minor issue\]
* Chopping of -2 from the length ... I dunno I'd check out -1 \[minor\]
* C-style matching. I don't like this. I would definitely use `re`. One solution uses `match` or `search`, but the solution above uses `finditer`, Reg-ex will allow to search for NGG and TTTN, which you are not doing
* I think you are looping a locus of 3 and 4 bp and I'm not sure how the C-style `==` would match this because thats 2 and 3 bp respectively. Just thoughts. You will have 'edge effects' at the beginning and end of sequences, e.g. TTT 3' will match but that isn't TTTN.
* I like your `count` method, but I think `finditer` is better because it gives you location information
* The dictionary is sort of fine. Personally, I would dump the dictionary into a `pandas` dataframe This way the results are output as e.g. CSV, spreadsheet flat-files for your colleagues to view.
* the `if` and `elif` would be placed in a subroutine, this is good practice to convert to an object orientated approach.
> 1 votes

---
Tags: crispr, python
---


# Question
Title: Are sequencing reads written differently when inserts and indexes are sequenced seperately?
I was going over different sequencing library protocols and noticed that the adapter sequence can vary between protocols. In some adapter sequences the index sits between the primer binding site and the insert. In this case during sequencing the index will be sequenced and then the insert that comes right after it. In other cases, the adapter sequence is built so that the index is outside the primer binding site, i.e. index - primer binding site - insert (see two leftmost parts of the attached figure). Here the insert is sequenced using the Read1 primer and afterwards the index is sequenced in a separate reaction using the 17 index primer.

In this case are the index reads just added on to the ends of the insert reads in the FastQ files? Is there a way to tell from the FastQ files where the index was positioned within the adapter?

# Answer
The order of sequencing in a dual index flow cell is:

1. Read 1
2. Index 1 (I7)
3. Index 2 (I5)
4. Read 2

Between each step there is a switch of the primer used (Read1-primer -\> Index1-primer (I7) -\> Index2-primer (I5) -\> Read2-primer)
There are differences in how the I5 is dealt with between dual- and single-read flow cells as explained in the Illumina Indexed Sequencing Overview. In this document there is a lot more detailed explanation of the sequencing workflow that may be helpful to you.

But as explained by the other answers, none of this will matter in most cases. Even if you decide to output the Index sequences they will be exported as separate FASTQ files.
> 1 votes

# Answer
This is definitely a bioinformatics question! Understanding how sequencing libraries are created in the wet-lab and how they are sequenced is fundamental to bioinformatics.

To answer the question, the index sequences are not included in the `fastq` file output from Illumina, or other machines. Illumina sequencers just use the index sequence to demultiplex the complex mixture of molecules in the sequencing pool. If a small "dot" on the glass sequencing slide has some sequence, but the indices don't match what the technician specified for a library, then that read is not included in the fastq output of that library since it could be contamination of some sort.

The sequences in the fasta file theoretically all should be the insert from the library that you prepared, but sometimes the molecule being sequenced is short, and there is read-through of the adapter sequence. This is the only circumstance under which the i7 or i5 index sequences might be in the normal `fastq` output sequence line. This is not ideal, and this is why people perform adapter trimming.
> 0 votes

---
Tags: rna-seq, ngs, adapter, library
---


# Question
Title: Get protein names corresponding to PDB ID
I have a list of about 4000 PDB IDs and would like to get the actual names of the proteins (e.g. lactate dehydrogenase, cytochrome c). I tried the batch header section at the Protein Databank Download page but it refused to accept my PDB IDs in formats (xxxx or xxxx.pdb, individually or space-separated) that worked in an interactive search for the protein structure.

Any suggestions?

# Answer
> 5 votes
You can use one of the UniProt Protein APIs.

As you said you have your pdb entries in a text file line by line you can, like this `example.txt` containing:

```
1brr
4lzm
2dyi

```

Using the commandline, you can use a little script like this to download the name, if it is available for the given pdb entry.

```
while read line;
do
 curl -X GET --header 'Accept:application/json' "https://www.ebi.ac.uk/proteins/api/proteins/pdb:$line" |
 jq -r '.[].protein.recommendedName.fullName.value'  |
 sed "s/^/$line\t/" >> pdb_names.txt;
done < example.txt;

```

You need to have `curl`, `sed` and `jq` installed on your system.

This gives you following output in `pdb_names.txt`

```
1brr    Bacteriorhodopsin
4lzm    Endolysin
2dyi    Ribosome maturation factor RimM

```

---

Update

if you want to speed it up, you can run it with `parallel`.

```
parallel -j 4 'curl -X GET --header "Accept:application/json" "https://www.ebi.ac.uk/proteins/api/proteins/pdb:{}" |   jq -r ".[]. .protein.recommendedName.fullName.value"  |  sed "s/^/{}\t/" >> pdb_names_parallel.txt' :::: example.txt

```

With the `-j` option you call how many jobs should run in parallel. The limit of the UniProt API is 200 request per second per user.

---

Update 7. Nov 2020

To get another info beside the protein name, you need to know how the `JSON` response from UniProt looks like.

To get also the scientific name, you can run following command:

```
parallel -j 4 'curl -X GET --header "Accept:application/json" "https://www.ebi.ac.uk/proteins/api/proteins/pdb:{}" | jq -r ".[] | .protein.recommendedName.fullName.value + \" - \" + .organism.names[0].value"  |  sed "s/^/{}\t/" >> pdb_names_parallel.txt' :::: example.txt

```

As result you get this:

```
1brr    Bacteriorhodopsin - Halobacterium salinarum (strain ATCC 700922 / JCM 11081 / NRC-1)
4lzm    Endolysin - Enterobacteria phage T4
2dyi    Ribosome maturation factor RimM - Thermus thermophilus (strain HB8 / ATCC 27634 / DSM 579)

```

# Answer
> 2 votes
Assuming you can use `R`, have you tried with `biomaRt`? For example, using `2bhl` (my PhD lover :D)

```
library(biomaRt)
ensembl <- useMart("ensembl",dataset="hsapiens_gene_ensembl")

# get list of all available info
filters <- listFilters(ensembl)
attributes <- listAttributes(ensembl)

getBM(attributes=c('hgnc_symbol','ensembl_gene_id','entrezgene_id',
                'protein_id','description',"superfamily"),
      filters = 'pdb',
      values = "2bhl",
      mart = ensembl)

```

Returns

```
hgnc_symbol ensembl_gene_id entrezgene_id protein_id description
1         G6PD ENSG00000160211          2539   ADO22353 glucose-6-phosphate dehydrogenase [Source:HGNC Symbol;Acc:HGNC:4057]
2         G6PD ENSG00000160211          2539   CAA27309 glucose-6-phosphate dehydrogenase [Source:HGNC Symbol;Acc:HGNC:4057]
3         G6PD ENSG00000160211          2539   AAA63175 glucose-6-phosphate dehydrogenase [Source:HGNC Symbol;Acc:HGNC:4057]
4         G6PD ENSG00000160211          2539   AAA52500 glucose-6-phosphate dehydrogenase [Source:HGNC Symbol;Acc:HGNC:4057]
5         G6PD ENSG00000160211          2539   AAA52501 glucose-6-phosphate dehydrogenase [Source:HGNC Symbol;Acc:HGNC:4057]
6         G6PD ENSG00000160211          2539   CAA39089 glucose-6-phosphate dehydrogenase [Source:HGNC Symbol;Acc:HGNC:4057]
7 ...

```

I am sure if you look at all the available filters and attributes, you can pinpoint the ID you are looking for.

---
Tags: protein-structure
---


# Question
Title: Excess average estimated identical-by-descent in genotype data
(Cross-post with Biostars: https://www.biostars.org/p/470788/)

We have some genotype data that we are putting through quality control in PLINK 1.9. As part of this QC, we have limited the data to subjects of self-reported and genomic European ancestries, and have run basic QC steps, removing call rates below 95% for individuals and variants, and removing variants excessively out of Hardy-Weinberg equilibrium (\< 10^-10).

Pruning for LD and limiting to the autosomes, we have calculated estimated IBD between our participants and the estimated the average pairwise IBD for each participant. We have found that an unusual number of our participants (233/13049) show a high mean pairwise IBD with others (\>3SD from the cohort mean, \> 0.028 in this instance, maximum 0.075). Examining examples of the individuals in question shows that this is driven by low-level relatedness (~0.03-0.1) with many other individuals, rather than a small subset of relatives in the data. In addition, the relatedness in each case stems entirely from IBD1 (one shared variant) and not at all from IBD2 (two shared variants). Individual average IBD correlates negatively (-0.69) with the genome-wide inbreeding coefficient.

I have one idea about what could be causing this (partial cross-contamination of the samples), but I am interested to know if others have other ideas what would be the cause?

# Answer
> 1 votes
A clear negative association within an inbreeding coefficient suggests convergent evolution driven by shared selective pressure. There is an outside chance of purifying selection, or else balancing selection, but selection is the key word. The phenotype of the locus needs investigating, i.e. outside nearly neutral theory what is constitutes positive selection?

Allopatry caused by genetic isolation and genetic drift is your null hypothesis. However Darwinian positive selection is the alternative hypothesis.

I would be wary against suggesting contamination unless you've a clear suspicion this is occuring. I do accept that in human genetics its often difficult to avoid, particularly in ancient DNA type work. In this case you would need to try and control for it. The argument I would put against contamination is that 'genetic relatedness' is localised within the genome, whereas contamination would be a genome-wide phenomena.

Convergence is common within for example HLA genotypes (T-cell) loci between global populations.

You could perform formal tests for postive selection either within pop gen (the locus will associate with significant homozygosity) or phylogeny (point mutation).

---
Tags: gwas, quality-control, genotyping
---


title: Get raw genome sequence

text:
# Question
Title: Get raw genome sequence
I'm currently working on testing some matching algorithms for strings. I would like to do some tests on raw genomic data as I expect different results from random strings given the lower entropy. How can I get a whole human genome sequence? (just the raw basis)

I have tried websites like https://www.ncbi.nlm.nih.gov/ it seems very confusing and full of unnecessary information given that I just want to have the raw basis without annotations.

# Answer
You want FASTA files from here.

If you want the entirety of the genome, get the dna-toplevel file. Files with rm and sm have repeat sequences masked, the former by replacing the repeats with Ns (hard mask) the latter by designating repeats as lower case letters (soft masked).

See the README for more info.
> 3 votes

---
Tags: fasta, genome, phylogenetics, human-genome
---


# Question
Title: Gene expression Table to Expression Matrix converstion
I have an RNAseq gene expression file (Count data) as follows, I need to conduct a Differential gene expression analysis between Patients, for that, I need to have this file as a Matrix of Rows as Genes and columns as counts (samples), how best can I get this as matrix file?

Thanks a lot for your help, Stay safe

# Answer
> 0 votes
Your question stated differently is as follows:

> I have data in long format and need to change it to wide format. For that, I want to keep the `gene_raw`, use `patient_id` for the column names and `gene_code` for the row names. How can I do this easily?

Assuming you know how to read that into R as a data.frame that I'll call `d`:

```
library(tidyr)
library(dplyr)
reformatted = d %>% select(patient_id, gene_code, gene_raw) %>% spread(patient_id, gene_raw) %>% column_to_rownames(gene_code)

```

Or something along those lines.

---
Tags: rna-seq, gene-expression, format-conversion
---


# Question
Title: What is the fastest way to get the reverse complement of a DNA sequence in python?
I am writing a python script that requires a reverse complement function to be called on DNA strings of length 1 through around length 30. Line profiling programs indicate that my functions spend a lot of time getting the reverse complements, so I am looking to optimize.

What is the fastest way to get the reverse complement of a sequence in python? I am posting my skeleton program to test different implementations below with DNA string size 17 as an example.

```
#!/usr/bin/env python
import random
import timeit

global complement
complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'}

DNAlength = 17

#randomly generate 100k bases
int_to_basemap = {1: 'A', 2: 'C', 3: 'G', 4: 'T'}
num_strings = 500000
random.seed(90210)
DNAstrings = ["".join([int_to_basemap[random.randint(1,4)] for i in range(DNAlength)])
              for j in range(num_strings)]
#get an idea of what the DNAstrings look like
print(DNAstrings[0:5])

def reverse_complement_naive(seq):
    this_complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'}
    return "".join(this_complement.get(base, base) for base in reversed(seq))

def reverse_complement(seq):
    return "".join(complement.get(base, base) for base in reversed(seq))

tic=timeit.default_timer()
rcs = [reverse_complement_naive(seq) for seq in DNAstrings]
toc=timeit.default_timer()
baseline = toc - tic
namefunc = {"naive implementation": reverse_complement_naive,
            "global dict implementation": reverse_complement}

for function_name in namefunc:
    func = namefunc[function_name]
    tic=timeit.default_timer()
    rcs = [func(seq) for seq in DNAstrings]
    toc=timeit.default_timer()
    walltime = toc-tic
    print("""{}
       {:.5f}s total,
       {:.1f} strings per second
       {:.1f}% increase over baseline""".format(
           function_name,
           walltime,
           num_strings/walltime,
           100- ((walltime/baseline)*100)   ))

```

By the way, I get output like this. It varies by the call, of course!

```
naive implementation
       1.83880s total,
       271916.7 strings per second
       -0.7% increase over baseline
global dict implementation
       1.74645s total,
       286294.3 strings per second
       4.3% increase over baseline

```

Edit: Great answers, everyone! When I get a chance in a day or two I will add all of these to a test file for the final run. When I asked the question, I had not considered whether I would allow for cython or c extensions when selecting the final answer. What do you all think?

Edit 2: Here are the results of the final simulation with everyone's implementations. I am going to accept the highest scoring pure python code with no Cython/C. For my own sake I ended up using user172818's c implementation. If you feel like contributing to this in the future, check out the github page I made for this question.

```
the runtime of reverse complement implementations.
10000 strings and 250 repetitions
╔══════════════════════════════════════════════════════╗
║ name                       %inc   s total  str per s ║
╠══════════════════════════════════════════════════════╣
║ user172818 seqpy.c        93.7%  0.002344  4266961.4 ║
║ alexreynolds Cython (v2)  93.4%  0.002468  4051583.1 ║
║ alexreynolds Cython (v1)  90.4%  0.003596  2780512.1 ║
║ devonryan string          86.1%  0.005204  1921515.6 ║
║ jackaidley bytes          84.7%  0.005716  1749622.2 ║
║ jackaidley bytesstring    83.0%  0.006352  1574240.6 ║
║ global dict                5.4%  0.035330   283046.7 ║
║ revcomp_translateSO       45.9%  0.020202   494999.4 ║
║ string_replace            37.5%  0.023345   428364.9 ║
║ revcom from SO            28.0%  0.026904   371694.5 ║
║ naive (baseline)           1.5%  0.036804   271711.5 ║
║ lambda from SO           -39.9%  0.052246   191401.3 ║
║ biopython seq then rc    -32.0%  0.049293   202869.7 ║
╚══════════════════════════════════════════════════════╝

```

# Answer
> 23 votes
I don't know if it's the fastest, but the following provides an approximately 10x speed up over your functions:

```
import string
tab = string.maketrans("ACTG", "TGAC")

def reverse_complement_table(seq):
    return seq.translate(tab)[::-1]

```

The thing with hashing is that it adds a good bit of overhead for a replacement set this small.

For what it's worth, I added that to your code as "with a translation table" and here is what I got on my workstation:

```
global dict implementation
       1.37599s total,
       363374.8 strings per second
       3.3% increase over baseline
naive implementation
       1.44126s total,
       346919.4 strings per second
       -1.3% increase over baseline
with a translation table
       0.16780s total,
       2979755.6 strings per second
       88.2% increase over baseline

```

If you need python 3 rather than python 2, then substitute `tab = str.maketrans("ACTG", "TGAC")` for `tab = string.maketrans("ACTG", "TGAC")`, since `maketrans` is now a static method on the `str` type.

For those wondering, using biopython is slower for this (~50% slower than the naive implementation), presumably due to the overhead of converting the strings to `Seq` objects. If one were already reading sequences in using biopython, though, I wouldn't be surprised if the performance was much different.

# Answer
> 9 votes
Here's a Cython approach that might suggest a generic approach to speeding up Python work.

If you're manipulating (ASCII) character strings and performance is a design consideration, then C or Perl are probably preferred options to Python.

In any case, this Cython test uses Python 3.6.3:

```
$ python --version
Python 3.6.3 :: Anaconda custom (64-bit)

```

To install Cython, e.g.:

```
$ conda install -c anaconda cython

```

The Cython code below seems to offer about the same speed bump as the translation table — perhaps similar code is run under the hood of that.

Two files are needed, starting with `setup.py`:

```
from distutils.core import setup
from Cython.Build import cythonize

setup(
  name = 'Reverse complement C test',
  ext_modules = cythonize("revcomp_c.pyx"),
)

```

And then a second file called `revcomp_c.pyx`:

```
from libc.stdlib cimport malloc

cdef int seq_len = 17
cdef char *seq_dest = <char *>malloc(seq_len + 1)
seq_dest[seq_len] = '\0'

def reverse_complement_c_v1(str seq):
    cdef bytes py_bytes = seq.encode('ascii')
    cdef char *seq_src = py_bytes
    cdef int i = 0
    cdef int d = 0
    for i in range(seq_len):
        d = seq_len - i - 1
        if seq_src[i] == 'A':
            seq_dest[d] = 'T'
        elif seq_src[i] == 'G':
            seq_dest[d] = 'C'
        elif seq_src[i] == 'C':
            seq_dest[d] = 'G'
        elif seq_src[i] == 'T':
            seq_dest[d] = 'A'
    return seq_dest[:seq_len].decode('ascii')

```

This can be compiled into a Python module like so:

```
$ python setup.py build_ext --inplace

```

Then we can modify the test bench to include this Cython module and the relevant test method:

```
#!/usr/bin/env python
import random
import timeit
from revcomp_c import reverse_complement_c_v1

global complement
complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'}

DNAlength = 17

#randomly generate 100k bases
int_to_basemap = {1: 'A', 2: 'C', 3: 'G', 4: 'T'}
num_strings = 500000
random.seed(90210)
DNAstrings = ["".join([int_to_basemap[random.randint(1,4)] for i in range(DNAlength)])
              for j in range(num_strings)]
#get an idea of what the DNAstrings look like
print(DNAstrings[0:5])

def reverse_complement_naive(seq):
    this_complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'}
    return "".join(this_complement.get(base, base) for base in reversed(seq))

def reverse_complement(seq):
    return "".join(complement.get(base, base) for base in reversed(seq))

tic=timeit.default_timer()
rcs = [reverse_complement_naive(seq) for seq in DNAstrings]
toc=timeit.default_timer()
baseline = toc - tic
namefunc = {"naive implementation": reverse_complement_naive,
            "global dict implementation": reverse_complement,
            "Cython implementation (v1)": reverse_complement_c_v1}

for function_name in namefunc:
    func = namefunc[function_name]
    tic=timeit.default_timer()
    rcs = [func(seq) for seq in DNAstrings]
    toc=timeit.default_timer()
    walltime = toc-tic
    print("""{}
       {:.5f}s total,
       {:.1f} strings per second
       {:.1f}% increase over baseline""".format(
           function_name,
           walltime,
           num_strings/walltime,
           100- ((walltime/baseline)*100)   ))

```

A sample run:

```
$ python ./revcomp_test.py
['ACTGCAATAGACCACGC', 'CAGCTTGAGCCATTAAC', 'GGCCCAAGAGTTCGAAC', 'CGACTGTCTCGAATTGT', 'ATCCGCTATGCGCCGTC']
naive implementation
       1.73295s total,
       288525.0 strings per second
       -2.5% increase over baseline
global dict implementation
       1.70896s total,
       292575.5 strings per second
       -1.0% increase over baseline
Cython implementation (v1)
       0.21458s total,
       2330158.7 strings per second
       87.3% increase over baseline

```

One easy way to speed this up is to use a `static const unsigned char` array as an ASCII lookup table, which maps a residue directly to its complement. This would replace the nest of `if` statements and probably give a nice little boost (and it appears it does, making it among the best performers so far!).

Here is my fast implementation of a reverse complement function in C:

https://gist.github.com/alexpreynolds/4f75cab4350e9d937f4a

You might be able to use this directly in Python via the `subprocess` library. Outsourcing the reverse complement step to a utility written in C will almost always beat the best that Python can do, and you can do nice and important things like bounds checking etc. without losing much speed.

It's unclear how "pure" the answer needs to be, but making a system call from Python seems fair if you're processing strings and your goal is performance.

Another direction to take may be to look at multithreading, if you don't need ordered output. If you have many thousands of sequences stored in memory, you could split an array of sequences up into smaller arrays by use of offsets or array indices. Each thread would work on "rc"-ing sequences in its own piece of the array.

# Answer
> 7 votes
The most reliable and simplest way is probably using `Biopython`:

```
from Bio.Seq import Seq

def my_reverse_complement(seq):
    return Seq(seq).reverse_complement()

print(my_reverse_complement('ATGCGTA'))
# TACGCAT

```

As Devon has already said here using `Biopython` isn't as fast as the naive Python solution, and I also tested that shown here with ipython. However, this is because `Biopython`'s implementation, although similar to the naive approach, includes other features; it can reverse complement RNA as well as DNA and it will tell you if you're mixing DNA and RNA

Note that if you really want a fast way you could look at `Cython` or another python extension. As I edit this now, there are several nice answers taking this approach from user172818 and Alex Reynolds.

# Answer
> 5 votes
Another python extension but without cython. The source code is available at the bottom of this answer or from this gist. On Mac with Python3:

```
string
   0.38166s total,
   1310067.4 strings per second
   84.4% increase over baseline
seqpy
   0.20524s total,
   2436182.7 strings per second
   91.6% increase over baseline

```

On Linux with Python2 (seqpy is the first):

```
seqpy
   0.10960s total,
   4561842.5 strings per second
   93.3% increase over baseline
string
   0.17811s total,
   2807330.9 strings per second
   89.1% increase over baseline

```

---

This is file `setup.py`:

```
from setuptools import setup, Extension

seqpy_module = Extension('seqpy', sources = ["seqpy.c"])

setup(
    name = "seqpy",
    description = "seqpy",
    ext_modules = [seqpy_module])

```

This is file `seqpy.c`:

```
#include <stdlib.h>
#include <stdint.h>
#include <Python.h>

static unsigned char seq_comp_table[256] = {
      0,   1,   2,   3,   4,   5,   6,   7,   8,   9,  10,  11,  12,  13,  14,  15,
     16,  17,  18,  19,  20,  21,  22,  23,  24,  25,  26,  27,  28,  29,  30,  31,
     32,  33,  34,  35,  36,  37,  38,  39,  40,  41,  42,  43,  44,  45,  46,  47,
     48,  49,  50,  51,  52,  53,  54,  55,  56,  57,  58,  59,  60,  61,  62,  63,
     64, 'T', 'V', 'G', 'H', 'E', 'F', 'C', 'D', 'I', 'J', 'M', 'L', 'K', 'N', 'O',
    'P', 'Q', 'Y', 'S', 'A', 'A', 'B', 'W', 'X', 'R', 'Z',  91,  92,  93,  94,  95,
     64, 't', 'v', 'g', 'h', 'e', 'f', 'c', 'd', 'i', 'j', 'm', 'l', 'k', 'n', 'o',
    'p', 'q', 'y', 's', 'a', 'a', 'b', 'w', 'x', 'r', 'z', 123, 124, 125, 126, 127,
    128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143,
    144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159,
    160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175,
    176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191,
    192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207,
    208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223,
    224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239,
    240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255
};

static PyObject *seqpy_revcomp(PyObject *self, PyObject *args)
{
    PyObject *r;
    char *seq, *rev;
    int i, len;
    PyArg_ParseTuple(args, "s#", &seq, &len);
    rev = (char*)malloc(len);
    for (i = 0; i < len; ++i)
        rev[len - i - 1] = seq_comp_table[(uint8_t)seq[i]];
    r = Py_BuildValue("s#", rev, len);
    free(rev);
    return r;
}

static PyMethodDef seqpy_methods[] = {
    {"revcomp", seqpy_revcomp, METH_VARARGS, "Reverse complement a DNA sequence"},
    {NULL, NULL, 0, NULL}
};

#if PY_MAJOR_VERSION >= 3
static struct PyModuleDef seqpy_module = { PyModuleDef_HEAD_INIT, "seqpy", NULL, -1, seqpy_methods };
PyMODINIT_FUNC PyInit_seqpy(void) { return PyModule_Create(&seqpy_module); }
#else
PyMODINIT_FUNC initseqpy(void) { Py_InitModule3("seqpy", seqpy_methods, NULL); }
#endif

```

To compile and test:

```
python setup.py build_ext -i
python -c 'import seqpy; print(seqpy.revcomp("GGGTT"))'

```

# Answer
> 3 votes
**Use a `bytearray` instead of a `string` and then employ `maketrans` to translate**

You do not need the more advanced string encoding capabilities of `string` to store a string of bases, but you're still paying for it in performance. Devon Ryan's suggestion of `maketrans` is the huge improvement, 10x faster than your naive implementation. Using the same approach, but swapping everything out for `byte`s allows a further 40% speed improvement, however:

```
# DNA sequence generation also needs to be swapped to bytes
tab_b = bytes.maketrans(b"ACTG", b"TGAC")

def reverse_complement_bytes(seq):
    return seq.translate(tab_b)[::-1]

```

In my comparison:

```
Table implementation
       0.24339s total,
       2054315.7 strings per second
       82.5% increase over baseline
Bytes table implementation
       0.14985s total,
       3336755.9 strings per second
       89.2% increase over table

```

# Answer
> 3 votes
Here is a revision of my original Cython answer which incorporates my suggestion to use a `char` lookup array:

```
from libc.stdlib cimport malloc

cdef int seq_len = 17
cdef char *seq_dest = <char *>malloc(seq_len + 1)
seq_dest[seq_len] = '\0'

cdef char *basemap = [ '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0',   '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0',  'T', '\0',  'G', '\0', '\0', '\0',  'C', '\0', '\0', '\0', '\0', '\0', '\0',  'N', '\0', '\0', '\0', '\0', '\0',  'A',  'A', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0',  't', '\0',  'g', '\0', '\0', '\0',  'c', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0',  'a',  'a' ]

def reverse_complement_c_v2(str seq):
    cdef bytes py_bytes = seq.encode('UTF-8')
    cdef char *seq_src = py_bytes
    cdef int i = 0
    for i in range(seq_len):
        seq_dest[seq_len - i - 1] = basemap[<int>seq_src[i]]
    return seq_dest[:seq_len].decode('UTF-8')

```

In testing under OS X, Python 3.6.3:

```
$ python ./revcomp.py
['ACTGCAATAGACCACGC', 'CAGCTTGAGCCATTAAC', 'GGCCCAAGAGTTCGAAC', 'CGACTGTCTCGAATTGT', 'ATCCGCTATGCGCCGTC']
...
Cython implementation (v1)
       0.20664s total,
       2419626.4 strings per second
       88.3% increase over baseline
Cython implementation (v2)
       0.14170s total,
       3528668.5 strings per second
       92.0% increase over baseline

```

Using my lookup array approach ("v2") adds a very decent performance bump over using `if` blocks ("v1"), and you can keep everything as a Python string.

# Answer
> 2 votes
Since at least version `1.71` of biopython you can use `Bio.Seq.reverse_complement`, which also works on plain strings natively (no conversion to Seq objects). On my mac I get 800k strings converted with that implementation ("biopython just rc") when using the benchmark.

```
10000 strings and 250 repetitions
                         percent increase over baseline  seconds total strings per second
name
alexreynolds Cython (v2)                          94.0%       0.002491          4015188.7
user172818 seqpy.c                                93.2%       0.002833          3529944.4
alexreynolds Cython (v1)                          90.7%       0.003885          2573863.7
devonryan string                                  86.0%       0.005804          1722983.5
jackaidley bytes                                  85.7%       0.005959          1678085.3
jackaidley bytesstring                            82.1%       0.007459          1340693.2
biopython just rc                                 70.5%       0.012275           814638.6
revcomp_translateSO                               49.9%       0.020822           480254.6
revcom from SO                                    41.7%       0.024253           412320.2
string_replace                                    35.2%       0.026946           371107.3
lambda from SO                                    10.7%       0.037111           269464.0
global dict                                        0.0%       0.041556           240641.7
biopython seq then rc                            -30.9%       0.054430           183721.9
naive (baseline)                                  -0.9%       0.041962           238313.6

```

# Answer
> 0 votes
```
def revcomplement(s):
    s = s.upper()
    c = { "A":"T","C":"G","T":"A","G":"C"}
    co = str()
    for i in s:
        co = c[i] + co
    return co

```

---
Tags: python, phylogenetics, benchmarking
---


# Question
Title: How to find all WGS assemblies accessions of a species
### Some background

Similar to the OP of https://www.biostars.org/p/377840/, I would like to programmatically BLAST a sequence to a local database of all WGS assemblies.
Since this isn't feasible for the average biology lab server (correct me if I am wrong), I plan to use ncbi-acc-download to download all WGS assemblies of the species of interest (not a popular species like E. coli, so it should be feasible). Then, I will create a BLAST database for the downloaded assemblies and BLAST the sequence to it.

### My question

How can I find all WGS assemblies accessions of a species?

---

My current plan is to search the NCBI Assembly database using Entrez and a search term such as `"wgs"[Properties] AND txid1337[orgn:exp]`.
**EDIT: IIUC, this approach might miss some WGS assemblies.** See my answer.

I am worried (and thus ask for your help) this isn't the right approach because there seem to be at least 3 other places in which assemblies can be found:

# Answer
There seem to be WGS assemblies that can't be found in NCBI Assembly database, e.g.: https://www.ncbi.nlm.nih.gov/nuccore/1779902990. I guess that such assemblies also cannot be found in The assembly\_summary.txt files that are described in https://ftp.ncbi.nlm.nih.gov/genomes/README\_assembly\_summary.txt.

My current best guess is that each WGS assembly has a "WGS master record" in NCBI Nuccore database. To find all WGS assemblies of a taxon whose uid in NCBI Taxonomy database is `1337`, search NCBI Nuccore database using Entrez and the search term `"wgs master"[Properties] AND txid1337[orgn:exp]`.
> 2 votes

---
Tags: assembly, phylogenetics, wgs
---


# Question
Title: Filtering out reads from a reference (e.g. rRNA) using bowtie2
As I understand it, bowtie2 can easily be used to split reads into one of two groups:

* reads for which both of a pair align well to a reference (using e.g. `--al-conc-gz`)
* reads for which one or both of a pair do not align a reference (using e.g `--un-conc-gz`)

But I really want to split this second group into reads for which neither of a pair align to the reference.

bowtie2's report give a nice breakdown, but it's not clear to me how to extract specific members of these groups:

```
1000000 reads; of these:
  1000000 (100.00%) were paired; of these:
    958118 (95.81%) aligned concordantly 0 times
    1329 (0.13%) aligned concordantly exactly 1 time
    40553 (4.06%) aligned concordantly >1 times
    ----
    958118 pairs aligned concordantly 0 times; of these:
      179 (0.02%) aligned discordantly 1 time
    ----
    957939 pairs aligned 0 times concordantly or discordantly; of these:
      1915878 mates make up the pairs; of these:
        711518 (37.14%) aligned 0 times
        32134 (1.68%) aligned exactly 1 time
        1172226 (61.18%) aligned >1 times
64.42% overall alignment rate

```

In this case, I really want those that didn't align at all (i.e. the 35.58% not included in the `64.42% overall alignment rate` which is the same as the `37.14%` of the pairs that did not align concordantly).

Since in my case the reference is rRNA and many of those sequences are shorter than my reads, I'm not surprised that so there are so few concordant reads. I'm really interested in reads that aren't concordant and aren't discordant discordant either.

So, if I'm sticking to bowtie2, it looks like my only option is to take the resulting SAM/BAM file and filter what I want based on flags and then convert that back into FASTQ files.

Is there another option that I'm missing? It seems that the remaining output options relate to unpaired reads or SAM files.

I'm looking into other aligners now...

My goal is to go from FASTQ files to FASTQ files without extra steps, but I can go through a SAM/BAM file stage if I need to and even build a solution around that if needed.

# Answer
> 2 votes
If you want to remove rRNA, use tools that are made for this. We usually go for SortMeRNA, which can do exactly what you describe your goal is, *i.e.* go from FASTQ to FASTQ, splitting your reads into rRNA and non-rRNA.

In our experience, it is very accurate, if a bit slow at times.

# Answer
> 0 votes
1. Bam flags will indicate which reads have both pairs not mapping.
2. Nothing mapped. You need to figure out why, and you can look at just about anything out figure that out.

---
Tags: sequence-alignment, bowtie2
---


# Question
Title: Why does Coronavirus genome contain poly A region at the 3' end in GenBank?
If you look at the Coronavirus genome below, it contains poly-A region at the end. I know that the viral RNA contains a poly-A tail. But, the genome in GenBank is stored in cDNA format. That means the end should have poly-T instead of A. Can anybody help me to understand this, please ?

```
>MN988713.1 Wuhan seafood market pneumonia virus isolate 2019-nCoV/USA-IL1/2020, complete genome
ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATCTGTTCTCTAAA
--------------------------------------------------------------------
GATCGAGTGTACAGTGAACAATGCTAGGGAGAGCTGCCTATATGGAAGAGCCCTAATGTGTAAAATTAAT
TTTAGTAGTGCTATCCCCATGTGATTTTAATAGCTTCTTAGGAGAATGACAAAAAAAAAAAA

```

# Answer
> 0 votes
I figured out that it has to do with second strand of cDNA. Usually genbank sequences are equivalent to 5’ to 3’ strand of cDNA. This strand is generated using second strand synthesis enzyme.

# Answer
> 0 votes
No its definitely a poly-A tail for both the *in situ* virus and the cDNA sequence. It couldn't be a poly-T *in situ* because,

1. the Ts are uracils (Us) because its an RNA virus. Thus if you substituted all the 'T's for 'U's you'd recover the *in situ* genome
2. specifically its a positive single-stranded RNA virus and the cDNA isn't a reverse complement but the positive strand.

If you look at the poly-A together with the 5' capping then you get the idea its to assist translocation and transcription to the endoplasmic reticulum. The genome pulls a few tricks which make it a bit more complicated than that, but thats the gist.

---
Tags: genome, covid-19, sars-cov-2, coronavirus
---


title: gprofiler output ENSG to gene name

text:
# Question
Title: gprofiler output ENSG to gene name
I would like to ask if there is any tool or R script that translates Ensembl IDs to gene names in a gprofiler GO output. An example of this file is the following one:

# Answer
This is usually done using AnnotationDbi in combination with the annotation database for human org.Hs.eg.db.

```
library(org.Hs.eg.db)
genes <- c("ENSG00000102055", "ENSG00000232605", "ENSG00000136828",
           "ENSG00000159197", "ENSG00000201394", "ENSG00000278580",
           "ENSG00000213440", "ENSG00000166816", "ENSG00000275895",
           "ENSG00000105723")
select(org.Hs.eg.db, genes, c("ENSEMBL", "ALIAS"), "ENSEMBL")

```

This returns:

```
           ENSEMBL     ALIAS
1  ENSG00000102055       I-4
2  ENSG00000102055  PPP1R2P9
3  ENSG00000102055   PPP1R2C
4  ENSG00000232605      <NA>
5  ENSG00000136828   RALGEF2
6  ENSG00000136828  RALGPS1A
7  ENSG00000136828   RALGPS1
8  ENSG00000159197     ATFB4
9  ENSG00000159197      LQT5
10 ENSG00000159197      LQT6
11 ENSG00000159197     MIRP1
12 ENSG00000159197     KCNE2
13 ENSG00000201394   RN5S140
14 ENSG00000201394 RNA5SP140
15 ENSG00000278580      <NA>
16 ENSG00000213440      <NA>
17 ENSG00000166816     DLACD
18 ENSG00000166816       DLD
19 ENSG00000166816      LDHD
20 ENSG00000275895   U2AF1L5
21 ENSG00000105723     GSK3A

```

If you want to apply a 1:1 mapping of gene IDs to common names, you can use the "SYMBOL" column. In a `data.frame`, you can do something like this:

```
ensembl_to_symbol <- function(e) {
  hits <- select(org.Hs.eg.db, e, "SYMBOL", "ENSEMBL")$SYMBOL
  # Gene IDs without symbols should just stay as IDs
  no_hit <- is.na(hits)
  hits[no_hit] <- e[no_hit]
  hits
}

recode_symbols <- function(s) {
  e <- ensembl_to_symbol(unlist(strsplit(s, split = ",")))
  paste(e, collapse = ",")
}

df <- data.frame(
  genes = c(paste(genes[1:5], collapse = ","),
            paste(genes[6:10], collapse = ",")
  )
)

df$genes <- sapply(df$genes, recode_symbols)

```

This isn't the most vectorized (fast) solution, but it's probably easier to understand this way.
> 1 votes

---
Tags: gene-ontology
---


# Question
Title: Visualizing FindMarkers result in Seurat using Heatmap
I used FindMarkers to find the DEGs between two clusters in my dataset using Seurat. I am trying to visualize the outcome using a heatmap but I failed to write the command in R.

```
dput(cluster4_3.markers)
              p_val  avg_logFC pct.1 pct.2    p_val_adj
Vim    1.803527e-44 -1.1689785 0.996 0.999 3.458985e-40
Lgals7 1.479278e-38  1.7615232 0.993 0.552 2.837108e-34
Anxa2  1.797016e-37 -1.3971764 0.629 0.991 3.446497e-33
Tmsb10 6.936830e-31  0.7937881 1.000 1.000 1.330415e-26
Ptma   7.518439e-27  0.5622808 1.000 1.000 1.441961e-22
Jund   2.044576e-25 -0.9881766 0.993 0.999 3.921293e-21

```

When I tried:

```
DoHeatmap(cluster4_3.markers, features = NULL) +
scale_fill_gradientn(colors = rev(RColorBrewer::brewer.pal(n = 7, name = "RdYlBu")))
Error in UseMethod(generic = "DefaultAssay", object = object) :
  no applicable method for 'DefaultAssay' applied to an object of class "data.frame"

```

It didn't work. I was trying to replicate this command:

```
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
DoHeatmap(pbmc, features = top10$gene) + NoLegend() +
  scale_fill_gradientn(colors = rev(RColorBrewer::brewer.pal(n = 7, name = "RdYlBu")))

```

I tried to extract some info from the csv file but I couldn't determine which parameter is important (p\_val avg\_logFC pct.1 pct.2 p\_val\_adj). Also, I didn't know which parameter to use to use the group\_by function. That's one of my attempts:

```
clus4vs3 <- cluster4_3.markers$p_val %>%
    group_by(pct.1) %>%
    top_n(n = 10, wt = avg_logFC)
Error in UseMethod("group_by_") :
  no applicable method for 'group_by_' applied to an object of class "c('double', 'numeric')"

```

Now that I understand 'group\_by' only accepts a categorical variable, how can I write a command that displays the top20 differentially expressed genes between cluster 4 and 3 in a heatmap ?

# Answer
In your `DoHeatmap()` call, you do not provide features so the function does not know which genes/features to use for the heatmap.

In your last function call, you are trying to group based on a continuous variable `pct.1` whereas `group_by` expects a categorical variable.
> 0 votes

---
Tags: seurat, heatmap
---


# Question
Title: How to predict with the pre-trained DNABERT model?
I was curious to give DNA BERT a try. This is a BERT (Bidirectional Encoder Representations from Transformers) model that was trained on short (`k=3,4,5, or 6`) k-mers of DNA. Overall, It is exciting work since it bridges NLP approaches and genetics data.

I span up a pytroch3.6 machine on SageMaker with a GPU and tried to assess the prediction of its pre-trained DNABERT model:

Preparing the data:

```
!git clone https://github.com/jerryji1993/DNABERT
%cd DNABERT
!python -m pip install --editable .
%cd examples
!python -m pip install -r requirements.txt
!curl -L https://northwestern.app.box.com/shared/static/g8m974tr86h0pvnpymxq84f1yxlhnvbi --output a6-new-12w-0.zip
!unzip a6-new-12w-0.zip
!cp ./sample_data/ft/prom-core/6/dev.tsv /home/ec2-user/SageMaker/DNABERT/examples/6-new-12w-0

```
2. Running prediction based on the pre-trained data:

```
!python run_finetune.py \
   --model_type=dna \
   --tokenizer_name=dna6 \
   --model_name_or_path=/home/ec2-user/SageMaker/DNABERT/examples/6-new-12w-0 \
   --task_name=dnaprom \
   --do_predict \
   --data_dir=/home/ec2-user/SageMaker/DNABERT/examples/sample_data/ft/prom-core/6  \
   --max_seq_length=75 \
   --per_gpu_pred_batch_size=128   \
   --output_dir=/home/ec2-user/SageMaker/DNABERT/examples/6-new-12w-0/out \
   --predict_dir=/home/ec2-user/SageMaker/DNABERT/examples/result \
   --n_process=4

```

Output:

```
11/09/2020 04:35:57 - WARNING - __main__ -   Process rank: -1, device: cuda, n_gpu: 1, distributed training: False, 16-bits training: False
11/09/2020 04:35:57 - INFO - transformers.configuration_utils -   loading configuration file /home/ec2-user/SageMaker/DNABERT/examples/6-new-12w-0/config.json
11/09/2020 04:35:57 - INFO - transformers.configuration_utils -   Model config BertConfig {
  "architectures": [
    "BertForMaskedLM"
  ],
  "attention_probs_dropout_prob": 0.1,
  "bos_token_id": 0,
  "do_sample": false,
  "eos_token_ids": 0,
  "finetuning_task": "dnaprom",
  "hidden_act": "gelu",
  "hidden_dropout_prob": 0.1,
  "hidden_size": 768,
  "id2label": {
    "0": "LABEL_0",
    "1": "LABEL_1"
  },
  "initializer_range": 0.02,
  "intermediate_size": 3072,
  "is_decoder": false,
  "label2id": {
    "LABEL_0": 0,
    "LABEL_1": 1
  },
  "layer_norm_eps": 1e-12,
  "length_penalty": 1.0,
  "max_length": 20,
  "max_position_embeddings": 512,
  "model_type": "bert",
  "num_attention_heads": 12,
  "num_beams": 1,
  "num_hidden_layers": 12,
  "num_labels": 2,
  "num_return_sequences": 1,
  "num_rnn_layer": 1,
  "output_attentions": false,
  "output_hidden_states": false,
  "output_past": true,
  "pad_token_id": 0,
  "pruned_heads": {},
  "repetition_penalty": 1.0,
  "rnn": "lstm",
  "rnn_dropout": 0.0,
  "rnn_hidden": 768,
  "split": 10,
  "temperature": 1.0,
  "top_k": 50,
  "top_p": 1.0,
  "torchscript": false,
  "type_vocab_size": 2,
  "use_bfloat16": false,
  "vocab_size": 4101
}

============================================================
<class 'transformers.tokenization_dna.DNATokenizer'>
11/09/2020 04:35:57 - INFO - filelock -   Lock 140334164036072 acquired on /home/ec2-user/.cache/torch/transformers/ea1474aad40c1c8ed4e1cb7c11345ddda6df27a857fb29e1d4c901d9b900d32d.26f8bd5a32e49c2a8271a46950754a4a767726709b7741c68723bc1db840a87e.lock
11/09/2020 04:35:57 - INFO - transformers.file_utils -   https://raw.githubusercontent.com/jerryji1993/DNABERT/master/src/transformers/dnabert-config/bert-config-6/vocab.txt not found in cache or force_download set to True, downloading to /home/ec2-user/.cache/torch/transformers/tmpa9hkn2q1
Downloading: 28.7kB [00:00, 12.5MB/s]
11/09/2020 04:35:57 - INFO - transformers.file_utils -   storing https://raw.githubusercontent.com/jerryji1993/DNABERT/master/src/transformers/dnabert-config/bert-config-6/vocab.txt in cache at /home/ec2-user/.cache/torch/transformers/ea1474aad40c1c8ed4e1cb7c11345ddda6df27a857fb29e1d4c901d9b900d32d.26f8bd5a32e49c2a8271a46950754a4a767726709b7741c68723bc1db840a87e
11/09/2020 04:35:57 - INFO - transformers.file_utils -   creating metadata file for /home/ec2-user/.cache/torch/transformers/ea1474aad40c1c8ed4e1cb7c11345ddda6df27a857fb29e1d4c901d9b900d32d.26f8bd5a32e49c2a8271a46950754a4a767726709b7741c68723bc1db840a87e
11/09/2020 04:35:57 - INFO - filelock -   Lock 140334164036072 released on /home/ec2-user/.cache/torch/transformers/ea1474aad40c1c8ed4e1cb7c11345ddda6df27a857fb29e1d4c901d9b900d32d.26f8bd5a32e49c2a8271a46950754a4a767726709b7741c68723bc1db840a87e.lock
11/09/2020 04:35:57 - INFO - transformers.tokenization_utils -   loading file https://raw.githubusercontent.com/jerryji1993/DNABERT/master/src/transformers/dnabert-config/bert-config-6/vocab.txt from cache at /home/ec2-user/.cache/torch/transformers/ea1474aad40c1c8ed4e1cb7c11345ddda6df27a857fb29e1d4c901d9b900d32d.26f8bd5a32e49c2a8271a46950754a4a767726709b7741c68723bc1db840a87e
11/09/2020 04:35:57 - INFO - transformers.modeling_utils -   loading weights file /home/ec2-user/SageMaker/DNABERT/examples/6-new-12w-0/pytorch_model.bin

```

However, the `result` directory is empty and I am not sure what did I miss?

```
!ls -lh /home/ec2-user/SageMaker/DNABERT/examples/result
total 0

!ls -l  /home/ec2-user/SageMaker/DNABERT/examples/6-new-12w-0/out
total 0

```

---

**EDIT:** I have tried also to run this on another ec2 machine and apparently it gets a segmentation fault, so my conclusion is that it's unstable. I will train a RoBERTAa model or similar by myself:

SageMaker doesn't propagate the segmentation fault to the notebook which is another weird thing about it.

```
11/12/2020 14:01:18 - INFO - filelock -   Lock 139866003860056 acquired on /root/.cache/torch/transformers/ea1474aad40c1c8ed4e1cb7c11345ddda6df27a857fb29e1d4c901d9b900d32d.26f8bd5a32e49c2a8271a46950754a4a767726709b7741c68723bc1db840a87e.lock
11/12/2020 14:01:18 - INFO - transformers.file_utils -   https://raw.githubusercontent.com/jerryji1993/DNABERT/master/src/transformers/dnabert-config/bert-config-6/vocab.txt not found in cache or force_download set to True, downloading to /root/.cache/torch/transformers/tmpwpxfkxlm
Downloading: 28.7kB [00:00, 20.2MB/s]
11/12/2020 14:01:18 - INFO - transformers.file_utils -   storing https://raw.githubusercontent.com/jerryji1993/DNABERT/master/src/transformers/dnabert-config/bert-config-6/vocab.txt in cache at /root/.cache/torch/transformers/ea1474aad40c1c8ed4e1cb7c11345ddda6df27a857fb29e1d4c901d9b900d32d.26f8bd5a32e49c2a8271a46950754a4a767726709b7741c68723bc1db840a87e
11/12/2020 14:01:18 - INFO - transformers.file_utils -   creating metadata file for /root/.cache/torch/transformers/ea1474aad40c1c8ed4e1cb7c11345ddda6df27a857fb29e1d4c901d9b900d32d.26f8bd5a32e49c2a8271a46950754a4a767726709b7741c68723bc1db840a87e
11/12/2020 14:01:18 - INFO - filelock -   Lock 139866003860056 released on /root/.cache/torch/transformers/ea1474aad40c1c8ed4e1cb7c11345ddda6df27a857fb29e1d4c901d9b900d32d.26f8bd5a32e49c2a8271a46950754a4a767726709b7741c68723bc1db840a87e.lock
11/12/2020 14:01:18 - INFO - transformers.tokenization_utils -   loading file https://raw.githubusercontent.com/jerryji1993/DNABERT/master/src/transformers/dnabert-config/bert-config-6/vocab.txt from cache at /root/.cache/torch/transformers/ea1474aad40c1c8ed4e1cb7c11345ddda6df27a857fb29e1d4c901d9b900d32d.26f8bd5a32e49c2a8271a46950754a4a767726709b7741c68723bc1db840a87e
11/12/2020 14:01:18 - INFO - transformers.modeling_utils -   loading weights file /home/ec2-user/DNABERT/examples/6-new-12w-0/pytorch_model.bin
Segmentation fault

```

# Answer
I don't why this RNN k-mer Pytorch model is not producing an output. I do know that,

```
 "rnn_dropout": 0.0,

```

... is bad news for any ANN (artificial neural network) let alone LSTM RNN, which is a pretty complex model.

Firstly, you have to assign a meaningful drop out rate and 0.5 is okay.

If dropout rate = 0 as an output then the model hasn't worked. A dropout layer is essential to present overfitting. Its a foundational part of deep learning. No-one emperically knows quite why it works, there's a few approximate good guesses, but a dropout is essential.

---

The other thing to note is what does "n\_process" mean? If you are parallelising a GPU across 4-cores, okay its got a lot of layers but even so thats alot of cores. It ain't a CPU and I'd get it working on a single GPU core first.
> 1 votes

---
Tags: deep-learning, pytorch, python
---


# Question
Title: I want to write python program to cut a DNA sequence at an EcoRI restriction site and print the two fragments after cutting
Seq = AAAAAAAAAAATTTTTTTTTTTGAATTCCCCCCCCCCCGGGGGGGGGGG

Restriction site = GAATTC

I want to cut the sequence at GA/ATTC

I tried split method in python but It does not cut at GA/ATTC. Please help

# Answer
```
import Bio
from Bio.Seq import Seq
my_seq = Seq("aaaaaaaaaaatttttttttttgaattccccccccccggggggggggg")
from Bio.Restriction import *
print(EcoRI.search(my_seq))
print(EcoRI.catalyse(my_seq))

```

> Actully the EcoRI cuts the DNA after G in restriction site GAATTC.
> 3 votes

---
Tags: python
---


# Question
Title: How can atomic RMS differences / atomic RMS distributions help us understand the structure of a protein?
For example, I have a paper that says: "the atomic RMS distribution about the mean coordinate positions is 0.41 +/- 0.04 A for the backbone atoms and 0.82 +/- 0.04 A for all atoms"

And I have no idea what any of this means, I found several papers with this kind of information but none of them explain why this is relevant in order to understand the structure of a protein.

# Answer
> 2 votes
The "RMS distribution about the mean coordinate positions" could be restated as "certainty of the structure". In other words, the main backbone that defines the large-scale shape in the structure is expected to be accurate to within about 0.41 angstrom, while the atoms more on the periphery (or otherwise not forming the backbone of the structure) are less accurately solved. How important this is will depend on what you want to do with the structure. If you want to get an idea for what it looks like then the accuracy is less important, if you want to do docking then the more accurate the atoms around any sort of pocket/interaction surface are the better.

---
Tags: protein-structure
---


# Question
Title: Metapath2vec on Drug-ADR Heterogeneous Graph
I am new in this field and I am having some problems regarding a new project.

I built a graph using Drugbank Data connected to SIDER Adverse Reactions. I used Organ- level Terms to classify the ADRs in 25 different groups. I aim to use the algorithm Metapath2vec in order to cluster drugs based on the group of ADRs they lead to. I am currently using the python package for Machine Learning on graphs Stellargraph and the related implementation of Metapath2vec.

I chose Clustering as downstream task with the resulting node embeddings, e.g. DBSCAN but the results are not promising. I would like to have some clusters of drugs related to different adverse drug reaction. Since I am new in this field everything I am trying is based on scientific literature but I don't know if this is the right approach for my objective.

This is the code related to the metapath2vec algorithm:

```
walk_length = 100  # maximum length of a random walk to use throughout this notebook

# specify the metapath schemas as a list of lists of node types.
metapaths = [
    ["drug", "adr", "drug"],
    ["drug", "adr", "drug", "drug"],
    ["drug", "drug"],
    ["drug", "adr", "group_adr", "adr", "drug"],
    ]

# Create the random walker
rw = UniformRandomMetaPathWalk(graph)

walks = rw.run(
    nodes=list(graph.nodes()),  # root nodes
    length=walk_length,  # maximum length of a random walk
    n=1,  # number of random walks per root node
    metapaths=metapaths,  # the metapaths
)

from gensim.models import Word2Vec

model = Word2Vec(walks, size=128, window=5, min_count=0, sg=1, workers=2, iter=1)

```

Should I change approach? what other algorithm for representation learning can I use in order to reach my goal? what can I improve in the presented approach to have better node embeddings?

# Answer
Clearly the additional references do put Metapath2vec into good use. Its a many to many ML relationship. Most catagorization/classification particularly in deep learning is many to one. So thats the good thing about the technique you are using, you've 25 catagories, which is quite large for machine learning(ML) and large for deep learning. The technique appears to be half way between supervised and unsupervised learning. The random walk establishes a zone (circle) of catagorisation.

**NLP Background**

The situation where it would work in NLP (natural language processing) is for example a spell checker, where you have one correct word and loads of possible typos, some of which you not observed in training. So what will happen is the random walk will pick up the local vacinity around the 'correct spelling' by the vectorization (converting a word into numbers) of the incorrect spelling will land it close to 'true spelling' in a 2D vectorization plot. The training following the random walk says 'yep you belong with me'.

The papers given talk about being able to identify the context of a word within a sentence which is a very NLP thing, one vectorisation being 'bags of words'. So here the sentence represents one of many categories of information, i.e. each is the central node within the graph, the vectorisation of the sentence places of meaning closer to one 'hub' than another and again the local vacinity is identified by the random walk. In training this 'random walk vacinity' is identified, so it will algorithmically recreate the 'random walk locality' so that sentences that have not been trained can have their meaning identified.

The central theory of deep learning, which not quite true of other ML, is any algorithm can be recreated through training. So like it will sort of recreate a 'random walk zone'

**Difference with drug discovery**

The papers you have map drugs to genes, but it means you have to have a "true drug to gene" relationship. The problem I see here is what is 'true' for a drug-protein interaction? With a spell-checker there is one 'true spelling', with a sentence there is pretty much one 'type of meaning', like a movie was great, average, terrible based on a ratings score and stuff like gender classification. If its not 'true' the random walk will get it wrong. So its a big assumption. If you have a circle and pick a point on the edge of the circle then try and draw a circle around the original circle, you can't unless the circle is so large it risks including stuff that is not part of the original circle. The problem the papers is how do they know what is the 'centre' of the circle? Its an a priori assumption and you need really good biological rationale.

With your stuff I agree its easier because adverse drug reaction, some drugs could be fatal or toxic (traditional cancer drugs or parasite drugs), so you'd pick the worst case adverse reaction. That then centres the definition around which everything else is relative. If however your 'true catagory' turned out to be an outlier its not going to work, because the random walk only works if its starts bang in the centre of the 'catagory'. Random walks form spheres, in case you are not aware. I can see this working as an advanced 'unsupervised learning' technique.

One very real advantage is if a drug doesn't have an adverse reaction it MIGHT (I'm guessing) not give a result, in which case you are good to go - that a good drug. If that is what happened that would be amazing because its not how ML, or particularly DL works. Traditional ML/DL will force a classification even if it doesn't fit. There are methods to accept some sort of "half-classification" but I don't know how to do it. This is one of the reaons you need lots of things to train.

I am assuming that vectorisation of your drug is fine because there is a huge history of this in drug-discovery going right back to the origins of the field.

**Summary**

So yeah I get the idea, I see the flaws and the advantages. The final issue is that most ML requires pretty large data sets and DL require huge data sets and when you've got 25 catagories. For DL approaches the data set size might need to be huge unless you can do 'transfer learning' (backing your training onto a pre-trained model). If your method resulted in a zero classification i.e. "doesn't fit anything I've got" - that would be cool and that is the key advantage I can see in the 'random walk' approach which is not available using traditional methods.

As long as you are aware, you do not have to use this algorithm to get a classification. This final point what appears missing in the papers is using a 'traditional' method by comparison. The ML before deep learning, stuff like random forests, routinely takes loads of catagories. If you are using a non-deep learning ML you would traditional use several different approaches to justify your result.

---

My answer is a bit rambly but I think you'll get the gist. Peronally I would complement with a standard ML, but all methods have their flaws. Deep learning (DL) has a big history in drug discovery (alot of work was done via Tensorflow), but I'm not sure you've data set size for this.
> 1 votes

---
Tags: networks, machine-learning, deep-learning, pytorch
---


title: Why SeqIO.parse method isn't working?

text:
# Question
Title: Why SeqIO.parse method isn't working?
I am trying to follow the tutorial from from the Biopython website here and I am right at the beginning (2.4.1), where I am trying a *Simple FASTA parsing example* which is

```
from Bio import SeqIO
for seq_record in SeqIO.parse("ls_orchid.fasta", "fasta"):
    print(seq_record.id)
    print(repr(seq_record.seq))
    print(len(seq_record))

```

For an unknown reason, I get the following error:

```
---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
<ipython-input-19-c64128d7dcfd> in <module>
      1 from Bio import SeqIO
----> 2 for seq_record in SeqIO.parse("ls_orchid.fasta", "fasta"):
      3     print(seq_record.id)
      4     print(repr(seq_record.seq))
      5     print(len(seq_record))

~\Anaconda3\lib\site-packages\Bio\SeqIO\__init__.py in parse(handle, format, alphabet)
    605     iterator_generator = _FormatToIterator.get(format)
    606     if iterator_generator:
--> 607         return iterator_generator(handle)
    608     if format in AlignIO._FormatToIterator:
    609         # Use Bio.AlignIO to read in the alignments

~\Anaconda3\lib\site-packages\Bio\SeqIO\FastaIO.py in __init__(self, source, alphabet, title2ids)
    181             raise ValueError("The alphabet argument is no longer supported")
    182         self.title2ids = title2ids
--> 183         super().__init__(source, mode="t", fmt="Fasta")
    184
    185     def parse(self, handle):

~\Anaconda3\lib\site-packages\Bio\SeqIO\Interfaces.py in __init__(self, source, alphabet, mode, fmt)
     45             raise ValueError("The alphabet argument is no longer supported")
     46         try:
---> 47             self.stream = open(source, "r" + mode)
     48             self.should_close_stream = True
     49         except TypeError:  # not a path, assume we received a stream

FileNotFoundError: [Errno 2] No such file or directory: 'ls_orchid.fasta'

```

I am using a Jupyter Notebook.

# Answer
> 2 votes
Personally I think its late and we're all tired.

> "FileNotFoundError: \[Errno 2\] No such file or directory: 'ls\_orchid.fasta'"

The traceback is pointing to line 2, the SeqIO as you point out.

In my opinion, the script should be in the same directory as the fasta file. You know how to identify the path in bash (`pwd` or `echo $PWD`), although you must be using Windows (Anaconda has \ rather than /). Alternatively the file names has a typo. You know the

```
for seq_record in SeqIO.parse("~\path\ls_orchid.fasta", "fasta"): # where path is the dir(s) leading to ls_orchid.fasta, but obviously use / if its Linux.

```

I would assume you can alternatively dump the .ipynb file in the location where the fasta file is.

As personal choice I prefer a formal IDE such as 'Visual Studio Code' or "PyCharm". In my personal opinion, I am aware some will disagree, Jupyter notebook good is for running code on a remote commercial server for data science work. However, the traceback - which you clearly have activated - looks good. I can see that you're Anaconda3 installation looks good, you must be using Windows (which I've never worked with) and you're running Jupiter from Anaconda3 ... so cool.

---
Tags: python, phylogenetics, phylogeny, seqio
---


# Question
Title: Error processing run file with bcl2fastq - Unable to open '.../RunInfo.xml' file for reading
I'm attempting to demux a nextseq run which I've successfully done several times in the past on other runs. I uploaded the run from my PC to the server (Ubuntu 18.04.4 LTS) without issue. Then I ran:

```bcl2fastq --runfolder-dir /data/201109_NB551126_0128_AH5G5KBGXF --sample-sheet ./SampleSheet.csv  --no-lane-splitting --barcode-mismatches 0```

From which I got the error:

>
'/data/201109_NB551126_0128_AH5G5KBGXF/RunInfo.xml' file for
> reading```
> </code></pre>

However, the RunInfo.xml file is present and readable:

```
$ ls -l
total 96
drwxrwxrwx 10 user user 90112 Nov 11 08:08 201109_NB551126_0128_AH5G5KBGXF
-rwxrwxr-x  1 user user 1929 Nov 11 07:37 SampleSheet.csv
$
$ ls -l 201109_NB551126_0128_AH5G5KBGXF/RunInfo.xml
-rwxrwxr-x 1 user user 28510 Nov 11 08:08 201109_NB551126_0128_AH5G5KBGXF/RunInfo.xml
$
$ head 201109_NB551126_0128_AH5G5KBGXF/RunInfo.xml
<?xml version="1.0"?>
<RunInfo xmlns:xsd="http://www.w3.org/2001/XMLSchema" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" Version="4">
  <Run Id="201109_NB551126_0128_AH5G5KBGXF" Number="128">
    <Flowcell>H5G5KBGXF</Flowcell>
    <Instrument>NB551126</Instrument>
    <Date>201109</Date>
    <Reads>
      <Read Number="1" NumCycles="75" IsIndexedRead="N" />
      <Read Number="2" NumCycles="6" IsIndexedRead="Y" />
    </Reads>

```

Has anyone had this issue or know of a fix? Thanks!

Update: contents of run folder

<pre><code>total 2452
-rwxrwxr-x 1 stewart stewart   977 Nov 11 07:59 '2020-11-09__14_05_58_94F12309 GetBestZ Camera 1 Focus Red top surface.snr.csv'
....
-rwxrwxr-x 1 stewart stewart   244 Nov 11 07:59 '2020-11-09__14_08_54_FE05AB9F SafeParabola called by CalculateSignalStd.snr.csv'
-rwxrwxr-x 1 stewart stewart   253 Nov 11 07:59 '2020-11-09__14_08_54_FE0620BC SafeParabola called by FindPeakUsingParabolaFit.snr.csv'
....
drwxrwxrwx 2 stewart stewart  4096 Nov 11 07:59  Config
drwxrwxrwx 3 stewart stewart  4096 Nov 11 07:59  Data
drwxrwxrwx 3 stewart stewart  4096 Nov 11 08:07  Images
drwxrwxrwx 2 stewart stewart 12288 Nov 11 08:07  InstrumentAnalyticsLogs
drwxrwxrwx 2 stewart stewart  4096 Nov 11 08:08  InterOp
drwxrwxrwx 2 stewart stewart  4096 Nov 11 08:08  Logs
drwxrwxrwx 2 stewart stewart  4096 Nov 11 08:08  Recipe
-rwxrwxr-x 1 stewart stewart    48 Nov 11 08:08  RTAComplete.txt
-rwxrwxr-x 1 stewart stewart  6200 Nov 11 08:08  RTAConfiguration.xml
drwxrwxrwx 2 stewart stewart  4096 Nov 11 08:08  RTALogs
-rwxrwxr-x 1 stewart stewart    36 Nov 11 08:08  RTARead1Complete.txt
-rwxrwxr-x 1 stewart stewart    36 Nov 11 08:08  RTARead2Complete.txt
-rwxrwxr-x 1 stewart stewart   925 Nov 11 08:08  RunCompletionStatus.xml
-rwxrwxr-x 1 stewart stewart 28510 Nov 11 08:08  RunInfo.xml
-rwxrwxr-x 1 stewart stewart 26802 Nov 11 08:08  RunParameters.xml
-rwxrwxr-x 1 stewart stewart  1929 Nov 11 08:08  SampleSheet.csv```

# Answer
> 0 votes
I repro the error.

```
ls /ml

```

> ls: cannot access '/ml': No such file or directory

If run

```
bcl2fastq --runfolder-dir ./data/201109_NB551126_0128_AH5G5KBGXF --sample-sheet ./SampleSheet.csv  --no-lane-splitting --barcode-mismatches 0

```

I suspect it will work.

At present you are attempting to call a directory `/data`, which doesn't exit. If you type `ls /etc` that will work. The difference between `./data` and `/data` in a file structure is alot, but on a command line it appears minimal.

---

Wait ... if that doesn't work .. because I suspect the calling directory is `data` then the following will,

```
bcl2fastq --runfolder-dir ./201109_NB551126_0128_AH5G5KBGXF --sample-sheet ./SampleSheet.csv  --no-lane-splitting --barcode-mismatches 0

```

OR

```
bcl2fastq --runfolder-dir ./data/201109_NB551126_0128_AH5G5KBGXF --sample-sheet ./data/SampleSheet.csv  --no-lane-splitting --barcode-mismatches 0

```

It depends on the path where you are making the call. You could sort this out using `$PWD` or just `echo $PWD` then copy and paste the path name.

---
Tags: illumina, demultiplexing
---


title: two aligners combined results

text:
# Question
Title: two aligners combined results
I would like to ask if someone uses two different aligners to produce count matrices and then run the same alogrithm for DEG analysis, would it make sense to find the intersection of the DEGs in order to conclude to the "real" DEGS?

# Answer
This sort of procedure is done in tool papers to benchmark results compared to different aligners/mappers. See the salmon or kallisto papers for examples. Using the intersect of DE results from different mappers will end up yielding somewhat conservative results, not unlike using different differential expression packages and taking the intersect of their outputs. In practice you could just take a more conservative p-value cut-off and save the hassle of using multiple tools, the results would end up being similar.
> 1 votes

---
Tags: rna-seq
---