text:
Can't install snakemake v5.15.0 from conda
I'm trying to install the latest version of snakemake from conda into a conda env, following the instructions here. I.e.
conda create -n snakemake -c conda-forge -c bioconda snakemake
However, this seems to only give me the option to install snakemake v5.3.0, while the latest no-arch version on conda is v5.15.0.
Any help greatly appreciated.
Solving environment: done

## Package Plan ##

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

  added / updated specs:
    - snakemake

The following packages will be downloaded:

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

The following NEW packages will be INSTALLED:

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

Proceed ([y]/n)?

Does it help if you include the specific version? conda create -n snakemake -c conda-forge -c bioconda snakemake=5.15

@mrhd - yes, it does, thanks for the suggestion

This is an example of the conda solver behaving incorrectly. In this case, it's choosing to pick a version that minimizes the number of dependencies required. Version 5.3.0 has 25 fewer dependencies, so if you let conda choose the version that's what it will pick. There's nothing you can really do about this other than specifying the version of snakemake that you want.
Update: It's very likely that mamba will correctly do what you want by default. Though I haven't used it, it's supposed to be a drop in replacement for conda here.

That works! Thanks Devon (I'm sure that's one of the things I tried, but perhaps when my channels were still in the wrong order...)

From next release on, installing via mamba is the default recommendation for Snakemake.

Actually, the issue does not seem to be the number of dependencies, but rather a lack of information in the stripped down current_repodata.json. If that is disabled and conda uses the full repodata, results are correct again, see https://github.com/conda/conda/issues/9905

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

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

anaconda is just miniconda with a bunch of additional packages installed.


text:
How to run a one-time 'setup' step for amrfinder in snakemake using conda
I'm trying to add amrfinder running within conda to a snakemake workflow, working on a slurm HPC.
The first time you run amrfinder, you have to run amrfinder -u to update/download the database of genes/proteins which it searches for. Since I'm using the --use-conda flag of snakemake, the first time I run the workflow, it creates a new env for the conda package, so this command has to be run before that copy of amrfinder will work.
My question is: how do I run amrfinder -u within snakemake? If I do it as part of the standard amrfinder command, it runs as many times as there are samples, and as all the jobs are looking at the same database, when it updates for one job, it breaks the other jobs.
I can't figure out how to run this command in a 'snakemake-y' way, as it has no input, it's just literally amrfinder -u.
Of course, I can make it work with a hack/workaround, but I'm interested in a snakemake way of doing it.
This is the main snakemake rule for running amr_finder_plus, the new rule, to run amrfinder -u needs to run before this
rule amr_finder_plus:
    input:
        assembly = 'my_assembly.fasta'
    output:
        amr_finder_plus_results = 'my_amr_finder_plus_results.tsv'
    conda:
        '../../envs/amrfinderplus.yaml'
    shell:
        'amrfinder -n {input.assembly} -O Salmonella --output {output.amr_finder_plus_results} --threads 4'

There are two approaches you can take:

There's a --conda-create-envs-only option that you can use to exit after creating the envs. You should get information about where the envs are so you can run amrfinder -u.
Add a new rule called "amr_setup" or something like that that runs:

.
amrfinder -u
touch amrfinder.setup

Then put amrfinder.setup into your snakemake workflow as an input file into the amr_finder_plus rule. Either of these will work

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
How do I tell that a MinION run has finished based on the files produced?
I want to automate some scripts which should run after a MinION run completes. Is there a file which is produced by MinKNOW which will enable me to tell that the run has finished?

Use rampart and just stop it once you've happy with coverage

You can just look for the presence of a file called "final_summary_*.txt" in the run folder. From Matt Loose.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
DNA methylation mediator analysis for environmental exposure and health outcome
We measured blood DNA methylation using epic array in 100 no-diabetic and 100 type 1 diatbetic patients matched for age, sex and smoking and also measured arsenic exposure using ICP-MS. Could you please suggest best way to find CpG sites associated with arsenic exposure and diabetes. We know that toxic element leads to health risk but we want to study this by using DNA methylation as mediator for disease risk.
Could you please suggest any R package for doing these types of analysis. I have used previously limma for toxic element exposure and DNA methylation association analysis but not sure how make connections with health risk phenotype?


text:
How to convert the .vcf (imputed) file with GT:GP format to GT:DS?
I have the genotyped data from impute2 output in .gen format (imputed to 1000G P3). The file has genotype posterior probabilities (GP:3 values per variant). I have converted .gen to .vcf using qctools and the .vcf file has GT:GP format. I need to convert the .vcf file with GT:GP format to GT:DS. Genotype dosages are recommended for use in qtltools/fastqtl analysis. However, I cannot find a tool that would keep the .vcf format and convert GP to DS. Any help much appreciated!

You can do this in Hail.
Here's the rough code to do it (0.1 versions).
Setup:
from hail import *
hc = HailContext()

Import the .gen file. VCF works too:
dataset = hc.import_gen(
    'src/test/resources/example.gen',
    'src/test/resources/example.sample')

Remap the genotype schema and export to VCF:
dataset.annotate_genotypes_expr('g = {GT: g.call(), DS: g.dosage()}')\
    .export_vcf('/tmp/out.vcf.bgz')

Take a look at the getting started page if you want to try it out!
I should note that you may be able to do QTL analyses in Hail, depending on the method you want to use. See blog post here.

Hm, I didn't know that the plugin existed so I wrote my own script to convert GP to minor allele dosage on github. Maybe someone else will find it useful :) https://github.com/7methylg/VCF-GP-to-DS

There's the dosage plugin for bcftools, but it only outputs tab separated values. It would not be too hard to extend the plugin to output a VCF with the DS tag instead, but it has not been done yet.  There's a good chance the bcftools devs would respond to a feature request...
In any case, this code:
curl https://raw.githubusercontent.com/samtools/bcftools/develop/test/convert.vcf > convert.vcf
bcftools +dosage convert.vcf > output.tsv
head -2 output.tsv

has the output:
#[1]CHROM       [2]POS  [3]REF  [4]ALT  [5]NA00001      [6]NA00002      [7]NA00003      [8]NA00004      [9]NA00005      [10]NA00006     [11]NA00007     [12]NA00008    [13]NA00009     [14]NA00010
X       2698560 G       A       0.1     0.0     0.1     0.2     0.3     0.2     0.2     0.2     0.2     0.1

This is using bcftools version 1.3.1.
Here is an excerpt from the bcftools manual for the dosage plugin:

dosage
print genotype dosage. By default the plugin searches for PL, GL and GT, in that order.

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/


text:
Workflow of metabarcoding analysis
For a diet analysis of an insect-eating animal, all species in a sample shall be identified. For this a sequencing of the metagenomic sample was done, where the COI/COX region was used as a barcode and amplified. This resulted in 20 million reads (after quality control).
While looking for an apropriate workflow, I am faced with these issues:

For many tools 20 million reads are to much to handle (BLAST for
instance)
For the majority of metagenomic tools the focus is on micro-organisms, not insects
The database to search against must contain COI/COX data, mainly from insects, however, the animal could also eat snails or frog eggs...
Since the reads are amplicons some analysis are not recommended

Leaving me with these questions:

What are the typical analysis steps in a workflow like this?
Which tools could I use for my purpose?
Where could I get a database of COI/COX data?

I'm working on something similar as a side project. I'll detail the steps that I'm using, but I'm curious if others can suggest a better way.

You don't mention if your data are paired-end, but in my case I have paired-end sequences that overlap in the middle. I used FLASh to merge the reads and seqtk to convert the reads to fasta format.
I used the cluster_fast algorithm in vsearch to cluster sequences at 100% identity, then filtered out rare sequence types (<10 sequences) to account for sequencing errors.
This resulted in a manageable number of unique sequences, which I searched against a custom database using vsearch, but which could also have been searched against the nr database using blast.

Regarding your third question ("Where could I get a database of COI/COX data?"), you can apparently search and download fasta sequences from BOLD.
However, I don't see a "select all" button to select all records found by one search in a single fasta file.
There seems to exist an R API: https://cran.r-project.org/web/packages/bold/bold.pdf which refers to this: http://www.boldsystems.org/index.php/resources/api#sequenceParameters>.
From this I seem to understand that one can directly query using urls:
wget -O Gastropoda.fa http://www.boldsystems.org/index.php/API_Public/sequence?taxon=Gastropoda&marker=COI

And I get a fasta file containing COI sequences for Gastropoda. There may be ways to refine the query geographically, using the geo parameter.

Yes, I already downloaded some of the BOLD data. I circumvented the missing "select all" button by using a firefox addon for selection boxes. I did not know that there is an R API, that will be helpful.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/


text:
Can Mauve export indels or is there an whole-genome aligner that can export indels?
I have aligned two whole genomes of the same species in Mauve using progressiveMauve. One of the genomes is 2 megabases longer than the other so I know it must contain some insertions. Looking at the alignment in the Mauve GUI, it seems that insertions are small and scattered throughout the genome. I have exported the SNPs from the GUI. Is there a similar option for indels?
I have also inspected the Mauve file written after the alignment. It contains the two original sequences in FASTA format but I noticed in contains gaps denoted by a dash (-). Would it be possible to create a VCF from that file?


text:
hgsql not found when running a script despite being accessible on command-line
I am working on macOS. I am following the steps described here to lift over an annotation from the reference genome of C. elegans to another assembly of C. elegans.
I have run this code:
time (doSameSpeciesLiftOver.pl -verbose=2 -buildDir=`pwd` \
      -ooc=`pwd`/${target}.ooc -fileServer=localhost -localTmp="/dev/shm" \
  -bigClusterHub=localhost -dbHost=localhost -workhorse=localhost \
  -target2Bit=`pwd`/${target}.2bit -targetSizes=`pwd`/${target}.chrom.sizes \
  -query2Bit=../${query}/${query}.2bit \
  -querySizes=../${query}/${query}.chrom.sizes ${target} ${query}) > do.log 2>&1

