# Question
Title: `TxDb.Hsapiens.UCSC.hg19.knownGene` and `EnsDb.Hsapiens.v75` can't be installed
I think there is a problem with `TxDb.Hsapiens.UCSC.hg19.knownGene` and `EnsDb.Hsapiens.v75`. I'm struggling to install all of these packages as they pause at `byte-compile and prepare package for lazy loading` . Can anyone try to re-install them in their local machine because am not the only one who is having the problem, please? They used to work fine until a week ago!

As usual, I install the packages by:

```
BiocManager::install(c("Homo.sapiens", "EnsDb.Hsapiens.v75", "TxDb.Hsapiens.UCSC.hg19.knownGene"))

```

I also tried installing them individually:

```
BiocManager::install("EnsDb.Hsapiens.v75")
BiocManager::install("TxDb.Hsapiens.UCSC.hg19.knownGene")

```

It pause at:

```
* installing *source* package ‘TxDb.Hsapiens.UCSC.hg19.knownGene’ ...
** using staged installation
** R
** inst
** byte-compile and prepare package for lazy loading

```

sessionInfo( )

```
R version 4.0.3 (2020-10-10)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur 10.16

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib

locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

loaded via a namespace (and not attached):
[1] BiocManager_1.30.10 compiler_4.0.3      tools_4.0.3

```

I've tried with and without RStudio. I'd be grateful if someone can try to replicate the error, please?

# Answer
I am able to successfully install these packages using:

```
BiocManager::install(c("Homo.sapiens", "EnsDb.Hsapiens.v75", "TxDb.Hsapiens.UCSC.hg19.knownGene"))

```

This completes without any errors / warnings (R v4.0.3 / macOS 10.16), so I don't believe the packages themselves/Bioconductor are 'broken', but I suspect you are having compiler issues due to the Big Sur update.

These are the steps I used to solve the problem:

1. Reinstall xcode command line tools (even if it says "up to date")

```
sudo rm -rf /Library/Developer/CommandLineTools
sudo xcode-select --install

```
2. Install gcc & llvm via Homebrew (instructions for installing Homebrew) or, if you already have gcc/llvm installed via Homebrew, skip to the next step

```
# This can take several hours
brew install gcc
brew install llvm

```
3. Once you have gcc & llvm installed via Homebrew:

```
brew cleanup
brew update
brew upgrade
brew reinstall gcc
brew reinstall llvm

```
4. Link some headers into /usr/local/include

```
sudo ln -s /Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include/* /usr/local/include/

# You can safely ignore warnings like this:
#ln: /usr/local/include//tcl.h: File exists
#ln: /usr/local/include//tclDecls.h: File exists
#ln: /usr/local/include//tclPlatDecls.h: File exists
#ln: /usr/local/include//tclTomMath.h: File exists
#ln: /usr/local/include//tclTomMathDecls.h: File exists
#ln: /usr/local/include//tk.h: File exists
#ln: /usr/local/include//tkDecls.h: File exists
#ln: /usr/local/include//tkPlatDecls.h: File exists

```
5. Create a new `~/.R/Makevars` file (if you already have a `~/.R/Makevars` file, save it in a different directory (away from `~/.R/`)) and include only these lines in the file:

```
FLIBS=-L/usr/local/gfortran/lib/gcc/x86_64-apple-darwin19/10.2.0 -L/usr/local/gfortran/lib -lgfortran -lquadmath -lm
CXX1X=/usr/local/gfortran/bin/g++
CXX98=/usr/local/gfortran/bin/g++
CXX11=/usr/local/gfortran/bin/g++
CXX14=/usr/local/gfortran/bin/g++
CXX17=/usr/local/gfortran/bin/g++

LLVM_LOC = /usr/local/opt/llvm
CC=/usr/local/gfortran/bin/gcc -fopenmp
CXX=/usr/local/gfortran/bin/g++ -fopenmp
CFLAGS=-g -O3 -Wall -pedantic -std=gnu99 -mtune=native -pipe
CXXFLAGS=-g -O3 -Wall -pedantic -std=c++11 -mtune=native -pipe
LDFLAGS=-L/usr/local/opt/gettext/lib -L$(LLVM_LOC)/lib -Wl,-rpath,$(LLVM_LOC)/lib
CPPFLAGS=-I/usr/local/opt/gettext/include -I$(LLVM_LOC)/include -I/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include

```
6. Compile a package from source in R/Rstudio

