# Question
Title: Missing ',' in line when Biopython reads a nexus tree
I want to edit a tree that I got from BEAST2 treeannotator in nexus-format.

Usually I use the module Phylo from Biopython for such work but `Phylo.read(r"filename.tree", "nexus")` gave me the next exception:

```
---------------------------------------------------------------------------
NexusError                                Traceback (most recent call last)
Input In [29], in <cell line: 1>()
----> 1 Phylo.read(r"filename.tree", "nexus")

File ~\miniconda3\lib\site-packages\Bio\Phylo\_io.py:60, in read(file, format, **kwargs)
     58 try:
     59     tree_gen = parse(file, format, **kwargs)
---> 60     tree = next(tree_gen)
     61 except StopIteration:
     62     raise ValueError("There are no trees in this file.") from None

File ~\miniconda3\lib\site-packages\Bio\Phylo\_io.py:49, in parse(file, format, **kwargs)
     34 """Parse a file iteratively, and yield each of the trees it contains.
     35
     36 If a file only contains one tree, this still returns an iterable object that
   (...)
     46
     47 """
     48 with File.as_handle(file) as fp:
---> 49     yield from getattr(supported_formats[format], "parse")(fp, **kwargs)

File ~\miniconda3\lib\site-packages\Bio\Phylo\NexusIO.py:40, in parse(handle)
     32 def parse(handle):
     33     """Parse the trees in a Nexus file.
     34
     35     Uses the old Nexus.Trees parser to extract the trees, converts them back to
   (...)
     38     eventually change Nexus to use the new NewickIO parser directly.)
     39     """
---> 40     nex = Nexus.Nexus(handle)
     42     # NB: Once Nexus.Trees is modified to use Tree.Newick objects, do this:
     43     # return iter(nex.trees)
     44     # Until then, convert the Nexus.Trees.Tree object hierarchy:
     45     def node2clade(nxtree, node):

File ~\miniconda3\lib\site-packages\Bio\Nexus\Nexus.py:668, in Nexus.__init__(self, input)
    665 self.options["gapmode"] = "missing"
    667 if input:
--> 668     self.read(input)
    669 else:
    670     self.read(DEFAULTNEXUS)

File ~\miniconda3\lib\site-packages\Bio\Nexus\Nexus.py:718, in Nexus.read(self, input)
    716     break
    717 if title in KNOWN_NEXUS_BLOCKS:
--> 718     self._parse_nexus_block(title, contents)
    719 else:
    720     self._unknown_nexus_block(title, contents)

File ~\miniconda3\lib\site-packages\Bio\Nexus\Nexus.py:759, in Nexus._parse_nexus_block(self, title, contents)
    757 for line in block.commandlines:
    758     try:
--> 759         getattr(self, "_" + line.command)(line.options)
    760     except AttributeError:
    761         raise NexusError("Unknown command: %s " % line.command) from None

File ~\miniconda3\lib\site-packages\Bio\Nexus\Nexus.py:1144, in Nexus._translate(self, options)
   1142         break
   1143     elif c != ",":
-> 1144         raise NexusError("Missing ',' in line %s." % options)
   1145 except NexusError:
   1146     raise

NexusError: Missing ',' in line 1 AB298157.1_2015_-7.9133750332192605_114.8086828279248, 2 AB298158.1_2007_-8.41698974207…

```

Using `Nexus.read(Nexus(), input=r"filename.tree")` gave the same result. Please could anyone help with this? I cannot understand the reason of this error because nexus file looks correct.

---

**Nexus format**

A NEXUS tree file is a rigid but very useful format which can be used to store a phylogenetic tree, along with associated metadata, for example, the alignment.

It will always comprise a:

> TAXA block

The TAXA block contains information about taxa. ... can comprise a,

> DATA block

The DATA block contains the data matrix (e.g. sequence alignment). ... in this context importantly will have a,

> TREES block

The TREES block contains trees described via the format, e.g. ((A,B),C);:

Each block will open with a BEGIN statement and terminate with an END statement. Again its versatility to store specific information within a block is valuable property of the format.

---

**NB!** I encountered some problems with my following work. So I would like to warn readers about them. It is allowed in the **newick** format to surround the label with *quotes*, & some programmes or scripts add them to those names that have ambiguous characters. But it can lead to exceptions during the following phylogenetic analysis, for instance, in **BEAST**. I wish you would be careful with this.

# Answer
> 1 votes
That's a complicated output.

Two solutions, the easiest is simply to use Andrew Rambaut's `figtree` http://tree.bio.ed.ac.uk/software/figtree/ . Figtree is a very, very good tree drawing program and it will work; it works on Mac, Windows and Linux.

You can then export the file from Figure - from memory in numerous formats, e.g. phylip format, and import it into Biopython if that is where it is needed. If you are parsing Beast output, Beast has A LOT of internal tools and you might be attempting to replace those via Biopython. I would avoid it. Also there is the package for processing Beast output (which I forget), I think its called Tracer.

With Biopython try parse, it might help if there are multiple trees:

```
from Bio import Phylo
trees = Phylo.parse("examples.nex", "nexus")
for tree in trees:
    print(tree.name)

```

I would avoid the 'r'. I agree it makes no difference, but really you'd need to be using unusual file names for 'r' to be needed. Its possible you are using it to get around \s, i.e. spaces in the file name. You could just think about simply not including \ and space type characters in file names. In regex it is used for certain but that's because its 'inescapable' (no pun intended, but I thought it was funny).

---

Just to wrap up and clarify, I recommended:

* To open the tree in Figtree;
* Find the export menu and output as either 'newick' or preferably 'nexus' or both;
* Re-run Biopython importing as 'nexus' or 'newick' (if nexus fails).

The OP reported ...

> I have just compared them (nexus-\>FigTree-\>nexus).

Good so it opened up in FigTree and exported, this was expected.

> Fig Tree compose all this actually in newick & stash it in nexus block decor. While in the previous case newick tree contains IDs, & there is a seclude block with ID-name pairs

I missed this bit, "newick stashed in a nexus block". That is potentially really cool because I suspect that is within the TREES block. If that's correct then that is nexus format. This should now open Biopython.

**Why the bug?**

One of the issues with Nexus is there are two variants regarding how the TAXA block relates to the TREES block. In one format they recode the taxa - I definitely understand why. Recoding is possibly (very likely) be being used in Beast to increase algorithm efficiency - it makes no sense slushing a load of taxa names around a convergence algorithm. Writing all that stuff out is not cool either.

However, when you are looking at a tree that's not what you want at all - efficiency is an irrelevance and readability is important. Thus recoding the taxa labels is definitely not wanted. This is a possible explanation for the read error.

Whatever the reason for the bug I am certain Biopython has a work around because in phylip alignment format this is certainly true. There is a custom module that handles the specific issues the format generates.

It was fun but time to move on.

---
Tags: phylogenetics, python, phylogeny, beast
---


# Question
Title: Constructing Two Triangle Heatmaps in One Square
I'm looking to make a two, two range colour scheme heatmap where each range corresponds to a separate triangle in the heatmap square.

I am looking for a relevant guide for this heatmap. Any recommended package/library or guides? Thanks!

# Answer
> 1 votes
This solution doesn't reproduce the plot, but it would do if the numerical ranges between the triangles do not overlap.

```
import matplotlib.pyplot as plt
import seaborn as sns

xx = plt.get_cmap('rocket').copy()
xx.set_under('none')
yy = plt.get_cmap('mako').copy()
yy.set_over('none')
axis = sns.heatmap(dataframe, vmin=-1, cmap=xx, cbar_kws={'pad': -0.04})
sns.heatmap(dataframe, vmin=0, vmax=1, cmap=yy, ax=axis)
plt.show()

```

My colour schemes are not great but they are distinct, you can search for more here: https://seaborn.pydata.org/tutorial/color\_palettes.html

The idea is this would address a situation where you compare the *r*-value OR a *p*-value against another set of values where are \> 1. That's what the `vmin` and `vmax` are doing.

This is a ONE np array or dataframe solution. In other words you will need to combine your data sets. The solution you want is a two dataframes in one heatmap square. You can put two triangles side-by-side, i.e. one dataframe per heatmap, but I don't know who to place two separate dataframes in one heatmap.

---
Tags: r, python, heatmap, complexheatmap
---


# Question
Title: How can I draw a diagram of hydrogen bonds only on the basis of C-alpha backbone?
We are working on an algorithm to assign proteins only on the basis of alpha-carbons.

My teacher told me to draw a diagram to represent hydrogen bonds only on the basis of alpha-carbons.

I showed him the following picture:

He said this picture doesn't count, as this diagram shows all the atoms.

I couldn't figure out how to do this without showing oxygen, non-alpha-carbon, and nitrogen atoms.

How can I draw a diagram of hydrogen bonds only on the basis of C-alpha backbone?

# Answer
There are three options:

1. it's a badly phrased Ramachandran plot question
2. it's a question aimed at the describing the properties of different SS
3. it's a trick question

φ & ψ can be used to classify SS (cf. Ramachandran plot), these are however proper dihedrals, so require 4 atoms, which you lack.

The sequential Cαs in a strand are increasingly distant by 3–4Å, while in a helix there's a periodicity every 3–4 residues. The images you were given are likely to remind you that between two strands in a sheet the Cαs are ~5Å away from the other strand, but the order is different —the indices will increase together in parallel strands, but not in anti-parallel strands.

A C-alpha trace has the issue that the residues are often not known (`UNK` is technically represented as homoalanine, but is an unknown) and the resolution is \> 4Å. Any SS allocation will be highly inaccurate making it a trick question.
> 2 votes

---
Tags: protein-structure, homework
---


# Question
Title: Why does GATK produce both 0/1 and 1/0 genotypes in the same file? Are the two not equivalent?
I have always thought that `1/0` and `0/1` in VCF genotype fields are equivalent. And yet, GATK uses both. For example, these are two variants called in the same sample and the same run of GATK 4.1.4.0:

```
chr7    117120317   .   ATTCATTGTTTTGAAAGAAAGATGGAAGAATGAACTGAAG    A   748.97  .   AC=1;AF=0.5;AN=2;DP=64;ExcessHet=3.0103;FS=0;MLEAC=1;MLEAF=0.5;MQ=60;QD=11.89;SOR=7.223 GT:AD:DP:GQ:PL:SB   1/0:0,36:63:99:2294,1042,933:0,0,0,36
chr7    117120306   .   GACGCGCT    G   516.6   .   AC=1;AF=0.5;AN=2;BaseQRankSum=0.826;DP=70;ExcessHet=3.0103;FS=0;MLEAC=1;MLEAF=0.5;MQ=60;MQRankSum=0;QD=8.2;ReadPosRankSum=-5.786;SOR=0.976  GT:AD:DP:GQ:PL:SB   0/1:37,26:63:99:524,0,1203:0,20,0,26

```

And here are the same lines redacted for clarity:

```
chr7    117120317   .   AT..AG  A   [...]   GT:AD:DP:GQ:PL:SB   1/0:0,36:63:99:2294,1042,933:0,0,0,36
chr7    117120306   .   GA..CT  G   [...]   GT:AD:DP:GQ:PL:SB   0/1:37,26:63:99:524,0,1203:0,20,0,26

```