However, the do.log file contains:
zsh:1: command not found: hgsql
HgStepManager: executing from step 'align' through step 'cleanup'.
HgStepManager: executing step 'align' Sun Nov 11 15:38:41 2018.
Using localhost, /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/GCA_000004515.3_Glycine_max_v2.0.2bit and ../GCF_000004515.3_V1.1/GCF_000004515.3_V1.1.2bit
# chmod a+x /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/run.blat/job.csh
# chmod a+x /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/run.blat/doAlign.csh
# ssh -x -o 'StrictHostKeyChecking = no' -o 'BatchMode = yes' localhost nice /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/run.blat/doAlign.csh
cd /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/run.blat
rm -rf tParts
/Users/username/Documents/projects/main/output/software/liftover/data/scripts/partitionSequence.pl 58018742 0 /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/GCA_000004515.3_Glycine_max_v2.0.2bit /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/GCA_000004515.3_Glycine_max_v2.0.chrom.sizes 2000 -lstDir=tParts
rm -rf qParts
/Users/username/Documents/projects/main/output/software/liftover/data/scripts/partitionSequence.pl 10000000 0 ../GCF_000004515.3_V1.1/GCF_000004515.3_V1.1.2bit ../GCF_000004515.3_V1.1/GCF_000004515.3_V1.1.chrom.sizes 1000 -lstDir=qParts
Can't open seq size file ../GCF_000004515.3_V1.1/GCF_000004515.3_V1.1.chrom.sizes: No such file or directory
Command failed:
ssh -x -o 'StrictHostKeyChecking = no' -o 'BatchMode = yes' localhost nice /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/run.blat/doAlign.csh

I think the biggest problem is that the hgsql command is not found.
If I run hgsql on the command line, I find it even though I get the following error:
dyld: Library not loaded: /usr/local/opt/mysql/lib/libmysqlclient.20.dylib
  Referenced from: /usr/local/bin/hgsql
  Reason: image not found
[1]    34059 abort      hgsql

Location of the hgsql installation:
which hgsql
/Users/username/Documents/projects/main/output/software/liftover/data/bin/hgsql

Can anyone help me with using hgsql, please?

It is possible that the shell being used in the process is different from your standard command line shell. I am not familiar with the tool, and it's hard to say more without knowing your OSX version, but it looks like it is using zsh. I think that zsh only fairly recently became the standard OSX shell. I don't know the tool or the shell very well, but it looks like you might be able to explicitly set up hgsql in your PATH in your .zshrc.
It's also possible that the liftover script was not built to run via zsh. I notice that it is trying to run csh scripts as part of the workflow. Then again the localhost may insulate it.
So I think that my first step would be to sort out which shell I am using, whether it works with the liftover pipeline, and how it's running on the localhost. I notice that the ssh call running the csh script on the localhost is the step that actually fails. What shell is executing the csh script? Does that shell know about hgsql? It could be that the way that the error bubbles up it is hiding which shell it is actually using.
I'd recommend trying that ssh call in isolation from the rest of the pipeline, and diagnosing the error that way.
Hope that helps.


text:
BUSCO does not find buscos in reference sequence of C. elegans (technical problem)
I am working on a GNU/Linux machine. I have installed BUSCO using conda. I have downloaded the nematode dataset nematoda_odb9.tar.gzhere. I downloaded the reference sequence of C. elegans from Ensembl. I then ran busco using the following command:
run_busco --in Caenorhabditis_elegans.WBcel235.dna.toplevel.fa --cpu 1 --out busco --lineage_path ./nematoda_odb9/ --mode genome --force > busco.o 2> busco.e
tblastn seems to have run properly. run_busco/blast_output/tblastn_busco.tsv contains lots of hits for many of the buscos.
However, I have 982 "Missing BUSCOs".
The file that contains logs and errors related to augustus (run_busco/augustus_output/augustus.log) contains:
Will create parameters for a EUKARYOTIC species!
creating directory ...anaconda/envs/busco3_0_2/config/species/BUSCO_busco_1395218643/ ...
creating ...anaconda/envs/busco3_0_2/config/species/BUSCO_busco_1395218643/BUSCO_busco_1395218643_parameters.cfg ...
creating ...anaconda/envs/busco3_0_2/config/species/BUSCO_busco_1395218643/BUSCO_busco_1395218643_weightmatrix.txt ...
creating ...anaconda/envs/busco3_0_2/config/species/BUSCO_busco_1395218643/BUSCO_busco_1395218643_metapars.cfg ...
The necessary files for training BUSCO_busco_1395218643 have been created.
Now, either run etraining or optimize_parameters.pl with --species=BUSCO_busco_1395218643.
etraining quickly estimates the parameters from a file with training genes.
optimize_augustus.pl alternates running etraining and augustus to find optimal metaparameters.

It does not seem to contain any errors.
I have the following error:
Failed to open sequence file run_busco/augustus_output/extracted_proteins/EOG091H07ZC.faa.1 In fact, the folder run_busco/augustus_output/extracted_proteins/ is empty.
Anyone would know why I have no extracted proteins?

Most likely your Augustus installation is not functioning properly. Did you install Augustus from source or using conda? I believe the conda version of Augustus is known to have some issues.
Something worth trying: Running a complete busco install in a VM. This image has all dependencies installed and properly configured.

It was installed via conda, the log output clearly shows the Anaconda being accessed for busco installation. Personally I'm surprised because Anaconda programs are very reliable, but that might just be the core installation (at least for  Mac OSX). Anyway its good advice to try a direct installation.
Could you please point out what "some issues" might be?

BUSCO worked with this Docker image: https://hub.docker.com/r/comics/busco I'll try again using the conda installation and point towards some issues but for the moment, I am happy I could run it in a Docker container.

@MichaelG. you are right that is quite reliable in general. Augustus however, is quite tedious in it's installation (e.g. the location of the config file and the location of it's species training-sets). For example, the manual of BRAKER(the automated gene annotation program) states: "Problems have been reported from users that tried to run BRAKER with AUGUSTUS releases maintained by third parties, i.e. Bioconda". They therefore recommend to install Augustus as distributed by the original developers

@charlesdarwin That's good to hear! Most likely the person that build that docker container did a correct installation of Augustus for you :)


text:
Why can I not install snakemake on my SLURM computer even though I can find it in the bioconda channel?
I searched for the package 'snakemake' on my SLURM cluster using:
conda search --channel bioconda snakemake
and I get many versions, up to 5.4.2.
I then try to install it using:
conda install --channel bioconda --yes snakemake=5.4.2
but it fails with the following message:
PackagesNotFoundError: The following packages are not available from current channels:

  - snakemake

Current channels:

  - https://repo.anaconda.com/pkgs/main/linux-64
  - https://repo.anaconda.com/pkgs/main/noarch
  - https://repo.anaconda.com/pkgs/free/linux-64
  - https://repo.anaconda.com/pkgs/free/noarch
  - https://repo.anaconda.com/pkgs/r/linux-64
  - https://repo.anaconda.com/pkgs/r/noarch
  - https://repo.anaconda.com/pkgs/pro/linux-64
  - https://repo.anaconda.com/pkgs/pro/noarch

Why?
Why is the bioconda channel not included in the current channels? I thought that was what the --channel bioconda did in my install command.
Note that conda install --channel bioconda --yes samtools=1.9 worked fine and installed samtools successfully.

I've heard some folks warn that conda is sensitive to the order in which channels are set/configured. I'm not exactly sure what the details are, but maybe that gives you somewhere to start troubleshooting?

I created a file .condarc in my home directory with the following content:
channels:
  - conda-forge
  - bioconda
  - defaults

and it works now. I don't understand why this is necessary though. I thought the channels in the .condarc file would be used if I did not specify them in my command.

This seems to be a bug in the version of conda you're using, it works in 4.6.4. Note, however, that you needed to add -c conda-forge to your original command anyway, since many of the snakemake dependencies come from there.

good point about the dependencies. I didn't know about that. I just always give the channel of the package I want to install but it sounds like I also need to give the channel of dependencies. So I'll keep using -c conda-forge -c bioconda -c r in all my conda install commands.


text:
GMAJ error: Exception in thread "main" java.lang.ClassCastException: java.lang.StringBuffer cannot be cast to java.lang.String
I have installed GMAJ using:
wget http://globin.bx.psu.edu/dist/gmaj/gmaj.zip  -O output/software/gmaj.zip
unzip -u output/software/gmaj.zip -d output/software/gmaj/
rm output/software/gmaj.zip

I ran lastz on two chromosomes from two different assemblies of C. elegans:
lastz c_elegans.PRJNA13758.WS263.genomicI.fa I.fa --format=maf > MAF.maf
and got a MAF file.
Then, I tried to visualise the file using:
java -jar output/software/gmaj/gmaj.jar MAF.maf
but I got the following error:
Exception in thread "main" java.lang.ClassCastException: java.lang.StringBuffer cannot be cast to java.lang.String
    at edu.psu.bx.gmaj.MajGui.setDefaults(MajGui.java:272)
    at edu.psu.bx.gmaj.Maj.<init>(Maj.java:53)
    at edu.psu.bx.gmaj.MajMain.main(MajMain.java:86)

Can someone help me resolve this error or suggest an alternative way to visualise my alignment?
I work on Mac OS X.
Cheers.

It's a programming error that should look like this:
StringBuffer sb = (...)
(...)
String s = (String)sb;

while it should be
StringBuffer sb = (...)
(...)
String s = sb.toString();

