title: Primer-Blast Custom Database [error]

text:
# Question
Title: Primer-Blast Custom Database [error]
I'm using Primer-Blast to compare a PCR template against a custom database of reference sequences, using Primer-Blast's *Custom* database functionality.

However whenever I try to upload a FASTA file as the custom database and then run Primer-Blast, I get the following error:

`Exception error: Fasta Reader: sequence id ends with 80 valid nucleotide characters. Was the sequence accidentally placed in the definition line?`

No matter how I try to format my FASTA file, I get this error.

Does anyone know how to properly format a FASTA file to have it accepted as a custom database of reference sequences to compare against the PCR template? Or is this a bug in the software?

Here's an example FASTA file, to upload both for the template and the custom database.

```
>test1
ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATCTGTTCTCTAAACGAACTTTAA
AATCTGTGTGGCTGTCACTCGGCTGCATGCTTAGTGCACTCACGCAGTATAATTAATAACTAATTACTGTCGTTGACAGG
ACACGAGTAACTCGTCTATCTTCTGCAGGCTGCTTACGGTTTCGTCCGTGTTGCAGCCGATCATCAGCACATCTAGGTTT
CGTCCGGGTGTGACCGAAAGGTAAGATGGAGAGCCTTGTCCCTGGTTTCAACGAGAAAACACACGTCCAACTCAGTTTGC

```

# Answer
> 1 votes
I just copied your FASTA sequence into a file (on Mac OS) and uploaded it to Primer BLAST, the search did not give an error.

If the error is caused by the line endings, maybe you can fix your downloaded files with a tool like dos2unix.

# Update:

On a whim, I tried adding a description to the FASTA header:

```
>test1 description
ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATCTGTTCTCTAAACGAACTTTAA
AATCTGTGTGGCTGTCACTCGGCTGCATGCTTAGTGCACTCACGCAGTATAATTAATAACTAATTACTGTCGTTGACAGG
ACACGAGTAACTCGTCTATCTTCTGCAGGCTGCTTACGGTTTCGTCCGTGTTGCAGCCGATCATCAGCACATCTAGGTTT
CGTCCGGGTGTGACCGAAAGGTAAGATGGAGAGCCTTGTCCCTGGTTTCAACGAGAAAACACACGTCCAACTCAGTTTG

```

This time, there was no error.

Does that work for you as well?

# Answer
> 0 votes
I did it by copy and paste. Here's the top line of the output:

```
    Primer pair 1
Sequence (5'->3')   Template strand Length  Start   Stop    Tm  GC% Self complementarity    Self 3' complementarity
Forward primer  TGTCGTTGACAGGACACGAG    Plus    20  148 167 59.97   55.00   5.00    2.00
Reverse primer  TTACCTTTCGGTCACACCCG    Minus   20  264 245 59.97   55.00   3.00    2.00
Product length  117
Primer pair 2
Sequence (5'->3')   Template strand Length  Start   Stop    Tm  GC% Self complementarity    Self 3' complementarity
Forward primer  ATCTGTGTGGCTGTCACTCG    Plus    20  82  101 60.04   55.00   5.00    3.00
Reverse primer  ACTCGTGTCCTGTCAACGAC    Minus   20  168 149 59.97   55.00   4.00    3.00
Product length  87
Primer pair 3
Sequence (5'->3')   Template strand Length  Start   Stop    Tm  GC% Self complementarity    Self 3' complementarity
Forward primer  TACTGTCGTTGACAGGACACG   Plus    21  145 165 60.00   52.38   7.00    3.00
Reverse primer  ACACCCGGACGAAACCTAGA    Minus   20  251 232 60.54   55.00   4.00    2.00
Product length  107
Primer pair 4
Sequence (5'->3')   Template strand Length  Start   Stop    Tm  GC% Self complementarity    Self 3' complementarity
Forward primer  GCTTACGGTTTCGTCCGTGT    Plus    20  192 211 60.94   55.00   4.00    0.00
Reverse primer  GTTGAAACCAGGGACAAGGC    Minus   20  292 273 59.33   55.00   4.00    2.00
Product length  101

```

A command line solution would work, based on the bug reported by the OP. However there's nothing in NCBIs e-utilities and Biopython doesn't do Primer-Blast, but will do Blast. You could do it Perl via the LWP Agent, e.g. here.

---
Tags: blast, phylogenetics, primer
---


# Question
Title: Python module for fetching NCBI id for a list of species
I have a list of scientific names of species. Is there a python module that can fetch NCBI taxonomy IDs?

# Answer
Biopython has an efetch frontend.

```
from Bio import Entrez
Entrez.email = "mg@example.edu"
handle = Entrez.efetch(db="nucleotide", id="KR052012.1", rettype="gb", retmode="text")
print(handle.read())

```

---

**OUT**

LOCUS KR052012 10714 bp RNA linear VRL 09-MAY-2015 DEFINITION Dengue virus strain Hb33/CHN/2014, complete genome. ACCESSION KR052012 VERSION KR052012.1 KEYWORDS
.... etc

You will likely need to perform an EFetch search and then extract out the accession numbers for import (as done above)
> 1 votes

# Answer
> *This isn't a Python library, but it's easily invoked from Python code, easy to install, and has a lot of useful features.*

I just discovered the amazing taxonkit library this week. In particular, the `taxonkit name2taxid` command seems to be what you're looking for.

```
$ cat NAMES
Polistes dominula
Sinorhizobium meliloti
Gossypium hirsutum
$ cat NAMES | taxonkit name2taxid
Polistes dominula       743375
Sinorhizobium meliloti  382
Gossypium hirsutum      3635

```

Easily installed using `conda install taxonkit`.
> 1 votes

# Answer
If you want to perform a species search, here's the raw efetch code. I use this within an alias and the complex quotes are needed for complicated reasons. The reg-ex may need to be adapted to your search. $virus here is your species of interest, e.g. `set virus='some species'` (csh or tcsh), \`virus='some species' (bash). You'll just need to translate this to Biopython.

```
alias esummary-acc \
'esearch -db protein -query "$virus" |\
efetch -format docsum |\
xtract -pattern DocumentSummary -element Title,AccessionVersion |\
grep "$virus" |\
perl -p -e '"'"'s/^.*\[.*\]\s+(\w+)/$1/'"'"''

```
> 0 votes

---
Tags: python, ngs, genomics, metagenome, taxonomy
---


title: Download sequences of isoforms

text:
# Question
Title: Download sequences of isoforms
I want to collect all isoforms of all genes (Fasta/Fastq file of nucleotide and protein sequences) Arabidopsis Col-0. I am wondering if there is a straightforward way to download the file from any database.

I am new to sequence analysis, so I apologize if asking a very basic question.

# Answer
> 1 votes
Col-0 is, to my knowledge, the standard reference Arabidopsis thaliana ecotype, so the easiest way to get the sequences for all genes is to just download them from Ensembl. The nucleotide sequences are available here and the protein sequences here.

# Answer
> 0 votes
For RefSeq data, go to NCBI homepage and search for `Arabidopsis Col-0 genome`. Click on the 'Get Data' button in the result box and you will be directed to the NCBI Datasets portal with the latest RefSeq assembly selected. Click on the blue 'Download' button, uncheck the genome and GFF3 files from list, check transcript and protein FASTA files from the list and download data.

---
Tags: rna-seq, phylogenetics, sequence-analysis, phylogeny, isoform
---


# Question
Title: Can you give me a detailed explanation about PrimeBlast?
I am trying to write a research paper on PRIMER-BLAST. I have read something about it on

https://www.ncbi.nlm.nih.gov/pubmed/22708584/

but still, it is not enough. So can you please tell me more about it or suggest some resources.

# Answer
I would start at the Primer3 website here.

Once you understand all the variables on the website via the manual here and the associated publications Koressaar et al (2007) and Untergasser et al (2012) you will understand the gist of the algorithm used in primer design. This gives you optimised primers for a single strand of DNA/RNA

Primer-Blast looks in addition at maximising sequence conservation within an alignment of different sequences. The alignment is created using a Blast (I think its Mega-Blast. Thus Primer-Blast looks at one more parameter accounting for mutation heterogeneity which might be encountered in the PCR/RT-PCR and the number of primers it recommends is therefore reduced over say PRimer3.

Prove it for yourself take a given sequence and run it through Primer3 and then Primer-Blast and assess the differences and similarities.

*References*

Untergasser A, Cutcutache I, Koressaar T, Ye J, Faircloth BC, Remm M, Rozen SG (2012) Primer3 - new capabilities and interfaces. Nucleic Acids Research 40(15):e115

Koressaar T, Remm M (2007) Enhancements and modifications of primer design program Primer3 Bioinformatics 23(10):1289-91
> 0 votes

---
Tags: blast, primer
---


# Question
Title: After completing the evaluation, how to predict PTM lysine or any other modification with unlabeled data for making a web server?
I have created a model with SVM classifier and evaluated it with 5-fold cross-validation. The dataset was sequence data containing lysine modified sites. Now I want to build a web server where anyone can input any sequence in FASTA format. Apparently, I have no label information at this point. So how can I will be able to extract features or use any feature selection method which requires label information? How can I extract or create the test set to predict? Suppose I have used f-score feature selection method in my model and selected 50 features from the given dataset. But for web server, having no label information, how can I select 50 features from the sequence that are provided by an user and then test it with my previously trained model?

As being a novice in the world of bioinformatics, I am seeking the answer which might be improper for others. But it will be very much helpful for me. I am giving some link of some papers and its corresponding web servers.

1. iPTM-mLys: identifying multiple lysine PTM sites and their different type (https://academic.oup.com/bioinformatics/article/32/20/3116/2196525)

Web server link : http://www.jci-bioinfo.cn/iPTM-mLys
2. iHyd-PseCp: Identify hydroxyproline and hydroxylysine in proteins by incorporating sequence-coupled effects into general PseAAC (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5190098/)

Web server link : http://www.jci-bioinfo.cn/iHyd-PseCp

# Answer
> 1 votes
Interesting application. I would suggest using Python's Flask system, which is described in a few pages in O'Reilly's Programming PyTorch, chapter 8. PyTorch in Production. Django is not really recommended (don't say that on Stackoverflow however). O'Reilly do an entire book on Flask, but thats to make the site pretty.

Obviously, the explanation here is Deep Learning rather than ML but the basic gist will hold. You should also think about building a docker for distribution which is easy (same chapter). We're in lockdown so I don't have this book with me.

The **difficult part** will be web architecture because you shouldn't have the CPU/GPU on the same server as the quiry interface, otherwise you are prone to a denial of service attack or just sytem overload. To do this you need to consider Representational state transfer (RESTful) and that is not easy at all (at least for me). Flask facilitates RESTful, but thats like systems admin stuff. I guess before getting into RESTful I'd see how popular your website is, but if its the same server as your institute website you'd need to be really careful because an external user could bring your institutional website down (you wouldn't be popular).

---
Tags: fasta, sequence-analysis, proteomics
---


# Question
Title: Methods to predict one connected structure from two known separate structures
Currently, I am facing such a question: how to predict one connected protein structures from two separate known structures?

It should be different from the protein-protein docking method, which is non-covalent interactions. My question is that the two proteins should "contently" link together.

Any suggestions are appreciated.

# Answer
> 1 votes
There are two approaches here.

## Physiological

If domain A and B of your protein interact with each other physiologically, then you need to protein-protein docking. This is not too straightforward as protein are flexible.

Several online servers exists such as Haddock and include ClusPro, which use more traditional approaches. While the online server DEMO (I-Tasser group) does a structural homology based search and a final MD based repacking to better minimise the structures (no protein-protein docking).