As you can see, one is genotyped as `1/0` and the other as `0/1`. The VCF specs say:

> GT (String): Genotype, encoded as allele values separated by either of / or |. The allele values are 0 for the reference allele (what is in the REF field), 1 for the first allele listed in ALT, 2 for the second allele list in ALT and so on. For diploid calls examples could be 0/1, 1 | 0, or 1/2, etc. Haploid calls, e.g. on Y, male non-pseudoautosomal X, or mitochondrion, are indicated by having only one allele value. A triploid call might look like 0/0/1. If a call cannot be made for a sample at a given locus, ‘.’ must be specified for each missing allele in the GT field (for example ‘./.’ for a diploid genotype and ‘.’ for haploid genotype). The meanings of the separators are as follows (see the PS field below for more details on incorporating phasing information into the genotypes):
>
> * / : genotype unphased
> * | : genotype phased

This suggests that although `1|0` is not the same as `0|1`, the *unphased* `0/1` and `1/0` genotypes are completely equivalent. Which makes perfect sense given that the `N/N` genotypes are unphased, so `n/N` must be equivalent to `N/n`.

So, my questions are: why does gatk produce some lines with `0/1` and others with `1/0`? Is it a bug? Am I missing something? I can't imagine a developer writing code to pick either `0/1` or `1/0` randomly, so presumably they indicate some difference in what gatk sees about a variant. So are `0/1` and `1/0` not 100% equivalent somehow?

---

<sub> I have seen How is the GT field in a VCF file defined? which suggests the two are equivalent, but this seems to be contradicted by GATK's own output which is why I am asking this again to get a clear and preferably referenced answer.</sub>

# Answer
> 2 votes
**Convention** The convention is it should be in ascending ordering, i.e. 0/1. The question is, why the convention broken?

**Phasing** I think, at a guess, the answer is it depends what reference genome you are using. I would suggest that this result would not occur if the reference genome is mum (or dad) and the alleles are occurring in the child. The output supports this idea because it is marked with '/', i.e. they are not phased, thus the reference genome is neither mum nor dad.

**Example provided** The thing I would like to point out is in the example given `117120306` and `117120317` are not randomly scattered across the genome but appear to be 11 bp apart. For a human genome those are very strongly linked.

**Linkage hypothesis** Thus I would suggest the answer is linkage in this example, given the reference genome isn't a parent. The mutations are occurring on different strands therefore better represented as 0/1 and then 1/0, i.e. one strand comprises 'A' allele and the opposite strand is the 'G' allele.

**Why is unphased linkage important** In this scenario the interpretation could be hugely important because it would be stating that each allele has come from different parents. This would be more important than maintaining convention. In an inherited disorder there is a risk the child has a dominant phenotype, when each parent could have been heterozygous. Clearly it depends whether the mutation is forcing an amino acid mutation, at 11bp it is certainly possible, i.e. mutation in the first codon upstream and mutation in the second codon downstream.

More generally in any population analysis, linkage is an enormously important issue and central to any calculation. The point is simply that most analysis infers independent segregation of each allele and in this example that remains a good assumption because the allele (I am suggesting) are on opposite strands. If they were on the same strand it would create a bias away from panmixia.

The hypothesis is formed on this example alone and if the alleles were 10 megabp apart it would be impossible to conclude they were linked, therefore the convention of 0/1 should hold. That's my guess, its not a bad guess.

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


# Question
Title: How do I use bash for loops to loop through a WDL Array?
Within the example below a WDL array of strings, `test` is defined as `"a","b","c","d"` the value `total` is a variable that signifies the length of the WDL array. Using the bash for loop, which will iterate from 1 up to length of the array, I seek to print out the values at each indicy in the array and output it to another array. For example when the loop is on 1 the value at indicy 1 within the array will be printed out in this case `"a"` and stored in the array out, if the loop is on 2 then the value `"b"` will be printed out and stored in the array out, etc. I tried the below example, however, I get an error when running. Can anyone help me with this problem?

How I tried to do this in WDL:

```
task A {
    input {
        Array[String] test = ["a","b","c","d"]
        Int total = length(test)
    }

    command {
       for (( c = 1; c <= ${total}; c++ ))
       do
           echo ${test[$c]}
       done
    }
    output {
        Array[String] out = read_lines(stdout())
    }
}

```

# Answer
> 3 votes
I suspect the issue here is that you are trying to use a bash variable inside of a WDL placeholder. The `${}` are interpreted by WDL as places where some value should should be inserted into the command. So when you say `echo ${test[$c]}`, `test[$c]` is seen a WDL expression. This isn't a valid WDL expression, causing a parsing error.

The choice to use `${}` as placeholder syntax in WDL was a rather inconvenient decision by the WDL team, because it clashes with the bash syntax. There is an alternative syntax available for the command section in WDL 1.0, which is preferred for this very reason. This alternative uses `~{}` as placeholder syntax, instead of `${}`, and wraps the command in `<<< >>>` (referred to as the "heredoc-like syntax"), rather than the usual `{ }`. In this heredoc-like syntax the `${}` isn't seen as a WDL placeholder and can, thus, be used in bash commands like normal. You can read about it in the WDL SPEC.

Using this syntax your WDL would look something like this:

```
version 1.0

workflow test {
    call A
}

task A {
    input {
        Array[String] test = ["a","b","c","d"]
        Int total = length(test)
    }

    command <<<
       ARRAY=(~{sep=" " test}) # Load array into bash variable
       for (( c = 0; c < ~{total}; c++ )) # bash array are 0-indexed ;)
       do
           echo ${ARRAY[$c]}
       done
    >>>

    output {
        Array[String] out = read_lines(stdout())
    }
}

```

Note that if the values in your array are more complex and might include characters with special meaning in bash (like spaces), it would probably be better to wrap the values in single quotes like so:
`ARRAY=('~{sep="' '" test}')`

And just in case you weren't aware, bash for-loops can also loop over arrays directly:

```
command <<<
       for c in ~{sep=" " test}
       do
           echo $c
       done
>>>

```

---
Tags: bash, wdl
---


title: BEAGLE do not see and/or not use OpenCL

text:
# Question
Title: BEAGLE do not see and/or not use OpenCL
I am conducting phylogeny analyses using BEAST 2 as software tool. It is said that it can be accelerated by BEAGLE-library. I have a personal laptop with Windows & a cluster with CentOS.

I installed both BEAST 2 & BEAGLE onto them. This gives me a boost by 3.6 times. It was said too that using OpenCL with BEAGLE gives me more speed-up. Nonetheless, my calculations takes from twelve hours to 3 days so I got really interested in this. Installing it on Windows were very simple. When I ran `beast -beagle_info`, I saw:

```
                        BEAST v2.6.6, 2002-2021
             Bayesian Evolutionary Analysis Sampling Trees
                       Designed and developed by
 Remco Bouckaert, Alexei J. Drummond, Andrew Rambaut & Marc A. Suchard

                   Centre for Computational Evolution
                         University of Auckland
                       r.bouckaert@auckland.ac.nz
                        alexei@cs.auckland.ac.nz

                   Institute of Evolutionary Biology
                        University of Edinburgh
                           a.rambaut@ed.ac.uk

                    David Geffen School of Medicine
                 University of California, Los Angeles
                           msuchard@ucla.edu

                      Downloads, Help & Resources:
                           http://beast2.org/

  Source code distributed under the GNU Lesser General Public License:
                   http://github.com/CompEvol/beast2

                           BEAST developers:
   Alex Alekseyenko, Trevor Bedford, Erik Bloomquist, Joseph Heled,
 Sebastian Hoehna, Denise Kuehnert, Philippe Lemey, Wai Lok Sibon Li,
Gerton Lunter, Sidney Markowitz, Vladimir Minin, Michael Defoin Platel,
          Oliver Pybus, Tim Vaughan, Chieh-Hsi Wu, Walter Xie

                               Thanks to:
          Roald Forsberg, Beth Shapiro and Korbinian Strimmer

--- BEAGLE RESOURCES ---

0 : CPU
    Flags: PRECISION_SINGLE PRECISION_DOUBLE COMPUTATION_SYNCH EIGEN_REAL EIGEN_COMPLEX SCALING_MANUAL SCALING_AUTO SCALING_ALWAYS SCALERS_RAW SCALERS_LOG VECTOR_SSE VECTOR_NONE THREADING_NONE PROCESSOR_CPU FRAMEWORK_CPU

1 : NVIDIA GeForce 940MX
    Global memory (MB): 2048
    Clock speed (Ghz): 1.19
    Number of cores: 384
    Flags: PRECISION_SINGLE PRECISION_DOUBLE COMPUTATION_SYNCH COMPUTATION_ASYNCH EIGEN_REAL EIGEN_COMPLEX SCALING_MANUAL SCALING_AUTO SCALING_ALWAYS SCALERS_RAW SCALERS_LOG VECTOR_NONE THREADING_NONE PROCESSOR_GPU FRAMEWORK_CUDA

2 : Intel(R) HD Graphics 620 (OpenCL 2.1 )
    Global memory (MB): 3219
    Clock speed (Ghz): 1.00
    Number of compute units: 24
    Flags: PRECISION_SINGLE PRECISION_DOUBLE COMPUTATION_SYNCH COMPUTATION_ASYNCH EIGEN_REAL EIGEN_COMPLEX SCALING_MANUAL SCALING_AUTO SCALING_ALWAYS SCALERS_RAW SCALERS_LOG VECTOR_NONE THREADING_NONE PROCESSOR_GPU FRAMEWORK_OPENCL

3 : Intel(R) Core(TM) i5-7200U CPU @ 2.50GHz (OpenCL 2.1 (Build 10))
    Global memory (MB): 8067
    Clock speed (Ghz): 2.50
    Number of compute units: 4
    Flags: PRECISION_SINGLE PRECISION_DOUBLE COMPUTATION_SYNCH COMPUTATION_ASYNCH EIGEN_REAL EIGEN_COMPLEX SCALING_MANUAL SCALING_AUTO SCALING_ALWAYS SCALERS_RAW SCALERS_LOG VECTOR_NONE THREADING_NONE PROCESSOR_CPU FRAMEWORK_OPENCL

```

However, it still uses only resourse #0, even through forced `beast -beagle_order 3 -threads 1 filename.xml`. Although, in comparison to Linux case, it was more successful.

I have downloaded intel\_sdk\_for\_opencl\_applications\_2020.3.494.tar.gz from here. Then I unzipped it & run **install.sh**. Installation was successful. After that, I have decided to go to **build** folder in **beagle-lib** & run `cmake -DCMAKE_INSTALL_PREFIX:PATH=$HOME ..` in order to go on to run `make install` but I get the next message:

