# Question
Title: Kallisto error: index input file could not be opened!
I am utilizing `Kallisto` in `Anaconda`/`miniconda` for RNA sequencing. I have successfully made an index; using said index to analyze RNA sequencing data has yielded an error of:

> "Error: index input file could not be opened!"

The implication is that the code used does in fact work but raises the question as to why the index file could not be opened. I have tried searching this error to see if others may have faced the same error message, but unfortunately no luck. Any recommendations and or solutions would be greatly appreciated. The code and output can be found below.

```
kallisto quant -i Homo_sapiens.GRCh38.cdna.all -o test -t 8 --single -l 250 -s 30 SRR8668755_1M_subsample.fastq

```

> \[quant\] fragment length distribution is truncated gaussian with mean = 250, sd = 30

> Error: index input file could not be opened!

# Answer
Usually errors like these occur when an input file doesn't exist or is corrupt. But you might also see these errors if you don't have permissions to access the file for whatever reason. Although unlikely, these can also occur if/when your system is out of memory.

Assuming the index is corrupt, I think just retry the command. For example:

```
curl -O ftp://ftp.ensembl.org/pub/release-108/fasta/homo_sapiens/cdna/Homo_sapiens.GRCh38.cdna.all.fa.gz

kallisto index \
    -i Homo_sapiens.GRCh38.cdna.all.release-108_k31.idx \
    Homo_sapiens.GRCh38.cdna.all.fa.gz

```

And test:

```
kallisto quant \
    -i Homo_sapiens.GRCh38.cdna.all.release-108_k31.idx \
    -o test \
    -t 8 \
    --single \
    -l 250 \
    -s 30 \
    SRR8668755_1M_subsample.fastq

```

---

To get this working using Windows/Anaconda Powershell, I just downloaded the cDNA sequences from https://ftp.ensembl.org/pub/release-108/fasta/homo\_sapiens/cdna/ and moved them to my working directory. I also downloaded and extracted the Windows Kallisto executable v0.46.1 from https://pachterlab.github.io/kallisto/download.html. Then, with the FASTQ file in my working directory:

```
 .\kallisto_windows-v0.46.1\kallisto\kallisto.exe index `
    -i Homo_sapiens.GRCh38.cdna.all.release-108_k31.idx `
    .\Homo_sapiens.GRCh38.cdna.all.fa.gz

```

And tested with:

```
.\kallisto_windows-v0.46.1\kallisto\kallisto.exe quant `
    -i .\Homo_sapiens.GRCh38.cdna.all.release-108_k31.idx `
    -o test `
    -t 8 `
    --single `
    -l 250 `
    -s 30 `
    .\SRR8668755_1M_subsample.fastq

```
> 5 votes

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


# Question
Title: CNVKit does not output all the accessible regions in the targets bed file
*This question was also asked on Biostars*

I am using CNVkit on my data using hg38 as reference. The command that I am using is the following:

```
cnvkit.py batch sample.bam  -n control.bam -m wgs -f reference.fasta  --target-avg-size 1000 --output-dir results/

```

So, on the CNVkit documentation page is reported that

> The `batch --method wgs` option uses the given reference genome’s sequencing-accessible regions (“access” BED) as the “targets” – these will be calculated on the fly if not provided.

Problem is, when it calculates accesibile regions from the reference, there are a lot of them listed in the log file that are not present in the final .bed output, and I cannot figure out why. Accessible regions are later divided in bins, and bins without enough coverage are discarded later in the process. These are the outputs I get on a short example:

**Log file:**

```
WGS protocol: recommend '--annotate' option (e.g. refFlat.txt) to help locate genes in output files.
NC_000001.10: Scanning for accessible regions
      Accessible region NC_000001.10:10000-177417 (size 167417)
      Accessible region NC_000001.10:227417-267719 (size 40302)
      Accessible region NC_000001.10:317719-471368 (size 153649)
      Accessible region NC_000001.10:521368-2634220 (size 2112852)
      Accessible region NC_000001.10:2684220-3845268 (size 1161048)
      Accessible region NC_000001.10:3995268-13052998 (size 9057730)
      Accessible region NC_000001.10:13102998-13219912 (size 116914)
      Accessible region NC_000001.10:13319912-13557162 (size 237250)
      Accessible region NC_000001.10:13607162-17125658 (size 3518496)
      Accessible region NC_000001.10:17175658-29878082 (size 12702424)
      Accessible region NC_000001.10:30028082-103863906 (size 73835824)
      Accessible region NC_000001.10:103913906-120697156 (size 16783250)
      Accessible region NC_000001.10:120747156-120936695 (size 189539)
      Accessible region NC_000001.10:121086695-121485434 (size 398739)
      Accessible region NC_000001.10:142535434-142731022 (size 195588)
      Accessible region NC_000001.10:142781022-142967761 (size 186739)
      Accessible region NC_000001.10:143117761-143292816 (size 175055)
      Accessible region NC_000001.10:143342816-143544525 (size 201709)
      Accessible region NC_000001.10:143644525-143771002 (size 126477)
      Accessible region NC_000001.10:143871002-144095783 (size 224781)
      Accessible region NC_000001.10:144145783-144224481 (size 78698)
      Accessible region NC_000001.10:144274481-144401744 (size 127263)
      Accessible region NC_000001.10:144451744-144622413 (size 170669)
      Accessible region NC_000001.10:144672413-144710724 (size 38311)
      Accessible region NC_000001.10:144810724-145833118 (size 1022394)
      Accessible region NC_000001.10:145883118-146164650 (size 281532)
      Accessible region NC_000001.10:146214650-146253299 (size 38649)
      Accessible region NC_000001.10:146303299-148026038 (size 1722739)
      Accessible region NC_000001.10:148176038-148361358 (size 185320)
      Accessible region NC_000001.10:148511358-148684147 (size 172789)
      Accessible region NC_000001.10:148734147-148954460 (size 220313)
      Accessible region NC_000001.10:149004460-149459645 (size 455185)
      Accessible region NC_000001.10:149509645-205922707 (size 56413062)
      Accessible region NC_000001.10:206072707-206332221 (size 259514)
      Accessible region NC_000001.10:206482221-223747846 (size 17265625)
      Accessible region NC_000001.10:223797846-235192211 (size 11394365)
      Accessible region NC_000001.10:235242211-248908210 (size 13665999)
      Accessible region NC_000001.10:249058210-249240621 (size 182411)
NT_113878.1: Scanning for accessible regions
      Accessible region NT_113878.1:0-106433 (size 106433)
NT_167207.1: Scanning for accessible regions
      Accessible region NT_167207.1:0-547496 (size 547496)
NC_000002.11: Scanning for accessible regions
      Accessible region NC_000002.11:10000-3529312 (size 3519312)
      Accessible region NC_000002.11:3579312-5018788 (size 1439476)
      Accessible region NC_000002.11:5118788-16279724 (size 11160936)
      Accessible region NC_000002.11:16329724-21153113 (size 4823389)
      Accessible region NC_000002.11:21178113-31705550 (size 10527437)
      Accessible region NC_000002.11:31705551-31725939 (size 20388)
      Accessible region NC_000002.11:31726790-31816827 (size 90037)
      Accessible region NC_000002.11:31816828-31816854 (size 26)
      Accessible region NC_000002.11:31816855-31816858 (size 3)
      Accessible region NC_000002.11:31816859-33092197 (size 1275338)
      Accessible region NC_000002.11:33093197-33141692 (size 48495)
      Accessible region NC_000002.11:33142692-87668206 (size 54525514)
      Accessible region NC_000002.11:87718206-89630436 (size 1912230)
      Accessible region NC_000002.11:89830436-90321525 (size 491089)
      Accessible region NC_000002.11:90371525-90545103 (size 173578)
      Accessible region NC_000002.11:91595103-92326171 (size 731068)
      Accessible region NC_000002.11:95326171-110109337 (size 14783166)
      Accessible region NC_000002.11:110251337-149690582 (size 39439245)
      Accessible region NC_000002.11:149790582-234003741 (size 84213159)
      Accessible region NC_000002.11:234053741-239801978 (size 5748237)
      Accessible region NC_000002.11:239831978-240784132 (size 952154)
      Accessible region NC_000002.11:240809132-243102476 (size 2293344)
      Accessible region NC_000002.11:243152476-243189373 (size 36897)
NC_000003.11: Scanning for accessible regions
      Accessible region NC_000003.11:60000-17137943 (size 17077943)
NT_113878.1: Joining over small gaps
NT_167207.1: Joining over small gaps
Wrote GCF_000001405.25_GRCh37.p13_genomic.bed with 2 regions

```

So, multiple accessible regions are present for sequences NC\_000001.10 , NC\_000002.11 and NC\_000003.11, but only regions for NT\_113878.1 and NT\_167207.1 are reported:

```
NT_113878.1     0       106433
NT_167207.1     0       547496

```

Why is that so? I cannot really figure out.

I am sure the missed regions have data, because I am also using another CNVs caller on the same exact data, and it shows high coverage for the same regions CNVkit is discarding. Also the coverage for control data is high in those discarder regions, so this shouldn't be the matter, but I have no ideas at this point.

**EDIT:** The question was answered on Biostars. Looking into the code, CNVkit is skipping chromosomes which don't have canonical names, like 'NCs' (so, the fact that it outputs the 'NTs' chromosomes is kind of incorrect, because they're also non-canonical and they should be skipped).

# Answer
The OP received an answer on their question on Biostars:

> The question was answered on Biostars. Looking into the code, CNVkit is skipping chromosomes which don't have canonical names, like 'NCs' (so, the fact that it outputs the 'NTs' chromosomes is kind of incorrect, because they're also non-canonical and they should be skipped).
> 0 votes

---
Tags: reference-genome, bed, wgs, cnv
---


title: Error installing VEP on mac

text:
# Question
Title: Error installing VEP on mac
I have tried multiple different ways to get Ensemble's Variant Effect Predictor (`VEP`) to properly install on my mac, but I keep getting same error.

```
In /Users/username/ensembl-vep/htslib
make: Nothing to be done for `all'.
 - unpacking ./Bio/tmp/biodbhts.zip to ./Bio/tmp/
./Bio/tmp/Bio-DB-HTS-2.11 - moving files to ./biodbhts
 - making Bio::DB:HTS
clang-4.0: error: invalid version number in '-mmacosx-version-min=12.4'
Warning: ExtUtils::CBuilder not installed or no compiler detected
Proceeding with configuration, but compilation may fail during Build

Created MYMETA.yml and MYMETA.json
Creating new 'Build' script for 'Bio-DB-HTS' version '2.11'
Building Bio-DB-HTS
Error: no compiler detected to compile 'lib/Bio/DB/HTS.c'.  Aborting
ERROR: Shared Bio::DB:HTS library not found

```

--\> When I check for the version of clang, I always get version 12, not sure where clang-4.0 is. I think error is specifically not correct path to current version of `C++`.

# Answer
> 2 votes
You can try and get Perl's ExUtils::Builder, described here correctly installed and working, for example

```
perl -MCPAN -e shell
force install Bundle::CPAN

```

You might want to use `brew` to make the build safe.

I don't do this stuff, I use OSX, I would install ensemble's VEP via `conda`

```
conda install -c bioconda ensembl-vep

```

You might want to (good idea) create a specific environment and obviously download `Anaconda` or `miniconda`.

```
conda create --name vep
conda activate vep

```

It will work, the output is below. I didn't proceed with the installation because of the shear amount of stuff downloaded into a given conda environment. There's everything in there but the kitchen sink, like there is `mafft`, `t-coffee`, `clustalw`, `muscle` and `clustalo` required in the download. I am honestly clueless because they are all alignment algorithms. Anyway I honestly wouldn't chase the MCPAN stuff, in this instance I definitely recommend `conda`. Finally, you honestly you don't want all this extra stuff on your hard drive when you've finished just delete the entire environment.

The final point I'd like to make is `muscle5` (which is also in there) has recently significantly been modified making all other aligners obsolete on the runtime issue, but has to be run under the `super5` mode.