( see also https://www.biostars.org/p/177992/ ) . You cannot do anything about it but asking the authors of the software.

Sent an email to Cathy Riemer asking for help.


text:
Mugsy error: Can't find species II dna at output/software/mugsy_x86-64-v1r2.3/mugsy line 501
I have installed mugsy in order to create a multiple genome alignment and a phylogenetic tree of several species of nematodes. The following command successfully pulls out the help:
output/software/mugsy_x86-64-v1r2.3/mugsy -h

To test that the software is installed correctly and that I am using it correctly, I ran the following alignment between the Wormbase sequence and the Ensembl sequence of the C. elegans genome:
output/software/mugsy_x86-64-v1r2.3/mugsy  --directory /.../project/output/software/mugsy_x86-64-v1r2.3/ --prefix mygenomes output/genome/celegans/ref/wormbase/seq/c_elegans.PRJNA13758.WS263.genomic.fa output/genome/celegans/ref/ensembl/seq/celegans.fa

I get the following output:
Parsing sequences for c_elegans  num_seqs:7
Parsing sequences for celegans  num_seqs:1
Can't find species II dna at output/software/mugsy_x86-64-v1r2.3/mugsy line 501.

The headers of the wormbase sequences are:
>I
>II
>III
>IV
>V
>X
>MtDNA

and the headers of the ensembl sequences are:
>I dna:chromosome chromosome:WBcel235:I:1:15072434:1 REF
>II dna:chromosome chromosome:WBcel235:II:1:15279421:1 REF
>III dna:chromosome chromosome:WBcel235:III:1:13783801:1 REF
>IV dna:chromosome chromosome:WBcel235:IV:1:17493829:1 REF
>MtDNA dna:chromosome chromosome:WBcel235:MtDNA:1:13794:1 REF
>V dna:chromosome chromosome:WBcel235:V:1:20924180:1 REF
>X dna:chromosome chromosome:WBcel235:X:1:17718942:1 REF

Maybe this line:
    Parsing sequences for celegans  num_seqs:1 should make me worry. Why does mugsy think there is only one sequence in that genome? Why do I get the following error and how can I fix it?:
Can't find species II dna at output/software/mugsy_x86-64-v1r2.3/mugsy line 501.

I am working on a GNU/Linux platform.


text:
Import a tab-separated file with differing numbers of elements in each row; prokka output
I am using prokka to annotate a bacterial genome:
prokka ecoli.fa
Prokka is outputting a tab-separated file (called PROKKA_12142017.tsv) with differing numbers of elements in each row:
locus_tag    ftype    gene    EC_number    product
BOFHCHGE_00001    CDS    hypothetical protein
BOFHCHGE_00002    CDS    rdgB_1    3.6.1.9    dITP/XTP pyrophosphatase
BOFHCHGE_00003    CDS    hypothetical protein
BOFHCHGE_00004    CDS    hypothetical protein
BOFHCHGE_00005    CDS    hypothetical protein

Does anybody know how I can get round this to import into a dataframe in R?
Is there a way to ask prokka to add something in the gene column?

Please make sure to send a bug report to the prokka authors, since this should be fixed on their end.

Author here - this was a bug which was fixed in August:
See the issue in Github.
You will need to use the latest Github HEAD version until the next release.

You don't really need to get around anything, your example can be loaded correctly with read.delim(). The rows with missing values are filled in with blanks. You may, however, prefer the readr package, which handles this a bit more elegantly (it'll tell you which values were missing and not do the annoying string->factor conversion).

Thanks, useful but it puts 'hypothetical protein' under the gene column when it actually belongs in the product column. The gene column should be empty for hypothetical proteins.

Oh, I thought that was intentional. In that case prokka has a bug and you'll need to test for "" as the product and just assign accordingly.

And I get 'Sec-independent protein translocase protein TatA' for EC number, so that's definitely wrong. I never see an empty element before a non-empty element.

You can use the stringsAsFactors = FALSE to prevent converting strings to factors in read.delim/read.csv...

@Llopis Quite true, but the readr package has that as the default :)

As Devon said, read.delim can deal with this perfectly well as long as the file is properly formatted (so missing fields still end with the field separator, like \t\t). For example, with this input:
Field1  Field2  Field3
Row1field1  Row1field2  Row1field3
Row2field1      Row2field3

Or, to show the fields more clearly:
Field1\tField2\tField3
Row1field1\tRow1field2\tRow1field3
Row2field1\t\tRow2field3

R can import:
> df<-read.delim("file",header=TRUE, sep="\t")
> df
      Field1     Field2     Field3
1 Row1field1 Row1field2 Row1field3
2 Row2field1            Row2field3

Alternatively, if you really want to add something to the missing columns, you can use:
$ awk -F"\t" '($3==""){$3="none" }1;' file
locus_tag   ftype   gene    EC_number   product
BOFHCHGE_00001  CDS hypothetical    protein
BOFHCHGE_00002  CDS rdgB_1  3.6.1.9 dITP/XTP    pyrophosphatase
BOFHCHGE_00003  CDS hypothetical    protein
BOFHCHGE_00004  CDS hypothetical    protein
BOFHCHGE_00005  CDS hypothetical    protein

Or, to add none to any empty field, not just the 3rd one:
awk -F"\t" '{for(i=1;i<=NF;i++){if($i==""){$i="none"; }}}1;' file > newfile

I fear that the .tsv file is mal-formatted then. df<-read.delim("file",header=TRUE, sep="\t") will misimport and EC column will end up having product names, gene having hypothetical protein, etc.

@charlesdarwin if it is indeed malformatted, yes. However, I would be very surprised if that were the case. Using empty fields as described above is standard practice. You can check easily enough though. FInd one of the lines with a missing field (using 12 here as an example), and then run awk 'NR==12' | od -c and make sure you see the \t\t where the empty field is. If you don't see that, this becomes a text parsing problem. I suggest you post a new question with some example lines and we'll see what we can do.

my 6th line has a missing EC number. My file is file.txt and my data frame is called df, what is the command I should run in the terminal with awk?

@charlesdarwin (sorry, I just saw this for some reason) it depends. If there is a \t\t for the empty field, you can do awk -F"\t" '($4==""){$4="none" }1;' file.txt > newfile.txt. If that doesn't work, you might need to ask a new question, or come into the site chatroom and ping me and I'll see if I can help.

newfile.txt is a mess. The element 1,1 is now BOFHCHGE_00001 CDS hypothetical protein none. I've written to the developer. I don't think it should be that difficult. I am using the GFF file instead. It is well formatted and contains all the information, too. Thanks for your help.

@charlesdarwin that's probably for the best. Sorry, but it's really hard to parse something without having the original file to be able to see exactly what's needed.

very easy to reproduce, download the E. coli genome from here: ftp://ftp.ensemblgenomes.org/pub/release-32/bacteria//fasta/bacteria_91_collection/escherichia_coli/dna/Escherichia_coli.HUSEC2011CHR1.dna.toplevel.fa.gz, run prokka on it (conda install prokka) and you'll get a malformed .tsv file.

@charlesdarwin OK. Sorry, but I don't use prokka and I don't want to install it for this. If you want to paste a few lines of the output file in chat or in your question, I can try to help.

Here is a link to the offending file: http://cristian-riccio.ch/wp-content/uploads/2017/12/PROKKA_12142017.tsv

@charlesdarwin oh wow. Yeah, that's definitely a bug, the fields are wrong and since different lines are missing different fields, I dont' see how you can change it unless you do so manually. That's really impressively bad design.

Since the missing values don't have field separators, the easiest option is to rearrange the columns so that the ones with missing values are at the end. They will then be the ones that are filled with missing values when you import them into R. You can use use fill=TRUE with read.table() or use readr, which will complain about missing columns, but fill them in with missing values automatically:
awk -F'\t' -v OFS='\t' '{print($1,$2,$5,$3,$4)}' file.txt |perl -pe 's/\t+/\t/g' >newfile.txt

With readr:
read_tsv("newfile.txt")

# A tibble: 5 x 5
       locus_tag ftype                  product   gene EC_number
           <chr> <chr>                    <chr>  <chr>     <chr>
1 BOFHCHGE_00001   CDS     hypothetical protein   <NA>      <NA>
2 BOFHCHGE_00002   CDS dITP/XTP pyrophosphatase rdgB_1   3.6.1.9
3 BOFHCHGE_00003   CDS     hypothetical protein   <NA>      <NA>
4 BOFHCHGE_00004   CDS     hypothetical protein   <NA>      <NA>
5 BOFHCHGE_00005   CDS     hypothetical protein   <NA>      <NA>

The perl regex is necessary because awk adds tabs after the missing values, which defeats the point of reordering the columns

I did as you said but the data frame is still wrong. On line 6, gene and product are interchanged.

Since you only provided 5 lines in the example, I can't comment on the results for the 6th line

all. I know this question was asked 2 years ago, but I wrote a script in R, with a little bash, to parse the tsv file in the old version of prokka. In the new version of prokka this problem is fixed.
First use bash to get the rows where ftype is CDS.
grep "CDS" all_prokka_annotations.tsv > all_prokka_annotations_CDS.tsv

Next we can parse the TSV with the CDS lines in R.
prokka_annotations = readLines("all_prokka_annotations_CDS.tsv")
prokka_annotations_split = strsplit(prokka_annotations,split="\t")
# first get the protein product name.
prokka_annotations_product = sapply(prokka_annotations_split, function(x) x[length(x)])
# now the I have all protein product names remove them
prokka_annotations_split = sapply(prokka_annotations_split, function(x) x[-length(x)])
#next get locus tag and ftype
locus_tag_cds = sapply(prokka_annotations_split, function(x) x[c(1,2)])
locus_tag_cds = t(locus_tag_cds)
# now remove locus tag and ftype
prokka_annotations_split = sapply(prokka_annotations_split, function(x) x[-c(1,2)])
# next get EC number
ec_number = lapply(prokka_annotations_split, function(x) x[grep(".",x,fixed=T)])
ec_number = as.character(ec_number)
ec_number[ec_number=="character(0)"] <- NA
## now get the gene names if available. to do this remove the EC numbers. whats left are gene names
prokka_gene_names = lapply(prokka_annotations_split, function(x) x[!grepl(".",x,fixed=T)])
prokka_gene_names = as.character(prokka_gene_names)
prokka_gene_names[prokka_gene_names=="character(0)"] <- NA
prokka_parsed_tsv_correct = cbind(locus_tag_cds,product=prokka_annotations_product,gene=prokka_gene_names,EC=ec_number)