```
-- JAVA_HOME=
-- JNI_INCLUDE_DIRS=/usr/lib/jvm/java/include;/usr/lib/jvm/java/include/linux;/usr/lib/jvm/java/include
-- JNI_LIBRARIES=/usr/lib/jvm/java/lib/libjawt.so;/usr/lib/jvm/java/lib/server/libjvm.so
-- Not using libtools for plugins
-- Could NOT find OpenCL (missing: OpenCL_LIBRARY OpenCL_INCLUDE_DIR)
CUDA_TOOLKIT_ROOT_DIR not found or specified
-- Could NOT find CUDA (missing: CUDA_TOOLKIT_ROOT_DIR CUDA_NVCC_EXECUTABLE CUDA_INCLUDE_DIRS CUDA_CUDART_LIBRARY)
-- Configuring done
-- Generating done
-- Build files have been written to: /export/home/rinkman/beagle-lib/build

```

It was solved through manual specification of directories `cmake -DOpenCL_INCLUDE_DIR=~/intel/system_studio_2020/opencl/SDK/include/ -DOpenCL_LIBRARY=~/intel/system_studio_2020/opencl/SDK/lib64/libOpenCL.so.1.2 -DCMAKE_INSTALL_PREFIX:PATH=$HOME ..`. The compiling was successful after that. `make test` was executed without any problems. I detected the file **libhmsbeagle-opencl.so.40.0.0** in **~/lib folder**. But `beast -beafle_info` did not show OpenCL. As there are older versions of BEAST & BEAGLE on the root, I tried to change very roughly startup script to run only my local BEAGLE but nothing changed:

```
#!/bin/sh

if [ -z "$BEAST" ]; then
    ## resolve links - $0 may be a link to application
    PRG="$0"

    # need this for relative symlinks
    while [ -h "$PRG" ] ; do
        ls=`ls -ld "$PRG"`
        link=`expr "$ls" : '.*-> \(.*\)$'`
        if expr "$link" : '/.*' > /dev/null; then
        PRG="$link"
        else
        PRG="`dirname "$PRG"`/$link"
        fi
    done

    # make it fully qualified
    saveddir=`pwd`
    BEAST0=`dirname "$PRG"`/..
    BEAST=`cd "$BEAST0" && pwd`
    cd "$saveddir"
fi

BEAST_LIB="$BEAST/lib"

if [ -z "$JAVA_HOME" ]; then
  JAVA=java
else
  JAVA="$JAVA_HOME"/bin/java
fi

# use BEAGLE_LIB if the BEAGLE library is not in a standard location
if [ -n "$BEAGLE_LIB" ]; then
    if [ -n "$BEAST_EXTRA_LIBS" ]; then
        BEAST_EXTRA_LIBS="$BEAST_EXTRA_LIBS:$BEAGLE_LIB"
    else
        BEAST_EXTRA_LIBS="$BEAGLE_LIB"
    fi
fi

# Explicitly add /usr/local/lib to library search path to ensure
# beast continues to find beagle when installed here.  (This is
# necessary due to known problems with certain JREs.)
if [ -z "$LD_LIBRARY_PATH" ]; then
    export LD_LIBRARY_PATH=/usr/local/lib
else
    export LD_LIBRARY_PATH="$LD_LIBRARY_PATH":/usr/local/lib
fi

BEAST_EXTRA_LIBS="$HOME/lib/"

# use BEAST_EXTRA_LIBS variable to load BEAGLE and other libraries from non-default locations
# this assumes that the library path contains all these libraries (or are set through LD_LIBRARY_PATH)
#if [ -n "$BEAST_EXTRA_LIBS" ]; then
#  "$JAVA" -Dlauncher.wait.for.exit=true -Xms256m -Xmx1g -Djava.library.path="$BEAST_EXTRA_LIBS" -Duser.language=en -cp "$BEAST_LIB/launcher.jar" beast.app.beastapp.BeastLauncher $*
#else
#  "$JAVA" -Dlauncher.wait.for.exit=true -Xms256m -Xmx1g -Duser.language=en -cp "$BEAST_LIB/launcher.jar" beast.app.beastapp.BeastLauncher $*
#fi
"$JAVA" -Dlauncher.wait.for.exit=true -Xms256m -Xmx1g -Djava.library.path="$BEAST_EXTRA_LIBS" -Duser.language=en -cp "$BEAST_LIB/launcher.jar" beast.app.beastapp.BeastLauncher $*

```

Please, could anyone help with this?

P.S. I am a noviciate in unix-systems.

# Answer
> 1 votes
I better explain this style of parameter rich MCMC takes ages and the issues are how speed this up.

---

Having given this some thought my suggestion is to install `miniconda` or `anaconda` and then install the port to OpenCL, then install beast and beagle.

Anaconda is here: https://www.anaconda.com/products/individual

Then install openCL via conda:

```
 conda create -n beasty python=3.9
 conda activate beasty
 conda install -c conda-forge pyopencl
 conda install -c bioconda beagle
 conda install -c bioconda beast2

```

My idea was you'd access OpenCL via python described here. This wouldn't work ... its complicated.

**Traditionally OpenCL and Beagle didn't work ... on an Intel Mac**

> 2 : Intel(R) HD Graphics 620 (OpenCL 2.1 )

However, if this is your chip set - I wouldn't give up. In any other scenario I wouldn't pursue this - personally, below.

One alternative is NVIDA e.g. Cloud computing used via CUDA and Beagle. CUDA will easily install and Beagle stacks straight on top.

Any heavy duty MCMC is very amenable to traditional CPU parallelisation - if you think at it - it would be and thats why GPU - CUDA will work. Personally, looking at your resources I would parallelise under you cluster using a MINIMUM of 8-cores per calculation outside CUDA/OpenCL. If that has not already been done, that is the solution in context in my opinion. Unless you're a Google developer few will have access to a NVIDA cluster in any case sufficient CPU parallelisation should be able to achieve the runtime close to that of a good GPU (its gonna be a lot of CPUs though).

**Conclusion** Honestly, my personal advise is that phylogenetics sump clusters no matter what the calculation is being used and that has always been the case. The best solution is to max out on available cores. When eukaryote genome guys complain about 1000 hour runtimes (obviously thats the total CPU parallelised), thats a routine calculation in phylogenetics. Buying in NVIDA doesn't necessarily solve this, because whilst MCMC is massively accelerated for maximum likelihood isn't accelerated to that extent by comparison.

**Extra information** In case the OP is not aware, the enormous advantage of Beast2 is the ability to resume an MCMC at the point it ended (or some checkpoint therein), you don't need to run the whole thing again. You should seek to exploit that, so a calculation started on a laptop can be finished on a cluster. I can't state sufficiently how powerful that is, it stands to reason because the random number stream is declared and every sample thereafter. Keeping all calculations is standard practice in any case, but especially with a Beast2 run and don't every discard it (hard disk is cheap right?). You might find you have to return to it and its a whole lot better than needing to do the whole thing again.

*Its time to move on.*

---
Tags: phylogenetics, linux, beast
---


# Question
Title: Create a NHX phylogeny format with r or python
I need help to add information in my phylogenetic tree on each branch.

I have two pieces of information at my disposal:

1. A phylogenetic tree in newick format :

(mammal:0.14,turtle:0.02,(rayfinfish:0.25,(frog:0.01,salamander:0.01)50:0.12)95:0.09);

2. A matrix with the edge number and information in the Information column:

> ```
> Edge    Inf
> 7        0
> 6        0
> 5        2
> 4        1
> 3        2
> 2        2
> 1        2
>
> ```

Does someone have an idea in order to add these informations by creating a `NHX format` in python (with **ete3** for instance) or with an **R package** ?

Here is the tree :

# Answer
> 1 votes
**Brief answer** Firstly this is not a final solution. The simple answer is to use TreeNode.add\_feature() in ETE3 which would enable NHX format to be inserted into the tree. There are two possible ways to implement the solution, but I'm not aware of a 'reader'. I'll return to this later under the section 'Question (again)'.

**Note** I assume that the OP made a typo with 'witg ETE3' and probably means 'with ETE3'?

**Interest** However as an approach to molecular phylogeny it is interesting theory and has clear application. Not in a simple tree, but when the volume of data overwhelms optical analysis. I am not aware of this being deployed in phylogeny - not done a literature search though, but would readily automate processing of complex trees.

**Theory** The OP uses the term 'edge' rather than branch, so this is referring to depth first search (DFS) type style calculations. In Phylogenetics there are 'branches' and 'nodes' and these are distinct, in this line of theory, I think, 'edges' can be nodes ... possibly, it is possible an 'edge' connects two nodes, in which an 'edge' is synonymous with a 'branch'.

The difficult I have with the information is that edge 0 is not defined. In Phylogenetics edge 0 would probably be the root. It is possibly the 'missing branch' in the trifurcation at the base of the tree. In this line of theory however 0 may simply be a null category, because an edge is not necessarily the same thing as a branch.

**Data structure** The information present is paired. thus ... salamander is 7,0 ... frog is 6,0 ... turtle is 4, 2... We can work this out because the data structure is 7,6,5,4,3,2,1 for odd numbers.

7-\>0, 6-\>0, 5-\>2 and so on ...

Moreover, 7 is at the top of the tree whilst, 1 is at the bottom and the 'edge' is an arbitrary designation. Phylogenetics wouldn't label nodes in this format, labelling would usually be with respect to the base of the tree (root). If I look at this from my world-view there's no 0-\> ... but 5 the ancestor of 7 and 6 goes to turtles.

Thus far the information isn't the DFS, its presented as the input, for example to calculate DFS, and can be used calculate all sorts of things in this style of graph theory.

**Question (again)** Back to the question, I am not aware of DFS style theory being used in trees, therefore I would not be aware of an interpreter that would read it. If the OP was just wanting to stash the data into a tree via NHX it could be stacked onto the root lineage. There isn't a root lineage in this tree so I've just added a root at '1:0.01', i.e. a nominal distance to provide the OP with a place to stash the data there. This will read as a tree and the additional 1:0.01 doesn't affect the integrity of the tree - especially in this given phylogeny where the location of the root is not in dispute.

**Idea 1**

As described above the Edge data gets stored in the 1:0.01 inserted as an additional bracket at the end of the tree.

```
((mammal:0.14,turtle:0.02,(rayfinfish:0.25,(frog:0.01,salamander:0.01)50:0.12)95:0.09)1:0.01);

```

**Idea 2**

This is pairwise data and NXH format could insert the pairwise data for each edge (via TreeNode.add\_feature()) for each 'edge' thus salamander would have 7-0 or 7,0 inserted using NHX. That is absolutely technically where the data belongs, but how this read later is more challenging.

**If I get chance** I'll update the question later with specific NHX output and possibly code to achieve it.

---

OP upvoted. Its a good question with interesting application.

# Answer
> 0 votes
The treedata book has a section on tree annotations that should work for you.

---
Tags: r, python, phylogenetics, phylogeny
---


title: Target genes for piRNA

text:
# Question
Title: Target genes for piRNA
Where I can find a database or tool to give me the `target genes of PIWI (piRNAs)` in `human`?

I found one but works for worm like `piRTarBase`

There are a lot for miRNAs around but not for piRNAs

# Answer
## Update

It does indeed appear that while piRBase etc. have human piRNAs, they do not have target information.

Looking at related questions, recommendations are made to map piRNA sequences to the genome to find the targets. One tool recommended is PiPipes, which might assume read input but otherwise might give some ideas.