```
## Package Plan ##

  environment location: /opt/anaconda3/envs/XXX # environment was edited  out

  added / updated specs:
    - ensembl-vep

The following packages will be downloaded:

    package                    |            build
    ---------------------------|-----------------
    argtable2-2.13             |    h1de35cc_1001         2.7 MB  conda-forge
    clustalo-1.2.4             |       h9722bc1_5         331 KB  bioconda
    clustalw-2.1               |       hcd10b59_7         326 KB  bioconda
    ensembl-vep-108.2          | pl5321h4a94de4_0         3.1 MB  bioconda
    expat-2.4.9                |       he9d5cce_0         127 KB
    htslib-1.14                |       hc38c3fb_0         2.0 MB  bioconda
    libdeflate-1.7             |       h9ed2024_5          52 KB
    libgfortran-5.0.0          |11_3_0_hecd8cb5_28         142 KB
    libgfortran5-11.3.0        |      h9dfd629_28         1.4 MB
    libzlib-1.2.12             |       hfd90126_3          65 KB  conda-forge
    mafft-7.508                |       ha5712d3_0         8.7 MB  bioconda
    muscle-5.1                 |       hb339e23_1         298 KB  bioconda
    mysql-connector-c-6.1.11   |       hccea1a4_1         1.2 MB
    paml-4.9                   |       hb4d813b_6         1.2 MB  bioconda
    perl-algorithm-diff-1.1903 | pl5321hdfd78af_3          36 KB  bioconda
    perl-base-2.23             | pl5321hdfd78af_2          12 KB  bioconda
    perl-bio-asn1-entrezgene-1.73| pl5321hdfd78af_3          24 KB  bioconda
    perl-bio-coordinate-1.007001| pl5321hdfd78af_3          36 KB  bioconda
    perl-bio-db-hts-3.01       | pl5321hc7384eb_4         116 KB  bioconda
    perl-bio-featureio-1.6.905 | pl5321hdfd78af_4          41 KB  bioconda
    perl-bio-samtools-1.43     | pl5321h1f540d2_3         161 KB  bioconda
    perl-bio-tools-phylo-paml-1.7.3| pl5321hdfd78af_3          58 KB  bioconda
    perl-bio-tools-run-alignment-clustalw-1.7.4| pl5321hdfd78af_3          25 KB  bioconda
    perl-bio-tools-run-alignment-tcoffee-1.7.4| pl5321hdfd78af_5          25 KB  bioconda
    perl-bioperl-1.7.8         |       hdfd78af_0           8 KB  bioconda
    perl-bioperl-core-1.7.8    | pl5321h527b516_0         1.8 MB  bioconda
    perl-bioperl-run-1.007003  | pl5321hdfd78af_0         340 KB  bioconda
    perl-business-isbn-3.007   | pl5321hdfd78af_0          18 KB  bioconda
    perl-business-isbn-data-20210112.006| pl5321hdfd78af_0          21 KB  bioconda
    perl-capture-tiny-0.48     | pl5321hdfd78af_2          19 KB  bioconda
    perl-class-data-inheritable-0.09| pl5321hdfd78af_0          10 KB  bioconda
    perl-constant-1.33         | pl5321hdfd78af_2          12 KB  bioconda
    perl-data-dumper-2.183     | pl5321ha5712d3_1          38 KB  bioconda
    perl-dbd-mysql-4.046       | pl5321hcd10b59_4          94 KB  bioconda
    perl-dbi-1.643             | pl5321ha5712d3_1         602 KB  bioconda
    perl-devel-stacktrace-2.04 | pl5321hdfd78af_1          17 KB  bioconda
    perl-digest-hmac-1.04      | pl5321hdfd78af_0          10 KB  bioconda
    perl-digest-md5-2.58       | pl5321ha5712d3_1          18 KB  bioconda
    perl-encode-locale-1.05    | pl5321hdfd78af_7          13 KB  bioconda
    perl-exception-class-1.45  | pl5321hdfd78af_0          21 KB  bioconda
    perl-file-listing-6.15     | pl5321hdfd78af_0          14 KB  bioconda
    perl-file-slurp-tiny-0.004 | pl5321hdfd78af_2          10 KB  bioconda
    perl-file-sort-1.01        | pl5321hdfd78af_3          20 KB  bioconda
    perl-file-spec-3.48_01     | pl5321hdfd78af_2           6 KB  bioconda
    perl-getopt-long-2.54      | pl5321hdfd78af_0          31 KB  bioconda
    perl-html-parser-3.80      | pl5321hb3ea9af_0          67 KB  bioconda
    perl-html-tagset-3.20      | pl5321hdfd78af_4          14 KB  bioconda
    perl-http-cookies-6.10     | pl5321hdfd78af_0          22 KB  bioconda
    perl-http-daemon-6.14      | pl5321hdfd78af_0          20 KB  bioconda
    perl-http-date-6.05        | pl5321hdfd78af_0          13 KB  bioconda
    perl-http-message-6.36     | pl5321hdfd78af_0          54 KB  bioconda
    perl-http-negotiate-6.01   | pl5321hdfd78af_4          14 KB  bioconda
    perl-io-html-1.004         | pl5321hdfd78af_0          16 KB  bioconda
    perl-io-socket-ssl-2.074   | pl5321hdfd78af_0         162 KB  bioconda
    perl-io-string-1.08        | pl5321hdfd78af_4          12 KB  bioconda
    perl-io-tty-1.16           | pl5321ha5712d3_1          30 KB  bioconda
    perl-ipc-run-20200505.0    | pl5321hdfd78af_0          76 KB  bioconda
    perl-libwww-perl-6.67      | pl5321hdfd78af_0          97 KB  bioconda
    perl-libxml-perl-0.08      | pl5321hdfd78af_3          40 KB  bioconda
    perl-lwp-mediatypes-6.04   | pl5321hdfd78af_1          24 KB  bioconda
    perl-mime-base64-3.16      | pl5321ha5712d3_2          16 KB  bioconda
    perl-net-http-6.22         | pl5321hdfd78af_0          21 KB  bioconda
    perl-net-ssleay-1.92       | pl5321h5d619eb_1         299 KB  bioconda
    perl-ntlm-1.09             | pl5321hdfd78af_5          16 KB  bioconda
    perl-perlio-gzip-0.20      | pl5321h1f540d2_3          16 KB  bioconda
    perl-sereal-4.019          | pl5321hdfd78af_0          11 KB  bioconda
    perl-sereal-decoder-4.025  | pl5321ha5712d3_1         156 KB  bioconda
    perl-sereal-encoder-4.025  | pl5321ha5712d3_1         365 KB  bioconda
    perl-set-intervaltree-0.12 | pl5321h9722bc1_2          31 KB  bioconda
    perl-socket-2.027          | pl5321ha5712d3_3          29 KB  bioconda
    perl-sub-uplevel-0.2800    | pl5321ha5712d3_4          14 KB  bioconda
    perl-test-deep-1.130       | pl5321hdfd78af_0          37 KB  bioconda
    perl-test-differences-0.69 | pl5321hdfd78af_0          19 KB  bioconda
    perl-test-exception-0.43   | pl5321hdfd78af_3          15 KB  bioconda
    perl-test-harness-3.44     | pl5321hdfd78af_0          85 KB  bioconda
    perl-test-most-0.38        | pl5321hdfd78af_0          20 KB  bioconda
    perl-test-warn-0.36        | pl5321hdfd78af_2          14 KB  bioconda
    perl-text-csv-2.01         | pl5321hdfd78af_0          70 KB  bioconda
    perl-text-diff-1.45        | pl5321hdfd78af_1          23 KB  bioconda
    perl-time-local-1.30       | pl5321hdfd78af_0          12 KB  bioconda
    perl-timedate-2.33         | pl5321hdfd78af_2          31 KB  bioconda
    perl-tree-dag_node-1.32    | pl5321hdfd78af_0          43 KB  bioconda
    perl-try-tiny-0.31         | pl5321hdfd78af_1          19 KB  bioconda
    perl-uri-5.12              | pl5321hdfd78af_0          59 KB  bioconda
    perl-url-encode-0.03       | pl5321h527b516_0          12 KB  bioconda
    perl-www-robotrules-6.02   | pl5321hdfd78af_4          14 KB  bioconda
    perl-xml-dom-1.46          | pl5321hdfd78af_1          59 KB  bioconda
    perl-xml-dom-xpath-0.14    | pl5321hdfd78af_2          12 KB  bioconda
    perl-xml-parser-2.44       | pl5321h50c80e2_8         162 KB  bioconda
    perl-xml-regexp-0.04       | pl5321hdfd78af_3          10 KB  bioconda
    perl-xml-xpathengine-0.14  | pl5321hdfd78af_3          28 KB  bioconda
    t-coffee-11.00.8cbe486     |       h26a2512_0         1.1 MB  bioconda
    unzip-6.0                  |       h3b3adb9_0          71 KB
    ------------------------------------------------------------
                                           Total:        28.7 MB

The following NEW packages will be INSTALLED:

  argtable2          conda-forge/osx-64::argtable2-2.13-h1de35cc_1001
  blast              bioconda/osx-64::blast-2.12.0-h0370960_3
  clustalo           bioconda/osx-64::clustalo-1.2.4-h9722bc1_5
  clustalw           bioconda/osx-64::clustalw-2.1-hcd10b59_7
  ensembl-vep        bioconda/noarch::ensembl-vep-108.2-pl5321h4a94de4_0
  entrez-direct      bioconda/osx-64::entrez-direct-16.2-h193322a_1
  expat              pkgs/main/osx-64::expat-2.4.9-he9d5cce_0
  htslib             bioconda/osx-64::htslib-1.14-hc38c3fb_0
  libdeflate         pkgs/main/osx-64::libdeflate-1.7-h9ed2024_5
  libgfortran        pkgs/main/osx-64::libgfortran-5.0.0-11_3_0_hecd8cb5_28
  libgfortran5       pkgs/main/osx-64::libgfortran5-11.3.0-h9dfd629_28
  libidn2            pkgs/main/osx-64::libidn2-2.3.2-h9ed2024_0
  libunistring       pkgs/main/osx-64::libunistring-0.9.10-h9ed2024_0
  libzlib            conda-forge/osx-64::libzlib-1.2.12-hfd90126_3
  llvm-openmp        pkgs/main/osx-64::llvm-openmp-14.0.6-h0dcd299_0
  mafft              bioconda/osx-64::mafft-7.508-ha5712d3_0
  muscle             bioconda/osx-64::muscle-5.1-hb339e23_1
  mysql-connector-c  pkgs/main/osx-64::mysql-connector-c-6.1.11-hccea1a4_1
  paml               bioconda/osx-64::paml-4.9-hb4d813b_6
  pcre               pkgs/main/osx-64::pcre-8.45-h23ab428_0
  perl               pkgs/main/osx-64::perl-5.32.1-0_h435f0c2_perl5
  perl-algorithm-di~ bioconda/noarch::perl-algorithm-diff-1.1903-pl5321hdfd78af_3
  perl-archive-tar   bioconda/noarch::perl-archive-tar-2.40-pl5321hdfd78af_0
  perl-base          bioconda/noarch::perl-base-2.23-pl5321hdfd78af_2
  perl-bio-asn1-ent~ bioconda/noarch::perl-bio-asn1-entrezgene-1.73-pl5321hdfd78af_3
  perl-bio-coordina~ bioconda/noarch::perl-bio-coordinate-1.007001-pl5321hdfd78af_3
  perl-bio-db-hts    bioconda/osx-64::perl-bio-db-hts-3.01-pl5321hc7384eb_4
  perl-bio-featureio bioconda/noarch::perl-bio-featureio-1.6.905-pl5321hdfd78af_4
  perl-bio-samtools  bioconda/osx-64::perl-bio-samtools-1.43-pl5321h1f540d2_3
  perl-bio-tools-ph~ bioconda/noarch::perl-bio-tools-phylo-paml-1.7.3-pl5321hdfd78af_3
  perl-bio-tools-ru~ bioconda/noarch::perl-bio-tools-run-alignment-clustalw-1.7.4-pl5321hdfd78af_3
  perl-bio-tools-ru~ bioconda/noarch::perl-bio-tools-run-alignment-tcoffee-1.7.4-pl5321hdfd78af_5
  perl-bioperl       bioconda/noarch::perl-bioperl-1.7.8-hdfd78af_0
  perl-bioperl-core  bioconda/osx-64::perl-bioperl-core-1.7.8-pl5321h527b516_0
  perl-bioperl-run   bioconda/noarch::perl-bioperl-run-1.007003-pl5321hdfd78af_0
  perl-business-isbn bioconda/noarch::perl-business-isbn-3.007-pl5321hdfd78af_0
  perl-business-isb~ bioconda/noarch::perl-business-isbn-data-20210112.006-pl5321hdfd78af_0
  perl-capture-tiny  bioconda/noarch::perl-capture-tiny-0.48-pl5321hdfd78af_2
  perl-carp          bioconda/noarch::perl-carp-1.38-pl5321hdfd78af_4
  perl-class-data-i~ bioconda/noarch::perl-class-data-inheritable-0.09-pl5321hdfd78af_0
  perl-common-sense  bioconda/noarch::perl-common-sense-3.75-pl5321hdfd78af_0
  perl-compress-raw~ bioconda/osx-64::perl-compress-raw-bzip2-2.201-pl5321h9722bc1_0
  perl-compress-raw~ bioconda/osx-64::perl-compress-raw-zlib-2.105-pl5321h9722bc1_0
  perl-constant      bioconda/noarch::perl-constant-1.33-pl5321hdfd78af_2
  perl-data-dumper   bioconda/osx-64::perl-data-dumper-2.183-pl5321ha5712d3_1
  perl-dbd-mysql     bioconda/osx-64::perl-dbd-mysql-4.046-pl5321hcd10b59_4
  perl-dbi           bioconda/osx-64::perl-dbi-1.643-pl5321ha5712d3_1
  perl-devel-stackt~ bioconda/noarch::perl-devel-stacktrace-2.04-pl5321hdfd78af_1
  perl-digest-hmac   bioconda/noarch::perl-digest-hmac-1.04-pl5321hdfd78af_0
  perl-digest-md5    bioconda/osx-64::perl-digest-md5-2.58-pl5321ha5712d3_1
  perl-encode        bioconda/osx-64::perl-encode-3.19-pl5321ha5712d3_1
  perl-encode-locale bioconda/noarch::perl-encode-locale-1.05-pl5321hdfd78af_7
  perl-exception-cl~ bioconda/noarch::perl-exception-class-1.45-pl5321hdfd78af_0
  perl-exporter      bioconda/noarch::perl-exporter-5.72-pl5321hdfd78af_2
  perl-exporter-tiny bioconda/noarch::perl-exporter-tiny-1.002002-pl5321hdfd78af_0
  perl-extutils-mak~ conda-forge/noarch::perl-extutils-makemaker-7.64-pl5321hd8ed1ab_0
  perl-file-listing  bioconda/noarch::perl-file-listing-6.15-pl5321hdfd78af_0
  perl-file-slurp-t~ bioconda/noarch::perl-file-slurp-tiny-0.004-pl5321hdfd78af_2
  perl-file-sort     bioconda/noarch::perl-file-sort-1.01-pl5321hdfd78af_3
  perl-file-spec     bioconda/noarch::perl-file-spec-3.48_01-pl5321hdfd78af_2
  perl-getopt-long   bioconda/noarch::perl-getopt-long-2.54-pl5321hdfd78af_0
  perl-html-parser   bioconda/osx-64::perl-html-parser-3.80-pl5321hb3ea9af_0
  perl-html-tagset   bioconda/noarch::perl-html-tagset-3.20-pl5321hdfd78af_4
  perl-http-cookies  bioconda/noarch::perl-http-cookies-6.10-pl5321hdfd78af_0
  perl-http-daemon   bioconda/noarch::perl-http-daemon-6.14-pl5321hdfd78af_0
  perl-http-date     bioconda/noarch::perl-http-date-6.05-pl5321hdfd78af_0
  perl-http-message  bioconda/noarch::perl-http-message-6.36-pl5321hdfd78af_0
  perl-http-negotia~ bioconda/noarch::perl-http-negotiate-6.01-pl5321hdfd78af_4
  perl-io-compress   bioconda/osx-64::perl-io-compress-2.201-pl5321h9722bc1_0
  perl-io-html       bioconda/noarch::perl-io-html-1.004-pl5321hdfd78af_0
  perl-io-socket-ssl bioconda/noarch::perl-io-socket-ssl-2.074-pl5321hdfd78af_0
  perl-io-string     bioconda/noarch::perl-io-string-1.08-pl5321hdfd78af_4
  perl-io-tty        bioconda/osx-64::perl-io-tty-1.16-pl5321ha5712d3_1
  perl-io-zlib       bioconda/noarch::perl-io-zlib-1.11-pl5321hdfd78af_0
  perl-ipc-run       bioconda/noarch::perl-ipc-run-20200505.0-pl5321hdfd78af_0
  perl-json          bioconda/noarch::perl-json-4.10-pl5321hdfd78af_0
  perl-json-xs       bioconda/osx-64::perl-json-xs-2.34-pl5321hcd10b59_5
  perl-libwww-perl   bioconda/noarch::perl-libwww-perl-6.67-pl5321hdfd78af_0
  perl-libxml-perl   bioconda/noarch::perl-libxml-perl-0.08-pl5321hdfd78af_3
  perl-list-moreuti~ bioconda/noarch::perl-list-moreutils-0.430-pl5321hdfd78af_0
  perl-list-moreuti~ bioconda/osx-64::perl-list-moreutils-xs-0.430-pl5321ha5712d3_1
  perl-lwp-mediatyp~ bioconda/noarch::perl-lwp-mediatypes-6.04-pl5321hdfd78af_1
  perl-mime-base64   bioconda/osx-64::perl-mime-base64-3.16-pl5321ha5712d3_2
  perl-net-http      bioconda/noarch::perl-net-http-6.22-pl5321hdfd78af_0
  perl-net-ssleay    bioconda/osx-64::perl-net-ssleay-1.92-pl5321h5d619eb_1
  perl-ntlm          bioconda/noarch::perl-ntlm-1.09-pl5321hdfd78af_5
  perl-parent        bioconda/noarch::perl-parent-0.236-pl5321hdfd78af_2
  perl-pathtools     bioconda/osx-64::perl-pathtools-3.75-pl5321ha5712d3_3
  perl-perlio-gzip   bioconda/osx-64::perl-perlio-gzip-0.20-pl5321h1f540d2_3
  perl-scalar-list-~ bioconda/osx-64::perl-scalar-list-utils-1.62-pl5321ha5712d3_0
  perl-sereal        bioconda/noarch::perl-sereal-4.019-pl5321hdfd78af_0
  perl-sereal-decod~ bioconda/osx-64::perl-sereal-decoder-4.025-pl5321ha5712d3_1
  perl-sereal-encod~ bioconda/osx-64::perl-sereal-encoder-4.025-pl5321ha5712d3_1
  perl-set-interval~ bioconda/osx-64::perl-set-intervaltree-0.12-pl5321h9722bc1_2
  perl-socket        bioconda/osx-64::perl-socket-2.027-pl5321ha5712d3_3
  perl-sub-uplevel   bioconda/osx-64::perl-sub-uplevel-0.2800-pl5321ha5712d3_4
  perl-test-deep     bioconda/noarch::perl-test-deep-1.130-pl5321hdfd78af_0
  perl-test-differe~ bioconda/noarch::perl-test-differences-0.69-pl5321hdfd78af_0
  perl-test-excepti~ bioconda/noarch::perl-test-exception-0.43-pl5321hdfd78af_3
  perl-test-harness  bioconda/noarch::perl-test-harness-3.44-pl5321hdfd78af_0
  perl-test-most     bioconda/noarch::perl-test-most-0.38-pl5321hdfd78af_0
  perl-test-warn     bioconda/noarch::perl-test-warn-0.36-pl5321hdfd78af_2
  perl-text-csv      bioconda/noarch::perl-text-csv-2.01-pl5321hdfd78af_0
  perl-text-diff     bioconda/noarch::perl-text-diff-1.45-pl5321hdfd78af_1
  perl-time-local    bioconda/noarch::perl-time-local-1.30-pl5321hdfd78af_0
  perl-timedate      bioconda/noarch::perl-timedate-2.33-pl5321hdfd78af_2
  perl-tree-dag_node bioconda/noarch::perl-tree-dag_node-1.32-pl5321hdfd78af_0
  perl-try-tiny      bioconda/noarch::perl-try-tiny-0.31-pl5321hdfd78af_1
  perl-types-serial~ bioconda/noarch::perl-types-serialiser-1.01-pl5321hdfd78af_0
  perl-uri           bioconda/noarch::perl-uri-5.12-pl5321hdfd78af_0
  perl-url-encode    bioconda/osx-64::perl-url-encode-0.03-pl5321h527b516_0
  perl-www-robotrul~ bioconda/noarch::perl-www-robotrules-6.02-pl5321hdfd78af_4
  perl-xml-dom       bioconda/noarch::perl-xml-dom-1.46-pl5321hdfd78af_1
  perl-xml-dom-xpath bioconda/noarch::perl-xml-dom-xpath-0.14-pl5321hdfd78af_2
  perl-xml-parser    bioconda/osx-64::perl-xml-parser-2.44-pl5321h50c80e2_8
  perl-xml-regexp    bioconda/noarch::perl-xml-regexp-0.04-pl5321hdfd78af_3
  perl-xml-xpatheng~ bioconda/noarch::perl-xml-xpathengine-0.14-pl5321hdfd78af_3
  t-coffee           bioconda/osx-64::t-coffee-11.00.8cbe486-h26a2512_0
  unzip              pkgs/main/osx-64::unzip-6.0-h3b3adb9_0
  wget               pkgs/main/osx-64::wget-1.21.3-h6dfd666_0

```