```


text:
Can the canu assembler output a fastq file of the final assembly just like HGAP4?
I have assembled some genome from Sequel PacBio data both with HGAP4 on the SMRT Link interface and using canu on the command line. The HGAP4 assembler outputs a fastq file of the final assembly such that I can see how reliable the sequence is along the genome. On the other hand, canu does no such thing, which means I do not know how much I can trust each nucleotide in the assembly? Is there a way canu can output the fastq file as well or at least some depth information at each base?

I've created an issue on the Canu github repository for this. I'm not aware of any existing functionality to output FASTQ files, but think that this would be a useful feature to have.
It would be possible to create something like this by aligning the original reads to the genome, generating a pileup from that alignment, then determining the sequence likelihood based on that pileup. Unfortunately, I don't have any scripts that I've written to do that... yet.


text:
What to do with the configuration test cell after configuring the MinION with it?
What to do with the configuration test cell after configuring the MinION with it?
I received a MinION with a configuration test cell in it. After running a configuration in the MinKNOW software while the CTC was in the MinION, I wonder what to do with the CTC? Am I going to need it later or can I throw it away? If I should keep it, in what conditions should I keep it (temperature, light, etc.)?
Also, do you know why we need to do a configuration of the MinION the first time?
I posted this question on Biology Stackexchange as well but maybe the community here will be more receptive?

I would keep it, in case you need to do another configuration later on. I have used mine for that purpose.
I think storage doesn't really matter, I just put it in a drawer somewhere. Keep it away from liquids and extreme temperatures :)


text:
Does MinKNOW work with Mac OSX high sierra 10.13.1?
I have done a compatibility check as described here on my computer to know if I can use MinKNOW on my mac laptop. It says that my OS is incompatible and that the OSX must be Yosemite or El Capitan. Did anyone manage to make MinKNOW work on High Sierra? I find it odd that they would not support more recent version of mac.

Have you contacted ONT support about this support@nanoporetech.com? They would be in a better place to deal with this, and make sure the issue doesn't happen for other people.

AFAIK Officially Ubuntu 16.04 is also not supported (but 14 is). But it works anyway on our system.

Some replies from ONT is:

OS X is supported and you can view the installation instructions on
  our Downloads page linked below. At the moment High Sierra is not
  officially supported as we have not yet fully tested it, but other
  users in the community have reported no issues with it so far.
The compatibility test only performs a simplistic check to see if you
  are using one of our officially supported OS X versions. Since you are
  on High Sierra, it will show a failure. That does not mean the
  software can't be installed and used. As I said, other users have
  reported no problems, but we can't officially support it if you
  encounter an issue.

So, I am going ahead and installing MinKNOW on my High Sierra.


text:
What are the pros and cons of the different basecallers in Oxford Nanopore Technology Sequencing?
What are the pros and cons of the different basecallers in Oxford Nanopore Technology Sequencing?
I am about to start a MinION run on my laptop. What should I consider when choosing my basecaller? Can I let MinION do its sequencing and generate fast5 files and call the bases later? I think I might use Albacore. Metrichor and MinKNOW can also do basecalling. Are there reasons to choose one over another?

First of all - yes, you can generate FAST5 files and basecall later. Basecalling during the sequencing run is useful if you want results more quickly. You can also recall your FAST5 files with multiple basecallers, if desired.
There are several ways to basecall currently:

MinKNOW
Albacore
Guppy

MinKNOW uses an embedded version of Albacore to perform its basecalling. However, the version in MinKNOW occasionally lags behind the released Albacore version.
Guppy is a GPU enhanced basecaller, and has not yet been released by ONT. It is used by default on the GridION. Results are similar, but not identical to Albacore, but Guppy is much faster if you have a supported GPU.
There used to be a cloud version of basecalling in Metrichor, but this is no longer available, and all basecalling must now be performed locally.
I'd suggest looking over Ryan Wick's comparison of the different available base callers, as he has done a comprehensive look (mainly from a qualitative perspective) at the different basecallers.

There are also third party free and open source basecallers that haven't been developed by Oxford Nanopore.
Of particular note is Chiron, which gave the best uncorrected assembly identity among the base callers that Ryan Wick tested. It's slower than Albacore, but appears to be more customisable, and could theoretically be used to model any DNA feature that is detectable in the electrical signal, rather than just the A/C/G/T base sequence.
There's a paper about Chiron on BiorXiv.

It's worth noting two things as of Dec 2018:

Albacore is being deprecated (but is still available from the Nanopore developer portal).
Guppy is under active development, so Ryan Wick's comparisons may not reflect the current state of things. (The claim is that the current, just-released version of Guppy uses a new "flip-flop" algorithm that improves accuracy over Albacore.)
MinKnow contains the current version of the production basecaller. As of the release today, this has switched from Albacore to Guppy. If you want, you can run MinKnow for your basecalling (which will run Guppy), but you will likely need a beefy computer with access to a good amount of storage to do that.

Depending on your local system I'd suggest to run the MinION without local basecalling, and perform basecalling later using Albacore on server/cluster. That will keep storage requirements on your nanopore-computer low. Have albacore output fastq files which will be okay for most typical applications. Keep the fast5 files to repeat basecalling with newer versions of basecallers and to call nucleotide modifications.


text:
Optimal design for an RNA-seq experiment with ERCC RNA Spike-In Control Mixes
I am planning an experiment on gene expression. I have 6 conditions and 3 replicates in each condition. I would like to use ERCC spike-in mixes for quality check: lower limit of detection, dynamic range and also check that the fold-changes are accurate. For example, by knowing that spike-in sequence xyz in condition 2 should be a 4-fold higher than in condition 1, does it actually get detected as a 4-fold change in my analysis? To assess the fold-changes, I need to use 2 different mixes. Only one mix should be used per sample. The manual suggests either adding mix 1 in a reference and mix 2 in all other samples or randomly assigning mixes 1 and 2 among samples. The latter is not very clear to me. Should I use the same mix for all replicates within a condition or should I randomise really among all the samples? Or can I randomise among conditions (for example 1, 3, 4 get mix 1 and 2, 5, 6 get mix 2)? So far, I am intending to do the latter, so using mix 1 in conditions 1, 3, 4 and mix 2 in conditions 2, 5, 6. This would also be cheaper because I can only buy equal volumes of mix 1 and 2. Do people have any advice on this?

I have 6 conditions and 3 replicates in each condition.

If including more biological replicates has a low experimental cost, you should be using 6 replicates per condition (and fewer reads per condition, i.e. maintaining the total number of reads) in order to pick up the majority of differentially-expressed genes. For that, I'd recommend using 3 lanes of a flow cell, with 12 dual-indexed barcodes (with the same index on each end) for each lane. Doing dual indexing will protect a bit against index switching and make the results more robust.

randomly assigning mixes 1 and 2 among samples. The latter is not very clear to me. Should I use the same mix for all replicates within a condition or should I randomise really among all the samples?

Randomly assign mixes among and within the conditions. Make sure that the mixes are equally represented (as much as possible) across the conditions, but that each condition has at least one of each mix. This can be done using something like the following in R:
numReplicates <- 7;
numEachMix <- ceiling(numReplicates / 2);
assignments <- sample(rep(c("mix1","mix2"),numEachMix))[1:numReplicates]

For working out the level of fold change this is reliably replicated (or "statistically significant"), you can split the samples into mix 1 and mix 2, then apply the same algorithms for ERCC differential expression as with the non-ERCC genes.
However, as ERCC is typically added in proportion to the amount of RNA (or, if you're lucky and have cell-sorted populations, possibly the number of cells), all it's really testing is the accuracy of sample prep.

Thank you for this very clear answer. If I randomise as you say, what happens if my pseudo-random number generator assigns the same mix to all the replicates of the same condition? Shall I just switch one over to the other mix in that case? In the analysis of the fold coverages, shall I just separate my samples into 2 groups: mix 1 and mix2 and do fold-change analysis between the two groups?

If all replicates of the same condition have the same mix, then switch one of them over -- "each condition has at least one of each mix", but it'd be better to distribute the numbers of mixes equally in each group.

so, with 3 replicates, shall I do mix 1, mix 2 and mix random?

Yes, but randomly assigned. In R, you could do this by something like the following:
c("mix1","mix2","mix1","mix2")[sample(4)][1:3]

that's quite clever, I will also set a seed with set.seed(date + 1, 2, 3, etc.)


text:
How to download the whole BLAST nt database into a specific folder?
I have successfully downloaded the whole nt BLAST database into the current folder using:
wget -b "ftp://ftp.ncbi.nlm.nih.gov/blast/db/nt.??.tar.gz"

However, I would like to download it in a specific folder, e.g. output/
I have tried:
wget -b "ftp://ftp.ncbi.nlm.nih.gov/blast/db/nt.??.tar.gz"  -O output/

but it does not work. Nothing happens. My 'wget-log' is empty.
I could move to that folder and download it there but I would like to remain in the same working directory.
C.

-O is to specify the name of the output file. You want the -P option.
So, from your example
wget -b "ftp://ftp.ncbi.nlm.nih.gov/blast/db/nt.??.tar.gz" -P output/

I recommend using NCBI's tool for this called update_blastdb.pl, shipped with the blast suite (I use the conda environment from bioconda).
cd to the directory you want nt DB in and run:
update_blastdb.pl --decompress nt

Pros: if a download fails for whatever reason it is retried, downloads are md5sum checked and re-running in a folder with an existing DB will (as the name suggests) only update the relevant DB portions. More information here!


text:
How to download the whole directory of an ensembl FTP page?
I am interested in downloading all the files at:
ftp://ftp.ensembl.org/pub/release-91/embl/bos_taurus/
I know how to download a single file from that page using wget (Unix) or download.file() (R) but is there an easy way of downloading all of them?

You can use the mirroring function of wget:
wget -m ftp://ftp.ensembl.org/pub/release-91/embl/bos_taurus/


text:
Problem installing RATT
I have followed the steps here to install RATT, i.e.:
svn co "https://svn.code.sf.net/p/ratt/code/" ratt-code

However, when I check my installation by running '$RATT_HOME/start.ratt.sh', it says:
zsh: no such file or directory: /start.ratt.sh

which makes sense once I check that the RATT_HOME variable is empty.
I downloaded their tuberculosis example, went into that directory and ran:
~/path/to/my/ratt/ratt-code/start.ratt.sh embl F11.fasta F11 SimStrain

and got the following output:
Please set the RATT_HOME variable.
At Sanger for bash it is RATT_HOME=/nfs/users/nfs_t/tdo/Bin/ratt; export RATT_HOME
At Sanger for tcsh setenv RATT_HOME /nfs/users/nfs_t/tdo/Bin/

What should I set my RATT_HOME variable to? The message only says what it should be set to at the Sanger Institute.

Usually the home of a program is where the executable is. In order to make the system aware you need to include it as a variable when the terminal is run.
So you need to add it to the file at ~/.bashrc of your computer (if you are using bash as shell [you can check it with echo $SHELL]). You need to add :
export RATT_HOME=path/where/RATT/is/installed/or/downloaded

As terdon says in a comment below:

Note that echo $SHELL just prints the default login shell of the
  user, and not necessarily the shell that is currently running.  Also,
  ~/.bashrc is not read by macOS bash which starts login shells by
  default. In any case, this sort of global variable should be in
  ~/.profile, not .bashrc. For more of this sort of detail see
  Sequence of scripts sourced upon login and How to correctly add
  a path to PATH.


text:
Why does cutadapt remove low quality bases from the ends of reads only?
I use cutadapt to remove low quality bases from my Illumina reads. The algorithm only removes low quality bases from the end until it reaches a good quality base. If there is a bad quality base beyond that, it is not trimmed. Why? Why doesn't the algorithm remove all low quality bases? It could replace low quality bases in the middle of a read by an N, for example. Or is it that the low quality base is still most probably correct and hence one does not want to lose the base information for mapping?

I am commenting on this part:

The algorithm only removes low quality bases from the end until it reaches a good quality base. If there is a bad quality base beyond that, it is not trimmed.

According to its user guide, cutadap is designed this way: it trims off bases from the 3'-end until it sees a base with quality higher than a threshold. This is not a good algorithm. For example, when we see a quality string 30,30,30,3,3,3,3,3,20,3,3, we would prefer to trim it down to 30,30,30 because a 3,3,3,3,3,20 tail is still not useful.
Some tools such as trimmomatic uses a sliding window to avoid this issue. In my view, a better algorithm is the so-called the "modified Mott algorithm" used by phred 20 years ago. seqtk among others implements this algorithm. Note that the Mott algorithm always trims from both ends, which is sometimes not preferred. BWA implements a variant of this algorithm, trimming from the 3'-end only.
In practice, though, different quality trimming algorithms probably work equally well because it is relatively rare to see a high-quality base in the middle of a low-quality tail. In addition, modern aligners can well handle low-quality tails; assemblers and error correctors can correct through such tails, too. It is usually not necessary to apply quality trimming.

It could replace low quality bases in the middle of a read by an N, for example. Or is it that the low quality base is still most probably correct and hence one does not want to lose the base information for mapping?

With a low-quality base, you might have an error/mismatch, but with an N, you always have a mismatch. Masking low-quality bases to N is not as good.
EDIT: responding to the following comment from OP:

What if the base is low-quality and then it matches the reference by mistake? Wouldn't it be better if it did not match and was penalized? I guess this should be modeled mathematically.

If a read could be mismapped due to a sequencing error, its true location is often in the genome. In this case, a competent mapper will give the alignment a low mapping quality (mapQ). A quality-aware mapper such as novoalign will penalize mapQ even more. In comparison, if you hard mask this sequencing error to "N", you will get a mapQ=0. You can see the difference between the two approaches mostly comes from mapQ.
In modern Illumina data, it is quite often to see a Q8 base in the middle of high-quality bases. >80% of them (in theory) are still correct. My hunch is that hard masking all of them would lead to considerable data loss.

Is it possible to expand on the last part? 'With a low-quality base, you might have an error/mismatch, but with an N, you always have a mismatch. Masking low-quality bases to N is not as good.' What if the base is low-quality and then it matches the reference by mistake? Wouldn't it be better if it did not match and was penalised? I guess this should be modelled mathematically.

If you check the read QC statistics of an Illumina run in e.g. fastQC, you will see that at the end of the read the quality decreases. This is because of exhaustion of chemicals at the end of the run. This is a general trend seen in all runs, therefore you can remove these low quality bases from the end of your run. If you have incidentally a bad quality base in the middle of a read this is not the same general trend but is sporadic. Which is hard to automatically remove with tools such as cutadapt. If you really want an N for each low quality bases within the read, you can make a simple program for that, but why would you need that? Most aligners can handle a few bad quality bases within a read, and for variant calling you'll need many reads (so the incidental error base will be overruled by the good ones).

It's also worth noting that many aligners penalize an N by default, so it's not like replacing super-low quality bases with an N is going to greatly boost alignment scores.

For Illumina (and 454) reads, the quality decreases with read length.  It's not linear and is run/library-dependent.  It has less to do with exhaustion of reagents and more to do with the strands on a spot being out of phase due to incomplete/missed base incorporation during sequencing.
It is common practice to trim off the low quality 3' ends as the first step in QC.  Most analyses and most can deal with a few low errors in the middle of the read and what you do about such sequencing errors may vary depending on what type and use of the data.  But in general, nobody wants the garbage at the end of the reads.
P.S. Even for the same library prep, the fwd/rev runs will have different quality curves.
P.P.S. Sometimes the first couple bases are also poor too.

"in general, nobody wants the garbage at the end of the reads." For RNA-seq, quality trimming may be ok, but for de novo assembly, trimming may reduce N50.

If you don't trim prior to assembly then a DBG assembler will have a lot of extra kmers, blowing up RAM.  It's been my experience that trimming with increase N50, not decrease it.  Although I don't think N50 is a great metric to evaluate assembly quality.  I suppose as an alternative to end-trimming, one could do error correction, e.g. with bfc.

Modern assemblers (e.g. allpaths, sga and spades) often come with an error correction step. They only trim an end when they can't correct it. Only a small fraction of rubbish k-mers can go into the graph.

Thanks, user172818, you're probably right.  But we have found it beneficial to use bfc prior to metaspades.


text:
Can I index a compressed FASTA file using STAR?
I am using STAR to align RNA-seq reads to a reference genome. Before the alignment, I need to generate an index of the reference genome. I use the following code to generate the index successfully:
STAR --runThreadN 8 --runMode genomeGenerate --genomeDir output/index/star --genomeFastaFiles ref.fa --sjdbGTFfile ref.gtf --sjdbOverhang 100
This works fine. However, I would like to keep my reference genome compressed to save disk space. So I am trying the following command:
STAR --runThreadN 8 --runMode genomeGenerate --genomeDir output/index/star --genomeFastaFiles ref.fa.gz --readFilesCommand "gunzip -c" --sjdbGTFfile ref.gtf --sjdbOverhang 100
but I get the following error:
EXITING because of INPUT ERROR: the file format of the genomeFastaFile: ref.fa.gz is not fasta: the first character is '' (31), not '>'.
 Solution: check formatting of the fasta file. Make sure the file is uncompressed (unzipped).
I am using the readFilesCommand successfully with compressed RNA-seq fastq files. Does anybody know if there is a similar option to use compressed references? Is there a workaround using Unix commands (maybe piping?) or do I need to decompress the reference, index it and then compress it again?

I can think of two possible workarounds. The simplest would be to try using zcat instead of gunzip -c, in case that works. Assuming it doesn't, you can use process substitution, if your shell supports it. Something like this:
STAR --runThreadN 8 --runMode genomeGenerate --genomeDir output/index/star \
    --genomeFastaFiles <(zcat ref.fa.gz) --sjdbGTFfile ref.gtf --sjdbOverhang 100

The <() creates a file descriptor and passes that to the command. It is very useful for passing streams to commands that expect files.
If your shell doesn't support this construct, assuming you are on a system that supports named pipes (Linux does), you can do the same thing a bit more manually:

Create a named pipe (FIFO)
mkfifo foo

Set it to read the uncompressed genome
zcat ref.fa.gz > foo &

The & is important since the command above will hang until something starts reading from foo. Alternatively, you can just run this in one terminal and the next command in another.
Read the input from the pipe
STAR --runThreadN 8 --runMode genomeGenerate --genomeDir output/index/star \
    --genomeFastaFiles foo --sjdbGTFfile ref.gtf --sjdbOverhang 100


text:
Can I use my computer while MinION is sequencing without negatively affecting the run?
I have a Mac Book Pro and I am about to start a sequencing run on the MinION. MinKNOW is going to be running for 2 days.
Can I use my computer during sequencing? Can I browse the web and use Excel, etc? I am not speaking about running a parallelised genome assembly on all my cores using my entire RAM but more mundane activities, like checking my emails or taking notes. Furthermore, can I lock my screen using Control + Shift + OFF without negatively affecting the sequencing?

By locking do you include suspending the machine or not?

Yes, that should be fine, especially if basecalling is disabled. Locking a screen should be fine as well.
I prefer to do basecalling after a run anyway, because there's a chance that it could chew up too much processor time and result in lots of skipped reads.

Since biological samples and reagents cost time and money, I wouldn't push it too much.  Checking your email would be fine, but I wouldn't be browsing the web (besides gmail/calendar) in case a plug-in (e.g. flash), bad javascript, etc. causes any trouble.
I'd disable any auto-updates too.
Locking the screen would be fine.
Personally, I'd buy a cheap chromebook and let the Mac work.


text:
Can I stop my nanopore sequencing run if there are no more reads being produced?
I am sequencing a whole genome on MinION. I have used this flow cell for 6 hours for the phage lambda control. That gave me around 300,000 reads. I started a 48h sequencing run on the same flow cell after washing it. The sequencing has been going on for 1 day and 21 hours. It is not producing any reads at the moment (Strand = 0). The number of single pores is also 0. Can I "Stop acquisition" without negatively affecting the results?
The image below is a screenshot taken from MinKNOW, the software used to control the MinION.

Yes, you can click "stop acquisition" and your run won't be negatively affected. All of the reads are saved as they are generated. I am not sure how this will impact live basecalling though if that is something that you do.

Live basecalling for the uncalled reads will continue until all have been called, or the application is quit.

Yes, you can stop the run and it should not affect the base calling. We normally do the runs and then do the base calling after the run has completed.


text:
Does the DNA or RNA of some ethnic groups map better than others' to the human reference genome sequence?
I believe that the human genome reference took DNA samples from different people and that some natural variation is included in extra contigs.
However, the human reference genome comes from a limited number of people.
Can the limited number of donors and the linear representation of the genome (as opposed to a graph representation that could potentially include variation) bias an analysis in favour of one ethnic group over another? Are there studies that have compared the alignment rates/mappabilities of reads of people coming from different geographic regions? Do some ethnic groups' reads align better to the reference?

For simple variants like SNPs it would not really be a problem to use the current genome assembly for other ethnic groups. But for more complex variants this could be indeed problematic, however not only for ethnic groups but also for individuals within the same population. Think of very complex regions such as HLA or KIR.
In studies where they compare different ethnic groups (e.g., 1000 G) they rather not discuss such matters, but act like the analysis is 100% perfect.
There is a nice comment in Genome Biology about this issue.


text:
How can I transform a mapped BAM file into an unmapped BAM file?
In order to use MergeBamAlignment (Picard), I need an unmapped BAM file and a mapped BAM file. I have two mapped BAM files:

one with reads mapped to the reference I want but without metadata such as read group
a second one with reads mapped to the incorrect reference but with metadata

When I try to use MergeBamAlignment with these two files, I get errors, such as "Second read from pair not found in unmapped bam". So I am trying to unmap the reads from my second BAM file. How can I do that? Is there a tool for it? Otherwise, I'll have to manually change the bitwise FLAGs, reference names, CIGARS, etc.

Looks like you need Picard's RevertSam operation (described in context here, manual here).
java -jar picard.jar RevertSam \
     I=input.bam \
     O=reverted.bam

What happens if you use the bam with the metadata as the input to your aligner?  Maybe it will keep all those fields.

bwa mem requires a fastq file, it does not accept a BAM

Okay, I didn't know you had to use that.  You see that bwa aln will take a bam?

good point, the RevertSam worked though


text:
How can I ask snakemake to produce a dag where each node represents a rule?
I am using snakemake to create workflows. It is very convenient to visualise my DAG using snakemake --dag target{sampleA,sampleB,sampleC}File | dot -Tpdf > dag.pdf. The resulting pdf shows all the rules' dependencies to get to the target files. However, there is one node per rule/sample combination in the pdf image. When I have many samples, this creates an unreadable image with little added value. Is there an option to merge all the nodes from the same rule together so that I have the workflow for just one sample?

I think you are looking for --rulegraph instead of --dag
  --rulegraph           Do not execute anything and print the dependency graph
                        of rules in the dot language. This will be less
                        crowded than above DAG of jobs, but also show less
                        information. Note that each rule is displayed once,
                        hence the displayed graph will be cyclic if a rule
                        appears in several steps of the workflow. Use this if
                        above option leads to a DAG that is too large.


text:
Can blat use more than one core/CPU to speed up the alignment?
I am using BLAT to align two versions of the genome of C. elegans. I can see in the Activity Monitor of my Mac Book Pro High Sierra that blat is using 100% of a CPU. However, is this programme able to use more than one core / CPU to speed up the alignment?
I cannot find such option in the help page of the tool here.
I installed blat using conda install -c bioconda blat.
I am using the latest version of blat, i.e. 36.

BLAT can only use one CPU. It is actually not the right tool for full-genome alignment. For "two versions" of the same species, MUMmer and minimap2 are orders of magnitude faster and probably give better alignment.
EDIT (moving comment to answer): OP comments that the purpose of this alignment is for lifting over annotations using the UCSC PSL-based pipeline. So far there are no converters for minimap2's SAM/PAF or MUMmer's delta formats. With the UCSC pipeline, you have to write a converter by yourself. An alternative approach is to use minimap2's liftover. The minimap2 cookbook gives an example:
minimap2 -cx asm5 --cs ecoli_ref.fa ecoli_canu.fa > ecoli_canu.paf
echo -e 'tig00000001\t200000\t300000' | paftools.js liftover ecoli_canu.paf -

Note that this pipeline is not as mature as UCSC liftOver and is not as scalable.

I need an alignment in PSL format for downstream tools needed for lifting over an annotation. Do you know if I can get an alignment in psl format with nucmer or minimap2?

@charlesdarwin LASTZ tool might also be worth a look for this, although it's quite old now...


text:
How can I specify an email address in the SLURM snakemake json configuration file?
I am working on a SLURM-based cluster. I am currently using the following slurm json file:
{
    "__default__" :
    {
        "A" : "overall",
        "time" : "96:00:00",
        "nodes": 1,
        "ntasks": 1,
        "cpus" : 24,
        "p" : "LONG",
        "mem": 16G,
        "output": "snakemake%A.out",
        "error": "snakemake%A.err"
    }
}

along with the following command:
snakemake --jobs 200 --use-conda --printshellcmds  --cluster-config slurm.json --cluster "sbatch -A {cluster.A} -p {cluster.p} -t {cluster.time} --output {cluster.output} --error {cluster.error} --ntasks {cluster.ntasks} --cpus-per-task {cluster.cpus} --mem {cluster.mem}"

I would like to be able to specify my email and receive an email when a job finishes. I have tried adding a line in the json file with "email": my@email.com and an option to the sbatch command --email {cluster.email} but that did not work.
Anyone knows how to do it?

After looking at some headers of files (#SBATCH --mail-user=email@comI have used --mail-user instead. I'll let you know if I receive an email.

I'd be surprised if --mail-user didn't end up working. Note also --mail-type.

so --mail-user on its own didn't work but by adding the --mail-type, it worked!

As noted by the comments and the SLURM manual (not so much bioinformatics related?)
You need the following terms in your batch/queuing command:
--mail-type=ALL
--mail-user=you@email.com

I used "ALL" which will send you all reports. If you just want the finish message then you would use "END" and "FAIL" like such:
--mail-type=END,FAIL

Here are the specifics from the manual if you wish to make a different selection:
--mail-type=

Notify user by email when certain event types occur. Valid type values
  are NONE, BEGIN, END, FAIL, REQUEUE, ALL (equivalent to BEGIN, END,
  FAIL, REQUEUE, and STAGE_OUT), STAGE_OUT (burst buffer stage out and
  teardown completed), TIME_LIMIT, TIME_LIMIT_90 (reached 90 percent of
  time limit), TIME_LIMIT_80 (reached 80 percent of time limit),
  TIME_LIMIT_50 (reached 50 percent of time limit) and ARRAY_TASKS (send
  emails for each array task). Multiple type values may be specified in
  a comma separated list. The user to be notified is indicated with
  --mail-user. Unless the ARRAY_TASKS option is specified, mail notifications on job BEGIN, END and FAIL apply to a job array as a
  whole rather than generating individual email messages for each task
  in the job array.

--mail-user=

User to receive email notification of state changes as defined by
  --mail-type. The default value is the submitting user.


text:
How do I tell featureCounts to ignore fragments that overlap with more than one feature?
I am using featureCounts from the subread package to count the number of fragments per gene in an RNA-seq experiment. I feed it a BAM file with paired-end reads aligned to a reference genome. In the help page of the executable, which I have installed using conda install subread, there is a section about overlaps, including for example the -O option which will assign the reads/fragments to all overlapping features. What is the default behaviour? Assigning to reads/fragments to the feature with largest overlap? not assigning multi-overlapping reads/fragments?
I would appreciate your help on this.
OS: macOS X
featureCounts version: 1.6.3

If a fragment overlaps multiple genomic features then it's unassignable, which is the appropriate way to handle ambiguous mappings when the resulting counts are then given to tools like DESeq2.

fair enough, so I will not use the -O option, thanks


text:
TRAL does not find "Phobos result file"
I want to use TRAL to annotate tandem repeats in the reference genome of Caenorhabditis elegans. For this, I need to install some external software, such as Phobos. I've downloaded Phobos and I am using it by typing the path to it's executable like this:
input/software/phobos-v3.3.12-mac/bin/phobos-mac-intel-64bit-start-from-command-line genome.fa output.txt
This works quite well.
In order for TRAL to find Phobos, I changed TRAL's configuration file as explained here.
I've replaced PHOBOS = phobos with PHOBOS = /Users/user/Documents/.../project/input/software/phobos-v3.3.12-mac/bin/phobos-mac-intel-64bit-start-from-command-line
Then, I followed the tutorial Run and parse de novo repeat detection software and typed the following commands in Python 3 in the Terminal:
import os
from tral.sequence import sequence
from tral.paths import PACKAGE_DIRECTORY

proteome_HIV = os.path.join(PACKAGE_DIRECTORY, "examples", "data", "HIV-1_388796.faa")
sequences_HIV = sequence.Sequence.create(file = proteome_HIV, input_format = 'fasta')

So far so good. But the following command creates an error:
tandem_repeats = sequences_HIV[0].detect(denovo = True, detection = {"detectors": ["Phobos"]})

Did not find Phobos result file in /var/folders/z2/rt6qxntd3tg5bjmqlqnm0mdm0000gn/T/tmpvc20y3i4/001/PHOBOS/phobos.o

Note that if I use PHOBOS instead of Phobos in the command, I get the following error:
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/user/anaconda/envs/python3env/lib/python3.6/site-packages/tral/sequence/sequence.py", line 162, in detect
    **kwargs['detection'])[0]
  File "/Users/user/anaconda/envs/python3env/lib/python3.6/site-packages/tral/sequence/repeat_detection_run.py", line 962, in run_detector
    Detectors(detectors, sequence_type)
  File "/Users/user/anaconda/envs/python3env/lib/python3.6/site-packages/tral/sequence/repeat_detection_run.py", line 872, in Detectors
    FINDER_LIST.keys()))
Exception: Unknown TR detector supplied (Supplied: ['PHOBOS']. Known TR detectors: dict_keys(['HHrepID', 'Phobos', 'TRED', 'T-REKS', 'TRF', 'TRUST', 'XSTREAM']))

I am working on Mac OS X. I installed TRAL with pip install tral.

I can reproduce this, so I think it might be a TRAL bug. See https://github.com/acg-team/tral/issues/12

Are you running on DNA or protein? In general, TRAL support for DNA is weaker than we'd like it. If you're searching exons, I might try performing the translation first and running in AA mode.

Running on DNA. I am interested in DNA repeats. Mini-and microsatellites are unlikely to code for anything. If I run it on the translation, then the repeats are not bound to appear at the DNA level because of the genetic code degeneracy.

OK. Just make sure you review the results carefully when using TRAL on DNA, as the parameters were trained for AA. We have a student working on improving DNA repeats now, so this should get easier.

TRAL 0.3.5 has bugs running and parsing Phobos. You can follow progress towards a fix on this issue: https://github.com/acg-team/tral/issues/16
Until the issue is fixed, you should remove Phobos from the repeat detector list in your tral config file (~/.tral/config.ini) and use a different repeat detector.
[sequence]
    [[repeat_detection]]
        DNA = TRED, T-REKS, TRF, XSTREAM


text:
How to transform a DNAStringSet from the Bioconductor package Biostrings to a data frame?
I am working on Mac OS X. I am using R version 3.5.1.
I have important a FASTA file into R using Biostrings::readDNAStringSet. This creates a DNAStringSet object in R.
I would like to transform it into a data frame so that I can then plot the widths with ggplot2.
as.data.frame on my DNAStringSet object returns a data frame with only one column. Is there a function that returns a data frame with columns width, seq and names present in the DNAStringSet object?

According to the documentation (?Biostrings::DNAStringSet):

width(x): A vector of non-negative integers containing the number of
            letters for each element in x.  Note that width(x) is
            also defined for a character vector with no NAs and is
            equivalent to nchar(x, type="bytes").
names(x): NULL or a character vector of the same length as x
            containing a short user-provided description or comment for
            each element in x.  These are the only data in an
            XStringSet object that can safely be changed by the user. All
            the other data are immutable!  As a general recommendation,
            the user should never try to modify an object by accessing
            its slots directly.
as.character(x, use.names=TRUE): Converts x to a character vector
            of the same length as x.  The use.names argument controls
            whether or not names(x) should be propagated to the names
            of the returned vector.

This should be enough for you to create a data.frame:
dss2df <- function(dss) data.frame(width=width(dss), seq=as.character(dss), names=names(dss))

Also, if you have a huge DNAStringSet, perhaps you do not want to convert all the sequences to characters; storing sequences as characters supposedly is less efficient. For plotting it might be sufficient to have just widths.
qplot(width(dss), geom='histogram'))


text:
Where can I find the DNA sequence of the C. elegans rRNA cluster?
Where can I find the DNA sequence of the C. elegans rRNA cluster?
Ideally, I would be able to do it using the command-line.
I've tried using NCBI but the results are unsatisfactory. I am interested in the DNA sequence that is in the current reference genome of WormBase.
I have tried using the Bioconductor biomaRt package and create a table containing the description but when I grep 'rRNA', I only get proteins doing things with rRNA. I used the following code:
ensembl <- biomaRt::useMart('ensembl', dataset = 'celegans_gene_ensembl')
wormGenes <- biomaRt::getBM(attributes = c('external_gene_name', 'chromosome_name', 'start_position', 'end_position', 'gene_biotype', 'description'), mart = ensembl)
wormGenes[grep('rRNA', wormGenes$description), 6]

[1] "rRNA 2'-O-methyltransferase fibrillarin  [Source:UniProtKB/Swiss-Prot;Acc:Q22053]"
[2] "Pre-rRNA-processing protein pro-1  [Source:UniProtKB/Swiss-Prot;Acc:Q22006]"
[3] "Pre-rRNA-processing protein TSR1 homolog  [Source:UniProtKB/Swiss-Prot;Acc:Q19329]"
[4] "rRNA methyltransferase 2, mitochondrial  [Source:UniProtKB/Swiss-Prot;Acc:O62251]"
[5] "Probable rRNA-processing protein EBP2 homolog  [Source:UniProtKB/Swiss-Prot;Acc:Q09958]"
[6] "Putative rRNA methyltransferase  [Source:UniProtKB/TrEMBL;Acc:Q5CCJ2]"

The info is actually in the gene_biotype column of my wormGenes data frame. I could select the rRNA genes using dplyr::filter(wormGenes, gene_biotype == 'rRNA'). Unfortunately, no description is available for the genes, so I will have to get it using the gene IDs.


text:
GFF file too big to load into R
I want to load a GFF (annotation) file into an R data frame in order to extract some information from it, e.g. the location of transposons in the genome.
I used the following code:
library(rtracklayer)
gff <- rtracklayer::import('c_elegans.PRJNA13758.WS263.annotations.gff3.gz')
Error: cannot allocate vector of size 210.2 Mb

It seems that the GFF3 file is too big for loading into R. Do you know if this is limited by R or by my computer's memory?
I ran this on a GNU/Linux computer.

It depends of course on your available RAM on your computer.
Under Windows you can increase the available RAM in R with memory limit:
> memory.limit()
[1] 8070
> memory.limit(size=10000)
[1] 10000
> memory.limit()
[1] 10000

But that is not necessary on linux/unix.

This only works (and is only required) on Windows. OP is using Linux. Not sure why he accepted this answer.

Indeed, memory.limit()
[1] Inf
Warning message:
'memory.limit()' is Windows-specific. I understood this as 'no memory limit' on Linux, so only limited by available RAM?


text:
How can I run the fdnadist command from the EMBOSS package?
My operating system is Mac OS X.
I would like to compute DNA distances between different sequences. For this, the dnadist package from the Phylip software would be good but hard to automate. EMBOSS seems to have an fdnadist command for this. I have installed emboss version 6.6.0 from the bioconda package but if I type fdnadist, the command is not found. How can I use dnadist or fdnadist on the command line?

To find emboss packages type wossname, then hit enter twice.  You will see all the programs.  There is one called distmat.  Then do tfm distmat, to get the fkn manual for distmat.

I found fdnadist in embassy-phylip package https://anaconda.org/bioconda/embassy-phylip


text:
Y Chromosome Aligned Reads in scATAC-seq data from a female-derived cell line?
I'm working with scATAC-Seq data on the K562 cell line, which is supposed to be derived from a female patient. While following the scATAC-seq data analysis pipeline, after performing bowtie alignment they recommend filtering out all reads aligned to mitochondrial DNA and on the Y chromosome.
Out of curiosity, I decided to count how many reads aligned to chrY and found that it can be quite high - chrY has approximately 10% as many reads as a similarly-sized chromosome (chr19 in hg19).
Is it normal to find reads mapping to chrY when the cell line is female-derived, or did I mess up somehow?

There are homologous regions between X an Y chromosomes: https://en.wikipedia.org/wiki/Pseudoautosomal_region
It is therefore normal to have some female-derived reads mapping in Y chromosome.
You should probably check what proportion of such reads fall in other parts of the Y chromosome than pseudoautosomal regions.

This is it -- I just noticed that most prominent peak corresponds to PRKY, a peudogene in the pseudoautosomal region that's similar to PRKX. Not all of my reads in chrY were pseudoautosomal, but all of my MACS2-called peaks were. Thanks!


text:
How to get a list of genes corresponding to the list of SNPs (rs ids)?
Is there a way for me to get a list of genes given a list of SNP rs ids? I found several questions asked with a similar goal years ago, and the answers are always about using multiple online tools, and/or R programming.
I wonder if there is a simpler solution recently? (and hopefully, I can solve it within the scope of python. )
The relevant question I found (that are asked and answered more than five years ago.)

https://support.bioconductor.org/p/39278/
https://www.biostars.org/p/10437/
https://www.biostars.org/p/413/

Actually functional annotation is a key step in SNP calling pipelines. I have used SnpEff and it worked perfectly. Even if you have to build your own data base given the reference genome and the gene models against the SNPs were called. I have a script if you want to use it as guide.

I wonder if there is a simpler solution recently? (and hopefully, I can solve it within the scope of python. )

A simpler solution, I don't know... but this is at least one Python solution using Biopython's ELink method via NCBI's Entrez E-utils.
The Biopython library is flexible enough (they have an in-depth tutorial worth reading) to modify the code below to to fit project-specific needs.  For example, the mapping between snp id and gene id is lost with this bulk request method.
from Bio import Entrez
Entrez.email = "youremail@gmail.com"

snp_ids = ["rs8179065",
           "rs12899845",
           "rs8179066",
           "rs80357332",
           "rs80356920"]

record = Entrez.read(Entrez.elink(dbfrom="snp",
                                  id=",".join(snp_ids).replace('rs', ''),
                                  db="gene"))

This returns a record object
[{'LinkSetDb': [{'DbTo': 'gene', 'LinkName': 'snp_gene', 'Link': [{'Id': '374654'}, {'Id': '90381'}, {'Id': '672'}]}], 'ERROR': [], 'LinkSetDbHistory': [], 'DbFrom': 'snp', 'IdList': ['8179065', '8179066', '12899845', '80356920', '80357332']}]

which contains several iterable fields -- 'Link' has a list of the Gene UID's that can be converted to names using ESummary.
for gene_id in record[0]['LinkSetDb'][0]['Link']:
    handle = Entrez.esummary(db="gene", id=gene_id['Id'])
    uid_record = Entrez.read(handle)
    handle.close()

    uid_summary = uid_record["DocumentSummarySet"]['DocumentSummary'][0]
    print (gene_id['Id'], uid_summary['Name'])
    # ...or something more interesting than print

374654 KIF7
90381 TICRR
672 BRCA1

An option to preserve discrete mapping, is to loop over each snp id and extract the gene(s):
although I don't endorse the method below, especially for long lists of snps, you may flood the server with too many requests, maybe someone more familiar can comment on request limits or best practices for Entrez
for snp_id in snp_ids:
    record = Entrez.read(Entrez.elink(dbfrom="snp",
                                  id=snp_id.replace('rs',''),
                                  db="gene"))
    results = record[0]['LinkSetDb'][0]['Link']
    for result in results:
        uid = result['Id']
        handle = Entrez.esummary(db="gene", id=uid)
        uid_record = Entrez.read(handle)
        handle.close()
        uid_summary = uid_record["DocumentSummarySet"]['DocumentSummary'][0]
        gene_name = uid_summary['Name']

        print(snp_id, uid, gene_name)

rs8179065 374654 KIF7
rs12899845 374654 KIF7
rs12899845 90381 TICRR
rs8179066 374654 KIF7
rs80357332 672 BRCA1
rs80356920 672 BRCA1

I can show you a simple way in R, using biomaRt. Let's say you have two snp ids you want to exmine.
SNPids <- c("rs431905509", "rs431905511")

library("biomaRt")

# To show which marts are available
listMarts()

# You need the SNP mart
mart <- useMart("ENSEMBL_MART_SNP")

# Find homo sapiens
listDatasets(mart)

# This will be the dataset we want to use
dataset <- useDataset("hsapiens_snp", mart=mart)

# Show available filters
listFilters(dataset)

# Now list all available attributes
listAttributes(dataset)

# To get the ensembl gene id belonging to the SNPs
getBM(attributes=c('refsnp_id', 'ensembl_gene_stable_id'),
      filters = 'snp_filter',
      values = SNPids,
      mart = dataset)

    refsnp_id ensembl_gene_stable_id
1 rs431905509        ENSG00000070371
2 rs431905509        ENSG00000100075
3 rs431905509        ENSG00000260924
4 rs431905511        ENSG00000145335

If you want a fast (binary) search with a locally-stored dataset, where you know exactly what SNPs and gene annotations are being queried, here's a way to set up files to do your own querying, which should scale up nicely if you have a lot of SNPs.
1) Get annotations of interest and write them into a sorted BED file.
For example, here's a way to get Gencode v26 gene annotations:
$ wget -qO- ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_26/gencode.v26.basic.annotation.gff3.gz \
    | gunzip -c \
    | convert2bed --input=gff - \
    > gencode.v26.bed

2) Get SNPs and write them into a text file sorted by SNP ID.
For hg38, for instance:
$ LC_ALL=C
$ wget -qO- http://hgdownload.cse.ucsc.edu/goldenpath/hg38/database/snp147.txt.gz \
   | gunzip -c \
   | awk -v OFS="\t" '{ print $5,$2,$3,($3+1) }' \
   | sort -k1,1 \
   > hg38.snp147.sortedByName.txt

3) Install pts-line-bisect:
$ git clone https://github.com/pts/pts-line-bisect.git
$ cd pts-line-bisect && make && cd ..

4) Run a query:
$ rs_of_interest=rs2814778
$ ./pts-line-bisect/pts_lbsearch -p hg38.snp147.sortedByName.txt ${rs_of_interest} \
   | head -1 \
   | cut -f2- \
   | bedmap --echo --echo-map-id-uniq --range 1000 - gencode.v26.bed \
   > answer.bed

The file answer.bed contains the genomic interval of the SNP of interest, along with the IDs of any Gencode v26 gene annotations within 1000 nt upstream and downstream of the SNP.
To repeat step 4 within a Python script, you can use subprocess.check_output():
...
pts_lbsearch = os.path.join(project_dir, "pts-line-bisect", "pts_lbsearch")
if not os.path.exists(pts_lbsearch):
  print "Content-type:application/json\r\n\r\n"
  print json.dumps({'error' : 'could not locate pts_lbsearch [%s]' % (pts_lbsearch)})
  sys.exit(errno.ENOENT)

if term.startswith('rs'):
  query_fn = os.path.join(annotations_dir, "hg38.snp147.sortedByName.txt")
  cmd = "%s -p %s %s" % (pts_lbsearch, query_fn, term)
  try:
    res = subprocess.check_output(cmd, shell=True)
    if res and len(res) > 0:
      (chr, start, stop) = res.strip().split('\n')[0].split('\t')
    else:
      pass
  except subprocess.CalledProcessError as err:
    sys.stderr.write("Query failed!\n")
    sys.exit(errno.EIO)

#
# map 'chr', 'start', and 'stop' variables to annotations with
# another subprocess call (not shown -- use example above)
#
...

Have you tried using SnpEff and SnpShift (http://snpeff.sourceforge.net) for annotation of the SNPs?
Even if your desired organism available in their annotation database you can create them easily with the gene bank files. Other than annotation the tools features a bunch of other amazing functionalities.

BiocManager::install('biomaRt')
require(biomaRt)
ensembl <- useMart("ENSEMBL_MART_SNP", dataset = "hsapiens_snp")
getBM(attributes=c("refsnp_id", "chr_name", "chrom_start", "chrom_end",
                   "allele", "mapweight", "validated", "allele_1", "minor_allele",
                   "minor_allele_freq", "minor_allele_count", "clinical_significance",
                   "synonym_name", "ensembl_gene_stable_id"),
      filters="snp_filter", values=c("rs6025","rs431905509"),
      mart=ensembl, uniqueRows=TRUE)

using ncbi-efetch and XSLT with the following stylesheet:
<?xml version='1.0' encoding="UTF-8"?>
<xsl:stylesheet
    xmlns:xsl='http://www.w3.org/1999/XSL/Transform'
    xmlns:r="https://www.ncbi.nlm.nih.gov/SNP/docsum"
    version='1.1'
    >
<xsl:output method="text"  encoding="UTF-8"/>

<xsl:template match="/">
<xsl:apply-templates select="r:ExchangeSet"/>
</xsl:template>
<xsl:template match="r:ExchangeSet">
<xsl:apply-templates select="r:Rs"/>
</xsl:template>
<xsl:template match="r:Rs">
<xsl:apply-templates select="r:Assembly/r:Component/r:MapLoc/r:FxnSet"/>
</xsl:template>

<xsl:template match="r:FxnSet">
<xsl:text>rs</xsl:text>
<xsl:value-of select="../../../../@rsId"/>
<xsl:text>  ncbi:</xsl:text>
<xsl:value-of select="@geneId"/>
<xsl:text>  </xsl:text>
<xsl:value-of select="@symbol"/>
<xsl:text>
</xsl:text>
</xsl:template>

</xsl:stylesheet>

usage (group the rs# by 3):
$ echo -e "rs8179065\nrs12899845\nrs8179066\nrs80357332\nrs80356920\n" | xargs  -L 3 echo  | tr " " "," | while read R; do curl -sL "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=snp&id=${R}&retmode=xml" | xsltproc --novalid  transform.xsl -; done | uniq
rs8179065   ncbi:374654 KIF7
rs12899845  ncbi:90381  TICRR
rs8179066   ncbi:374654 KIF7
rs80357332  ncbi:672    BRCA1
rs80356920  ncbi:672    BRCA1


text:
Where can I find a list of SNP rs IDs that belong to the X chromosome?
Where can I find a list of SNP rs IDs that belong to the X chromosome?
I'm trying to exclude the SNPs on X-chromosome in a data set, but I don't know where these SNPs come from.
I tried to use dbSNP to search, it returns 26933249 SNPs, but the webpage seems to refuse me to download the information.

Have you tried the send to option in the top right corner, just below the search bar?

In what species? Your search didn't specify a species, are you really interested in rsIDs from any species?

You just searched for X[Chromosome], you didn't specify a species. Presumably, your dataset comes from a specific species, so you should limit your search to that species. For instance, for human, you can search dbSNP for X[Chromosome] human[organism]. Then, click on the send to link at the top left and download the data:

This is veeeeeery slow, however. I suggest that you instead download the data for chrX from dbSNP's FTP site (this assumes your species of interest is human): ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606/chr_rpts/chr_X.txt.gz
Once you have this file, you can use it to filter your dataset:

Collect all rsIDs from the file:
zcat chr_X.txt.gz | awk '/^[0-9]/{print "rs"$1}' > chrX.rsIDs

Filter your dataset.
You haven't shown the format of your dataset, so if it is just a list of rsIDs like this:
$ cat dataset
rs77669712
rs199706910
rs2799940
rs12205639
rs10958205
rs200778574

You can remove any that are on chrX with:
$ grep -vxFf chrX.rsIDs dataset
rs77669712
rs199706910
rs12205639
rs10958205
rs200778574

That will print all rsIDs from the text file dataset that are not also in the text file chrX.rsIDs. Note the -x option for grep which tells it to only consider matches that cover the entire line. This ensures that rs1234 is not removed when you're looking for rs123.


text:
Is there a way to tell which chromosome a gene is on, by looking at the "Chromosome/scaffold name"
I recently got a data set, from which I need to figure out which chromosome a gene is from, but the head of the data reads like:
Gene ID Description     Gene type       Gene End (bp)   Gene Start (bp) Strand  Associated Gene Name    Chromosome/scaffold name

while one entry of the data is like:
ENSG00000252303 RNA, U6 small nuclear 280, pseudogene [Source:HGNC Symbol;Acc:HGNC:47243]       snRNA   67546754        67546651        1       RNU6-280P       CHR_HG2128_PATCH

It looks like the last column is the information about which chromosome the gene comes from, but how to interpret the "CHR_HG2128_PATCH"?

What you are looking at is a patch, not an actual chromosome location. HG-2128 is an issue ID, something that the Genome Reference Consortium uses to track issue with references. Information regarding that specific patch is here. The issue is described as a "Possible misassembly or indel variation in GRCh38 within AL049860.8."
During the alignment step before you received your data, someone used a top level assembly version of the reference genome. This includes extra, accessioned scaffolds tacked on to the end of the file that contain the fix. These scaffolds represent what will be found in the next major release of the genome. The fixes aren't applied directly to the reference as soon as they are found due to the inevitable disruption of coordinates. So reads mapped to this FIX patch rather than the primary original scaffold in the primary assembly, on which the RNU6-280P gene resides.
The coordinates are slightly different. In the current GRCh38.p12
 assembly the RNU6-280P gene is located at chr6:67546651-67546754 (source). You can see that the coordinates are very close to those indicated in the table you provided (Gene Start, Gene End) but are indeed different.
Here is some more information on patches.