## Original answer

piRBase has this. See "target mRNA" field.
> 1 votes

---
Tags: rna-seq, gene, genome, non-coding
---


# Question
Title: What can be the reason of getting negative branches lengths after BEAST analysis?
BEAST2 is currently being used for tree reconstruction prior phylogeographic analysis. The sample size and loci are described below.

I thought that BEAST/BEAST2 does not allow negative lengths of branches on prior distribution level so I was surprised when I have looked at my consensus tree using **FigTree**.

**Question** What could be the reason of this and how to wrestle with this phenomenon?

I am analysing about 1500 taxa of HAV (subtype: IA, region: VP1) under the next parameters:

* GTR+Г,
* Strict Clock,
* Coalescent constant population.

The tree was generated by **TreeAnnotator** through

```
~/beast/bin/treeannotator -heights median -burnin 20 -limit 0.5 VP1_test_bt_ExpPop.trees VP1_test_bt_ExpPop.tree

```

# Answer
I found the answer on ResearchGate from Santiago Sanchez-Ramirez & let me quote it here:

> Try using the option `-heights ca` when you run TreeAnnotator. This stands for "common ancestor trees" and aims at summarizing clade ages across all posterior trees and not only the values for subset of trees that have that clade. If that subset of trees is low in frequency the average node height might end up being older than the direct ancestor. Here is the explanation for negative branches: "MCC trees produced by TreeAnnotator can have a descendent node that is older than its direct ancestor (a negative branch length). This may seem like an error but is actually the correct behaviour. The MCC tree is, by default, generated with average node heights across all trees in the sample which contain that clade. The negative branch lengths result when a clade is at low frequency and tends not to occur in those trees that have the MCC tree's ancestral clade (or vice versa). This means the average heights are for the adjacent nodes are derived from different sets of trees and may not have any direct ancestor-descendent relationship." In this paper you will find more information about the `ca` option: Looking for trees in the forest: Summary tree from posterior samples
> 2 votes

# Answer
I am avoiding directly answering the question because in summary, I feel, this given issue is akin to project supervision and is simply too much time.

One issue raised however is something called 'burnin'.

**Burnin is foundational to MCMC**.

In the comments I stated the burnin is too low. The burnin used ...

```
~/beast/bin/treeannotator -heights median -burnin 20 ...

```

The OP believes the burnin represents a percentage, i.e. 20%.

Citing the documentation on TreeAnnotator at https://beast.community/treeannotator

> BurnIn (as trees) In this case you just specify the actual number of trees you wish to discard rather than the number of states.

In this given case, considering this is 1500 taxa, I would estimate the actual burnin used in comparison to an initial assessment of completion time (technically called 'convergence') to be 0.4% maximum, i.e. well over 0.5 log smaller than the burnin believed to be used.

It is easy to assess empirically e.g. comparing `wc infile` and `wc outfile`, where infile and outfield are MCMC treefiles.

A standard burnin is 10% of the total MCMC. The threshold would be independently assessed and if inadequate a higher burnin would be used, of which 20% is reasonable.

To explain the documentation, in an MCMC only a sample of the trees are written to disk (often 1 in 1000, but it can be anything) and the burnin is 10% of the total (unsampled).

Thus factoring the sampling is important in calculating the precise burnin in any MCMC output analysis. Otherwise the burnin that was intended, isn't the burnin used in practice.

Perhaps you can now see my point? Not the specific technical matter above, but the issue concerning complexity and time commitment.

---

To address the matter concerning Beast1 vs Beast2 raised as a comment below, i.e. the proposal the original MCMC was performed in Beast2 and Beast2 is exempt. The issue under discussion is TreeAnnotator.

The TreeAnnotator section of the Beast2 documentation is here https://www.beast2.org/treeannotator/

> TREEANNOTATOR ....

> BurnIn This option allows you to select the amount of burn-in, i.e., the number of samples that will be discarded at the start of the run, so that you are only analysing the part of the trace that is in equilibrium.

Again the Beast2 documentation is clear, 'the number of samples', not the percentage or proportion. Furthermore note it states 'samples', not 'states', i.e. the number of trees resulting from the sample frequency that was assigned.

A burnin can be performed prior to TreeAnnotator analysis and that is good practice, albeit not essential. This might be where the report output the OP cites was obtained. Furthermore, if this was done here upstream, then no burnin is required for TreeAnnotator. Moreover if a correct upstream burnin was performed, burnin should either be; avoided for TreeAnnotator, or else the total combined burnin calculated and correctly reported. The final output ain't just a tree, there's a lot of important analysis/ 'meta data'.

**Note** I agree with the OP's second comment below this answer using information from the comment below the question: a value of -0.31 can't be ignored, i.e. its a notable value for any branch length and very large in the given context.

*Its definitely time to move on*

OP upvoted, good luck with the project and all the best.
> 0 votes

---
Tags: phylogenetics, phylogeny, beast
---


# Question
Title: How to sort legends alphabetically in plotly dash with category_orders function
I am using Plotly dash for creating a scatter plot. The legends I am getting in the final figure are randomly placed. I want to sort labels alphabetically from A to Z with `category_orders`.

```
fig = px.box(
        df,
        x=selected_x,
        y=selected_y,
        points='all',
        hover_data=hover_data,
        color=colour_by,
        width=800,
        height=600,
        labels={
            selected_y: "{} {}".format(selected_y_gene, selected_y),
            selected_x: "{} {}".format(selected_x_gene, selected_x),
        },

```

# Answer
Before you construct the figure could you try sorting the values in the dataframe in the order you require? Thus set the True/False flag according to what you want.

```
orderedDF = df.sort_values(
  by=col_label,
  ascending=True)

```

then add the following into the plot denoting the labels:

```
sort=False

```

Thus added below.

```
fig = px.box(
    orderedDF,
    x=selected_x,
    y=selected_y,
    points='all',
    hover_data=hover_data,
    color=colour_by,
    width=800,
    height=600,
    labels={
        selected_y: "{} {}".format(selected_y_gene, selected_y),
        selected_x: "{} {}".format(selected_x_gene, selected_x),
# I think it should go here
        sort=False
    }
    )

```

The gist is you are wanting to stop the automated re-ordering and force the dataframe order onto it.
> 2 votes

---
Tags: python
---


# Question
Title: How to determine NCBI's SRS google cloud bucket or AWS bucket
It appears that AWS and GCP host SRA data and it's beneficial to grab the data from that source when running on GCP for example.

Given an SRR accession like SRR1929796 https://trace.ncbi.nlm.nih.gov/Traces/sra/?run=SRR1929796

Is there a metadata API / source that can tell me which google cloud bucket to retrieve the data from like they show on the website?

```
gs://sra-pub-crun-4/SRR1929796/SRR1929796.1

```

It also sounds like fastq-dump and fasterq-dump can be configured with Google Cloud credentials to download that for me.

I'm just wondering is there a way to determine the google cloud storage bucket path it's in `gs://sra-pub-crun-4`?

# Answer
It sounds like it to me too. I have never directly accessed NCBIs stuff via Google buckets and I don't know what permissions they are setting for direct access.

So ... I briefly checked it out

```
gsutil ls gs://sra-pub-crun-4

```

and 'yes' its possible but ...

> BadRequestException: 400 Bucket is a requester pays bucket but no user project provided.

'Requester pay' luckily I wasn't logged in :-)

What the message implies given access to your (my) GCP account and they'll bill you (me). I don't think it will be much, particularly for a 'list' request, however I ain't gonna find out for a data set I don't need. I do doubt I'm going to get billed for 'list' request, however I've never paid 'list' request I prefer to keep it that way. In contrast, NCBI's front-end is free and thats what they want, because they could of course set the bucket so they pay.

So in answer the question, if you wanna pay go ahead ... dial up a GCP account link it to pay.

---

Here's my first answer

Please note I'm not really certain what information the OP seeks here, e.g. gs:// locations type stuff. *If the question is where's the SRS database on every gs:// path ...* I don't know.

I think the OP wants an API ... Google Cloud certainly use APIs by the bucket load (no pun intended but I thought it was funny). I don't use them because they are excessively verbose and I prefer SDK and gcloud tools, of which `gsutil` is one.

> It also sounds like fastq-dump and fasterq-dump can be configured with Google Cloud credentials to download that for me.

If I have understood correctly and that is all you want then sure

```
gsutil cp gs://sra-pub-crun-4/SRR1929796/SRR1929796.1 .

```

will do it, if permitted (or in hindsight prepared to pay) and there will be an API version of this.

---

**Background**

Installation via Google's `Cloud SDK`, via https://cloud.google.com/storage/docs/gsutil\_install#install
> 0 votes

---
Tags: fastq, sratoolkit, sra, cloud
---


title: Multiple file Python script

text:
# Question
Title: Multiple file Python script
I am currently doing an analysis requiring a multiple-input script. I have a directory filled with 900+ samples.

Each sample is comprised of two-files in a fastq format. This is an example of a sample:

```
MIP_REV_BC_01_01A_S1_R1_001.fastq.gz
MIP_REV_BC_01_01A_S1_R2_001.fastq.gz

```

The script has a Read1 & Read2 variable, and I need to make it so the R1 file is inputted for Read1 and the R2 file is inputted for Read2 for every sample. To be clear, the script is supposed to be executed on every sample separately.

This was an attempt:

```
with open(filename) as file:
    lines=file.readlines()
    for line in lines:
        parts[]=line.split("_")
        if parts[7]="R1":
            r1[].append(line)
        else if parts[7]="R2":
             r2[].append(line)

```

We were bogged down by syntax errors, but either way, this script doesn't seem to do what we wanted it to. However, it captures the idea of what we're trying to achieve.

**Question** I want to input a directory into the function so it automatically runs the script on the samples in the directory.

Any idea what I should do?

# Answer
> 2 votes
Edited with explanation below. I don't really understand the objective, because I initially thought you were editing the content of fastaq file, but you seem to be editing a file which contains fastaq filenames. Is that correct? I know the syntax but I don't know really what the goal is.

**Note** I'm using the pathlib which has a glob function, do consider `rglob` under pathlib.

```
import re
import pathlib as Path

# pathlib is super, good to habit to use it
mypath = Path('/Users/fewbugs/workingd')
allfilenames = mypath.glob("*.fastaq.gz")
r1, r2 = [], []
x = re.compile('_')
for filename in allfilenames:
    # the 'r' is super important
    with open(filename, 'r') as file:
    # there is no need for a readlines()
        for line in file:
    # split via re, it is better to compile (faster)
        parts=x.split(line)
    # you must use '=='. REMEMBER lists start at 0 not 1
           if parts[6]=="R1":
    # if you want to append to a list, declare the list first
            r1.append(line)
    # there is no such thing as 'else if' REMEMBER lists start at 0 not 1
            elif parts[6]=="R2":
                 r2.append(line)
    # when testing good idea to trap the stuff you're missing
            else:
                 print ("This one escaped %s" % (line))

```

Please be CAREFUL with my spacing, everything needs to be multiples of 4. I've not checked the code.

---