Many of these work by magic as far as the user is concerned. But understanding what is going on is really important, even if it is overly laborious. When you don't know where a protein goes, it's called blind or global docking, which is computational expensive and so imprecise that the effort required to prove it is correct for peer-review is too steep for it to be useful. If you are doing local docking it is okay, which may be the case here depending on the linker. Specifically, you need to generate more 'decoys' (potential poses) the more unsure you are: tens of thousands or more. Then you have find clusters and see how each scores following more precise refinement. As an example of a tediously and manually set up run, Rosetta has several docking protocols, which explain the steps required —of note are the many judgement calls, such as in the local refinement, where one has to judge whether setting `-dock_pert 3 8` (3 Å and 8° perturbation) is too strict etc.

## Pretty picture

The protein-protein docking approach will result in a globular protein, which may be incorrect and majorly is not great for pictures.

If you want a pretty picture of two or more domains most likely you would want them placed like beads on a string. If you have a linker that is not solved the recommended distance to space them is $(3.5 N\_{missing})^{0.5}$ Å —unfortunately I do not know where this lore comes from. If it helps I wrote a script that does this in automatic but manual placement works best.

## Linkers

In the situation you have a good confidence about how the two parts fit together (or have placed them like beads on a string), you may want to add linkers. If you absolutely don't care about correct bond topology here is a blog post on how to join bits (shoddily) with PyMOL and here is one to do it properly with Rosetta or Pyrosetta.

If one wanted to do a better job, one tool to use is Rosetta remodel, which has a bit of a learning curve, but is very powerful for protein design. This requires the PDB to be on a single chain (pymol `alter` command, followed by `create` will do the trick: more here) and actually the amino acids to not have gaps —PDB number being the same as the pose numbering. Rosetta gives better results if the poses are energy minimised, but for this it does not matter too much. Then a "blueprint file" formatted as discussed in the documentation. Briefly, in a way that residues left alone are in the form `1 M .` (where residue 1 is a methionine), while the residues before and after the insertion are like `100 D L PIKAA D` (meaning D100 is changed to a Asp in loop SS; not actual numbers or residue), while between entries are added `0 X L PIKAA G` etc. To model loops, `-generic_aa G` is a good option as this tells the default rough fitting (centroid mode) to use glycine and not valine as generic residue. The top pose will actually not have a superb score, so a 3 (regular) or more cycles (15 thorough) with Rosetta Relax application either of the whole protein or with a movemap specific to the linker, may be done.

There actually is an option in Rosetta remodel to have a domain insertion, i.e. it would swing around like a medieval flail. However, it will not repack and therefore properly dock so this option is a bad idea for this specific application.

---
Tags: proteins, protein-structure, assembly
---


# Question
Title: Matching files by chromosome during wdl scatter gather
I'm trying to implement a polygenic score pipeline in WDL and I'm quite new to this pipeline management.

The first step of the pipeline takes a file of GWAS summary statistics and splits it by chromosome.

The corresponding wdl task:

```
task split {
    input {
        File gwas
        String output_prefix
    }

    command {
        ./splitGwas -i ${gwas} -o ${output_prefix}
    }

    output {
        Array [File] gwas_by_chr = glob("${output_prefix}_*.assoc")
    }
}

```

The next part of the process is to compute posterior effects of SNPs by chromosome. The inputs to this task are going to be the files split by chromosome in the previous step and an LD matrix for the chromosome.

The corresponding wdl call would be something like this:

```
Array [File] ld_matrices = read_lines(file_of_ld_matrices_by_chr)

scatter(pair in zip(split.gwas_by_chr, ld_matrices)) {
    call sbayes.run {
        input:
            GCTB = gctb_executable_path,
            gwas = pair.left,
            ld_matrix = pair.right,
            output_prefix = out
    }
}

```

Since I'm grabbing the split files using glob in the previous step, I do not particularly know which element in the array corresponds to which chromosome. So I would be making a mistake if I zip the summary stats for one chromosome with the LD matrix corresponding to a different chromosome.

The programmatic way would be to iterate the scatter over chromosome name but then I would lose the ability to specify the dependency between the two tasks.

Is there a better way to do this?

