# Question
Title: Filtering reads greater than 5 from HT-seq count files
I have some raw counts from HTSeq after aligning with hg38 human reference genome. I want to do filtering in a way that the filtered count files should have the same number of lines. The reason behind this is that I want to use DESeq2 for differential gene expression which requires all the count files to be of the same length to carry out the analysis.

I have written a python program to select reads greater than 5 from the second column.

My count files look like this

```
ENSG00000000003 20
ENSG00000000005 0
ENSG00000000419 123
ENSG00000000457 35
ENSG00000000460 56
ENSG00000000938 0
ENSG00000000971 0
ENSG00000001036 164
ENSG00000001084 83
ENSG00000001167 68
ENSG00000001460 8
ENSG00000001461 13
ENSG00000001497 80
ENSG00000001561 0
ENSG00000001617 39
ENSG00000001626 0
ENSG00000001629 146
ENSG00000001630 4

```

Python program I have written is as follows

```
import sys

inputlist = [sys.argv[1], sys.argv[2], sys.argv[3], sys.argv[4], sys.argv[5], sys.argv[6], sys.argv[7], sys.argv[8]]

for i in range(len(inputlist)):
    inputfile = open(inputlist[i], 'r')
    for line in inputfile:
        splitline = line.strip().split('\t')
        if int(splitline[1]) >= 5:
            if not splitline[0] in genelist:
                genelist.append(splitline[0])
    inputfile.close()

outputlist = open(sys.argv[9], 'w')
for gene in genelist:
    outputlist.write(gene + '\n')
outputlist.close()

for i in range(len(inputlist)):
    inputfile = open(inputlist[i], 'r')
    outputfile = open(sys.argv[i+10], 'w')
        for line in inputfile:
        splitline = line.strip().split('\t')
        if splitline[0] in genelist:
            outputfile.write(line)
    inputfile.close()
    outputfile.close()

```

`python filtercounts.py 'mir99-1-L1count.txt' 'mir99-1-L2count.txt' 'mir99-2-L1count.txt' 'mir99-2-L2count.txt' 'mir99-3-L1count.txt' 'mir99-3-L2count.txt' 'genelist.txt' 'mir99-1-L1filtered.txt' 'mir99-1-L2filtered.txt' 'mir99-2-L1filtered.txt' 'mir99-2-L2filtered.txt' 'mir99-3-L1filtered.txt' 'mir99-3-L2filtered.txt'`

This program has worked well for me in past , but now I see that although we see many such reads removed but there are many zero reads or reads less than 5 in the dataset.

Some of my output looks like this

```
ENSG00000000003 0
ENSG00000000419 0
ENSG00000000457 15
ENSG00000000460 12
ENSG00000001036 0
ENSG00000001084 0
ENSG00000001167 0
ENSG00000001460 1
ENSG00000001461 0
ENSG00000001497 1
ENSG00000001617 0
ENSG00000001629 0
ENSG00000001631 0
ENSG00000002016 2
ENSG00000002330 2
ENSG00000002549 0
ENSG00000002822 0
ENSG00000002834 0
ENSG00000002919 0
ENSG00000003056 4
ENSG00000003147 0
ENSG00000003249 1

```

Any help regarding the problematic behavior of the problem will be appreciated.

# Answer
After seeing that you want to use DESeq2, I'd recommend you do the count filtering within R, after generating the count matrix. Note that DESeq2 includes functions for loading from HTSeq-count output files:

https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#htseq-count-input

It is not necessary to filter out low count data when processing with DESeq2, but if desired this can be done after converting to a DESeq2 object:

> While it is not necessary to pre-filter low count genes before running the DESeq2 functions, there are two reasons which make pre-filtering useful: by removing rows in which there are very few reads, we reduce the memory size of the dds data object, and we increase the speed of the transformation and testing functions within DESeq2. It can also improve visualizations, as features with no information for differential expression are not plotted.

```
keep <- apply(counts(dds), 1, function(x){all(x >= 5)} # for all gene counts >= 5
keep <- apply(counts(dds), 1, function(x){any(x >= 5)} # for any gene counts >= 5
dds <- dds[keep,]

```

https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#pre-filtering

Note that DESeq2 cares about replicates, and will almost always pin an adjusted p-value near 1 if any conditions have 1 or 2 replicates (three is the minimum number required to identify outlier values). There's not really any way to work around that other than including more replicates for the low-replicate conditions.

Another program might give different results, but that won't get around the fundamental issue: two replicates is not sufficient to determine sample variation. Either accept that p-values cannot be used as a reliable estimate of model fit this case, or add more replicates.
> 1 votes

---
Tags: python, deseq2, genomics
---


# Question
Title: How can I assemble my genome from raw files?
I've had my whole genome sequenced (at 30x average coverage) by a lab, and they have provided the raw files to me (BAM, FASTQ, and VCF).

How can I assemble it?

And does assembly provide any further insights that I can't get from the separated fragments?

I'd like to assemble and annotate it, and find insights, but I don't know where to start.

# Answer
Given you've tagged this as 'human-genome' there is likely little value in performing genome assembly. De-novo assembly of WGS is generally for non-model organisms without a reference genome available, and would require much higher than 30X coverage (and ideally long reads). By providing the BAM and VCF, they'e already saved you some upstream work unless there are needs to re-process the data from FASTQ
> 2 votes

# Answer
No, assembling won't tell you anything new. When working with organisms whose reference genome is known, all you care about is where *your* sequences differ from the reference. If you know where your sequences differ, then you can infer the rest, they will be whatever is in the reference.

All you want to do is identify *variants*, those sites where your genome differs from the reference. That's what the VCF file contains. Bear in mind that the reference isn't some sort of platonic ideal of the perfect human, it's just the collected, pooled DNA of 13 anonymous, apparently healthy individuals. You can expect anywhere from about 3 to 5 million differences in your own genome, the overwhelming majority of which will be completely innocuous.

Since you already have your VCF file, your list of variants, that's all you need to investigate. If you really want to, you can make your assembly by taking the reference genome and then adding the variants in your VCF file, but it won't be of much practical use.
> 2 votes

# Answer
To add the "do not assemble your genome" answers. There might be specific circumstances when you would actually like to do that. For example, if the 30x you got are HiFi long reads and have a reason to suspect there are some interesting large rearrangements in your genome (rare disease?), it might be informative to actually assemble the genome and compare it to the reference to find these rearrangements.

Of course there is also plenty of possibilities in between (still assuming long reads), for example to map your sequencing data to a pan-genome graph or simply use a mapper that is aware of genomic rearrangements (like ngmlr) and call the structural variants.

But in the end it boils down to two aspects you did not mention in your question: What would you like to learn about your genome?
> 2 votes

---
Tags: assembly, human-genome, sequence-annotation
---


# Question
Title: include a glimma interface in a shiny app
I am trying to code a shiny app for RNA-Seq data analysis.

I would like to include glimma interactive plots in it. However, in my current interface, clicking the action button `Glimma plot` opens the glimma interface in a new page of my web browser and not in the shiny UI.

I use `shiny`, `glimma` and `EdgeR`. Here is my code:

```
ui <- fluidPage(
  actionButton(inputId="run.glimma", label="Glimma plot"),
  htmlOutput("glimma")

server <- function(input, output, session){
  observeEvent(input$run.glimma, {
    glMDSPlot(dgeObj.norm(), labels=rownames(dgeObj.norm()$samples), groups=dgeObj.norm()$samples[,as.character(input$mds.grouping.feature)])
  }),

  output$glimma <- renderUI({
      includeHTML("glimma-plots/MDS-Plot.html")
    })
}

```

I only indicated the part of the code that I use to generate the `glimma` MDS plot.

Is there a way to display the `glimma` plots in the shiny app UI and not in a new webpage ?

Thank you

# Answer
I have come across this page today with same query. After visiting here and there, finally following approach helped to resolve this issue. I must confess that I am just a beginner to R Shiny.

Lets say we have ui.R, server.R and global.R.

**Inside global.R, add following code:**

```
glimmaMDSOutput <- function(outputId, width = "100%", height = "400px") {
        htmlwidgets::shinyWidgetOutput(outputId, 'glimmaMDS', width, height, package = 'Glimma')
}

renderglimmaMDS <- function(expr, env = parent.frame(), quoted = FALSE) {
        if (!quoted) { expr <- substitute(expr) }
        htmlwidgets::shinyRenderWidget(expr, glimmaMDSOutput, env, quoted = TRUE)
}

```

**Inside ui.R, add following:**

```
glimmaMDSOutput(outputId = 'MDSplot')

**Inside server.R, add following:**
output$MDSplot <- renderglimmaMDS({
                Glimma::glimmaMDS(x = .., labels = .., groups = .., top = .., gene.selection = "pairwise", launch = FALSE, html = NULL, width = 1000, height = 500);
})

```
> 2 votes

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


# Question
Title: Getting the licenses of all conda installed packages
I would like to get a list of the licenses of all packages I installed using Conda. pip has what I need in the form of pip-licenses, is there a similar thing for conda?

# Answer
To my knowledge there's no single listing of the licenses, but you can:

```
conda list | awk '{if(NR>3) printf("%s=%s", $1, $2)}' | xargs conda info | grep license

```

Allegedly `conda info --license` should work, but at least in conda 4.6.3 it doesn't show any license information. It unfortunate that there's no convenient way to get a tabular list of package:version:channel:license information from the command line.
> 7 votes