Forget this bit I thought you were processing the content of the fastaq file.

```
outfile1 = Path(mypath.parent, 'outfile1')
outfile2 = Path(mypath.parent, 'outfile2')

with open (outfile1, 'w') as f:
    for item1 in r1:
        f.write("%s\n" % item1)
with open (outfile2, 'w') as f2:
    for item2 in r2:
        f2.write("%s\n" % item2)

```

---

## A simpler way?

... it looks like you're just processing filenames...

Might I suggest just globing no need for all this stuff.

```
mypath = Path('/Users/fewbugs/workingd')
r1filenames = mypath.glob("*_R1_*.fastaq.gz")
r2filenames = mypath.glob("*_R2_*.fastaq.gz")

```

If that is the aim, thats it! Just be careful about the `glob` versus `rglob`.

**Idea** If you want two paired lists this might be better done in a nested comprehension, e.g. a counter could be inserted at the central glob (\*).

That may very well be what you're after. i.e.

```
R1 = ['stuff_R1_001.fastaq.gz', 'stuff_R1_002.fastaq.gz' 'stuff_R1_003.fastaq.gz']
R2 = ['stuff_R2_001.fastaq.gz', 'stuff_R2_002.fastaq.gz' 'stuff_R2_003.fastaq.gz']

```

---
Tags: python, fastq
---


# Question
Title: Suggestion for python graphics library to plot aligned protein sequences
I'm working towards developing some visualization method for multiple protein sequences and subsequences using PYTHON. I want to show the alignment of many fragment sequences to the original sequences as well as the position of fragments and the amino acids. The amino acids are to be colored by some physico-chemical properties, something similar but not exactly as this:

Can somebody suggest some Python libraries suitable for this type of work? The output can be graphics file like png or interactive graphics or html, doesn't matter really for me

Thanks

# Answer
> 3 votes
You might have a look at the biotite package, and more specifically to its `biotite.sequence.graphics` module: https://www.biotite-python.org/apidoc/biotite.sequence.graphics.html#module-biotite.sequence.graphics.

Disclaimer: I haven't used these functionalities myself, apart from using one of their colour schemes for DNA.

---
Tags: python, visualization
---


# Question
Title: Regression plot of a continuous trait - is there a binary equivalent?
Within the GWAS wikipedia page, you can view the following plot:

> "Illustration of a simulated genotype by phenotype regression for a single SNP. Each dot represents an individual. A GWAS of a continuous trait essentially consists of repeating this analysis at each SNP."

This is a regression plot showing on the y-axis the phenotype value for a continuous trait. Is it possible to create a similar plot with a binary trait - I can't think what this would look like, or if there is an equivalent plot I could create.

# Answer
I have not really seen this type of plot used when a binary outcome is being analyzed. I think in part this is due to the fact that the plot will be rather boring and difficult to interpret as there are only a limited set of possible values (4 if the response and independent variable are both binary), and without introducing some level of random noise, all of the data observations with the same X-Y value pair would overlap.

The situation is a little better when the independent variable is continuous: then you can show more of a distribution. In the figure below, pulled from Figure 2 of this paper explaining how to visualize the discrimination abilities of a logistic regression model, the authors have plotted the risk score produced by the logistic regression model for each patient (x-axis) with the patient's actual outcome on the y-axis. Note that some random noise has been introduced to each observation's y-value so that they do not cover up each other. This plot shows that the risk score is able to somewhat separate the two outcomes.

If you are dealing with a categorical independent variable like the figure included in your question, I think it is better to visualize the proportion of the outcome for each independent variable category like this:

Here the outcome is whether the borrower defaulted on their loan (1 = yes, 0 = no) and the independent variable is whether they are a student (yes/no). For the y-axis, I have calculated the proportion of loan recipients in each group that defaulted. This could easily be done with GWAS data such as genotype vs whether the patient died of heart disease.

I also suggest you check out this answer on how to visualize a contingency table, as the data you are describing (assuming the independent variable stays categorical) is a contingency table.
> 1 votes

---
Tags: snp, gwas, logistic-regression
---


# Question
Title: RNAseq small sample size in control group
I am doing an RNAseq experiment with affected (n=6) and control groups (n=6) in cow. However it turned out that 4 of my control samples have very low quality so I have to discard them. My problem is that there is no way to repeat the sampling and sequencing. Are there any solutions to perform DE,GSEA with these conditions (somehow generating bootstrapped samples from control groups or something similar)? Thanks in advance!

# Answer
> 3 votes
There is no obstacle in performing a DE analysis and all downstream on a 6 vs 2 comparison. The only thing to keep in mind is that the representation of the control is not particularily good, but you can technically perform the exact same analysis as if it was 6 vs 6, e.g. with DESeq2, edgeR or limma-voom/trend.

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


title: Low mapping reads for RNA sequence

text:
# Question
Title: Low mapping reads for RNA sequence
I am new to RNA sequencing analysis. I have RNA samples for bacteria for different treatment. I did contaminant filtration, remove adapters and checked fastQC report: attached below for one sample (fwd and rev) .  After mapping with the reference genome using bbsplit but the mapping is only 0.62%.

Any help/suggestions highly appreciated.

Thanks

# Answer
> 1 votes
In addition to swbarnes2 comment, I'd suggest a more systematic approach.

Kraken2 is a tool for metagenomic analysis of reads, though it can be hard to set up.

mash screen is another method for metagenomic analysis, much easier to set up but noisier.

If you are dealing with contamination, treating your data as metagenomic might help.

Also, as always, ensure that your workflow and reference assembly are correct!!!

---
Tags: rna-seq, read-mapping
---


# Question
Title: LD Score Regression Derivation hard to follow
I am trying to understand the derivations from Sullivan et al. (2015) in the Supplementary Material. There, it is mentioned in the first page that the least squares estimate of the j-th SNP effect, considering the polygenicity linear model $φ=Xβ + ε$, is $\hat{β\_j} := X^T\_j φ/N$ . Normally, the least squares solution of $\hat{\beta}$ is $(X^TX)^{-1}X^T\phi$. By comparing the two equations with each other, this means that, in this case, $X^TX=N I\_{M\times M}$. However, LD is assumed, therefore $E(X\_{ij}X\_{ik})\neq 0$, so the covariance matrix cannot be diagonal. Could someone explain to me how the derivation for the beta estimate was performed, please? What am I missing? I checked that the same thing is mentioned in more than one sources such as here.

# Answer
> 2 votes
Not entirely clear on what you're asking but the beta estimates in LDSC come from GWAS summary statistics. As I understand it, they aren't "derived" so much as just set equivalent to what they linear model defines them as, which is their additive effects on the phenotype as implied by the model. So for each SNP j, it's effect size (beta) is whatever difference one copy of the allele makes on the phenotype (phi), divided by N, the number of SNPs.

---
Tags: statistics, gwas, phylogenetics
---


title: Processing proteomics data

text:
# Question
Title: Processing proteomics data
*\[this question was also asked on BioStars\]*

In relevance to my research, I've been looking for proteomics data (control vs. diabetic) and I found a dataset in the article "Diabetes causes marked inhibition of mitochondrial metabolism in pancreatic β-cells".