```
# To check whether openmp is enabled, compile data.table:
install.packages("data.table", type = "source")

```
7. If your package fails to compile, a couple of SO users have had to install a fresh gfortran (re: https://stackoverflow.com/a/65334247/12957340), which you can download from https://github.com/fxcoudert/gfortran-for-macOS/releases/tag/10.2-bigsur-intel

Once you can install packages from source without error/issue, update Bioconductor

```
BiocManager::install(version = "3.12")

```

Then, finally, you should be able to install your packages using

```
BiocManager::install(c("Homo.sapiens", "EnsDb.Hsapiens.v75", "TxDb.Hsapiens.UCSC.hg19.knownGene"))

```
> 1 votes

---
Tags: bioconductor, software-installation
---


title: VlnPlot problem

text:
# Question
Title: VlnPlot problem
I was confused about the VlnPlot that it shows like the screenshot. It only shows the dot but no distribution of the data, is that normal?

# Answer
> 6 votes
I assume you are referring to `VlnPlot()` of Seurat. The reason that you are getting such a plot is because your distribution is highly skewed, most of your points are zeros, not surprising for scRNA-seq. And the plot is actually showing your distributions, it is just the points at zero stacking on top of each other.

You can double check this by plotting a histogram for each of the clusters.

An example violin plot:

Cluster 0, which is mostly zeros, kind of looks like yours.

Cluster 6, which I believe is what you would expect to see, contains a lot of cells that express the gene and much less skewed, hence the shape denoting the density.

So what you are seeing is normal, your gene is not expressed/not detected is all.

# Answer
> 0 votes
Try using ridgeplots instead of violin plots. Ridge plots will show you that bump at zero

---
Tags: scrnaseq
---


# Question
Title: What changed wwith Python in WDL change from draft 2 to version 1?
I'm updating my WDL from draft 2 (whichs works fine) to version 1. I'm stuck with a Py snippet

```
task getReadGroup {

    input {
        String inputFile
    }

    command {echo input = ${inputFile}}

    # command <<<
    # python <<CODE
    # print ("${inputFile}".strip().split('/'))
    # CODE
    # >>>

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

```

This command gives the expected result:

```
input = /mnt/lustre/working/genomeinfo/data/20191015_LongRead/20191015_LongRead_1/data/LXBAB201271/LXBAB201271/LXBAB201271/20190417_0612_1-A7-D7_PAD44174_42dcde1a/fast5/test/

```

But the commented out python code:

```
print ("${inputFile}".strip().split('/'))

```

Gives an empty list - "${inputFile}" is empty in the code command for some reason!

Any ideas? I've tried using a File instead of String to make it just like the example in the docs (below) other than that I see no diff:

```
task heredoc {
  input {
    File in
  }

  command<<<
  python <<CODE
    with open("${in}") as fp:
      for line in fp:
        if not line.startswith('#'):
          print(line.strip())
  CODE
  >>>
}

```

# Answer
Complements of Phil Shapiro;

Use ~{} instead of ${}.

So this example from the current docs is wrong:

```
task heredoc {
  input {
    File in
  }

  command<<<
  python <<CODE
    with open("${in}") as fp:
      for line in fp:
        if not line.startswith('#'):
          print(line.strip())
  CODE
  >>>
}

```

It should be :

```
task heredoc {
  input {
    File in
  }

  command<<<
  python <<CODE
    with open("~{in}") as fp:
      for line in fp:
        if not line.startswith('#'):
          print(line.strip())
  CODE
  >>>
}

```
> 2 votes

---
Tags: python, wdl
---


# Question
Title: How to generate a consensus sequence from a multi-reference BAM file?
I am trying to generate a consensus sequence from a BAM file that was generated by mapping reads to a reference FASTA containing multiple sequences.

Usually, I generate consensus sequences from BAM files using samtools and bcftools:

```
samtools mpileup -vf reference.fasta filename.sorted.bam | bcftools call -m -O z - > filename.vcf.gz

bcftools index filename.vcf.gz

bcftools consensus -f reference.fasta filename.vcf.gz > filename.consensus.fasta

```

However, this doesn't work when the reference FASTA contains more than 1 sequence.

Is there a way to generate a consensus sequence from a multi-reference BAM file?

It would also be awesome if the method filled in gaps in coverage with Ns or Xs in the generated consensus sequence but this isn't essential.

Thank you for your help!

# Answer
This post on biostars has several suggestions.

THe one that seems most straightforward, starting where you are, is to use vcftools:

```
cat ref.fa | vcf-consensus file.vcf.gz > out.fa

```

Based on that I think that you should also be able to solve this problem using alfred and GATK.

If all of these suffer the same problem of needing a single sequence (which I very much doubt that they do), you could split the bam / ref files into multiple files, run whatever consensus workflow, and then combine them again. But that doesn't seem like it should be necessary.

See also this post and this post.
> 0 votes

---
Tags: sam, samtools, ngs, mpileup, bcftools
---


# Question
Title: MSA (protein) with biopython or something else?
I am very new to bioinformatics (and python in general), but I would like to use python to more efficiently analyse enzymes both in terms of structure and functio, using Jupiter notebook. I would like to ask what is the best program/source code for multiple sequence alignments (amino acids) to identify conserved binding sites etc. I see that biopython has a few ways of creating alignments, but I have the impression that it is more focused on nucleotide sequences?

# Answer
In my personal experience, MUSCLE is the easiest program to use in conjunction with Biopython. Biopython features a command line wrapper for this program, which makes it very easy to use. Make sure to download the appropriate MUSCLE program from drive5 and save it somewhere. E.g., if you are using Jupyter in Linux:

```
!wget https://www.drive5.com/muscle/downloads3.8.31/muscle3.8.31_i86linux64.tar.gz
!tar -xzvf muscle3.8.31_i86linux64.tar.gz
!cp muscle3.8.31_i86linux64 /usr/local/bin
!chmod 755 /usr/local/bin/muscle3.8.31_i86linux64

```

Then, you can run MUSCLE like so:

```
from Bio.Align.Applications import MuscleCommandline
def runMUSCLE(infile, outfile):
    muscle_exe = r"/usr/local/bin/muscle3.8.31_i86linux64" #Here is where we installed MUSCLE
    muscle_cline = MuscleCommandline(muscle_exe,
                                     input=infile,
                                     out=outfile,
                                     clwstrict=True #Output in clustal format (more visually pleasing), otherwise the output is in FASTA. Whichever you need.
                                    )
    muscle_cline()

```

You don't need to specify that your sequence is amino acidic, however keep in mind that the input file must be in FASTA format.
> 1 votes

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


# Question
Title: How to calculate various properties of a protein structure per atom using PDB2PQR and Rosetta tools (or any other tools)?
Currently, I am implementing a descriptor of protein structure and would like to calculate properties per atoms like charges, hydrophobicity, hydrogen bond donors / acceptors, hydrophobicity, aromaticity etc.

I tried using PDB2PQR and succeeded with protein structures but it crashed while calculating nucleic acid bound proteins and choosing other force fields did not help. Also, as far as I know, Rosetta tool can calculate hydrophobicity however, it assigns the property per residue but not per atom (correct me if I am wrong).

Any other suitable tool suggestions would be appreciated.

# Answer
## Per Atom

You are correct. Rosetta scorefunction does not store any per atom data. The scoring operates at the per residue level. Whereas each atom has its coordinates and properties in full atom mode, in Pyrosetta it is clear that an atom is a just part of a residue and every operation is applied at the residue level. It's a team effort: the functional group of the side chain works as one and in most cases two atoms have the same role (e.g. `OE1` and `OE2` in glutamate behave identically).

So, the real question is whether you really care about the atom level in a classical mechanical setting?

## Formal charges

A wee parenthesis, we need to make a difference between formal charge and Marsilli-Gasteiger partial charge —one is an integer and the other is a float. Both require a stated protonation state in most applications —say `aspartate` vs. `aspartic acid`.

`PDB2PQR` formats to a PQR file for the ligand docking application `autodock`. It does not assign how a residue is protonated. In the case of Rosetta, the partial charge depends on the 4 column of the topology (params file) of the amino acid, i.e. it is fixed for all amino acids —Rosetta does not implement a Drude particle to make something polarisable. I very strongly suspect this is the same in PDB2PQR.

Therefore, something like the Advance Poisson-Boltzmann Solver may give more interesting data in terms of partial charges, but also requires the uses to specify the protonation state.

A small parenthesis is histidine, which can adopt four forms. Two tautomers at neutral pH: single proton on either the epsilon (`HIE` in Amber) or delta nitrogen (`HID`), while it can be deprotonated (anion) or diprotonated (`HID`, cation). Rosetta switches between these two automatically.

## Properties

You list a few properties of interest, which all seem to be based on the residue a priori of structure.

For example:

* hydrophobicity, unless you have a polarisable model (Drude particle) this is structure independent. RSA, relative surface area, is a must better indication (per residue), which then can get classified as core or surface based on this. Per atom would require just a wee tinker.
* aromaticity, this is also structure independent, but π-π stacking or π-cation interactions are noted at the residue level.
* hydrogen bond acceptor/donor. Whether an atom is a donor or acceptor is actually defined in the atom type used in the topology file. Whether that atom fulfils that bonding is a different matter. If you really really wanted to get per atom values of this, in Rosetta you could exclude the backbone atoms you could mutate to alanine without minimising and see if it's maintained. So differentiating between backbone and sidechain is fine, but per atom is a different matter.

A structure-dependent property is how flexible the residue is, as measured as an RMS. This can be calculated per atom. This is akin to B-factors, which is how big is the "blur" of the atom. However, if you get the `NZ` atom on a lysine that is switching between two or more positions you get a partial occupancy. If this atom is basically hand-waving along an axis in an arc you get an anisotropic B-factor. In other words it does not vibrate in all 3 cartesian axes equally. So a B-factor is a terrible approximation. An RSMD of an ensemble is even worse. In Rosetta, you can a quick vibration using Backrub, which is very quick at telling you how much motion does an atom experience.

## Here be dragons

> Warning: As mentioned, it takes two atoms to resulting in a force, hence why per atom scores are a bad idea.

However if for some *exceptional* reason this is required, say figuring out which group on a ligand is causing a clash, in pyrosetta it is possible:

```
import pandas as pd
import pyrosetta
pyrosetta.init()
pose = pyrosetta.pose_from_sequence("ELVISISALIVE")
target_resi = 9
residue = pose.residue(target_resi)
scores = []
for i in range(1, residue.natoms() + 1):
    iname = residue.atom_name(i)
    for r in range(1, pose.total_residue() + 1):
        other = pose.residue(r)
        for o in range(1, other.natoms() + 1):
            oname = other.atom_name(o)
            score = pyrosetta.toolbox.atom_pair_energy.etable_atom_pair_energies(residue,
                                                                                i,
                                                                                other,
                                                                                o,
                                                                                scorefxn)
        scores.append(dict(zip(['residue_index_A', 'atom_name_A', 'residue_index_A', 'atom_name_B', 'lj_atr', 'lj_rep', 'fa_solv', 'fa_elec'],
                         [target_resi, residue.atom_name(i), r, other.atom_name(o), *score])))

df = pd.DataFrame(scores)
df.pivot_table(index='atom_name_A',
               aggfunc='sum')[['lj_atr','lj_rep','fa_elec','fa_solv']]

```
> 2 votes

---
Tags: pdb, mutations, rosetta
---


# Question
Title: How to do `bedtools intersection` using pandas alone?
I have two pandas Dataframes, using python3.x:

```
import pandas as pd

dict1 = {0:['chr1','chr1','chr1','chr1','chr2'],
    1:[1, 100, 150, 900, 1], 2:[100, 200, 500, 950, 100],
    3:['feature1', 'feature2', 'feature3', 'feature4', 'feature4'],
    4:[0, 0, 0, 0, 0], 5:['+','+','-','+','+']}

df1 = pd.DataFrame(dict1)

print(df1)

##       0    1    2         3  4  5
## 0  chr1    1  100  feature1  0  +
## 1  chr1  100  200  feature2  0  +
## 2  chr1  150  500  feature3  0  -
## 3  chr1  900  950  feature4  0  +
## 4  chr2    1  100  feature4  0  +

dict2 = {0:['chr1','chr1'], 1:[155, 800], 2:[200, 901],
    3:['feature5', 'feature6'], 4:[0, 0], 5:['-','+']}

df2 = pd.DataFrame(dict2)
print(df2)
##       0    1    2         3  4  5
## 0  chr1  155  200  feature5  0  -
## 1  chr1  800  901  feature6  0  +

```

The columns to focus on in these dataframes are the first three columns: location, start, and end. Each start:end value represents a distance on location (e.g. `chr1`, `chr2`, `chr3`).

I would like to output the intersection of `df1` against `df2`. Here is the correct output:

```
chr1    155 200 feature2    0   +
chr1    155 200 feature3    0   -
chr1    900 901 feature4    0   +

```

I was hoping there was a pandas only solution...

# Answer
Your setup:

```
import pandas as pd

dict1 = {0:['chr1','chr1','chr1','chr1','chr2'],
    1:[1, 100, 150, 900, 1], 2:[100, 200, 500, 950, 100],
    3:['feature1', 'feature2', 'feature3', 'feature4', 'feature4'],
    4:[0, 0, 0, 0, 0], 5:['+','+','-','+','+']}

df1 = pd.DataFrame(dict1)

print(df1)

##       0    1    2         3  4  5
## 0  chr1    1  100  feature1  0  +
## 1  chr1  100  200  feature2  0  +
## 2  chr1  150  500  feature3  0  -
## 3  chr1  900  950  feature4  0  +
## 4  chr2    1  100  feature4  0  +

dict2 = {0:['chr1','chr1'], 1:[155, 800], 2:[200, 901],
    3:['feature5', 'feature6'], 4:[0, 0], 5:['-','+']}

df2 = pd.DataFrame(dict2)

```

Now use pyranges (https://github.com/biocore-ntnu/pyranges):

```
# pip install pyranges
# conda install -c bioconda pyranges

import pyranges as pr
cols = "Chromosome Start End Name Score Strand".split()
# Add column names to the dfs so PyRanges knows which
# represent Chromosomes, Starts, Ends and optionally Strand

df1.columns = cols
df2.columns = cols

# create PyRanges-objects from the dfs
gr1, gr2 = pr.PyRanges(df1), pr.PyRanges(df2)

# intersect the two
gr = gr1.intersect(gr2)

# get the result as a   df
df = gr.df

```

Results:

```
print(gr)
# +--------------+-----------+-----------+------------+-----------+--------------+
# | Chromosome   |     Start |       End | Name       |     Score | Strand       |
# | (category)   |   (int32) |   (int32) | (object)   |   (int64) | (category)   |
# |--------------+-----------+-----------+------------+-----------+--------------|
# | chr1         |       155 |       200 | feature2   |         0 | +            |
# | chr1         |       900 |       901 | feature4   |         0 | +            |
# | chr1         |       155 |       200 | feature3   |         0 | -            |
# +--------------+-----------+-----------+------------+-----------+--------------+
# Stranded PyRanges object has 3 rows and 6 columns from 1 chromosomes.
# For printing, the PyRanges was sorted on Chromosome and Strand.
print(df)
#   Chromosome  Start  End      Name  Score Strand
# 0       chr1    155  200  feature2      0      +
# 1       chr1    900  901  feature4      0      +
# 2       chr1    155  200  feature3      0      -

```

See the docs for more info: https://biocore-ntnu.github.io/pyranges/
> 8 votes

# Answer
I don't think Pandas has this implemented functionality out-of-the-box. Even if it did, solutions not designed specifically for bioinformatics probably rarely handle intervals on different chromosomes correctly unless you split the intervals by chromosome first.

Pandas does *handle* intervals (see docs for the Interval and IntervalIndex classes), but I've never used these. I *have* used the excellent intervaltree package. Whatever approach and tools you decide to use, there's a good chance you'll have to write some code. Computing interval intersections with `intervaltree` would look something like this.

```
>>> from collections import defaultdict
>>> from intervaltree import IntervalTree
>>> from pandas import DataFrame
>>>
>>> dict1 = {
...     0:['chr1','chr1','chr1','chr1','chr2'],
...     1:[1, 100, 150, 900, 1],
...     2:[100, 200, 500, 950, 100],
...     3:['feature1', 'feature2', 'feature3', 'feature4', 'feature4'],
...     4:[0, 0, 0, 0, 0],
...     5:['+','+','-','+','+']
... }
>>> dict2 = {
...     0:['chr1','chr1'],
...     1:[155, 800],
...     2:[200, 901],
...     3:['feature5', 'feature6'],
...     4:[0, 0],
...     5:['-','+']
... }
>>> df1 = DataFrame(dict1)
>>> df2 = DataFrame(dict2)
>>>
>>>
>>> # Construct an interval index from one of the dataframes,
>>> #   one tree per chromosome
>>> index = defaultdict(IntervalTree)
>>> for n, row in df2.iterrows():
...     index[row[0]].addi(row[1], row[2], row)
...
>>>
>>> # Query the index from the other dataframe to build a 3rd
>>> #   dataframe representing the intersection
>>> dict3 = { i: list() for i in range(6) }
>>> for n, row in df1.iterrows():
...     overlapping = index[row[0]].overlap(row[1], row[2])
...     for itvl in overlapping:
...         itsct_start = max(row[1], itvl.begin)
...         itsct_end = min(row[2], itvl.end)
...         dict3[0].append(row[0])
...         dict3[1].append(itsct_start)
...         dict3[2].append(itsct_end)
...         dict3[3].append(row[3])
...         dict3[4].append(row[4])
...         dict3[5].append(row[5])
...
>>> df3 = DataFrame(dict3)
>>> df3
      0    1    2         3  4  5
0  chr1  155  200  feature2  0  +
1  chr1  155  200  feature3  0  -
2  chr1  900  901  feature4  0  +
>>>

```
> 4 votes

# Answer
As mentioned by OP, another option is to use pybedtools, which in my opinion is pretty convenient for people already familiar with BedTools.

Let's even say df1's format is slightly different than df2:

```
import pandas as pd

dict1 = {0: ['chr1', 'chr1', 'chr1', 'chr1', 'chr2'], 1: [1, 100, 150, 900, 1],
         2: [100, 200, 500, 950, 100],
         3: ['feature1', 'feature2', 'feature3', 'feature4', 'feature4'],
         4: [0, 0, 0, 0, 0], 5: ['+', '+', '-', '+', '+'],
         6: ["remember", "the", "5th", "of", "november"]}

df1 = pd.DataFrame(dict1)

print(df1)

      0    1    2         3  4  5         6
0  chr1    1  100  feature1  0  +  remember
1  chr1  100  200  feature2  0  +       the
2  chr1  150  500  feature3  0  -       5th
3  chr1  900  950  feature4  0  +        of
4  chr2    1  100  feature4  0  +  November

dict2 = {0: ['chr1', 'chr1'], 1: [155, 800], 2: [200, 901],
         3: ['feature5', 'feature6'], 4: [0, 0], 5: ['-', '+']}

df2 = pd.DataFrame(dict2)

print(df2)

      0    1    2         3  4  5
0  chr1  155  200  feature5  0  -
1  chr1  800  901  feature6  0  +

```

Using pybedtools we get:

```
from pybedtools import BedTool

bed1 = BedTool.from_dataframe(df1)
bed2 = BedTool.from_dataframe(df2)

intersected_bed = bed1.intersect(bed2)

print(intersected_bed)

chr1    155 200 feature2    0   +   the
chr1    155 200 feature3    0   -   5th
chr1    900 901 feature4    0   +   of

```

We can then keep using the intersected BedTool object or convert it back into a DataFrame:

```
# create a pandas.DataFrame, passing args and kwargs to pandas.read_table
intersected_df = intersected_bed.to_dataframe(header=None,
                                              names=["chrom", "start", "end", "name",
                                                     "score", "strand", "whatever"])

print(intersected_df)

  chrom  start  end      name  score strand whatever
0  chr1    155  200  feature2      0      +      the
1  chr1    155  200  feature3      0      -      5th
2  chr1    900  901  feature4      0      +       of

```
> 3 votes

# Answer
I implemented pandas "Intervals" and ... it should be a few lines, clearly there are limitations. For non-overlapping data it is very cool however. It will work for overlapping data, BUT if the data you are using as the interval data is overlapping, it falls over. It could work if an independent (non-overlapping) interval was constructed.

Anyway, the point of the excercise was to keep it all in pandas, "start\_stop" should not be needed. This can't be done via Intervals, you've got to drop out, and that will slow it down and harder to parse.

The query contig `df3.iloc[1,]`, is still a bug (hence not printed), should be \[x,\] where x is the value\_counts hit, but this is just a pilot.

The script would be fine if it the df was wrangled using `groupby` rather than `value_counts` because `value_counts` trashes important data, which is what I'm trying to recover via `df3.iloc[1,]`. I like pandas elegant solutions, but you know another day another dollar.

Note: I changed 200 to 199 otherwise there's too much overlap.

**Overlap Output** (not labeled very well)

```
                 start  stop
index
(1.0, 100.0]        0     0
(100.0, 199.0]      1     0
(150.0, 500.0]      1     1
(900.0, 950.0]      0     1

```

---

```
import pandas as pd

def tableit (dict):
    df = pd.DataFrame(dict)
    table = pd.pivot_table(df, values=['start', 'stop'], columns=['chromosome'], index=['feature', 'misc', 'misc2'], fill_value = '-')
    return table

def start_stop(point, table2, bin):
    df3 = pd.cut(table2[point,'chr1'], bin[0]).value_counts()
    df3 = df3.to_frame(point).reset_index()
    for x in bin:
            if not (x == bin[0]):
                df3_tmp = pd.cut(table2[point,'chr1'], x).value_counts()
                df3_tmp = df3_tmp.to_frame(point).reset_index()
                df3 = pd.concat([df3, df3_tmp]
            else:
                continue
    df3.set_index('index', inplace=True)
    return df3, df3.iloc[1,] # [1,] is a bug and the query locus needs correctly reading from df3_tmp

def main():
    dict1 = {'chromosome':['chr1','chr1','chr1','chr1','chr2'], 'start':[1, 100, 150, 900, 1], 'stop':[100, 199, 500, 950, 100], 'feature':['feature1', 'feature2', 'feature3', 'feature4', 'feature4'], 'misc':[0, 0, 0, 0, 0], 'misc2':['+','+','-','+','+']}
    dict2 = {'chromosome':['chr1','chr1'], 'start':[155, 800], 'stop':[200, 901], 'feature':['feature5', 'feature6'], 'misc':[0, 0], 'misc2':['-','+']}

    table1 = tableit(dict1)
    table2 = tableit(dict2)

    bin = [[i, j] for i, j in zip(table1['start', 'chr1'], table1['stop', 'chr1'])]

    df_start, locus = start_stop('start', table2, bin)
# locus.name gets to the data
    df_stop, _ = start_stop('stop', table2, bin)
    df_merge = pd.merge(df_start, df_stop, on="index", how='outer')
    print (df_merge)

if __name__ == "__main__":
    main()

```

---

The following would have been better,

```
table_join = pd.concat([table1, table2]).fillna('-')

```

and then sort them.. perhaps including a new column to specify the original dataframe and screen for the overlap with table\_join.iloc\[i, column\] - table\_join.iloc\[1-i, other\_column\].

The key step in all approaches is to convert from "long format" to "wide format" and thats the most important thing I've done in all code and makes pandas work.

**Output Table\_join**

```
                     start        stop
chromosome            chr1 chr2   chr1 chr2
feature  misc misc2
feature1 0    +        1.0    -  100.0    -
feature2 0    +      100.0    -  199.0    -
feature3 0    -      150.0    -  500.0    -
feature4 0    +      900.0    1  950.0  100
feature5 0    -      155.0    -  200.0    -
feature6 0    +      800.0    -  901.0    -

```
> 2 votes

---
Tags: python, bedtools, interval, pybedtools, pandas
---


# Question
Title: Why GAG protein HIV-1 (virus) is in Staphylococcus sp. SS21 and in Klebsiella quasipneumoniae (Bacteria)?
I am one ignorant...

Why GAG protein HIV-1 (virus) is in Staphylococcus sp. SS21 and in Klebsiella quasipneumoniae (Bacteria)? These are the links in GenBank:

```
Staphylococcus sp. SS21
LOCUS       WP_193573829             503 aa            linear   BCT 29-OCT-2020
DEFINITION  hypothetical protein [Staphylococcus sp. SS21].

HIV-1
LOCUS       KU749412                9188 bp    RNA     linear   VRL 08-DEC-2016
DEFINITION  HIV-1 isolate DEMC14PK009 from Pakistan, partial genome.

```

protein align:

```
mgarasllsggqldrwerirlrpggkkkyrlkhivwasrelerfavnpgllessegcrkileqlqpslqtgseelrslyntvavlycvqqridvkdtkealekieeeqnkskkkaqqtaadkgnnspvsqnfpivqnlqgqmvhqpisprtlnawvkvveekafspevipmfsalsegatpqdlntmlntvgghqaamqmlketineeaaewdrlhpvhagpvapgqmreprgsdiagttstlqeqigwmtnnpaipvgeiykrwiilglnkivrmysptsildirqgpkepfrdyvdrfyktlraeqasqdvknwmtetllvqnanpdcktilkalgpaatleemmtacqgvggpshkarilaeamsqvtqptnilmqrgnfrtprktvkcfncgkeghiakncraprkkgcwkcgkeghqmkdcserqanflgkiwpshkgrpgnflqnrpeptappappeesfrfgegtttpaqkqepidkdlyplaslkslfgndpssq
|||||| | |  || |||||||||||| | ||| ||||||| ||| ||||||  |||  |  |||| |||| |||||| |||| ||||   | |||||||| ||||||||   || |  |        | | ||||||||||||| |||||||||||| ||||||||||||| ||||||||||||||||||||||||||||| |||||||||||||||||||||||||||||||||||||||| ||| | | ||  ||||||||||||||||||||||| |||||||||||||||||||| |||||||| |||||||| ||||||||||||||| |||| |||||||||||||||| |||| ||||||| |      ||| ||     | |||||||||| | ||||||||||||||||||||||| ||||||||||||||||||||||| |||||||||   ||||| | || || |||  ||   || ||||||| || ||
MGARASILRGEKLDTWERIRLRPGGKKHYMLKHLVWASRELDRFALNPGLLETLEGCKQIIKQLQPALQTGTEELRSLFNTVATLYCVHAGIPVKDTKEALDKIEEEQNKCQQKAKQ--AEAAAAGKISLNYPIVQNLQGQMVHQAISPRTLNAWVKVIEEKAFSPEVIPMFTALSEGATPQDLNTMLNTVGGHQAAMQMLKDTINEEAAEWDRLHPVHAGPVAPGQMREPRGSDIAGTTSTLAEQIAWITGNPPVPVGEIYKRWIILGLNKIVRMYSPXSILDIRQGPKEPFRDYVDRFFKTLRAEQATQDVKNWMTDTLLVQNANPDCKTILRALGPGATLEEMMTACQGVGGPGHKARVLAEAMSQTT--NAIMMQRSNFKGPKRTIKCFNCGKEGHLARNCRAPRKKGCWKCGKEGHQMKDCNERQANFLGKIWPSHKGRPGNFLQSRPEPTAPPA---ESFRF-EETT-PAPKQEQKDK--EPLTSLKSLFGSDPLSQ

```

The links in GenBank of Klebsiella quasipneumoniae :

```
LOCUS       WP_199660587             458 aa            linear   BCT 06-JAN-2021
DEFINITION  hypothetical protein, partial [Klebsiella quasipneumoniae].
ACCESSION   WP_199660587
VERSION     WP_199660587.1

```

protein align:

```
-----------------------grkkyk-lkhlvwasrelerfalnpgllesaegcqqimeqlqsalgtgseelkslfntvatlwcvhqrieikdtkealdkveevkskqkaqqaqrtqqaaaatgsssqnypivqnaqgqmvhqsmsprtlnawvkvveekafspevipmftalsegatpqdlnmmlnivgghqaamqmlkdtineeaaewdrthpvhagpvapgqireprgsdiagttstlqeqigwmtnnppipvgeiykrwivlglnkivrmyspvsildirqgpkepfrdyvdrffktlraeqatqevknwmtesllvqnanpdcktilrglgpgatleemmtacqgvggpghkarvlaeamsqaqqsnimmqrgnfrg-qrtikcfncgkeghiarnckaprkrgcwkcgkeghqmkdcnerqanflgriwpshkgrpgnfpqsrpeptappaesf-geevtpsprqe-----pvd--kel-----ypp
                       | ||   ||||||||||| ||||||||||  ||| ||  ||| || || ||| ||||||||| |||  |  ||||||||  ||   |   | |   |  ||| |  | ||||||| |||||||  ||||||||||| |||||||||||||||||||||||||| ||| |||||||||||||||||||||||| |||||||||||| ||||||||||||||| ||| | | ||| |||||||||| |||||||||||| ||||||||||||||||||||||||||||||||||||||| ||||||||||||||| |||||||||||||||||||||||||||||||||     ||||| || |  ||||||||||||| |||| |||| |||||||||||||||||||||||| |||||||||||| ||||||||||||||  || || | ||     |    | |
MGARASILRGEKLDTWERIRLRPGGKKHYMLKHLVWASRELDRFALNPGLLETLEGCKQIIKQLQPALQTGTEELRSLFNTVATLYCVHAGIPVKDTKEALDKIEEEQNK-CQQKAK--QAEAAAAGKISLNYPIVQNLQGQMVHQAISPRTLNAWVKVIEEKAFSPEVIPMFTALSEGATPQDLNTMLNTVGGHQAAMQMLKDTINEEAAEWDRLHPVHAGPVAPGQMREPRGSDIAGTTSTLAEQIAWITGNPPVPVGEIYKRWIILGLNKIVRMYSPXSILDIRQGPKEPFRDYVDRFFKTLRAEQATQDVKNWMTDTLLVQNANPDCKTILRALGPGATLEEMMTACQGVGGPGHKARVLAEAMSQ-TTNAIMMQRSNFKGPKRTIKCFNCGKEGHLARNCRAPRKKGCWKCGKEGHQMKDCNERQANFLGKIWPSHKGRPGNFLQSRPEPTAPPAESFRFEETTPAPKQEQKDKEPLTSLKSLFGSDPLSQ

```

The links in GenBank of hypothetical protein HMV63\_24500 \[Escherichia coli\] and protein align with HIV-1 GAG:

```
LOCUS       HAJ5901352               246 aa            linear   BCT 11-MAY-2020
DEFINITION  TPA: hypothetical protein HMV63_24500 [Escherichia coli].
ACCESSION   HAJ5901352
VERSION     HAJ5901352.1

-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------mtsnppipvgdiykrwiilglnkivrmyspvsildirqgpkepfrdyvdrffktlraeqatqdvknwmtetllvqnanpdcktilralgpgasleemmtacqgvggpshkarvlaeamsqannanimmqksnfkgprrnvkcfncgkeghiarncraprkrgcwkcgkeghqmkdcterqanflgkiwpshkgrpgnflqsrpeptappaesfkfeettptprqeprdrepltslkslfgsdplsq
                                                                                | ||| ||| ||||||||||||||||||| |||||||||||||||||||||||||||||||||||||| |||||||||||||||||||||| |||||||||||||| ||||||||||||  || |||| |||||| |  |||||||||||||||||||| ||||||||||||||| |||||||||||||||||||||||||||||||||||| |||||| | ||| | |||||||||||||||||
MGARASILRGEKLDTWERIRLRPGGKKHYMLKHLVWASRELDRFALNPGLLETLEGCKQIIKQLQPALQTGTEELRSLFNTVATLYCVHAGIPVKDTKEALDKIEEEQNKCQQKAKQAEAAAAGKISLNYPIVQNLQGQMVHQAISPRTLNAWVKVIEEKAFSPEVIPMFTALSEGATPQDLNTMLNTVGGHQAAMQMLKDTINEEAAEWDRLHPVHAGPVAPGQMREPRGSDIAGTTSTLAEQIAWITGNPPVPVGEIYKRWIILGLNKIVRMYSPXSILDIRQGPKEPFRDYVDRFFKTLRAEQATQDVKNWMTDTLLVQNANPDCKTILRALGPGATLEEMMTACQGVGGPGHKARVLAEAMSQTTNA-IMMQRSNFKGPKRTIKCFNCGKEGHLARNCRAPRKKGCWKCGKEGHQMKDCNERQANFLGKIWPSHKGRPGNFLQSRPEPTAPPAESFRFEETTPAPKQEQKDKEPLTSLKSLFGSDPLSQ

```

I remember Blast find result in **plant** bacteria or virus too

# Answer
> -3 votes
for me could be one lab escape...

it seems the page is for me this

```
  https://science.sciencemag.org/content/266/5193/1981

```

there is a pay wall so I can not read that... The title is

```
Crystal structure of the catalytic domain
of HIV-1 integrase:
similarity to other polynucleotidyl transferases

```

i know nothing of argument, but it seems someone unknow for unknow reason infect one bacteria of Klebsiella pneumonia with the virus HIV-1. Someone know what happen if one common bacteria infect with the HIV-1 is diffused? what happen if there is something activate virus replication in the bacteria when the bacteria is in some human?

the genbank file would be:

```
LOCUS WP_142766724 415 aa linear BCT 07-JAN-2020
DEFINITION DDE-type integrase/transposase/recombinase, partial
 [Klebsiella pneumoniae].

```

The general speech

```
We, You don't have the right (or the intelligence) to change
one iota of any DNA or RNA. Genetic material, if we still
want to live a few years, must always be read-only,
never re-written from us because consequences are unknows.

```

# Answer
> 2 votes
The interesting numbers of an alignment are coverage (91%) and identity (82.68%) of API70200.1 (the protein encoded in KU749412) against WP\_199660587.

There are two cases:

* technical error: samples were contaminated, which happens not infrequently.
* biology did something cool: horizontal gene transfer

The former is likely to be present when you have a single sequence where this gene is present and it is not from a complete genome.

The latter, unless it's an ancient transfer, would likely result in flanking regions and genes.

So blasting the gene to the Klebsiella quasipneumoniae proteome, shows there are three hits (two poor). There is no strain name in list, whereas there are three complete genomes of three strains. So the sequence WP\_199660587 is a one off and not from a complete genome. Therefore, the clues strongly point towards it being a technical error.

---
Tags: sequence-alignment, proteins
---


title: ambiguous calls in nanopolish

text:
# Question
Title: ambiguous calls in nanopolish
I am converting scores into frequencies (from nanopolish results) using Nanopore scripts. The script filters out lots of rows on the premise that reads where `LLR < T * N` are ambiguous.

where:

```
LLR = log-likelihood ratio (log-likelihood methylated / log-likelihood unmethylated)
T   = arbitrary threshold. Above that value the C is methylated.
N   = number of Cs (CpG) in the fragment analyzed.

```

What? Why? Where?

# Answer
Assuming that it is the postprocessing script `calculate_methylation_frequency.py` that you are using, user Wouter de Coster's comment is correct that the methylation event is filtered out at the frequency estimation stage if it has a low LLR. It will have a low LLR if it is considered not confident enough to report. This means not only `LLR > T * N`, but also the statistical threshold `LLR > 0` (see below).

Specifically, here is the code in the script:

```

def update_call_stats(key, num_called_cpg_sites, is_methylated, sequence):
    if key not in sites:
        sites[key] = SiteStats(num_called_cpg_sites, sequence)

    sites[key].num_reads += 1
    sites[key].called_sites += num_called_cpg_sites
    if is_methylated > 0:
        sites[key].called_sites_methylated += num_called_cpg_sites

```

That is the function that decides whether to add a site to output, which depends on `is_methylated`, which is set here:

```
    # [skip several lines of code]
        num_sites = int(record['num_motifs'])
        llr = float(record['log_lik_ratio'])

        # Skip ambiguous call
        if abs(llr) < args.call_threshold * num_sites:
            # [MP NOTE: SKIP if LLR < T * N]
            continue
        sequence = record['sequence']

        # [MP NOTE: SKIP IF LLR <= 0]
        is_methylated = llr > 0

        # [skip several lines of code]

        update_call_stats(key, num_sites, is_methylated, sequence)

```

It's a little messy to read because the `sites` dictionary is a global outside of `update_call_stats()`.

So overall, it seems like **if you want sites included that are below this threshold, you will have to either:**

1. use an earlier stage of processing (output of `nanopolish call-methylation` I believe according to the docs) or
2. alter the script to change the heuristic thresholds.
> 1 votes

---
Tags: nanopore, methylation
---


# Question
Title: How do I create a VCF file of all known pathogenic mutations in a gene of interest?
I would like to create a list (probably .vcf format would be good) of all known pathogenic missense mutations in a human gene of interest and then add other variants that could lead to the same pathogenic amino acid substitutions. I intend to search targeted sequencing data of my gene of interest for variants that would cause these particular substitions.

I'm very new to bioinformatics so I'm not so familiar with the appropriate databases and tools to use, but this is a rough idea of how it could be done:

1. extract all mutations listed in Ensembl within my gene of interest that are labelled with `HGMD_MUTATION` (is there an option to download a vcf file with all HGMD mutations within my gene? I've only managed to do so for individual mutations.)
2. annotate list with reference and mutant amino acids for every mutation in the list
3. add additional mutations (which may not be present in HGMD) which could lead to the same substitution. E.g. I'm interested in the substitution Proline (reference codon: CCG) to Leucine (TTR or CTN) at AA position 102, so I would like to list all hypothetical mutations that would create a leucine codon in this position, not only the one given in Ensembl (which is CCG -\> CTG).

Could you point me towards the best tools to use for these purposes? Unfortunately I haven't got a HGMD subscription. Many thanks!

# Answer
> 1 votes
I'm not sure how to generate the additional mutations, but I would say that HGMD is not the way to find all the pathogenic variants. I would probably filter this table by either Clinical significance or Evidence-\>Phenotype association.

---
Tags: vcf, genomics, ensembl, mutations
---


# Question
Title: Getting two class of promoters : ATAC seq data
In this paper they classify the promoters into two categories, highlighted in italics below.

> **Two Types of Promoters**
>
> Mammalian promoters can initiate transcription across broad or narrow genomic regions, and these promoter shapes, broad or narrow (sharp), correlate with distinct transcriptional regulatory mechanisms.33 We evaluated cardiac promoters predicted from CAGE clusters for these 2 major types of promoters by calculating the interquartile range of promoter CAGE clusters by determining the base pair distance between 10% and 90% of a promoter’s total signal. *We observed the expected 2 distinct populations, defined as sharp (interquartile range \<10 bp) and broad promoters (interquartile range ≥10 bp; Figure 2A). Broad promoters were those associated many different cellular functions, including housekeeping functions. In contrast, genes with sharp (narrow) promoters were those encoding proteins critical for tissue-specific functions seen by the presence myofilament, muscle tissue, and muscle contractions genes (Figure 2B). Thus, tissue-specific genes important for LV specification and function were more likely to have sharp promoters.*

Now My question is can I do this only based on ATAC seq data. If yes how do i get that information of narrow and broad promoters.

# Answer
Its highly unlikely that you will be able to categorize promoters on this basis using only ATAC data. This is a consequence of people using the terms "promoters" and "transcription start sites" interchangeably, when they are actaully different things. CAGE data measures transcription start sites which is the specific location where transcription starts. CAGE does not identify promoters, which are the complete region around the start of a gene that contains, at a minimum, all the sequences necessary to assemble the RNA polymerase II pre-initiation complex. Some definitions of promoter are even broader, including all regulatory sequence proximal the the gene start.

ATAC identifies regions of open chromatin, and I would expect the open chromatin regions around the starts of genes to be more like the broadest definition of promoter, and therefore to be significantly larger than the transcription starts sites they contain, whether they be broad or narrow TSSs.

You could of course look at a distribution of the peak sizes over TSSs and see if you found a bimodal distribution. If you did you could correlate the genes with broad/narrow TSS peaks with their known functions so see if it fitted the idea that narrows TSSs
> 1 votes

---
Tags: atac-seq
---


# Question
Title: Any way to quantify the variation of genes that expressed in Affymetrix expression data?
am experimenting preprocessed Affymetrix microarrays expression data matrix (Affymetrix probe-sets in rows (32830 probesets), and RNA samples in columns (735 samples)) for my downstream analysis. Here is how the data looks like:

```
> dim(eset_HTA20)
[1] 32830   735
## load assayData
HTA20_rma <- load("HTA20_RMA.RData")

> eset_HTA20[1:10,1:3]
             Tarca_001_P1A01 Tarca_003_P1A03 Tarca_004_P1A04
1_at                6.062215        6.125023        5.875502
10_at               3.796484        3.805305        3.450245
100_at              5.849338        6.191562        6.550525
1000_at             3.567779        3.452524        3.316134
10000_at            6.166815        5.678373        6.185059
100009613_at        4.443027        4.773199        4.393488
100009676_at        5.836522        6.143398        5.898364
10001_at            6.330018        5.601745        6.137984
10002_at            4.922339        4.711765        4.628124
10003_at            2.689344        2.771010        2.556756

```

**objective**:

I am interested in to reduce the dimension of the dataset by tossing off low expressed in genes in the experiment. To do so, I need to quantify the variation of the gene that expressed. I checked the possible approach and using the coefficient of variation (cv) could be one of the approaches I could try. It is not intuitive for me how to quantify the variation of genes that expressed.

Can anyone point me out how to lay out this procedure on the above expression data? Is there any Bioconductor package to do this kind of task? any idea to make this happen? thanks

# Answer
You can use the following code to calculate the coefficient of variation:

```
# expr is your expression matrix.
SD <- apply(expr, 1, sd)
CV <- sqrt(exp(SD^2) - 1)

```

It might be implemented in some package but it is so brief that you can write it again yourself. Then you can filter out those that are below certain percentage of the distribution of CV (like 0.1)

```
expr[CV > quantile(CV, probs = 0.1), ]

```
> 2 votes

# Answer
Another way to reduce dimensionality of your dataset is via principal component analysis or PCA. It is implemented as `pcaMethods` R package, for instance.
> 0 votes

---
Tags: r, microarray, gene-expression, limma
---


# Question
Title: RNA_Seq Analysis in R, propmapped function issue
I'm currently trying to learn RNA-seq data analysis and differential expression in R. I've been using the Combine lab tutorial.

Everything was running well until the tutorial asked to use the `propmapped` function (on this page) on my BAM files and then view the promapped proportion for each sample. In this tutorial we should be getting promapped proportion ranging from 0.05 to 0.150, however in my case I have no reads that align to my mapped genome (promapped = 0 in my case, for each sample). I don't know why as I've been using the data they provided for this tutorial.

I was wondering if someone could follow the first few steps of this tutorial until the propmapped step and tell me if they eventually get the right proportion of reads or not after the propmapped function. This could tell me if the tutorial is the culprit or not.

# Answer
> 1 votes
There must be a problem in the BAM files you are using. Here are steps I would do to figure where exactly is the problem.

1. Check the data you have downloaded. The tutorial is using 12 fastq.gz files (SRR1552444 - SRR1552455). You can check on sra how big they are supposed to be (For example SRR1552444 is supposed to be an ~1.7Gb file) and if your downloaded files are about the same. Check if the files are fine (In the terminal you can type `zcat data/SRR1552444.fastq.gz | head` to check the first few reads). If the files are 0 sized or somehow corrupted, you need to redownload them.
2. Check if your mapping reference looks alright (size wise, look at the first few lines) and if the reference is indexed (there is a `.fai` in the directory).
3. Check the size of the BAM files, they also need to be in the range of a few hundred Mb to a few Gb. I expect these will be corrupted for sure, otherwise, you would get your numbers. You might also want to check if they are indexed (although this should not be a problem, I suspect they will get indexed automatically). If there is nothing wrong the bamfiles, you might want to check another tool to summarise a bamfile (like `samtools flagstat`) to check if at the end the problem really is in `propmapped`. If that's the case.

This should at least help you to zoom to the core of your problem if not resolve it.

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


# Question
Title: What is "unmapped read segments" in the output of samtools idxstats?
`samtools idxstats` produces a four column output (see here)

```
The output is TAB-delimited with each line consisting of reference sequence name, sequence length,
# mapped read-segments and # unmapped read-segments. It is written to stdout. Note this may count reads
multiple times if they are mapped more than once or in multiple fragments.

```

What is the meaning of *unmapped reads* here? Naively, I would take it to mean all the reads that didn't map to a specific feature... however in a sample with millions of reads most of them should be counted as such, which is not the case in my calculations, where the number of unmapped reads is of the same order as the number of mapped reads.

I will very much appreciate explanations from people experience with `samtools`.

# Answer
An unmapped read segment is any entry with bit 4 of the flag field set (i.e., it's unmapped). Assuming the reference to which you aligned reasonably matches the sample from which the reads were generated, less than ~10% of to total reads should be unmapped.

Please note that BAM files have no idea what a "feature" is. They only contain information about reads and the chromosomes/contigs to which they align.
> 2 votes

---
Tags: samtools, read-mapping, reads, featurecounts
---


# Question
Title: How to extract reads with INDELs > a given size?
I'm trying to modify this https://www.biostars.org/p/253774/

To get reads with deletions \> 20bp

I think this gives reads with exactly 20bp dels:

```
samtools view -h in.bam chr1:12340-200000 \ | awk '$1 ~ "^@" || $6 ! "20D"' \ | samtools view -b - > out.bam

```

is it possibe to do \>?

Thanks!

# Answer
> 2 votes
There's no simple way to do that with awk, since you would need to parse the CIGAR string and iterate over the operations. In python using pysam you would iterate over the reads and filter as follows:

```
keep = False
for OP, LEN in b.cigartuples:
    if OP == 2 and LEN > 20:
        keep = True
        break

if keep:
    ...do stuff...

```

---
Tags: sam, samtools, reads, awk, indel
---


# Question
Title: Multi Omics comparison RNA seq : Condition specific gene signature
The authors of an article are comparing different cell type and the organ type and they show the ontology enrichment of the RNA-seq:

> Cell-type-specific and organ-specific marker genes were identified using a two-step procedure. First, we performed pairwise differential expression analysis for each of the three type of structural cells, comparing them across organs. Second, based on the resulting pairwise comparisons, we identified for each organ those genes that were upregulated in each cell type compared to at least five other organs. Pairwise comparisons between organs were performed using the limma package in R, separately for each of the three structural cell types. Significantly differential genes were selected based on statistical significance (adjusted P \< 0.05), average expression (log2(CPM) \> −1) and sequencing coverage (median number of reads greater than 10 in the group with stronger signal). On the basis of these pairwise comparisons, we counted the total number of times each gene was upregulated in a specific organ compared to all other organs. Genes that were upregulated in comparison to five or more other organs were selected as marker genes of the corresponding organ.

Now how is this part done ? **"Genes that were upregulated in comparison to five or more other organs were selected as marker genes of the corresponding organ.\*"**

The figure for this results is this:

# Answer
> 1 votes
From their methods section:

> On the basis of these lists of marker genes and receptor–ligand pairs, we inferred potential cell–cell interactions for all pairs of one structural cell type and one haematopoietic immune cell type, quantifying the enrichment for known receptor–ligand pairs among all pairs of marker genes between the structural cell type and the haematopoietic cell type (Fig. 2a). First, we counted all pairs of marker genes for each pair of structural and immune cell types. Second, we calculated the fraction of these gene pairs that were annotated as receptor–ligand pairs. Third, we tested whether this fraction was greater than the fraction of annotated receptor–ligand pairs across all pairs of genes. Fisher’s exact test was used to obtain P values and odds ratios, as implemented by the function ‘fisher.test’ in R. P values were adjusted for multiple testing using the Benjamini–Hochberg method, as implemented the function ‘p.adjust’ in R. Finally, significantly enriched pairs (adjusted P \< 0.05) of structural and immune cell types were connected by edges to generate a graph of cell–cell interactions.

> In addition to testing for significant enrichment (as described in the previous paragraph), we calculated aggregated receptor and ligand gene signatures (Fig. 2c) based on a manually curated list of immune-related receptors and ligands (Supplementary Table 4). These gene signatures were derived by normalizing expression values for each gene and averaging across all genes in a given set of biologically related receptors and ligands.

The actual values plotted then seem to be from EnrichR, since they mention that elsewhere in the methods.

---
Tags: rna-seq, figure-explanation
---


title: Reads count in metagenomics

text:
# Question
Title: Reads count in metagenomics
**Background:** I am developing a pipeline for metagenomic studies of human gut microbiote. In particular, I am mapping the reads data originating from shotgun whole genome sequencing onto a gene catalogue (similar to IGC) and counting reads mapped on each feature/gene. Till now I have been using NGLess, which is a rather convenient tool. However, we have some questions about how exactly it counts reads and it might be not the best option for scaling up our pipeline.

**Question:**
I am looking for alternative tools for gene count. A brief search has shown that `htseq-count` and `featureCount` are two frequently used tools. However they mainly appear in RNAseq context, and I am not certain whether they could be used in my case. I would also appreciate comments on the possibility of using `samtools mpileup`.

# Answer
> 2 votes
I think that there is a good argument for using technologies adapted from RNA-seq for metagenomic abundance quantification. The tool that I am most familiar with (that I use for this purpose) is kallisto, which uses pseudoalignment instead and is therefore less compute intensive.

For some background on the conceptual linkages between metagenomic abundance and transcriptomic abundance you can see the kallisto metagenome abundance paper or the explanatory blog post.

I think that it is preferable to use a tool that is specifically adapted for metagenomics, as there are always some numerical issues that might have been accounted for by the authors. I think that for this purpose it is maybe best not to naively use tools like `mpileup`, though I have no specific objection to it. However, if you look around you can see various papers using e.g. Cufflinks and Cuffdiff for effectively the same purpose, and Cufflinks is definitely intended for RNA-seq.

Depending on your specific application, you might want to use a different tool like mOTUS2, if you are interested in estimating abundance of specific genomes or lineages in a metagenome mixture. This method is probably more accurate as it is based on averaging across the whole genome, but it is not very useful if what you are interested in is more e.g. functional profiling of the metagenome.

---
Tags: read-mapping, metagenome, featurecounts
---


# Question
Title: Getting System.IO.FileNotFoundException with MaxQuant running on Linux with mono
I am getting this error with MaxQuant 1.6.14 running with mono (5.14.0.177=he1b5a44\_1001).

> System.IO.FileNotFoundException: Could not find file "/home/ubuntu/run-dir/combined/combinedRunInfo"

And I not sure why I am getting this error. The process has read and write permissions to that location.

# Answer
...seems the mqpar file was not compatible with the maxquant version.

Make sure the `mqpar.xml` file shows the same version number as the MaxQuant binary. `MaxQuant` is not really made for the command line. So, the errors shown are sometimes quite misleading.
> 1 votes

---
Tags: maxquant
---


title: FASTQ file reads extraction paired ends

text:
# Question
Title: FASTQ file reads extraction paired ends
I have paired end reads and want to extract only those which have this sequence within them:

```
TGTATGTAAACTTCCGACTTCAACTGTA

```

I tried using `grep -B1 -A2 "TGTATGTAAACTTCCGACTTCAACTGTA" file_1.fq | grep -v "^--$"` and the same from the mate pair. After this I tried aligning them , but they dont align because the reads headers are different. Can someone help in extracting the mate pair of one read from the other file (2.fq), so that they can be aligned ? regards.

# Answer
> 3 votes
Your question can be rewritten as, "how do I get the union of reads in two files?" Since you have already subset the mates, the steps are:

1. Process each file, storing the list of reads names in each.
2. Take the union of these.
3. Iterate through the original/non-subset files and write only the reads with names in the union to the output.

The files will be both subset correctly and synced properly such that bowtie2 or any other aligner can correctly process them.

# Answer
> 1 votes
It's a little bit unclear what you're asking, but if you want to repair paired-end files that have become disordered (or the linking has been lost), then the repair function of bbtools should work. This will pair up any coincident read IDs, and create a separate file for any reads that have no pairs in the other file:

https://jgi.doe.gov/data-and-tools/bbtools/bb-tools-user-guide/repair-guide/

```
repair.sh in1=broken1.fq in2=broken2 out1=fixed1.fq out2=fixed2.fq outs=singletons.fq repair

```

If grep is used on one file, and repair used on that file in combination with the original file from the matching read, then the relevant read pairs should appear in the output:

```
repair.sh in1=filtered_1.fq in2=file_2 out1=fixed_1.fq out2=fixed_2.fq outs=unpaired_2.fq repair

```

---
Tags: grep
---


# Question
Title: RNASeq: Normalization, stabilization, gene length and rlog
I was thinking about the best method for normalization, which takes gene length into account (in order to compare genes)...

Do you think I can do that? : \- taking raw counts and dividing each gene by its length \- using the function rlog (DESeq2) on these counts divided by gene length (I would modify the rlog function to allow it to be used on decimal data).

I was wondering if it would be ok to do that.

Thank you very much for your help

# Answer
I am not sure what you mean by "decimal data", however, if you really want to use "gene length" information for normalization, take a look at `EDASeq`. Unlike `DESeq` or `EdgeR`, its normalization step takes gene length into account to produce "normalized counts", which you can feed into `DESeq` or `EdgeR` according to their manual.

Below is from the said manual:

> The normalized counts (or the original counts and the offset) obtained using the EDASeq package can be supplied to packages such as edgeR (Robinson, McCarthy, and Smyth 2010) or DESeq (Anders and Huber 2010) to find differentially expressed genes.
> 1 votes

# Answer
I've found that normalisation of variance-stabilised counts from short-read sequencing based on the length of the longest transcript isoform seems to work well. Here's a bit more detail on that from the methods of our paper:

> A variance-stabilizing transformation (VST) was applied to count data using the DESeq2 function varianceStabilizingTransformation, producing gene expression values that are adjusted for variance across sample conditions. To compare different genes to each other, this matrix was further transformed by dividing by transcript length in kbp to generate VSTPk values.

\[FWIW, single-molecule long reads don't need any length-normalisation carried out on them\]
> 1 votes

# Answer
tximport is recommended by the DESeq author, that can account for gene lengths.

Of course, you don't really want to account for gene length, so much as transcript length. And since you might have different transcripts in there, you will have to use something like RSEM or kallisto which will intelligently figure out the right transcript length.
> 0 votes

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


# Question
Title: how to repeat records in fastq n times efficiently?
How to iterate/repeat a record n times in a fastq file using bioawk? I wrote a python code using biopython, but it is very very slow. So, I am wondering if I can get some help by using bioawk. Thank you so much!

for example, in the file.fq file, I'd like to repeat read\_id1 for 100times, read\_id2 for 100times, ..., and write them to a new file.

# Answer
> 0 votes
Not sure why you would want to do this, but if the aim is to repeat **all** records in the fastq/fasta files `n` times, the most efficient way I can think of is just to concatenate the orginal file `n` times, similar to this thread.

i.e adapting the accepted answer would look something like this (I have not tested this):

```
for i in {1..n};do cat file.fq >> file_ntimes.fq ; done

```

**Update based on comments:**

To replicate each fastq entry back to back, maybe this hacky unix command can be adjusted (say for n=4):

```
n=10
paste - - - - <file.fq | awk '{for(i=0;i<$n;i++)print}' | tr '\t' '\n' > file_10times.fq

```

The idea is to `paste` blocks of 4 lines (i.e. 1 fastq entry) into one line and replicate this line n times. `tr` is reversing the paste action (since the separation by paste is a tab character). This won't work if fastq file for some reason has tabs in it which I don't expect but didn't test either. Note, for fasta file this would be `paste - -` rather than `paste - - - -`

I haven't tested this, but if it works for you I could polish it up if needed.

**Update2**

Based on Xioafeis comment this is solution seems to work as well:

```
bioawk -cfastx '{for(i=1; i<=10; i++) {print "@"$name"\n"$seq"\n+\n"$qual}}' file.fq > file_10times.fq

```

---
Tags: fastq, shell, perl, bioawk, repeat
---


# Question
Title: How to get matching pattern along with ID in a single command in grep?
I have a file containing multiple fast sequneces. For a specific consensus pattern as input, I extracted all the matching patterns from target fasta sequences with

```
grep -o -E "CC[GT]AAA[GC][AC]TT[GC]" input.fasta

```

However, the above command will retrieve just the matching sequences from fasta sequences and **I also wanted to get the corresponding fasta header of each match**.

For example, if **input.fasta** file is something like this,

```
>Gene 1
TGATGAAAAATGATAGAT
ATTGGGGGAAAAAAAAAT

>Gene 2
TTTCCTAAAGATTGT
AAATTTAAAAATGTTTTT

```

(Gene 2 has matching subsequence CCTAAAGATTG)

Output:

```
CCTAAAGATTG   Gene2

```

I prefer a solution with grep. But other possible solutions also helpful.

# Answer
> 1 votes
Use this Perl one-liner:

```
perl -lne '$id = $1 if /^>(.+)/; ($m) = /(CC[GT]AAA[GC][AC]TT[GC])/; print join "\t", $id, $m if $m;' input.fasta

```

The Perl one-liner uses these command line flags:
`-e` : Tells Perl to look for code in-line, instead of in a file.
`-n` : Loop over the input one line at a time, assigning it to `$_` by default.
`-l` : Strip the input line separator (`"\n"` on \*NIX by default) before executing the code in-line, and append it when printing.

**SEE ALSO:**
`perldoc perlrun`: how to execute the Perl interpreter: command line switches

---
Tags: fasta, motifs, grep
---


title: freebayes segmentation fault

text:
# Question
Title: freebayes segmentation fault
I have a draft, Illumina-based, genome and bunch of reseq samples. I generated bam files using `bowtie2` and marked duplicates using `picard`.

Now, when I tried to call SNPs using freebayes, I got cryptic `Segmentation fault (core dumped)` even when I tried to run it in verbose debugging mode (`-dd`). Do you have an idea where might be a problem?

freebayes -f data/Afus1/genome.fa -b reseq1.sorted.rmdup.bam reseq2.sorted.rmdup.bam -dd
  loading fasta reference data/Afus1/genome.fa
                                                                         Opening 2 BAM format alignment input files

     done
Number of ref seqs: 821466
Number of target regions: 0
no sample list file given, reading sample names from bam file
found 2 samples in BAM file
Segmentation fault (core dumped)
```
</code></pre>

# Answer
> 2 votes
Oh you silly sausage!

This obviously means that you are running out of stack memory. Apparently...

> freebayes is creating a single string containing the entire vcf header, which gets quite big when the reference has a lot of contigs with long names, causing stack space depletion.

Try increasing the maximum stack size (in this example to 128M) and running your command again.

```
ulimit -S -s 131072

```

---
Tags: snp
---


title: Can SBOL encode DNA Assembly linkers?

text:
# Question
Title: Can SBOL encode DNA Assembly linkers?
Where in an SBOL file could you specify the linkers that would be placed as prefixes / suffixes between parts for assembly? Lets say you have a final construct of your design and want to assemble the component parts - could you encode what linkers you used between each component, ie as an annotation? Would this be different between SBOL versions?

# Answer
> 2 votes
Every version of SBOL allows you to specify prefixes/suffixes as annotations on a sequence.

* In SBOL 3, the construct would be a `Component` of type `DNA` (https://identifiers.org/SBO:0000251), and the linker a `SequenceFeature` on that `Component`. The `role` of the sequence feature might depend on the particular assembly you are describing, and whether you're talking before or after assembly (in those cases where the process changes the sequence), but a good one in many cases would be Nucleotide Binding Site (https://identifiers.org/SO:0001655)
* In SBOL 2, the representation is the same, except the construct would be a `ComponentDefinition` and the linker a `SequenceAnnotation`.
* In SBOL 1, the representation is still the same, except the construct would be a `DnaComponent` and the linker an SBOL 1 `SequenceAnnotation` (which has different fields than in SBOL 2).

If the assembly method does not change the sequence, then all three can represent this pretty much equally well. If you want to represent a sequence that is modified by assembly or the plasmid carriers from which parts are extracted, then SBOL 3 is the best because it cleanly represents deletions during assembly via the `sourceLocation` field.

---
Tags: assembly, sbol
---


# Question
Title: Regarding canonical analysis of principal coordinates (CAP) for linking microbial data with elements
I used the below code to do CAP analysis to link microbial data with elements, but not sure how to do statistical analysis to show which element variation is responsible most for the difference in microbial composition. I can only able to see two axis in figure so how I can get complete information about the variation based on these mineral element and individually like which mineral responsible most for microbial diversity?

```
    # CAP ordinate
    cap_ord <- ordinate(
        physeq = phyloseq_prop,
        method = "CAP",
        distance = bray_not_na,
        formula = ~ ~ P + Mn + Cu +Zn + Rb + Mg )

        # CAP plot
    cap_plot <- plot_ordination(
      physeq =phyloseq_prop ,
      ordination = cap_ord,
        color = "SampleType",
        axes = c(1,2)
    ) +
        aes(shape = Region) +
        geom_point(aes(colour = SampleType), alpha = 0.4, size = 4) +
        geom_point(colour = "grey90", size = 1.5) +
        scale_color_manual(values = c("#CBD588", "#5F7FC7", "orange","#DA5724", "#508578", "#CD9BCD",
       "#AD6F3B", "#673770","#D14285","#652926", "#C84248",
      "#8569D5")
        )

        # Now add the environmental variables as arrows
    arrowmat <- vegan::scores(cap_ord, display = "bp")

    # Add labels, make a data.frame
    arrowdf <- data.frame(labels = rownames(arrowmat), arrowmat)

    # Define the arrow aesthetic mapping
    arrow_map <- aes(xend = CAP1,
        yend = CAP2,
        x = 0,
        y = 0,
        shape = NULL,
        color = NULL,
        label = labels)

    label_map <- aes(x = 1.3 * CAP1,
        y = 1.3 * CAP2,
        shape = NULL,
        color = NULL,
        label = labels)

    arrowhead = arrow(length = unit(0.02, "npc"))

    pdf("cap2.pdf", height=8, width=8)
    # Make a new graphic
    cap_plot +
      geom_segment(
        mapping = arrow_map,
        size = .5,
        data = arrowdf,
        color = "gray",
        arrow = arrowhead
      ) +
      geom_text(
        mapping = label_map,
        size = 4,
        data = arrowdf,
        show.legend = FALSE
      )

      dev.off();

     cap_ord

     Call: capscale(formula = distance ~  P + Mn + Cu +Zn + Rb + Mg, data = data)

                                  Inertia Proportion Rank
Total          2.069e+01  1.000e+00
Constrained    5.240e+00  2.533e-01   10
Unconstrained  1.545e+01  7.468e-01   97
Imaginary     -1.090e-03 -5.268e-05    1
Inertia is squared Bray distance

Eigenvalues for constrained axes:
  CAP1   CAP2   CAP3   CAP4   CAP5   CAP6   CAP7   CAP8   CAP9  CAP10
2.5733 0.8789 0.5647 0.3797 0.3037 0.1873 0.1243 0.1030 0.0800 0.0453

Eigenvalues for unconstrained axes:
 MDS1  MDS2  MDS3  MDS4  MDS5  MDS6  MDS7  MDS8
4.790 1.343 1.138 0.916 0.784 0.649 0.518 0.501
(Showing 8 of 97 unconstrained eigenvalues)

```

Many thanks

# Answer
I don't think that you can say that the 25% can be said to relate to the chemical variable specifically; rather, that 25% is the variance that is explained by all constrained axes. While these axes are drawn in such a way that they maximize variance associated with your independent variables, I do not think that they are directly interpretable as being 100% due to the elements.

The number of constrained axes is a parameter chosen as an input to the CAPscale method:

> An extremely important point is to determine how many PCO axes should be retained; i.e., what should be the choice for $m$? One might be concerned that the method could suffer from a poor choice for $m$. If $m$ is too small, then there may be some ecologically important information in the data that will not be included in the canonical analysis. On the other hand, if $m$ is too large relative to $N$, then a misleading canonical plot could result. For example, we cannot set $m = N - 1$, for then the squared canonical correlations would all be equal to 1.0, because we have $N$ points!

So the value of the variance there will vary depending on how many constrained axes you use. If I understand right, you could obtain 100% distance variance explained by choosing as many constrained axes as you have observations.

I think that it is worth looking into how others have solved this problem. For example, one has to first decide whether or not this 25% is an unexpectedly large number, i.e. compute p-values (I am not crazy about p-values in most contexts but here they could serve as a helpful heuristic), as suggested by the CAPscale paper:

> As a further option, one may test the hypothesis of either (1) no significant differences in multivariate location among groups (CDA), or (2) no significant relationship with quantitative environmental or other variables (CCorA). This is done by using the trace statistic (sum of canonical eigenvalues (sum of squared canonical correlations; see the Appendix) and obtaining a P value by permutation (e.g., Anderson 2001b).

This seems to be directly implemented as the `anova.cca()` function in the vegan package (I believe that they have decided that CAP is a flavor of CCA for that package). I understand the motivation to get % variance as it's intuitive, but in complex statistical analysis of distances it is a little bit less clear what it really means. Some kind of comparison to a null hypothesis is probably necessary to say whether it is meaningful.

For how others have approached similar questions you can see this discussion. For other implementations that may be more user-friendly, you might consider the vegan package, which is a richer toolbox for ordination in ecology (even if it is less focused on metagenomics). Here is the CAPscale function in vegan. The phyloseq authors have correctly decided to build large pieces of phyloseq out of vegan, so I think that you can pass phyloseq back into vegan functions and e.g. use those analyses directly on your `cap_ord` object. But see that link for more info.

For a little more easy-to-understand background on statistics of distance methods (particularly in microbiomes), you might see here.
> 1 votes

---
Tags: r, statistics, microbial-genomics, phyloseq
---


# Question
Title: How to check if a given gene is expressed in a group of microarray samples if I do not have control group to compare with?
I have a microarray gene expression dataset consisting of placenta samples. I want to check whether these placenta samples are mixed with maternal decidua tissue. I have marker genes for decidua samples SCARA5, PIGF, SLPI. How to check whether these genes are expressed in placenta samples? In other words, how to check if a gene is up or down regulated without control group?

# Answer
> 3 votes
You ask about which "genes are expressed" and then you mention "if a gene is up or down regulated". These are different, and given your application I think what you actually want to know if these marker genes are expressed. This is not a question of up/down-regulation relative to a control, so you do not need control samples for this.

It seems most of your datasets are Affymetrix, and there are R/Bioconductor packages that will give you the probability a gene is expressed. One such I found is frma (https://bioconductor.org/packages/release/bioc/html/frma.html), which returns an object "weight" that is "a vector of weights which roughly correspond to the probability of expression for each gene." I haven't worked with microarrays in a while, but I recall many other R/Bioconductor packages offering this functionality.

A detection p-value is calculated for Illumina array probes in Illumina's native BeadStudio software.

---
Tags: gene-expression, transcriptome, microarray
---


# Question
Title: How to tell which channel is the reference channel in two channel (red green) array data?
I'm sorry if this is a stupid question - I have not analysed two channel array data before.

I have a GenePix data file, which does not seem to indicate which channel (red or green) is used as reference channel.

If it possible to tell this from the data itself?

Where raw is my data object, and G is a dataframe of the green channel data and R is the data frame of the red channel data, the ranges are as follows:

```
#range(raw$G)
[1]    45 21939
>
> range(raw$R)
[1]  39 332

```

# Answer
> 2 votes
No, you can't get this from the data itself. If there is a control/reference sample on one channel, you need someone to annotate which channel it's on.

---
Tags: microarray
---


# Question
Title: Retrieve overlap widths with Bioconductor
I want to retrieve all overlaps between two sets of genomic intervals (GRanges) using Bioconductor. I am interested in the widths of the resulting overlapping ranges.

In cases with discontinuous ranges, `intersect()` works, but I am interested in the case where both sets contain contiguous ranges (no gaps).

I am using the findOverlaps function as in the example below:

```
gr1 <- GRanges(seqnames="chr2", ranges=IRanges(c(1, 11, 21, 31), c(10, 20, 30, 40)))
gr2 <- GRanges(seqnames="chr2", ranges=IRanges(c(1, 5, 15), c(4, 14, 30)))
hits <- findOverlaps(gr1, gr2)

```

This code yields the following intersections in a "Hits" object:

```
Hits object with 5 hits and 0 metadata columns:
      queryHits subjectHits
      <integer>   <integer>
  [1]         1           1
  [2]         1           2
  [3]         2           2
  [4]         2           3
  [5]         3           3
  -------
  queryLength: 4 / subjectLength: 3

```

But I could not find a way to extract the intersections widths from the resulting object.

It this the wrong approach ?

> Note: Using `intersect(gr1, gr2)` in this case yields a single contiguous overlap range spanning the whole chromosome, but I want each overlap to be returned independently, just like findOverlaps does.

# Answer
Whilst `intersect` is sort of doing what you want, it merges adjacent overlaps together which is not what you want I suppose. There doesn't seem to be dedicated `Grange` function to do this out of the box, but you could go through all overlapping pairs resulting from `findOverlap` and calculate it's widths using pintersect:

> punion, pintersect, psetdiff and pgap are generic functions that compute the element-wise (aka "parallel") union, intersection, (asymmetric!) difference and gap between each element in x and its corresponding element in y.

> "pintersect"(x, y, resolve.empty=c("none", "max.start", "start.x"), ...) # x, y IRanges objects.

```
hits <- findOverlaps(gr1, gr2)

width(pintersect(gr1[queryHits(hits)],
                 gr2[subjectHits(hits)]))
[1]  2  6  4  6 10

```
> 1 votes

---
Tags: r, bioconductor, genomicranges
---


# Question
Title: How do I infer someone's haplotype of a gene from data on several of the gene's SNPs?
For a a few SNPs from given gene, I know there happen to be several haplotypes that most of the population has. But from a person's data on those SNPs, we have two alleles at each SNP, so it's not trivial to match the SNP data to infer the person's haplotypes. When I search this I find tools that seem to be for haplotype reconstruction and SNP imputation, such as SHAPEIT, but this is not what I'm looking for.

Might someone know of a tool for this inference of known haplotypes from SNP data? It would be helpful even if someone could give me a name for this problem, as I'm not sure of the right term to search.

# Answer
It depends on the scale of the data (i.e. number of individuals, number of SNPs). Generally this can be done by matching up individual SNP genotypes to haplotype tag SNPs. The specific combination of genotyped SNPs is used to determine the underlying haplotype, based on tag SNP frequencies for each haplotype.

More details of the process of identifying tag SNPs can be found here:

https://www.broadinstitute.org/haploview/tagger
> 1 votes

---
Tags: snp, haplotypes
---


# Question
Title: How does one convert a log-likelihood substitution matrix into a probability matrix?
I am trying to calculate the substitution probabilities for amino acids (i.e. the probability of one amino acid given another). In theory, you should be able to derive them from a log-likelihood matrix (also called log-odds), such as Blosum62.

Does anyone know how to calculate probabilities from log-likelihoods or have a good modern source for a substitution probability matrix?

# Answer
**UPDATE2**

@chubs wants to compare different clades. This is possible within recent releases of ConSurf:

It is additionally quite possible to use ConSurf's default outputs to get per-position conservation scores that can be used for post-processing in any fashion, as well as prefab scripts that can be used with PyMol to make custom visualization:

> For all cases, ConSurf creates the following outputs:
>
> 1. The sequence and MSA colored by ConSurf conservation scores.
> 2. A text file that summarizes for each position the normalized score calculated, the assigned grade (1-through-9), the reliability estimation (for the Bayesian method) and the amino acids/nucleotides observed in the respective MSA column.
> 3. The sequences selected for the MSA and the MSA constructed (unless those files were uploaded by the user).
> 4. A file with the frequency of each amino acid/nucleotide observed in each column of the MSA.
> 5. The evolutionary tree, which was calculated by the server or uploaded by the user, together with the MSA are shown using the WASABI platform (Veidenberg et al., 2015). Moreover, using WASABI, a user can select a subtree containing a fraction of the homologous sequences and conduct a follow-up ConSurf analysis with these selected sequences. To refine a ConSurf analysis to a selected subtree, a user can choose any internal node on the phylogenetic tree and open a WASABI menu by using a right mouse click. Selecting the option 'run ConSurf on subtree' will issue a new window with the new ConSurf run for the selected subtree sequences (see example in Figure 3).

**UPDATE**

It seems that what @chubs really wants to do is map protein conservation onto a structure. There are various tools for doing this (here, here). ConSurf is the older and more established tool, using evolutionary conservation to rule out trivial apparent associations due to phylogenetic non-independence.

here is the somewhat related 3DPatch workflow, which uses profile HMM information as a measure of conservation:

Here is the output that they get:

Likewise ConSurf gets similar looking outputs:

**ORIGINAL ANSWER**

I think that more information is needed here to answer this question properly (what is your data? what is your ultimate goal? what are the things you've already tried?), but we can say a few things right away.

I would strongly suggest looking into substitution models, which are literally substitution probability matrices. Note that these are probabilities *per unit time*, which is to say that they are calibrated to e.g. evolutionary rates along a specific phylogenetic tree. I would therefore strongly suggest incorporating a phylogenetic model of some kind.

If you have e.g. a multiple sequence alignment, you can estimate a phylogenetic tree by one of the many tree estimation methods, which will give you a substitution probability matrix "for free" as this estimation is taken care of in the process of tree estimation. These will again be substitution probabilities per unit time; you will honestly not be able to get rid of the time-dependence here, as that's how evolution works.

I am a little behind the times, but some tools you might look into for tree estimation are PhyML and FastTree. If you only have sequences and haven't aligned them yet, you can use a tool like MUSCLE, ClustalO, or MAFFT.

BLOSUM and related models are derived from specific slices of protein sequence identity (evolutionary history) calibrated to certain scopes, it's unlikely that they are appropriate to your case, unless your data fits inside those rather narrow slices.
> 1 votes

---
Tags: sequence-alignment, phylogenetics
---


# Question
Title: Why is the UniProt REST API returning multiple results, when I am only providing one ID to be checked?
Can some explain why when I run a REST API query on UniProt for one ID ('Q9UJL9'), the results for three IDs are returned?

The query I'm using is:

https://www.uniprot.org/uniprot/?query=Q9UJL9&sort=score&columns=id,go(molecular%20function)&format=tab

(I didn't hyperlink it since you need to see the query to see where I'm going wrong).

When I run the query above, the output is:

```
Entry   Gene ontology (molecular function)
Q9UJL9  DNA-binding transcription factor activity, RNA polymerase II-specific [GO:0000981]; DNA-binding transcription repressor activity, RNA polymerase II-specific [GO:0001227]; metal ion binding [GO:0046872]; RNA polymerase II transcription regulatory region sequence-specific DNA binding [GO:0000977]
Q8WWY3  identical protein binding [GO:0042802]; ribonucleoprotein complex binding [GO:0043021]; RNA binding [GO:0003723]; snRNP binding [GO:0070990]; U4atac snRNA binding [GO:0030622]; U4 snRNA binding [GO:0030621]
Q6A162  structural molecule activity [GO:0005198]

```

I checked the uniprot page of the query I want (Q9UJL9), the other items are interactors (along with many others) but not synonyms/linked (it seems, e.g my query is a zinc finger, one of the linked records, Q6A162, is a keratin protein).

# Answer
it is because you do a text search and the "word" Q9UJL9 appears in multiple entries. If you want an id search, specifically say so and for that the query syntax is `id:Q9UJL9`. You will want to follow a redirect, and there are some rare cases where that might still lead to more than one entry if the id was made secondary and attached to multiple entries. e.g. P00001.
> 3 votes

---
Tags: uniprot, api
---


# Question
Title: installing multiple bioconductor packages at once
I was wondering if there is more elegant way of installing and loading multiple packages in Bioconductor similar to `pacman` with CRAN packages.

I tried:

```

# install and load  the package  manager
 if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

# install and load bioconductor packages
bio_pkgs = c("BiocManager", "ggbio", "biomaRt", "EnsDb.Hsapiens.v75", "cBioPortalData", "GenomicRanges", "GenomicFeatures", "TxDb.Hsapiens.UCSC.hg19.knownGene", "Homo.sapiens")

 if (!requireNamespace(bio_pkgs, quietly = TRUE))
    BiocManager::install(bio_pkgs)

library(easypackages)
libraries("ggbio", "biomaRt", "EnsDb.Hsapiens.v75", "renvcBioPortalData", "GenomicRanges", "GenomicFeatures", "TxDb.Hsapiens.UCSC.hg19.knownGene", "Homo.sapiens")

```

# Answer
> 7 votes
You can do something like this:

```
bio_pkgs <- c("ggbio", "biomaRt", "EnsDb.Hsapiens.v75", "renvcBioPortalData",
          "GenomicRanges", "GenomicFeatures", "TxDb.Hsapiens.UCSC.hg19.knownGene", "Homo.sapiens")

# install:
BiocManager::install(bio_pkgs)

# load all at once
invisible(lapply(pkgs, function(x) library(x, character.only=TRUE)))

```

Be aware though that if you version-control your packages with something like `renv` then it is encouraged to explicitely load every package as suggested here so `renv` can properly recognize the packages to snapshot.

---
Tags: bioconductor
---


# Question
Title: how to easily make complex Euler/Venn diagrams in R?
From my microarray data I have selected three subsets of probes. I would like to display these subsets as a Euler diagram to visualize the overlaps, but am struggling with inputting the combinations into Eulerr. All three subsets need to be nested within the full set of probes; one subset is fully contained within another, and the third overlaps with the other two. Is there a package into which I could feed the vectors of probe names and let it determine the overlaps for me? Many thanks in advance!

# Answer
Since you are in R, consider using the UpsetR library. It doesn't make Venn diagrams but it helps to visualize overlaps between any number of groups.

http://gehlenborglab.org/research/projects/upsetr/#:~:text=UpSetR%20is%20an%20R%20package,based%20on%20groupings%20and%20queries.
> 2 votes

# Answer
You might consider different ways of using the `venneuler` package. I'm not super thrilled with its usability (or with the accuracy of the overlaps), but you can get it to make custom Venn diagrams.

Here is a thread with some suggestions along these lines.

Below is some sample R code that I slightly modified from an old project that I used to make a 3-category Venn diagram.

I didn't run this to make sure that it works because of an irritating Java dependency issue on my OSX laptop that I don't feel like fixing right now, but as long as you have a JDK on your machine and can install the packages it *should* work.

```
# note irritating rJava dependency of venneuler
library(plotrix)
library(venneuler)

#A-B sharing: 200414
#B-C sharing: 17561
#A-C sharing: 16764
#All share 14201

# get coordinates for circles- replace as needed
# this function is supposed to plot stuff, but doesn't actually
ven = venneuler(c(C=48147-17561-16764+14201,
                  B=347181-17561-200414+14201,
                  A=325971-200414-16764+14201,
                  "C&B"=17561, "B&A"=200414, "C&A"=16764, "A&B&C"=14201))

cats = c('A','B','C')

# note- this following only shows up in the .pdf printer, not in the X11 utility plotter
pdf('venn.pdf')
plot(c(-.5,1.5), c(-.5,1.5), type='n', xaxt='n', yaxt='n', bty='n', xlab='', ylab='')
draw.circle(ven$center[1,1], ven$center[1,2], ven$diameters[1],
        col = hsv(1,1,1,.5))
draw.circle(ven$center[2,1], ven$center[2,2], ven$diameters[2],
            col = hsv(.5,1,1,.5))
draw.circle(ven$center[3,1], ven$center[3,2], ven$diameters[3],
            col = hsv(.1,1,1,.5))

text(ven$center[1,1]+.065, ven$center[1,2], labels=cats[3], cex=1.3)
text(ven$center[2,1], ven$center[2,2]-.2, labels=cats[2], cex=1.3)
text(ven$center[3,1], ven$center[3,2]+.2, labels=cats[1], cex=1.3)

dev.off()

```
> 0 votes

---
Tags: r
---


# Question
Title: How can I extract information from .sam files?
I have 10 `.sam` files after my `bowtie2` alignment on ten single-pair sequences. I would like to build a graph based on that output data, however I don't see any other solution but just writing the alignment percentage manually after each iteration.

Is there a tool that can help me to build a `.csv` file so I can manipulate the output data more easily?

# Answer
Use `samtools flagstat` with option `-O tsv`:

> Using **-O tsv** selects a tab-separated values format that can easily be imported into spreadsheet software. In this format the first column contains the values for QC-passed reads, the second column has the values for QC-failed reads and the third contains the category names.

**SEE ALSO:**
What Does Samtools Flagstat Results Mean?

**INSTALLATION:**
Some old versions of `samtools` do not have option `-O tsv`. To install a recent versions of `samtools`, use either:

`conda` (macOS, Linux) (but see also: Tutorial: installing samtools with conda - overcoming two common errors):

```
conda create --name samtools samtools

```

or `brew` (macOS):

```
brew install samtools

```
> 2 votes

---
Tags: samtools, sam, bowtie2
---


# Question
Title: Find pattern that is present twice and allow <=2 mismatches on each
I have a fastq file of 400,000 reads (so speed is important). In the sequences there are barcodes integrated that should be present twice. Given a barcode, I want to find the sequences that have the barcode present twice with \<= 2 mismatches. So, with a barcode 'ATTCGACCGATAGG', I would like to retrieve all of the following sequences-

> TATCTTGTGGAAAGGACGAAACACCGAACACAAAGCATAGATGCGTTTAAGAGCTATGCTGGAAAACAGCATAGCAAGTTTAAATAAGGCTAGTCCGTTATCAACTTGAAAAAGTGGCACCGAGTCGGTGCTTTTTTT**ATTCGACCGATAGG**GGTGGCAGGGGAGGCCGAGGAGGAAGAAGGGGAGGTGGCAG**ATTCGACCGATAGG**TGGCGTAACTAGATCTTGAGACAAA TATCTTGTGGAAAGGACGAAACACCGGTCCGAGCAGAAGAAGAAGTTTAAGAGCTATGCTGGAAACAGCATAGCAAGTTTAAATAAGGCTAGTCCGTTATCAACTTGAAAAAGTGGCACCGAGTCGGTGCTTTTTTT**ATTCGACCGATAGG**GGTGGCAGGGGAGGCCGAGGAGGAAGAAGGGGAGGTGGCAG**ATTCGACCGATAGG**TGGCGTAACTAGATCTTGAGACAAA TATCTTGTGGAAAGGACGAAACACCGAGTCCGAGCAGAAGAAGAAGTTTAAGAGCTATGCTGGAAACAGCATAGCAAGTTTAAATAAGGCTAGTCCGTTATCAACTTGAAAAAGTGGCACCGAGTCGGTGCTTTTTTT**ATTCGACCGATAGG**GGTGGCAGGGGAGGCCGAGGAGGAAGAAGGGGAGGTGGCAG**ATTCGACCGATAGG**TGGCGTAACTAGATCTTGAGACAAA TATCTTGTGGAAAGGACGAAACACCGAGTCCGAGCAGAAGAAGAAGTTTAAGAGCTATGCTGGAAACAGCATAGCAAGTTTAAATAAGGCTAGTCCGTTATCAACTTGAAAAAGTGGCACCGAGTCGGTGCTTTTTTT**ATTCGACGATAGG**GGTGGCAGGGGAGGCCGAGGAGGAAGAAGGGGAGGTGGCAG**ATTCGACCGATAGG**TGGCGTAACTAGATCTTGAGACAAA

Note that the first barcode in the fourth sequence is short of one character. I have tried with biopython and regex but it's just too slow given I have 5k barcodes. I am wondering if there is a fast solution available in python or in something like grep, awk or anything else. Thanks.

# Answer
> 1 votes
This isn't perfect but you can build on it to deal with any issues (ie overalpping hits). Also if you need it faster just chunk it and run in parallel threads.

I just copy pasted your examples seqs into files to read (called them seqs and code).

```
#need blast installed
import subprocess
import pandas as pd

path = '/Users/ur/dir/'
#make a blast db
with open (path + 'seqs.txt', 'r') as fin:
    with open (path + 'seqs.fa', 'w') as fout:
        for i, line in enumerate(fin):
            fout.write('>'+str(i)+'\n')
            fout.write(line)

subprocess.run([
    'makeblastdb',
    '-in', path + 'seqs.fa',
    '-dbtype', 'nucl'
])

#use blast
subprocess.run([
    'blastn',
    '-task', 'blastn-short',
    '-query', path + 'code.txt',
    '-db', path + 'seqs.fa',
    '-out', path + 'out.txt',
    '-outfmt', '6'
])

#grab seqs with 2 hits (ie ~2x len(query)
df = pd.read_csv(
    path + 'out.txt',
    sep='\t',
    header=None,
    names=['query', 'db_name', 'percent','len'],
    usecols = [0,1,2,3]
)
df = df.groupby('db_name').agg('sum')
print(*list(df.query('len > 27').index))

```

I this case they are all hits. If you look at the blast output you'll see one hit was fragmented (blast isn't perfect with small seqs). But total len works to resolve this.

---
Tags: sequence-alignment, fastq, long-reads, bioawk, barcode
---


# Question
Title: Quast duplication ratio and mismatches percent
When analyzing Quast results it seems that it doesn't calculate mismatches and indels in a useful way if the "Duplication ratio" is over 1. For example, that's what I get for an assembly with two contigs covering almost the full genome

And that's what I get if I remove one of the contigs:

See how the mismatches and indels per 100kbp are almost double in the first case but that does not correspond to the reality.

Any better way to get these metrics? Would it be safe to divide error rates by the duplication ratio?

# Answer
I took a quick look at the quast code that I believe is responsible here:

```
def analyze_coverage(ref_aligns, reference_chromosomes, ns_by_chromosomes, used_snps_fpath):
    indels_info = IndelsInfo()
    genome_mapping = {}
    for chr_name, chr_len in reference_chromosomes.items():
        genome_mapping[chr_name] = [0] * (chr_len + 1)
    with open(used_snps_fpath, 'w') as used_snps_f:
        for chr_name, aligns in ref_aligns.items():
            for align in aligns:
                ref_pos, ctg_pos = align.s1, align.s2
                strand_direction = 1 if align.s2 < align.e2 else -1
                for op in parse_cs_tag(align.cigar):
                    if op.startswith(':'):
                        n_bases = int(op[1:])
                    else:
                        n_bases = len(op) - 1
                    if op.startswith('*'):
                        ref_nucl, ctg_nucl = op[1].upper(), op[2].upper()
                        if ctg_nucl != 'N' and ref_nucl != 'N':
                            # MISMATCHES INCREMENTED HERE
                            indels_info.mismatches += 1
                            if qconfig.show_snps:
                                used_snps_f.write('%s\t%s\t%d\t%s\t%s\t%d\n' % (chr_name, align.contig, ref_pos, ref_nucl, ctg_nucl, ctg_pos))
                        ref_pos += 1
                        ctg_pos += 1 * strand_direction
                    elif op.startswith('+'):
                        # INDELS INCREMENTED HERE
                        indels_info.indels_list.append(n_bases)
                        indels_info.insertions += n_bases
                        if qconfig.show_snps and n_bases < qconfig.MAX_INDEL_LENGTH:
                            ref_nucl, ctg_nucl = '.', op[1:].upper()
                            used_snps_f.write('%s\t%s\t%d\t%s\t%s\t%d\n' % (chr_name, align.contig, ref_pos, ref_nucl, ctg_nucl, ctg_pos))
                        ctg_pos += n_bases * strand_direction
                    elif op.startswith('-'):
                        indels_info.indels_list.append(n_bases)
                        indels_info.deletions += n_bases
                        if qconfig.show_snps and n_bases < qconfig.MAX_INDEL_LENGTH:
                            ref_nucl, ctg_nucl = op[1:].upper(), '.'
                            used_snps_f.write('%s\t%s\t%d\t%s\t%s\t%d\n' % (chr_name, align.contig, ref_pos, ref_nucl, ctg_nucl, ctg_pos))
                        ref_pos += n_bases
                    else:
                        ref_pos += n_bases
                        ctg_pos += n_bases * strand_direction
                if align.s1 < align.e1:
                    for pos in range(align.s1, align.e1 + 1):
                        genome_mapping[align.ref][pos] = 1
                else:
                    for pos in range(align.s1, len(genome_mapping[align.ref])):
                        genome_mapping[align.ref][pos] = 1
                    for pos in range(1, align.e1 + 1):
                        genome_mapping[align.ref][pos] = 1
            for i in ns_by_chromosomes[align.ref]:
                genome_mapping[align.ref][i] = 0

    covered_ref_bases = sum([sum(genome_mapping[chrom]) for chrom in genome_mapping])
    return covered_ref_bases, indels_info

```

This indicates to me that it simply iterates over contigs (or rather, called SNPs from contig alignments to a reference) and counts up the number of mismatches/indels that occur in total across all aligned contigs.

This suggests to me that if haplotigs are present, the areas that they cover will be double-counted (or n-counted where there are n haplotigs, assuming that all haplotigs covering a reference position have roughly equal divergence).

From this point of view, I think that it is heuristically reasonable to divide by genome fraction. However, if you want something really rigorous you might consider digging into the quast `used_snps` file that this code iterates over to decide which SNPs/indels/whatever you think are worth keeping.
> 1 votes

---
Tags: ngs, assembly
---


# Question
Title: Differences in webtool and standalone Clustal Omega Percent Identity Matrix
I'm trying to align sequences using MUSCLE in my own computer. To compute Percent Identity Matrix I'm trying to use Clustal2.1 like it's done in https://www.ebi.ac.uk/Tools/msa/muscle/. Even though I'm doing the same thing, I get rounded up values in my computer.

For example, this is the webtool output:

```
#
#
#  Percent Identity  Matrix - created by Clustal2.1
#
#

     1: JX310755.1  100.00   92.79   65.29   63.98
     2: M97264.1     92.79  100.00   65.66   63.92
     3: AJ292230.1   65.29   65.66  100.00   66.88
     4: KM379144.1   63.98   63.92   66.88  100.00

```

And this is the output in my computer.

```
#
#
#  Percent Identity  Matrix - created by Clustal2.1
#
#

     1: JX310755.1     100      93      65      64
     2: M97264.1        93     100      66      64
     3: AJ292230.1      65      66     100      67
     4: KM379144.1      64      64      67     100

```

I have both installed clustalw using `sudo apt install clustalw` and by downloading the tarballs from http://www.clustal.org/download/current/clustalw-2.1.tar.gz, with exact results.

The commmand I'm using is `clustalw2 alignment -tree -pim`.

Any thoughts on why this could be? Thanks in advance.

# Answer
So, after a year I came out with the solution - as M\__ kindly suggested, using Clustal Omega.

First, align usign muscle:

`muscle -in fasta -quiet -out alignment`

Then, generate percent distance matrix using clustal omega:

`clustalo --infile alignment --distmat-out output --percent-id --full`
> 1 votes

# Answer
Try clustal omega, instead of clustalw because it is really old. It should report to 2 dp. Variations do occur depending on how indels are treated, but this is just rounding.

You could try

```
alignment -tree -pim -score

```

The idea is to convert to absolute rather than percentage and might overcome the rounding error
> 1 votes

---
Tags: sequence-alignment
---