# Answer
This one uses JQ for parsing and querying conda generated json output (https://stedolan.github.io/jq/)

```
conda list | awk '{if(NR>3) printf("%s=%s\n", $1, $2)}' | xargs -I{} conda search --info --json "{}" | jq --raw-output '.[][0] | "\(.name)\t\(.version)\t\(.license)"'

```

builds the list, asks conda search for all info on package as json, feeds json to JQ, which selects the first item and outputs its name, version and license as tsv.
> 1 votes

# Answer
I can't comment, but here's an improvement for the top answer to list more information:

```
conda list | awk '{if(NR>3) printf("%s=%s\n", $1, $2)}' | xargs conda info | grep 'name\|version\|license'

```
> 1 votes

---
Tags: bioconda, conda
---


# Question
Title: Does cutadapt trim trailing N's first and then use max_n to filter reads?
**Background** I want to trim leading and (likely just) trailing N's from my WES (Illumina NextSeq500) reads with `cutadapt` (`--trim-n`).

I also want to filter for reads with \> 40%? N's (`--max-n`).

**Question** Does `cutadapt` first remove the trailing/leading N's before removing reads with \>X% N's? This is because if, hypothetically, 40% of the N's are on the end of the reads, I want to keep them.

Cheers and thanks in advance.

# Answer
> 3 votes
Yes, it behaves just like you expect, where a read with \>40% Ns but only at the ends will still be kept (after trimming).

---

Cutadapt runs in two stages, with read modification first and read filtering after. So first, `--trim-n` removes Ns at the start and end of reads are removed, and only then is a rule like `--max-n` applied.

Here's an example where the first 5 sequences have Ns at the ends and the next 5 have Ns interspersed:

```
>seq01  0%
CACATCAAGAGGTCCTTGGCGACTATTAGCGCCCGGCGGACCCAATCGTA
>seq02 20%
NNNNNCAAGAGGTCCTTGGCGACTATTAGCGCCCGGCGGACCCAANNNNN
>seq03 40%
NNNNNNNNNNGGTCCTTGGCGACTATTAGCGCCCGGCGGANNNNNNNNNN
>seq04 60%
NNNNNNNNNNNNNNNTTGGCGACTATTAGCGCCCGNNNNNNNNNNNNNNN
>seq05 80%
NNNNNNNNNNNNNNNNNNNNGACTATTAGCNNNNNNNNNNNNNNNNNNNN
>seq06 10%
CACATCANGAGGTCCTTNGCGACTATTNGCGCCCGGCNGACCCAATCNTA
>seq07 20%
CANATCANGAGGNCCTTNGCGANTATTNGCGCNCGGCNGACCNAATCNTA
>seq08 40%
CNNATCNNGAGNNCCTNNGCGNNTATNNGCGNNCGGNNGACNNAATNNTA
>seq09 60%
CNNANCNNGNGNNCNTNNGNGNNTNTNNGNGNNCNGNNGNCNNANTNNNA
>seq10 80%
CNNNNNNNGNNNNCNNNNGNNNNTNNNNGNNNNCNNNNGNNNNANNNNNA

```

With just read filtering:

```
$ cutadapt --quiet --max-n 0.4 seqs.fa -o -
>seq01  0%
CACATCAAGAGGTCCTTGGCGACTATTAGCGCCCGGCGGACCCAATCGTA
>seq02 20%
NNNNNCAAGAGGTCCTTGGCGACTATTAGCGCCCGGCGGACCCAANNNNN
>seq03 40%
NNNNNNNNNNGGTCCTTGGCGACTATTAGCGCCCGGCGGANNNNNNNNNN
>seq06 10%
CACATCANGAGGTCCTTNGCGACTATTNGCGCCCGGCNGACCCAATCNTA
>seq07 20%
CANATCANGAGGNCCTTNGCGANTATTNGCGCNCGGCNGACCNAATCNTA
>seq08 40%
CNNATCNNGAGNNCCTNNGCGNNTATNNGCGNNCGGNNGACNNAATNNTA

```

With trimming also:

```
$ cutadapt --quiet --max-n 0.4 --trim-n seqs.fa -o -
>seq01  0%
CACATCAAGAGGTCCTTGGCGACTATTAGCGCCCGGCGGACCCAATCGTA
>seq02 20%
CAAGAGGTCCTTGGCGACTATTAGCGCCCGGCGGACCCAA
>seq03 40%
GGTCCTTGGCGACTATTAGCGCCCGGCGGA
>seq04 60%
TTGGCGACTATTAGCGCCCG
>seq05 80%
GACTATTAGC
>seq06 10%
CACATCANGAGGTCCTTNGCGACTATTNGCGCCCGGCNGACCCAATCNTA
>seq07 20%
CANATCANGAGGNCCTTNGCGANTATTNGCGCNCGGCNGACCNAATCNTA
>seq08 40%
CNNATCNNGAGNNCCTNNGCGNNTATNNGCGNNCGGNNGACNNAATNNTA

```

---
Tags: quality-control, trimming, cutadapt
---


# Question
Title: pip installation within conda env: Error downgrading scikit-learn - WARNING: No metadata
I am trying to get the bioinformatics tool `deepbgc` up and running. Now I have created a `conda` environment and have currently installed `deepbgc` onto it. Now when I run the tool, the program works but it outputs the following warning message as a result:

```
 /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/sklearn/base.py:306:
 UserWarning: Trying to unpickle estimator RandomForestClassifier from version 0.18.2
 when using version 0.21.3. This might lead to breaking code or invalid results. Use at your own risk.

```

This led me to believe that the program is trying to utilize an older version of the `scikit-learn` package (`version 0.18.2` instead of the currently installed 0.21.3 version) for this `RandomForestClassifer` estimator. This pointed me to the next logical solution of downgrading my current version of `scikit-learn`. So, I ran a simple `pip` install of the older package version:

```
pip install scikit-learn==0.18.2

```

For some reason, the package downgrade failed and outputted the following error:

```
    error: Command "g++ -pthread -B /EFS/tools/miniconda/envs/deepbgc2/compiler_compat -
Wl,--sysroot=/ -Wsign-compare -DNDEBUG -g -fwrapv -O3 -Wall -fPIC -
I/EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/numpy/core/include -
I/EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/numpy/core/include -
I/EFS/tools/miniconda/envs/deepbgc2/include/python3.7m -c
sklearn/cluster/_dbscan_inner.cpp -o build/temp.linux-x86_64-cpython-
37/sklearn/cluster/_dbscan_inner.o -MMD -MF build/temp.linux-x86_64-cpython-
37/sklearn/cluster/_dbscan_inner.o.d" failed with exit status 1

      [end of output]

  note: This error originates from a subprocess, and is likely not a problem with pip.
  WARNING: No metadata found in /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages
  Rolling back uninstall of scikit-learn
  Moving to /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/scikit_learn-0.21.3.dist-info/
   from /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/~cikit_learn-0.21.3.dist-info
  Moving to /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/sklearn/
   from /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/~klearn
error: legacy-install-failure

× Encountered error while trying to install package.
╰─> scikit-learn

note: This is an issue with the package mentioned above, not `pip`.

```

Currently, I have `Python` version `3.7` installed on the environment and version 22.1.2 of `pip`. Any guidance appreciated, thanks!

# Answer
`pip` and `conda` don't mix and good advice is don't ever mix them. Thus to avoid this delete your old environment, create a new `conda` environment, or use a `conda` environment where `pip` has not been used ...

```
conda install -c anaconda scikit-learn=0.18.2

```

If you are using `pip` it is recommended to use it only with `virtualenv` `venv`. If there is a dependency issue delete the environment and start over is much easier.

The rationale is for not mixing them is:

* `conda` installs lots of stuff not just Python and it gets the `env` (environment) right for this purpose
* `pip` only installs Python within a Python environment

---

The simple answer is

```
conda install -c bioconda deepbgc

```

conda will then install the correct version of scikit-learn automatically. If doing this it important it is in its own (new) `conda` environment however.

If the above doesn't work further stuff can be done, but the 'simple answer' should work.

---

@Joe\_G says its hanging ... on `solving environment` just leave it ... you will then see why,

```
The following NEW packages will be INSTALLED:

  _tflow_select      pkgs/main/osx-64::_tflow_select-2.3.0-mkl
  absl-py            pkgs/main/noarch::absl-py-0.15.0-pyhd3eb1b0_0
  appdirs            pkgs/main/noarch::appdirs-1.4.4-pyhd3eb1b0_0
  astor              pkgs/main/osx-64::astor-0.8.1-py37hecd8cb5_0
  biopython          pkgs/main/osx-64::biopython-1.76-py37h1de35cc_0
  blas               pkgs/main/osx-64::blas-1.0-mkl
  bottleneck         pkgs/main/osx-64::bottleneck-1.3.5-py37h67323c0_0
  brotli             pkgs/main/osx-64::brotli-1.0.9-hca72f7f_7
  brotli-bin         pkgs/main/osx-64::brotli-bin-1.0.9-hca72f7f_7
  c-ares             pkgs/main/osx-64::c-ares-1.18.1-hca72f7f_0
  ca-certificates    pkgs/main/osx-64::ca-certificates-2022.07.19-hecd8cb5_0
  certifi            pkgs/main/osx-64::certifi-2022.9.14-py37hecd8cb5_0
  cycler             pkgs/main/noarch::cycler-0.11.0-pyhd3eb1b0_0
  dataclasses        pkgs/main/noarch::dataclasses-0.8-pyh6d0b6a4_7
  deepbgc            bioconda/noarch::deepbgc-0.1.23-py_0
  fftw               pkgs/main/osx-64::fftw-3.3.9-h9ed2024_1
  fonttools          pkgs/main/noarch::fonttools-4.25.0-pyhd3eb1b0_0
  freetype           pkgs/main/osx-64::freetype-2.11.0-hd8bbffd_0
  gast               pkgs/main/noarch::gast-0.5.3-pyhd3eb1b0_0
  giflib             pkgs/main/osx-64::giflib-5.2.1-haf1e3a3_0
  grpcio             pkgs/main/osx-64::grpcio-1.42.0-py37ha29bfda_0
  h5py               pkgs/main/osx-64::h5py-3.7.0-py37h4a1dd59_0
  hdf5               pkgs/main/osx-64::hdf5-1.10.6-h10fe05b_1
  hmmer              bioconda/osx-64::hmmer-3.3.2-h9722bc1_2
  importlib-metadata pkgs/main/osx-64::importlib-metadata-4.11.3-py37hecd8cb5_0
  intel-openmp       pkgs/main/osx-64::intel-openmp-2021.4.0-hecd8cb5_3538
  joblib             pkgs/main/noarch::joblib-1.1.0-pyhd3eb1b0_0
  jpeg               pkgs/main/osx-64::jpeg-9e-hca72f7f_0
  keras              pkgs/main/osx-64::keras-2.2.4-0
  keras-applications pkgs/main/noarch::keras-applications-1.0.8-py_1
  keras-base         pkgs/main/osx-64::keras-base-2.2.4-py37_0
  keras-preprocessi~ pkgs/main/noarch::keras-preprocessing-1.1.2-pyhd3eb1b0_0
  kiwisolver         pkgs/main/osx-64::kiwisolver-1.4.2-py37he9d5cce_0
  lcms2              pkgs/main/osx-64::lcms2-2.12-hf1fd2bf_0
  lerc               pkgs/main/osx-64::lerc-3.0-he9d5cce_0
  libbrotlicommon    pkgs/main/osx-64::libbrotlicommon-1.0.9-hca72f7f_7
  libbrotlidec       pkgs/main/osx-64::libbrotlidec-1.0.9-hca72f7f_7
  libbrotlienc       pkgs/main/osx-64::libbrotlienc-1.0.9-hca72f7f_7
  libcxx             pkgs/main/osx-64::libcxx-14.0.6-h9765a3e_0
  libdeflate         pkgs/main/osx-64::libdeflate-1.8-h9ed2024_5
  libffi             bioconda/osx-64::libffi-3.2.1-1
  libgfortran        pkgs/main/osx-64::libgfortran-5.0.0-11_2_0_h246ff09_26
  libgfortran5       pkgs/main/osx-64::libgfortran5-11.2.0-h246ff09_26
  libpng             pkgs/main/osx-64::libpng-1.6.37-ha441bb4_0
  libprotobuf        pkgs/main/osx-64::libprotobuf-3.20.1-h8346a28_0
  libtiff            pkgs/main/osx-64::libtiff-4.4.0-h2ef1027_0
  libwebp            pkgs/main/osx-64::libwebp-1.2.2-h56c3ce4_0
  libwebp-base       pkgs/main/osx-64::libwebp-base-1.2.2-hca72f7f_0
  llvm-openmp        pkgs/main/osx-64::llvm-openmp-14.0.6-h0dcd299_0
  lz4-c              pkgs/main/osx-64::lz4-c-1.9.3-h23ab428_1
  markdown           pkgs/main/osx-64::markdown-3.3.4-py37hecd8cb5_0
  matplotlib-base    pkgs/main/osx-64::matplotlib-base-3.5.2-py37hfb0c5b7_0
  mkl                pkgs/main/osx-64::mkl-2021.4.0-hecd8cb5_637
  mkl-service        pkgs/main/osx-64::mkl-service-2.4.0-py37h9ed2024_0
  mkl_fft            pkgs/main/osx-64::mkl_fft-1.3.1-py37h4ab4a9b_0
  mkl_random         pkgs/main/osx-64::mkl_random-1.2.2-py37hb2f4e1b_0
  munkres            bioconda/noarch::munkres-1.0.7-py_1
  ncurses            pkgs/main/osx-64::ncurses-6.3-hca72f7f_3
  numexpr            pkgs/main/osx-64::numexpr-2.8.3-py37h2e5f0a9_0
  numpy              pkgs/main/osx-64::numpy-1.21.5-py37h2e5f0a9_3
  numpy-base         pkgs/main/osx-64::numpy-base-1.21.5-py37h3b1a694_3
  openssl            pkgs/main/osx-64::openssl-1.1.1q-hca72f7f_0
  packaging          pkgs/main/noarch::packaging-21.3-pyhd3eb1b0_0
  pandas             pkgs/main/osx-64::pandas-1.3.5-py37h743cdd8_0
  pillow             pkgs/main/osx-64::pillow-9.2.0-py37hde71d04_1
  pip                pkgs/main/osx-64::pip-22.2.2-py37hecd8cb5_0
  prodigal           bioconda/osx-64::prodigal-2.6.3-ha5712d3_4
  protobuf           pkgs/main/osx-64::protobuf-3.20.1-py37he9d5cce_0
  pyparsing          pkgs/main/osx-64::pyparsing-3.0.9-py37hecd8cb5_0
  python             pkgs/main/osx-64::python-3.7.7-hc70fcce_0_cpython
  python-dateutil    pkgs/main/noarch::python-dateutil-2.8.2-pyhd3eb1b0_0
  pytz               pkgs/main/osx-64::pytz-2022.1-py37hecd8cb5_0
  pyyaml             pkgs/main/osx-64::pyyaml-6.0-py37hca72f7f_1
  readline           pkgs/main/osx-64::readline-8.1.2-hca72f7f_1
  scikit-learn       pkgs/main/osx-64::scikit-learn-1.0.2-py37hae1ba45_1
  scipy              pkgs/main/osx-64::scipy-1.7.3-py37h214d14d_2
  setuptools         pkgs/main/osx-64::setuptools-63.4.1-py37hecd8cb5_0
  six                pkgs/main/noarch::six-1.16.0-pyhd3eb1b0_1
  sqlite             pkgs/main/osx-64::sqlite-3.39.3-h707629a_0
  tensorboard        pkgs/main/osx-64::tensorboard-1.14.0-py37h80053f4_0
  tensorflow         pkgs/main/osx-64::tensorflow-1.14.0-mkl_py37h085be34_0
  tensorflow-base    pkgs/main/osx-64::tensorflow-base-1.14.0-mkl_py37h5a24fda_0
  tensorflow-estima~ pkgs/main/noarch::tensorflow-estimator-1.14.0-py_0
  termcolor          pkgs/main/osx-64::termcolor-1.1.0-py37hecd8cb5_1
  threadpoolctl      pkgs/main/noarch::threadpoolctl-2.2.0-pyh0d69192_0
  tk                 pkgs/main/osx-64::tk-8.6.12-h5d9f67b_0
  typing_extensions  pkgs/main/osx-64::typing_extensions-4.3.0-py37hecd8cb5_0
  werkzeug           pkgs/main/noarch::werkzeug-2.0.3-pyhd3eb1b0_0
  wheel              pkgs/main/noarch::wheel-0.37.1-pyhd3eb1b0_0
  wrapt              pkgs/main/osx-64::wrapt-1.14.1-py37hca72f7f_0
  xz                 pkgs/main/osx-64::xz-5.2.6-hca72f7f_0
  yaml               pkgs/main/osx-64::yaml-0.2.5-haf1e3a3_0
  zipp               pkgs/main/osx-64::zipp-3.8.0-py37hecd8cb5_0
  zlib               pkgs/main/osx-64::zlib-1.2.12-h4dc903c_3
  zstd               pkgs/main/osx-64::zstd-1.5.2-hcb37349_0

Proceed ([y]/n)?

```

Its an enormous data pipeline. Everything is there, thats why it takes so long.
> 1 votes

---
Tags: python, deeptools, conda, scipy
---


# Question
Title: How to get started with louvain/leiden algorithm with UMAP in R
A collegue of mine recently suggested to try the louvain algorithm for clustering multiplex cytometry data. However, implementations of louvain are kind of rare in R. To my knowledge the only stand-alone implementation is included in the `igraph` package. So far i tried to use the adjacency matrix that the `umap` function returns (`uwot` package) to create a graph using `igraphs` `graph_from_adjacency_matrix` function. I can run the louvain algorithm on the graph, but the result is always a few thousand clusters with a hand-full if cells. changing the resolution parameter does not change anything.

If i remember correctly, `Seurats` `findClusters` function uses `louvain`, however i don't want to use PCA reduction before clustering, which is requiered in `Seurat` to find clusters.

Does anybody know of a different implementation of louvain or leiden algorithm that or an easier way to use `igraphs`functions? I obviously don't have to use the matrix umap returned, but what would be a better approach to make my data "louvain-able"?

# Answer
Both conda and PyPI have leiden clustering in Python which operates via `iGraph`

```
conda install -c conda-forge leidenalg

pip install leiden-clustering

```

Used via

```
import leidenalg as la
import igraph as ig

```

Example output

The docs are here. You will not need much Python to use it.

There is an entire Leiden package in R-cran here

Python implements Rs Leiden package here,

```
conda install -c conda-forge r-leiden

```
> 1 votes

---
Tags: r, seurat, clustering, umap
---


title: Bulk protein annotation with Prokka

text:
# Question
Title: Bulk protein annotation with Prokka
Can prokka be used to annotate a set of protein sequences (so bulk protein annotation), or does it require an entire genome for the annotation?

I've been looking through the documentation but haven't found a clear answer to this. I tested it, and the annotation of my protein sequences worked, but I wanted to check if the program expected a genome and that I can trust the output.

Thanks!

# Answer
> 1 votes
Prokka doesn't need whole genomes, I've definitely used it for smaller contigs before so I think you're fine!

---
Tags: functional-annotation
---


# Question
Title: Why are genetic algorithm chromosomes represented using 8 bit binary?
I am learning about genetic algorithms.

Why chromosomes are represented using 8 bit binary.

This arose when I apply crossover techniques to the bits.

Please help me to find the answer.Thanks in advance !

# Answer
> 1 votes
Not always. The representation of chromosomes in genetic algorithms is often dependent on the specific problem to be solved and number of variables to be found through the optimization process.

Not only you can chose a 16-bit representation, but integer-based representations are also used as mentioned here.

At the computation level, using a binary representation usually facilitates cross over operations on the bit arrays. 8 bits are compact while already give you a range of 0-255 different values for a given variable, thus provides a first option to go unless you need a larger value space.

---
Tags: gene, gene-expression, algorithms, genetics, chromosomes
---


# Question
Title: Does the number of RNA reads per cell obtained from the 10X scRNA experiment depend on amount of mRNA in given cell?
As we know, the amount of RNA reads per cell obtained from 10X scRNA experiment vary between cells. I wonder if this is effect of technical issues or does the number of RNA reads per cell obtained from the 10X scRNA experiment depend on the amount of mRNA in that cell? For the latter to be the case, we would have to read almost all the mRNA fragments in a cell, which I don't think is likely. Thank you for your replies in advance!

# Answer
> 1 votes
I am inclined to agree with you that the reads-per-cell is dependent on the cell type and cell cycle state. In the datasets I have analyzed, I observed that immune cells tend to have far fewer reads than fibroblasts, which have fewer reads than larger epithelial cell types.

EDIT: I have also observed that macrophage populations generally have more RNA than lymphocytes. Granulocytes have even less if you even capture them in your sample.

# Answer
> 1 votes
The number of reads per cell for any given experiment is largely technical. It depends on sequencing depth. 10X recommend 20K reads per cell for the v3 chemistry (and 50K reads per cell for v2)

---
Tags: ngs, scrnaseq, seurat, single-cell, 10x
---


# Question
Title: Nested variable expansion in WDL code block
I am developing a gatk pipeline for variant calling and want to include trimmomatic as a preprocessing step. I would like the user to specify a program file of parameters to feed into trimmomatic, and a file of adapter sequences to provide to the ILLUMINACLIP parameter. My goal is to include a variable (“~{adapter\_file}” below) in the program file that will expand to include the adapter filename, like so:

Adapter file:

```
>IL_Multiplex_five_prime_fwd
AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT
>IL_Multiplex_three_prime_rev
AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC

```

Program file:

```
ILLUMINACLIP:~{adapter_file}:2:30:10:7:keepBothReads
MINLEN:36
SLIDINGWINDOW:4:15
LEADING:3
TRAILING:3

```

In this way, the user can change the adaptor file without having to alter the program file each time. However, I can’t seem to get the ~{adapter\_file} variable to expand inside the ~{trim\_program\_array} variable. Here’s a simplified version of my trimmomatic WDL file:

Sample workflow:

```
version 1.0

Workflow trim {
    input {
        File fastq1 = “R1.fastq.gz”
        File fastq2 = “R2.fastq.gz”
        File program_file = “trim_program.txt”
        File adapter_file = “adapters.fa”
    }

    Array[String] trim_program_array = read_lines(program_file)

    call trimmomatic {
        input:
            fastq1 = fastq1,
            fastq2 = fastq2,
            trim_program_array = trim_program_array,
            adapter_file = adaptor_file
    }
}

Task trimmomatic {
    input {
        File fastq1
        File fastq2
        File outfile1 = “trimmed_R1.fastq.gz”
        File outfile2 = “trimmed_R2.fastq.gz”
        Array[String] trim_program_array
        File adapter_file
     }

     command {
        TrimmomaticPE \
        ~{fastq1} ~{fastq2} \
        ~{outfile1} ~{outfile2} \
        ~{sep = “ “ trim_program_array}
     }

     output {
        File out1 = outfile1
        File out2 = outfile2
     }
}

```

My goal is for the first element of ~{sep = " " trim\_program\_array} (aka the first line of trim\_program.txt) to evaluate to “ILLUMINACLIP:adapters.fa:2:30:10:7:keepBothReads” rather than “ILLUMINACLIP:~{adapter\_file}:2:30:10:7:keepBothReads”. I know this is because ~{adapter\_file} is nested inside trim\_program\_array\[0\], so it’s only expanding the “outer” variable. How can I convince WDL to expand the variable nested inside the first element of that array before executing the bash command?

# Answer
Putting the placeholder inside of an input file won't work, since the contents will just be turned into a string when read with `read_lines` and won't get treated like WDL code. Since the lines from the `program_file` are simply getting pasted into the command, you can move the contents into the command section. In a placeholder (`~{}`), adding a string to an optional value will result in an empty string if that optional value is not supplied (see the spec here). So you can use optional inputs to make the different lines non-mandatory:

```
task trimmomatic {
    input {
        File fastq1
        File fastq2
        String outfile1 = “trimmed_R1.fastq.gz”
        String outfile2 = “trimmed_R2.fastq.gz”
        File? adapter_file
        Int? minlen
        String? slidingwindow
        Int? leading
        Int? trailing
     }

     command {
        TrimmomaticPE \
        ~{fastq1} ~{fastq2} \
        ~{outfile1} ~{outfile2} \
        ~{"ILLUMINACLIP:" + adapter_file + ":2:30:10:7:keepBothReads"} \
        ~{"MINLEN:" + minlen} \
        ~{"SLIDINGWINDOW:" + slidingwindow} \
        ~{"LEADING:" + leading} \
        ~{"TRAILING:" + trailing}
     }

     output {
        File out1 = outfile1
        File out2 = outfile2
     }
}

```

If, say, `adapter_file` were to be omitted in the inputs, then the produced command should look something like this:

```
TrimmomaticPE \
R1.fq R2.fq \
trimmed_R1.fastq.gz trimmed_R2.fastq.gz \
 \
MINLEN:36 \
SLIDINGWINDOW:4:15 \
LEADING:3 \
TRAILING:3

```

Also note that the output file paths in the input section should be `String`s. Since they don't exist yet, the execution engine will complain that it can't find them if they are `File`s

Ps. I'm not familiar with the actual trimmomatic command options, so please do check that this produces the desired command. I suspect that the `:2:30:10:7:keepBothReads` bit could be split up into additional inputs as well.
> 1 votes

---
Tags: wdl
---


# Question
Title: How to map PDB chains to Uniprot IDs using API services
I have a lot of PDB IDs and I need to get uniprot fasta sequences of these PDB IDs special chains by API services. For example, imagine that I need to get fasta sequence of '1kf6' 'A' chain. The uniprot entry (Accession) for this '1kf6' chain is 'P00363'. If I find a way to get this entry, I can simply retrieve its fasta sequence using bioservices package:

```
from bioservices import *
u = UniProt()
sequence = u.retrieve("P00363","fasta")
print(sequence)

```

The main problem is that I don't know how can I find this uniprot entry("P00363" here) for each chains in PDB IDs(here I imagine that I don't know this "P00363" is the uniprot entry for '1kf6' 'A' chain and I'm trying to find it using API services). I tried to use ID mapping. In ID mapping cross-references, the chain name for each entry is specified and if I can add a cross-reference column to u.search('1kf6', columns="id,entry name,length, genes,pdb"), I can simply find the entry related to my favourite chain. The following code returns some entries each related to one of '1kf6' chains. The problem is that I don't know exactly which entry is related to 'A' chain that here is my favorite chain. In the following code, id searches for entry and other parameters search for sth else related to them:

```
from bioservices import *
u = UniProt()
res = u.search("1kf6",columns="id,entry name,length, genes")
print(res)

```

Each one of the following words can be used and added to columns part:

```
'id', 'entry name', 'genes', 'genes(PREFERRED)', 'genes(ALTERNATIVE)', 'genes(OLN)', 'genes(ORF)', 'organism', 'organism-id', 'protein names', 'proteome', 'lineage(ALL)', 'lineage-id', 'virus hosts', 'fragement', 'sequence', 'length', 'mass', 'encodedon', 'comment(ALTERNATIVE PRODUCTS)', 'comment(ERRONEOUS GENE MODEL PREDICTION)', 'comment(ERRONEOUS INITIATION)', 'comment(ERRONEOUS TERMINATION)', 'comment(ERRONEOUS TRANSLATION)', 'comment(FRAMESHIFT)', 'comment(MASS SPECTROMETRY)', 'comment(POLYMORPHISM)', 'comment(RNA EDITING)', 'comment(SEQUENCE CAUTION)', 'feature(ALTERNATIVE SEQUENCE)', 'feature(NATURAL VARIANT)', 'feature(NON ADJACENT RESIDUES)', 'feature(NON STANDARD RESIDUE)', 'feature(NON TERMINAL RESIDUE)', 'feature(SEQUENCE CONFLICT)', 'feature(SEQUENCE UNCERTAINTY)', 'version(sequence)', 'domains', 'domain', 'comment(DOMAIN)', 'comment(SIMILARITY)', 'feature(COILED COIL)', 'feature(COMPOSITIONAL BIAS)', 'feature(DOMAIN EXTENT)', 'feature(MOTIF)', 'feature(REGION)', 'feature(REPEAT)', 'feature(ZINC FINGER)', 'ec', 'comment(ABSORPTION)', 'comment(CATALYTIC ACTIVITY)', 'comment(COFACTOR)', 'comment(ENZYME REGULATION)', 'comment(FUNCTION)', 'comment(KINETICS)', 'comment(PATHWAY)', 'comment(REDOX POTENTIAL)', 'comment(TEMPERATURE DEPENDENCE)', 'comment(PH DEPENDENCE)', 'feature(ACTIVE SITE)', 'feature(BINDING SITE)', 'feature(DNA BINDING)', 'feature(METAL BINDING)', 'feature(NP BIND)', 'feature(SITE)', 'go', 'go(biological process)', 'go(molecular function)', 'go(cellular component)', 'go-id', 'interpro', 'interactor', 'comment(SUBUNIT)', 'citation', 'citationmapping', 'created', 'last-modified', 'sequence-modified', 'version(entry)', '3d', 'feature(BETA STRAND)', 'feature(HELIX)', 'feature(TURN)', 'comment(SUBCELLULAR LOCATION)', 'feature(INTRAMEMBRANE)', 'feature(TOPOLOGICAL DOMAIN)', 'feature(TRANSMEMBRANE)', 'annotation score', 'score', 'features', 'comment(CAUTION)', 'comment(TISSUE SPECIFICITY)', 'comment(GENERAL)', 'keywords', 'context', 'existence', 'tools', 'reviewed', 'feature', 'families', 'subcellular locations', 'taxonomy', 'version', 'clusters', 'comments', 'database', 'keyword-id', 'pathway', 'score'

```

I used some of the above words in columns part but they couldn't retrieve chain name. Anyone know which one of the above parameters can give me the chain name for this PDB ID? or anyone has any other idea for finding uniprot entry(UniprotKB AC) related to each PDB ID chain by API services?

# Answer
I found the answer and I thought it may be useful for others. we can use pypdb package for that purpose (it works in unix systems):

```
import pypdb

all_info = pypdb.get_all_info('1kf6')
print(all_info)

```

the output is as follows:

```
{'polymer': [{'@entityNr': '1', '@length': '602', '@type': 'protein', '@weight': '66057.6', 'chain': [{'@id': 'A'}, {'@id': 'M'}], 'Taxonomy': {'@name': 'Escherichia coli', '@id': '562'}, 'macroMolecule': {'@name': 'Fumarate reductase flavoprotein subunit', 'accession': {'@id': 'P00363'}}, 'polymerDescription': {'@description': 'FUMARATE REDUCTASE FLAVOPROTEIN'}, 'enzClass': {'@ec': '1.3.5.4'}}, {'@entityNr': '2', '@length': '243', '@type': 'protein', '@weight': '27021.9', 'chain': [{'@id': 'B'}, {'@id': 'N'}], 'Taxonomy': {'@name': 'Escherichia coli', '@id': '562'}, 'macroMolecule': {'@name': 'Fumarate reductase iron-sulfur subunit', 'accession': {'@id': 'P0AC47'}}, 'polymerDescription': {'@description': 'FUMARATE REDUCTASE IRON-SULFUR PROTEIN'}, 'enzClass': {'@ec': '1.3.5.1'}}, {'@entityNr': '3', '@length': '130', '@type': 'protein', '@weight': '14898.8', 'chain': [{'@id': 'C'}, {'@id': 'O'}], 'Taxonomy': {'@name': 'Escherichia coli', '@id': '562'}, 'macroMolecule': {'@name': 'Fumarate reductase subunit C', 'accession': {'@id': 'P0A8Q0'}}, 'polymerDescription': {'@description': 'FUMARATE REDUCTASE 15 KDA HYDROPHOBIC PROTEIN'}}, {'@entityNr': '4', '@length': '119', '@type': 'protein', '@weight': '13118.9', 'chain': [{'@id': 'D'}, {'@id': 'P'}], 'Taxonomy': {'@name': 'Escherichia coli', '@id': '562'}, 'macroMolecule': {'@name': 'Fumarate reductase subunit D', 'accession': {'@id': 'P0A8Q3'}}, 'polymerDescription': {'@description': 'FUMARATE REDUCTASE 13 KDA HYDROPHOBIC PROTEIN'}}], 'id': '1KF6'}

```

This output gives all the information related to all chains and we can simply find uniprotkb AC of each chain in this output. I also found another answer for this question. This answer can be found here and the following code:

```
import requests
from xml.etree.ElementTree import fromstring

pdb_id = '1kf6.A'
pdb_mapping_url = 'http://www.rcsb.org/pdb/rest/das/pdb_uniprot_mapping/alignment'
uniprot_url = 'http://www.uniprot.org/uniprot/{}.xml'

def get_uniprot_accession_id(response_xml):
    root = fromstring(response_xml)
    return next(
        el for el in root.getchildren()[0].getchildren()
        if el.attrib['dbSource'] == 'UniProt'
    ).attrib['dbAccessionId']

def get_uniprot_protein_name(uniport_id):
    uinprot_response = requests.get(
        uniprot_url.format(uniport_id)
    ).text
    return fromstring(uinprot_response).find('.//{http://uniprot.org/uniprot}recommendedName/{http://uniprot.org/uniprot}fullName'
    ).text

def map_pdb_to_uniprot(pdb_id):
    pdb_mapping_response = requests.get(
        pdb_mapping_url, params={'query': pdb_id}
    ).text
    uniprot_id = get_uniprot_accession_id(pdb_mapping_response)
    uniprot_name = get_uniprot_protein_name(uniprot_id)
    return {
        'pdb_id': pdb_id,
        'uniprot_id': uniprot_id,
        'uniprot_name': uniprot_name
    }

print (map_pdb_to_uniprot(pdb_id))

```

The output of this code is as follows:

```
{'pdb_id': '1kf6.A', 'uniprot_id': 'P00363', 'uniprot_name': 'Fumarate reductase flavoprotein subunit'}

```
> 3 votes

# Answer
An alternative solution involves the use of the SIFTS resource from the EBI. I wrote a python parser for their xml format, which will give you a mapping between a pdb and a uniprot entry at the single residue level.

```
$ curl --silent "ftp://ftp.ebi.ac.uk/pub/databases/msd/sifts/xml/1kf6.xml.gz" | gunzip | python parse_sifts.py
1kf6    A       M       0       P00363  M       1
1kf6    A       Q       1       P00363  Q       2
1kf6    A       T       2       P00363  T       3
1kf6    A       F       3       P00363  F       4
1kf6    A       Q       4       P00363  Q       5
1kf6    A       A       5       P00363  A       6
1kf6    A       D       6       P00363  D       7
1kf6    A       L       7       P00363  L       8
1kf6    A       A       8       P00363  A       9
1kf6    A       I       9       P00363  I       10
[...]

```

The returned fields are as follows:

* pdb ID
* chain
* pdb residue
* pdb position
* uniprot ID
* uniprot residue
* uniprot position

In the above snippet, just put your PDB ID of choice to retrieve the data you need.
> 2 votes

# Answer
While not as nice as using SIFT. You can use the uniprot sparql endpoint as well to get mappings.

```
PREFIX up: <http://purl.uniprot.org/core/>
SELECT
  ?pdbId ?chain ?pdbChain ?uniprot
WHERE
{
  # A space separated list of pairs of PDB identifiers and the chain code.
  VALUES(?pdbId ?pdbChain) { ('6VXC' 'A') ('1BG3' 'B') }

  # Make an IRI out of the pdbId
  BIND(iri(concat('http://rdf.wwpdb.org/pdb/', ?pdbId)) AS ?pdb)

  # Map to UniProt entries
  ?uniprot rdfs:seeAlso ?pdb .
  ?pdb up:database <http://purl.uniprot.org/database/PDB> ;
   up:chainSequenceMapping ?chainSm .
  ?chainSm up:chain ?chainsPlusRange .

  # Extract the list of chains from the text representation.
  BIND(STRBEFORE(?chainsPlusRange, '=') AS ?chain)

  # Filter those that match.
  FILTER(CONTAINS(?chain, ?pdbChain))
}

```
> 2 votes

# Answer
Check out this file on SIFT. ftp.ebi.ac.uk/pub/databases/msd/sifts/csv/pdb\_chain\_uniprot.csv It contains residue-wise mapping of PDB chains with uniprot ids.

```
# 2020/04/12 - 21:03 | PDB: 15.20 | UniProt: 2020.02
PDB,CHAIN,SP_PRIMARY,RES_BEG,RES_END,PDB_BEG,PDB_END,SP_BEG,SP_END
101m,A,P02185,1,154,0,153,1,154
102l,A,P00720,1,40,1,40,1,40
102l,A,P00720,42,165,41,None,41,164
.
.

```
> 1 votes

---
Tags: fasta, python, pdb, api
---


# Question
Title: How to promote assemblies into genomes in NCBI?
Note: I've never submitted an assembly/genome to NCBI, so excuse if my perspective is flawed.

I'm working with *Drosophila subobscura*. (spring fruit fly)

I see here https://www.ncbi.nlm.nih.gov/data-hub/genome/?taxon=7241 that there are at least 2 assemblies of better quality than the only available genome (I checked checked N50 and number of scaffolds).

However, these assemblies weren't promoted to genomes.

**Questions**

1. What is required to consider them genomes?
2. Is it the scaffold to chromosome notation that is missing? Or also the annotations?
3. How would I promote the better assemblies into genomes?
4. And also, if I wanted to annotate one of the better assemblies, how would I proceed?

# Answer
> 1 votes
The reason is simple, insect genomes are often heavily fragmented and historically relied on wet-lab experimentation - polytene chromosome hybridisation to correcty assemble the chromosomes. Thus you probe the various contigs, prep the contigs often from salivary glands, stick them under a microscope and map the location of the probes. Hence, the reluctance to assign a 'genome' based on contigs because the chromosomes cannot be assembled. This is particularly true for a very close relative to the premier model of the Insecta *D. melanogaster*. A lot of work will have gone into the correct chromosome annotation of this genome. Thus, you can't simply 'promote' a genome without further genomics analysis and it will be a community decision based on the analysis (i.e. reviewers).

Bioinformatic assemblies often get it right. What I would do is use a reference genome assembly based on the nearest ancestor, which you'll know more than me, I'm pretty sure its not *D. melanogaster*. This is within the context of a comparative genomic analysis between the species/subspecies associated with *D. subobscura* and then say to NCBI 'thats a genome'. You then have to publish the analysis and submit your assembled genome to NCBI. Its not trivial, but if you 'live next door' to one of the genetic models of eukaryotes (well metazoa) the standards of annotation and assembly will be high.

---
Tags: genome, assembly, phylogenetics
---


# Question
Title: Warning message related to tensorflow package when trying to run deepbgc bioinformatics tool
I have come across some related issues pertaining to this tensorflow warning but none related specifically to deepbgc

I have created a conda environment and have installed onto it the bioinformatics analysis tool deepbgc. Now the program runs but it does output some warning messages related to the tensorflow package. When I run the "deepbgc info" command through terminal, I get the following output messages related to the tensor flow package:

```
Using TensorFlow backend.
WARNING 03/10 12:55:05   From /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:74: The name tf.get_default_graph is deprecated. Please use tf.compat.v1.get_default_graph instead.

WARNING 03/10 12:55:05   From /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:517: The name tf.placeholder is deprecated. Please use tf.compat.v1.placeholder instead.

WARNING 03/10 12:55:05   From /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:4138: The name tf.random_uniform is deprecated. Please use tf.random.uniform instead.

WARNING 03/10 12:55:06   From /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:133: The name tf.placeholder_with_default is deprecated. Please use tf.compat.v1.placeholder_with_default instead.

WARNING 03/10 12:55:06   From /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:3445: calling dropout (from tensorflow.python.ops.nn_ops) with keep_prob is deprecated and will be removed in a future version.
Instructions for updating:
Please use `rate` instead of `keep_prob`. Rate should be set to `rate = 1 - keep_prob`.
WARNING 03/10 12:55:06   From /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:174: The name tf.get_default_session is deprecated. Please use tf.compat.v1.get_default_session instead.

WARNING 03/10 12:55:06   From /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:181: The name tf.ConfigProto is deprecated. Please use tf.compat.v1.ConfigProto instead.

WARNING 03/10 12:55:06   From /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:186: The name tf.Session is deprecated. Please use tf.compat.v1.Session instead.

WARNING 03/10 12:55:06   From /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:190: The name tf.global_variables is deprecated. Please use tf.compat.v1.global_variables instead.

WARNING 03/10 12:55:06   From /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:199: The name tf.is_variable_initialized is deprecated. Please use tf.compat.v1.is_variable_initialized instead.

WARNING 03/10 12:55:06   From /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:206: The name tf.variables_initializer is deprecated. Please use tf.compat.v1.variables_initializer instead.

```

Could this be related to the version of tensorflow? Currently, I have version 1.15.4 installed onto the conda environment. Are these error messages significant and could they impact the data running through deepbgc? Or can they mostly be ignored? Any guidance appreciated, thanks!

# Answer
1. The 'blunderbus' way to do this is as follows (not recommended in this instance :-) ),