# Answer
> 1 votes
VEP requires `gcc` and `g++` for installation. Before running the build script, try setting the `CC` and `CXX` environment variables to the paths of your `gcc` and `g++` compilers, respectively:

```
export CC=gcc
export CXX=g++

```

You'll need to install these compilers if they don't already exist on your system.

---
Tags: ensembl, software-installation, vep
---


# Question
Title: Ubiquitous regulation of highly specific marker genes
I am fairly new to scRNAseq analysis and keep running into the same problem in the two datasets I am currently working with. We work with kidney inflammation and have two new mouse models, which we are now analysing in scRNAseq. I followed the usual Seurat workflow and ended up getting some nice UMAPs. With `FindAllMarkers` I was able to identify the clusters pretty well.

But once I do differential expression analysis (either with `FindMarkers` or `pseudobulking`) I got weird results.

**Background** There is a highly specialized epithelial celltype in the kidney (the podocyte) that specifically expresses a host of markers (like e.g. *Nphs1* and *Nphs2*). They get downregulated when the kidney gets damaged (as it does in our disease models).

**Problem**

1. I do not expect to see those markers to be present or regulated in any other celltype. But when I do the **DE**, those podocyte-specific marker genes pop up in almost all other clusters and are described to be highly regulated (downregulated, as I would expect in the pdoocytes).

* So in my eyes, this points to "doublets of podocytes" and other celltypes being involved here. Maybe some damaged podocytes got into the same droplet as other celltypes and thus distort the analysis.

2. I then used somewhat stringent filtering via `scDblFinder`* to remove a lot of doublets (~ 10 % of my dataset).

* However, I still get the same problem.

3. Then I analysed each batch (we have 3 in total) individually, but in each one this problem comes up.

Overall, I see 6.4% of doublets (that is 8,710 doublets in a total of 135,844 cells). I have two podocyte clusters with 1.7% doublets (that is 462 doublets in a total of 26,693 podocytes).

Some help / advice would be highly appreciated as I am becoming a little frustrated with this issue.

Thanks a lot!

\*, Bioconductor package

# Answer
***An explanation in 3 parts***

**Final Edit:** It turns out there was contamination in the data which you can read about in the comments, without being able to see it myself I could not know this. If you see this in an experiment and you're clustering from totally raw unfiltered data this maybe the case for you. Nevertheless as you subcluster within cell types of particular interest to get more refined marker genes the advice I provide below may be of use.

**EDIT:** The note above about looking at the violin plots may be useful, but again what your downstream analysis or experiment should impact your decision making.

**TLDR:** The presence of marker genes in other cell types is not surprising nor necessarily a problem. A marker gene that seurat finds are highly expressed in a cell type relative to the average in all cells not a gene that is uniquely expressed, if you are looking for uniquely expressed genes you will probably not find them.