I am actually looking for the data that was used for creating figure 1. I checked the PRIDE dataset shared in the paper (https://www.ebi.ac.uk/pride/archive/projects/PXD012979) and I could find the raw files (RAW file, proprietary format) but I was not sure where to look for the processed data (does the .mzid file contain precessed data) and how to generate figure 1 (added below).

I checked a few workflows for processing the raw data (e.g. Progenesis and Knime- OpenMS). But the experimental conditions are missing in the PRIDE link. I am not sure how to get started with this. Could someone please give some suggestions?

EDIT: I understand from the answer posted below, one could use graph capture to digitize the values of relative abundance and create the heatmap. I am actually interested mainly in the enzymes involved in the glycolytic pathway Fig. 1 c,d of the paper (image shown above). Unfortunately, I couldn't find the relative abundance reported for proteins involved in this pathway.

# Answer
You asked two questions

1. Is it possible to find from the dataset (i.e does the raw file contain information on the sample whether it is a diabetic sample or a control sample?).
2. Should the diabetic and control dataset be normalized separately?

Before I continue further I'd just like to point out the graph in the Supplementary information of the Nature research article of the paper in question. This graphic has not been edited in anyway.

I will not comment further but it would be very unusual for an analyst to produce a y-axis as done here. Just to labour/labor the point this is the appendix of a Nature research article (not 'even' a Nature letter).

The answer to the first question is you can't and will need to contact the corresponding author. There is a possible work around. The reason you need to contact the authors is:

1. The authors have not labeled the x-axis in any heatmap figure. This is not standard procedure and in this analysis you cannot link the heatmap values for that given mouse to the data base - because the authors have not declared which mouse is which. In addition, one cannot link the RNA-seq analysis to the proteomics data and they are paired samples, because only 3 diabetic mice are presented for the RNA-seq data but 4 mice for the proteomics data, therefore again it is not possible to link the two data sets (I understand this last point is not the precise question, but is important in the 'conclusion' below).
2. The link you provided had 12 samples - not 8 as in the paper. It is very important to understand why there is a different number of samples (where one 'sample' is one mouse) in the proteomics data base compared with the apparent lower number samples presented in the paper.

To reiterate, only the authorship knows this information therefore they need to be contacted. The corresponding author is different from the data base curator who has been previously contacted (with no response). The corresponding author has to uphold the terms of publication as described at the end of the paper.

However, you still cannot link these values to the publicly available data sets, therefore until x-axis sampling is declared it is not a singular solution.

**Question 2**

Normalisation, I would like to point out in the appendix the authors transformed the RNA-seq data (log2) but did not do this for the proteomics data. For example aquaporins were upregulated 94 fold for proteomics data, but 4.61 fold using log2 RNA-seq data. The reason they did this is not clear and understanding this is quite important to the analysis. Why present an untransformed ratio in one analysis, but for the same experiment transform the values for precisely paired data?

The Appendix describes this information in the Supplementary Information in the following table.

**Conclusion** This is a Nature paper and the authors are therefore under obligation, as per terms of publishing, to firstly supply all information on request and any information present in the published article must be explained in the Supplement Information. Those terms are stated at the end of the paper. It is clear the absence of labelling the x-axis in the heatmap in anyway (even if it is done in the legend), or the inability to link pairwise RNA-seq to proteomics data (which is very important), or the discrepancy in the samples on the database verses the samples in the paper, is not explained in Supplementary Information (Appendix).

Under these technical circumstances the authors are obliged to supply the missing information on request. You may require a chain of request escalation within your organisation because they may not respond to the initial requests, but you don't know unless you try. However, at some level in the organisation, e.g. Director, the terms of publication would be upheld by authorship in question.

---

First answer

> But the experimental conditions are missing in the PRIDE link

The experimental conditions of the heatmap are as follows: each box is a separate animal and the heat map colors (colours) are a log2 representation of Figure A (below). In other words to understand the heat map using Figure A is the starting point. You can obtain the values from this analysis easily using a graph capture package using Figure A (not Figure B). Its a long while since I've used graph capture, so I forget the packages. If you feel there are anomalies you can then contact the authors for the precise values.

log2 is a bit weird, natural log would have been preferred:

FIGURE A "Abundance of the indicated proteins involved in lipid synthesis, measured by proteomics, in islets isolated from control (black, Ctrl, n=4) and 2-week diabetic βV59M (white, Diab, n=4) mice. Each data point indicates a separate mouse. Mean±S.E.M. \**p\<0.01, \***p\<0.001. HMGCS1, 3-Hydroxy-3-methylglutaryl-CoA synthase 1; HMGCR, Hydroxy-3-methylglutaryl- CoA reductase; DHCR7, 7-dehydrocholesterol reductase; ACLY, ATP citrate lyase; GPD1, glycerol phosphate dehydrogenase; ACC1, acetyl-CoA carboxylase 1; AACS, acetoacetyl- CoA synthetase; FASN, Fatty acid synthase."

FIGURE B "Heat maps of mRNA and protein expression of the indicated lipid metabolism genes in islets isolated from control and 2-week diabetic bV59M mice. Each box corresponds to a different animal. Colour indicates log2 fold-change."

https://static-content.springer.com/esm/art%3A10.1038%2Fs41467-019-10189-x/MediaObjects/41467\_2019\_10189\_MOESM1\_ESM.pdf

Data availability RNA data https://www.ebi.ac.uk/ena/browser/view/PRJEB31793?show=reads Checks out ok.

Proteomics I checked one link only ... http://proteomecentral.proteomexchange.org/cgi/GetDataset?ID=PXD012979 I agree something doesn't look right here just in the way it is reported, because the paper should highly polished and they did miss the link out. This could be word limit issues (Data availability might be included). The curator is listed as roman.fischer@ndm.ox.ac.uk and that would be the person to contact.
> 1 votes

---
Tags: proteins, proteomics, mass-spectrometry
---


# Question
Title: De novo antibody sequencing fromMS/MS Ion Trap proteomics raw signal
I've read about different proteomics tools that can do de novo assembly of MS/MS Ion Trap proteomics signatures for de novo antibody sequencing. Given the way CDR3s of antibodies are created, together with Somatic Hypermutation, can generate a huge diversity in the CDR regions of the heavy and light chains, the reference-based methods can fall short of determining the full sequence of an unknown antibody without a set of reference sequences including those with its CDRs.

Are there any other tools that can determine antibody sequences from MS/MS proteomics raw signals?

# Answer
> 2 votes
A few notes on data generation first:

There are currently two broad approaches for the generation of bottom-up or "shotgun" MS proteomics: data-dependent acquisition (DDA) and data-independent acquisition (DIA).1 In tandem MS (MS/MS), the DDA approach only puts forward certain peptides generated during the first cycle of MS for fragmentation during the second cycle, while with the DIA approach, all peptides generated during the first MS cycle can be fragmented in the second round.

The recently-introduced Orbitrap Fusion mass spectrometry permits various types of MS2 acquisition methods. There are four main MS2 strategies:

* HCD-OT (higher-energy-collisional-dissociation with Orbitrap detection),
* HCD-IT (HCD with ion trap, IT),
* CID-IT (collision-induced-dissociation with IT) and
* CID-OT on Orbitrap Fusion

Reading from the biorxiv preprint of Stitch:

https://www.biorxiv.org/content/10.1101/2022.03.07.483237v1

"With the aid of specialized software packages like DiPS, Supernovo, or ALPS, full and accurate sequences of the heavy and light chains can be reconstructed with the MS/MS spectra of the digested peptides."

Before running Stitch or similar tools, a peptide profiling tool needs to be executed on the .raw or .mgf data alone.

STITCH: https://www.biorxiv.org/content/10.1101/2022.03.07.483237v2

This currently requires the output of analysing .raw data via PEAKS, which is not free unfortunately. Authors currently looking into using other free de novo sequencing software for input into Stitch. The software also allows for FASTA format, which can then be executed on the output of PepNovo, Novor or similar.

DiPS: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5461544/

ALPS: https://www.nature.com/articles/srep31730

"The results show that ALPS was able to assemble three complete monoclonal antibody sequences of length 216–441 AA, at 100% coverage and 96.64–100% accuracy."

DeepNovo: https://www.pnas.org/doi/10.1073/pnas.1705691114

"This field has been actively studied over the past 20 y, and many de novo sequencing tools have been proposed, such as PepNovo, PEAKS, NovoHMM, MSNovo, pNovo, UniNovo, and Novor among others (6–19). The recent "gold rush" into mAbs has undoubtedly elevated the application of de novo sequencing to a new horizon (20–23)"

"DeepNovo considerably outperformed state of the art methods, achieving 7.7–22.9% higher accuracy at the amino acid level and 38.1–64.0% higher accuracy at the peptide level"

Available for download, requires activation via sending email to author, which then points to the company behind PEAKS, which will give you a 15 day evaluation license.

NovoHMM: the download page leads to a missing URL

DeepResolution2: https://www.sciencedirect.com/science/article/abs/pii/S0039914022002119

Capable of analyzing GC-MS data rather than MS/MS Ion Trap

DePS: https://arxiv.org/abs/2203.08820

DeepNovoV2: https://github.com/volpato30/DeepNovoV2

DPST: https://arxiv.org/abs/2203.13132

Link does not work https://github.com/Yan98/DPST

PepNovo: http://proteomics.ucsd.edu/software-tools/531-2/

The new version is from 2012, includes models for CID and LTQ, but not HCD models.

pNovo: http://pfind.org/

Includes models for CID, HCD, ETD and HCD+ETD.

Novor: https://www.rapidnovor.com/

Filling in a form, someone will contact you. Doesn't seem to be free.

METHODOLOGY used so far:

1 - ThermoRawFileParser: Convert .raw data to .gcf

2 - PepNovo: produce list of peptides, one-liner to convert to .fasta

3 - Stitch: analyse given a list of V(D)J+C references from templates and the .fasta file from above.

---
Tags: machine-learning, proteomics
---


# Question
Title: Number of collinear genes are zero while running MCScanX
I am trying to do a synteny analysis using MCScanX. As per requirement I have done all by all blast between two sets of protein from two closely related species. And also prepared the gff file. The file formats are as follows:

cyp\_dan.blast

NP_001347641.1  XP_042591522.1  78.31   166     32      1       1       166     1       162     5e-91   267
NP_001347641.1  XP_042577882.1  60.24   166     62      2       1       166     1       162     4e-63   196
.
.
XP_704272.2     XP_018953536.1  31.03   87      59      1       275     360     90      176     8e-07   50.4
XP_704272.2     NP_001005958.2  33.33   63      38      1       273     331     38      100     1e-06   48.9
XP_704272.2     XP_018969808.1  31.33   83      56      1       275     356     91      173     2e-06   49.3```

cyp_dan.gff
NC_056572.1     gene-LOC1221465.1       2966    7093
NC_056572.1     id-LOC122146575.1       2966    3166
NC_056572.1     id-LOC122146575.1       3246    3357
.
.
NC_002333.2     gene-trnP;Dbxre.1       16527   16596
NC_002333.2     rna-trnP;Parent.1       16527   16596
NC_002333.2     exon-trnP-1;Par.1       16527   16596

cmdline: /home/software/MCScanX/MCScanX cyp_dan
Output:
############### Parameters ###############
# MATCH_SCORE: 50
# MATCH_SIZE: 5
# GAP_PENALTY: -1
# OVERLAP_WINDOW: 5
# E_VALUE: 1e-05
# MAX GAPS: 25
############### Statistics ###############
# Number of collinear genes: 0, Percentage: 0.00
# Number of all genes: 108556
##########################################

Here, both of the species are from same family since there is no chance to have "zero collinear gene".

# Answer
I'm not familiar with the MCScanX tool, but I see two potential issues that are worth investigating.

1. The `cyp_dan.gff` file does not look like GFF or any of its common variants (GFF3, GTF, etc.). Again, I'm unfamiliar with MCScanX, so I don't know precisely how its wants the input data formatted. But if it wants a GFF file, I would double-check to make sure that what you have matches what the tools expects/requires.
2. The gene/protein IDs in the `cyp_dan.gff` file do not seem to match the query or subject IDs in the `cyp_dan.blast` file. No tool is going to be able to find collinear adjacent genes if the IDs in the BLAST results don't match the IDs in the gene annotations.
> 1 votes

---
Tags: blastp, gff
---


title: Seurat Integration

text:
# Question
Title: Seurat Integration
I'm following the instructions for integration: https://satijalab.org/seurat/articles/integration\_introduction.html

And it's taking a while to run:

```
immune.anchors <- FindIntegrationAnchors(object.list = ifnb.list, normalization.method = "SCT",
    anchor.features = features)

```

I'm integrating 6 datasets.

Is this function designed well for large amount of data?

# Answer
> 1 votes
Are you using the same data from the integration tutorial or using your own data? How many cells/genes are in your dataset?

Seurat's default integration method (CCA) is known to be runtime/memory intensive. It can handle large datasets but may require lots of CPU cores/memory. I have found that the amount of RAM also seems to increase as more CPU cores are used. That said, I have been able to run the integration tutorial with the included dataset on an 8GB laptop.

The easiest solution is to simply use more resources when integrating: increase the amount of RAM and consider using multiple cores (there is a separate Seurat tutorial on parallelization). If you are encountering runtime/memory issues and you cannot scale up your resources, I would suggest looking at alternative integration methods:

1. Seurat includes another integration approach, RPCA which is supposed to be less memory intensive at the cost of being more conservative with integration
2. There are several other integration methods that provide similar results. I have used Harmony and found that it requires significantly less resources than the Seurat integration approaches (both CCA and RPCA). A quick google search of single cell RNA-seq integration methods will turn up other popular methods with benchmarks to compare and contrasts different approaches.

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


title: How to create drug synergy plot?

text:
# Question
Title: How to create drug synergy plot?
I want to replicate the synergy drug response curve with the Python library given in the following paper. Could you kindly advice on:

* Suitable graphics libraries for Python
* Code or pseudo code

Apologies because I don't currently have a test code.

I would really appreciate if anyone could guide me in this. Thank you in advance.

Fig 1a

https://www.nature.com/articles/s41586-022-04437-2/figures/1

# Answer
> 1 votes
To get you started. Complex editing in matplotlib isn't intuitive.

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

from numpy import linspace, sin
# you don't need math its just to get the value of pi for this data
import math
import matplotlib.pyplot as plt

xaxis = linspace(1,1.75*math.pi,10)
drug2 = linspace(1,1.45*math.pi,10)

drug1 = sin(xaxis)
drug2b = sin(drug2)
plt.plot(xaxis, drug1, xaxis, drug2b)
plt.legend(['Drug 1', 'Drug 2'])
plt.show()

```

The authors look to have curve fitted and for that you would use `import scipy.optimize as opt`

---
Tags: python, drugs
---


# Question
Title: Help using BioPython's IC_Chain functionality for getting protein 3D coordinates from torsion angles
I have been successfully using BioPython's functionality for extracting protein torsion angles using the `IC_Chain` module. However, I am currently looking to use the inverse functionality as follows:

```
from Bio.PDB import PDBParser, internal_coords

# Read PDB entry.
parser = PDBParser(PERMISSIVE=1, QUIET=1)
structure = parser.get_structure('', 'example.pdb')

# Generate protein internal coordinates (angles).
ic_chain_bound = internal_coords.IC_Chain(structure)
ic_chain_bound.atom_to_internal_coordinates()

# Generate back 3D coordinates from angles.
ic_chain_bound.internal_to_atom_coordinates()

```

but this fails with the following error:

```
Traceback (most recent call last):
  File "test_biopython_internal_coords.py", line 97, in <module>
    main(args)
  File "test_biopython_internal_coords.py", line 89, in main
    ic_chain_bound.internal_to_atom_coordinates()
  File "/software/anaconda/installation/envs/alphafold/lib/python3.7/site-packages/Bio/PDB/internal_coords.py", line 661, in internal_to_atom_coordinates
    verbose=verbose, start=start, fin=fin
  File "software/anaconda/installation/envs/alphafold/lib/python3.7/site-packages/Bio/PDB/internal_coords.py", line 415, in assemble_residues
    Dict[AtomKey, numpy.array], ric.assemble(verbose=verbose)
  File "/software/anaconda/installation/envs/alphafold/lib/python3.7/site-packages/Bio/PDB/internal_coords.py", line 1732, in assemble
    atomCoords = self.get_startpos()
  File "/software/anaconda/installation/envs/alphafold/lib/python3.7/site-packages/Bio/PDB/internal_coords.py", line 1632, in get_startpos
    sp = self.residue.parent.internal_coord.initNCaC.get(self.rbase, None)
AttributeError: 'NoneType' object has no attribute 'initNCaC'

```

Is this not the intended use? Could someone help me figure out why this is failing?

I'm using this PDB file in my example, but the code previously failed for all other PDB files I tried.

Thanks!

# Answer
> 2 votes
It turns out that in order to generate the internal coordinates for a protein structure, one shouldn't manually create an `IC_Chain` object, but rather just use the `atom_to_internal_coordinates()` method of the structure and the `internal_to_atom_coordinates()` one for the inverse functionality.

```
structure.atom_to_internal_coordinates()

# Now you have access to the bond angles and lengths
# of the first chain of structure as follows:
chain = list(structure.get_chains())[0]
ic_chain = chain.internal_coord

# Dictionary mapping 3 AtomKeys to a Hedron
# (with len12, len23, angle fields):
# ic_chain.hedra

# Dictionary mapping 4 AtomKeys to a Dihedraon
# (with angle field for torsion angles):
# ic_chain.dihedra

structure.internal_to_atom_coordinates()

```

# Answer
> 2 votes
not able to grasp the question, tried to figure out a way to show the internal\_coordinates:

```
from Bio.PDB import PDBParser, PICIO, PDBIO

# Read PDB entry.
parser = PDBParser(PERMISSIVE=1, QUIET=1)
structure = parser.get_structure('', 'example.pdb')

print(dir(structure)) # 'atom_to_internal_coordinates' and
                      # 'internal_to_atom_coordinates' seems
                      # to be structure methods

# Generate protein internal coordinates (angles).

structure.atom_to_internal_coordinates(verbose = True)

PICIO.write_PIC(structure[0]['A'], 'internal')

# Generate back 3D coordinates from angles.

pippo = PICIO.read_PIC('internal')

pippo.internal_to_atom_coordinates(verbose = True)

print(pippo)

io = PDBIO()

io.set_structure(pippo)

io.save('atom_coord.pdb',  preserve_atom_numbering=True)

```

compare `example.pdb` to `internal` and `atom_coord.pdb`;

have a look at `internal` file contents if you comment out

`structure.atom_to_internal_coordinates(verbose = True)` (around line 17 of

code, # Generate protein internal coordinates (angles). block-code )

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


# Question
Title: How reliable is Alphafold2 for folding random-coil peptide
I have a peptide X of length 50aa. Alphafold2 predicted the structure which takes a helix form and with high average pLDDT score (\>90).

But during actual experiment (e.g. crystallography, spectroscopy), peptide X is determined to have a *random coil* structure.

My question is how can I resolve the difference between AF2 prediction and crystallography? To what extend we can use AF2 peptide structure prediction for downstream analysis (e.g. docking, etc)?

# Answer
## Possibility 1: technical error

AlphaFold2 uses MSA to infer residue co-evolution (due to epistasis) and thus proximity. A MSA can become contaminated with incorrect homologues, especially shorter protein. *Generally* these will have low scores and be filtered out, but if they don't they will come in like a wrecking ball.

## Possibility 2: conformational switching

Even though on circular dichroism a random coil and helix could not be more opposites, several proteins adopt both conformations and switch. Amyloids are a classic example of this, switching between a beta sheet and a random coil. Synuclein is another example.

A technically problematic (vide infra) alternative option is that the helical conformation is only found when the protein is bound to another.

## Solution: AlphaFold2 testing

Whereas EBI-AF2 provides a single model, the AF2 algorithm can provide multiple possible solutions, which in many cases are actually bona fide alternate conformations. On a protein I tried, the phosphorylations in PhosphoSitePlus (high-throughput) explained wonderfully a second model. The model with the best score is not the "most normal" conformer: not only such a statement makes no sense, but is the one that best obeys the evolutionary signals. I would suggest looking at https://github.com/sokrypton/ColabFold Set the target number of decoys to more than 5 and if all are helices then it's a technical issue. In particular the MSA data would need to be dissected if one really cared. But if nothing comes up... it would be a low reward area of investigation as one would need to discount other less and less likely possibilities.

If the decoys *do* show alternative folds, i.e. coils, then a few things spring to mind that could be tested *in silico* (applicable only to some cases):

* different pH sidechain protonation
* ion binding
* phosphorylations
* binding to nucleic acids
* binding to known protein candidates —not worth exploring if there's no very likely binding partners. If different paralogues bind differ targets it's tricky
* binding to membranes

Looking in PyMOL at the location of the high pLDDT scores and after running APBS (Adv. Poisson-Boltzmann Solver) within PyMOL may reveal something in the helical AlphaFold2 model. Once the conditions which explain the helical model (e.g. bond to DNA as a heterotetramer in the presence of zinc at pH 10 has a very low ∆G), those conditions could be tested in vitro.
> 2 votes

---
Tags: protein-structure, pdb
---


# Question
Title: Biopython - extracting query coverage from XML Blast output
I need to extract (or to calculate precisely) the query coverage from Blast output. I am using NCBIWWW.qblast to query Blast in XML format.

Is it possible to know the query coverage parsing the XML output? If yes, how?

```
result_handle_xml = NCBIWWW.qblast("blastp","pdb",my_query.seq,format_type='XML')

```

here it suggest to use '-6 qcovs', but only for BLAST command line version and only for tabular format

This is the fasta file that I open:

```
>sp|P02489|CRYAA_HUMAN Alpha-crystallin A chain OS=Homo sapiens
GN=CRYAA PE=1 SV=2
MDVTIQHPWFKRTLGPFYPSRLFDQFFGEGLFEYDLLPFLSSTISPYYRQSLFRTVLDSG
ISEVRSDRDKFVIFLDVKHFSPEDLTVKVQDDFVEIHGKHNERQDDHGYISREFHRRYRL
PSNVDQSALSCSLSADGMLTFCGPKIQTGLDATHAERAIPVSREEKPTSAPSS

```

# Answer
> 2 votes
I had the same problem. I'm pretty sure you can calculate it yourself using the following:

1. The length of your sequence (seq\_len)
2. The value from the `<Hsp_query-from>` tag from the hit (`<Hsp_query-from>`)
3. The value from the `<Hsp_query-to>` tag from the hit (`<Hsp_query-to>`)

It can be calculated by:

```
query_cover = (<Hsp_query-to> - <Hsp_query-from> + 1) / seq_len * 100

```

Please note that I did not get this from any documentation. I figured it out by staring at the xml output returned from the Biopython `NCBIWWW.qblast()` function. This solution seems to match with the results I get by blasting the same sequence on the interactive web portal.

# Answer
> 1 votes
```
record = NCBIXML.read(result_handle_xml)
query_len = record.query_length
for alignment in record.alignments:
    _ = alignment.hsps.__next__()
    query_cover = (_.query_end - _.query_start + 1) / query_len * 100

```

---
Tags: python, blast, parsing
---


# Question
Title: Can someone help me estimating the runtime of the pipeline applied by the vertebrate genome project?
# The vertebrate genome project (VGP)

has a lot of interesting publications such as this one.

The rough pipeline is outlined below:

Here the pipeline in more detail:

While the paper describes all the steps of the iterative assembly pipeline, I cannot find information on the **runtime** of their **assembly-pipeline**.

I was excited to check out their github page to find out more in the documentation. Unfortunately did I not find any information on the runtime in the https://github.com/VGP/vgp-assembly repo.

I'have written some e-mails asking collaborators of VGP but I have still not gotten any answers.

Did any of you apply this assembly-pipeline? Could anyone tell me the runtime of the assembly-pipeline?

I hope this question is in the appropriate stackexchange, if not please point me to what you think fits best. Thanks for your suggestions, please be gentle I'm new to all of this.

For more info on the VGP check out this link

# Answer
Note to OP that your tag of @user172818 was not successful in your last comment. That's probably why they didn't respond.

Because no one has gone there so far, I will simply link to HiFiAsm. See linked Nature Methods paper. This tool accepts HiFi reads, which are now standard, and gives you a GFA that you can quite easily convert to FA.

Polishing with Illumina would actively harm the assembly, as it would remove real heterozygosity or add errors where Illumina reads map ambiguously.

Run "out of the box"t with default options on HiFi it will probably give you much better results than VGP, haplotype-resolved.

If you also add Hi-C PE reads with no additional options, it will give you something much much better than VGP.

It's great. See the docs for practical info on running it.

There is no reason to follow VGP workflows, as they are very very outdated, as pointed out in comments.
> 2 votes

---
Tags: assembly, genome-sequencing
---


title: Sequence alignment with MUSCLE

text:
# Question
Title: Sequence alignment with MUSCLE
For protein multiple sequence alignment (MSA), I am using MUSCLE (v5.1). The default setting without any specifications gives a 'bad' output even though my proteins are very similar (UniProt sequence and the corresponding PDB SEQRES sequences). 'Bad' is defined as unnecessary gaps are introduced into the sequences.

However when the super5 algorithm is used the MSA output is correct, which is described here https://drive5.com/muscle5/manual/cmd\_super5.html

In addition, when I increase the gap penalty this issue can be solved but as far as I know, super4 v5.1 doesn't give you the opportunity to adjust gap penalty.

Documentation describes super5 quite similarly than the default algorithm, it "just" reduces runtime for long sequences.

Can someone please give me a hint what am I missing and/or why super5 outperforms default in terms of accuracy in such "easy" cases?

# Answer
**HMM** Super5 is a hidden Markov model (HMM) where gap penalties and the substitution matrix are the parameters at the centre of the model. That is the answer, thats why you can't change it - it does it itself and thats its purpose. Thus one algorithm (super5) is HMM for the gap-penalty and the other one I suspect isn't.

---

**Explanation** To explain that further HMM uses observed outcomes to estimate/optimise the hidden parameters, the gap penalty being one of the two key hidden parameters. Thus the algorithm is exploring a parameter landscape through reiteration generating multi-multiple alignments (if that makes sense) using different gap penalties and considering the error between the alignments (highlighted below). The search strategy is key to the success of the output. Muscle 5 describes 'introducing variations (perturbations) into HMM parameters', I assume this is singularly with reference to super5. The parameter landscape might be quite small because it doesn't seem to be a directional search strategy in any way (there's no memory in it).

The optimisation algorithm, i.e. what is being observed is a metric called the 'pair-wise error (PWE) and pair-wise dispersion (PWD)', so what appears to be happening is the algorithm is attempting to minimise the number of gaps in an alignment and this would be compatible with your results.

**Speed** Super5 is breaking the alignment up into smaller alignment and solving each 'centroid' before bringing together the 'centroids'. This is an obvious strategy particularly in the era of multiple processors because it overcomes RAM bottlenecks, we are all doing that in any algorithmic context.

**Overall** The algorithm was not absolutely clear where PWE and the feedback loop is being deployed whether it was in the generation of centroids or bringing together the centroids and its relationship to the guide tree. There was a large reliance on BLOSSUM62 - which didn't quite make sense in the given context because the background was clear about the different substitution matrices and the impact on gap penalties.

**Summary** Static gap penalty parameterisation is antithesis of what super5 is about, that is why you cannot use them. The speed is clear its using parallel processing and subalignments - there are no surprises there, thats a common strategy and easy to code these days.

**Final points**

* I am clear the ms has not been accepted and is available here https://www.biorxiv.org/content/10.1101/2021.06.20.449169v1.full#ref-10 . I don't want to comment further on the manuscript for several reasons, which I prefer to avoid explaining.
* It is better to wait until the paper is published because it will likely be different.
* The algorithm itself will unquestionably be subject to simulation studies by other authors.
* Simulation studies were absent in the BioRXIV manuscript.
* Muscle historically has been at the forefront of methods to increase sensitivity. If you read the Clustal Omega publication they describe the implementation of historical Muscle sensitivity improvements, in other words they must have been approved over the course of time.

Other authors are likely to investigate whether simulation conditions will result in super5 causing a 'false peak' scenario. That is the Achilles heal in any HMM.

There are certainly situations where super5 is useful in automation, that is without question. In some pipelines there is no option of checking the alignments within the pipeline because there are just to many and it stops the pipeline (human eyes the opposite of a data pipeline).

In practice optical annotation is used to overcome limitations post alignment, i.e. to resolve 'errors' in gap penalties and this is desirable because it prevents a phenomena of 'circular mathematics' which is a caution in any alignment method. The author does discuss this phenomena without calling this directly 'circular mathematics'. Any alignment algorithm can erroneously impose itself on the alignment and the metric of 'error' is extremely important in this regard. Optical analysis breaks circularisation. There is a very, very obvious solution to all this without using HMM at one extreme or 'optical analysis' at the other extreme, but one day at a time.

**Example of HMM** 'Circularisation' is the Achilles heal of my explanations, because I explain stuff, then reexplain it which more likely increases confusion and intersperse everything with waffle. However, I do follow a HMM style algorithm because next time I explain this it will be more honed and punchier. The hidden parameter being 'clarity' and the observed phenomena being the frequency of blank expressions, which is sought to minimise on each iteration. I do hope my search of the parameter landscape is a bit more directional than that used in Super5.

---

parameters looking at the documentation

```
-perturb SEED

```

It is important to do the same alignment using different random number streams. A given seed should result in the same gap pattern. If there are different gap patterns between different 'SEEDS', there's something suboptimal in the HMM and that is not good at all. If you think the Muscle 5 is implicit perfection, you don't need to check it - but its a trust thing. Its good the author made it explicit, i.e. you can set it.

The other parameters are:

```
-perm PERM
-consiters N
-refineiters N

```

The second most important parameter in my view is `-perm` \- this is permutations on the UPGMA guide tree. If different `-perm` values give different alignments thats not good. It must work reasonably well because the author will have checked it.
> 1 votes

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


title: Accessing .bam/.cram files from AWS S3?

text:
# Question
Title: Accessing .bam/.cram files from AWS S3?
What's the best/easiest way for accessing .bam/.cram files from S3?

I have .bam/.bai .cram/.crai and .bed + .gff files that I want to make available easily to others so they can browse in IGV and other solutions from their computers. There seem to be different solutions after cursory Googling, what would you recommend?

Here is what I've found so far:

IGV:

Samtools: samtools directly supports AWS

```
The full URL format is:

    s3[+SCHEME]://[ID[:SECRET[:TOKEN]]@]BUCKET/PATH

```

# Answer
> 3 votes
You could set up a public read-only bucket in an S3 account and put your files there. URLs would look like:

```
https://${username}.s3.amazonaws.com/public/${filename}

```

For samtools, the request for a BAM file's index should be built off that URL. For IGV, I think you can specify separate, custom URLs for the data and the index resources.

Public buckets are available to all, without need for credentials. So that might be an issue to think about, with respect to billing for web requests and data transfer.

I'm still not 100% certain that samtools/tabix support https requests reliably (I have a prior question on here about that), but that might be one approach.

# Answer
> 1 votes
Doo you have your heart set on AWS, or are you open to using Google Cloud Platform (GCP)?

If you're open to using GCP, there are some advantages. For example GATK can read and write directly to and from Google storage buckets, and IGV can use your login credentials to access data in Google storage as described here:

https://googlegenomics.readthedocs.io/en/latest/use\_cases/browse\_genomic\_data/igv.html

Note that if you make the data public, by default you're going to get billed if people download it. Not sure how that's handled on AWS, but on GCP you can set the bucket (or just the file) to be "requester-pays", ie the person who downloads is the person who pays for it:

https://cloud.google.com/storage/docs/requester-pays

In general there are several ways to interact with Google buckets.

*Disclosure: I work at the Broad Institute, which uses a lot of GCP infrastructure; I do not get any compensation from Google Cloud, but I recently wrote a book about Genomics in the Cloud that uses GCP.*

# Answer
> 1 votes
I have recently open sourced a Java NIO SPI for S3 that allows direct reads from S3 using Java applications that use the NIO model. This allows applications like GATK and some parts of Picard to understand `s3://` URIs.

It is used as a drop in so as long as the jar is on your classpath or in your JRE `ext` folder then Java will delegate nio operations to this lib when it "sees" an `s3://` URI. It will use the credentials and configuration of your current profile so no additional config is require and no recompilation of an application is needed.

The github repo is https://github.com/awslabs/aws-java-nio-spi-for-s3

*Disclosure: I am the author of the library, employed by the AWS Health AI Genomics Team, and have financial interests in Amazon*

---
Tags: sam, samtools, igv, cram, aws
---


# Question
Title: Running into an error when loading Flexaid on PyMol 3.7
I have attached a screenshot from terminal, which is the result of me attempting to run Flexaid on my Mac.

I was wondering if anyone knows why this isn't running?

Thanks

# Answer
> 1 votes
*Answer from @maximilian-press, converted from comment. Please improve where possible:*

It seems that many of the files are not being found, including both run executables and your pdb files. That's what all the `No such file or directory` stuff is about.

I would suggest confirming that you have completely installed PyMol and any of the relevant dependencies such as `Process_ligand` that are not being found. Here is one resource that might be relevant for that tool. Is `run_flexaid.sh` a script that you have written or a third-party tool? I can't find it by google. Can you show us how you ran it?

I am guessing that `sh ./run_flexaid.sh` is the terminal command that you are using to run the script, but I still don't know where the script comes from. It's possible that running within the script is not using your environment properly somehow. What happens if you try to run `Process_ligand` on its own (if that's possible?). Have you confirmed that your input pdb files are in the right directory? Basically, it's hard to know what's going on without seeing the script and knowing where your input files live on your computer.

*Note: Comments are **temporary** "Post-It" notes left on a question or answer. They should be used for clarification or constructive criticism of questions or answers, not for answering (no matter how minor). More details on this StackExchange policy here.*

---
Tags: proteins, docking, pymol
---


# Question
Title: How to identify [hyper]thermophilic proteins by BLAST search?
I need to find homologs of certain proteins by BLAST search or similar. But I want filter (select) homologs from organisms that are classified as thermophilic or hyperthermophilic. I am wondering if there are suggestions for me how approach to this. Thank you.

# Answer
It seems like the main issue is finding a list of species to filter your blast results on. I will therefore address this list of species issue.

There are some resources such as Microbe Directory that contain organismal characteristics similar to what you describe. Their main site appears nonfunctional but their github still works and appears to have the full database. This database includes an `extreme_environment` field, which appears to be what you are going for.

I feel like PATRIC used to have something similar but I can't locate it now.

For more ideas you might look at the Microbe Directory paper to see if they cite other similar databases.
> 3 votes

# Answer
> Warning: Hearsay

As the others have said, filtering by a list of taxonomy ids is the only one. However, there is one extra consideration if this is for wet lab applications: expression yields even when codon optimised may be affected by phylum choice.

Archaeal protein encoded by codon-optimised genes for *E. coli* generally perform less well than the equivalent bacterial thermophilic protein. I have not read a study on it, but I have heard its tale over and over and had issues myself. GC-content and the amino acid composition do correlate, but I have heard of good results in *E. coli* with *Geobacillus* (high-GC thermophilic firmicute) as with *Bacillus* (high GC mesophilic firmicute) after all. I don't know if its isoelectric points close to pH 7 or what. But mid GC-content bacterial thermophiles, such as *Aquifex* or some clades of *Betaproteobacteria*, may be a better call than *Archaea*.

So as you need to filter your blast search by a series of phyla, you may want to be selective if for expression in a different host...
> 1 votes

# Answer
I recommend use the local BLAST search. You can form the request for desired taxa in Entrez or through NCBI website & download fasta-file onto your disk.
> 0 votes

---
Tags: proteins, blast, sequence-homology
---


title: Bowtie2 mapping with multiple indexes

text:
# Question
Title: Bowtie2 mapping with multiple indexes
**Background:**

We performed NGS using cells collected from mice in a xenotransplantation study.

As such, the FASTQ files contain reads of DNA from both mice and human cells.

I expect ~30% of reads are contaminated by mice cells.

**Question:**

Is it wrong to directly align using the indexed human reference with bowtie2?

Is there a way to use two indexes at once with bowtie2? The goal here would be to get the % of reads that map to each genome.

# Answer
I am not aware of a method using two indices in `bowtie2` but here is a simple workaround:

Get human reference genome as fasta and suffix all fasta names with `_human`. Do the same with the mouse genome using `_mouse`. `cat` both together and build an index. Then you can later track back whether the alignment was done to human or mouse.

Edit: Be aware though that for those parts of the genomes that have high sequence similarity you will get multimapping reads with MAPQ of 0 (or at least very low).
> 3 votes

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