Cross posted on biostars (https://www.biostars.org/p/434582/) so pardon me if that is discouraged.

# Answer
> 1 votes
You could use the `sub` function to obtain the chromosome from the globbed file names and then use that to retrieve the correct LD matrices from a Map (see `read_map` or `read_json`).

eg.:

```
Map[String, String] ld_matrices = read_json(file_of_ld_matrices_by_chr)

scatter(gwas_for_chr in split.gwas_by_chr) {
    String chr = sub(sub(gwas_for_chr, prefix, ""), "\.assoc$", "")

    call sbayes.run {
        input:
            GCTB = gctb_executable_path,
            gwas = gwas_for_chr,
            ld_matrix = ld_matrices[chr],
            output_prefix = out
    }
}

```

The `file_of_ld_matrices_by_chr` would have to look like this:

```
{
  "chr1": "/path/to/file1",
  "chr2": "/path/to/file2"
}

```

---
Tags: wdl
---


# Question
Title: Model the effect of treatment on mutation burden
I am struggling a bit to model the following problem. Basically, I would like to model tumor mutation burden (mutations per megabase, a continuous) as a function of treatment (categorical). The following is an example of how data looks like

```
patient-id therapy-1 therapy-2 therapy-3 cancer-type TMB
PID-1      1         NA        1         CT-1        1.3
PID-2      0         0         1         CT-1        0.14
PID-3      NA        0         1         CT-2        5
PID-4      1.        1         NA        CT-3        2.3
PID-5      1         1         NA        CT-3        10
...
...

```

1 - treated; 0 - not treated; NA - No information

I was trying multiple linear regression, but in most cases one patient got multiple therapies. This is where I am confused on how to interpret the results of linear regression. So far, I haven't included the `cancer-type` information into the model. However, IMO, including cancer-type into model will yield specific treatment effects in specific cancer type? Or is it over-expectation of what we can achieve from linear model?

Also, is there a better way to address the problem to get reasonable associations between treatment and TMB in this case?

Thank you for the help. Stay safe!

# Answer
> 0 votes
Depends how you consider "linear regression", if you mean Pearsons, Spearman's analysis then general linear modelling (GLM) is a far better approach if you have continuous data. I don't really understand TMB, but that is a critical output then GLM will work. Logistic regression is an alternative approach data type depending. You have to perform a Q-Q plot and the residual against against any transformations you make (you might need to transform TMB). The lower the residual the better, the straight the Q-Q plot the better.

To explain how GLM works long hair and height correlate in a straight regression analysis. The implication is if you want to be tall - cut your hair! Of once you account for gender and send it through a GLM then gender and hairlength will correlate far stronger than hair length and height.

Mac OSX has a very easy to use pointy clicky, but highly professional, GLM package on its AppStore called Wizzard (bit expensive). It doesn't do general linear mixed modelling however. Everything does GLM of course, its just this package is really easy to use and you can figure out GLM by trial and error.

---
Tags: r, python, mutations, linear-regression
---


# Question
Title: How can I bypass 'out of memory' problem using a HISAT2 for human genome indexing?
I am trying to perform a `HISAT2` indexing using a `GCF_000001405.39_GRCh38.p13_genomic.gtf` and `GCF_000001405.39_GRCh38.p13_genomic.fna` files downloaded from https://www.ncbi.nlm.nih.gov/assembly/GCF\_000001405.39 using the procedure described in Pertea et al. (https://www.ncbi.nlm.nih.gov/pubmed/27560171), e.g. first I create `extract_splice_sites.py` and `extract_exons.py` and then my code looks like this

```
#!/bin/bash
#SBATCH --time=02:30:00
#SBATCH --account=def-myusername
#SBATCH --mem=160000M
#SBATCH --mail-user=mymail@gmail.com
#SBATCH --mail-type=All

module load nixpkgs/16.09  intel/2018.3
module load hisat2/2.1.0

hisat2-build --ss imam.ss --exon imam.exon GCF_000001405.39_GRCh38.p13_genomic.fna imam002

```

I am trying to receive a set of `imam002*.ht2` files but every time the procedure ends with a `Ran out of memory` message. From here https://github.com/griffithlab/rnaseq\_tutorial/wiki/Indexing I learned that I must allocate 160Gb of RAM if I want to index the entire genome. Even when I do that and submit my script via sbatch on a HPC I still got this

```
slurmstepd: error: Detected 1 oom-kill event(s) in step 30827017.batch cgroup.

```

Just in case here's the window that I have when I don't sbatch my script, I still have an error

```
Inactive Modules:
  1) openmpi/2.1.1

The following have been reloaded with a version change:
  1) gcccore/.5.4.0 => gcccore/.7.3.0     2) icc/.2016.4.258 => icc/.2018.3.222     3) ifort/.2016.4.258 => ifort/.2018.3.222     4) imkl/11.3.4.258 => imkl/2018.3.222     5) intel/2016.4 => intel/2018.3

Settings:
  Output files: "imam002.*.ht2"
  Line rate: 7 (line is 128 bytes)
  Lines per side: 1 (side is 128 bytes)
  Offset rate: 4 (one in 16)
  FTable chars: 10
  Strings: unpacked
  Local offset rate: 3 (one in 8)
  Local fTable chars: 6
  Local sequence length: 57344
  Local sequence overlap between two consecutive indexes: 1024
  Endianness: little
  Actual local endianness: little
  Sanity checking: disabled
  Assertions: disabled
  Random seed: 0
  Sizeofs: void*:8, int:4, long:8, size_t:8
Input files DNA, FASTA:
  GCF_000001405.39_GRCh38.p13_genomic.fna
Reading reference sizes
  Time reading reference sizes: 00:00:23
Calculating joined length
Writing header
Reserving space for joined string
Joining reference sequences
  Time to join reference sequences: 00:00:13
  Time to read SNPs and splice sites: 00:00:15
  Ran out of memory; automatically trying more memory-economical parameters.
./hisat2_index.sh: line 24:   943 Killed                  hisat2-build --ss imam.ss --exon imam.exon GCF_000001405.39_GRCh38.p13_genomic.fna imam002

```

When I use the same script for another genome (not human) - I don't have this problem.

Is it really a problem of RAM and, if it is, how is usually human genome gets indexed by HISAT2? Is it even possible?

# Answer
> 1 votes
You haven't allocated the requisite 160GB of RAM required. Or rather, you allocated 160GB of memory, which is what you were told to do, but they most likely meant 160**GiB**.

Computers use binary, not decimal, so a gigabyte is 1024 megabytes. This is the same annoying mix-up people run into when their hard drives don't have the space the marketing department says it does. By requesting 160000M in your job script, you actually requested a few gigabytes less than what your source says you need.

$$ 160000 \text{ megabytes} \times \frac{1 \text{ gigabyte}}{1024 \text{ megabytes}} = 156.25 \text{ gigabytes} $$

$$ 160 \text{ gigabytes} \times \frac{1024 \text{ megabytes}}{1 \text{ gigabyte}} = 163840 \text{ megabytes} $$

If "160GB" is what you need, then `163840M` is the correct value for the job script.

I would honestly suggest cutting the Gordian knot and not dealing with this "gigabytes vs. gibibytes" ordeal by simply using the `G` suffix.

```
#SBATCH --mem=160G

```

---
Tags: rna-seq, human-genome, genome-index
---


# Question
Title: Running snakemake with multiple conda environments without rebuilding environments every time
I have some snakemake pipelines that have rules which run in different environments. Is it possible to run the pipeline without having to rebuild the environments from scratch every time when run in a different directory?

This is the way I would normally run snakemake:

`snakemake --directory <my_output_dir> --use-conda`

My snakefile might have rules that need to use different environments:

```
rule A:
   ....
   conda: envA.yml
   ....

rule B:
   ....
   conda: envB.yml
   ....

```

So if I wanted to do two different analyses with the same pipeline, is there anyway that I can point them to the same environments without having to install everything separately in each output directory?

The snakemake documentation makes it seem like they don't currently support this. You can either use 1) --conda-prefix when running snakemake and then point to a single env or 2) use the --use-conda directive and the env.yml files by rule, in which case it will create new envs any time you run it in a new output directory. So my case doesn't seem to be covered, but does anyone have a workaround to suggest? I often find conflicts between some of the R packages I use and python packages, so I tend to keep them in their own environments to avoid headache.

I'm starting to get neck-deep in conda environments at this point which is taking up space from my disk quota, but also I am bothered by the fact that I may actually want to run analysis of data from different projects with the same environments. If the environments are rebuilt at each execution, I do not know that they are the same, short of specifying the exact version info for everything I use in the env.yaml files (and even then, they could be slightly different if the repo has changed).

I have looked at this similar question: Running Snakemake in one single conda env in which it is suggested there is a way to do such a thing, but I'm afraid I don't understand how to implement it without an example.

Thanks

# Answer
> 3 votes
Each YAML file will always create its own environment. The packages in each aren't downloaded and installed each time, if they're already installed elsewhere they're just hard-linked in. You'll want to use both `--use-conda` and `--conda-prefix`. The latter will point to the base directory in which all conda environments will be stored. This has the benefit that running the same workflow in different places won't result in a new env being created (by default snakemake will do so inside a .snakemake directory otherwise).

# Answer
> 1 votes
@vantom, while I am not very familiar with snakemake on a whole, I do have extensive knowledge writing workflows in WDL as well as CWL. The approach both of these languages take to overcome a similar issue that you are describing, is to use containerization through ie docker. This approach allows you to bundle the specific tools you need for a specific task, making different images for different tasks. There are several large repositories of images online ie: Docker Hub, or more biologically focused: Biocontainers

Looking through the documentation, it does appear that snakemake supports running your rule from within the context of a docker image. Without knowing how your compute environment is setup, I am not sure if there is any additional configuration you would need, however the documentation seems to suggest that it should be as easy as describing the name of the container. As long as you can run docker, (or an alternative like singularity) you should be able to avoid the current headache you are facing.

```
rule NAME:
    input:
        "table.txt"
    output:
        "plots/myplot.pdf"
    containers:
        "docker://joseespinosa/docker-r-ggplot2"
    script:
        "scripts/plot-stuff.R"

```

---
Tags: snakemake
---


# Question
Title: "megahit: command not found" after correct installation with miniconda
I want to use `MEGAHIT` (git: https://github.com/voutcn/megahit) to assemble metagenome shotgun sequencing reads (illumina).

I install the package with the command `conda install -c bioconda megahit` as described in on the github page which yields:

```
Solving environment: done

## Package Plan ##

  environment location: /home/datura-workstation/miniconda3

  added / updated specs:
    - megahit

The following NEW packages will be INSTALLED:

  megahit            bioconda/linux-64::megahit-1.2.9-h8b12597_0

Proceed ([y]/n)?

```

Where I click yes and the packages is installed. However, after installation I am not able to use the command `megahit` as this results in the error: `megahit: command not found`.

I use Ubuntu 18.04.3 LTS and have installed `miniconda3`.

I reproduced the installation of `megahit` on a windows computer with a virtual machine using the same Ubuntu version. On my virtual Ubuntu I CAN use the `megahit` command.

Could anyone please help me solve this problem? I searched the internet but cannot find a good solution on the issue (also not on the github issue page of this package).

# Answer
> 2 votes
type `conda info --envs`. Look at your list of environments. Say it's `base` that has `megahit` in it, do `conda activate base`, then just do `megahit`.

---
Tags: assembly
---


# Question
Title: How to run Cox proportional Hazard model for each group in R?
My population is divided into 10 groups. I want to run the cox-proportional hazard model separately for each of these 10 groups so that in the end i can compare which groups has highest hazard ratio depending uopn the covariates. How can i do this in R? Earlier i was using this formula.

```
Model_1 <- coxph(serve_object ~ GROUP + age+ Smoking+ Diabetes, data = DF
summary(Model_1)

```

But the output is very confusing. it gives hazard ratio of 9 out of 10 groups. I think it dealt group as covariate.

# Answer
> 2 votes
The fact that you supplied `GROUP` as a covariate, you are getting 9 coefficients as one of the ten levels is used as the reference level.

If you would like to fit one model per group (and in this scenario you would not be correcting for `GROUP`), you should first split your data frame and fit your models on the resulting list's elements:

```
lapply(split(DF, DF$GROUP),
       function(x) coxph(serve_object ~ age + Smoking + Diabetes, data = DF))

```

Example with the `mtcars` data set when grouped by number of cylinders, `cyl` having values of 4, 6 or 8:

```
> lapply(split(mtcars, mtcars$cyl),
+        function(x) lm(disp ~ hp, data = x))
$`4`

Call:
lm(formula = disp ~ hp, data = x)

Coefficients:
(Intercept)           hp
    59.0369       0.5579

$`6`

Call:
lm(formula = disp ~ hp, data = x)

Coefficients:
(Intercept)           hp
   290.9178      -0.8799

$`8`

Call:
lm(formula = disp ~ hp, data = x)

Coefficients:
(Intercept)           hp
   320.2083       0.1572

```

---
Tags: statistics, statsmodels
---


# Question
Title: Is it possible to go from a digital sequence to an actual nucleotide?
I'm trying to learn about what the state of technology is at presently. It seems like we clearly can go from nucleotide to digitally stored sequence, but can we transcribe something from the digital sequence?

# Answer
> 3 votes
Oligonucleotide synthesis is a service offered by many providers such as IDT.

# Answer
> 0 votes
To add to @Maximilian Press' answer.

There is quite a lot of research on how to store arbitrary digital data in DNA. Nick Goldman describes his research on the topic in his TEDx talk. The last progress in the field I have noticed was DNA fountain, a tranformation of bites in DNA with maximal efficiency of data storage in DNA (information / nt) and reliable recovery when sequenced afterwards.

---
Tags: fasta, sequence-analysis, phylogenetics
---


# Question
Title: Intersecting a data frame based on the multiple columns
I have 2 data frames

```
> head(df1)
   Chromosome Start_Position End_Position       Gene
1:      chr12       46254706     46254706       ARID2
2:      chr17        7579849      7579849        TP53
3:       chr9       21974777     21974780      CDKN2A
4:       chr1       27105930     27105930      ARID1A
5:       chr1      201982973    201982973        ELF3
6:       chr3       89480459     89480459       EPHA3
> head(df2)
  Chromosome Start_Position End_Position
1          1          10583       863511
2          1       12841835     12854479
3          1       17231299     17232877
4          1      120380739    142540462
5          1      142540484    142684781
6          1      142685823    143462598
>

```

How I can intersect genes from df1 which are in the range of df2? Because both data frames have chromosome number, start and end but I don't know how to identify genes from df1 based on df2

I have tried

```
df1[df1$Start_Position%in%df2$Start_Position,]

```

But I obtain one gene

I also tried

```
> new_dataset <- df1 %>% right_join(df2, by=c("Chromosome","Start_Position"))
Error: Can't join on 'Start_Position' x 'Start_Position' because of incompatible types (character / integer)
>

```

# Answer
> 1 votes
reformat df2 so that you have `chr1` is instead of `1`, then you basically have two bed files. Then use `bedtools intersect`. https://bedtools.readthedocs.io/en/latest/content/tools/intersect.html

# Answer
> 0 votes
In R:

```
> write.table(df1, 'df1.txt', quote=F, sep="\t", row.names=F, col.names=F)
> write.table(df2, 'df2.txt', quote=F, sep="\t", row.names=F, col.names=F)

```

Convert both files to sorted BED:

```
$ sort-bed df1.txt > df1.bed
$ awk -v FS="\t" -v OFS="\t" '{ print "chr"$1 }' df2.txt | sort-bed - > df2.bed

```

Then map `df1` gene names to `df2` coordinates:

```
$ bedmap --echo --echo-map-id-uniq --delim '\t' df2.bed df1.bed > answer.bed

```

To bring back into R, e.g.:

```
> dfMap <- read.table('answer.bed', sep='\t', col.names=c('Chromosome', 'Start_Position', 'End_Position', 'Genes'), colClasses=c("factor", "numeric", "numeric", "character"))

```

---
Tags: r
---


# Question
Title: Protein PDB/CIF text file reader that can fold a same section into one line
Is there a software that can recognise the protein PDB/CIF files, and fold multiple lines of a same section into one line, so that it is easier to handle? I need to understand the format in CIF so as to create my own CIF file.

For example, when `Notepad++` is used to read a `.py` file, it can fold the same section into one line.

So for a PDB/CIF file, e.g. all the atom lines for a same chain can be folded into one line, and I can expand it if needed.

Like this picture, assuming "This PC" is the whole protein, while the other objects are the different chains.

# Answer
Programs that work with PDB and mmCIF files may present the structure as a tree, for example:

(this is a Find dialog from ccp4MG)

but I'd be surprised if a text editor had such a feature. Normally, you don't edit PDB files in a text editor.
> 1 votes

---
Tags: pdb, text-processing
---


# Question
Title: Can't install snakemake v5.15.0 from conda
I'm trying to install the latest version of snakemake from conda into a conda env, following the instructions here. I.e.

`conda create -n snakemake -c conda-forge -c bioconda snakemake`

However, this seems to only give me the option to install snakemake v5.3.0, while the latest `no-arch` version on conda is v5.15.0.

Any help greatly appreciated.

```
Solving environment: done

## Package Plan ##

  environment location: /home/ubuntu/anaconda3/envs/snakemake

  added / updated specs:
    - snakemake

The following packages will be downloaded:

    package                    |            build
    ---------------------------|-----------------
    aioeasywebdav-2.4.0        |        py36_1000          18 KB  conda-forge
    aiohttp-3.6.2              |   py36h7b6447c_0         544 KB
    appdirs-1.4.3              |   py36h28b3542_0          15 KB
    async-timeout-3.0.1        |           py36_0          12 KB
    attrs-19.3.0               |             py_0          39 KB
    bcrypt-3.1.7               |   py36h7b6447c_0          40 KB
    blas-2.14                  |         openblas          10 KB  conda-forge
    boto3-1.9.66               |           py36_0         105 KB
    botocore-1.12.189          |             py_0         3.4 MB
    datrie-0.8                 |   py36h7b6447c_0         147 KB
    docutils-0.16              |           py36_0         668 KB
    dropbox-9.4.0              |           py36_0         756 KB
    filechunkio-1.6            |           py36_0           7 KB  bioconda
    ftputil-3.2                |           py36_0          85 KB  bioconda
    google-api-core-1.16.0     |           py36_1          86 KB
    googleapis-common-protos-1.51.0|           py36_2          72 KB
    idna_ssl-1.1.0             |           py36_0           7 KB
    importlib_metadata-1.5.0   |           py36_0          48 KB
    jsonschema-3.2.0           |           py36_0          95 KB
    liblapacke-3.8.0           |      14_openblas          10 KB  conda-forge
    multidict-4.7.3            |   py36h7b6447c_0          68 KB
    networkx-2.4               |             py_0         1.2 MB
    numpy-1.17.0               |   py36h99e49ec_0          24 KB  r
    numpy-base-1.17.0          |   py36h2f8d375_0         5.2 MB  r
    pandas-1.0.3               |   py36h0573a6f_0         8.6 MB
    protobuf-3.11.4            |   py36he6710b0_0         635 KB
    psutil-5.7.0               |   py36h7b6447c_0         315 KB
    pygraphviz-1.3.1           |           py36_0         205 KB  bioconda
    pynacl-1.3.0               |   py36h7b6447c_0         1.1 MB
    pyrsistent-0.16.0          |   py36h7b6447c_0          93 KB
    pysftp-0.2.9               |           py36_0          31 KB  bioconda
    ratelimiter-1.2.0          |        py36_1000          12 KB  conda-forge
    rsa-3.1.4                  |           py36_0          87 KB  bioconda
    s3transfer-0.1.13          |           py36_0          79 KB
    snakemake-5.3.0            |           py36_1           4 KB  bioconda
    snakemake-minimal-5.3.0    |           py36_1         280 KB  bioconda
    sqlite-3.31.1              |       h62c20be_1         2.0 MB
    typing_extensions-3.7.4.1  |           py36_0          40 KB
    wrapt-1.12.1               |   py36h7b6447c_1          49 KB
    yarl-1.4.2                 |   py36h7b6447c_0         132 KB
    zipp-2.2.0                 |             py_0          12 KB
    ------------------------------------------------------------
                                           Total:        26.2 MB

The following NEW packages will be INSTALLED:

  _libgcc_mutex      pkgs/main/linux-64::_libgcc_mutex-0.1-main
  aioeasywebdav      conda-forge/linux-64::aioeasywebdav-2.4.0-py36_1000
  aiohttp            pkgs/main/linux-64::aiohttp-3.6.2-py36h7b6447c_0
  appdirs            pkgs/main/linux-64::appdirs-1.4.3-py36h28b3542_0
  asn1crypto         pkgs/main/linux-64::asn1crypto-1.3.0-py36_0
  async-timeout      pkgs/main/linux-64::async-timeout-3.0.1-py36_0
  attrs              pkgs/main/noarch::attrs-19.3.0-py_0
  bcrypt             pkgs/main/linux-64::bcrypt-3.1.7-py36h7b6447c_0
  blas               conda-forge/linux-64::blas-2.14-openblas
  boto3              pkgs/main/linux-64::boto3-1.9.66-py36_0
  botocore           pkgs/main/noarch::botocore-1.12.189-py_0
  ca-certificates    pkgs/main/linux-64::ca-certificates-2020.1.1-0
  cachetools         pkgs/main/noarch::cachetools-3.1.1-py_0
  cairo              pkgs/main/linux-64::cairo-1.14.12-h8948797_3
  certifi            pkgs/main/linux-64::certifi-2020.4.5.1-py36_0
  cffi               pkgs/main/linux-64::cffi-1.14.0-py36h2e261b9_0
  chardet            pkgs/main/linux-64::chardet-3.0.4-py36_1003
  configargparse     pkgs/main/noarch::configargparse-1.1-py_0
  cryptography       pkgs/main/linux-64::cryptography-2.8-py36h1ba5d50_0
  datrie             pkgs/main/linux-64::datrie-0.8-py36h7b6447c_0
  decorator          pkgs/main/noarch::decorator-4.4.2-py_0
  docutils           pkgs/main/linux-64::docutils-0.16-py36_0
  dropbox            pkgs/main/linux-64::dropbox-9.4.0-py36_0
  expat              pkgs/main/linux-64::expat-2.2.6-he6710b0_0
  filechunkio        bioconda/linux-64::filechunkio-1.6-py36_0
  fontconfig         pkgs/main/linux-64::fontconfig-2.13.0-h9420a91_0
  freetype           pkgs/main/linux-64::freetype-2.9.1-h8a8886c_1
  fribidi            pkgs/main/linux-64::fribidi-1.0.5-h7b6447c_0
  ftputil            bioconda/linux-64::ftputil-3.2-py36_0
  gitdb              pkgs/main/noarch::gitdb-4.0.2-py_0
  gitpython          pkgs/main/noarch::gitpython-3.1.1-py_1
  glib               pkgs/main/linux-64::glib-2.63.1-h5a9c865_0
  google-api-core    pkgs/main/linux-64::google-api-core-1.16.0-py36_1
  google-auth        pkgs/main/noarch::google-auth-1.13.1-py_0
  google-cloud-core  pkgs/main/noarch::google-cloud-core-1.3.0-py_0
  google-cloud-stor~ pkgs/main/noarch::google-cloud-storage-1.27.0-py_0
  google-resumable-~ pkgs/main/noarch::google-resumable-media-0.5.0-py_1
  googleapis-common~ pkgs/main/linux-64::googleapis-common-protos-1.51.0-py36_2
  graphite2          pkgs/main/linux-64::graphite2-1.3.13-h23475e2_0
  graphviz           pkgs/main/linux-64::graphviz-2.40.1-h21bd128_2
  harfbuzz           pkgs/main/linux-64::harfbuzz-1.8.8-hffaf4a1_0
  icu                pkgs/main/linux-64::icu-58.2-h9c2bf20_1
  idna               pkgs/main/noarch::idna-2.9-py_1
  idna_ssl           pkgs/main/linux-64::idna_ssl-1.1.0-py36_0
  importlib_metadata pkgs/main/linux-64::importlib_metadata-1.5.0-py36_0
  jinja2             pkgs/main/noarch::jinja2-2.11.1-py_0
  jmespath           pkgs/main/noarch::jmespath-0.9.4-py_0
  jpeg               pkgs/main/linux-64::jpeg-9b-h024ee3a_2
  jsonschema         pkgs/main/linux-64::jsonschema-3.2.0-py36_0
  ld_impl_linux-64   pkgs/main/linux-64::ld_impl_linux-64-2.33.1-h53a641e_7
  libblas            conda-forge/linux-64::libblas-3.8.0-14_openblas
  libcblas           conda-forge/linux-64::libcblas-3.8.0-14_openblas
  libedit            pkgs/main/linux-64::libedit-3.1.20181209-hc058e9b_0
  libffi             pkgs/main/linux-64::libffi-3.2.1-hd88cf55_4
  libgcc-ng          pkgs/main/linux-64::libgcc-ng-9.1.0-hdf63c60_0
  libgfortran-ng     pkgs/main/linux-64::libgfortran-ng-7.3.0-hdf63c60_0
  liblapack          conda-forge/linux-64::liblapack-3.8.0-14_openblas
  liblapacke         conda-forge/linux-64::liblapacke-3.8.0-14_openblas
  libopenblas        conda-forge/linux-64::libopenblas-0.3.7-h5ec1e0e_6
  libpng             pkgs/main/linux-64::libpng-1.6.37-hbc83047_0
  libprotobuf        pkgs/main/linux-64::libprotobuf-3.11.4-hd408876_0
  libsodium          pkgs/main/linux-64::libsodium-1.0.16-h1bed415_0
  libstdcxx-ng       pkgs/main/linux-64::libstdcxx-ng-9.1.0-hdf63c60_0
  libtiff            pkgs/main/linux-64::libtiff-4.1.0-h2733197_0
  libuuid            pkgs/main/linux-64::libuuid-1.0.3-h1bed415_2
  libxcb             pkgs/main/linux-64::libxcb-1.13-h1bed415_1
  libxml2            pkgs/main/linux-64::libxml2-2.9.9-hea5a465_1
  markupsafe         pkgs/main/linux-64::markupsafe-1.1.1-py36h7b6447c_0
  multidict          pkgs/main/linux-64::multidict-4.7.3-py36h7b6447c_0
  ncurses            pkgs/main/linux-64::ncurses-6.2-he6710b0_0
  networkx           pkgs/main/noarch::networkx-2.4-py_0
  numpy              r/linux-64::numpy-1.17.0-py36h99e49ec_0
  numpy-base         r/linux-64::numpy-base-1.17.0-py36h2f8d375_0
  openssl            pkgs/main/linux-64::openssl-1.1.1g-h7b6447c_0
  pandas             pkgs/main/linux-64::pandas-1.0.3-py36h0573a6f_0
  pango              pkgs/main/linux-64::pango-1.42.4-h049681c_0
  paramiko           pkgs/main/noarch::paramiko-2.7.1-py_0
  pcre               pkgs/main/linux-64::pcre-8.43-he6710b0_0
  pip                pkgs/main/linux-64::pip-20.0.2-py36_1
  pixman             pkgs/main/linux-64::pixman-0.38.0-h7b6447c_0
  prettytable        conda-forge/noarch::prettytable-0.7.2-py_3
  protobuf           pkgs/main/linux-64::protobuf-3.11.4-py36he6710b0_0
  psutil             pkgs/main/linux-64::psutil-5.7.0-py36h7b6447c_0
  pyasn1             pkgs/main/noarch::pyasn1-0.4.8-py_0
  pyasn1-modules     pkgs/main/noarch::pyasn1-modules-0.2.7-py_0
  pycparser          pkgs/main/noarch::pycparser-2.20-py_0
  pygraphviz         bioconda/linux-64::pygraphviz-1.3.1-py36_0
  pynacl             pkgs/main/linux-64::pynacl-1.3.0-py36h7b6447c_0
  pyopenssl          pkgs/main/linux-64::pyopenssl-19.1.0-py36_0
  pyrsistent         pkgs/main/linux-64::pyrsistent-0.16.0-py36h7b6447c_0
  pysftp             bioconda/linux-64::pysftp-0.2.9-py36_0
  pysocks            pkgs/main/linux-64::pysocks-1.7.1-py36_0
  python             pkgs/main/linux-64::python-3.6.10-hcf32534_1
  python-dateutil    pkgs/main/noarch::python-dateutil-2.8.1-py_0
  python-irodsclient conda-forge/noarch::python-irodsclient-0.8.2-py_0
  pytz               pkgs/main/noarch::pytz-2019.3-py_0
  pyyaml             pkgs/main/linux-64::pyyaml-5.3.1-py36h7b6447c_0
  ratelimiter        conda-forge/linux-64::ratelimiter-1.2.0-py36_1000
  readline           pkgs/main/linux-64::readline-8.0-h7b6447c_0
  requests           pkgs/main/linux-64::requests-2.23.0-py36_0
  rsa                bioconda/linux-64::rsa-3.1.4-py36_0
  s3transfer         pkgs/main/linux-64::s3transfer-0.1.13-py36_0
  setuptools         pkgs/main/linux-64::setuptools-46.1.3-py36_0
  six                pkgs/main/linux-64::six-1.14.0-py36_0
  smmap              pkgs/main/noarch::smmap-3.0.2-py_0
  snakemake          bioconda/linux-64::snakemake-5.3.0-py36_1
  snakemake-minimal  bioconda/linux-64::snakemake-minimal-5.3.0-py36_1
  sqlite             pkgs/main/linux-64::sqlite-3.31.1-h62c20be_1
  tk                 pkgs/main/linux-64::tk-8.6.8-hbc83047_0
  typing_extensions  pkgs/main/linux-64::typing_extensions-3.7.4.1-py36_0
  urllib3            pkgs/main/linux-64::urllib3-1.25.8-py36_0
  wheel              pkgs/main/linux-64::wheel-0.34.2-py36_0
  wrapt              pkgs/main/linux-64::wrapt-1.12.1-py36h7b6447c_1
  xmlrunner          conda-forge/noarch::xmlrunner-1.7.7-py_0
  xz                 pkgs/main/linux-64::xz-5.2.5-h7b6447c_0
  yaml               pkgs/main/linux-64::yaml-0.1.7-had09818_2
  yarl               pkgs/main/linux-64::yarl-1.4.2-py36h7b6447c_0
  zipp               pkgs/main/noarch::zipp-2.2.0-py_0
  zlib               pkgs/main/linux-64::zlib-1.2.11-h7b6447c_3
  zstd               pkgs/main/linux-64::zstd-1.3.7-h0b5b093_0

Proceed ([y]/n)?

```

# Answer
This is an example of the conda solver behaving incorrectly. In this case, it's choosing to pick a version that minimizes the number of dependencies required. Version 5.3.0 has 25 fewer dependencies, so if you let conda choose the version that's what it will pick. There's nothing you can really do about this other than specifying the version of snakemake that you want.

Update: It's very likely that mamba will correctly do what you want by default. Though I haven't used it, it's supposed to be a drop in replacement for conda here.
> 4 votes

# Answer
I don't know if it means anything but...

* It seems you are using anaconda instead of miniconda as required by the snakemake docs. Maybe reinstall conda using miniconda3.
* conda is installing python 3.6 instead of the more recent 3.8... odd
> 0 votes

---
Tags: snakemake, conda
---


# Question
Title: remove sequences from fasta file matching a string in the header
I have a file with 16S sequences. some headers contain species information. For my purposes I would like to exclude a number of species from the file, therefore I would like to do a pattern matching within a 16S header (that contains other characters beyond the species name) to then remove that header and the following sequence.

So, I have:

```
>NCBI_ID1;high_rank_taxonomy;...;species_name1
ATCCTT
>NCBI_ID2;high_rank_taxonomy;...;species_name2
ATTTTCCC
>NCBI_ID3;high_rank_taxonomy;...;species_name3
TTTCCCCAAAA
>NCBI_ID4;high_rank_taxonomy;...;species_name4
CCCTTTAAA

```

I would like to exclude, say, species 1,2 and 4 to retain:

```
>NCBI_ID3;high_rank_taxonomy;...;species_name3
TTTCCCCAAAA

```

I found this solution in an adjacent question:

```
awk 'BEGIN{while((getline<"ids.txt")>0)l[">"$1]=1}/^>/{f=!l[$1]}f' seq.fa

```

where I have my species names in the file ids.txt and sequences in file seq.fa

but that would work if the whole of my header was just the species name (contents of each line in ids.txt), which is not the case. I thought I could be smart and use `*species_name1*` instead of `species_name1` as contents of ids.txt but that apparently doesn't work. Could someone help me with modifying the current solution to accommodate the existence of additinal strings in the headers?

Many thanks! A

# Answer
> 2 votes
you may use simple R commands, using the 'seqinr' package

```
#load package
library(seqinr)
#load file containing sequences
data<-read.fasta("filename.fa")
#create a vector containing species names: these are the last part of the string
vec.names<-unlist(lapply(strsplit(names(data), ";"), function(x)x[length(x)]))
#find names to keep: indices which are not in the species to remove
species.to.remove<-c("species_name1", "species_name2", "species_name4")
vec.tokeep<-which(! vec.names %in%  species.to.remove)
#write the final output
write.fasta(sequences=data[vec.tokeep], names=names(data)[vec.tokeep], file.out="output.fa")

```

---
Tags: fasta, bash, shell, awk
---


# Question
Title: Gene information from chromosome location
I have a file that contains species name, chromosome, start and end locations like the following format `hg19,chr8,77778733,77779026`. I would like to know some information about the regions `77778733-77779026`. How can I quickly find out if it lies in an intron or exon of a gene? Or it is in the upstream or downstream regions of a gene? Is there a tool or web service from where I can get this information?

# Answer
> *How can I quickly find out if it lies in an intron or exon of a gene?*

I have an answer on Biostars that might help: https://www.biostars.org/p/124515/#124522

> *Or it is in the upstream or downstream regions of a gene?*

One option is to build a table of annotated bins and then map your regions-of-interest to those annotations.

These annotations would consist of whether the bin is inside a gene, or, after picking the nearest gene by strand, locate the bin to an upstream region of the gene TSS or to an downstream region of the gene TES.

To start, one could download stranded gene annotations:

```
$ GENCODE_URL=ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.basic.annotation.gff3.gz
$ GENCODE_GENES=${PWD}/../data/gencode.v28.genes.bed
$ wget -qO- "${GENCODE_URL}" | gunzip -c | awk '($3=="gene")' | gff2bed - | awk -v FS="\t" -v OFS="\t" '{ n=split($10,a,";"); p=a[4]; n=split(p,a,"="); $4=a[2]; print $0; }' | cut -f1-6 > ${GENCODE_GENES}

```

Chop up nuclear chromosomes of the genome of interest into 1kb bins (or whatever gradation of bin size you prefer):

```
$ ASSEMBLY=hg38
$ GENOME=${PWD}/../data/${ASSEMBLY}.bed
$ CHOP_SIZE=1000
$ CHOPPED_GENOME=${PWD}/../data/chopped.${ASSEMBLY}.bed
$ fetchChromSizes ${ASSEMBLY} | awk -v FS="\t" -v OFS="\t" '{ print $1, "0", $2 }' | grep -ve "\.*_\.*" | grep -v "chrM" | sort-bed - > ${GENOME}
$ bedops --chop ${CHOP_SIZE} ${GENOME} > ${CHOPPED_GENOME}

```

Then assign the midpoints of bins inside a gene, or to the nearest of any particular set of intergenic categories.

Here's a script I wrote which uses the output of `closest-features` to do that. This could be easily modified for any set of distance categories:

```
#!/usr/bin/env python

import sys
import numpy as np

'''
Input should be bin midpoints run through closest-elements, against some set of BED6+ (stranded) gene annotations
'''

def assign_to_category(dist):
  if dist < -50000:
    return '-50kb_to_-Inf'
  elif dist < -20000 and dist >= -50000:
    return '-20kb_to_-50kb'
  elif dist < -10000 and dist >= -20000:
    return '-10kb_to_-20kb'
  elif dist < -5000 and dist >= -10000:
    return '-5kb_to_-10kb'
  elif dist < -1000 and dist >= -5000:
    return '-1kb_to_-5kb'
  elif dist < 0 and dist >= -1000:
    return 'TSS_to_-1kb'
  elif dist == 0:
    return 'Gene'
  elif dist > 0 and dist <= 1000:
    return 'TES_to_+1kb'
  elif dist > 1000 and dist <= 5000:
    return '+1kb_to_+5kb'
  elif dist > 5000 and dist <= 10000:
    return '+5kb_to_+10kb'
  elif dist > 10000 and dist <= 20000:
    return '+10kb_to_+20kb'
  elif dist > 20000 and dist <= 50000:
    return '+20kb_to_+50kb'
  elif dist > 50000:
    return '+50kb_to_+Inf'

for in_line in sys.stdin:
  (bin_str, upstream_gene_str, upstream_gene_dist, downstream_gene_str, downstream_gene_dist) = in_line.rstrip().split('|')
  try:
    upstream_gene_dist = int(upstream_gene_dist)
  except ValueError:
    upstream_gene_dist = np.nan
  try:
    downstream_gene_dist = int(downstream_gene_dist)
  except ValueError:
    downstream_gene_dist = np.nan
  upstream_gene_elems = upstream_gene_str.split('\t')
  downstream_gene_elems = downstream_gene_str.split('\t')

  '''
  We have six possible scenarios, depending on the minimum distance and strandedness of closest elements
  '''

  '''
  1) If either distance value is NaN, the bin is assigned to the non-NaN distance element
  '''
  if upstream_gene_dist is np.nan and downstream_gene_dist is np.nan:
    raise SystemError("Bad NaN input on:\n[{}]\n".format(in_line.rstrip()))
  elif upstream_gene_dist is np.nan:
    gene_id = downstream_gene_elems[3]
    if downstream_gene_elems[5] == '+':
      # bin is closer to downstream TSS
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(-(abs(downstream_gene_dist))), gene_id))
    else:
      # bin is closer to downstream TES
      sys.stdout.write('{}\t{}\t\n'.format(bin_str, assign_to_category(abs(downstream_gene_dist)), gene_id))
    continue
  elif downstream_gene_dist is np.nan:
    gene_id = upstream_gene_elems[3]
    if upstream_gene_elems[5] == '+':
      # bin is closer to upstream TES
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(abs(upstream_gene_dist)), gene_id))
    else:
      # bin is closer to upstream TSS
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(-(abs(upstream_gene_dist))), gene_id))
    continue

  '''
  2) If either distance value is zero, the bin is assigned to "Gene" to state overlap with a gene
  '''
  if upstream_gene_dist == 0:
    gene_id = upstream_gene_elems[3]
    sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(0), gene_id))
    continue

  if downstream_gene_dist == 0:
    gene_id = downstream_gene_elems[3]
    sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(0), gene_id))
    continue

  '''
  3) If the upstream strand is positive, if the downstream strand is positive, then:

   -- the upstream distance is a signed (-) distance to the upstream gene's TES
   -- the downstream distance is a signed (+) distance to the downstream gene's TSS

  '''
  if upstream_gene_elems[5] == '+' and downstream_gene_elems[5] == '+':
    if abs(upstream_gene_dist) < abs(downstream_gene_dist):
      # bin is closer to upstream TES
      gene_id = upstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(abs(upstream_gene_dist)), gene_id))
    else:
      # bin is closer to downstream TSS
      gene_id = downstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(-(abs(downstream_gene_dist))), gene_id))
    continue

  '''
  4) If the upstream strand is negative, if the downstream strand is positive, then:

   -- the upstream distance is a signed (-) distance to the upstream gene's TSS
   -- the downstream distance is a signed (+) distance to the downstream gene's TSS

  '''
  if upstream_gene_elems[5] == '-' and downstream_gene_elems[5] == '+':
    if abs(upstream_gene_dist) < abs(downstream_gene_dist):
      # bin is closer to upstream TSS
      gene_id = upstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(-(abs(upstream_gene_dist))), gene_id))
    else:
      # bin is closer to downstream TSS
      gene_id = downstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(-(abs(downstream_gene_dist))), gene_id))
    continue

  '''
  5) If the upstream strand is positive, if the downstream strand is negative, then:

   -- the upstream distance is a signed (-) distance to the upstream gene's TES
   -- the downstream distance is a signed (+) distance to the downstream gene's TES

  '''
  if upstream_gene_elems[5] == '+' and downstream_gene_elems[5] == '-':
    if abs(upstream_gene_dist) < abs(downstream_gene_dist):
      # bin is closer to upstream TES
      gene_id = upstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(abs(upstream_gene_dist)), gene_id))
    else:
      # bin is closer to downstream TES
      gene_id = downstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(abs(downstream_gene_dist)), gene_id))
    continue

  '''
  6) If the upstream strand is negative, if the downstream strand is negative, then:

   -- the upstream distance is a signed (-) distance to the upstream gene's TES
   -- the downstream distance is a signed (+) distance to the downstream gene's TES

  '''
  if upstream_gene_elems[5] == '-' and downstream_gene_elems[5] == '-':
    if abs(upstream_gene_dist) < abs(downstream_gene_dist):
      # bin is closer to upstream TSS
      gene_id = upstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(-(abs(upstream_gene_dist))), gene_id))
    else:
      # bin is closer to downstream TES
      gene_id = downstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(abs(downstream_gene_dist)), gene_id))
    continue

  '''
  Sanity check: If we get here, something went wrong
  '''
  raise SystemError("Bad logic on:\n[{}]\n".format(in_line.rstrip()))

```

To use it, e.g.:

```
$ MP_LEFT=500
$ MP_RIGHT=-499
$ MP_LEFT_UNDO=-500
$ MP_RIGHT_UNDO=499
$ CATEGORIZED_BINS=${PWD}/../data/categorized_bins.${ASSEMBLY}.bed
$ ./assign_chops.py < <(bedops --range ${MP_LEFT}:${MP_RIGHT} --everything ${CHOPPED_GENOME} | closest-features --dist - ${GENCODE_GENES}) | sort-bed - | bedops --range ${MP_LEFT_UNDO}:${MP_RIGHT_UNDO} --everything - > ${CATEGORIZED_BINS}

```

The output from this file annotates bins with the closest feature to a gene. It looks something like this:

```
$ more categorized_bins.hg38.bed
...
chr1  32000 33000 TES_to_+1kb RP11-34P13.3
chr1  33000 34000 TES_to_+1kb FAM138A
chr1  34000 35000 Gene  FAM138A
chr1  35000 36000 Gene  FAM138A
chr1  36000 37000 Gene  FAM138A
chr1  37000 38000 TSS_to_-1kb FAM138A
chr1  38000 39000 -1kb_to_-5kb  FAM138A
...

```

For example, the region `chr1:33000-34000` is closest to the 1kb window downstream of the TES of `FAM138A`.

Once you have this file, you can quickly map it to your sorted regions-of-interest file:

```
$ bedmap --echo --echo-map --delim '\t' regions-of-interest.bed categorized_bins.hg38.bed > answer.bed

```

The file `answer.bed` will contain your regions of interest and the annotated bins which overlap them.
> 1 votes

# Answer
> How can I quickly find out if it lies in an intron or exon of a gene? Or it is in the upstream or downstream regions of a gene? Is there a tool or web service from where I can get this information?

Since you emphasise speed and suggest a web service, I'm guessing you just want to check a few individual loci. In that case, just use a genome browser like those from UCSC or Ensembl. If you look up those coordinates in UCSC you can see the sequence falls within the last exon of the ZFHX4 protein coding gene and other information.
> 0 votes

# Answer
There is already some very useful answers out there but in my experience Homer does a very efficient job at annotating bed files. This gives you a lot of details about nearest transcription start sites or exonic regions and more.

Be sure to convert your data to bed format first
> 0 votes

---
Tags: chromosomes, hg19
---


# Question
Title: Generate an appropriate output from DNaSP6 to Arlequin program
I am performing a haplotype and nucleotide diversity analysis of my Sanger sequences (from Watermelon mosaic virus coat protein region) using the Arlequin program. My aim is to see the probability that two randomly sampled alleles are different, and the expansion and genetic variation of the virus respectively.

Some of my sequences have gaps (-) because the PCR reaction could not cover all my region of interest. After making the alignment in Mega and export in a Nexus file to DnaSp6 I saw that there are many options to export to the Arlequin program after DnaSP6. I do not know what is the best option to generate the output from DNAsp6 to Arlequin because I must choose one option for the sites with the alignment gaps option and another option for the invariable sites option.

These are the options:

```
Sites with Alignment Gaps options:

Not considered: These sites are ignored (complete deletion).
Considered: Gaps are considered just like another nucleotide variant (fifth state).
Only gaps are considered: Only gaps information is considered to build haplotypes.

Invariable Sites options:

Removed: Invariable (monomorphic) sites will not be included in the output file.
Included: Invariable (monomorphic) sites will be included in the output file.

```

From my point of view I think that it is more certain for my analysis (above explained) to use the Not consider gaps option (because they can buffer the information of the other sequences) and include invariable sites option because I want to analyse these parameters of my dataset.

Is it sensible? I would like to know what is the best option for my dataset.

# Answer
Yes, your question makes alot of sense and its a good question because hits at the heart of differences between population genetic analysis and phylogenetic analysis. However, when you speak about gaps you are referring to 'indels' (insertion/deletion).

**Indels** The conservative approach is to simply delete sites with gaps. Arlequin uses a lot of maximum likelihood solutions and gaps would automatically be deleted in an Arlequin max. likelihood analysis, because it is the output per alignment per matrix. If the matrix doesn't have a 'gap'.

**Indels and pop gen** Thats tricky, it is better to remove them. HOWEVER the indel may contain underlying population signal, if its a maximum likelihood analysis (pop gen does this too) I suspect they will be ignored. If its pairwise, such Fst including indels might give you a better answer because there is more polymorphisms to assess population structure. So if you're doing a pairwise analysis, I'd do it both ways and then report the methodology correctly if there are differences in results and simply explain the theory.

**Invariable sites** This is different from an indel. An invariable site is a site which carries no phylogenetic signal, it carries information but a phylogenetic analysis could draw any conclusions because all the nucleotides/amino acids are identical at a give alignment position.

The output you obtain in allele frequency will be different and here's the rub ...

**Phylogenetics analysis** Arlequin is really pop gen, but it will contain some phylogenetics programs. If you see anything with 'Jukes Cantor correction', 'Kimura 2 (or 3) parameter model', 'HGW model', '"REV" model' (eek I can't remember its name any its a reverible 6-parameter model of nucleotides), you are dealing with phylogenetics. Keeping invariable sites is normal here because the rate matrix is calculation "per alignment site".

**Population genetics analysis** Fst, Fis, linkage disequilibrium statistics, REMOVE THEM. Why? Because you want SNPs (single nucleotide polymorphisms). The invariable sites will get in the way of the statistics. For example is the locus recombines and you want to assess a linkage disequilibrium statistic (recombination), invariable sites will recombine, but the algorithm can't see this because the nucleotides are all the same. The only thing the algorithm can see and analyse are SNPs, which 'representative' of all the surrounding invariant sites. Thus removing invariants is essential. Arelquin might do this automatically (but I don't know the preparatory algorithms), but if you submit multiple loci Arlquin will want to know which are different loci. It will then recode the loci automatically if a single locus comprises more than one SNP. We are going into the underlying theory of linkage disequilibrium here and I'm not sure if say Fst/Fis is particularly bothered with this or is just concerned with the populations of genetic diversity.

**Summary** So gaps ... if they are genuine indels and not just 'N', ie indeterminant nucleotide calls, remove them for phylogenetics. If its pop gen do the calculation both ways and assess the output, if there is no difference report the 'gap/indel removal' result.

Invariant sites .. for pop gen you remove them but Arlequin might do this automatically to prepare a data set for analysis (I don't know). For phylogenetics you never remove invariant sites.

Goodluck
> 1 votes

---
Tags: phylogeny, phylogenetics, haplotypes
---


# Question
Title: Can I use BioEdit to quantify the sequence similarity?
I am using `BioEdit`. It can perform an align of multiple sequences. Can this sequence similarity be outputted as a value, e.g. 70%?

# Answer
> 1 votes
No Bioedit is an alignment editor only. It is not in anyway an analysis program.

Dump the format in fasta or mega (its the option is available) and import it into Mega X here. You can quantify the diversity as a pairwise matrix, population analysis (if you have a clear a priori that you can subdivide your sequences into different groups) or phylogeny ... which provides diversity as mutations per alignment position. Its pointy and clicky and easy to use.

---
Tags: sequence-alignment, phylogeny, bioedit
---


title: efetch nucleotide -> protein ids

text:
# Question
Title: efetch nucleotide -> protein ids
I'm downloading translated viral genomes initially via Blast - which had shortcomings - and now by `efetch`.

What I want to do is obtain a complete list of all translated protein ids from a nucleotide Genbank entry.

```
efetch -db nucleotide -format gpc -id MT039887 | xtract -insd ....

```

The protein id and the translated gene sequences are within the `efetch -db nucleotide -format gpc -id MT039887` output ... but how do I use `xtract` to parse them?

# Answer
If you want the cds in fasta format, try:

```
esearch -db nucleotide -query MT039887 | efetch -format fasta_cds_na

```

And for the proteins sequence in fasta format:

```
elink -db nucleotide -target protein -id MT039887 | efetch -format fasta

```
> 1 votes

# Answer
In case you are still interested in just getting the protein ids and not sequence, you can use `xtract` like this:

```
efetch -db nucleotide -format gpc -id MT039887 | xtract  -insd CDS gene product protein_id
MT039887.1  orf1ab  orf1ab polyprotein           QHZ00388.1
MT039887.1  S       surface glycoprotein         QHZ00389.1
MT039887.1  orf3a   orf3a protein                QHZ00390.1
MT039887.1  E       envelope protein             QHZ00391.1
MT039887.1  M       membrane glycoprotein        QHZ00392.1
MT039887.1  orf6    orf6 protein                 QHZ00393.1
MT039887.1  orf7a   orf7a protein                QHZ00394.1
MT039887.1  orf8    orf8 protein                 QHZ00395.1
MT039887.1  N       nucleocapsid phosphoprotein  QHZ00396.1
MT039887.1  orf10   orf10 protein                QHZ00397.1

```
> 1 votes

---
Tags: genome, phylogeny, phylogenetics, efetch
---


# Question
Title: Can someone recommend software for reference-based genotype-estimation?
**Problem**: I have low-coverage sequencing data from humans in vcf format. For each sample, I have PL scores / preliminary genotype calls from the software ATLAS (it's a bit like GATK for those who don't know). Because the samples are low coverage (mean 1x, with lots of missing data), I am looking to enhance the genotype calls by using a reference panel and genetic map, before go on to phase the samples. Briefly, the reference panel and genetic map can be used as a kind of prior to inform one about genotypes when there is little data present.

**What I've done so far**: So far I have been using Beagle4 with a command:

```
beagle4 gl=lowCoverageSamples.vcf ref=highCoverageReference.vcf map=geneticMap out=lowCoverageSamples.newCalls.vcf

```

However, Beagle4 is extremely slow, very memory inefficient and quite old. There have been many advances in phasing methods, which are now much more accurate, faster and memory efficient than the beagle4 model. However, I haven't come across any similar improved methods for the genotype re-estimation step.

Does anyone knew of any newer pieces of software which do such a thing, which may more faster and more accurate than Beagle4. It needs to have 2 things in particular:

1) Accept data input in PL score format, as to allow the appropriate representation of uncertainty in low coverage samples.

2) Produce a genotype call / genotype probability for each SNP in each individual as output.

# Answer
> 1 votes
In case anyone was interested, a new piece of software for rapidly imputing and phasing low coverage sequence data from genotype likelihoods has just been released. It is called GLIMSPE and is exactly what I was after.

It looks to be much faster (1200x speedup) compared to Beagle 4.0, and is more accurate at imputing low frequency variants. It is also more amenable to use with larger reference panels (such as the Haplotype Reference Consortium panel).

---
Tags: vcf, genotyping, phasing
---


# Question
Title: Is there a single cell RNAseq equivalent of GTEx or TCGA?
Or do I need to find individual studies and obtain data the long way!

# Answer
I think the most well-known scRNA dataset that is similar to GTEx in the Human Cell Atlas (https://www.humancellatlas.org). It is pretty much still ongoing but there are some preliminary data.
> 2 votes

---
Tags: rna-seq, scrnaseq, single-cell, gtex
---


# Question
Title: BERT Language Model and Gene Sequences - How Do I Relate Clusters of Sequences?
I hope you'll indulge a question from a computer scientist with limited bioinformatics knowledge. I've been working with the Google tool for language modeling called BERT. It's generally regarded as state of the art when encoding language. It allows you to say that the sentences "I went to the bank to cash my check." and "My wife is going to the Credit Union to pay the mortgage" are more similar than if compared to "he was fishing and fell off the bank." In short, BERT takes sentences and encodes them so that related **concepts** are *near* each other and unrelated concepts are *farther* apart.

BERT is an amazing tool and I wondered if it could successfully model the language of the genome if it were encoded into "sentences" (I chose 7-mer words and sentences of length 64, arbitrarily). Below is the result of encoding the full genomes of human, bonobo, gorilla, chimp, and orangutan. There are about 270K "sentences" here. You can see there are four major blobs, with indications of small tight clusters. BERT has decided that the points indicated in red are somehow semantically related. I think is is significant, or at the very least interesting.

My first thought was that clustering was trivial degenerate cases. Maybe all the "AAAAAAA..." were put near each other, and all the "TTTTTTTTT...." etc. I printed out samples from the clusters and visually that is not the case. Also I've made histograms of the average Levenshtein distance in groups. Within a tight group, the LD is statistically much lower (but not trivially low) as compared to unclustered regions (the big blob is all unclustered).

If there is a biological basis for the clustering, then new sequences could be similarly transformed, and classified according to cluster membership (and related biological function). This exhausts my biological knowledge. I've taken sequences from clusters and used BLAST to try and get some insight but I'm **way** out of my comfort zone with that.

My primary question is, given a set of sequences (in my case len=448), how do you figure out what they all have in common (given my hypothesis is that BERT says they are related based on proximity)? And my backup question is, do you think this is interesting or a waste of time?

# Answer
> 2 votes
The UMAP results look very pretty :)

If I were to try to make sense of this blob, I'd try to do the following.

1. Pull out Gene Ontology (GO) terms for a list of candidate genes (you can try to subset Uniprot / Swissprot for some examples). Then I feed through some genes through your fragmented BERT kmer approach to see where they land. I'd be particularly interested in seeing where transmembrane proteins and enzymes land here.
2. I'd try to color this by introns vs exons (you should be able to extract this sort of information from your genome files). There is a good chance that the difference between introns / exons could explain some of these differences.
3. It maybe worth looking into genes that only encode RNAs vs genes that encode proteins.

# Answer
> 1 votes
In truth deep neural networks quite mainstream now, albeit they are a pain because of their reliance on very expensive graphics cards to process the data (GPU) rather than traditional CPU.

There are two issues,

* The algorhithms need an application to a biological phenomena and this is dependent on the model. DNA applications (usually) use a different model to language applications (below)
* The model is only as good what it is trained against. In languages this would be people saying stuff, like in your example. In biology it is a phenotype. Thus you need to feed it the right answer and all variants of the right answer. It will then spot the pattern and recognise this pattern in de novo inputs.

What we are talking is deep learning using a neural network yes we use them, but again the models are (usually) different. Convolution neural networks (CNN) are common in DNA/genomics because it assess the physical linkage of DNA. CNN is also used in image recognition, so like when you type 'dog' into google and see dog images, chances are the search is a consequence of CNN classification.

In languages recurrent neural networks particularly long short term memory (LSTM) models are common because they need to learn the patterns and idioms of a language. Google of course will only use their development called TensorFlow to build their models, which is very geeky, for example they have a programming construct called a "placeholder". In my view a placeholder is something found in Adobe PhotoShop rather than a programming language etc... TensorFlow exponents outside Google tend to use the Keras interface. In Python there is Natural Language Toolkit — NLTK 3.5 to do this.

**Summary** So, to answer your question, by dumping DNA sequence into a language RNN its fun, but it wouldn't mean very much because the evolution of language, or construction of a given language (i.e. the training set) is a very different data to mechanism(s) involved in the evolution of DNA.

Fundamentally, we don't speak DNA.

I think it would be fun dumping whale songs and dolphin squeaks into these programs to see whether there is any convergence with our own speech.

---

My advice is to use an unsupervised training approach to identify your clusters. This would be model-free and training-free and would help further analysis in "supervised-learning" OR just stick with the clusters it generates.

---
Tags: genome, blast, machine-learning
---


# Question
Title: Map gene IDs to protein structure if exist
I want to map some gene ids to their proteins structures if they exist. The genes are rows of this CSV. There are other information that may help also such as Product description but that may lead to anbiguities. I am looking for a programmatic way, biopython or pdb API?

# Answer
> 2 votes
Not sure if you want the structure exactly or a handle of some sort, but the best way i know to get ALL your 5000 identifiers to map to uniprot/pdb is uniprot's mapping service. You can either copy your ids straight from excel or parse with python ex.:

```
import pandas as pd
import os

file_name = 'aap7847_Zhang_SM_Table-S5.xlsx'
sheetname = 'Table S5'
df = pd.read_excel(file_name, sheet_name=sheetname, skiprows=1)
ids = list(df['Gene_ID'])

name_of_file_to_save_to = "geneIds.txt"
with open(name_of_file_to_save_to, "w") as newf:
    for id in ids:
        print("{}\n".format(id), file=newf)
    newf.close()

```

You can then save the Uniprot's response and either add it as a column to your file or map to other parts of Elixir ecosystem with SIFTS, for example, as @Matteo Ferla says.

---

In the case that you maybe need a single id's information, you could hit their API too, although this wouldn't work for bulk or at least didn't work for me. You can map the uniprot accession you get from this to PDB or whatever else you need!

```
import requests
def get_uniprot_accession(geneid, format='list'):
    resp = requests.get(
        "https://www.uniprot.org/uniprot/?query={}&sort=score&format={}".format(geneid, format)).text
    if format == 'list':
        resp = resp.strip()
    return resp

```

---
Tags: gene, protein-structure
---


title: how to evaluate neural networks

text:
# Question
Title: how to evaluate neural networks
I have made some architectures of neural networks in order to classify proteins in three categories. I have calculated true positive, true negative, false negative, false positive, Mathews, sensitivity, specificity, accuracy for each category. My question is what is the right way to evaluate the architecture. Except for the c-index, could I calculate a cumulative sensitivity, accuracy etc?or should I evaluate separately for each category? Thank you!

# Answer
> 0 votes
Wow, I would just go for accuarcy as the singular most important measure, but be careful of sampling bias.

If the classification training had an unbiased sample\* and you have an accuracy of \>90%, in general ANN that wouldn't be great, but in bioinformatcs thats really good. If your training comprises sampling bias that gets tricky and you either need to accomodate that in your accuarcy ROC etc ... or else retrain to eliminate the sample bias.

You can account for all of them of course by wrapping everything into a single analysis.

\*, by this I mean that each of the 3 classifications comprises 33% in each training set. If this is badly skewed, this requires consideration and gets complicated. If it is precisely 1/3 each then 'accurary' is the key stat.

---
Tags: proteins, networks
---


title: load the phenotype data for ballgown

text:
# Question
Title: load the phenotype data for ballgown
I am trying to reproduce the work of this paper \[1\], and I have run StringTie successfully, but after that I have to run Ballgown but could not understand this command:

```
>pheno_data = read.csv("geuvadis_phenodata.csv")

```

Actually, I dont know how to make the .csv file of my samples. What phenotype information should this csv contain?

I have checked the supplementary data of this research article but didn't find any .csv file.

How I can make the .csv file for my samples to load the phenotype data for ballgown?

> 8. Load the phenotype data for the samples. An example file called geuvadis\_phenodata.csv is included with the data files for this protocol (ChrX\_data). In general, you will have to create this file yourself. It contains information about your RNA sequencing samples, formatted as illustrated in this csv (comma-separated values) file. Each sample should be described on one row of the file and each column should contain one variable. To read this file into R we use the command read.csv. In this file the values are separated by commas, but if the file were tabdelimited you could use the function read.table.
>
> > pheno\_data = read.csv("geuvadis\_phenodata.csv")

---

1. Pertea *et al*, Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown, *Nature Protocols* volume 11, pages 1650–1667 (2016)

# Answer
> 1 votes
You can easily do it via Notepad++ or even Excel. Just don't forget to save the output file as a `.csv` file (comma-separated value). The info inside should look like (my working example):

```
id,state
lar_control141,control
lar_control142,control
lar_control143,control
lar_probiotic142,probiotic
lar_probiotic143,probiotic
lar_probiotic144,probiotic
lar_probiotic145,probiotic
lar_probiotic146,probiotic

```

That should work

---
Tags: rna-seq, bioconductor, differential-expression
---


# Question
Title: Why this ballgown error (first column of pData does not match the names of the folders) comes out
I am trying to create a ballgown object, and I would like to set a phenodata information. For that I do

```
phenodata <- read.csv('l002.csv')

```

and my l002.csv file contains the following info

```
id,state
RNA-CCC_8N_S5_L002,control
RNA-CCC_8T1_S6_L002,tube
RNA-CCC_12N_S7_L002,control
RNA-CCC_12T_S8_L002,tube
RNA-CCC_13N_S105_L002,control
RNA-CCC_13T_S106_L002,tube
RNA-CCC_23N_S107_L002,control
RNA-CCC_23T_S108_L002,tube

```

A ballgown folder contains folders that have exactly the same names, e.g. RNA-CCC\_8N\_S5\_L002, RNA-CCC\_8T1\_S6\_L002, etc... Then, I use this piece of code for bg object creation

```
bg = ballgown(dataDir = "ballgown", samplePattern = "RNA-CCC_", pData = phenodata )

```

but every time I run it, I have the same error every time:

```
Error in ballgown(dataDir = "ballgown", samplePattern = "RNA-CCC_", pData = phenodata) :
  first column of pData does not match the names of the folders containing the ballgown data.

```

What can be a problem since first column of pData matches exactly the names of the folders containing the ballgown data..

If I get rid of phenodata ballgown object creation works, but I need this phenodata for further manipulations.

# Answer
Well, it looks like if you're using Windows for managing files - that's not a good idea - it sorts names of your ballgown folders in alphabetical order while you need it differently.

That is why the one should use Linux - when you're `ls -ltr` your ballgown folder in Linux it gives you the "right" sequence of folder names to be put inside a ballgown object using pData.

This is what helped me to get rid of this error:

In Linux enter your ballgown folder and do `ls -ltr | awk '{print$9}'` and you will have the **right** sequence that you have to put inside you `.csv` file using Notepad++, indicating what's what. Below's my example compared to the above mentioned initial example:

```
id,state
RNA-CCC_12N_S7_L002,control
RNA-CCC_12T_S8_L002,tube
RNA-CCC_13N_S105_L002,control
RNA-CCC_13T_S106_L002,tube
RNA-CCC_23N_S107_L002,control
RNA-CCC_23T_S108_L002,tube
RNA-CCC_8N_S5_L002,control
RNA-CCC_8T1_S6_L002,tube

```

I don't have my error anymore and after

```
phenodata <- read.csv('my_phenodata_file.csv')
bg = ballgown(dataDir = "ballgown", samplePattern = "RNA-CCC_", pData = phenodata)

```

my results are nicely wrapped

```
Sat May 02 22:11:57 2020
Sat May 02 22:11:57 2020: Reading linking tables
Sat May 02 22:11:59 2020: Reading intron data files
Sat May 02 22:12:04 2020: Merging intron data
Sat May 02 22:12:06 2020: Reading exon data files
Sat May 02 22:12:14 2020: Merging exon data
Sat May 02 22:12:16 2020: Reading transcript data files
Sat May 02 22:12:20 2020: Merging transcript data
Wrapping up the results
Sat May 02 22:12:22 2020

```

The problem is explained here as well

https://www.biostars.org/p/228389/

Hope this will help anyone
> 2 votes

---
Tags: rna-seq, stringtie
---


title: awk shell script

text:
# Question
Title: awk shell script
I have a folder named quant. Inside this folder I have 27 folders named SRR8068516\_quant to SRR8068543\_quant. I want to edit a file inside these folders titled quant.sf and produce an output. I attempted writing a shell script to do so:

```
    #!/bin/bash
for fn in /home/usr/Downloads/salmon-1.2.1_linux_x86_64/quants/SRR/quant.sf{8068516..8068543};
do
samp=`basename ${fn}`
echo "Processing sample ${samp}"
awk 'BEGIN { OFS=FS="\t" } { sub("\\..*", "", $1); print }' quant.sf > equant.sf
done

```

When I run the script it gives me an error saying : **Processing sample quant.sf8068516 awk: cannot open quant.sf (No such file or directory).** What is wrong that I am doing. Kindly help.

# Answer
> 2 votes
You are running the command on the same file, `quant.sf` each time. You need to use your `$fn` variable instead:

```
#!/bin/bash
for fn in /home/usr/Downloads/salmon-1.2.1_linux_x86_64/quants/SRR/quant.sf{8068516..8068543};
do
    samp=`basename ${fn}`
    echo "Processing sample ${samp}"
    awk 'BEGIN { OFS=FS="\t" } { sub("\\..*", "", $1); print }' "$fn" > equant.sf
done

```

However, that would also fail since it will be simply overwriting a file called `equant.sf` in your current directory with the results of processing each file.

I can't be sure, since you haven't actually shown us your file structure, but I suspect what you're looking for is actually this:

```
#!/bin/bash
for fn in /home/usr/Downloads/salmon-1.2.1_linux_x86_64/quants/SRR/quant.sf{8068516..8068543};
do
    samp=`basename ${fn}`
    outName="${fn/quant.sf/equant.sf}"
    echo "Processing sample ${samp}"
    awk 'BEGIN { OFS=FS="\t" } { sub("\\..*", "", $1); print }' "$fn" > "$outName"
done

```

The `outName="${fn/quant.sf/equant.sf}"` replaces the first occurrence of `quant.sf` in the variable `$fn` with `equant.sf`, keeping the path and thereby ensuring that the output of each file is stored in the same directory with the name `equant.sfNNNNNN`. For example:

```
$ fn=/home/usr/Downloads/salmon-1.2.1_linux_x86_64/quants/SRR/quant.sf8068516
$ echo ${fn/quant.sf/equant.sf}
/home/usr/Downloads/salmon-1.2.1_linux_x86_64/quants/SRR/equant.sf8068516

```

This means that your awk command will become:

```
awk 'BEGIN { OFS=FS="\t" } { sub("\\..*", "", $1); print }' /home/usr/Downloads/salmon-1.2.1_linux_x86_64/quants/SRR/quant.sf8068516 > /home/usr/Downloads/salmon-1.2.1_linux_x86_64/quants/SRR/equant.sf8068516

```

---
Tags: shell, awk
---


# Question
Title: merge the matrix from computematrix of deeptools
I have several matrix from computematrix of deeptools. I need to merge two of them using "computeMatrixOperations cbind -m input1.mat.gz input2.mat.gz -o output.mat.gz" but I am running to error "/computeMatrixOperations.py", line 521, in cbindMatrices hm.matrix.matrix = np.hstack((hm.matrix.matrix, np.empty(hm2.matrix.matrix.shape)))"

Is there any other way to merge the matrixes?

# Answer
> 1 votes
You're trying to cbind two matrices of different sizes, so computeMatrix is throwing an error since it's unclear which row goes with which. You need to ensure that matrices have exactly the same number of rows and that the Nth row in each match (i.e., they should be bound together with cbind).

---
Tags: sequence-alignment, ngs, chip-seq, deeptools
---


# Question
Title: What does the number mean in an HGVSp annotation?
Let's take the example `p.Arg452Pro` that I got from an annotated VCF file that contains an `HGVSp` column.

What does `452` mean in this case ?

Theoretically, is it possible to have the same `HGVSp` on the same gene but on a different transcript level ?

For example:

```
Chromosome      Position        Ensembl_Gene        Ensembl_transcript      HGVSp
4               4613173         ENSG00052342        ENST00141456            p.Arg452Pro
4               4613176         ENSG00052342        ENST00141458            p.Arg452Pro

```

Thank you in advance.

# Answer
It's the position in the protein. `p.Arg452Pro` means that there is a variant that changes an **Arg**inine residue at position 452 of the protein into a **Pro**line.

And yes, it is absolutely possible to have the same HGVS p. notation on the same gene but on different transcripts. More than possible, it is quite common. Your particular example isn't actually a real variant since that position on chromosome 4 is intronic in both hg19 and hg38, but TP53:p.R175L1 is a real example of what you describe:

As you can see in the image or by visiting the link, this variant on TP53 has the same HGVS p. notation on 6 transcripts. It also has a different p. notation on two more transcripts that are not shown in the image as well as no p. notation for a few other transcripts in which it doesn't affect the protein product.

---

1Disclaimer: The link goes to https://varsome.com, a free human variant search engine that is developed by the company I work for.
> 4 votes

---
Tags: vcf, sequence-annotation
---


# Question
Title: Is there any value in scaffolding the output contigs of MEGAHIT assembler given a metagenomic dataset?
As far as I understood, for most assembly programs, the scaffolding step takes into consideration paired-end information in order to get from contigs (contiguous sequences) to scaffolds (longer sequences that might have some N-filled gaps).

My assembly software of choice, MEGAHIT, uses paired-end information to build the contigs, but it does not output a standard scaffold. So I am wondering the following, related things:

-is it meaningful to run a scaffold program on the output of MEGAHIT? I imagine there might be some instances in which paired-end information could span a gap.

-Which software would you recommend for it? (I've tried soapdenovo2 and SSPACE, but they appear not to be actively maintained so I have the issue of 'it doesn't work and I can't do anything about it')

-Could relevant information (eg. two contigs being connected by pair-end information) be recovered by alternative and perhaps more user-friendly means, such as exploring the assembly graph with Bandage?

Thank you for your time!

# Answer
> 1 votes
**Update 2: It looks like your approach has actually been suggested here as one way to use the PE information. I guess megahit may not be really using the PE information anyways. I still believe that it's kind of weird but other people do suggest it, so maybe it's worth trying what Torsten suggests.**

I think that scaffolding with PE information from the reads used as input to Megahit is somewhat sketchy. You could still try it, but I'd be worried about artifacts, just because de novo assemblers are so heuristic.

However, I think that it is perfectly ok to scaffold using *orthogonal* data. Here are some examples of orthogonal data:

1. A reference assembly (weird for metagenomics but still possibly ok).
2. long reads, e.g. ONT or PacBio (low coverage probably fine).
3. Proximity ligation data, e.g. Hi-C reads.
4. Mate-pair reads (also lead to artifacts at a somewhat high rate, but technically valid).
5. Linked-reads (e.g. 10X).
6. Probably other methods as well.

The tools used in each case would be somewhat different. For a little more information about how you might use this, here is a recent review of tech.

Full disclosure: I work for a company that sells Hi-C kits for such applications.

Update: Realized that I missed one part of the Q. I think that visually exploring the assembly using e.g. bandage is always a good idea. Quite possibly you can make some scaffolding decisions that way, but to me it sounds somewhat painful to do, especially in a metagenome where there are going to be a lot of multiple-branching collapsed regions.

# Answer
> 0 votes
You can use SOAPdenovo-Fusion to scaffold contigs produced by MEGAHIT as suggested by one of the developers: https://github.com/aquaskyline/SOAPdenovo2

---
Tags: assembly, metagenome, scaffold
---


# Question
Title: Downloading all COI sequences from BOLD fails
I have metabarcoding sequence data (COI) from bulk animal samples (including arthropoda, nematoda, annelida, mollusca) and I want to BLAST all of these sequences. I used following command to do this: `blastn -remote -db nt -query COI_all.fasta -num_alignments 2 -out COI_blasted.txt`. However this results in errors similar to this post: https://www.biostars.org/p/359971/ .

These errors probably appear due to the number of sequences in my file (around 700) and the remote connection is thus interrupted.

I found that a solution would be to use `blastn` with a local database and since the samples are so diverse, I would like to download ALL animal COI sequences from BOLD (or gen bank). It would not be a problem if non-animal (e.g. plant) sequences would also be included.

I think the BOLD database would be great to BLAST my sequences to. However, I'm currently struggling to find a good way to download all animal `COI` sequences from BOLD.

When entering `COI-5P` as search term on http://v4.boldsystems.org/index.php/Public\_SearchTerms I receive error: `Your search terms resulted in too many matching terms. Please try again with more specific search criteria.`. I could likely download the sequences from all the `phyla` etc seperately and merge them, but I'd rather just download 1 file.

I also tried to use the API by running: `wget http://v4.boldsystems.org/index.php/API_Public/sequence?marker=COI-5P`. A download starts but around 3.7 MB download, it is stuck and the file I receive only contains ~5000 sequences.

**UPDATE: I've contacted BOLD about the stalling behavior and this is their reply: "This issue is because of the large API request that retrieves millions of records, which our system does not handle. Please break up the search by smaller groups, such as classes."**

**Does anyone have a solution to download all COI sequences from BOLD in one file?**

I could also download COI sequences from gen bank using the ftp://ftp.ncbi.nlm.nih.gov/blast/db/ URL, but I'm not sure which exact files I need. For `16S, 18S,..` it is obvious, but not for COI. Any suggestions?

Thanks for the help.

# Answer
I couldn't help but looking a bit deeper into this. There is actually an R package that uses the BOLD API.

However, they also say that there will be timeouts when downloading large sets of data. Luckily, they suggest a workaround (downloading sequences separately for each taxonomic Class).

So you should be able to download the marker sequences using the suggested steps:

```
# install packages (if necessary)
install.packages("bold")
install.packages("taxize")
install.packages("seqinr")
# load packages
library(tidyverse)
library(bold)    # API interface to BOLD
library(taxize)  # for NCBI taxonomy lookup
library(seqinr)  # for FASTA output

# get class-level taxa within "Animalia" from NCBI taxonomy
taxa <- downstream("Animalia", db = "ncbi", downto = "class")
# check if taxa present in BOLD
checks <- bold_tax_name(taxa$Animalia$childtaxa_name)
taxa_bold <- checks[!is.na(checks$taxon),]$taxon

# download sequences from BOLD for each class-level taxon
sequences <- map(taxa_bold[1:3], bold_seq, marker = 'COI-5P') %>%
  flatten() %>%
  bind_rows()

# write sequences to file
write.fasta(
  sequences = as.list(sequences$sequence),
  names = as.list(sequences$id),
  nbchar = 80,
  file.out = 'coi5p.fasta')

```

Note that I selected the first tree taxa (`taxa_bold[1:3]`) for testing, so that it runs quicker. If you want to get all taxa, just remove the `[1:3]`.

One last importand caveat: In the help for `bold_seq` it says that

> Notes from BOLD on the `marker` param: "All markers for a specimen matching the search stringwill be returned. ie. A record with COI-5P and ITS will return sequence data for both markers evenif only COI-5P was specified."You will likely end up with data with markers that you did not request - just be sure to filter thoseout as needed.

So you might want to double-check that the downloaded sequences are actually COI-5P. But if you are using them as a BLAST database, you might also be happy with COI-5P **and** ITS markers being present in the database (it affects the e-values, though).
> 2 votes

# Answer
A straightforward way of getting the sequences from GenBank is from their website using this URL: https://www.ncbi.nlm.nih.gov/nuccore?term=COI\[Gene Name\]

There you can click on "Send to" -\> select "Gene Features" and "FASTA Nucleotide" -\> Click "Create File".

There are \>3 mio hits, so the download could take a while.
> 1 votes

---
Tags: public-databases, data-download, blastn
---


title: What is PyMOLs search path for scripts?

text:
# Question
Title: What is PyMOLs search path for scripts?
I know I can run my PyMOL scripts if I specify a valid path from the current working directory and also using absolute paths, but can I have a central script repository, so I don't always need to specify the full path or copy my scripts around from project to project?

Environment info:

```
Kernel: Linux 5.6.8-1-MANJARO x86_64
Python: 3.8.2
Python2: 2.7.18
PyMOL: 2.3.0 Open-Source

```

# Answer
You can use pymolrc and create aliases to your scripts.

If your script is located e.g. in

```
/Users/tillyboy/pymolscripts/myscript1.pml

```

Just enter this line to your pymolrc:

```
alias myscript1, run /Users/tillyboy/pymolscripts/script1.pml

```

Now myscript1 is available for all your Pymol sessions independent of your current working directory.
> 1 votes

---
Tags: pymol
---


# Question
Title: Is there a way to get the code from "github.mit.edu"?
everyone. I am a postdoc who just begins to analyze single-cell RNA-seq data. Nowadays, I found a really interesting paper (https://doi.org/10.1016/j.celrep.2018.10.047, PMID: 30404002). So, I decided to study bioinformatics by drawing the same graphs presented in the paper during the corona outbreak season. However, it was turned out that I can't get any access to the website in which code is deposited. Here is the address for the website.

Do you guys know how to get the code to learn bioinformatics and make graphs?

Thank you.

# Answer
> 0 votes
It is very strange that they put their code somewhere cannot be easily accessed. If they just put it under a normal Github repo, there would not be such a problem. One way to get the code is simply to email the corresponding author. People are usually responsive if you are interested in their work.

---
Tags: scrnaseq, sequence-analysis
---