**1) A note on seurat**

So most genes are expressed in most cells there is nothing wrong with your data or your analysis. Gene expression is very messy and the term "marker gene" is a bit of a misnomer. Note that other analysis pipelines like Monocle have different methods for determining what constitutes a marker gene.

Please consider this meta-analysis of gabanergic cerebral cortex neurons (i.e. inhibitory cortical neurons, interneurons etc etc).
https://ieeexplore.ieee.org/abstract/document/9669888

The fact is that most cell populations express a gradient of genes and that each gene likely has multiple roles in each cell. If your goal is to simply publish on expressed marker genes taking the top highest uniquely expressed genes which achieve significance in each cluster will be sufficient, just be sure to report ALL your settings and provide the meta data when you publish. If you are going on to do MERFISH or stain with the results you may want to be more selective.

FindaAllMarkers() uses a simple DEG analysis where it takes the average expression of each feature (the genes you selected, which I hope you did because it's not wroth looking at every gene for clustering just the top most variable) in the cluster against the average expression for these genes in all other clusters. Note that this means you are doing A TON of statistical tests so they use the bonferroni correction which would normally be considered to conservative but in this case is probably still too lenient. You should only be reporting on the very top 5 maybe 10 expressed genes. If your lucky one of these will be noticeably more prominent than the others.

**2) Alternatives to seurat** Monocle package has a different definition which arguably is better: http://cole-trapnell-lab.github.io/monocle-release/docs/#classifying-and-counting-cells though this package is less widely used. The monocle score accounts for how widely expressed across all clusters a marker is and the more unique to a cluster a gene is the greater its score in addition to its actual expression strength.
There is also a package called pandora, but I know next to nothing about it. Read up on them and decide which has features you like best.

Lastly I think part of your problem is just wanting biology to be a cleaner science than it is. The reality is that things we call cell types are simply transient states a population of cells passes through during development and various environmental perturbations, the reality is each cell is a dynamic system and you are sampling a portion of these graded changes. Some states are more stable than others and we have taken to calling these cell types but remember you can always subcluster within a cluster and pull out more clusters which can themselves be subclustered and so on and so forth. It calls into question the definition of a cell type even more and to be honest at this point finding cell types is somewhat of an art, the science is yet to be settled.

I encourage you to do more reading on gene expression within single cells and the definition of a cell type (or if such a thing even exists). Here is some reading to get you started.

Cembrowski, Mark S., and Vilas Menon. "Continuous variation within cell types of the nervous system." Trends in Neurosciences 41, no. 6 (2018): 337-348. https://www.sciencedirect.com/science/article/pii/S0166223618300596

Muñoz-Manchado, A. B., Gonzales, C. B., Zeisel, A., Munguba, H., Bekkouche, B., Skene, N. G., ... & Hjerling-Leffler, J. (2018). Diversity of interneurons in the dorsal striatum revealed by single-cell RNA sequencing and PatchSeq. Cell reports, 24(8), 2179-2190. https://www.sciencedirect.com/science/article/pii/S2211124718311562

If you really want to lean into this you should start digging into the concept of meta-cells, there is a package called Metacell2 which may be of use, be sure to read the documentation and methods section for Metacell (the first version) as Metacell2 won't make much sense without it:

Harris, Benjamin D., Megan Crow, Stephan Fischer, and Jesse Gillis. "Single-cell co-expression analysis reveals that transcriptional modules are shared across cell types in the brain." Cell Systems 12, no. 7 (2021): 748-756. -this is not metacell but is a meta analysis of cell markers using metacells with some pretty surprising and counter intuitive results about marker genes

Baran, Y., Bercovich, A., Sebe-Pedros, A., Lubling, Y., Giladi, A., Chomsky, E., ... & Tanay, A. (2019). MetaCell: analysis of single-cell RNA-seq data using K-nn graph partitions. Genome biology, 20(1), 1-19.

*Ensure you read all documentation carefully and understand the stats they are using. The methods use are a specific definition. Also a real marker gene will be robust to multiple methods of identifying them, if they disappear based on a change in normalization or in other packages you may have your answer there. Bear in mind if you publish your results should be reproducible and people will need to know exactly how you got your results so share your code and metadata please.*

**3) What do you want to get from this analysis?**

As I mentioned earlier if your goal is to get published, simply finding some markers in seurat may suffice. If your goal is to do further experiments based on these results consider your hypothesis carefully. If you want more certainty pull some similar kidney scRNA seq data from GEO and compare your results to their or integrate your data with theirs seurat has a tutorial you can follow along with to do this https://satijalab.org/seurat/articles/integration\_introduction.html#performing-integration-on-datasets-normalized-with-sctransform-1 Note you can also use this pipeline to do batch correction, that may help you. The sct transform give pretty good results.

Sorry nature isn't like a biology textbook =(, I'm not saying that to be facetious, I just want to remind you that living things are very complex. Lean into that complexity and make some interesting discoveries or just accept it and dig out some results you can publish. Talk with your supervisor and explain the issue. Hope this helps!

**4) I lied sorry this is a response in 4 parts, make sure you report exactly what you did and share your metadata and pipeline**

Did I mention your results should be reproducible? How can one do this? Here are some tips. Make sure your keeping your analysis on github, so you can keep track of all the changes it goes through and revert to previous versions if you change things. Keep all your code, its worth having a file that is full of old code you threw away cause it was messy or you didn't end up using it. Not broken code but stuff you don't think is useful. Report your final results, share your data and make sure a monkey could follow you instructions to reproduce your results as well as being able to double check they get the same thing. Use a conda environment to ensure software updates will not destroy your code for future lab members or people who read your future paper. Save all file paths as variables at the top of your code in your scripts so people can easily change the file structure without having to go through your code and struggle finding which read.csv() functions they need to edit.
> 2 votes

# Answer
The violin plots show only little expression in non-podocyte clusters, but the podocyte marker transcripts are only the most obvious manifestation of this problem. Also, small but significant amounts of endothelial and other celltype transcripts show up in most of the other clusters. Identification of the celltype works fine but these contaminating counts throw off most of the downstream analyses.

By now, we did manage to find a solution for the problem, though: It turns out the problem was due to high amounts of ambient RNA contamination. I tried three different packages to remove the contaminating counts (DIEM, DecontX and SoupX) and found that DecontX (and SoupX) worked best for my dataset.
> 2 votes

# Answer
You should make violin plots of the genes of interest (x-axis is your different clusters (i.e. cell-types) and y-axis is expression). Maybe they are differentially expressed but at the same time not expressed very highly.

Also why does their expression matter in other cell types? Do you have strong evidence to expect that these genes would only be expressed in podocytes?
> 1 votes

---
Tags: r, rna-seq, scrnaseq, seurat
---


# Question
Title: Need to convert .bed to .vcf. Can the reference build (and necessary .fasta) be determined for a .bed file?
I have some .bed files (and .bim and .fam files) containing data for a number of different samples, and I need to convert them to .vcf.

I've found `bed2vcf` which is from the `bedr` R package, and also available as a stand-alone python program of the same name, but it requires a .fasta file corresponding to the reference genome used to create the .bed file. So we need to know what genome reference build was used to get the right `.fasta`.

**Question** Is there any tool or way to inspect the files I have and determine what genome reference build is required to do a correct conversion to .vcf, or is there any other method?

There is a `bed2vcf` related post here.

# Answer
There are two different kinds of `.bed` files in bioinformatics; one related to storing genomic coordinate data and one which is part of the plink binary dataset, so it's important to check which kind of file you are actually converting. If you have `.bed`, `.bim` and `.fam` files, you could use PLINK 2.0 to export them as VCF:

```
plink2 --bfile [prefix] --export vcf --out <prefix>

```

The relevant options from the docs:

```
<...>

  --bfile  <prefix> ['vzs'] : Specify .bed + .bim[.zst] + .fam prefix

<...>

  --out <prefix>     : Specify prefix for output files.

<...>

  --export <output format(s)...> [{01 | 12}] ['bgz'] ['id-delim='<char>]
           ['id-paste='<column set descriptor>] ['include-alt']
           ['omit-nonmale-y'] ['spaces'] ['vcf-dosage='<field>] ['ref-first']
           ['bits='<#>] ['sample-v2']
    Create a new fileset with all filters applied.  The following output
    formats are supported:

...

    * 'vcf',     : VCF (default version 4.3).  If PAR1 and PAR2 are present,
      'vcf-4.2',   they are automatically merged with chrX, with proper
      'bcf',       handling of chromosome codes and male ploidy.
      'bcf-4.2'    When the 'bgz' modifier is present, the VCF file is
                   block-gzipped.  (This always happens with BCF output.)
                   The 'id-paste' modifier controls which .psam columns are
                   used to construct sample IDs (choices are maybefid, fid,
                   iid, maybesid, and sid; default is maybefid,iid,maybesid),
                   while the 'id-delim' modifier sets the character between the
                   ID pieces (default '_').

...

```
> 2 votes

---
Tags: ngs, fasta, vcf, bed
---


# Question
Title: split fast reference fast file by 22 chromosome
HG38 human reference is the collection of reference from all 22 chromosomes,x,y chromosomes, and some others. I was targeting to split the reference file into 22 chromosomes fasta files. I am following this link https://crashcourse.housegordon.org/split-fasta-files.html to split the fasta files.More specifically I want to split the fasta file into all possible (placed,unlocalized and unplaced) sequence and the merge placed and localized sequences. However, I am unsure where to put the HLA part of the reference. Should I put it in Chromosome 6? Also, should I dump the unknown chromosomes as I am specifically interested in chromosomes 1-22?

Any suggestion on splitting the fasta file into 22 chromosomes would be nice.

# Answer
The HLA region is already a part of chromosome 6. When you split the file, it will just be included with the rest of the chromosome. Whether or not you include the remaining contigs depends on your specific objectives. Inclusion of the unplaced and unlocalized scaffolds usually helps to reduce mapping errors, for example. To split/extract **only** chr1-22 from the UCSC hg38.fa.gz, we can just use . With the following in a file called `script.awk`:

```
BEGIN {
    for(i=1;i<=22;i++) {
        arr["chr" i]
    }
}

/^>/ {
    f=substr($0,2)
}

f in arr {
    print > f ".fa"
}

```

Run using:

```
awk -f script.awk <(zcat hg38.fa.gz)

```
> 2 votes

---
Tags: fasta, reference-genome
---


title: Breseq error (code 137). Any ideas?

text:
# Question
Title: Breseq error (code 137). Any ideas?
The code I ran was here, so nothing fancy.

```
breseq -l 110 -o RifR_align -r Big_burk_assembly.fasta RifRNano_nanopore.fastq.gz

```

This was the output.

```
NOW PROCESSING Read alignment to reference genome
[system] bowtie2-build -q "RifR_align/data/reference.fasta" "RifR_align/02_reference_alignment/reference"
    [system] bowtie2 -t --no-unal -p 1 -L 31 --ma 1 --mp 3 --np 0 --rdg 2,3 --rfg 2,3 --ignore-   quals --local -i S,1,0.25 --score-min L,1,0.9 -k 2000 --reorder -x "RifR_align/02_reference_alignment/reference" -U "RifR_align/01_sequence_conversion/RifRNano_nanopore.converted.fastq" -S "RifR_align/02_reference_alignment/RifRNano_nanopore.stage1.sam" --un "RifR_align/02_reference_alignment/RifRNano_nanopore.stage1.unmatched.fastq"
    Time loading reference: 00:00:00
    Time loading forward index: 00:00:00
    Time loading mirror index: 00:00:00
    Killed
    Use of uninitialized value $l in scalar chomp at /home/lts/miniconda3/envs/assemblers/bin/bowtie2 line 588, <BT> line 4626.
Use of uninitialized value $l in substitution (s///) at /home/lts/miniconda3/envs/assemblers/bin/bowtie2 line 589, <BT> line 4626.
Use of uninitialized value $l in print at /home/lts/miniconda3/envs/assemblers/bin/bowtie2 line 593, <BT> line 4626.
(ERR): bowtie2-align exited with value 137
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!> FATAL ERROR <!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

Error running command:
[system] bowtie2 -t --no-unal -p 1 -L 31 --ma 1 --mp 3 --np 0 --rdg 2,3 --rfg 2,3 --ignore-quals --local -i S,1,0.25 --score-min L,1,0.9 -k 2000 --reorder -x "RifR_align/02_reference_alignment/reference" -U "RifR_align/01_sequence_conversion/RifRNano_nanopore.converted.fastq" -S "RifR_align/02_reference_alignment/RifRNano_nanopore.stage1.sam" --un "RifR_align/02_reference_alignment/RifRNano_nanopore.stage1.unmatched.fastq"
Result code: 35072
FILE: libbreseq/common.h   LINE: 1411
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!> STACK TRACE <!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Backtrace with 6 stack frames.
breseq(+0x44cf4) [0x56102c3c3cf4]
breseq(+0x4a4f1) [0x56102c3c94f1]
breseq(+0x6d142) [0x56102c3ec142]
breseq(+0x36cc4) [0x56102c3b5cc4]
/lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0xf3) [0x7f71ba61e0b3]
breseq(+0x40559) [0x56102c3bf559]
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

```