```
    conda update --all
    conda update -n base conda
    conda update -n myenv --all

```
2. The refined way to do this is

```
    conda update -n myenv tensorflow # see `Best guess, notes

```

In theory `conda` was supposed to sort out the deprecation, because all the versions were supposed to be harmonious. I'm a bit puzzled (see below). What might have happened is `tensorflow` has been upgraded (or `keras`) against the version the developer had written it for.

Here `myenv` is the name you gave that given conda environment.

In your case, I would look at point 2, because the environment is very specific to your project and ensures reproducibility.

*Note* A purist would specify the exact channel `-c` that `tensorflow` is found in, or whatever dependency. This actually can be very important because there can be name conflicts, obviously not for `tensorflow` because I'm pretty sure Google have the name registered.

The reason I do point 1 is the problem might a dependency of `tensorflow` and you'd have to figure it out, but again see 'notes'.

---

**Best guess**

In this particular example,

```
conda update -n myenv keras

```

Might be the ticket, because what I think is going on is its not actually `tensorflow` but the `keras` sitting above it.

The alternative (probably better) guess is ....

```
conda update -n myenv -c bioconda deepbgc

```

This is probably what I would try first.

**Note** I'm surprised that `conda`/`mamba` has not resolved this issue on installation because thats exactly what is it is supposed to do. In other words its supposed to make a perfect installation of `deepbgc` precisely as the developer had intended it. There is a risk that updating might not be the solution ... and solution is a new `conda`/`mamba` installation (see note 2).

**Note 2** Watch for the `Python` version, you are using 3.7 ... that might be an issue. If the developer had written this for \>`3.7` that could cause problems. Different Python versions have different methods and if one of the 'new' methods was being used that would be an issue. Ultimately your package/pipeline is really complex and the more stuff the more chance stuff doesn't fit together. We're currently on `3.10` and `3.9` is really stable.

**BTW** `3.11` will be a major upgrade of `Python3` and will be the 'goto' version.

---

I looked at `deepbgc` its built on `Python3.7`, I don't know why because it was released in 2019. My advice is start again with a new install and follow the instructions exactly as on Github. The order of the build is really important (omitting pip use `conda`/`mamba` instead). It doesn't seem to be very `conda` but thats how the developer did it. In theory the whole pipeline could have been placed on `conda`. Instructions here https://github.com/Merck/deepbgc

The key issue is that it is built on an old build of Python, with a lot of dependencies which have substancially changed. There is concern about ...

```
conda create -n deepbgc python=3.7 hmmer prodigal

