# Question
Title: Parsing gtf file for transcript ID and transcript name
I used this previously to get ensembl ID and its gene name from the gencode gtf annotation file

This was the script I used

```
awk '{
    for (i = 1; i <= NF; i++) {
        if ($i ~ /gene_id|gene_name/) {
        printf "%s ", $(i+1)
        }
    }
    print ""
}' Homo_sapiens.GRCh37.70.gtf | sed -e 's/"//g' -e 's/;//g' -e 's/ /\t/' | sort -k1,1 | uniq > Homo_sapiens.GRCh37.70.txt

```

And the input file:

```
##description: evidence-based annotation of the human genome (GRCh38), version 21 (Ensembl 77)
##provider: GENCODE
##contact: gencode@sanger.ac.uk
##format: gtf
##date: 2014-09-29
chr1    HAVANA  gene    11869   14409   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENSG00000223972.5"; gene_type "transcribed_unprocessed_pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "transcribed_unprocessed_pseudogene"; transcript_status "KNOWN"; transcript_name "DDX11L1"; level 2; havana_gene "OTTHUMG00000000961.2";
chr1    HAVANA  transcript  11869   14409   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENST00000456328.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "processed_transcript"; transcript_status "KNOWN"; transcript_name "DDX11L1-002"; level 2; tag "basic"; transcript_support_level "1"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000362751.1";
chr1    HAVANA  exon    11869   12227   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENST00000456328.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "processed_transcript"; transcript_status "KNOWN"; transcript_name "DDX11L1-002"; exon_number 1; exon_id "ENSE00002234944.1"; level 2; tag "basic"; transcript_support_level "1"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000362751.1";
chr1    HAVANA  exon    12613   12721   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENST00000456328.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "processed_transcript"; transcript_status "KNOWN"; transcript_name "DDX11L1-002"; exon_number 2; exon_id "ENSE00003582793.1"; level 2; tag "basic"; transcript_support_level "1"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000362751.1";
chr1    HAVANA  exon    13221   14409   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENST00000456328.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "processed_transcript"; transcript_status "KNOWN"; transcript_name "DDX11L1-002"; exon_number 3; exon_id "ENSE00002312635.1"; level 2; tag "basic"; transcript_support_level "1"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000362751.1";
chr1    HAVANA  transcript  12010   13670   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENST00000450305.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "transcribed_unprocessed_pseudogene"; transcript_status "KNOWN"; transcript_name "DDX11L1-001"; level 2; ont "PGO:0000005"; ont "PGO:0000019"; tag "basic"; transcript_support_level "NA"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000002844.2";
chr1    HAVANA  exon    12010   12057   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENST00000450305.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "transcribed_unprocessed_pseudogene"; transcript_status "KNOWN"; transcript_name "DDX11L1-001"; exon_number 1; exon_id "ENSE00001948541.1"; level 2; ont "PGO:0000005"; ont "PGO:0000019"; tag "basic"; transcript_support_level "NA"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000002844.2";
chr1    HAVANA  exon    12179   12227   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENST00000450305.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "transcribed_unprocessed_pseudogene"; transcript_status "KNOWN"; transcript_name "DDX11L1-001"; exon_number 2; exon_id "ENSE00001671638.2"; level 2; ont "PGO:0000005"; ont "PGO:0000019"; tag "basic"; transcript_support_level "NA"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000002844.2";
chr1    HAVANA  exon    12613   12697   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENST00000450305.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "transcribed_unprocessed_pseudogene"; transcript_status "KNOWN"; transcript_name "DDX11L1-001"; exon_number 3; exon_id "ENSE00001758273.2"; level 2; ont "PGO:0000005"; ont "PGO:0000019"; tag "basic"; transcript_support_level "NA"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000002844.2";
chr1    HAVANA  exon    12975   13052   .   +   .   gene_id "ENSG00000223972.5"; transcript_id "ENST00000450305.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "transcribed_unprocessed_pseudogene"; transcript_status "KNOWN"; transcript_name "DDX11L1-001"; exon_number 4; exon_id "ENSE00001799933.2"; level 2; ont "PGO:0000005"; ont "PGO:0000019"; tag "basic"; transcript_support_level "NA"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000002844.2";

```

What I'm trying to do now is replace the gene\_ID with transcript\_ID and the gene\_name with transcript\_name but I'm doing something wrong as I'm not getting the transcript\_ID but the gene Ensembl ID instead. I'm not sure if I'm choosing the right column for parsing the transcript\_ID and name.

So I'm expecting to get this

```
transcript_id "ENST00000450305.2";  gene_name "DDX11L1";transcript_name "DDX11L1-001"

```

Any help or suggestion would be highly appreciated

# Answer
> 3 votes
I don't see how your `awk` command would work since you're using whitespace as a field delimiter. In any case, you can use a much shorter perl command for this:

```
$ perl -lne 'print "@m" if @m=(/((?:transcript_id|gene_name|transcript_name)\s+\S+)/g);' Homo_sapiens.GRCh37.70.gtf
transcript_id "ENSG00000223972.5"; gene_name "DDX11L1"; transcript_name "DDX11L1";
transcript_id "ENST00000456328.2"; gene_name "DDX11L1"; transcript_name "DDX11L1-002";
transcript_id "ENST00000456328.2"; gene_name "DDX11L1"; transcript_name "DDX11L1-002";
transcript_id "ENST00000456328.2"; gene_name "DDX11L1"; transcript_name "DDX11L1-002";
transcript_id "ENST00000456328.2"; gene_name "DDX11L1"; transcript_name "DDX11L1-002";
transcript_id "ENST00000450305.2"; gene_name "DDX11L1"; transcript_name "DDX11L1-001";
transcript_id "ENST00000450305.2"; gene_name "DDX11L1"; transcript_name "DDX11L1-001";
transcript_id "ENST00000450305.2"; gene_name "DDX11L1"; transcript_name "DDX11L1-001";
transcript_id "ENST00000450305.2"; gene_name "DDX11L1"; transcript_name "DDX11L1-001";
transcript_id "ENST00000450305.2"; gene_name "DDX11L1"; transcript_name "DDX11L1-001";

```

The scriptlet will print the contents of the `@m` array if this line matches what you need. The regex magic is:

* `(?: )`: this is a non-capturing group. We need this so that `@m` only holds what is inside the *capturing* parentheses (see below).
* `( (?:transcript_id|gene_name|transcript_name)\s+\S+)` : the outer parentheses will capture everything matched: one of `transcript_id|gene_name|transcript_name` followed by 1 or more whitespace characters (`\s+`) and one or more non-whitespace characters.
* `print @m if @m=...` : only print if the match was successful.

---
Tags: gtf, text-processing
---


title: Sampling haplotypes

text:
# Question
Title: Sampling haplotypes
I am trying to simulate different genome of peoples, I have data (VCF files) of various genes from the 1000K Gene project.

I want to simulate different whole genomes i.e generate a new population by combining real haplotypes I have. I am wondering what is the best way to approach the problem. That is what is an efficient method of generating realistic (Not just randomly select parts and recombining them) new genotypes based on the real genotypes I already have.

I am using Bioconductor packages VariantAnnotation to read and manipulate the VCF files and TxDb.Hsapiens.UCSC.hg19.knownGene to determine the Genes positions.

The data looks as the following:

```
> gene58@rowRanges
GRanges object with 91 ranges and 1 metadata column:
              seqnames    ranges strand | paramRangeID
                 <Rle> <IRanges>  <Rle> |     <factor>
  rs551585351        1 229566998      * |         <NA>
  rs528384854        1 229567027      * |         <NA>
  rs542093083        1 229567063      * |         <NA>
  rs561849701        1 229567128      * |         <NA>
  rs531042647        1 229567160      * |         <NA>
          ...      ...       ...    ... .          ...
  rs565479298        1 229569784      * |         <NA>
  rs572772527        1 229569785      * |         <NA>
     rs605430        1 229569803      * |         <NA>
     rs605428        1 229569804      * |         <NA>
  rs368699658        1 229569810      * |         <NA>
  -------
  seqinfo: 86 sequences from "hg19" genome

```

The genotype matrix of the gene:

```
> gene.mat[1:5, 1:5]
            HG00867 HG02371 HG00759 HG00766 HG00844
rs551585351 "0|0"   "0|0"   "0|0"   "0|0"   "0|0"
rs528384854 "0|0"   "0|0"   "0|0"   "0|0"   "0|0"
rs542093083 "0|0"   "0|0"   "0|0"   "0|0"   "0|0"
rs561849701 "0|0"   "0|0"   "0|0"   "0|0"   "0|0"
rs531042647 "0|0"   "0|0"   "0|0"   "0|0"   "0|0"

```

# Answer
For simple haplotype simulation of unrelated individuals, there is Montana's venerable HapSim. The model used in that tool is simple and may be good enough for your purposes. You will need to acquire a recombination rate map from for example the HapMap consortium or from 1000G. A cursory google search also revealed this tool and walkthrough that may be suited to your problem: https://adimitromanolakis.github.io/sim1000G/inst/doc/SimulatingFamilyData.html

There are lots of other tools with more advanced models, but it sounds like those would be overkill here.
> 2 votes