The programs seemed to have no trouble aligning the reference genome, however the Bowtie2 package is throwing error code 137 and terminating the program.

# Answer
\[Not quite an answer but too long for a comment\]

I am guessing that your assembly `Big_burk_assembly.fasta` is a rather large assembly. Assembly indexing can be very memory intensive for complex assemblies. If it is eating enough memory to kill the bowtie2 indexing process (as noted in comment by @terdon), I'd consider filtering it, for example if there are lots of very short contigs that are unlikely to be meaningful, consider removing those.

Just looking at the way you've named the assembly file, it looks like it might be a somewhat ad hoc assembly, so maybe there is a meaningful way to filter it. We might be able to suggest more disciplined ways to do this if you tell us more about it.

Alternately, just get more memory as suggested in comments by @terdon. This will be hard if you are running on a personal computer, but on a cluster you can usually specify this in a job parameter file. Here (section "submit scripts") is how to set it on SLURM, for example.
> 1 votes

---
Tags: reference-genome, bowtie2
---


title: aligner for 1D^2 oxford nanopore data

text:
# Question
Title: aligner for 1D^2 oxford nanopore data
I am wondering if anyone has tried aligning 1d^2 RNA-seq data. I am trying that on minimap2. It is giving me very low ~50% of mapped reads. it gives ~95% mapped reads for 1D.

I am using the following command for 1d and 1d^2:

```
minimap2 -ax splice -uf -k14 ref.fa direct-rna.fq > aln.sam

```

Update:

I also tried to play around the command using :

```
minimap2 -ax splice -k14 -a ref.fa direct-rna.fq > aln.sam
minimap2 -ax splice -ur -k14 ref.fa direct-rna.fq > aln.sam
minimap2 -ax splice -ub -k14 ref.fa direct-rna.fq > aln.sam

```

It did not work better.

FYI: I am trying to find answer for this from few days. I have posted this question on Biostars before. I did not get any response. Hence, I am trying to find the answer here.

Thank you in advance. I really appreciate any help.

# Answer
> 1 votes
Since this is an old unanswered question where a lot of method development has since happened, going to put some answers here:

* This rather old tool also mentioned in comments (graphmap) is supposedly good for this, though it may not be relevant anymore.
* This paper uses minimap (not minimap2) alignments as a heuristic for read quality, suggesting that minimap2 should be good enough.

Note that depending on how you have processed (specifically basecalled) the data, 1d2 "second" reads may have low accuracy. That matches your ~50% unaligned proportion. You may need to do additional work during basecalling to fully use 1d2, as discussed in this paper.

Note: I am not an expert on nanopore and I put this answer together using google. Others are welcome to edit/editorialize.

---
Tags: sequence-alignment, nanopore
---


# Question
Title: dog coordinates (canFam3) to human coordinates (hg19)
I've converted dog coordinates to human using UCSC LiftOver. These are 200bp intergenic regions that are differentially methylated from normal dogs to cancer dogs. I've converted these to human coordinates and found that a lot of them overlap with differentially methylated regions we found in the human model.

1. Is this okay to do?
2. What stringencies should I check or modify, the min ratio of bases needed to overlap?
3. Lastly, is LiftOver taking the 200bp dog sequence and saying that at least 10% of the sequence has to align with a region within the human genome?

I could not find an answer online.

# Answer
Based on this old answer, this is not recommended for liftOver. For between-species liftovers, it is claimed that you are likely to want a different tool such as pslMap. For a lot of rather detailed discussion of how liftOver operates in this context and how to parameterize properly, see here.

See also this old Biostars question addressing the issue of mapping not between versions of the same assemblies but different genomes entirely. This thread also notes that using liftOver can be dangerous. That isn't to say that you can't learn something interesting, just that you might make mistakes.

pslMap apparently drops non-syntenic regions, which may lose a lot of data, but it is supposedly making rather fewer assumptions about assembly relationships.

orthoMap is apparently a similar tool to pslMap, but I can't find a way to run it except for a very old and restricted web interface. There are other tools with very similar names ("orthomap") but somewhat different uses, which may increase confusion.

A counterpoint to all the foregoing is that this rather detailed sum-up of liftover tools does not seem to be very concerned about the issues of using tools like liftOver between species. **It specifically uses dog to human liftover as the example, just as you want to**. It's notable that the assumption is made in that page that a high-quality chain file exists, which does seem to be true for you.

To go through your questions:

1. depends, see above.
2. For suggested parameters for running blat/chain generation/liftOver for this application, see here.
3. This is derived in part from the \[chain file\] I believe rather than liftOver parameters (http://genomewiki.ucsc.edu/index.php/LiftOver\_Howto). However overall this is set by the parameters you used in generating a chain file, which probably boils down to parameters of the workflow here. liftOver itself is just relying on the chain alignments for each region with the best score.
> 1 votes

---
Tags: ngs, methylation, sequence-annotation, liftover
---


# Question
Title: Where to assign a read pair when a mate spans a HiC bin boundary?
Chromosome capture experiments, like Hi-C, typically bin the genome into distinct bins. A contact matrix, $M$, is given by $M\_{i,j} = |\\{pairs | mate\_1 \in bin\_i, mate\_2 \in bin\_j\\}|$.

If one of the read mates spans a bin boundary, where should that pair be assigned?

For example, if $mate\_1$ overlaps $bin\_i$ and $mate\_2$ overlaps both $bin\_j$ and $bin\_{j+1}$, should this read pair be assigned to $M\_{i,j}$ or $M\_{i,j+1}$? Should it be randomly assigned to either bin with equal probability? Should it be assigned to both?

I know this may be a rare event, but if Hi-C datasets are generated using billions of reads (such as Rao *et al.*, Cell, 2014), this is bound to happen.

# Answer
> 1 votes
\[Old question, but possibly still relevant to answer.\]

This will not necessarily be rare. One or the other of the mates could span the ligation junction itself, in which case your read start and end could be entirely different chromosomes. In fact, some tools specifically measure the quality of Hi-C libraries by measuring the number of such ligations directly observable within reads. This will be a much larger number of reads than bin-crossing reads.

The safest thing to do in all such cases (IMO) is to take the read start position as the location of interest for each mate. That way, if a junction does occur, you still have an unambiguous position which is your best guess of where each end of the ligation product actually originates. This gets around both the minor issue of the bin-crossing problem and the major issue of the ligation junction problem. There is probably some imprecision here and you throw out some information, but sorting out the complex possible situations otherwise seems like a lot of hassle for minimal gain.

---
Tags: hi-c, chromosome-capture
---


# Question
Title: How can I resolve "missing name of population" message in PGDSpider?
I want to use Arlequin for AMOVA. I have my FASTA files which I want to convert into .arp files which are recognizable by arlequin. When I try to convert it using Population Genetics Data (format) Spider, or PGDSpider by Lisher and Excoffier (2018).

The error message says the name of my populations is missing.

Can anyone tell me what's wrong with my FASTA files? Check the following example:

```
>CcHr1
ATAAAGATATTGGTACGTTGTATATATTATTGGGTGTGTGGTGTGGTATAGTTGGTACAGGATTATCATT
ACTAATTCGTCTTGAGTTAGGTACTTCTGGTGTCTTGACTGATGAGCATTTTTATAATGTAATTGTAACA
GCCCACGCTTTTGTAATGATTTTTTTTATGGTTATACCGATTATAATTGGTGGTTTTGGAAATTGAATAG
TTCCTTTATTAATTGGTGCTCCAGACATAAGGTTCCCTCGGATAAATAATATAAGATTTTGATTATTACC
TCCTGCTTTTGTACTACTGATTAGTAGAAGCCTTGTTGAGGGAGGGGCTGGGACTGGCTGAACTGTTTAT
CCTCCACTTAGGAGGGTAACAGGACATAGTGGGGCTTCTGTAGACTTGGCTATTTTTTCCTTACATCTCG
CGGGCATATCATCAATTTTGGGGGCAATTAATTTTATTACAACGATCTTCAATATGCGGTCCCCTGGGTT
AACACTCGAGCGGATAAGGCTATTTGTATGGTCTATTTTAGTAACAGTTTTTTTATTGTTACTATCACTT
CCAGTTCTTGCAGGGGCAATTACGATGCTCTTAACTGATCGTAATTTTAGTACTTCATTTTTTGATCCAG
>CcMath1
ATAAAGATATTGGTACGTTATATATATTATTGGGTGTGTGGTGTGGAATAGTTGGCACAGGGTTGTCATT
ATTAATTCGTCTTGAGTTAGGTACTTCTGGTATCTTGACTGATGAGCATTTTTATAACGTAATTGTAACA
GCCCACGCTTTTGTAATAATTTTTTTTATGGTTATGCCAATTATAATTGGTGGTTTTGGAAATTGAATAG
TCCCTTTATTAATTGGTGCTCCAGACATAAGGTTCCCTCGAATAAATAATATAAGATTTTGACTATTACC
TCCTGCGTTTGTACTGCTGATTAGTAGAAGCCTTGTTGAGGGAGGGGCTGGAACTGGCTGAACTGTTTAT
CCTCCACTTAGAAGGGTAACAGGGCATAGTGGGGCTTCTGTAGACTTGGCGATCTTTTCCTTGCATCTCG
CGGGTATATCATCAATTTTGGGGGCAATTAATTTTATTACAACGATCTTCAATATGCGGTCCCCTGGGTT
AACACTCGAACGGATAAGGTTATTTGTATGGTCTATTTTAGTAACAGTTTTTTTATTGTTATTGTCACTT
CCAGTTCTTGCAGGGGCTATTACGATACTCTTAACTGATCGTAATTTTAGTACTTCGTTTTTTGACCCAG

```

As you can see I have to populations but it only names it as pop1.

# Answer
> 2 votes
You could try specifying the `population` keyword in the record header. From the PGDSpider manual (version 2.1.1.5 (May 2018)), at the bottom of page 67:

> Population definition in header line: String after “population:” until the first space is parsed as population name e.g.: `>sample1 population:pop1 gene:HMA4`

# Answer
> 1 votes
It's important to understand what Arlequin is doing for MANOVA. It will compare two or more populations to assess whether there is a difference in the allele frequency between the populations. It will also assess "interaction" (interpreting that in "gene flow" needs some thinking).

The basic theory is that if two or more populations cannot be distinguished under a formal parametric test there must be very strong gene flow supporting the lack of population structure between these populations.

**Geographic foci** Population designation is usually done between geographic foci, thus population A is usually geographic location (focus) A and population B is a geographically separated region (focus). Thus the fasta header needs to know the geographic region of the sample. Thats it.

**Population info** The issue is that information is required within the fasta header. Thus population delineation is therefore essential. This is confirmed in the PGDSpider manual below the line for Arlequin: which you don't need to do because thats what the conversion does.

Arlequin therefore needs to group each member of the population into a (usually geographically seperated) population to perform MANOVA. It will require this information specifying on the fasta header line. At the moment the data is saying there is a single population, it thinks all the sequences are from the same geographic location.

**Procedure** Thus it's likely simply a case of assigning two or more broad geographic foci for your samples. Then going through the fasta file and designating each sequence into one of the geographic focus under the keyword "population".

Just modifying @Steve's reading of the manual

```
> CcHr1 population:geoFocusA gene:X
ATAAAGATATTGGTACGTTGTATATATTATTGGGTGTGTGGTGTGGTATAGTTGGTACAGGATTATCATT

>CcMath1 population:geoFocusB gene:X
 ATAAAGATATTGGTACGTTATATATATTATTGGGTGTGTGGTGTGGAATAGTTGGCACAGGGTTGTCATT

```

Then run PGDSpider

```
java -Xmx1024m -Xms512m -jar PGDSpider2-cli.jar -inputfile examples\example_Structure.fasta -inputformat FASTA -outputfile examples\output_Arlequin.arp -outputformat ARLEQUIN -spid examples\Structure_Arlequin.spid

```

You did attain the critical anomaly in the sequence format,

* *each line of a sequence should have fewer than 80 characters*

Your rows have 70 residues.

---

> Genetic structure (only required for AMOVA):

> The genetic structure specifies the hierarchical genetic structure of the samples. It is possible to define groups of populations. o start of the subsection: \[\[Structure\]\] o name for the genetic structure (string within ” ”): StructureName = “name” o number of groups defined in the structure (int value): NbGroups = 5 o group definitions (list containing the names of the samples belonging to the group, entered within braces “{ }”): Id1 2 ACTCGGGTTCGCGCGC # the first pseudo-haplotype ACTCGGGCTCACGCGC # the second pseudo-haplotype

```
NbGroups=2
Group ={
  population1
}
Group ={
  population2
  population3
}

```

---

**Advantage** The advantage of MANOVA is you can have LOADS of populations and it will assess which have high levels of gene flow between them. Statistically it is fantastic for loads of different reasons (quite complicated reasons). The disadvantage of output of MANOVA is there will likely be significant differences between many of the populations - possibly all of them. This is a result of the power of parametric statistics, but may not be particularly helpful for the analysis. If that is the outcome switching to *Fst* is the next part of the analysis (and a different question). Goodluck and you can always post back if you get stuck.

---
Tags: phylogenetics, fasta, phylogeny, populations
---


title: How to analyse Orthofinder results

text:
# Question
Title: How to analyse Orthofinder results
I read this paper (*Colletotrichum shisoi sp. nov., an anthracnose pathogen of Perilla frutescens in Japan: molecular phylogenetic, morphological and genomic evidence*). We ran Orthofinder and also created an Upset plot. Unfortunately, I do not understand how they calculated the following percentages in the *"Conservation of proteins amongst Colletotrichum species"* section:

> A total of 11,914 orthogroups with two or more proteins were identified (Fig. 4b). Of these, 7,950 groups (74.0% proteins from C. shisoi, 63.0% from C. higginsianum, 78.3% from C. incanum, 74.2% from C. tofieldiae and 73.7% from C. graminicola) were conserved in all five species (Supplementary Tables S3 and S4)

How is it possible to calculate the values?

Thank you in advance

# Answer
I think they simply calculate the percentage of proteins in each species that belongs to an orthogroup conserved in the five species.

There should be a tab separated text file called Orthogroups.GeneCount.tsv (at least in the Orthofinder version used in the paper). This file contains counts of the number of genes for each species in each orthogroup. You can use this file to get the total number of genes in each species associated with those 7950 orthogroups and then calculate its proportion compared to the total number of genes.

Hope this helps
> 4 votes

---
Tags: clustering, orthologues, clusterprofiler, orthofinder
---


# Question
Title: File I/O error using nglview.show_biopython(structure)
So I have been trying to get into visualizing proteins in python, so after some research I ended up on a tutorial that was teaching you how to visualize a protein from the COVID-19 virus, so I went and setup anaconda, got jupyter notebook working in vscode, and downloaded the necessary files from the PDB database, and made sure they were in the same directory as my notebook but when I run the the `nglview.show_biopython(structure)` function I get `ValueError: I/O opertation on a closed file`. I'm stymied this is my first time using jupyter notebook so maybe there is something I'm missing, I don't know.

This what the code looks like

```
from Bio.PDB import *
import nglview as nv

parser = PDBParser()
structure = parser.get_structure("6YYT", "6YYT.pdb")
view = nv.show_biopython(structure)

```

The error looks like this

```
Output exceeds the size limit. Open the full output data in a text editor
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
~\AppData\Local\Temp\ipykernel_8776\2743687014.py in <module>
----> 1 view = nv.show_biopython(structure)

c:\Users\jerem\anaconda3\lib\site-packages\nglview\show.py in show_biopython(entity, **kwargs)
    450     '''
    451     entity = BiopythonStructure(entity)
--> 452     return NGLWidget(entity, **kwargs)
    453
    454

c:\Users\jerem\anaconda3\lib\site-packages\nglview\widget.py in __init__(self, structure, representations, parameters, **kwargs)
    243         else:
    244             if structure is not None:
--> 245                 self.add_structure(structure, **kwargs)
    246
    247         if representations:

c:\Users\jerem\anaconda3\lib\site-packages\nglview\widget.py in add_structure(self, structure, **kwargs)
   1111         if not isinstance(structure, Structure):
   1112             raise ValueError(f'{structure} is not an instance of Structure')
-> 1113         self._load_data(structure, **kwargs)
   1114         self._ngl_component_ids.append(structure.id)
   1115         if self.n_components > 1:
...
--> 200         return io_str.getvalue()
    201
    202

ValueError: I/O operation on closed file

```

I only get this error when using `nglview.show_biopython`, when I run the `get_structure()` function it appears to read the file just fine. I can visualize other molucles just fine, or maybe that's because I was using the ASE library instead of a file. I don't know, that's why I'm here.

Update: Recently I found out that I can visualize the protein using `nglview.show_file()` instead of using `nglview.show_biopython()`. Even though I can visualize proteins now and technically my problem has been solved I would still like to know why the `show_biopython()` function isn't working properly.

# Answer
> 4 votes
I also figured out another way to fix this problem. After going back to the tutorial I was talking about I saw that it was made back in 2021. After seeing this I wonder if we were using the same verions of each package, turns out we were not. I'm not sure what version of nglview they were using, but they were using biopython 1.79 which was the latest verion back in 2021 and I was using biopython 1.80. While using biopython 1.80 I was getting the error seen above. But now that I'm using biopython 1.79 I get this

So I guess there is something going on with biopython 1.80, so I'm going to stick with 1.79

# Answer
> 3 votes
I had a similar problem with:

```
from Bio.PDB import *
import nglview as nv

parser = PDBParser(QUIET = True)
structure = parser.get_structure("2ms2", "2ms2.pdb")

save_pdb = PDBIO()
save_pdb.set_structure(structure)
save_pdb.save('pdb_out.pdb')

view = nv.show_biopython(structure)
view

```

error was like in question:

```
.................site-packages/nglview/adaptor.py:201, in BiopythonStructure.get_structure_string(self)
    199 io_str = StringIO()
    200 io_pdb.save(io_str)
--> 201 return io_str.getvalue()

ValueError: I/O operation on closed file

```

I modified `site-packages/nglview/adaptor.py:201, in BiopythonStructure.get_structure_string(self)`:

```
def get_structure_string(self):
        from Bio.PDB import PDBIO
        from io import StringIO
        io_pdb = PDBIO()
        io_pdb.set_structure(self._entity)
        io_str = StringIO()
        io_pdb.save(io_str)
        return io_str.getvalue()

```

with :

```
def get_structure_string(self):
        from Bio.PDB import PDBIO

        import mmap

        io_pdb = PDBIO()

        io_pdb.set_structure(self._entity)

        mo = mmap_str()

        io_pdb.save(mo)

        return mo.read()

```

and added this new class `mmap_str()` , in same file:

```
import mmap
import copy

class mmap_str():

    import mmap #added import at top

    instance = None

    def __init__(self):

        self.mm = mmap.mmap(-1, 2)

        self.a = ''

        b = '\n'

        self.mm.write(b.encode(encoding = 'utf-8'))

        self.mm.seek(0)

        #print('self.mm.read().decode() ',self.mm.read().decode(encoding = 'utf-8'))

        self.mm.seek(0)

    def __new__(cls, *args, **kwargs):
        if not isinstance(cls.instance, cls):
            cls.instance = object.__new__(cls)
        return cls.instance

    def write(self, string):

        self.a = str(copy.deepcopy(self.mm.read().decode(encoding = 'utf-8'))).lstrip('\n')

        self.mm.seek(0)

        #print('a -> ', self.a)

        len_a = len(self.a)

        self.mm = mmap.mmap(-1, len(self.a)+len(string))

        #print('a :', self.a)

        #print('len self.mm ', len(self.mm))

        #print('lenght string : ', len(string))

        #print(bytes((self.a+string).encode()))

        self.mm.write(bytes((self.a+string).encode()))

        self.mm.seek(0)

        #print('written once ')

        #self.mm.seek(0)

    def read(self):

        self.mm.seek(0)

        a = self.mm.read().decode().lstrip('\n')

        self.mm.seek(0)

        return a

    def __enter__(self):

        return self

    def __exit__(self, *args):

        pass

```

if I uncomment the print statements I'll get the :

```
IOPub data rate exceeded.
The notebook server will temporarily stop sending output
to the client in order to avoid crashing it.

```

error , but commenting them out I get:

while using the`nglview.show_file(filename)` I get:

tha's because, as could be seen looking at the `pdb_out.pdb` file

outputted by my code, `Biopytho.PDB.PDBParser.get_structure(name , filename)` doesnt retrieve the pdb header responsible for generate full CRYSTALLOGRAPHIC SYMMETRY/or biopython can't handle it (not sure about this, help if you know better), but just the coordinates.

Still don't understand what is going on with the :

```
--> 201 return io_str.getvalue()

ValueError: I/O operation on closed file

```

it could be something related to jupiter ipykernal ? hope somebody could shed more light into this, dont know how the framework runs, but is definitely different from a normal python interpreter. As an example:

Same code in one of my Python virtualenv, will run forever, so it could be ipykernel dont like `StringIO()`s or do something strange to them ?

OK thanks to the hint in the answer below, I went inspecting PDBIO.py in github repo for version Biopython 1.80 and compared the save method of PDBIO : def save(self, file, select=\_select, write\_end=True, preserve\_atom\_numbering=False): with the one in Biopython 1.79,

see first bit:

and last bit:

so apparently the big difference is the `with fhandle:` block in version 1.80.

So I realized that changing `adaptor.py` with adding a subclass of `StringIO` that looks like:

```
from io import StringIO
class StringIO(StringIO):

    def __exit__(self, *args, **kwargs):

        print('exiting from subclassed StringIO !!!!!')

        pass

```

and modifying `def get_structure_string(self):` like this:

```
def get_structure_string(self):
        from Bio.PDB import PDBIO
        #from io import StringIO
        io_pdb = PDBIO()
        io_pdb.set_structure(self._entity)
        io_str = StringIO()
        io_pdb.save(io_str)
        return io_str.getvalue()

```

was enough to get my Biopython 1.80 work in jupiter with nglview.

That told I am not sure what are the pitfalls of not closing the StringIO object we use for the visualization, but apparently its what Biopython 1.79 was doing like my first answer using a mmap object was doing too (not closing the mmap\_str)

# Answer
> 3 votes
This is an issue with Biopython 1.80, in particular introduced by a commit that automatically closes IO objects sent to `PDBIO.save`.

I made a PR to change the behaviour which fixes this issue and should hopefully be in the next Biopython release. Until then, avoiding using Biopython 1.80 should work.

Thanks @jeremiah-wade and @pippo1980 for digging into this.

# Answer
> 1 votes
tried to understand git, I ended up with this, seems more coherent with the previous habits in the biopython project, but cant push it.

Anyone could pass it along ? \[of course needs to be checked out, my tests are OK, but I am a novice\].

---
Tags: python, structural-biology, biopython
---


# Question
Title: How can I fix "Ill-defined genetic structure" Error message in Arlequin?
I am performing an AMOVA analysis via Arlequin.

At this point, I have:

1. collected sequences of mtDNA from the *cox*1 gene,
2. aligned the different genes using clustalX
3. converted the files into .arp format.

I then altered the structure of my populations in Arlequin so that all the different populations are in one group. When I tried to run AMOVA I get the error message

> "Ill-Defined genetic structure".

Any clues as to what I have to adjust?

# Answer
> 0 votes
If one group is one population this will not work with Arlequin AMOVA. There must be more than one population (I might have misread the question however).

There could several reasons for this output. The most likely is there are insufficient SNPs within the cox1 locus for MANOVA to identify a population structure with sufficient statistical robustness. What it requires are a minimum of 3 independent alleles (SNPs), really 3 genotypes. If the SNPs are linked on the same bit of DNA that doesn't count - thats just one genotype. What I suspect is happening is there is a dearth (not enough) genetic information within cox1 to perform MANOVA.

There could be other explanations, but this is the first issue to address. This isn't just Arlequin it is standard theory for any allele frequency based calculations and why this is important is more complicated to explain.

---

Please note, it is the run-up to Christmas .. okay put Christmas on hold. To respond to the comments ...

Firstly, if you hold the data set it is very easy to identify the number of genotypes. There are past posts on this herein. I am not inclined to describe alternative approaches to an algorithm I've still not yet published :-) ... however it is really easy and needs an empirical test. It is good to assess number of genotypes *per se* Arlequin will do this analysis ... finding that within the output I dunno.

Extending loci is always good and for mtDNA I would not recommend with 16S nor cox2. rRNA should be avoided in allele frequency calculations, its okay for taxonomic assignment (bacteria rather than mtDNA). cox2 is a small locus. I would recommend for mtDNA

* cobB, formerly called cytB
* intergenic region of mtDNA

cobB is a good target, it has a lot more selection on it for within species analysis than cox1. The intergenic region is rarely used, its a good locus because of the higher rates of mutation therein.

---
Tags: phylogenetics, phylogeny, allele-frequency, populations
---


# Question
Title: RNA folding at specific temperature with ViennaRNA in python
I am trying to get dot-bracket-notations of single stranded RNA via the ViennaRNA python package (https://pypi.org/project/ViennaRNA/) at different temperatures. I have read in the docs (https://www.tbi.univie.ac.at/RNA/ViennaRNA/doc/html/examples\_c.html) that the folding prediction at different temperatures is possible with C..but nowhere it is described how to specify the temperature when using python. Does anyone know how to specify the temperature at which the RNA folding is predicted with the ViennaRNA package in python?

Thanks for the help!

# Answer
> 1 votes
I think that the "Global Variables" section of the docs describes this:

> For the Python interface(s) SWIG places global variables of the C-library into an additional namespace cvar. For instance, changing the global temperature variable thus becomes

```
RNA.cvar.temperature = 25

```

Globals are somewhat surprising to me as a way of handling this. It seems like you should be able to specify this differently. Maybe I'm missing something?

---
Tags: python, rna, secondary-structure
---


# Question
Title: How can I use sequences with different lengths in arlequin?
I have mtDNA sequences of the *cox1* gene from various populations of one species. I want to perform an AMOVA allele frequency analysis via Arlequin. This allele frequency package however requires my sequences to have the same length.

* Is it possible to use my sequences?
* Do I have to align them first maybe?

The length difference is most likely due to indels/deletions. Do I have to use specific parts of the sequence for each gene of the subunit?

# Answer
All that is needed is to align the sequences within a formal aligner like MAFFT or muscle. Then pass that through the PDGSpider converter. Arlequin should just take it from there.

Manual annotation of the alignment is normal e.g. `seqoTron` and is traditionally regarded as the best approach rather than do this algorithmically. Thus I'm suggesting (strongly suggesting) do this prior the PGDSpider conversion. If an alignment has not been performed Arlequin should throw an error.

If there are trailing gaps (indels or '-') at the beginning and end of the alignment Arlequin will/should ignore these. You can of course simply cut them out with the alignment editor (Mega X will also work for alignment editing). The important point is that Arlequin must correctly identify homologues its an allele frequency package and nothing more. If there are errors in the alignment Arlequin will believe there is population structure when it could be absent.

I can't honestly believe that Arlequin will count an indel '-' as an informative site. It should just ignore those sites, because it is a very respectable package that has been around for years and years.

Okay, I am now signing off for Christmas... with occasional exceptions.
> 0 votes

---
Tags: phylogenetics, genetics, allele-frequency, populations
---


title: bsub with -M and -hl a good idea?

text:
# Question
Title: bsub with -M and -hl a good idea?
I have been submitting my snakemake pipeline to a node-shared high performance computer (HPC) via bsub the following way:

```
bsub -W {resources.runtime_min} -n {threads} -M {resources.mem_mb} -R 'rusage[mem={resources.mem_mb}]' -o run_logs/{name}_%J.out -e run_logs/{name}_%J.err"

```

Works great, no problem.

However, I noticed that the memory limit was `tasks*memory` given. I am not aware of any memory issues on the HPC. The "issue" was not noticed until I asked why my memory limit was so high compared to what I gave with the `bsub` flag `-M` (I wanted to understand how everything worked).

I was recommend to use a `-hl` flag to hard cap the usage. It was actually an administator that recommended it, after I asked him about how it worked. They now want to implement it for other users. I found this in the documentation too, so I understand this.

But I see no one else is using this and it is also not suggested when you want to use the `bsub` outfiles to understand how much resources you need (how much memory should I use).

Does anyone have thoughts on this?

Somehow, it makes sense to use `-hl`, but I am hesitant since no one seems to use `-hl` for some reason, and I feel like I do not fully understand it.

I'm asking here, because I felt I kept bugging the administrator with the same questions :'-). Thank you in advance for any advice.

# Answer
> 1 votes
What this question appears to be about is capping the RAM usage on a node within a shared HPC.

Being 'nice' on an HPC is very honourable indeed, it is however the system administrators jurisdiction. It simply requires an email to them and they will provide the authoritative answer. Thus if there were RAM usage issues within the nodes then they will be aware of it. Limiting RAM on a job is risky because the biggest 'fall-over' in bioinformatics is a job running into a RAM bottleneck. If that happens the process just freezes and literally will take 1 log more run-time. If the machine doesn't have enough RAM for the job, thats a separate issue and easily solvable via cloud computing.

The individual who made the suggestion might be running into RAM bottleneck issues on the machine. However, they should inform the HPC administrator who could then provide an assessment whether all other users should cap their RAM usage. There are alternatives - the user experiencing RAM bottlenecks can parallelise their job and this will overcome RAM issues. This can be a bit tricky on a course-grain parallelisation with a shared HPC.

So in summary, I wouldn't do this personally and if I was concerned I would seek the advice of the administrator who will give the definitive response. If admin says 'do it' then do it, but it your jobs take way-way longer than they should I'd relax the constraint and inform admin.

Admin could enforce this constraint without the user base being aware of what was happening if it was a major issue on the server.

These days if there are HPC issues its easier to shift on a cloud computer and is not worth the arguments when everyone is trying to get their jobs done before the vacation.

Merry Christmas!

---

From the comments ...

* System admin recommended this ... okay they should now inform the rest of the user base (which is what you are implying), unless there are been a specific issue. Thats a bit different and it is polite to comply with their request;
* What you can do is justify to the system admin why you job requires extra RAM;
* Systems admin might have automatically restricted the RAM access on the account - you can ask them to lift that;

You said,

> "1 log more run-time"?

That is extremely important to understand. Bioinformatic jobs can use loads of RAM, if the job exceeds the RAM capacity of the node essentially the job "freezes". What it really does is slow the time to completion of the job dramatically, so I'm saying in my experience it's about 10-fold longer to complete than normal. That is what I am calling a RAM-bottleneck.

**Final point** On the `bsub` parallelisation appears to be efficiently implemented and if RAM issues are occurring parallelising the job is the best solution. This would be a work around if the systems admin wants to restrict RAM usage and as a user you are wanting quick run-time (fast completion). Not all algorithms will parallelise but most will these days.

---
Tags: linux, hpc
---


# Question
Title: DESeq2 for large number of samples takes too much RAM
I am trying to run a very large number of transposase-accessible chromatin (ATAC)-seq samples through DESeq2 to find peaks of differential chromatin accessible across my study genome.

I have about 100 samples, and 1,000,000 "genes" which are really peaks.

Doing this takes a significant amount of RAM; right now I am running on a 1 TB bigmem node and still getting an out-of-memory error. The code works perfectly fine for a small subset of samples.

Does anyone have an alternative strategy to use DESeq2 on a very large dataset? Should I look into doing a Wilcox-rank test myself? Thanks.

# Answer
> 4 votes
I'll preface this by saying that I don't think DESeq2 is the right tool to use for ATAC-Seq data. My own study of ATAC-Seq patterns \[admittedly only a couple of runs that were our first exploration of the technique\] indicates that it has a fractal expression which is difficult to properly summarise at any scale, especially at a bulk sample level. I found this poster from a quick search on Twitter, which recommends AIAP for ATAC-Seq analysis (but don't read too much into that; I have no experience with doing ATAC-Seq analysis). What I *do* have a lot of experience with is DESeq2.

DESeq2 has certain assumptions about the expression patterns of genes which are built into its differential expression model, and the further you get from gene-level analysis, the less appropriate DESeq2 is for data analysis (for example, Swish is a better approach for isoform level and allelic analysis).

However, if you're absolutely set on using DESeq2... the DESeq2 data analysis documentation has a few different suggestions for speeding up analysis:

> The above steps should take less than 30 seconds for most analyses. For experiments with complex designs and many samples (e.g. dozens of coefficients, ~100s of samples), one may want to have faster computation than provided by the default run of DESeq. We have two recommendations:
>
> 1. By using the argument fitType="glmGamPoi", one can leverage the faster NB GLM engine written by Constantin Ahlmann-Eltze. Note that glmGamPoi’s interface in DESeq2 requires use of test="LRT" and specification of a reduced design.
> 2. One can take advantage of parallelized computation. Parallelizing DESeq, results, and lfcShrink can be easily accomplished by loading the BiocParallel package, and then setting the following arguments: parallel=TRUE and BPPARAM=MulticoreParam(4), for example, splitting the job over 4 cores. However, some words of advice on parallelization: first, it is recommend to filter genes where all samples have low counts, to avoid sending data unnecessarily to child processes, when those genes have low power and will be independently filtered anyway; secondly, there is often diminishing returns for adding more cores due to overhead of sending data to child processes, therefore I recommend first starting with small number of additional cores. Note that obtaining results for coefficients or contrasts listed in resultsNames(dds) is fast and will not need parallelization. As an alternative to BPPARAM, one can register cores at the beginning of an analysis, and then just specify parallel=TRUE to the functions when called.

Also, applying a low count filter can help a lot.

> **Pre-filtering**
>
> While it is not necessary to pre-filter low count genes before running the DESeq2 functions, there are two reasons which make pre-filtering useful: by removing rows in which there are very few reads, ***we reduce the memory size of the dds data object, and we increase the speed of the transformation and testing functions within DESeq2.*** It can also improve visualizations, as features with no information for differential expression are not plotted.
>
> Here we perform a minimal pre-filtering to keep only rows that have at least 10 reads total. Note that more strict filtering to increase power is automatically applied via independent filtering on the mean of normalized counts within the results function.

```
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]

```

> Alternatively, a popular filter is to ensure at least X samples with a count of 10 or more, where X can be chosen as the sample size of the smallest group of samples:

```
# unevaluated chunk...
keep <- rowSums(counts(dds) >= 10) >= X
dds <- dds[keep,]

```

# Answer
> 1 votes
Absolutely, look at Wilcox-rank, not purely because of the computational overhead, but the parametric assumption that DESeq2 uses are not correct on large data sets. This might describe your data sets too:

> **Abstract** DESeq2 and edgeR, have unexpectedly high false discovery rates.
>
> Expanding the analysis to limma-voom, NOISeq, dearseq, and **Wilcoxon rank-sum test**, we found that \[false discovery rates\] FDR control is often failed **except for the Wilcoxon rank-sum test**
>
> Particularly, the actual **FDRs of DESeq2** and edgeR sometimes exceed **20%** when the target FDR is 5%. Based on these results, for population-level RNA-seq studies with large sample sizes, **we recommend the Wilcoxon rank-sum test**.
>
> **Results and Discussion** DESeq2 and edgeR had only an 8% overlap in the DEGs they identified

---

---
Tags: r, deseq2, differential-expression, edger, atac-seq
---


# Question
Title: What reasons are there to choose Illumina if PacBio provides longer and better reads?
PacBio provides longer read length than Illumina's short-length reads. Longer reads offer better opportunity for genome assembly, structural variant calling. It is not worse than short reads for calling SNP/indels, quantifying transcripts. Sounds like PacBio can do whatever Illumina platform can offer.

Would that be any reason for going with Illumina if PacBio's long reads can do everything and more?

# Answer
There are **so** many reasons why one might want to prefer Illumina over PacBio (also note that it's a false dichotomy, at least Oxford Nanopore is a competitive sequencing platform):

* The first (IMHO and the most common reason) is still the cost of both sequencing and the instruments. Illumina can sequence a Gbp of data for \\$7 - \\$93. PacBio sequencing is according the same webpage \\$115 per Gbp, however at our sequencing center it's ~$200. Though ONT might have already a cheaper solution. Edit, I just found a google sheat with prices that seems to be frequently updated, seems that the ratio still holds Illumina short reads ~10x cheaper than PacBio.
* RNA-seq (i.e. analysis of a gene expression) is not possible with PacBio due to preferential sequencing of smaller fragments; shorter genes would always be shown to be more expressed. To be clear, it's possible to sequence RNA with PacBio (the keyword is iso-seq), but the analysis of gene expression is problematic.
* It's way easier to extract fragmented DNA (concerns small non-model organisms; although recently a single mosquito was sequenced, so we can expect a further improvement)
* other sequencing techniques as RAD-seq that allow genotyping with very little effort and cost, I have never seen anybody even considering using long reads for such genotyping
* Genome profiling (assembly-less genome evaluation) based on kmer spectra analysis is not possible with PacBio data due to higher error rates. Conflict of interest: I am a developer of one of the tools for genome profiling (smudgeplot)
* DNA in Formalin-Fixed Paraffin-Embedded samples is fragmented anyway (~300bp), therefore, there is no point in sequencing them with more expensive long-read technology (contributed by @mRoten)
* I bet there will be a plenty other applications I am not aware of
> 16 votes

# Answer
Many analyses performed on Illumina machines these days require large numbers of reads. For example, most analyses in ChIP-seq, RNA-seq, ATAC-seq etc, need 10s or even 100s of millions of reads for the statistics to work out properly.

But this isn't limited to just sequencing as an assay experiments. High depths are important for things like somatic variant calling as well. Or simply sequencing of medical gene pannels, where you might only have 20 amplicons of a 500bp each, but need to do it on 100s of patients as cheaply as possible.
> 5 votes

# Answer
Three reasons for Illumina:

* Much better for a large number of samples (easily handle 96 samples).
* SNP calling is much better - much greater depth
* Hardware costs, an Illumina MiSeq machine is cheap

PacBio SNP calling has been improved towards the standard of Illumina by DNA modification to the inputs to be sequenced but it depends on the application.
> 3 votes

# Answer
As others have mentioned, cost is probably one of the main issue.

Other than that, PacBio sequencing also requires more input DNA material (~3-10 ug for human) while Illumina sequencing works with just 100ng. This makes PacBio sequencing more challenging when samples are rare. I'm sure people at PacBio are trying to resolve some of these concerns already though.
> 0 votes

---
Tags: phylogenetics, ngs, phylogeny, long-reads, pacbio
---


title: Grep peculiar pattern inside a file

text:
# Question
Title: Grep peculiar pattern inside a file
I would like to grep everything inside the second square bracket of a file:

```
>AB1311658.1 unnamed protein product [C.sp. balc.]

>evm.model.pilon.8 [actin][D.r.]

>NP_0232343.1 [fucose mutarotase] [E.l.]

```

Obtaining:

```
C.sp. balc.

D.r.

E.l.

```

How can I solve it?

Marco

# Answer
> 4 votes
One way using :

```
sed -n '/^>/s/.*\[\([^]]*\)\]$/\1/p' input.fa

```

Results:

```
C.sp. balc.
D.r.
E.l.

```

Explanation:

* Turn off automatic printing of the pattern space using `-n`
* Address lines starting with a header using `/^>/`
* Perform a substitution using a capture group to select everything between the last set of brackets on these lines, with `s/.*\[\([^]]*\)\]$/\1/`
* Print the pattern space with `p`

# Answer
> 3 votes
With GNU `grep` (the default on Linux):

```
$ grep -oP '.*\[\K[^]]+' file
C.sp. balc.
D.r.
E.l.

```

The `-o` means "print only the matched portion of the line" and `-P` enables Perl Compatible Regular Expressions which give us `\K` for "ignore everything matched up to this point". So we then look for the longest possible string ending in a `[` and discard it (`.*\[\K`), and then look for the longest possible stretch of non-`]` characters after that.

For very large files, you can speed things up a little by focusing only on the header lines:

```
grep -oP '^>.*\[\K[^]]+' file

```

# Answer
> 2 votes
For OSX, I don't have a `-P` option (Perl pattern matching) so the best I can do is finish the one-liner by piping into `perl`

```
grep -o '\[.*\]\[[^]]*' input.fa | perl -p -e 's/[][]/ /g'

```

The `-o` makes `grep` match the substring, rather than the entire line.

With GNU `grep` (the default on Linux) and the `-P` option, it's a single line of grep where `\K` or `(?=pattern)`, `(?<=pattern)` are used for non-capture characters and these are placed around the `[`

```
grep -oP '(?<=\[).*(?=\]).*(?<=\[)[^]]*' input.fa

```

The `?<=` is a "positive lookbehind" and the `?=` is a "positive lookahead".

It's a lot easier to use a `perl` one-liner.

# Answer
> 1 votes
Some of the solutions suggested above utilize GNU `grep.` If you are a Mac user, you can easily install GNU version of the commands into your Mac. You can do this by installing the GNU Command Line Tools. The easiest way is to use Homebrew, then `brew install coreutils.`

The default commands will still be MacOS version (unless you intentionally set otherwise), but when you want to use GNU version, you could just add a `g` before the command. For example, `sort` vs `gsort`, `grep` vs `ggrep`. You can also check manual of the command the same way, e.g. `man gsort`. Trust me this will make your life easier as a Mac user.

Also, in recent years default MacOS Terminal shell is `zsh.` If you desire you can switch it back to `bash`.

---
Tags: fasta, genome, grep
---


# Question
Title: How does one account for both single end and paired end reads as input in a Nextflow pipeline?
Suppose you have a pipeline that takes reads as input from the user and trims them.

```
// Declare syntax version
nextflow.enable.dsl=2

/*
 * pipeline input parameters
 */

params.reads = "$projectDir/data/*_{1,2}.fastq"

/*
 * define the `TRIM` process that creates trimmed.reads.fastq files
 * given the reads.fastq file.
 */

process TRIM {

    tag "$pair_id"

    input:
    tuple val(pair_id), path(reads)

    output:
    path(pair_id)

    script:
    """
    atria -r ${reads[0]} -R ${reads[1]} -t 1 -l 12 -q 30 -o $pair_id
    """
}

workflow {

    read_pairs_ch = Channel.fromFilePairs( params.reads, checkIfExists:true )

    trim_ch=TRIM(read_pairs_ch)

}

```

This particular one takes paired end reads as input. What if the user wanted to trim both single end and paired end reads, how to account for single end reads as the input?

Best wishes.

# Answer
> 2 votes
There are many ways to accomplish this. Assuming you don't need to process single-end and paired-end reads simultaneously (i.e. within the same run), one way using the `fromFilePairs` factory method is to just use an additional flag to define the number of files each emitted item is expected to hold. And then adjust your process definition accordingly:

```
params.reads = "${projectDir}/data/*_{1,2}.fastq.gz"
params.singleEnd = false

process TRIM {

    tag "${pair_id}"

    cpus 1

    input:
    tuple val(pair_id), path(reads)

    output:
    tuple val(pair_id), path("${pair_id}/*.atria.fastq{,.gz}")

    script:
    def single = reads instanceof Path

    def read1 = !single ? /-r "${reads[0]}"/ : /-r "${reads}"/
    def read2 = !single ? /-R "${reads[1]}"/ : ''

    """
    atria \\
        -t ${task.cpus} \\
        ${read1} \\
        ${read2} \\
        -l 12 \\
        -q 30 \\
        -o "${pair_id}"
    """
}

```

```
workflow {

    Channel
        .fromFilePairs( params.reads, size: params.singleEnd ? 1 : 2 )
        .ifEmpty { error "Could not find any reads matching pattern: ${params.reads}" }
        .set { reads }

    TRIM( reads ).view()
}

```

Although we've used an additional flag, it does let the pipeline fail immediately if it can't find the files it expects. I think this approach is more explicit and provides a better user experience:

```
$ nextflow run main.nf --reads './data/*.fastq'
N E X T F L O W  ~  version 22.10.0
Launching `main.nf` [nasty_allen] DSL2 - revision: 2cbc311974
Could not find any reads matching pattern: ./data/*.fastq

 -- Check script 'main.nf' at line: 38 or see '.nextflow.log' file for more details

```

Doh! I forgot that these are single ended reads:

```
$ nextflow run main.nf --reads './data/*.fastq' --singleEnd
N E X T F L O W  ~  version 22.10.0
Launching `main.nf` [clever_gutenberg] DSL2 - revision: 2cbc311974
executor >  local (3)
[8e/f90c90] process > TRIM (X) [100%] 3 of 3 ✔
[Y, /path/to/work/e6/48df1485bb7de8355ed8ce0f1b53da/Y/Y.atria.fastq]
[Z, /path/to/work/eb/aca3389344ac03ff283ded0c46a7a6/Z/Z.atria.fastq]
[X, /path/to/work/8e/f90c90210103c03a3fbde697e01227/X/X.atria.fastq]

```

# Answer
> 2 votes
My personal strategy is to always read data from a samplesheet which then auto-detects whether a sample is paired-end or single-end and stores this information in a so-called meta map (inspired by nf-core) which then passes this information to the modules.

Example:

Say you have a samplesheet that is a three-column csv file with a header:

```
$ cat samplesheet.csv
sample,r1,r2
sample1,/Users/atpoint/sample1_R1.fq.gz,/Users/atpoint/sample1_R2.fq.gz
sample2,/Users/atpoint/sample2_R1.fq.gz,

```

Read this into Nextflow and detect whether r2 has an entry. If so it is paired-end, else it is single-end:

```

# main.nf

nextflow.enable.dsl = 2

// Channel for the samplesheet
ch_samplesheet = Channel.fromPath("/Users/atpoint/samplesheet.csv")

// Parse it line by line
ch_reads = ch_samplesheet.splitCsv(header:true).map {

    // This is the read1 and read2 entry
    r1 = it['r1']
    r2 = it['r2']

    // Detect wiether single-end or paired-end
    is_singleEnd = r2.toString()=='' ? true : false

    // The "meta" map, which is a Nextflow/Groovy map with id (the sample name) and a single_end logical entry
    meta = [id: it['sample'], single_end: is_singleEnd]

    // We return a nested map, the first entry is the meta map, the second one is the read(s)
    r2.toString()=='' ? [meta, [r1]] : [meta, [r1, r2]]

}

```

If you `.view()` the `ch_reads` you get:

```
[[id:sample1, single_end:false], [/Users/atpoint/sample1_R1.fq.gz, /Users/atpoint/sample1_R2.fq.gz]]
[[id:sample2, single_end:true], [/Users/atpoint/sample2_R1.fq.gz]]

```

A process that can use that could be:

```

process A {

    tag "$meta.id"

    input:
    tuple val(meta), path(reads)

    script:

    // The meta map holds the single/paired information that can now easily be parsed in the module in a standardized way
    if(meta.single_end){

        """
        do_single_end_things with ${reads}
        """

    } else {

        """
        do_paired_end_things ${reads}[0] ${reads}[1]
        """

    }

}

A(ch_reads)

```

The meta map idea (cudos nf-core) is simple, yet versatile and powerful as it allows standardized parsing of such metadata-like information (single\_end, libtype, id/sample\_name, any other value you can think of) and most importantly, it allows to process single-end paired end data at the same time. All the decision making is based on the meta map and therefore automated if the modules are coded appropriately as above.

Hope that helps.

---
Tags: nextflow
---