```

because the author is not specifying the versions of hmmer or prodigal. These should have been built into the `conda` upload.
> 1 votes

---
Tags: python, deeptools, conda, tensorflow
---


# Question
Title: How to add N bases to PyRange object (i.e., shifting PyRanges coordinates)
I would like to add `n` (e.g., 1k) bases to a PyRange object (i.e., shifting PyRanges coordinates). How can I do that?

**Input:**

```
import pyranges as pr
gr = pr.random(3)

# gr:
# +--------------+-----------+-----------+--------------+
# | Chromosome   |     Start |       End | Strand       |
# | (category)   |   (int32) |   (int32) | (category)   |
# |--------------+-----------+-----------+--------------|
# | chr1         |  10000000 |  20000000 | -            |
# | chr2         |  20000000 |  30000000 | +            |
# | chr22        |  30000000 |  40000000 | -            |
# +--------------+-----------+-----------+--------------+

```

**Desired output:**

```
# magical command[gr + 1000]:
# +--------------+-----------+-----------+--------------+
# | Chromosome   |     Start |       End | Strand       |
# | (category)   |   (int32) |   (int32) | (category)   |
# |--------------+-----------+-----------+--------------|
# | chr1         |  10001000 |  20001000 | -            |
# | chr2         |  20001000 |  30001000 | +            |
# | chr22        |  30001000 |  40001000 | -            |
# +--------------+-----------+-----------+--------------+