# Answer
The data in your matrix appears to be from a VCF file (For reference, see https://github.com/samtools/hts-specs/blob/master/VCFv4.3.pdf).

In the VCF specification, a heterozygous genotype can be specified as `0|1`. The zero and one indicate the allele numbers (0=ref, 1=first alt) that make up the genotype. The "pipe" (`|`) is used to indicate that the alleles are in phase with the previous entry in the file. Therefore, if you observe data like this:

```
"0|1"    "0|0"
"0|1"    "1|1"
"0|1"    "1|1"

```

Then this means that the first sample has the haplotypes \[0,0,0\] and \[1,1,1\]; and that the second sample has the haplotypes \[0,1,1\] and \[0,1,1\].
> 2 votes

# Answer
Ensembl have haplotypes from 1000 Genomes available by transcript.
> 0 votes

---
Tags: bioconductor, vcf, snp, haplotypes, simulated-data
---


# Question
Title: How to get the longest fasta sequence including all possible switching isoforms of a gene out of isoforms
1. I am working with RNA-seq data without a reference genome/transcriptome and am instead using a **Trinity *de novo* transcriptome assembly**.
2. I analyzed both isoform and gene expression abundances using RSEM.
3. And got differentially expressed GENES via DESeq2. In the next step I want to get fasta sequences of differentially expressed genes to use them in the subsequent BlastX/Blast2Go analyses.

Considering the fact that my Trinity created assembly is based on isoform id:

How can I extract the gene sequence out of it? Is that a right thing to combine all the isoforms of a single gene in order to generate the longest possible isoform assembly?

If so, how can I do it to have a union-exon (combined form of all the exons present in the isoforms of an individual gene) sequence used in my Blast analysis?

If not, how should I proceed with the data type available? Trinity assembled file includes the isoform column as:

```
TRINITY_DN1780_c0_g1_i2
TRINITY_DN1780_c0_g1_i1
TRINITY_DN1780_c0_g1_i3
TRINITY_DN1780_c0_g1_i4
TRINITY_DN1780_c0_g1_i5

```

and the DESeq2 produces the output (gene differentials) as:

```
TRINITY_DN1780_c1_g1

```

# Answer
> 2 votes
From your question it isn't clear if you are trying to identify the isoforms that have a very high or very low differential expression value, or if you are trying to correctly find an orthlog for your top hits using BLAST.

If your question is about isoform-level quantification, identifying DE at the isoform level is challenging and DESeq2 is only designed for gene-level analyses. That is why it is only outputting genes rather than a specific isoform.

# Answer
> 2 votes
There's no clear way to go from a trinity "gene" (this is just a cluster of isoforms with shared sequence) to a linear gene sequence that will be blastable and not cause a number of headaches (e.g., due to swapping the order of exons).

What I would suggest that you do is either Blast2GO the isoforms within the Trinity-defined "genes" or take the highest expressed isoform within a "gene" as the sequence you feed into Blast2GO for that gene. Alternatively just use isoforms for everything and avoid this whole issue.

# Answer
> 0 votes
There are a lot of questions in this post, and the main ones have been answered. Here is an answer to the one that has not been answered yet:

**How can I do it to have a union-exon (combined form of all the exons present in the isoforms of an individual gene) sequence used in my Blast analysis?**

I don't think you should do this, but if you want one of the usually many plausible linearized representations of transcript isoforms, then Lace would do the trick.

# Answer
> 0 votes
It appears from the comment in another question that this question is also about transcript-level differential expression analysis.

If instead of using RSEM you use kallisto for transcript quantification, then you could use sleuth to perform transcript-level differential expression analysis.

---
Tags: rna-seq, blast
---


title: Volcano plot in R

text:
# Question
Title: Volcano plot in R
*This question has also been asked on biostars*

How can I reproduce this volcano plot?

I'm only able to do the traditional one, I'm kind knew too these field.

# Answer
This plot is clearly done using core R functions. There are smoother alternatives how to make a pretty volcano plot (like `ggplot` with example here), but if you really wish to, here is my attempt to reproduce it :

I obviously had to generate data since I do not have the expression data from the figure, but the procedure will be about the same with the real data. Here is the R script :

```
### define interesting genes
interesting_genes <- c('Il17s', 'Casz1', 'Rorc', 'Il22')
fc_interesting_genes <- c(-1.8, -1.3, -0.9, -0.85)
pval_interesting_genes <- log10(c(0.0005, 0.01, 1e-4, 4e-05))

### define colours (grey, blue, red)
palette <- c(rgb(0, 0, 0, max = 255, alpha = 65),
    rgb(0, 0, 255, max = 255, alpha = 125),
    rgb(255, 0, 0, max = 255, alpha = 125))

### generate some expression like data, fold change using normal distr, and tranformed uniform for pvalues
fold_changes <- c(rnorm(1000, 0, 2), fc_interesting_genes)
pvalues <- c(log10(runif(1000) / 3^(1 + abs(fold_changes[1:1000]))), pval_interesting_genes)

# find deferentially expressed genes to colour
left_biased <- fold_changes < -1
right_biased <- fold_changes > 1

point_type = 20
point_size = 1.5

### plot corpus
plot(fold_changes, pvalues,
     pch = point_type, cex = point_size, col = palette[1],
     xlab = "Fold Change (log2)", ylab = "P value",
     xlim = c(-6,6), ylim = c(-5, 0), yaxt="n")
aty <- axTicks(2)
labels <- sapply(aty,function(i)
            as.expression(bquote(10^ .(i)))
          )
labels[6] <- 1
axis(2,at=aty,labels=labels)

### plot deferentially expressed genes in colour + genes on interest with black
points(fold_changes[left_biased], pvalues[left_biased],
       pch = point_type, cex = point_size, col = palette[2])
points(fold_changes[right_biased], pvalues[right_biased],
       pch = point_type, cex = point_size, col = palette[3])
points(fc_interesting_genes, pval_interesting_genes, pch = point_type, cex = point_size)

### add text
labelx <- c(-4, -4, 4, 4)
labely <- c(-1.5, -0.5, -0.5, -1.5)
labels <- c('Il17s', 'Casz1', 'Rorc', 'Il22')
text(labelx, labely, labels)
text(-4, -4, "down\n194 genes")
text(4, -4, "up\n156 genes")

### add lines between text and data points
for(i in 1:4){
    lines(c(labelx[i], fc_interesting_genes[i]), c(labely[i], pval_interesting_genes[i]))
}

```
> 13 votes

# Answer
Without a minimal reproducible example (MWE), I can't reproduce the plot but I would suggest using existing `volcanoplot` functions such as the `limma` package on Bioconductor. Here is a typical workflow for a differential expression analysis that produces a violinplot:

```
#install package
source("https://bioconductor.org/biocLite.R")
biocLite("limma")

#load package
library("limma")

#define phenotype/experimental conditions for samples
group = rep(c("Treatment","Control"),
c(number_of_treated_samples,number_of_controls))
design = model.matrix(~group);
colnames(design) = c("Control","TreatmentvsControl")

#fit empirical Bayes model
fit <- lmFit(NormalisedExpression, design)
fit <- eBayes(fit)

#summarise significant results
tt <- topTable(fit, coef="TreatmentvsControl",
                adjust="BH",n=nrow(NormalisedExpression))

#plot results
volcanoplot(fit, coef="TreatmentvsControl")

```

This plot can be customised in a similar manner to base R plots by passing the relevant arguments as shown in the limma documentation. This is the recommended plot format that readers in the field will be familiar with. While I do not recommend deviating from this standard, it is possible to reverse axes in R and the source code for the plot function can be retrieved from R by entering `volcanoplot` (without any arguments) if you wish to modify it.
> 3 votes

# Answer
Easy to visualize volcano plot using Python. Check this tutorial for inverted volcano plot
> 0 votes

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


# Question
Title: Snakemake and SLURM -unable to use it with with --cluster
I'm trying to automate some of my previous bash scripts with snakemake. Below is one of the recipe I created:

```
import os
shell.executable("/bin/bash")

workdir: "/proj/uppstore2018136/publication/jupyterNotebooks"

ENSEMBL_RELEASE='94'
ENSEMBL_GRCm38_BASE='ftp://ftp.ensembl.org/pub/release-'+ENSEMBL_RELEASE+'/fasta/mus_musculus/dna'
ENSEMBL_GRCm38_GTF_BASE='ftp://ftp.ensembl.org/pub/release-'+ENSEMBL_RELEASE+'/gtf/mus_musculus'
GTF_FILE='Mus_musculus.GRCm38.'+ENSEMBL_RELEASE+'.gtf'

rule all:
    input:
        "genome/genome.fa",
        "genome/"+GTF_FILE,
        "indexes/bowtie2/reference.1.bt2",
        "indexes/bowtie2/reference.2.bt2",
        "indexes/bowtie2/reference.3.bt2",
        "indexes/bowtie2/reference.4.bt2",
        "indexes/bowtie2/reference.rev1.bt2",
        "indexes/bowtie2/reference.rev2.bt2",
        "indexes/hisat2/reference.1.ht2",
        "indexes/hisat2/reference.2.ht2",
        "indexes/hisat2/reference.3.ht2",
        "indexes/hisat2/reference.4.ht2",
        "indexes/hisat2/reference.5.ht2",
        "indexes/hisat2/reference.6.ht2",
        "indexes/hisat2/reference.7.ht2",
        "indexes/hisat2/reference.8.ht2"

rule getGenome:
    output:
        output1="genome/genome.fa",
        output2="genome/"+GTF_FILE
    shell:
           '''
           curl -o genome/Mus_musculus.GRCm38.dna.primary_assembly.fa.gz {ENSEMBL_GRCm38_BASE}/Mus_musculus.GRCm38.dna.primary_assembly.fa.gz
           gunzip genome/Mus_musculus.GRCm38.dna.primary_assembly.fa.gz
           mv genome/Mus_musculus.GRCm38.dna.primary_assembly.fa genome/genome.fa
           curl -o genome/{GTF_FILE}.gz {ENSEMBL_GRCm38_GTF_BASE}/{GTF_FILE}.gz
           gunzip genome/{GTF_FILE}.gz
           '''

rule bowtie2Build:
    input: "genome/genome.fa"
    params:
        basename="indexes/bowtie2/reference"
    output:
        output1="indexes/bowtie2/reference.1.bt2",
        output2="indexes/bowtie2/reference.2.bt2",
        output3="indexes/bowtie2/reference.3.bt2",
        output4="indexes/bowtie2/reference.4.bt2",
        outputrev1="indexes/bowtie2/reference.rev1.bt2",
        outputrev2="indexes/bowtie2/reference.rev2.bt2"
    threads: 20
    shell:
        '''
        module load bowtie2
        bowtie2-build --threads {threads} {input} {params.basename}
        '''

rule hisat2Build:
    input: "genome/genome.fa",
        "genome/"+GTF_FILE
    params:
        basename="indexes/hisat2/reference"
    output:
        output1="indexes/hisat2/reference.1.ht2",
        output2="indexes/hisat2/reference.2.ht2",
        output3="indexes/hisat2/reference.3.ht2",
        output4="indexes/hisat2/reference.4.ht2",
        output5="indexes/hisat2/reference.5.ht2",
        output6="indexes/hisat2/reference.6.ht2",
        output7="indexes/hisat2/reference.7.ht2",
        output8="indexes/hisat2/reference.8.ht2",
    threads: 20
    shell:
        '''
        module load HISAT2
        hisat2_extract_splice_sites.py genome/{GTF_FILE} > genome/genome.ss
        hisat2_extract_exons.py genome/{GTF_FILE} > genome/genome.exon
        hisat2-build -f -p {threads} {input} --ss genome/genome.ss --exon genome/genome.exon {params.basename}
        rm genome/genome.ss
        rm genome/genome.exon
        '''

rule clean:
    shell:
        '''
        rm genome/*
        rm indexes/bowtie2/*
        rm indexes/hisat2/*
        '''

```

And this is my driver script which I'm submitting as a job itself:

```
#!/bin/bash

module load snakemake

snakemake --jobs 1 --local-cores 20 --cluster-config /crex/proj/uppstore2018136/publication/jupyterNotebooks/recipes/cluster.json --latency-wait 120 --max-jobs-per-second 1 -c "sbatch -A {cluster.account} -p {cluster.partition} -n {-cluster.n} -t {cluster.time} -C {cluster.C}" -s indexMake

```

But I'm getting this error:

```
Building DAG of jobs...
Using shell: /bin/bash
Provided cluster nodes: 1
Job counts:
    count   jobs
    1   all
    1   bowtie2Build
    1   getGenome
    1   hisat2Build
    4

rule getGenome:
    output: genome/genome.fa, genome/Mus_musculus.GRCm38.94.gtf
    jobid: 2

RuleException in line 30 of /crex/proj/uppstore2018136/publication/jupyterNotebooks/recipes/indexMake:
NameError: The name '-cluster' is unknown in this context. Please make sure that you defined that variable. Also note that braces not used for variable access have to be escaped by repeating them, i.e. {{print $1}

}

```

I really am not sure what I'm doing wrong here ? Is there a conflict with variables or there is something wrong with the snakemake version I'm using (4.5.0).

# Answer
In your submission script you have `-n {-cluster.n}`, hence the error from snakemake I guess. Are you sure it is not a typo and it should be `-n {cluster.n}` instead?
> 5 votes

---
Tags: snakemake
---


# Question
Title: Rdkit.Chem.Draw.savefig() creates empty figures
I can draw a molecule in a jupyter notebook with RDKit. Though, when I use the savefig() method only a white canvas is saved. No molecule is saved.

```
%pylab inline
from rdkit import Chem
from rdkit.Chem import Draw

smiles = 'CCC'
m = Chem.MolFromSmiles(smiles)
fig = Draw.MolToMPL(m)
title('Test')
fig.savefig('mol.jpeg')

```

Any idea what goes wrong here?

# Answer
> 1 votes
For some reason you have to set `bbox_inches` to tight. Without that rdkit does not fit the molecule into the canvas and might plot an empty figure.

```
fig.savefig('mol.jpeg', bbox_inches='tight')

```

---
Tags: python
---


title: How to remove frame on RDKit figures?

text:
# Question
Title: How to remove frame on RDKit figures?
RDKit can plot molecules, thought the structures are surrounded by an ugly frame with ticks. How can I turn that off and plot the molecules without that frame?

Note, the `%pylab inline` portion of this code is from the Jupyter environment.

```
%pylab inline
from rdkit import Chem
from rdkit.Chem import Draw

smiles = 'CCC'
m = Chem.MolFromSmiles(smiles)
fig = Draw.MolToMPL(m)
title('Test')

```

# Answer
It looks like `rdkit` and `rdkit`'s `Chem` just use the `matplotlib`. So, if you import `matplotlib` you can edit your figure like any other `matplotlib` figure.

I produced this from the command line, so it may work differently in the Jupyter environment.

It looks like `Chem` doesn't properly scale the drawing based on how large the chemical is, so I added `plt.xlim` and `plt.ylim` to make the molecule larger in the frame.

```
from rdkit import Chem
from rdkit.Chem import Draw
import matplotlib.pyplot as plt
smiles = 'CCC'
m = Chem.MolFromSmiles(smiles)
fig = Draw.MolToMPL(m)
plt.title('Test')
plt.axis('off')
coord_min = 0.35
coord_max = 0.65
plt.xlim((coord_min, coord_max))
plt.ylim((coord_min, coord_max))
fig.savefig('mol.jpeg', bbox_inches='tight')

```
> 2 votes

---
Tags: python, visualization
---


# Question
Title: How do I interpret the function O(f^-2 log n) in Shekhar et al. (2017)?
The results of Shekhar et al. (2017) "indicate that ASTRAL requires

```
$\mathcal{O}(f^{-2} \log{n})$

```

gene trees to reconstruct the species tree correctly with high probability where `n` is the number of species and `f` is the length of the shortest branch in the species tree.

How do I interpret this function, specifically the big `O` notation, when applying it to my data to estimate whether I have enough loci?

# Answer
Big O notation means "of the order of". Thus $\mathcal{O}(f^{-2} \log{n})$ means that the number of trees required goes with $\log{n}$ for number of species (i.e. if you square the number of genes, you should double the number of trees) and with the inverse square of the shortest branch length (i.e. if you double the shortest branch length, you quarter the number of trees needed).

Big O notation generally allows to to compare two input sizes, but cannot really tell you anything about a single point. I.e. it is hard (impossible?) to tell from this function whether you have sufficient loci.
> 2 votes

---
Tags: phylogenetics, phylogeny, algorithms
---


# Question
Title: How to count the number of mapped read in 100-bp window from a BAM/SAM file
Although I know how to get total number of mapped read using samtools flagstat (`samtools flagstat file_sorted.bam`) but I want to count total number of mapped read in a **non-overlapping sliding window of fixed size** (let's say 500 bp) only with respect to reference (hg19) (assigning each read only once by its start position).

An example of the desired output would be something like:

```
Chromosome    Region_Start     Region_Stop    Num_Unique_Reads
chr2          1                500            200
chr2          501              1000           202
chr2          1001             1501           230
chr2          .                 .              .
chr2          .                 .              .
chr2          .                 .              .
chr2          8001              8500           20

```

Can anyone suggest me how to perform that window operation?

# Answer
## one-liner

Here's a gritty one-liner to count the number of reads in a region if you have just one region that you want to investigate. Change the `1` in `($4 >=1)` and the `500` in `($4 <=500)` to set your window. Change `"hg19"` to your target sequence.

Note, this one-liner does not double-count reads because of `uniq`.

```
samtools view file_sorted.bam | \
    awk '{ if (($3 == "hg19") && ($4 >=1) && ($4 <=500)) { \
  print($1) }}' | sort | uniq | wc -l

```

## python 3

If you need something a little more programmatic to get a known region you should consider trying the `pysam` package.

```
import pysam
samfile = pysam.AlignmentFile("file_sorted.bam", "rb")
region_set = set()
start = 100
stop =  500
counter = 0
for read in samfile.fetch('interesting_contig', 100, 500):
    region_set.add(read.query_name)
    counter += 1
print("{} unique reads in [{}, {}]".format(len(region_set), start, stop))
print("{} alignments in [{}, {}]".format(counter, start, stop))

```

For my `bam` file, this program outputs:

```
1618 unique reads in [100, 500]
2487 alignments in [100, 500]

```

## python 3, non-overlapping windows

It is important to note that the `samfile.fetch()` function conforms to 1-based indexing, consistent with SAM files, unlike pysam's otherwise default 0-based indexing that matches the BAM-indexing.

```
import pysam
samfile = pysam.AlignmentFile("file_sorted.bam", "rb")
print(samfile.references)
print(samfile.lengths)
window = 500
print("ref\tstart\tstop\tuniq_reads")
for i in range(len(samfile.references)):
    refname = samfile.references[i]
    seqlen = samfile.lengths[i]
    for j in range(1, seqlen, window):
        stop = j+window-1 if j+window-1 < samfile.lengths[i] else samfile.lengths[i]
        region_set = set()
        for read in samfile.fetch(refname, j, stop):
            region_set.add(read.query_name)
        print("{}\t{}\t{}\t{}".format(refname, j, stop, len(region_set)))

```

For the small bam file I looked at the output was:

```
ref         start  stop  uniq_reads
OdonB       0      500   1668
OdonB       501    1000  2647
OdonB       1001   1500  129400
OdonB       1501   2000  598
OdonB       2001   2500  538
OdonB       2501   3000  872
OdonB       3001   3500  561
OdonB       3501   4000  49120
OdonB       4001   4500  48433
OdonB       4501   4806  28121

```
> 10 votes

# Answer
using a streaming process with bioalcidaejdk rather than running many random-access jumps.

The script: **the.code**

```
final int WINDOW_SIZE=500;

// map position to count number of reads
final Map<Integer,Long> pos2count = new TreeMap<>();
String currChrom=null;

for(;;)
    {
    // get next read in stream
    final SAMRecord rec = iter.hasNext()?iter.next():null;
    //ignore 'bad' reads
    if(rec!=null ) {
        if(rec.getReadUnmappedFlag()) continue;
        if(rec.getDuplicateReadFlag()) continue;
        if(rec.isSecondaryOrSupplementary()) continue;
        }
    // loop over the windows before the current
    final Iterator<Integer> posIter = pos2count.keySet().iterator();
    while(posIter.hasNext())
        {
        final int start = posIter.next();
        final int end = start + WINDOW_SIZE;
        if(rec!=null && rec.getContig().equals(currChrom) && end >= rec.getStart()) break;
        final long count = pos2count.get(start);
        out.println(currChrom+"\t"+start+"\t"+end+"\t"+count);
        //remove this position from pos2count
        posIter.remove();
        }
    //no more read, exit
    if(rec==null) break;

    //retrieve all the windows for the current read
    for(Integer pos : IntStream.rangeClosed(rec.getStart(), rec.getEnd()).
            mapToObj(POS-> WINDOW_SIZE*(int)(POS/(1.0*WINDOW_SIZE))).
            collect(Collectors.toSet()))
        {
        //increment the count of reads
        final Long count = pos2count.get(pos);
        pos2count.put(pos,1L + (count==null?0L:count));
        }

    currChrom = rec.getContig();
    }

```

execute:

```
$java -jar dist/bioalcidaejdk.jar -f the.code  src/test/resources/S1.bam

RF01    0   500 38
RF01    500 1000    61
RF01    1000    1500    69
RF01    1500    2000    80
RF01    2000    2500    73
RF01    2500    3000    56
RF01    3000    3500    20
RF02    0   500 39
RF02    500 1000    76
RF02    1000    1500    73
RF02    1500    2000    72
RF02    2000    2500    59
RF02    2500    3000    13
RF03    0   500 34
RF03    500 1000    74
RF03    1000    1500    85
RF03    1500    2000    67
RF03    2000    2500    51
RF03    2500    3000    7
RF04    0   500 31
RF04    500 1000    68
RF04    1000    1500    90
RF04    1500    2000    75
RF04    2000    2500    29
RF05    0   500 47
RF05    500 1000    100
RF05    1000    1500    53
RF05    1500    2000    5
RF06    0   500 49
RF06    500 1000    89
RF06    1000    1500    29
RF07    0   500 58
RF07    500 1000    65
RF07    1000    1500    7
RF08    0   500 58
RF08    500 1000    69
RF08    1000    1500    3
RF09    0   500 58
RF09    500 1000    69
RF09    1000    1500    2
RF10    0   500 53
RF10    500 1000    39
RF11    0   500 50
RF11    500 1000    32

```
> 2 votes

---
Tags: sequence-alignment, sam, ngs, samtools, read-mapping
---


# Question
Title: Genomics: Downloading genomic data from the NCBI Genbank database
I'm taking a genomics course at school, and we have a final project due in a few weeks. Now, I can write bash scripts to sort and analyze the data, etc. However, I'm not quite sure how to actually *get* the data from the database. I was hoping to get some help there.

I'd like to get the assemblies of genomes for the species in the *Borrelia* genus. Is there any way that I can do this using commands like "wget", for example. I'm not even sure what URL I'd have to give the computer to go fetch these. Then once I give it the correct URL, how should I tell it what file in that directory to download?

I'd greatly appreciate some help on this issue.

# Answer
For downloading data from commandline use the E-utilities of NCBI.
> 3 votes

---
Tags: genomics, phylogenetics, genbank
---


# Question
Title: How to find all variable-length seqs with an exact 5' and 3' match in a FASTA file
## Context

I am interested in finding all of the promotors specific to a particular sigma factor. I have identified the -35 and -10 sites from the literature, bold denotes -10, -35, binding sites:

```
  -35 site               -10 site
5' TTTACAtatttatttcagacaacGTCTTT 3'
   ^^^^^^                 ^^^^^^

```

The `CAPITALIZED` and `^emphasized^` 5' and 3' sequence components are the conserved sequence that must have an exact match. The lower-case bases are the variable 10-25bp region for which the sequence does not matter.

The sigma factor is from *C. difficile* and I am particularly interested in searching *C. acetobutylicum*, however, it would also be interesting to see where else this sequence crops up. Obviously such a short sequence will arise often by random chance, so it is important to define the spacer region, which might be anywhere between 15-25 nucleotides and may be any sequence.

In particular, I would like to find these binding regions from this sequence and this one.

## Desired output

A text file that contains the locations and sequences of all of the variable-width sequences as defined above. Something like:

```
seq   start  stop  gap  dir  seq
seq1  1      22    10   F    TTTACAagtcagatctGTCTTT
seq2  2      23    11   F    TTTACAagtcagactctGTCTTT
seq2  122    101   10   R    AAAGACatgctagccaTGTAAA

```

Thank you.

# Answer
> 2 votes
## python 3 solution

This python program will get you what I think you want. It is not elegant given that it (A) stores the whole sequence in memory and (B) spends most of its time slicing strings rather than working with a buffer. This program is not recommended for FASTA files for which the largest sequence will not fit well in RAM.

All that you need to do to modify it for other genomes or other binding sites is to change the three lines that begin with `forseq`, `revseq`, and `filename`.

```
from Bio import SeqIO
from Bio.Seq import Seq

forseq = Seq("tttaca".upper())
revseq = Seq("gtcttt".upper())
filename = "NC_003030.fasta"

forcom = forseq.reverse_complement()
revcom = revseq.reverse_complement()
print("seq\tstart\tstop\tgap\tdir\tseq")
for record in SeqIO.parse(filename, "fasta"):
    for i in range(len(record.seq)):
        for v1, v2, direc in [ (forseq, revseq, "F"), (revcom, forcom, "R") ]:
            if record.seq[i:i+len(v1)] == v1:
                for j in range(10,25):
                    if record.seq[i+len(v1)+j:i+len(v2)+j+len(v2)] == v2:
                        print("{}\t{}\t{}\t{}\t{}\t{}".format(
                            record.id, i,
                            i + len(v1)+j+len(v2)-1,
                            j, direc,
                            record.seq[i: i+len(v1) + j + len(v2)] ))

```

Here are the first few lines of example output from the first fasta file that you linked.

```
seq          start    stop     gap  dir  seq
NC_003030.1  97512    97543    20   F    TTTACAAAAGAAAAGATATGCATAATGTCTTT
NC_003030.1  283173   283205   21   F    TTTACACTTATGTTTATATAGTTCAACGTCTTT
NC_003030.1  303104   303125   10   R    AAAGACGCAAGCATAATGTAAA
NC_003030.1  1127691  1127725  23   F    TTTACAAGGTGTTACTGATAGAACAAAGGGTCTTT
NC_003030.1  1301958  1301988  19   F    TTTACATAAAAAGATGATAGTTTTTGTCTTT

```

## python 3 with redundancy

If you need sequence redundancy, this version works. I used the redundancy codes available from www.reverse-complement.com You may use the redundant characters in \[ R, Y, S, W, K, M, B, V, D, H, N \] or lowercase variants. I'm not sure how to handle the '-' character at the moment but am open to suggestions. Code and results below.

Again, the algorithm isn't great given all the for loops so other answers are appreciated!

```
from Bio.Seq import Seq

forseq = "tttaca" # Just use python strings here. No Seq() as above
revseq = "ntnntn"
filename = "NC_003030.fasta"

revcomp = {'A':'T', 'C':'G', 'G':'C', 'T':'A',
           'R':'Y', 'Y':'R', 'S':'S', 'W':'W',
           'K':'M', 'M':'K', 'B':'V', 'V':'B',
           'D':'H', 'H':'D', 'N':'N'}
redund = {'A':['A'], 'C':['C'], 'G':['G'], 'T':['T'],
          'R':['A','G'], 'Y':['C','T'],
          'S':['C','G'], 'W':['A','T'],
          'K':['G','T'], 'M':['A','C'],
          'B':['C','G','T'], 'V':['A','C','G'],
          'D':['A','G','T'], 'H':['A','C','T'],
          'N':['A','C','G','T']}

def revcom(seq):
    return "".join(revcomp.get(base, base) for base in reversed(seq))

def seqmatch(query, reference):
    """Checks if the query sequence (perhaps a promoter with redundant bases)
    matches the reference. Both query and reference must be the same length.

    Returns true if there is a match (redundancy allowed), returns false
    otherwise.
    """
    assert len(query) == len(reference)
    for i in range(len(query)):
        if reference[i] in redund[query[i]]:
            pass
        else:
            return False
    return True

forseq = forseq.upper()
revseq = revseq.upper()
forcom = revcom(forseq)
revcom = revcom(revseq)
print("seq\tstart\tstop\tgap\tdir\tseq")
for record in SeqIO.parse(filename, "fasta"):
    for i in range(len(record.seq)):
        for v1, v2, direc in [ (forseq, revseq, "F"), (revcom, forcom, "R") ]:
            if seqmatch( v1, record.seq[i:i+len(v1)] ):
                for j in range(10,25):
                    if seqmatch( v2, record.seq[i+len(v1)+j:i+len(v2)+j+len(v2)] ):
                        print("{}\t{}\t{}\t{}\t{}\t{}".format(
                            record.id, i,
                            i + len(v1)+j+len(v2)-1,
                            j, direc,
                            record.seq[i: i+len(v1) + j + len(v2)] ))

```

Here is some sample output:

```
seq          start  stop  gap  dir  seq
NC_003030.1  733    767   23   R    GATGAAGATCAAGAAACCGATACAAACAATGTAAA
NC_003030.1  736    767   20   R    GAAGATCAAGAAACCGATACAAACAATGTAAA
NC_003030.1  739    767   17   R    GATCAAGAAACCGATACAAACAATGTAAA
NC_003030.1  742    767   14   R    CAAGAAACCGATACAAACAATGTAAA
NC_003030.1  743    767   13   R    AAGAAACCGATACAAACAATGTAAA
NC_003030.1  2665   2697  21   R    TATTAAAATTTTAGAGGATGATGGTGATGTAAA
NC_003030.1  2668   2697  18   R    TAAAATTTTAGAGGATGATGGTGATGTAAA
NC_003030.1  3798   3822  13   R    CACCAGAAGACTTAAAAATTGTAAA
NC_003030.1  3801   3822  10   R    CAGAAGACTTAAAAATTGTAAA

```

---
Tags: fasta, python, genomics
---


# Question
Title: Derive a GTF containing protein coding genes from a GTF file with Exons and CDS
### Why I need a compatible file

I’m trying to run velocyto with the R package to analyse RNA velocity (cell trajectories) with single cell RNASeq data. I have performed single cell analysis from 10x Genomics data using cellranger.

I have successfully aligned the reads to get loom files and imported these into R. I can get the velocity from these files by following the vignettes. However, I cannot reproduce the RNA velocity analysis based on “gene structure”.

I’m working with a different organism to the example (not human or mouse) so the annotation data provided in the example does not work. I have a GTF file for the latest annotation of this organism. However, it only contains “exon” and “CDS” as features. This appears to be the source of the problem. The “find.ip.sites” function requires a GTF with “features” = “gene” and one of the “attributes” to be “protein\_coding”. These requirements are hard-coded into the velocyto.R function.

I have the following GTF files from the AtRTD2 dataset. The chromosome labels match the loom files in R.

```
Chr1    TAIR10  exon    3631    3913    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  exon    3996    4276    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  exon    4486    4605    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  exon    4706    5095    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  exon    5174    5326    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  exon    5439    5899    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  CDS 3760    3913    .   +   0   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  CDS 3996    4276    .   +   2   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  CDS 4486    4605    .   +   0   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  CDS 4706    5095    .   +   0   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  CDS 5174    5326    .   +   0   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  CDS 5439    5630    .   +   0   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    Araport11   exon    6788    7069    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   exon    7157    7450    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   exon    7564    7649    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   exon    7762    7835    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   exon    7942    7987    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   exon    8236    8325    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   exon    8417    8464    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   exon    8571    8737    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   CDS 7315    7450    .   -   1   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   CDS 7564    7649    .   -   0   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   CDS 7762    7835    .   -   2   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   CDS 7942    7987    .   -   0   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   CDS 8236    8325    .   -   0   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";

```

The specifications for a GTF or GFF2 are:

> ### Fields
>
> Fields must be tab-separated. Also, all but the final field in each feature line must contain a value; "empty" columns should be denoted with a '.'
>
> 1. seqname - name of the chromosome or scaffold; chromosome names can be given with or without the 'chr' prefix. Important note: the seqname must be one used within Ensembl, i.e. a standard chromosome name or an Ensembl identifier such as a scaffold ID, without any additional content such as species or assembly. See the example GFF output below.
> 2. source - name of the program that generated this feature, or the data source (database or project name)
> 3. feature - feature type name, e.g. Gene, Variation, Similarity
> 4. start - Start position of the feature, with sequence numbering starting at 1.
> 5. end - End position of the feature, with sequence numbering starting at 1.
> 6. score - A floating point value.
> 7. strand - defined as + (forward) or - (reverse).
> 8. frame - One of '0', '1' or '2'. '0' indicates that the first base of the feature is the first base of a codon, '1' that the second base is the first base of a codon, and so on..
> 9. attribute - A semicolon-separated list of tag-value pairs, providing additional information about each feature.

### What I’m looking for

Is there a way to derive a compatible GTF file containing protein coding genes from the exons and CDS? I want to produce a GTF which contains genes in the features and protein\_coding in the attributes. Is it possible to do this with existing tools or scripts?

---

### What I’ve tried so far

I can modify the source code of the “find.ip.sites” function to run on my GTF file with these features missing. However, this requires running internal functions from the package written in Rcpp and means my workflow override future updates to the package. Running the rest of the vignette returns errors as no introns long enough have been identified (setting the thresholds lower is also incompatible with the general linear models called). Therefore I think it is better to generate a compatible GTF or GFF3 file rather than alerting the source code of the functions.

While intended for GTF files, the package functions can import GFF3 files, despite being described for GTF in the documentation. I’ve tried generating a GFF3 file using gffread from Cufflinks and replacing “feature” = “mRNA” with “gene” and adding “protein\_coding”. This also returns errors when running the velocity algorithm. It does not work for either the GTF file used as an input for cellranger or the file generated by it. Is there a way to annotate protein coding genes and intron/exon boundaries based on a GTF file containing only exon and CDS annotations?

With cufflinks 2.2.1, a GFF3 was generated

```
gffread - E AtRTD2_19April2016.gtf -o- > AtRTD2_19April2016.gff
sed -i '/mRNA/s/gene_name=AT/gene_type=protein_coding;gene_name=AT/g' AtRTD2_19April2016.gff'
sed -i 's/mRNA/gene/g' AtRTD2_19April2016.gff

```

This is the GFF3 files that I've tried:

```
# gffread - E AtRTD2_19April2016.gtf -o-
##gff-version3
Chr1    TAIR10  gene    3631    5899    .   +   .   ID=AT1G01010.1;geneID=AT1G01010;gene_type=protein_coding;gene_name=AT1G01010
Chr1    TAIR10  exon    3631    3913    .   +   .   Parent=AT1G01010.1
Chr1    TAIR10  exon    3996    4276    .   +   .   Parent=AT1G01010.1
Chr1    TAIR10  exon    4486    4605    .   +   .   Parent=AT1G01010.1
Chr1    TAIR10  exon    4706    5095    .   +   .   Parent=AT1G01010.1
Chr1    TAIR10  exon    5174    5326    .   +   .   Parent=AT1G01010.1
Chr1    TAIR10  exon    5439    5899    .   +   .   Parent=AT1G01010.1
Chr1    TAIR10  CDS 3760    3913    .   +   0   Parent=AT1G01010.1
Chr1    TAIR10  CDS 3996    4276    .   +   2   Parent=AT1G01010.1
Chr1    TAIR10  CDS 4486    4605    .   +   0   Parent=AT1G01010.1
Chr1    TAIR10  CDS 4706    5095    .   +   0   Parent=AT1G01010.1
Chr1    TAIR10  CDS 5174    5326    .   +   0   Parent=AT1G01010.1
Chr1    TAIR10  CDS 5439    5630    .   +   0   Parent=AT1G01010.1
Chr1    Araport11   gene    6788    8737    .   +   .   ID=AT1G01020_P2;geneID=AT1G01020;gene_type=protein_coding;gene_name=AT1G01020
Chr1    Araport11   exon    6788    7069    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   exon    7157    7450    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   exon    7564    7649    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   exon    7762    7835    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   exon    7942    7987    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   exon    8236    8325    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   exon    8417    8464    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   exon    8571    8737    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   CDS 7315    7450    .   -   1   Parent=AT1G01020_P2
Chr1    Araport11   CDS 7564    7649    .   -   0   Parent=AT1G01020_P2
Chr1    Araport11   CDS 7762    7835    .   -   2   Parent=AT1G01020_P2
Chr1    Araport11   CDS 7942    7987    .   -   0   Parent=AT1G01020_P2
Chr1    Araport11   CDS 8236    8325    .   -   0   Parent=AT1G01020_P2

```

The specifications for a GTF or GFF2 are:

The specifications for a GF3 are:

> ### Fields
>
> Fields must be tab-separated. Also, all but the final field in each feature line must contain a value; "empty" columns should be denoted with a '.'
>
> 1. seqid - name of the chromosome or scaffold; chromosome names can be given with or without the 'chr' prefix. Important note: the seq ID must be one used within Ensembl, i.e. a standard chromosome name or an Ensembl identifier such as a scaffold ID, without any additional content such as species or assembly. See the example GFF output below.
> 2. source - name of the program that generated this feature, or the data source (database or project name)
> 3. type - type of feature. Must be a term or accession from the SOFA sequence ontology
> 4. start - Start position of the feature, with sequence numbering starting at 1.
> 5. end - End position of the feature, with sequence numbering starting at 1.
> 6. score - A floating point value.
> 7. strand - defined as + (forward) or - (reverse).
> 8. phase - One of '0', '1' or '2'. '0' indicates that the first base of the feature is the first base of a codon, '1' that the second base is the first base of a codon, and so on..
> 9. attributes - A semicolon-separated list of tag-value pairs, providing additional information about each feature. Some of these tags are predefined, e.g. ID, Name, Alias, Parent - see the GFF documentation for more details.

# Answer
In your case, I would definitely suggest following @Emily\_Ensembl's advice and using the Arabidopsis GTF from Ensembl. But for future reference, if an Ensembl GTF wasn't available, you could build something like this using the `gtf` class from `cgat`

The cgat module and dependancies can be installed by following the instructions here. Specifically, you can download and run their installation script. This will set up a dedicated conda environment to run the versions of dependancies needed. Further documentation can be found here.

```
# download installation script:
curl -O https://raw.githubusercontent.com/cgat-developers/cgat-apps/master/install.sh

# install the development version (recommended, no production version yet):
bash install.sh --devel

```

This requires python3 and anaconda to install. You may need to run this script in the conda environment:

```
# show available environments
conda info --envs
# activate by path
source activate /home/user/local/bin/conda-install/envs/cgat-a

```

Then run the follow as a python3 script. Copy the contents in a file such as `convert_gtf.py`. Then run it in the terminal on your gtf file:

```
python convert_gtf.py genes.gtf > genes_new.gtf

```

```
#python 3.6.4

import sys

from cgat import GTF
from cgatcore import iotools

for gene in GTF.flat_gene_iterator(
     GTF.iterator(iotools.open_file(sys.argv[1])), strict = False):
    gene_start = min (e.start for e in gene)
    gene_end = max(e.end for e in gene)

    if any(e.feature == "CDS" for e in gene):
        gene_type = "protein_coding"
    else:
        gene_type = "non_coding"

    for exon in gene:
        exon_line.setAttribute("gene_biotype", gene_type)

    gene_line = GTF.Entry().fromGTF(gene[0])
    gene_line.feature = "gene"
    gene_line.start = gene_start
    gene_line.end = gene_end

    # Its not clear what the transcript_id for a gene line is.
    # Technically speaking, Ensembl GTFs do not follow the GTF
    # standard as GTF entries must include transcript_ids, but
    # the ensembl ones don't.

    if hasattr(gene_line, "gene_name"):
        gene_line.attributes["gene_name"] = gene[0].gene_name
    else:
        gene_line.attributes["gene_name"] = gene[0].gene_id

    gene_line.attributes["gene_biotype"] = gene_type

    print (str(gene_line)+"\n")
    for exon in gene:
        print(str(exon_line)+"\n")

```
> 4 votes

---
Tags: scrnaseq, sequence-annotation, single-cell, gtf, 10x-genomics
---


title: Snakemake: Cannot find first rule?

text:
# Question
Title: Snakemake: Cannot find first rule?
Back again with another snakemake query. This time I decided to port my read cleaning and alignment pipeline to snakemake. Repeating the steps from previous question and trying not to make any typo here, I'm now receiving this "error" while trying a dry run using `snakemake snakemake -nr -s alignMake`:

```
Building DAG of jobs...
MissingRuleException:
No rule to produce snakemake (if you use input functions make sure that they don't raise unexpected exceptions).

```

Below is my snakefile:

```
import os
shell.executable("/bin/bash")

workdir: '/proj/uppstore2018136/publication/jupyterNotebooks'

pathToPEReads = 'raw_data/reads/BL6/totalRNA'
pathToSMReads = 'raw_data/reads/BL6/smallRNA'

pathToHISAT2index = 'indexes/hisat2/reference'
pathToBOWTIE2index = 'indexes/bowtie2/reference'

pathToTemp = 'temp/'
pathToBAM = 'processed_data/BAM'

readsPE, = glob_wildcards('raw_data/reads/BL6/totalRNA/{sample}_1.fastq.gz')
readsSE, = glob_wildcards('raw_data/reads/BL6/smallRNA/{sample}.fq.gz')

samToolsThreads = 6

rule all:
    input:
        expand('processed_data/BAM/{sample}.bam', sample=readsPE),
        expand('processed_data/BAM/{sample}.bam', sample=readsSE)

rule trimgalorePE:
    input:
        fastq1 = pathToPEReads + '{sample}_1.fastq.gz',
        fastq2 = pathToPEReads + '{sample}_2.fastq.gz'
    output:
        output1 = pathToTemp + '{sample}_1_val_1.fq.gz',
        output2 = pathToTemp + '{sample}_2_val_2.fq.gz'
    shell:
        '''
        module load bioinfo-tools
        module load cutadapt
        module load TrimGalore
        trim_galore --quality 20 --paired --length 80 --illumina --stringency 2 --clip_R1 11 --clip_R2 11 -o {pathToTemp} --paired {input.fastq1} {input.fastq2}
        '''

rule trimgaloreSEandDust:
    input:
        fastq = pathToSMReads + '{sample}.fq.gz',
    output:
        output = pathToTemp + '{sample}_trimmed.fq.gz',
    shell:
        '''
        module load bioinfo-tools
        module load cutadapt
        module load TrimGalore
        trim_galore --quality 20 --length 16 --max_length 32 --small_rna --fastqc --stringency 2 -o {pathToTemp} {input.fastq}
        '''

rule align_hisat2:
    input:
        fastq1 = pathToTemp + '{sample}_1_val_1.fq.gz',
        fastq2 = pathToTemp + '{sample}_2_val_2.fq.gz'
    output:
        output = 'processed_data/BAM/{sample}.bam'
    threads: 10
    shell:
        '''
        module load bioinfo-tools
        module load HISAT2
        module load samtools
        hisat2 -q -p {threads} --no-discordant --no-mixed --rna-strandness RF --dta -x {pathToHISAT2index} -1 {input.fastq1} -2 {input.fastq2} | samtools sort -@ {samToolsThreads} -o {output}
        '''

rule align_bowtie2:
    input:
        fastq = pathToTemp + '{sample}_trimmed.fq.gz'
    output:
        output = 'processed_data/BAM/{sample}.bam'
    threads: 10
    shell:
        '''
        module load bioinfo-tools
        module load bowtie2
        module load samtools
        bowtie2 --no-unal --end-to-end -D 20 -R 3 -N 1 -L 20 -i S,1,0.50 -p {threads} -x {pathToBOWTIE2index} -U {fastq} | samtools sort -@ {samToolsThreads} -o {output}
        '''

rule clean:
    shell:
        '''
        rm temp/*
        '''

```

I believe it is due snakemake being unable to find the first rule ? Maybe there is something wrong with parsing ? I tested the `glob_wildcards` and `expand` function independently and they appear to work properly.

Edit: Here is the directory tree:

```
.
├── citations
├── doc
├── genome
│   ├── genome.fa
│   └── Mus_musculus.GRCm38.94.gtf
├── indexes
│   ├── bowtie2
│   │   ├── reference.1.bt2
│   │   ├── reference.2.bt2
│   │   ├── reference.3.bt2
│   │   ├── reference.4.bt2
│   │   ├── reference.rev.1.bt2
│   │   └── reference.rev.2.bt2
│   └── hisat2
│       ├── reference.1.ht2
│       ├── reference.2.ht2
│       ├── reference.3.ht2
│       ├── reference.4.ht2
│       ├── reference.5.ht2
│       ├── reference.6.ht2
│       ├── reference.7.ht2
│       └── reference.8.ht2
├── LICENSE
├── notebooks
├── plots
├── processed_data
│   └── BAM
├── raw_data
│   └── reads
│       └── BL6
│           ├── smallRNA
│           │   ├── smBL6_1.fq.gz
│           │   ├── smBL6_2.fq.gz
│           │   └── smBL6_3.fq.gz
│           └── totalRNA
│               ├── EKIxxxx012_1.fastq.gz
│               ├── EKIxxxx012_2.fastq.gz
│               ├── EKIxxxx013_1.fastq.gz
│               ├── EKIxxxx013_2.fastq.gz
│               ├── EKIxxxx014_1.fastq.gz
│               ├── EKIxxxx014_2.fastq.gz
│               ├── ESSxxxx077_1.fastq.gz
│               └── ESSxxxx077_2.fastq.gz
├── README.md
├── recipes
│   ├── alignMake
│   ├── cluster.json
│   ├── indexMake
│   ├── indexMakeScript.out
│   └── indexMakeScript.sh
├── results
├── scripts
├── slurmLogs
│   ├── slurm-6575541.out
│   ├── slurm-6575543.out
│   ├── slurm-6575815.out
│   └── slurm-6580867.out
├── source
└── temp

```

# Answer
you are starting `snakemake` the wrong way.

```
snakemake snakemake

```

means "Running the pipeline in `Snakefile` until you reach a rule named `snakemake` or a file name `snakemake` is created."

If your snakefile is called `snakemake` and this is the pipeline you want to start, the correct syntax would be:

```
snakemake -s snakemake

```
> 9 votes

---
Tags: python, snakemake
---


# Question
Title: Error while performing "lift-over" of genomic coordinates between two species using Ensembl Compara API
I'm trying to convert coordinates between the reference mouse genome and CAST/EiJ mice strain. Since Ensebml already has alignment available in the browser, I decided to use Ensembl perl API to lift over the alignment for the regions I'm interested in. The problem is that it doesn't seem to be working. I'm getting following error message:

```
-------------------- EXCEPTION --------------------
MSG: The database do not contain any CACTUS_HAL data for Bio::EnsEMBL::Compara::GenomeDB=HASH(0x557db0d661e8)!
STACK toplevel coord.pl:42
Date (localtime)    = Mon Oct 15 07:42:17 2018
Ensembl API version = 94
---------------------------------------------------

```

which I assume means that I'm either searching for the wrong alignment record or the alignment itself is not available. Following is the perl script which I modified from the Ensembl Compara API tutorial:

```
use strict;
use Bio::EnsEMBL::Registry;
use Bio::EnsEMBL::Utils::Exception qw(throw);
use Bio::AlignIO;
use Getopt::Long;

my $species = "murinae";
my $targetSpecies = "mus_musculus_casteij";
my $alignment_type = "CACTUS_HAL";
my $output_file = undef;
my $coord_system = "chromosome";
my $seq_region = "14";
my $seq_region_start = 5000000;
my $seq_region_end = 5010000;
my $output_format = "clustalw";
my $genome_dbs;

use Bio::EnsEMBL::Registry;
Bio::EnsEMBL::Registry->load_registry_from_db(
    -host => 'ensembldb.ensembl.org',
    -user => 'anonymous',
    -port => 5306);

my $genome_db_adaptor = Bio::EnsEMBL::Registry->get_adaptor(
    'Multi', 'compara', 'GenomeDB');

throw("Cannot connect to Compara") if (!$genome_db_adaptor);

my $this_genome_db = $genome_db_adaptor->fetch_by_name_assembly($species);
my $target_genome_db = $genome_db_adaptor->fetch_by_name_assembly($targetSpecies);
push(@$genome_dbs, $this_genome_db);
push(@$genome_dbs, $target_genome_db);

my $method_link_species_set_adaptor = Bio::EnsEMBL::Registry->get_adaptor(
    'Multi', 'compara', 'MethodLinkSpeciesSet');

my $method_link_species_set =
$method_link_species_set_adaptor->fetch_by_method_link_type_GenomeDBs(
      $alignment_type,
  $genome_dbs);

throw("The database do not contain any $alignment_type data for $target_genome_db!")
    if (!$method_link_species_set);

my $slice_adaptor = Bio::EnsEMBL::Registry->get_adaptor($species, 'core', 'Slice');

throw("Registry configuration file has no data for connecting to <$species>")
if (!$slice_adaptor);

my $query_slice = $slice_adaptor->fetch_by_region('toplevel', $seq_region, $seq_region_start, $seq_region_end);

throw("No Slice can be created with coordinates $seq_region:$seq_region_start-".
    "$seq_region_end") if (!$query_slice);

# Fetching all the GenomicAlignBlock corresponding to this Slice:
my $genomic_align_block_adaptor = Bio::EnsEMBL::Registry->get_adaptor(
    'Multi', 'compara', 'GenomicAlignBlock');

my $genomic_align_blocks =
$genomic_align_block_adaptor->fetch_all_by_MethodLinkSpeciesSet_Slice(
      $method_link_species_set,
  $query_slice);

my $all_aligns;

# Get a Bio::SimpleAlign object from every GenomicAlignBlock
foreach my $this_genomic_align_block (@$genomic_align_blocks) {
    my $simple_align = $this_genomic_align_block->get_SimpleAlign;
    push(@$all_aligns, $simple_align);
}

# print all the genomic alignments using a Bio::AlignIO object
my $alignIO = Bio::AlignIO->newFh(
-interleaved => 0,
-fh => \*STDOUT,
-format => $output_format,
    -idlength => 10
);

foreach my $this_align (@$all_aligns) {
    print $alignIO $this_align;
}

exit;

```

Any pointers on what I'm doing wrong ?

# Answer
While I wouldn't consider this as the best answer or solution - I decided to use Ensembl REST api to extract the coordinates of alignment blocks.

The quick fix - in python looks like this:

```
#This is the main function which uses Ensembl REST API to access the Enesembl Compara database
#And retrive the alignment blocks which correspond to the MM10 standard.
#The function accepts chromosome, start and array and returns the alignment blocks
#reported by the API. Since the blocks can be on different chromosomes/and/or/contigs, it
#returns the whole output in form of dataframe - containing chromosome, start and stop of the
#alignment blocks corresponding to input and all the other alignment blocks corresponding the
#16 mice species, if available.
#!!!!!!!----------Caution start----------!!!!!!!
#The API is a bit buggy, so it is advised not to request for specific species.
#!!!!!!!-----------Caution end-----------!!!!!!!
def liftover(chrom, start, end):
    mmStart, mmEnd, mmchr = [], [], []
    caStart, caEnd, cachr = [], [], []
    #ryStart, ryEnd, rychr = [], [], []
    server = "https://rest.ensembl.org"
    ext = "/alignment/region/mus_musculus/" + str(chrom) + ":" + str(start) + "-" + str(end) + ":1?species_set_group=murinae;method=CACTUS_HAL;content-type=application/json"
    headers={ "Content-Type" : "application/json"}
    r = requests.get(server+ext, headers={ "Content-Type" : "application/json"})
    if(r.ok):
        json_string = json.dumps(r.json())
        jData = json.loads(r.content)
        results = r.json()
        print("The response contains {0} properties".format(len(jData)))
        print("\n")
    else:
        r.raise_for_status()

    for element in jData:
        element.pop('tree', None)

    for element in jData:
        for key, value in element.items():
            for z in value:
                if (z["species"]=="mus_musculus"):
                    mmStart.append(z["start"])
                    mmEnd.append(z["end"])
                    mmchr.append(z["seq_region"])
                if (z["species"]=="mus_musculus_casteij"):
                    caStart.append(z["start"])
                    caEnd.append(z["end"])
                    cachr.append(z["seq_region"])
    if (len(mmStart) != len(caStart)):
        print("Length mismatch between species blocks - proceed with caution!")
    return [min(mmStart), max(mmEnd), min(caStart), max(caEnd), len(mmStart)]

```

Again, it depends on how alignment blocks are reported - and requires further digging to resolve which mm10 alignment block resembles CAST alignment block.
> 2 votes

---
Tags: ensembl, perl, api
---


# Question
Title: How to plot p-values in a circular barplot?
Im trying to plot reactome pathway in for of circular bar plot using ggplot2

my data

```
    dput(df2)
structure(list(ID = 1:10, Pathway_names = structure(c(5L, 1L,
7L, 6L, 2L, 3L, 9L, 8L, 10L, 4L), .Label = c("Antigen Presentation: Folding assembly and peptide loading of class I MHC",
"Antigen processing-Cross presentation", "Class I MHC mediated antigen processing & presentation",
"Cytokine Signaling in Immune system", "Endosomal/Vacuolar pathway",
"ER-Phagosome pathway", "Immunoregulatory interactions between a Lymphoid and a non-Lymphoid cell",
"Interferon alpha/beta signaling", "Interferon gamma signaling",
"Interferon Signaling"), class = "factor"), Entities_found = structure(c(1L,
2L, 3L, 4L, 4L, 5L, 6L, 7L, 8L, 9L), .Label = c("51", "54", "55",
"57", "75", "77", "96", "117", "176"), class = "factor"), Entities_pValue = c(1.11e-16,
1.11e-16, 0.003908431, 2.54e-12, 2.77e-10, 0.0053068, 1.64e-13,
1.11e-16, 1.11e-16, 3.68e-06)), row.names = c(NA, -10L), class = "data.frame")

```

my code

```
ggplot(df2, aes(x = factor(Pathway_names), y = Entities_found)) + ######################################
  geom_bar(aes(fill = Entities_pValue), stat = "identity", position = "stack") +
  #scale_fill_manual(values = c(blank = "white", dark = "black")) +
  coord_polar() +
  theme_minimal() +
  guides(fill = FALSE)

```

So far im not sure how to include both p value and entities in the plot so that the pathway that are coming up can be labeled accordingly

New figure

# Answer
> 3 votes
The `adjustcolor` function from `grDevices` can be used to generate colours with transparency. What you want is colours based on the pathway and transparency to form a colour gradient for the p-value.

`adjustcolor( "red", alpha.f = 0.2)`

For the ggplot example, the `aes` function has an `alpha` parameter.

```
ggplot(df2, aes(x = factor(Pathway_names), y = Entities_found), alpha = Entities_pValue) + ######################################
geom_bar(aes(fill = Entities_found, alpha = Entities_pValue), stat = "identity", position = "stack") +
coord_polar() +
theme_minimal() +
 guides(fill = FALSE)

```

---
Tags: r, visualization, ggplot2
---


# Question
Title: How can I identify a recessive and dominant gene?
For example, I have two allelic genes. How can I identify a recessive and dominant gene? Are there any databases with this information? Or the answer to my question is to study the concept of "haplosufficiency"?

# Answer
I'd probably look up the genes in OMIM, which provides detailed description of the functions and phenotypes associated with genes, including inheritance mode.

For example, this is what you would get for TP53:
> 5 votes

---
Tags: genome
---


# Question
Title: Why is a PacBio read length larger than the aligned reference region?
I recently had some Iso-Seq sequencing done on my organism catfish on the new Sequel platform and got weird alignments for a size selected 4 Kilobase and up fraction after running the isoseq3 pipeline. When I investigated the exon intron junctions using data from (Pac Bio Iso-Seq, Illumina RNA-seq) against our private reference genome, I get PacBio Iso-Seq reads where the sequence length is greater than the seq length on the reference for my bam alignments:

**Seq length:** 6751

**Seq length on ref:** 5267

When blasting the sequence in NCBI it shows 5% similarity to the publicly-available genome. I'll attach a JBrowse image with the same RNA sequenced with Illumina paired RNA-seq. I also included a image and circled the pacbio contig and Illumina is above.

## Question 1:

Is this a problem with JBrowse?

## Question 2:

Does this mean that some DNA was sequenced with the catfish cDNA?

## Question 3:

Why is this aligned sequence's length longer than the length of the reference genome region to which it is aligned?

Thanks

# Answer
### Q1

Probably not - this is probably just displaying the contents of the input (I presume a `bam` file) exactly as they are encoded.

### Q2

The splicing looks a little weird, but **the exon-mapping pattern in the 5' part of the gene makes it seem like RNA rather than DNA**. I think that the large intron is what makes you think that it is DNA, correct?

* Explanations for weird splicing
    + The weird read is a rare splice variant but a true mRNA.
    + The weird read is not completely matured and has an intron that has not yet been spliced. The IsoSeq library prep protocol picked this up somehow and it ended up being sequenced.

## Q3

* Explanations for why **Seq length: 6751** and **Seq length on ref: 5267**
    + There is probably a large amount of clipping on the 5' or 3' end of the read. This could be polyA or some other transcribed-but-not-yet-completely-matured seq. **Or**, the extra sequence could be sequencing error, PacBio adapters, or a chimeric molecule.

### things to try

How to check what is going on:

I'm guessing that you had to upload a `bam` file to JBrowse to make the IsoSeq track. Use the sequence name that you got from Jbrowse and run the following command to extract the alignment info from the bam file:

```
samtools view /path/to/your/isoseq/file.bam | grep 'weird_read_id'

```

Look at the cigar string to see what it says about 5' or 3' clipping and corroborate with the sequence data.
> 5 votes

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


title: Importing GFF file with Biopython

text:
# Question
Title: Importing GFF file with Biopython
Is there a way to import a GFF file for an organism with Biopython in the same way you can for a Genbank file?

For example,

```
from Bio import Entrez as ez
ez.email = '...'
handle = ez.efetch(db='gene', id='AE015451.2', rettype='genbank', retmode='xml')
h = handle.read()
print(h)

```

This will import a Gebank xml file for Pseudomonas putida. Is there a way to do this for a GFF file?

# Answer
No, there is currently no GFF support in biopython.

However, you can read in GFF files into python using this package, gffutils. There are also a few other packages to read/write GFF files, like `gff3`.
> 10 votes

# Answer
I would use `BCBio` for gff handling as it is written to directly interface with BioPython’s object model.

The only downside is that I believe it is no longer actively supported. It is however the package that the BioPython docs use generally. There are plans to properly incorporate GFF/GTF parsers in to BP in the not too distant future according to that link though.
> 2 votes

# Answer
Most GFF parsers handle the work of reading annotation data into objects for convenient data access, but most do not handle the important task of resolving relationships between features. The **tag** Python library does both, grouping related features together for inspection, traversal, and feature-by-feature processing.

**CAVEAT**: The tag library only supports GFF3.

**DISCLAIMER**: I'm the developer of the tag library, so this is a shameless plug. :-)
> 1 votes

---
Tags: python, phylogenetics, genbank, gff
---


# Question
Title: Transcript Coordinate Ranges to Genomic Coordinates
I have 2 GFF3 files:

1. Features using transcript IDs as the landmarks. i.e. "CDS" feature types using coordinates from transcript space.
2. Features using chromosome IDs as the landmarks. i.e. "exon" feature types using coordinates from chromosome space.

I would like to transform the coordinates of the features in file 1 into the coordinate space of file 2. i.e. transcript-based coordinates into genomic-based coodinates.

Here is an example for file1:

```
cat transcript_orfs.gff3
##gff-version 3
##sequence-region Tx.1 1 4000
Tx.1    ORF_finder  gene    1   4000    .   +   .   ID=1
Tx.1    ORF_finder  mRNA    1   4000    .   +   .   ID=2;Parent=1
Tx.1    ORF_finder  exon    1501    2500    .   +   0   ID=3;Parent=2
Tx.1    ORF_finder  CDS 1501    2500    .   +   0   ID=4;Parent=2

gt sketch -addintrons transcript_orfs.png transcript_orfs.gff3

```

Here is an example for file1:

```
cat genome.gff3
##gff-version 3
##sequence-region chr3A_part1 1 454103970
chr3A_part1 genome_assembly gene    1001    6000    .   +   .   ID=Tx
chr3A_part1 genome_assembly mRNA    1001    6000    .   +   .   ID=Tx.1;Parent=Tx
chr3A_part1 genome_assembly exon    1001    3000    .   +   .   ID=Tx.1.exon1;Parent=Tx.1
chr3A_part1 genome_assembly five_prime_UTR  1001    2000    .   +   .   ID=Tx.1.utr5;Parent=Tx.1
chr3A_part1 genome_assembly CDS 2001    3000    .   +   0   ID=Tx.1.cds1;Parent=Tx.1
chr3A_part1 genome_assembly exon    4001    6000    .   +   .   ID=Tx.1.exon2;Parent=Tx.1
chr3A_part1 genome_assembly CDS 4001    5000    .   +   2   ID=Tx.1.cds2;Parent=Tx.1
chr3A_part1 genome_assembly three_prime_UTR 5001    6000    .   +   .   ID=Tx.1.utr3;Parent=Tx.1

gt sketch -addintrons genome.png genome.gff3

```

I want to convert the coordinates of the features from file1 into chromosome genomic coordinates. I expect to get something like the following output:

```
cat output.gff3
##gff-version 3
##sequence-region chr3A_part1 1 454103970
chr3A_part1 ORF_finder  gene    1001    6000    .   +   .   ID=1
chr3A_part1 ORF_finder  mRNA    1001    6000    .   +   .   ID=2;Parent=1
chr3A_part1 ORF_finder  exon    2501    3000    .   +   0   ID=3.1;Parent=2
chr3A_part1 ORF_finder  CDS 2501    3000    .   +   0   ID=4.1;Parent=2
chr3A_part1 ORF_finder  exon    4001    4500    .   +   0   ID=3.2;Parent=2
chr3A_part1 ORF_finder  CDS 4001    4500    .   +   1   ID=4.2;Parent=2

gt sketch -addintrons output.png output.gff3

```

I have looked into using `mapFromTranscripts()` from Bioconductor's `GenomicRanges` library but I have made little progress trying to decipher the manual.

# Answer
This task has the same "flavor" as many I've done before, but each case is so subtly different that it's impossible to write a generalized tool that will work correctly for all circumstances.

R isn't my wheelhouse, but I was able to throw something together pretty quickly with standard Python. This should get you at least 95% of the way to where you're going. :-)

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

from argparse import ArgumentParser, FileType
from collections import defaultdict
from re import search, sub

cli = ArgumentParser()
cli.add_argument('txspace', type=FileType('r'))
cli.add_argument('chrspace', type=FileType('r'))
args = cli.parse_args()

exons = defaultdict(list)
for line in args.txspace:
    if line.startswith('#') or line.strip() == '':
        continue
    values = line.split('\t')
    featuretype = values[2]
    if featuretype not in ('CDS', 'exon'):
        continue
    txid = values[0]
    exons[txid].append(values)

for line in args.chrspace:
    values = line.split('\t')
    if len(values) != 9:
        print(line, end='')
        continue
    featuretype = values[2]
    if featuretype not in ('gene', 'mRNA'):
        continue
    print(line, end='')
    if featuretype == 'mRNA':
        txid = search('ID=([^;\n]+)', values[8]).group(1)
        start = int(values[3])
        txexons = exons[txid]
        for exon in txexons:
            exon[0] = values[0]
            exon[3] = str(int(exon[3]) + start)
            exon[4] = str(int(exon[4]) + start)
            exon[8] = sub(r'Parent=[^;\n]+', 'Parent=' + txid, exon[8])
            print('\t'.join(exon), end='')

```
> 2 votes

# Answer
The linked gist contains a class I wrote for doing this on GTF files:

https://gist.github.com/IanSudery/d8349c22823a475ceb489c3e8aeb448e

It uses the GTF class from cgat, which can be found here: https://github.com/cgat-developers/cgat-apps

You would use it to do this task like so:

```
from cgat import gtf
from cgat import iotools
import TranscriptCoordsInterconverter

trans_coord_file = gtf.iterator(iotools.open_file(args[1])
transcripts = dict()

for transcript in gtf.transcript_iterator(trans_coords_file):
    transcript_id = transcript[0].transcript_id
    transcripts[transcript_id] = transcript

trans_coord_file.close()

genome_coord_file = gtf.iterator(iotools.open_file(args[2]))
for transcript in gtf.transcript_iterator(genome_coord_file):

    transcript_id = transcript[0].transcript_id
    contig = transcript[0].contig
    strand = transcript[0].strand

    if transcript_id not in transcripts:
        continue

    converter = TranscriptCoordsInterconverter(transcript)
    for exon in transcripts[transcript_id]:
        genome_intervals = converter.transcript_interval2genome_intervals((exon.start,exon.end))
        for interval in genome_intervals:
            new_entry = GTF.Entry().fromGTF(exon)
            new_entry.contig = contig
            new_entry.strand = strand
            new_entry.start = interval[0]
            new_entry.end = interval[1]
            print (str(new_entry)+"\n")

```
> 1 votes

---
Tags: bioconductor, sequence-annotation, gff3, interval
---


# Question
Title: MrBayes runs as single-core, but hangs when multi-core
I compiled MrBayes v3.2.7 from github, with the configuration options `./configure --without-beagle --enable-mpi`. The program is called mb.

I tried one of the included example files "sceloporus.nex". I find that if I run `mb sceloporus.nex`, then it appears to make progress, and writes files to disk.

I used MPI with 1 core, and it behaves as before, making progress, writing files to disk.

But if I run with more than 1 core, using MPI, with `mpirun -np 2 mb sceloporus.nex`, then it only makes it through some of the analysis setup, doesn't write any files to disk, and just hangs. If I look at `top`, I see activity on the number of cores I specified. I finally had to manually terminate the program.

Is there something I should be doing differently to have MrBayes properly run in parallel?

Here is a copy of my terminal from these 3 attempts:

```
[ptitle@iu mbTest]$ ls
sceloporus.nex
[ptitle@iu mbTest]$ ~/MrBayes/src/mb sceloporus.nex

                            MrBayes pre-3.2.7 x86_64

                      (Bayesian Analysis of Phylogeny)

                             (Parallel version)
                         (1 processors available)

              Distributed under the GNU General Public License

               Type "help" or "help <command>" for information
                     on the commands that are available.

                   Type "about" for authorship and general
                       information about the program.

   Executing file "sceloporus.nex"
   DOS line termination
   Longest line length = 1632
   Parsing file
   Expecting NEXUS formatted file
   Reading data block
      Allocated taxon set
      Allocated matrix
      Defining new matrix with 123 taxa and 1606 characters
      Data is Dna
      Missing data coded as ?
      Gaps coded as -
      Taxon   1 -> AZYuJAS289
      Taxon   2 -> AZYuJAS290
      Taxon   3 -> AZlaJAS294
      Taxon   4 -> AZcoTBP271
      Taxon   5 -> AZcoGM1054
      Taxon   6 -> AZcoDGM994
      Taxon   7 -> AZgiP26438
      Taxon   8 -> AZgrGM1012
      Taxon   9 -> AZgrGM1013
      Taxon  10 -> AZgrGM1024
      Taxon  11 -> AZlaDGM856
      Taxon  12 -> AZmaDGM984
      Taxon  13 -> AZmaDGM992
      Taxon  14 -> AZmaRM4401
      Taxon  15 -> AZmaRM4402
      Taxon  16 -> AZmaRM4403
      Taxon  17 -> AZmoGM1045
      Taxon  18 -> AZmoDGM492
      Taxon  19 -> AZmoDGM699
      Taxon  20 -> AZmoDGM704
      Taxon  21 -> AZmoRM4394
      Taxon  22 -> AZpiKWS238
      Taxon  23 -> AZyuDGM826
      Taxon  24 -> AZyuDGM838
      Taxon  25 -> AZyuDGM852
      Taxon  26 -> AZyuDGM858
      Taxon  27 -> AZyuDGM865
      Taxon  28 -> CAfrP25750
      Taxon  29 -> CAimDGM534
      Taxon  30 -> CAimDGM558
      Taxon  31 -> CAimDGM694
      Taxon  32 -> CAimKWS170
      Taxon  33 -> CAinCRF530
      Taxon  34 -> CAinoDGM74
      Taxon  35 -> CAinCRF529
      Taxon  36 -> CAkeDS2SP7
      Taxon  37 -> CAkeP26556
      Taxon  38 -> CAla78323
      Taxon  39 -> CAlaDAW220
      Taxon  40 -> CAlaM14610
      Taxon  41 -> CAlaM23289
      Taxon  42 -> CAmoJQR132
      Taxon  43 -> CArvJOS138
      Taxon  44 -> CArvESA441
      Taxon  45 -> CAsaBUR167
      Taxon  46 -> CAsaMCFMAG
      Taxon  47 -> CAsarnMCC
      Taxon  48 -> CAsa150079
      Taxon  49 -> CAsa150088
      Taxon  50 -> CAsa179886
      Taxon  51 -> CAsaM19838
      Taxon  52 -> CAsdDGM691
      Taxon  53 -> CAsdDGM748
      Taxon  54 -> CAsdSWT701
      Taxon  55 -> MXBCDGM474
      Taxon  56 -> MXBCDGM475
      Taxon  57 -> MXBCDGM486
      Taxon  58 -> MXchRM4860
      Taxon  59 -> MXduA83134
      Taxon  60 -> MXsoM14488
      Taxon  61 -> MXsoP26449
      Taxon  62 -> MXso210345
      Taxon  63 -> NMdoTWR230
      Taxon  64 -> NMedDGM938
      Taxon  65 -> NMhiU48819
      Taxon  66 -> NMluGM947
      Taxon  67 -> NMluDGM948
      Taxon  68 -> NMsiDGM891
      Taxon  69 -> NMsoDGM900
      Taxon  70 -> NMsoDGM904
      Taxon  71 -> NMsoDGM905
      Taxon  72 -> NMvaDGM924
      Taxon  73 -> NVchDGM612
      Taxon  74 -> NVch162077
      Taxon  75 -> NVclDGM296
      Taxon  76 -> NVclDGM678
      Taxon  77 -> NVesDGM210
      Taxon  78 -> NVlyDGM610
      Taxon  79 -> NVwaDGM608
      Taxon  80 -> UTemGM1058
      Taxon  81 -> UTgrDGM156
      Taxon  82 -> UTgrDGM169
      Taxon  83 -> UTgrDGM175
      Taxon  84 -> UTkaGM1034
      Taxon  85 -> UTsnDM142
      Taxon  86 -> UTsnDG043
      Taxon  87 -> UTsnnDGM645
      Taxon  88 -> UTwsGM261
      Taxon  89 -> UTwsGM263
      Taxon  90 -> UTwsnDGM751
      Taxon  91 -> UTwsJRM4437
      Taxon  92 -> hunsiM376
      Taxon  93 -> hunsiM377
      Taxon  94 -> hunsiRWM619
      Taxon  95 -> lickGM378
      Taxon  96 -> lineatulus
      Taxon  97 -> orcuOM14377
      Taxon  98 -> orcuOM14453
      Taxon  99 -> orcuOM14673
      Taxon 100 -> orcuRWM798
      Taxon 101 -> orcuS201108
      Taxon 102 -> orcuS201124
      Taxon 103 -> zostZ16290
      Taxon 104 -> zosV161292
      Taxon 105 -> zosV161293
      Taxon 106 -> zosV161311
      Taxon 107 -> zosWM1890
      Taxon 108 -> zosGM364
      Taxon 109 -> zosGM365
      Taxon 110 -> zosGM387
      Taxon 111 -> zosGM393
      Taxon 112 -> zosRO332
      Taxon 113 -> zosOM37006
      Taxon 114 -> zosWM1601
      Taxon 115 -> clarkii
      Taxon 116 -> edwardti80
      Taxon 117 -> smaragdnus
      Taxon 118 -> graciosus
      Taxon 119 -> grammicus
      Taxon 120 -> maculosus
      Taxon 121 -> bicanthlis
      Taxon 122 -> torquatus
      Taxon 123 -> variabilis
      Successfully read matrix
      Setting default partition (does not divide up characters)
      Setting model defaults
      Seed (for generating default start values) = 1544128156
      Setting output file names to "sceloporus.nex.run<i>.<p|t>"
   Exiting data block
   Reading mrbayes block
      Excluding character(s)
      Setting Nst to 6
      Setting Rates to Gamma
      Successfully set likelihood model parameters
      Setting Brlenspr to Unconstrained:GammaDir(1.00,1.0000,0.70,0.50)
      Successfully set prior model parameters
      Setting number of generations to 5000000
      Running Markov chain
      MCMC stamp = 9284864458
      Seed = 753696747
      Swapseed = 1544128156
      Model settings:

         Data not partitioned --
            Datatype  = DNA
            Nucmodel  = 4by4
            Nst       = 6
                        Substitution rates, expressed as proportions
                        of the rate sum, have a Dirichlet prior
                        (1.00,1.00,1.00,1.00,1.00,1.00)
            Covarion  = No
            # States  = 4
                        State frequencies have a Dirichlet prior
                        (1.00,1.00,1.00,1.00)
            Rates     = Gamma
                        The distribution is approximated using 4 categories.
                        Shape parameter is exponentially
                        distributed with parameter (1.00).

      Active parameters:

         Parameters
         ---------------------
         Revmat              1
         Statefreq           2
         Shape               3
         Ratemultiplier      4
         Topology            5
         Brlens              6
         ---------------------

         1 --  Parameter  = Revmat
               Type       = Rates of reversible rate matrix
               Prior      = Dirichlet(1.00,1.00,1.00,1.00,1.00,1.00)

         2 --  Parameter  = Pi
               Type       = Stationary state frequencies
               Prior      = Dirichlet

         3 --  Parameter  = Alpha
               Type       = Shape of scaled gamma distribution of site rates
               Prior      = Exponential(1.00)

         4 --  Parameter  = Ratemultiplier
               Type       = Partition-specific rate multiplier
               Prior      = Fixed(1.0)

         5 --  Parameter  = Tau
               Type       = Topology
               Prior      = All topologies equally probable a priori
               Subparam.  = V

         6 --  Parameter  = V
               Type       = Branch lengths
               Prior      = Unconstrained:GammaDir(1.0,1.0000,0.7,0.5)

      Number of chains per MPI processor = 8

      The MCMC sampler will use the following moves:
         With prob.  Chain will use move
            0.94 %   Dirichlet(Revmat)
            0.94 %   Slider(Revmat)
            0.94 %   Dirichlet(Pi)
            0.94 %   Slider(Pi)
            1.89 %   Multiplier(Alpha)
            9.43 %   ExtSPR(Tau,V)
            9.43 %   ExtTBR(Tau,V)
            9.43 %   NNI(Tau,V)
            9.43 %   ParsSPR(Tau,V)
           37.74 %   Multiplier(V)
           13.21 %   Nodeslider(V)
            5.66 %   TLMultiplier(V)

      Division 1 has 624 unique site patterns
      Initializing conditional likelihoods
      Using standard FMA likelihood calculator for division 1 (single-precision)

      Initial log likelihoods and log prior probs for run 1:
         Chain 1 -- -38239.042144 -- 53.554896
         Chain 2 -- -37743.352685 -- 53.554896
         Chain 3 -- -38049.750737 -- 53.554896
         Chain 4 -- -38759.582835 -- 53.554896

      Initial log likelihoods and log prior probs for run 2:
         Chain 1 -- -38348.591921 -- 53.554896
         Chain 2 -- -39602.594468 -- 53.554896
         Chain 3 -- -38874.065643 -- 53.554896
         Chain 4 -- -39227.723762 -- 53.554896

      Using a relative burnin of 25.0 % for diagnostics

      Chain results (5000000 generations requested):

         0 -- [-38239.042] (-37743.353) (-38049.751) (-38759.583) * [-38348.592] (-39602.594) (-38874.066) (-39227.724) (...0 remote chains...)
^C
[ptitle@iu mbTest]$ ls
sceloporus.nex  sceloporus.nex.mcmc  sceloporus.nex.run1.p  sceloporus.nex.run1.t  sceloporus.nex.run2.p  sceloporus.nex.run2.t
[ptitle@iu mbTest]$ mpirun -np 1 ~/MrBayes/src/mb sceloporus.nex

                            MrBayes pre-3.2.7 x86_64

                      (Bayesian Analysis of Phylogeny)

                             (Parallel version)
                         (1 processors available)

              Distributed under the GNU General Public License

               Type "help" or "help <command>" for information
                     on the commands that are available.

                   Type "about" for authorship and general
                       information about the program.

   Executing file "sceloporus.nex"
   DOS line termination
   Longest line length = 1632
   Parsing file
   Expecting NEXUS formatted file
   Reading data block
      Allocated taxon set
      Allocated matrix
      Defining new matrix with 123 taxa and 1606 characters
      Data is Dna
      Missing data coded as ?
      Gaps coded as -
      Taxon   1 -> AZYuJAS289
      Taxon   2 -> AZYuJAS290
      Taxon   3 -> AZlaJAS294
      Taxon   4 -> AZcoTBP271
      Taxon   5 -> AZcoGM1054
      Taxon   6 -> AZcoDGM994
      Taxon   7 -> AZgiP26438
      Taxon   8 -> AZgrGM1012
      Taxon   9 -> AZgrGM1013
      Taxon  10 -> AZgrGM1024
      Taxon  11 -> AZlaDGM856
      Taxon  12 -> AZmaDGM984
      Taxon  13 -> AZmaDGM992
      Taxon  14 -> AZmaRM4401
      Taxon  15 -> AZmaRM4402
      Taxon  16 -> AZmaRM4403
      Taxon  17 -> AZmoGM1045
      Taxon  18 -> AZmoDGM492
      Taxon  19 -> AZmoDGM699
      Taxon  20 -> AZmoDGM704
      Taxon  21 -> AZmoRM4394
      Taxon  22 -> AZpiKWS238
      Taxon  23 -> AZyuDGM826
      Taxon  24 -> AZyuDGM838
      Taxon  25 -> AZyuDGM852
      Taxon  26 -> AZyuDGM858
      Taxon  27 -> AZyuDGM865
      Taxon  28 -> CAfrP25750
      Taxon  29 -> CAimDGM534
      Taxon  30 -> CAimDGM558
      Taxon  31 -> CAimDGM694
      Taxon  32 -> CAimKWS170
      Taxon  33 -> CAinCRF530
      Taxon  34 -> CAinoDGM74
      Taxon  35 -> CAinCRF529
      Taxon  36 -> CAkeDS2SP7
      Taxon  37 -> CAkeP26556
      Taxon  38 -> CAla78323
      Taxon  39 -> CAlaDAW220
      Taxon  40 -> CAlaM14610
      Taxon  41 -> CAlaM23289
      Taxon  42 -> CAmoJQR132
      Taxon  43 -> CArvJOS138
      Taxon  44 -> CArvESA441
      Taxon  45 -> CAsaBUR167
      Taxon  46 -> CAsaMCFMAG
      Taxon  47 -> CAsarnMCC
      Taxon  48 -> CAsa150079
      Taxon  49 -> CAsa150088
      Taxon  50 -> CAsa179886
      Taxon  51 -> CAsaM19838
      Taxon  52 -> CAsdDGM691
      Taxon  53 -> CAsdDGM748
      Taxon  54 -> CAsdSWT701
      Taxon  55 -> MXBCDGM474
      Taxon  56 -> MXBCDGM475
      Taxon  57 -> MXBCDGM486
      Taxon  58 -> MXchRM4860
      Taxon  59 -> MXduA83134
      Taxon  60 -> MXsoM14488
      Taxon  61 -> MXsoP26449
      Taxon  62 -> MXso210345
      Taxon  63 -> NMdoTWR230
      Taxon  64 -> NMedDGM938
      Taxon  65 -> NMhiU48819
      Taxon  66 -> NMluGM947
      Taxon  67 -> NMluDGM948
      Taxon  68 -> NMsiDGM891
      Taxon  69 -> NMsoDGM900
      Taxon  70 -> NMsoDGM904
      Taxon  71 -> NMsoDGM905
      Taxon  72 -> NMvaDGM924
      Taxon  73 -> NVchDGM612
      Taxon  74 -> NVch162077
      Taxon  75 -> NVclDGM296
      Taxon  76 -> NVclDGM678
      Taxon  77 -> NVesDGM210
      Taxon  78 -> NVlyDGM610
      Taxon  79 -> NVwaDGM608
      Taxon  80 -> UTemGM1058
      Taxon  81 -> UTgrDGM156
      Taxon  82 -> UTgrDGM169
      Taxon  83 -> UTgrDGM175
      Taxon  84 -> UTkaGM1034
      Taxon  85 -> UTsnDM142
      Taxon  86 -> UTsnDG043
      Taxon  87 -> UTsnnDGM645
      Taxon  88 -> UTwsGM261
      Taxon  89 -> UTwsGM263
      Taxon  90 -> UTwsnDGM751
      Taxon  91 -> UTwsJRM4437
      Taxon  92 -> hunsiM376
      Taxon  93 -> hunsiM377
      Taxon  94 -> hunsiRWM619
      Taxon  95 -> lickGM378
      Taxon  96 -> lineatulus
      Taxon  97 -> orcuOM14377
      Taxon  98 -> orcuOM14453
      Taxon  99 -> orcuOM14673
      Taxon 100 -> orcuRWM798
      Taxon 101 -> orcuS201108
      Taxon 102 -> orcuS201124
      Taxon 103 -> zostZ16290
      Taxon 104 -> zosV161292
      Taxon 105 -> zosV161293
      Taxon 106 -> zosV161311
      Taxon 107 -> zosWM1890
      Taxon 108 -> zosGM364
      Taxon 109 -> zosGM365
      Taxon 110 -> zosGM387
      Taxon 111 -> zosGM393
      Taxon 112 -> zosRO332
      Taxon 113 -> zosOM37006
      Taxon 114 -> zosWM1601
      Taxon 115 -> clarkii
      Taxon 116 -> edwardti80
      Taxon 117 -> smaragdnus
      Taxon 118 -> graciosus
      Taxon 119 -> grammicus
      Taxon 120 -> maculosus
      Taxon 121 -> bicanthlis
      Taxon 122 -> torquatus
      Taxon 123 -> variabilis
      Successfully read matrix
      Setting default partition (does not divide up characters)
      Setting model defaults
      Seed (for generating default start values) = 1544128192
      Setting output file names to "sceloporus.nex.run<i>.<p|t>"
   Exiting data block
   Reading mrbayes block
      Excluding character(s)
      Setting Nst to 6
      Setting Rates to Gamma
      Successfully set likelihood model parameters
      Setting Brlenspr to Unconstrained:GammaDir(1.00,1.0000,0.70,0.50)
      Successfully set prior model parameters
      Setting number of generations to 5000000
      Running Markov chain
      MCMC stamp = 9909051894
      Seed = 1261289135
      Swapseed = 1544128192
      Model settings:

         Data not partitioned --
            Datatype  = DNA
            Nucmodel  = 4by4
            Nst       = 6
                        Substitution rates, expressed as proportions
                        of the rate sum, have a Dirichlet prior
                        (1.00,1.00,1.00,1.00,1.00,1.00)
            Covarion  = No
            # States  = 4
                        State frequencies have a Dirichlet prior
                        (1.00,1.00,1.00,1.00)
            Rates     = Gamma
                        The distribution is approximated using 4 categories.
                        Shape parameter is exponentially
                        distributed with parameter (1.00).

      Active parameters:

         Parameters
         ---------------------
         Revmat              1
         Statefreq           2
         Shape               3
         Ratemultiplier      4
         Topology            5
         Brlens              6
         ---------------------

         1 --  Parameter  = Revmat
               Type       = Rates of reversible rate matrix
               Prior      = Dirichlet(1.00,1.00,1.00,1.00,1.00,1.00)

         2 --  Parameter  = Pi
               Type       = Stationary state frequencies
               Prior      = Dirichlet

         3 --  Parameter  = Alpha
               Type       = Shape of scaled gamma distribution of site rates
               Prior      = Exponential(1.00)

         4 --  Parameter  = Ratemultiplier
               Type       = Partition-specific rate multiplier
               Prior      = Fixed(1.0)

         5 --  Parameter  = Tau
               Type       = Topology
               Prior      = All topologies equally probable a priori
               Subparam.  = V

         6 --  Parameter  = V
               Type       = Branch lengths
               Prior      = Unconstrained:GammaDir(1.0,1.0000,0.7,0.5)

      Number of chains per MPI processor = 8

      The MCMC sampler will use the following moves:
         With prob.  Chain will use move
            0.94 %   Dirichlet(Revmat)
            0.94 %   Slider(Revmat)
            0.94 %   Dirichlet(Pi)
            0.94 %   Slider(Pi)
            1.89 %   Multiplier(Alpha)
            9.43 %   ExtSPR(Tau,V)
            9.43 %   ExtTBR(Tau,V)
            9.43 %   NNI(Tau,V)
            9.43 %   ParsSPR(Tau,V)
           37.74 %   Multiplier(V)
           13.21 %   Nodeslider(V)
            5.66 %   TLMultiplier(V)

      Division 1 has 624 unique site patterns
      Initializing conditional likelihoods
      Using standard FMA likelihood calculator for division 1 (single-precision)

      Initial log likelihoods and log prior probs for run 1:
         Chain 1 -- -39098.978535 -- 53.554896
         Chain 2 -- -38054.044640 -- 53.554896
         Chain 3 -- -39425.015629 -- 53.554896
         Chain 4 -- -40157.756230 -- 53.554896

      Initial log likelihoods and log prior probs for run 2:
         Chain 1 -- -40084.553712 -- 53.554896
         Chain 2 -- -39002.725314 -- 53.554896
         Chain 3 -- -38719.752582 -- 53.554896
         Chain 4 -- -39156.129424 -- 53.554896

      Using a relative burnin of 25.0 % for diagnostics

      Chain results (5000000 generations requested):

         0 -- [-39098.979] (-38054.045) (-39425.016) (-40157.756) * [-40084.554] (-39002.725) (-38719.753) (-39156.129) (...0 remote chains...)
^C[ptitle@iu mbTest]$ mpirun -np 2 ~/MrBayes/src/mb sceloporus.nex

                            MrBayes pre-3.2.7 x86_64

                      (Bayesian Analysis of Phylogeny)

                             (Parallel version)
                         (2 processors available)

              Distributed under the GNU General Public License

               Type "help" or "help <command>" for information
                     on the commands that are available.

                   Type "about" for authorship and general
                       information about the program.

   Executing file "sceloporus.nex"
   DOS line termination
   Longest line length = 1632
   Parsing file
   Expecting NEXUS formatted file
   Reading data block
      Allocated taxon set
      Allocated matrix
      Defining new matrix with 123 taxa and 1606 characters
      Data is Dna
      Missing data coded as ?
      Gaps coded as -
      Taxon   1 -> AZYuJAS289
      Taxon   2 -> AZYuJAS290
      Taxon   3 -> AZlaJAS294
      Taxon   4 -> AZcoTBP271
      Taxon   5 -> AZcoGM1054
      Taxon   6 -> AZcoDGM994
      Taxon   7 -> AZgiP26438
      Taxon   8 -> AZgrGM1012
      Taxon   9 -> AZgrGM1013
      Taxon  10 -> AZgrGM1024
^C[ptitle@iu mbTest]$

```

# Answer
> 1 votes
It is not trivial getting MB to run in parallel because from my experience it required granting privileges to regular users. I used this via MPI (MPICH) on a Slurm/SGE cluster.

Essentially, this is coarse grain parallel processing.

It is certainly worth parallelizing MB because each chain can be run on a different nodes and avoids RAM bottlenecks on complex models and big data sets. For example, 8 to 10+ chains can be run, ensuring convergence quickly by running chains across nodes.

---

**Update** I later found it was much easier to run via EasyBuild with OpenMPI using the command `runmpi`.

---
Tags: phylogenetics, phylogeny, molecular-genetics, mrbayes
---