```

# Answer
You can use the `.extend()` method. Starting with your data above, `gr.extend({"3": -10000, "5": 10000})` would generate the following:

```
+--------------+-----------+-----------+--------------+
| Chromosome   |     Start |       End | Strand       |
| (category)   |   (int32) |   (int32) | (category)   |
|--------------+-----------+-----------+--------------|
| chr1         |  10010000 |  20010000 | -            |
| chr2         |  19990000 |  29990000 | +            |
| chr22        |  30010000 |  40010000 | -            |
+--------------+-----------+-----------+--------------+

```

I don't know a way of to apply this "range extension" in a strand specific way (to be able to generate the desired output exactly) but if this is really needed, you should be able to generate 2 `gr` objects, one for each strand, do the manipulation and bring these together.
> 2 votes

---
Tags: python, pyranges
---


# Question
Title: How is an X chromosome encoded into a fasta string?
A human "X" chromosome has a centromere and two "identical" chromatids. If the chromatids are not identical, this fact is not assembled, correct?

The fasta string for a chromomsome is actually representing just one of the chromatids, correct?

So from my understanding, for diploids, the fasta chromosome is actually a mix of 4 chromatids: 2chromatids * 2 pairs of chromosomes into one string. Which for some heterozygous nucleotides (same percentage of each allele) will just be randomly selected from the reads as the correct one.

(unless we are talking about some methods like **10x** that can distinguish between chromosomes pairs)

I was thinking which side of the X is read first into the fasta string, and I came to the above conclusion.

# Answer
So, the two chromatids you see only during meta-phase of replication. In more "regular" state of a cell, there are only the two haplotypes (in the case of a diploid organism)

Well, and regardless of the ploidy, we usually aim for a haploid (monoploid) representation, so it is usually exactly as you say - an arbitrary one of the two states on each heterozygous site is chosen.

However, there are also approaches to phase genomes (to sort which alleles are on which haplotypes). That is possible not just with 10x, but also with long reads (to certain extend) or using trio-binning \- a strategy to sequence two parents and the offspring.

However, even phased genomes still contain the reference (fasta) in a haploid state and usually provide the alternative haplotype as a second fasta file.

For the record, these apply to all chromosomes, not just the X.
> 2 votes

---
Tags: assembly
---


# Question
Title: Remote blastn - what is breaking my (bash) loop?
I am trying to search the NCBI non-redundant database for sequences similar to a few other (~40) sequences that I already have.

So I've tried running blastn remotely, and looping through multiple queries. I've tried it both as a job (via slurm) and from the command line.

In addition, for every query, I am trying to run 3 searches with slightly different parameters (in WebBlast, I noticed that I recover more sequences of interest this way).

I am trying to dump the output data in a single file (stuff like duplicates can be sorted out later).

This is my code:

```
for i in /directory/*.fa;do #this is the directory containing the queries
j=${i#"/directory/"}
k=${j%".fa"} #getting a short version of the name for the log file
echo "$k" >> "log.log"
blastn -query "$i" -db "nt" -remote -entrez_query 700:1000[slen] -perc_identity 94 -outfmt "6 sacc sseqid sseq evalue bitscore length" >> "results.tsv"
blastn -query "$i" -db "nt" -remote -entrez_query 1000:1400[slen] -perc_identity 94 -outfmt "6 sacc sseqid sseq evalue bitscore length" >> "results.tsv"
blastn -query "$i" -db "nt" -remote -perc_identity 94 -outfmt "6 sacc sseqid sseq evalue bitscore length" >> "results.tsv"
done

```

However, it...doesn't seem to work? Either only the first search for the first query is completed, and then the code does nothing else, or nothing happens. (in a rare case, the first two searches worked) There seems to be no clear distinction between 'it works a tiny bit' and 'it doesn't work at all', which leads me to think it might be something outside the code, e.g. connection issues.

Has anyone faced this before, and would you know how to fix it?

**Update** I tried to run a single search as well. As well as adding the waiting period. Right now it doesn't work (same issue, no error message, just sort of appearing to run without results). In the link that Scot provided below, from the NCBI, it was mentioned 'We will move searches of users who submit more than 100 searches in a 24 hour period to a slower queue, or, in extreme cases, will block the requests.' I haven't run \>100 searches yesterday (was planning to before I knew this, but that's a different issue).

I think I'll limit the searches as recommended and try to run it from a different computer later on.

# Answer
> 2 votes
I gave it a whirl using a random bit of SARS-CoV-2

```
blastn -query out.fa -db "nt" -remote  -outfmt "6 sacc sseqid sseq evalue bitscore length" >> "results.tsv"

blastn -query out.fa -db "nt" -remote -perc_identity 94 -outfmt "6 sacc sseqid sseq evalue bitscore length" >> "results2.tsv"

```

Both gave a 500 line output - the limit of the blast specified.

One possibility is to place a `sleep 10` line between each `blastn` command (3 additional lines).

The rationale is specifications NCBI set:

> Do not contact the server more often than **once every 10 seconds**. Do not poll for any single RID more often than once a minute. Use the URL parameter email and tool, so that the NCBI can contact you if there is a problem. Run scripts weekends or between 9 pm and 5 am Eastern time on weekdays if more than 50 searches will be submitted.

---

The suggestions haven't worked, so I'd recommend trying `sleep 60` between each `blastn`command on the basis this might be an issue with `RID`. If that fails the fall back is to run the job overnight US time via `cron tab`, which I'll post the code if it is required. Thus,

> Run scripts weekends or between 9 pm and 5 am Eastern time on weekdays if more than 50 searches will be submitted.

The final thing to add is that it is a fairly long job so connectivity might be an issue. If this was the case then it might also be resolved by running the script overnight in the OPs timezone, i.e. via `cron`

---

`sleep 60` works but is very slow.

* Agreed one option is to download the database and perform the Blast locally.
* Alternative use mmseqs2
* You could perform the Blast on the Cloud, like GCP (much faster).
* Personally I would look at performing over night runs when the network was not in use.

# Answer
> 2 votes
This is likely a throttling issue - see here for some details, specifically the one about contacting the NCBI server more than once every 10 seconds.

You may be able to submit a multiFASTA of all your queries to get around this issue, rather than individual FASTA sequences, or alternatively, add in a `sleep 10` command between your BLAST commands and at the end of each `for` loop iteration to throttle the requests.

---
Tags: blast, bash, loop
---


title: Mapping PDB ID + chain ID to UniProt ID

text:
# Question
Title: Mapping PDB ID + chain ID to UniProt ID
This question was addressed here.

Unfortunately, including `bioservices`, the solutions `pypdb` and `map_pdb_to_uniprot` do not work for now.

Besides, I do not need residue level mapping as one of the suggested solution.

Does the failure of known methods have anything to deal with the recent update of the UniProt database?

Please let me know if you know a workaround.

Thank you!

# Answer
> 0 votes
Please find the solution here: https://www.biostars.org/p/9540519/#9540826

---
Tags: pdb, uniprot
---


# Question
Title: What can be the bias of aligning paired-reads in a single-end mode?
I don't know if I'm in the right place but I have a technical problem to fix. I would have to align **paired-end reads** from Illumina sequencing to compare a normal genome with a tumor one.

When I align with **bwa-mem (default parameters) with paired-end** mode I get a bigWig file where the coverage of the two samples is quite similar. When I align with **bwa-mem (default parameters) with single-end** mode I get coverage holes (definitely does not align any read) in some specific regions. But this only happens in tumor genomes.

Is it possible that something is really going on in that region? Or is it just a bias? And if it's a bias, why in the normal condition, in that specific region, do they align?

I think I understand that by treating paired-end reads as single-end I lose the directionality of mate2 if I align mate 1. But I still can't understand why I have this alignment defect only in the tumor genome and not in the normal genome (I observe this phenomenon in all the samples analyzed, not just in one).

Could someone give me some explanation or recommend any checks to do?

**spoiler** \- these holes fall exactly in highly repetitive regions

# Answer
In principle one should simply not align paired-end reads as if they were single end reads, not only because you are not taking advantage of the forward/reverse information contained in the reads, but also you are not taking into account the expected distance between the two pairs, both of which result in better mapping rates/accuracy.

Having said that, and not recalling from the top of my head right now the relevant BWA-mem algorithmic implementation details, it is generally difficult to map (even) paired-reads to repetitive regions. Still there will be some rate of success. In tumor samples however, we could reasonably imagine that chromosomal rearrangements, INDELs, CNVs, etc. will make it even harder to confidently assign these reads to a location in the reference, since these extra variation will increase e.g. multiple mapping locations and mapping qualities for some fraction of the reads, thus resulting in regions with low or no coverage compared to non-tumor samples.

Repetitive regions might exacerbate the problem and get affected even further compared to other regions.
> 3 votes

# Answer
I tried to give me this explanation: with single-end alignment of paired-end reads, losing the directionality of the reads and also the expected fragment length, the R2 no longer maps if it maps the R1, but can map as a single reads anywhere, or don't map at all. I thought that therefore I don't see alignment (in tumors) in repetitive regions because of the high mutation rate I have in the tumor sample. While in normal I don't have this situation. However, aligning with paired-end reads, I see coverage in these regions again despite the high mutation rate because R1 and R2 help each other map in the best position, and so even if I have mutations I still see the alignment. Could it be appropriate to make a variant call despite being on repetitive regions?
> 1 votes

---
Tags: sequence-alignment, sam, fastq, bwa, somatic
---


# Question
Title: Gene expression changes in one condition
I want to look at an RNA-seq data-set on 1 condition to look for expressed genes. As there is only one condition, I can’t perform differential analysis.

The only thing I can think to do is set a threshold (i.e. all replicates must have \> 10 gene counts) then say these are expressed. But I don’t think this is a very good way of doing it.

Do you know of any better ways that we can do this?

# Answer
> 2 votes
If you've got enough replicates (ideally at least 6, but can be done with 3), you should be able to do a Variance-stabilised Transformation using DESeq2:

```
vsd <- vst(dds, blind=FALSE)
rld <- rlog(dds, blind=FALSE)
head(assay(vsd), 3)

```

More details on this can be found in the DESeq2 Vignette:

https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#extracting-transformed-values

# Answer
> 1 votes
*Answer from @james-hawley, converted from comment:*

In any high throughput sequencing experiment, there will be lowly expressed RNA molecules that are not captured by the library prep + sequencing. So defining what "is expressed" is somewhat arbitrary from RNA-seq data, regardless of how many replicates you have and what method you use. That said, you can control some of the noise in your RNA-seq data using a variance-stabilized transformation, like @gringer suggests, and look at the distribution of genes/transcripts. Some threshold like $rlog \ge x$ in at least $n$ samples is usually a good starting point (for some arbitrary $x$ and $n$).

# Answer
> 1 votes
I don't think there is a very good way to do what you want to do. I'd think some kind of RPM would be better than strictly raw counts; you have to account for overall library size. Or TPM might be better.

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


# Question
Title: Why does the Smith-Waterman algorithm have a complexity of $O(m^2 \times n)$?
I recently read that original the Smith-Waterman algorithm to do a local pairwise alignment of sequences has a complexity of $O(m^2 \times n)$, where $m$ and $n$ would be the lengths of the sequences with $m$ the potentially longer one. What I am wondering is where the additional complexity arises in the algorithm compared to the Needelmann-Wunsch algorithm of global alignment, which has complexity $O(m \times n)$.

Am I correct to assume the $m\times n$ for the Needelmann-Wunsch algorithm comes from building up the scoring matrix - whereas the traceback is negligible in complexity i.e. $O(m\times n + m)$? Hence we arrive at $O(m\times n)$.

Why is there a $m^2$ term for the Smith-Watermann algorithm?

# Answer
I think your source is wrong. You are correct to assume the brunt of its complexity comes from building the $m \*n$ matrix. The main differences are in the calculations per cell and the traceback, but these don't make $m ^2$. Wikipedia also gives an identical $O(m * n)$ as worst case for both algorithms.

edit: Aparently in this paper the algorithms were improved to go from $m^2n$ to $m\*n$
> 2 votes

---
Tags: sequence-alignment
---


title: Name this alternative INDEL format

text:
# Question
Title: Name this alternative INDEL format
I have recently been working with variants which are represented in a non-standard format. This alternative format uses a period instead of prefixing variants with a reference base. What is this format called?

There is no specific output format associated with it. It's part of some calling code I've inherited. I figured that this existed elsewhere. But, I could be wrong.

**Expected Format:**

Insertion:

```
pos=4
alt=ATTT
ref=A

```

Deletion:

```
pos=4
alt=A
ref=ATTT

```

**Alternative Format:**

Insertion:

```
pos=5
alt=TTT
ref=.

```

Deletion:

```
pos=5
alt=.
ref=TTT

```

# Answer
> 1 votes
The alternative deletion representation that you are describing, i.e. using `.` to represent alternate alleles that are not present, is allowed within the definition of the VCF format (in other words, this representation can be found in a standard format). An example of this is provided on the Wikipedia page for VCF:

```
#CHROM POS      ID         REF   ALT    QUAL  FILTER
20     1230237  .          T     .      47    PASS

```

While missing values are allowed for *ALT*, they are not typically used because they lead to ambiguous variant definitions.

For missing *REF* values, the VCF format specification has the following explanation:

> Each base must be one of A,C,G,T,N (case insensitive). Multiple bases are permitted. The value in the POS field refers to the position of the first base in the String. For simple insertions and deletions in which either the REF or one of the ALT alleles would otherwise be null/empty, the REF and ALT Strings must include the base before the event (which must be reflected in the POS field), unless the event occurs at position 1 on the contig in which case it must include the base after the event....

To work around this ambiguity, VCF files are commonly left-normalised so that the specified sequences and variants can be precisely placed within the reference genome. Left normalisation can be carried out using `bcftools norm`.

More details on left normalisation can be found here:

https://genome.sph.umich.edu/wiki/Variant\_Normalization#Left\_alignment

> * A variant is normalized if and only if it is parsimonious and left aligned.
> * A variant is parsimonious if and only if it is represented in as few nucleotides as possible without an allele of length 0.
> * A variant is left aligned if and only if it is no longer possible to shift its position to the left while keeping the length of all its alleles constant.

In other words, these are ***non-parsimonious variants***.

---
Tags: variant-calling, variants
---


# Question
Title: How to plot stacked bar chart using R showing mean with range and labelled values?
I am trying to plot a stacked bar graph with mean values from multiple repeats, also showing the range of the data and values labelled.

My data frame is below:

```
structure(list(Average.DB.per.FA = c(0, 0.5, 1, 1.5, 2, 0, 0.5,
1, 1.5, 2, 0, 0.5, 1, 1.5, 2, 0, 0.5, 1, 1.5, 2, 0, 0.5, 1, 1.5,
2, 0, 0.5, 1, 1.5, 2), condition = c("CDD", "CDD", "CDD", "CDD",
"CDD", "CDD", "CDD", "CDD", "CDD", "CDD", "CDD", "CDD", "CDD",
"CDD", "CDD", "STD", "STD", "STD", "STD", "STD", "STD", "STD",
"STD", "STD", "STD", "STD", "STD", "STD", "STD", "STD"), percentage = c(0,
0.807218011, 0.192781989, 0, 0, 0, 0.737968015, 0.262031985,
0, 0, 0, 0.739096101, 0.260903899, 0, 0, 0, 0.466751289, 0.346156739,
0.170617491, 0.016474482, 0, 0.462997344, 0.331833715, 0.183359127,
0.021809813, 0.012053568, 0.47623097, 0.307175237, 0.182459607,
0.022080618)), row.names = c(NA, -30L), class = "data.frame")

```

I tried to use the following code:

```
    ggplot(long_DF, aes(fill=Average.DB.per.FA, y=percentage, x=condition)) +
  geom_bar(position=position_dodge(), stat="identity")

```

But the graph I obtained does not have bars add up to 1:

I would like to have graphs like below (which I created by Prism using the same set of data but Prism does not allow me to label bars with numbers):

Is there any way to reproduce the graph by Prism using R and add values to each stacked bar?

Appreciate your help and thanks in advance.

# Answer
> 1 votes
It seems like this is missing a grouping variable in the aesthetics definition. The R help page for `position_dodge` mentions the following:

> Position\_dodge() requires the grouping variable to be be specified in the global or geom\_\* layer. Unlike position\_dodge(), position\_dodge2() works without a grouping variable in a layer. position\_dodge2() works with bars and rectangles, but is particulary useful for arranging box plots, which can have variable widths.

Bearing this in mind, setting `group=Average.DB.per.FA` seems to give the result you want:

```
long_DF %>%
  ggplot() +
  aes(fill=Average.DB.per.FA,
      y=percentage, x=condition, group=Average.DB.per.FA) +
  geom_col()

```

... although this is adding up to 3 rather than 1, as would be expected from looking at the category counts:

```
> long_DF %>%
+   group_by(condition) %>%
+   summarise(sumPct = sum(percentage))
# A tibble: 2 × 2
  condition sumPct
  <chr>      <dbl>
1 CDD            3
2 STD            3

```

Note: as seen above, I'm using `geom_col`instead of `geom_bar` with a different statistical summary method. If you're using *stat=%something%*, there's probably a better way to approach the task. The following is from the first paragraph of the R help page for `geom_bar`/`geom_col`:

> There are two types of bar charts: geom\_bar() and geom\_col(). geom\_bar() makes the height of the bar proportional to the number of cases in each group (or if the weight aesthetic is supplied, the sum of the weights). If you want the heights of the bars to represent values in the data, use geom\_col() instead. geom\_bar() uses stat\_count() by default: it counts the number of cases at each x position. geom\_col() uses stat\_identity(): it leaves the data as is.

---
Tags: r, ggplot2
---


# Question
Title: Running a Python Script from a bash script on multiple input files across different directory paths
I am in the process of writing up a bash script that runs three bioinformatics tools and then executes a python script. This python script compiles some of the data from the informatics tools output into a single merged output file. Now I need to run this on a large number of samples at a time so I am trying to set up the script to read through the directory path to where all the output from the tools is stored.

These are the directory paths for the python script input files:

```
EFS/Analyses/BGC.test/Sample_Name/antismash2/sample.LargeContigs.gbk
EFS/Analyses/BGC.test/Sample_Name/Bagel4/00.OverviewGeneTables.json
EFS/Analyses/BGC.test/Sample_Name/Deepbgc/Deepbgc.bgc.tsv

```

Now there are multiple Sample\_Name directories containing identical structures in the EFS/Analyses/BGC.test path. I would like to run the script on all of them

To do this, I have written a `for` loop within the bash script as follows:

```
sample_paths=$(find $EFS/Analyses/BGC.test -mindepth 1 -maxdepth 1 -type d \ )
while read sample_paths;
do
   sample=$(basename $sample_paths)
   GenBank_Input="$sample/antismash2/$sample.LargeContigs.gbk";
   TsvFile_Input="$sample/Deepbgc/Deepbgc.bgc.tsv";
   Json_Input="$sample/Bagel4/00.OverviewGeneTables.json";
done <<< “$sample_paths”

export GenBank_Input
export TsvFile_Input
export Json_Input

```

The goal here being to define the sample path for every sample and then for every sample, looping through to the specific directories of each input file. I then exported the variables and used os.environ from the os package to use them in the python script as such:

```
    #parse gbk input file
    for record in SeqIO.parse(os.environ["GenBank_Input"], 'genbank'):
    #parse tsv input file
    file_handle = open(os.environ["TsvFile_Input"])
    #parse json input file
    json_data = json.load(open(os.environ["Json_Input"]))

```

Now the error I get when I run the script is a `FileNotFoundError: [Errno 2] No such file or directory: '“”/antismash2/“”.LargeContigs.gbk'` so the error seems pretty clear cut. The issue seems to be with the path that I defined within the while loop variables but as of right now, I have been unable to resolve. I am trying to familiarize myself more with bash scripting as I am still fairly new to it. Any guidance on how I might adjust the code to resolve this and get the script working is much appreciated, thanks!

# Answer
> 2 votes
Your loop is more complicated than you need and it won't work properly anyway since you're using the same variable (`sample_paths`) to iterate and save. Try this instead:

```
for dir in EFS/Analyses/BGC.test/*/; do ... done

```

Next, all you really need is the sample names. The rest can be inferred:

```
sample_names=""
for dir in EFS/Analyses/BGC.test/*/; do
    sample_names+="$(basename "$dir") "
done
## remove the extra space added in the end
sample_names=$(sed 's/  *$//' <<< "$sample_names")
export sample_paths

```

Then, in python, you would do:

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

import os

base_dir="EFS/Analyses/BGC.test"
for sample_name in os.environ["sample_names"].split(' '):
    gbk_file = """%s/%s/antismash2/%s.LargeContigs.gbk""" % (
        base_dir, sample_name, sample_name
    )
    tsv_input = """%s/%s/Deepbgc/Deepbgc.bgc.tsv""" % (base_dir, sample_name)
    json_input = """%s/%s//Bagel4/00.OverviewGeneTables.json""" % (base_dir, sample_name)
    print("============== Sample Name:",sample_name, "=============")
    print("GenBank file:", gbk_file)
    print("Tsv file:", tsv_input)
    print("Tsv file:", json_input)

```

I created the following test directories:

```
$ tree
.
└── EFS
    └── Analyses
        └── BGC.test
            ├── Sample_1
            │   ├── antismash2
            │   │   └── sample.LargeContigs.gbk
            │   ├── Bagel4
            │   │   └── 00.OverviewGeneTables.json
            │   └── Deepbgc
            │       └── Deepbgc.bgc.tsv
            ├── Sample_2
            │   ├── antismash2
            │   │   └── sample.LargeContigs.gbk
            │   ├── Bagel4
            │   │   └── 00.OverviewGeneTables.json
            │   └── Deepbgc
            │       └── Deepbgc.bgc.tsv
            └── Sample_3
                ├── antismash2
                │   └── sample.LargeContigs.gbk
                ├── Bagel4
                │   └── 00.OverviewGeneTables.json
                └── Deepbgc
                    └── Deepbgc.bgc.tsv

```

And running the two snippets of code above produces:

```
$ foo.py
============== Sample Name: Sample_1 =============
GenBank file: EFS/Analyses/BGC.test/Sample_1/antismash2/Sample_1.LargeContigs.gbk
Tsv file: EFS/Analyses/BGC.test/Sample_1/Deepbgc/Deepbgc.bgc.tsv
Tsv file: EFS/Analyses/BGC.test/Sample_1//Bagel4/00.OverviewGeneTables.json
============== Sample Name: Sample_2 =============
GenBank file: EFS/Analyses/BGC.test/Sample_2/antismash2/Sample_2.LargeContigs.gbk
Tsv file: EFS/Analyses/BGC.test/Sample_2/Deepbgc/Deepbgc.bgc.tsv
Tsv file: EFS/Analyses/BGC.test/Sample_2//Bagel4/00.OverviewGeneTables.json
============== Sample Name: Sample_3 =============
GenBank file: EFS/Analyses/BGC.test/Sample_3/antismash2/Sample_3.LargeContigs.gbk
Tsv file: EFS/Analyses/BGC.test/Sample_3/Deepbgc/Deepbgc.bgc.tsv
Tsv file: EFS/Analyses/BGC.test/Sample_3//Bagel4/00.OverviewGeneTables.json

```

All this said, I don't really understand why you are using the shell at all. You can implement all of this directly in python:

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

import os

base_dir="EFS/Analyses/BGC.test"
sample_names = os.listdir('EFS/Analyses/BGC.test/')

for sample_name in sample_names:
    gbk_file = """%s/%s/antismash2/%s.LargeContigs.gbk""" % (
        base_dir, sample_name, sample_name
    )
    tsv_input = """%s/%s/Deepbgc/Deepbgc.bgc.tsv""" % (base_dir, sample_name)
    json_input = """%s/%s//Bagel4/00.OverviewGeneTables.json""" % (base_dir, sample_name)

    # do something with these files here

```

---
Tags: python, bash
---


title: Summing up data effectively in excel

text:
# Question
Title: Summing up data effectively in excel
I have some tsRNA count data in excel which I need to add selectively in blocks based on the count

```
refID   Count
1) tRNA-Ala-AGC-1   1
tRNA-Ala-AGC-1  2
tRNA-Ala-AGC-1  2
tRNA-Ala-AGC-1  1
tRNA-Ala-AGC-1  3
tRNA-Ala-AGC-1  300
tRNA-Ala-AGC-1  51
tRNA-Ala-AGC-1  86
tRNA-Ala-AGC-1  6
tRNA-Ala-AGC-1  2
tRNA-Ala-AGC-1  1
tRNA-Ala-AGC-1  5

2)tRNA-Ala-AGC-11 1
tRNA-Ala-AGC-11 2
tRNA-Ala-AGC-11 1
tRNA-Ala-AGC-11 64
tRNA-Ala-AGC-11 1
tRNA-Ala-AGC-11 1
tRNA-Ala-AGC-11 2
tRNA-Ala-AGC-11 3
tRNA-Ala-AGC-11 1
tRNA-Ala-AGC-11 1
tRNA-Ala-AGC-11 1

3) tRNA-Ala-AGC-3
tRNA-Ala-AGC-3
tRNA-Ala-AGC-3
tRNA-Ala-AGC-3
tRNA-Ala-AGC-3
tRNA-Ala-AGC-3
tRNA-Ala-AGC-3
tRNA-Ala-AGC-3
tRNA-Ala-AGC-3
tRNA-Ala-AGC-3
tRNA-Ala-AGC-3
tRNA-Ala-AGC-3
tRNA-Ala-AGC-3

```

I am using sum and subtotal functions in excel to get the total counts of each of the blocks. I have about 22,000 values and several hundred groups like this. Is there a way to sum up all such blocks all at once based on the groups or at least a shorter way to sum things work other than autosum and subtotal?

Any help will be useful.

# Answer
Suppose the header is on row 1 and starts at column 'A'. Then, for the above example, just point to a cell at say `C1` to `C3` and write ...

```
=sum(B2:B12)
=sum(B14:B25)
=sum(B27:B38)

```

This gives the subtotals.

For the total sum point the mouse at `D1` and write:

```
=sum(B2:B38)

```

... except of course the real value is ...

```
=sum(B2:B22000)

```

In all cases you must put the `=` to start the cell.

Its better to try out `R`, or Python `pandas` or just a Python list `print (sum([1,2,3,4,5])` but Excel's okay.

---

The issue is its slow, its taking hours to calculate the sum of 22 000 cells - thats not a lot even for Excel. In Python that wouldn't even exceed 1 second, literally. What you are facing is **a RAM bottleneck** \- *there is literally no other explanation*.

Two suggestions firstly, look at Windows Excel, particularly if you are using OSX. This is because historically Windows was much more RAM efficient than OSX. Now however, I don't know. Secondly, simply switch to a desktop with 4 Gig RAM or more. Literally thats all it will take.

The drop back position is to do the calculation in two halves on the low RAM machine, i.e. split the data set in two. Not a cool solution but it will probably work. The final solution is to look at the stuff below the line.

---

The solutions are to shift to `VB` which none of us code in these days, but its the code underneath Excel, or the alternative is to write a custom script in a different language that will do the adding such as `Python`, `R` ... even `bash`.
> 1 votes

# Answer
It looks like what you want to do is a pivot table, which provides summary statistics for fields, grouped by specified categories. As a prerequisite, this requires an input dataset that is well formatted; more generally:

* Put each observation or sample in its own row; don't leave row gaps between samples or groups. Leaving gaps makes bulk analysis more difficult and time consuming. It can still be done, it just means that the data needs to be pre-cleaned, or alternatively analysed.
* Put all your variables in columns - the thing that vary between samples, like ‘strain’ or ‘DNA-concentration’
* Have column names be explanatory, but without spaces
* Do not combine multiple pieces of information in one cell

For more information, see this Data Carpentries Workshop on the recommended way to structure data in spreadsheets:

https://datacarpentry.org/organization-genomics/01-tidiness/index.html#structuring-data-in-spreadsheets

If the data is prepared in this way, then selecting the whole table and creating a pivot table should work. Unfortunately, the programs keep changing where in the menu that can be found. Currently in LibreOffice, it's `Data` -\> `Pivot Table` -\> `Insert or Edit...`, wheareas on my version of Excel on Mac OSX, it's found at `Data` -\> `Summarise With Pivot Table`.

This will load up a visual interface that allows you to drag and drop fields to summarise. For this data, you'd probably want *RefID* in rows, and *sum(Count)* as the data field:

After setting this up and clicking `Ok`, you should be greeted with a smaller table containing summarised variables:
> 0 votes

---
Tags: sequence-analysis
---


# Question
Title: Calculating average coverage for .bam files (sequence data)
*(Full discolosure that this is my first time working with sequence data, and with the bash scripting.)*

I need to calculate the average coverage for any .bam file.

After some searching I wrote the following script:

```
# Script to calculate the average coverage of a genome sample

# compute the total length of the sample
tot=$(samtools view -H $1 | grep -P '^@SQ' | cut -f 3 -d ':' | awk '{sum+=$1} END {print sum}')
echo $tot

#compute the coverage at each point
samtools depth -ao cov_$1 $1 #your .bam file has depth calculated for each point

#calculate the average coverage
sum=$(awk '{sum+=$3} END {print sum}' cov_$1)
echo $sum
echo

avg=$(echo "$sum/$tot" | bc -l)
echo ''The average coverage of the sample $1 is $avg x.''

```

The logic behind the calculation of tot is as follows:

`samtools view -H`  gets the header of the bam file

`grep -P` gets all lines in this header with SQ

`cut -f 3 -d ':'` removes fields that are lists of more than 3 letters while -d uses delimiter ':' instead of tab and this leaves only the numbers i.e. lengths of the chromosomes

`awk '{sum+=$1} END {tot=sum}'` calculates the length/total bases in the chromosome

However, my output seems incorrect. Any advice on where I might be making an error (in code or thought process), or a better method would be appreciated.

*One thing that worries me is that my tot variable is the total number of bases, so if the samtools depth function does not calculate coverage for each base then my sum and tot will not be in the same units and my calculated coverage will be wrong...*

# Answer
> 9 votes
For a quick estimate you’re making it more complicated than necessary.

1. The theoretical average coverage is $\frac{n \cdot \hat l}{N}$ where $n$ is the number of reads, $\hat l$ is the average read length and $N$ is the genome size.
2. `samtools idxstats` gives you the chromosome lengths and number of mapped reads in one convenient list.

Putting this together, we get

```
samtools idxstats data.bam \
| awk -vreadlen=120 '
    {
        len += $2
        nreads += $3
    }
    END {
        print nreads * readlen / len
    }
'

```

(In this case I’ve used 120 as the read length.)

However, note that this over-estimates *empirical* sequencing depth, mainly due to the fact that a given read may not be fully mapped; thus the the average read length is an upper bound on the average *mapped* length. If you have an estimate for the average percentage of mapped bases per read, you can multiple that by the above coverage estimate for a more accurate number.

---

If you need a more accurate estimate, your approach is fundamentally sound (and the result is probably correct!). For low-coverage data you can speed up your analysis by dropping the `-a` argument to `samtools depth`: you don’t need to record zero-coverage bases, they don’t contribute to the total. You also don’t need to store the coverage into an intermediate file, thus reducing IO cost:

```
sum=$(samtools depth "$1" | awk '{sum+=$3} END {print sum}')

```

Some more remarks on your code:

* The `-P` flag for `grep` is not supported everywhere and thus limits portability, so you want to avoid using it if possible. It turns out that you don’t need it here: `grep '^@SQ'` works just as well.
* While the philosophy of UNIX is based on solving problems by combining individual tools that each solve one part, it is sometimes convenient to merge individual steps:
    + You can often combine `cut` and `awk`, in this case by setting the field separator:
    ```
    awk -vFS=: '{sum+=$3} END {print sum}'

    ```
    + You can often combine `grep` and `awk`:
    ```
    awk -vFS=: '/^@SQ/ {sum+=$3} END {print sum}'

    ```
* *Always* double-quote your shell variables, otherwise your script will break if the values contain unexpected spaces or other special characters. This is especially important for filenames, because spaces and other special characters *are valid characters in paths!*
* Your code is clean and readable. Well done.

# Answer
> 2 votes
With the current versions of samtools, you can also use:

```
samtools coverage sample.bam

```

The documents are here,

> samtools coverage – produces a histogram or table of coverage per chromosome

A further example from the site,

```
samtools coverage -r chr1:1M-12M input.bam

```

**Output**

```
#rname  startpos  endpos    numreads  covbases  coverage  meandepth  meanbaseq  meanmapq
chr1    1000000   12000000  528695    1069995   9.72723   3.50281    34.4       55.8

```

The histogram option is

```
samtools coverage -A -w 32 -r chr1:1M-12M input.bam

```

... and output available on the website link.

---
Tags: sequence-analysis, bash, coverage
---


title: Compare VCF files

text:
# Question
Title: Compare VCF files
I have two strains, **X** and **Y**, and I would like to compare the variants found in both. I have my own parent strain, **Z**, that was sequenced and assembled to be used as the reference genome.

I first aligned **X** and **Y** to **Z** using minimap2 and then I used FreeBayes as the variant caller. Once I had the two separate VCF files I looked online to see what was the best way to filter out the unique and common variants in **X** and **Y**.

I found the `bcftools isec` tool to be my best bet so I went ahead and ran this line of code:

`./bcftools isec -p out/path/ X.vcf.gz Y.vcf.gz`

Now I have three new VCF files:

* `0000.vcf` (variants unique to X)
* `0001.vcf` (variants unique to Y)
* `0002.vcf` (variants in both X and Y)

However, it's still in VCF format. I need a way to parse through the VCF file to list out the SNPs with high quality with their position on the genome. I tried using `bcftools stats` but it just lists out how many SNPs there are and not their position.

**TLDR:** Can someone help me find a way to list out the locations of SNPs from a VCF file in human readable format?

# Answer
> 1 votes
If I understood you correctly you are searching for a way to represent the genomic locations in your vcf file. However, the vcf file is already a human readable file (vcf specification).

The first lines are the header, but afterwards you can find your information in a tab delimited fashion:

```
#CHROM POS ID REF ALT QUAL FILTER INFO

```

So the information on genomic positions is directly readable from the vcf file (first 2 columns)

---
Tags: vcf, snp, variant-calling, bcftools, variants
---


title: Regarding chi square?

text:
# Question
Title: Regarding chi square?
I need some suggestions regarding an analysis. I might be wrong as I have no stats background.

I have list of exon with differential usage across samples. I can find lots of denovo mutations from both proband and siblings in those exons. Now I am trying to check was there any enrichment of proband specific denovo mutations in those exons. I can think of doing a chi square test. But I am not sure what will be expected value in chi square test. Imagine I am having 700 denovo proband mutation in 650 samples out of 1500 probands and 300 denovo mutations in 320 unaffected siblings, out of 1200. Any suggestions? Thanks

|  | Denovo | WT | Total |
| --- | --- | --- | --- |
| Proband | 553 | 1202 | 1755 |
| Unaffected sibs | 505 | 813 | 1318 |
| **Total** | 1058 | 2015 | 3073 |

# Answer
> 1 votes
I'm a bioinformatician, but not a statistician, and it's been a while since I did any $\chi^2$ tests, so take this answer with a grain of salt.

From what I remember, the most important thing is to be clear about what your observed and expected values are, and to be consistent. In other words, ***it doesn't matter what you pick, as long as you make it clear what you've picked***. For your situation there are at least four possible "expected" values that can be used:

1. Average counts based on your available data (this is probably the most commonly used approach)
2. Equal proportions based on your available data (this is less common, because it alters table totals)
3. Unaffected sib counts (also less common, because it alters totals)
4. Population average counts based on public data

Unfortunately, there's no answer that will work for all situations. It's a good idea to think about which of these (or any others) makes sense in your particular case, what you're actually trying to compare, and how much you want to generalise to other populations. It's almost always the case that you know your own data better than anybody else, so are in the best position to make this decision. For example, if this is a very rare variant in the general population, you might choose to look for general population data for the variant, but if you're only interested in the impact of people who are aware of a family history, then unaffected sibs might be a better approach.

In any case, the calculations for generating the expected values are similar: work out the proportion of the variant in the expected population, and apply that proportion to the total count in the observed population.

In the first case, the expected counts would be something like this \[note that the table totals are preserved\]:

|  | Denovo | WT | Total |
| --- | --- | --- | --- |
| Proband | \*604.2 | 1150.8 | 1755 |
| Unaffected sibs | 453.8 | 864.2 | 1318 |
| **Total** | 1058 | 2015 | 3073 |

* calculation: row total (1755) as a proportion of full total (3073) = 0.5711, which is then multiplied by the column total (1058) = 604.2

In the second case, it would be this \[note that the row totals are different\]:

|  | Denovo | WT | Total |
| --- | --- | --- | --- |
| Proband | 529 | 1007.5 | 1536.5 |
| Unaffected sibs | 529 | 1007.5 | 1536.5 |
| **Total** | 1058 | 2015 | 3073 |

In the third case, it's simply comparing Proband vs Unaffected counts, normalising based on the row total:

|  | Denovo | WT | Total |
| --- | --- | --- | --- |
| Proband | 672.4 | 1082.6 | 1755 |
| Unaffected sibs | 505 | 813 | 1318 |
| **Total** | 1177.4 | 1895.6 | 3073 |

Note that the second row can be omitted from $\chi^2$ calculations here, because it's *identical* to the expected values, so adding that difference will not change the calculated $\chi^2$ sum.

\[this could also be done based on the column total, if you wanted to compare Denovo vs WT, instead of Proband vs Unaffected\]

# Answer
> 1 votes
It seems reasonable, it depends how you do the Chi2.

The Chi2 distribution depends on a good estimate of the expected (E). Here your observed (O) is the "affecteds" and you want to assess whether this has a suspected elevated mutation rate, then that looks a meaningful measure. There's two approaches

1. Your E is the unaffected siblings and E is calculated by averaging this value. You could also perform the reverse calculation where E was the average for the affecteds and O is the unaffected. There will be an issue about the degrees of freedom, which would require some thought.
2. There is an alternative way to perform Chi2, which I forget because O and E are matched (both siblings), here there is no averaging. You'd need to look it up, the numerator is the summed difference for each sibling group ... but I can't remember.

Point 2 would be a good test, but you'd need to check the test.

A phylogeneticist would also use the Chi2 in a very different test. They would draw a tree and designate two rates, one for affected and unaffected. If there was elevated mutation rate in the affected it could be tested on a likelihood ratio test, i.e. the rate of evolution was increased. This can be done because the test tree would be one where the rate was the same between affecteds and unaffected progeny. The way the test works is the likelihood difference between the two trees is assessed and because the number of increased parameters (degrees of freedom) is known and log probabilities can be subtracted from each other and directly assessed (I think its degrees of freedom + 1 ... 1 additional parameter so its 2) - and therefore the difference in likelihood is tested directly against the Chi2. Its kinda neat and exploits the behaviour of log probabilities.

---
Tags: statistics
---


# Question
Title: How to get the external (cluster) job id and not internal snakemake job id in cluster log files?
**Background** I am using `snakemake` and I want to add the external (cluster) jobid to my cluster output and error log files. Below is my shell script.

`snakemake --configfiles config/config_resources.yml config/config_parameters.yml --latency-wait 60 --use-conda --printshellcmds --jobname PipeName_{name}.jobid{jobid} --jobs 300 --cluster "bsub -W 480 -n {threads} -M {resources.mem_gb}G -R 'rusage[mem={resources.mem_gb}G]' -o logs/{name}_{jobid}.out -e logs/{name}_{jobid}.err"`.

But I only get the internal one in de log name from snakemake (e.g. 1 to 5), instead of the one from the cluster.

**Question** How to I get the external job id in there and not the local one made by snakemake? Or where in my shell script did I go wrong?

Cheers and thanks in advance.

# Answer
bsub supports some filename placeholders using a percent sign, like `%J` for the LSF job ID, so you can use that:

`snakemake ... --cluster "bsub ... -o logs/{name}_%J.out -e logs/{name}_%J.err"`

That way snakemake will format the cluster command string with its wildcards (like `name`) as it constructs the bsub command, and then bsub will put the job ID in place in those filenames. Since Snakemake just sees that as one big command to run for the cluster and doesn't know the details other than formatting the text, I think the two shouldn't clash.
> 1 votes

---
Tags: snakemake, hpc
---