title: How to install blasr on Mac OS X?

text:
# Question
Title: How to install blasr on Mac OS X?
I am trying to install blasr on mac osx.

I followed the steps here

I get the error below after using the 'make' command:

```
make PREFIX=/tmp/pitchfork blasr
Cloning into 'pitchfork'...
remote: Counting objects: 10203, done.
remote: Compressing objects: 100% (51/51), done.
remote: Total 10203 (delta 25), reused 51 (delta 18), pack-reused 10128
Receiving objects: 100% (10203/10203), 1.16 MiB | 207.00 KiB/s, done.
Resolving deltas: 100% (5413/5413), done.
/Applications/Xcode.app/Contents/Developer/usr/bin/make -f mk/init.mk init
/Users/cr517/Documents/pitchfork/bin/checkSystem
/Users/cr517/Documents/pitchfork/bin/checkCC gcc
mkdir -p "/Users/cr517/Documents/pitchfork/workspace"
mkdir -p "/Users/cr517/Documents/pitchfork/distfiles"
mkdir -p "/Users/cr517/Documents/pitchfork/staging"
mkdir -p "/tmp/pitchfork/bin"
mkdir -p "/tmp/pitchfork/etc"
mkdir -p "/tmp/pitchfork/include"
mkdir -p "/tmp/pitchfork/lib"
mkdir -p "/tmp/pitchfork/lib/pkgconfig"
mkdir -p "/tmp/pitchfork/share"
mkdir -p "/tmp/pitchfork/var/pkg"
echo "export DYLD_LIBRARY_PATH=/tmp/pitchfork/lib:\$DYLD_LIBRARY_PATH"|sed -e 's/::*/:/g' > "/tmp/pitchfork/setup-env.sh"
echo "export PATH=/tmp/pitchfork/bin:\$PATH"|sed -e 's/::*/:/g' >> "/tmp/pitchfork/setup-env.sh"
echo "OLDPS1=\$PS1" >> "/tmp/pitchfork/setup-env.sh"
echo "[[ \"\$PS1\" ]] && PS1=\"pitchfork($(/Users/cr517/Documents/pitchfork/bin/pitchfork short-branch)) \$PS1\" || PS1=\($(/Users/cr517/Documents/pitchfork/bin/pitchfork short-branch)\)" >> "/tmp/pitchfork/setup-env.sh"
touch mk//../initialized.o
/Applications/Xcode.app/Contents/Developer/usr/bin/make -C ports/thirdparty/zlib do-install
make[1]: Nothing to be done for `do-install'.
/Applications/Xcode.app/Contents/Developer/usr/bin/make -C ports/thirdparty/libbzip2 do-install
rm -rf /Users/cr517/Documents/pitchfork/staging/libbzip2-1.0.6
mkdir -p /Users/cr517/Documents/pitchfork/staging/libbzip2-1.0.6/lib/pkgconfig /Users/cr517/Documents/pitchfork/staging/libbzip2-1.0.6/include /Users/cr517/Documents/pitchfork/staging/libbzip2-1.0.6/bin
cp -aL /usr/bin/bzip2 /Users/cr517/Documents/pitchfork/staging/libbzip2-1.0.6/bin/
cp -aL /usr/include/bzlib.h /Users/cr517/Documents/pitchfork/staging/libbzip2-1.0.6/include
cp -aL /usr/lib/libbz2.* /Users/cr517/Documents/pitchfork/staging/libbzip2-1.0.6/lib
ln -sfn libbz2.dylib.1.0 /Users/cr517/Documents/pitchfork/staging/libbzip2-1.0.6/lib/libbz2.dylib.1
ln -sfn libbz2.dylib.1   /Users/cr517/Documents/pitchfork/staging/libbzip2-1.0.6/lib/libbz2.dylib
sed -e "s@^prefix=.*@prefix=/tmp/pitchfork@" bzip2.pc > /Users/cr517/Documents/pitchfork/staging/libbzip2-1.0.6/lib/pkgconfig/bzip2.pc
cd /Users/cr517/Documents/pitchfork/staging/libbzip2-1.0.6 && tar cf - * | tar xf - -C /tmp/pitchfork
find /Users/cr517/Documents/pitchfork/staging/libbzip2-1.0.6 ! -type d|awk -F '/Users/cr517/Documents/pitchfork/staging/libbzip2-1.0.6/' '{print $2}' > /tmp/pitchfork/var/pkg/libbzip2-1.0.6
/Applications/Xcode.app/Contents/Developer/usr/bin/make -C ports/thirdparty/boost do-install
/Users/cr517/Documents/pitchfork/bin/pitchfork fetch --url https://prdownloads.sourceforge.net/boost/boost_1_60_0.tar.gz
fetching https://prdownloads.sourceforge.net/boost/boost_1_60_0.tar.gz
{'url': 'https://prdownloads.sourceforge.net/boost/boost_1_60_0.tar.gz', 'subcommand': 'fetch'}
Traceback (most recent call last):
  File "/Users/cr517/Documents/pitchfork/bin/pitchfork", line 525, in <module>
    sys.exit(main())
  File "/Users/cr517/Documents/pitchfork/bin/pitchfork", line 518, in main
    pitchfork(pargs)
  File "/Users/cr517/Documents/pitchfork/bin/pitchfork", line 444, in pitchfork
    raise Exception("[ERROR] %s" % ''.join(_err))
Exception: [ERROR] --2018-01-24 00:35:30--  https://prdownloads.sourceforge.net/boost/boost_1_60_0.tar.gz
Resolving prdownloads.sourceforge.net... 216.34.181.59
Connecting to prdownloads.sourceforge.net|216.34.181.59|:443... connected.
ERROR: cannot verify prdownloads.sourceforge.net's certificate, issued by ‘CN=COMODO RSA Domain Validation Secure Server CA,O=COMODO CA Limited,L=Salford,ST=Greater Manchester,C=GB’:
  Unable to locally verify the issuer's authority.
To connect to prdownloads.sourceforge.net insecurely, use `--no-check-certificate'.
make[1]: *** [boost_1_60_0.tar.gz] Error 1
make: *** [boost] Error 2

```

I think the important info is here:

```
raise Exception("[ERROR] %s" % ''.join(_err))
    Exception: [ERROR] --2018-01-24 00:35:30--  https://prdownloads.sourceforge.net/boost/boost_1_60_0.tar.gz
    Resolving prdownloads.sourceforge.net... 216.34.181.59
    Connecting to prdownloads.sourceforge.net|216.34.181.59|:443... connected.
    ERROR: cannot verify prdownloads.sourceforge.net's certificate, issued by ‘CN=COMODO RSA Domain Validation Secure Server CA,O=COMODO CA Limited,L=Salford,ST=Greater Manchester,C=GB’:
      Unable to locally verify the issuer's authority.
    To connect to prdownloads.sourceforge.net insecurely, use `--no-check-certificate'.

```

but how can I specify the --no-check-certificate option? That's not a make option, I've tried it.

I tried homebrew but got this message:

```
Updating Homebrew...
==> Auto-updated Homebrew!
Updated 2 taps (caskroom/cask, caskroom/versions).
Warning: Calling 'depends_on :mysql' is deprecated!
Use 'depends_on "mysql"' instead.
/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/blat.rb:19:in `<class:Blat>'
Please report this to the brewsci/science tap!

Warning: Calling 'depends_on :python3' is deprecated!
Use 'depends_on "python3"' instead.
/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/matplotlib.rb:56:in `<class:Matplotlib>'
Please report this to the brewsci/science tap!

Warning: Calling <<-EOS.undent is deprecated!
Use <<~EOS instead.
/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/plink.rb:23:in `block in <class:Plink>'
Please report this to the brewsci/science tap!

Warning: Calling <<-EOS.undent is deprecated!
Use <<~EOS instead.
/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/plink.rb:23:in `block in <class:Plink>'
Please report this to the brewsci/science tap!

Warning: Calling <<-EOS.undent is deprecated!
Use <<~EOS instead.
/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/plink.rb:23:in `block in <class:Plink>'
Please report this to the brewsci/science tap!

Warning: Calling <<-EOS.undent is deprecated!
Use <<~EOS instead.
/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/blasr.rb:27:in `block in <class:Blasr>'
Please report this to the brewsci/science tap!

Warning: Calling <<-EOS.undent is deprecated!
Use <<~EOS instead.
/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/blasr.rb:27:in `block in <class:Blasr>'
Please report this to the brewsci/science tap!

Warning: Calling <<-EOS.undent is deprecated!
Use <<~EOS instead.
/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/blasr.rb:27:in `block in <class:Blasr>'
Please report this to the brewsci/science tap!

Warning: Calling <<-EOS.undent is deprecated!
Use <<~EOS instead.
/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/blasr.rb:36:in `block in <class:Blasr>'
Please report this to the brewsci/science tap!

Warning: Calling <<-EOS.undent is deprecated!
Use <<~EOS instead.
/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/blasr.rb:36:in `block in <class:Blasr>'
Please report this to the brewsci/science tap!

Warning: Calling <<-EOS.undent is deprecated!
Use <<~EOS instead.
/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/blasr.rb:36:in `block in <class:Blasr>'
Please report this to the brewsci/science tap!

==> Installing blasr from brewsci/science
Warning: Calling <<-EOS.undent is deprecated!
Use <<~EOS instead.
/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/blasr.rb:27:in `block in <class:Blasr>'
Please report this to the brewsci/science tap!

Warning: Calling <<-EOS.undent is deprecated!
Use <<~EOS instead.
/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/blasr.rb:27:in `block in <class:Blasr>'
Please report this to the brewsci/science tap!

Warning: Calling <<-EOS.undent is deprecated!
Use <<~EOS instead.
/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/blasr.rb:27:in `block in <class:Blasr>'
Please report this to the brewsci/science tap!

Warning: Calling <<-EOS.undent is deprecated!
Use <<~EOS instead.
/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/blasr.rb:36:in `block in <class:Blasr>'
Please report this to the brewsci/science tap!

Warning: Calling <<-EOS.undent is deprecated!
Use <<~EOS instead.
/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/blasr.rb:36:in `block in <class:Blasr>'
Please report this to the brewsci/science tap!

Warning: Calling <<-EOS.undent is deprecated!
Use <<~EOS instead.
/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science/Formula/blasr.rb:36:in `block in <class:Blasr>'
Please report this to the brewsci/science tap!

==> Downloading https://github.com/PacificBiosciences/blasr/archive/smrtanalysis
Already downloaded: /Users/cr517/Library/Caches/Homebrew/blasr-2.2.tar.gz
==> make STATIC= HDF5INCLUDEDIR=/usr/local/opt/hdf5/include HDF5LIBDIR=/usr/loca
Last 15 lines from /Users/cr517/Library/Logs/Homebrew/blasr/01.make:
/Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/include/c++/v1/memory:4912:1: note: candidate template ignored: could not match 'shared_ptr<type-parameter-0-0>' against 'const RegionAnnotation'
operator<(const shared_ptr<_Tp>& __x, const shared_ptr<_Up>& __y) _NOEXCEPT
^
/Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/include/c++/v1/memory:4977:1: note: candidate template ignored: could not match 'shared_ptr<type-parameter-0-0>' against 'const RegionAnnotation'
operator<(const shared_ptr<_Tp>& __x, nullptr_t) _NOEXCEPT
^
/Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/include/c++/v1/memory:4985:1: note: candidate template ignored: could not match 'shared_ptr<type-parameter-0-0>' against 'const int'
operator<(nullptr_t, const shared_ptr<_Tp>& __x) _NOEXCEPT
^
1 error generated.
make[1]: *** [bin/KBandMatcher.o] Error 1
make: *** [alignment-build] Error 2
make: *** [sequtils-build] Error 2
clang++  -O3 -Wno-div-by-zero -I ..//common -I /usr/local/opt/hdf5/include bin/PrintPBBWT.o -o bin/bwtpbPrint -L/usr/local/opt/hdf5/lib -lz -lpthread
clang++  -O3 -Wno-div-by-zero -I ..//common -I /usr/local/opt/hdf5/include bin/BWTQuery.o -o bin/bwtquery -L/usr/local/opt/hdf5/lib -lz -lpthread

If reporting this issue please do so at (not Homebrew/brew or Homebrew/core):
https://github.com/brewsci/homebrew-science/issues

```

I still need to report the issue at the provided address.

I tried Anaconda 'conda install -c bioconda blasr ':

```
PackageNotFoundError: Packages missing in current channels:

  - blasr

We have searched for the packages in the following channels:

  - https://conda.anaconda.org/bioconda/osx-64
  - https://conda.anaconda.org/bioconda/noarch
  - https://conda.anaconda.org/conda-forge/osx-64
  - https://conda.anaconda.org/conda-forge/noarch
  - https://conda.anaconda.org/r/osx-64
  - https://conda.anaconda.org/r/noarch
  - https://repo.continuum.io/pkgs/main/osx-64
  - https://repo.continuum.io/pkgs/main/noarch
  - https://repo.continuum.io/pkgs/free/osx-64
  - https://repo.continuum.io/pkgs/free/noarch
  - https://repo.continuum.io/pkgs/r/osx-64
  - https://repo.continuum.io/pkgs/r/noarch
  - https://repo.continuum.io/pkgs/pro/osx-64
  - https://repo.continuum.io/pkgs/pro/noarch

```

# Answer
> 3 votes
The path of least resistance is probably to run the software inside a container. If you install docker, you can just do the following:

```
git clone https://github.com/muccg/docker-blasr.git
docker-blasr/bin/blasr

```

(this worked for me on OSX 10.11.6)

# Answer
> 2 votes
I've had many issues with installing blasr on my own, but it's a breeze using (bio)conda, see https://anaconda.org/bioconda/blasr

```
conda install -c bioconda blasr

```

# Answer
> 1 votes
is there reason why you want to compile it by your own? You can install blasr though homebrew and brewsci/science tap :

```
brew tap brewsci/science
brew install blasr

```

---
Tags: errors, pacbio, software-installation, blasr
---


title: Problem installing RATT

text:
# Question
Title: Problem installing RATT
I have followed the steps here to install RATT, i.e.:

```
svn co "https://svn.code.sf.net/p/ratt/code/" ratt-code

```

However, when I check my installation by running '$RATT\_HOME/start.ratt.sh', it says:

```
zsh: no such file or directory: /start.ratt.sh

```

which makes sense once I check that the RATT\_HOME variable is empty.

I downloaded their tuberculosis example, went into that directory and ran:

```
~/path/to/my/ratt/ratt-code/start.ratt.sh embl F11.fasta F11 SimStrain

```

and got the following output:

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

```

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

# Answer
> 3 votes
Usually the home of a program is where the executable is. In order to make the system aware you need to include it as a variable when the terminal is run.

So you need to add it to the file at `~/.bashrc` of your computer (if you are using bash as shell \[you can check it with `echo $SHELL`\]). You need to add :

```
export RATT_HOME=path/where/RATT/is/installed/or/downloaded

```

As terdon says in a comment below:

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

---
Tags: software-installation, sequence-annotation, functional-annotation, ratt
---


# Question
Title: What is the difference between the entity_src_genCategory tag and the entity_src_natCategory tag in pdb-xmls?
I want to read out the taxonomy id of the source organism of a pdb entry using the xml-file of the entry that is provided under `http://files.rcsb.org/view/{pdb_accession}-noatom.xml`.

I have so far found the taxonomy id in the tag **PDBx:entity\_src\_genCategory** (example entry 1UKV) but also in the tag **PDBx:entity\_src\_natCategory** (example entry 2BWI).

XML for 1UKV with PDBx:entity\_src\_genCategory

```
<PDBx:entity_src_genCategory>
  <PDBx:entity_src_gen entity_id="1" pdbx_src_id="1">
     <PDBx:gene_src_common_name>baker's yeast</PDBx:gene_src_common_name>
     <PDBx:gene_src_genus>Saccharomyces</PDBx:gene_src_genus>
     <PDBx:host_org_genus>Escherichia</PDBx:host_org_genus>
     <PDBx:host_org_species>Escherichia coli</PDBx:host_org_species>
     <PDBx:pdbx_alt_source_flag>sample</PDBx:pdbx_alt_source_flag>
     <PDBx:pdbx_gene_src_ncbi_taxonomy_id>4932</PDBx:pdbx_gene_src_ncbi_taxonomy_id>
     <PDBx:pdbx_gene_src_scientific_name>Saccharomyces cerevisiae</PDBx:pdbx_gene_src_scientific_name>
     <PDBx:pdbx_host_org_ncbi_taxonomy_id>469008</PDBx:pdbx_host_org_ncbi_taxonomy_id>
     <PDBx:pdbx_host_org_scientific_name>Escherichia coli BL21(DE3)</PDBx:pdbx_host_org_scientific_name>
     <PDBx:pdbx_host_org_strain>BL21(DE3)</PDBx:pdbx_host_org_strain>
     <PDBx:pdbx_host_org_vector_type>PLASMID</PDBx:pdbx_host_org_vector_type>
     <PDBx:plasmid_name>pET19tev</PDBx:plasmid_name>
  </PDBx:entity_src_gen>

```

XML for 2BWI with PDBx:entity\_src\_natCategory

```
<PDBx:entity_src_natCategory>
  <PDBx:entity_src_nat entity_id="1" pdbx_src_id="1">
     <PDBx:pdbx_alt_source_flag>sample</PDBx:pdbx_alt_source_flag>
     <PDBx:pdbx_ncbi_taxonomy_id>223</PDBx:pdbx_ncbi_taxonomy_id>
     <PDBx:pdbx_organism_scientific>ACHROMOBACTER CYCLOCLASTES</PDBx:pdbx_organism_scientific>
  </PDBx:entity_src_nat>

```

What is the difference between them? Are there other tags that contain information about the source organism of the entry? Can they occur together in one file? If so, is one of them more valid?

# Answer
> 1 votes
`src_nat` = isolated directly from a natural tissue

`src_gen` = the source was genetically manipulated

The files that you use (PDBML/XML) are derived from the primary PDB data format which is PDBx/mmCIF. The names in XML seem to correspond to mmCIF names, so you can read documentation on http://mmcif.wwpdb.org/, and here are direct links:

Can they occur together in one file? I suppose they can for complexes: one "entity" could be natural and another genetically manipulated.

---
Tags: pdb, taxonomy
---


# Question
Title: Database of position weight matrices for protein motifs?
I am trying to identify proteins that carry a consensus target sequence for a kinase. Usually if I am working with, for e.g. transcription factors binding to DNA I would use position weight matrices (PWMs) from one of a number of databases of PWMs to search the sequences with something like FIMO.

However, I don't know of any databases that contain protein motif PWMs. Specially I am looking for one for yeast Cdc28, but it seems like it would be a generally useful resource that I can't believe doesn't exist. For example, while the MEME suite was originally developed to identify protein motifs, and ouputs PWMs, all of its databases refer to DNA or RNA.

# Answer
PFam has HMM models for protein domains, this is almost certainly what you're after. You can download all PFam-A models from ftp://ftp.ebi.ac.uk/pub/databases/Pfam/current\_release/Pfam-A.hmm.gz.

Then, you're looking for the "Pkinase", PF00069 family (I found this simply by searching the PFam website for "kinase"), so extract it from the file:

```
zcat Pfam-A.hmm.gz | awk -vRS='//' -vORS='//' '/PF00069/' > PF00069.pfam

```

Then, format it for HMMER format (you will need to install HMMER)

```
hmmpress PF00069.hmm

```

And now you can use it to search your sequences (fasta or flat file, probably others; see the documentation for details):

```
hmmscan PF00069.hmm prots.pep

```

I tried this on a fasta file with the sequences of CBK1\_YEAST, CSK21\_YEAST (kinases) and FLNA\_HUMAN (non-kinase control) and got:

```
# hmmscan :: search sequence(s) against a profile database
# HMMER 3.1b2 (February 2015); http://hmmer.org/
# Copyright (C) 2015 Howard Hughes Medical Institute.
# Freely distributed under the GNU General Public License (GPLv3).
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# query sequence file:             prots.pep
# target HMM database:             PF00069.hmm
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

Query:       sp|P53894|CBK1_YEAST  [L=756]
Description: Serine/threonine-protein kinase CBK1 OS=Saccharomyces cerevisiae (strain ATCC 204508 / S288c) GN=CBK1 PE=1 SV=1
Scores for complete sequence (score includes all domains):
   --- full sequence ---   --- best 1 domain ---    -#dom-
    E-value  score  bias    E-value  score  bias    exp  N  Model    Description
    ------- ------ -----    ------- ------ -----   ---- --  -------- -----------
    2.7e-64  203.2   9.0    3.5e-46  143.8   0.2    3.1  3  Pkinase   Protein kinase domain

Domain annotation for each model (and alignments):
>> Pkinase  Protein kinase domain
   #    score  bias  c-Evalue  i-Evalue hmmfrom  hmm to    alifrom  ali to    envfrom  env to     acc
 ---   ------ ----- --------- --------- ------- -------    ------- -------    ------- -------    ----
   1 ?   -2.9   3.2      0.19      0.19     198     226 ..     222     243 ..     191     277 .. 0.57
   2 !  143.8   0.2   3.5e-46   3.5e-46       1     153 [.     352     505 ..     352     522 .. 0.95
   3 !   63.6   0.0   9.8e-22   9.8e-22     154     264 .]     568     672 ..     551     672 .. 0.83

  Alignments for each domain:
  == domain 1  score: -2.9 bits;  conditional E-value: 0.19
                           TTSSCHHHHHHHHHTHHHHHHHCTSHH.. CS
               Pkinase 198 sgekgkekvekeldqlekilkilgetkek 226
                           + +++++++++      + ++i ++++++
  sp|P53894|CBK1_YEAST 222 Q-QQQQQQQQQ------QHMQIQQQQQQQ 243
                           1.111111111......122222221222 PP

  == domain 2  score: 143.8 bits;  conditional E-value: 3.5e-46
                           EEEEEEEEEESSEEEEEEEETTTTEEEEEEEEEHHCCCHHHH.HHHHHHHHHHHHHSSTTB--EEEEEEESSEEEEEEE--TTEBHHHH CS
               Pkinase   1 yekleklGeGsfGkVykaveketgkivAvKkikkekakkkke.kkvlrEikilkklkhpnivrlyevfeekdelylvleyveggeladl 88
                           +++++ +G+G+fG+V  + +k+tgki+A+K++ k+++ kk++  +v +E  +l+ +++p +v ly++f++ ++lyl++e+++gg+l+++
  sp|P53894|CBK1_YEAST 352 FHTVKVIGKGAFGEVRLVQKKDTGKIYAMKTLLKSEMYKKDQlAHVKAERDVLAGSDSPWVVSLYYSFQDAQYLYLIMEFLPGGDLMTM 440
                           67899********************************999999********************************************** PP

                           HHHHSS--HHHHHHHHHHHHHHHHHHHHTTEE-SS-SGGGEEEETTTEEEE-SGTTSEECCSSCS CS
               Pkinase  89 lsrkgslseeeakkiakqilegleylHskgiiHrDlKpeNiLidekgelKitDFGlakelesssk 153
                           l r + ++e+ +++++++ + ++e +H+ g+iHrD+Kp+NiLid +g++K++DFGl++ ++++++
  sp|P53894|CBK1_YEAST 441 LIRWQLFTEDVTRFYMAECILAIETIHKLGFIHRDIKPDNILIDIRGHIKLSDFGLSTGFHKTHD 505
                           **999*****************************************************9999876 PP

  == domain 3  score: 63.6 bits;  conditional E-value: 9.8e-22
                           BCTCCSCGGGS-HHHHTTSTCSHHHHHHHHHHHHHHHHHSS-STTTSSCHHHHHHHHHTHHHHHHHCTSHH......TTS-HHHHHHHH CS
               Pkinase 154 ltsfvgtreYlAPEvlkeneyskkvDvWslGvilyelltgkppfsgekgkekvekeldqlekilkilgetkeklpeaselseeakdllk 242
                           + s vgt++Y+APE+  +++y++++D+WslG i+ye l g ppf +e+ +e+++k+++          e+++++p + ++s ea dl++
  sp|P53894|CBK1_YEAST 568 AYSTVGTPDYIAPEIFLYQGYGQECDWWSLGAIMYECLIGWPPFCSETPQETYRKIMNF---------EQTLQFPDDIHISYEAEDLIR 647
                           6789*****************************************66665533333333.........557799*************** PP

                           HHT-SSGCCSTT....HHHHHTSGGG CS
               Pkinase 243 kllkkdpkkRlt....aeellqhpyl 264
                           +ll+    +Rl+    a+e+++hp++
  sp|P53894|CBK1_YEAST 648 RLLTHA-DQRLGrhggADEIKSHPFF 672
                           ***875.678877779********98 PP

Internal pipeline statistics summary:
-------------------------------------
Query sequence(s):                         1  (756 residues searched)
Target model(s):                           1  (264 nodes)
Passed MSV filter:                         1  (1); expected 0.0 (0.02)
Passed bias filter:                        1  (1); expected 0.0 (0.02)
Passed Vit filter:                         1  (1); expected 0.0 (0.001)
Passed Fwd filter:                         1  (1); expected 0.0 (1e-05)
Initial search space (Z):                  1  [actual number of targets]
Domain search space  (domZ):               1  [number of targets reported over threshold]
# CPU time: 0.01u 0.00s 00:00:00.01 Elapsed: 00:00:00.00
# Mc/sec: inf
//
Query:       sp|P15790|CSK21_YEAST  [L=372]
Description: Casein kinase II subunit alpha OS=Saccharomyces cerevisiae (strain ATCC 204508 / S288c) GN=CKA1 PE=1 SV=1
Scores for complete sequence (score includes all domains):
   --- full sequence ---   --- best 1 domain ---    -#dom-
    E-value  score  bias    E-value  score  bias    exp  N  Model    Description
    ------- ------ -----    ------- ------ -----   ---- --  -------- -----------
    7.5e-62  195.1   7.2    3.2e-55  173.4   0.7    2.1  2  Pkinase   Protein kinase domain

Domain annotation for each model (and alignments):
>> Pkinase  Protein kinase domain
   #    score  bias  c-Evalue  i-Evalue hmmfrom  hmm to    alifrom  ali to    envfrom  env to     acc
 ---   ------ ----- --------- --------- ------- -------    ------- -------    ------- -------    ----
   1 !   25.5   1.1   4.3e-10   4.3e-10       1      60 [.      40      94 ..      40     115 .. 0.87
   2 !  173.4   0.7   3.2e-55   3.2e-55      56     264 .]     128     363 ..     101     363 .. 0.92

  Alignments for each domain:
  == domain 1  score: 25.5 bits;  conditional E-value: 4.3e-10
                           EEEEEEEEEESSEEEEEEEETTTTEEEEEEEEEHHCCCHHHHHHHHHHHHHHHHHSSTTB CS
                Pkinase  1 yekleklGeGsfGkVykaveketgkivAvKkikkekakkkkekkvlrEikilkklkhpni 60
                           ye ++k+G+G++++V+++v+ +++ ++++K++k  k+kk k     rEikil  l+ +++
  sp|P15790|CSK21_YEAST 40 YEIENKVGRGKYSEVFQGVKLDSKVKIVIKMLKPVKKKKIK-----REIKILTDLSNEKV 94
                           78999****************************99888855.....7******9998877 PP

  == domain 2  score: 173.4 bits;  conditional E-value: 3.2e-55
                            S.STTB--EEEEEEE..SSEEEEEEE--TTEBHHHHHHHHSS--HHHHHHHHHHHHHHHHHHHHTTEE-SS-SGGGEEEE.TTTEEEE CS
                Pkinase  56 k.hpnivrlyevfee..kdelylvleyveggeladllsrkgslseeeakkiakqilegleylHskgiiHrDlKpeNiLid.ekgelKi 139
                            + h ni++l++++++  +++  lv+eyv++ ++  l+    +l++ e+++++ ++l++l+y+Hs+gi+HrD+Kp+N++id ++++l++
  sp|P15790|CSK21_YEAST 128 NgHANIIHLFDIIKDpiSKTPALVFEYVDNVDFRILYP---KLTDLEIRFYMFELLKALDYCHSMGIMHRDVKPHNVMIDhKNKKLRL 212
                            35*************8878888*********9988887...5**************************************66677*** PP

                            -SGTTSEECCSSCSBCTCCSCGGGS-HHHH.TTSTCSHHHHHHHHHHHHHHHHHSS-STTTSSCHHHHHHHHHTHHHHHHHCTSHH.. CS
                Pkinase 140 tDFGlakelessskltsfvgtreYlAPEvl.keneyskkvDvWslGvilyelltgkppfsgekgkekvekeldqlekilkilgetkek 226
                            +D+Gla+ ++ + +++ +v++r +  PE+l ++++y+ + D+Ws G +l++++++++pf   +g++    ++dql+ki+k+lg+++++
  sp|P15790|CSK21_YEAST 213 IDWGLAEFYHVNMEYNVRVASRFFKGPELLvDYRMYDYSLDLWSFGTMLASMIFKREPFF--HGTS----NTDQLVKIVKVLGTSDFE 294
                            ******************************999***************************..3443....789**********77766 PP

                            ....TTS...............................-HHHHHHHHHHT-SSGCCSTTHHHHHTSGGG CS
                Pkinase 227 lpeasel...............................seeakdllkkllkkdpkkRltaeellqhpyl 264
                                +                                 ++e++dl+++ll++d+++Rlta+e++ hp++
  sp|P15790|CSK21_YEAST 295 KYLLKYEitlprefydmdqyirkpwhrfindgnkhlsgNDEIIDLIDNLLRYDHQERLTAKEAMGHPWF 363
                            5554444588999999****************************************************9 PP

Internal pipeline statistics summary:
-------------------------------------
Query sequence(s):                         1  (372 residues searched)
Target model(s):                           1  (264 nodes)
Passed MSV filter:                         1  (1); expected 0.0 (0.02)
Passed bias filter:                        1  (1); expected 0.0 (0.02)
Passed Vit filter:                         1  (1); expected 0.0 (0.001)
Passed Fwd filter:                         1  (1); expected 0.0 (1e-05)
Initial search space (Z):                  1  [actual number of targets]
Domain search space  (domZ):               1  [number of targets reported over threshold]
# CPU time: 0.01u 0.00s 00:00:00.01 Elapsed: 00:00:00.02
# Mc/sec: 4.91
//
Query:       sp|P21333|FLNA_HUMAN  [L=2647]
Description: Filamin-A OS=Homo sapiens GN=FLNA PE=1 SV=4
Scores for complete sequence (score includes all domains):
   --- full sequence ---   --- best 1 domain ---    -#dom-
    E-value  score  bias    E-value  score  bias    exp  N  Model    Description
    ------- ------ -----    ------- ------ -----   ---- --  -------- -----------

   [No hits detected that satisfy reporting thresholds]

Domain annotation for each model (and alignments):

   [No targets detected that satisfy reporting thresholds]

Internal pipeline statistics summary:
-------------------------------------
Query sequence(s):                         1  (2647 residues searched)
Target model(s):                           1  (264 nodes)
Passed MSV filter:                         0  (0); expected 0.0 (0.02)
Passed bias filter:                        0  (0); expected 0.0 (0.02)
Passed Vit filter:                         0  (0); expected 0.0 (0.001)
Passed Fwd filter:                         0  (0); expected 0.0 (1e-05)
Initial search space (Z):                  1  [actual number of targets]
Domain search space  (domZ):               0  [number of targets reported over threshold]
# CPU time: 0.00u 0.00s 00:00:00.00 Elapsed: 00:00:00.00
# Mc/sec: inf
//

```

---

HMMER is a great tool and can also analyse multiple sequence alignments to build the hmms if that's what you need. I strongly urge you to reads its documentation which is both informative and a surpisingly fun read (it includes such phrases as "*Their format is “proprietary”, which is an open source term of art that means both “I haven’t found time to document them yet” and “I still might decide to change them arbitrarily without telling you”*").
> 2 votes

---
Tags: public-databases, motifs
---


# Question
Title: Different outputs of computeMatrix depending on binsize
I'd like to output a figure in +/-1kb around a center using computeMatrix as following command:

```
computeMatrix reference-point -S data.bw -R centers.bed -a 1000 -b 1000 -bs=10 -out newdata10.gz --referencePoint center

```

The bed file is like:

```
chr1 100 101

```

If I change `-bs=1`, two outputs of the same input newdata10.gz are quite different. Which's correct with my commands? Maybe I don't understand the meaning of binsize option.

# Answer
Please don't send things like this to the deepTools mailing list and post them here.

As I mentioned on the list, what you're likely being surprised by is that computeMatrix is not summing values, but taking the average in each bin. You can change this with the `--averageTypeBins` option (the default is `mean`, you presumably want `sum`).

**Update**: For those wondering the cause of the confusion was due to not using `--missingDataAsZero`. In other words, including or not including NAs will lead to different results, since different numbers of values will be included.
> 1 votes

---
Tags: deeptools
---


# Question
Title: Downloading an open-source MAF with tumor and normal BAMs?
I'm writing an R package that manipulates MAF files. Normally, I generate a small example of a file and use this for testing. (e.g. a subsetted BAM or MAF).

In this case, I would need a MAF and a tumor/normal BAM pair. My problem is I don't know where I could download an open-source MAF with said BAMs.

Where could researchers download files like this, and not violate privacy laws (i.e. they *must* be publicly-available BAMs)?

# Answer
There's the NCI Genomic Data Commons:

> For all submitted sequence data, including BAM alignment files, the GDC generates new alignments in BAM format using the latest human reference genome GRCh38 with standard alignment pipelines. Using these standard alignments, the GDC generates high level derived data, including normal and tumor variant and mutation calls in VCF and MAF formats, and gene and miRNA expression and splice junction quantification data in TSV formats.

But that needs an NIH account for individual-level data, which might be too restrictive for you.

There's also the Texas Cancer Research Biobank (7 human cancer patients with matched normal):

> These tumor and normal specimen sequence files (FASTQ and BAM), somatic variant calls (VCF and MAF) and germline MAF files are available for each consented patient as described on the TCRB Open Access Privacy page.

That also requires account registration, but the registration process is automated and much simpler (requires just a name, and acceptance of usage conditions). There's a published paper about that dataset.

I had some vague feeling that another repository existed (created within the last year), but I can't seem to find it from quick web searches.
> 1 votes

---
Tags: sam, public-databases, maf
---


title: What are the count units in DEseq2?

text:
# Question
Title: What are the count units in DEseq2?
This is a pretty silly/simple question, I suppose.

What units are DESeq2 counts? Or are the units arbitrary but internally consistent estimate of actual reads?

# Answer
> 7 votes
The counts are "reads" for single-end datasets and "fragments" for paired-end datasets. In other words, they're what featureCounts and htseq-count produce. The "normalized counts" that you'll be able to access in DESeq2 are, as aptly named, simply normalized versions such that values are comparable across samples.

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


# Question
Title: How to download the whole directory of an ensembl FTP page?
I am interested in downloading all the files at:

ftp://ftp.ensembl.org/pub/release-91/embl/bos\_taurus/

I know how to download a single file from that page using wget (Unix) or download.file() (R) but is there an easy way of downloading all of them?

# Answer
> 5 votes
You can use the mirroring function of wget:

```
wget -m ftp://ftp.ensembl.org/pub/release-91/embl/bos_taurus/

```

---
Tags: ensembl, data-download
---


# Question
Title: How to read the "SEQRES" section from a PDB file, using R
I'm using R and its Rpdb package to work with PDB files. I use the `read.pdb` function but, nevertheless, I couldn't find the argument which will allow me to view the `SEQRES` section and the full amino-acid sequence of the protein.

For example, I want to find which amino-acid is in the 17th position of the sequence and I want to extract it from the 'SRQRES' section in the PDB file (the original length of the sequence was reduced for the sake of the example):

```
SEQRES   1 A  390  GLU PRO GLU TRP THR TYR PRO ARG LEU SER CYS GLN GLY
SEQRES   2 A  390  SER THR PHE GLN LYS ALA LEU LEU ILE SER PRO HIS ARG
SEQRES   3 A  390  PHE GLY GLU ALA ARG GLY ASN SER ALA PRO LEU ILE ILE
SEQRES   4 A  390  ARG GLU PRO PHE ILE ALA CYS GLY PRO LYS GLU CYS LYS
SEQRES   5 A  390  HIS PHE ALA LEU THR HIS TYR ALA ALA GLN PRO GLY GLY
SEQRES   6 A  390  TYR TYR ASN GLY THR ARG GLU ASP ARG ASN LYS LEU ARG
SEQRES   7 A  390  HIS LEU ILE SER VAL LYS LEU GLY LYS ILE PRO THR VAL
SEQRES   8 A  390  GLU ASN SER ILE PHE HIS MET ALA ALA TRP SER GLY SER
SEQRES   9 A  390  ALA CYS HIS ASP GLY ARG GLU TRP THR TYR ILE GLY VAL
SEQRES  10 A  390  ASP GLY PRO ASP SER ASN ALA LEU ILE LYS ILE LYS TYR

```

...

The original file can be found here 1a4g.pdb.

Any ideas on how to accomplish that?

# Answer
I have not used `Rpdb` but I have used, and would recommend, bio3d.

To read the PDB file 1A4G:

```
library(bio3d)
pdb1a4g <- read.pdb("1a4g")

```

To get the SEQRES information, there are two options. Get it from the `pdb` object directly:

```
seq1a4g <- pdb1a4g$seqres

```

Or call `pdbseq`:

```
seq1a4g <- pdbseq(pdb1a4g)

```

Using the second method, the 17th position (note this is not the same as the 17th position of the original protein sequence):

```
seq1a4g[17]
 92
"Q"

```
> 4 votes

---
Tags: r, pdb
---


# Question
Title: factoextra: Error in if (xlab == "x") xlab <- "x value" : argument is of length zero
I just followed the tutorial trying to do clustering on my own data, but I can't visualize it because `fviz_cluster` function is giving me an error:

> Error in if (xlab == "x") xlab \<- "x value" : argument is of length zero

Here is the code I am running:

```
library("fpc")
db <- fpc::dbscan(rtsne$Y, eps = 0.15, MinPts = 5)
# Plot DBSCAN results
library("factoextra")
fviz_cluster(db, data = rtsne$Y, stand = FALSE,
        ellipse = FALSE, show.clust.cent = FALSE,
         geom = "point",palette = "jco", ggtheme = theme_classic())

```

How can I make `fviz_cluster` function working?

> Update

Filed a bug: https://github.com/kassambara/factoextra/issues/62

# Answer
I actually found an answer just accidentally. It is very unfortunate that `factoextra` documentation does not explicitly say that `data` parameter should be of `data.frame` type. In my case `rtsne$Y` was a `matrix`. The error that `factoextra` is giving is very uninformative.

Here is the code that ultimately worked:

```
fviz_cluster(db,
----->  data = as.data.frame(rtsne$Y),
        stand = FALSE,
        ellipse = FALSE, show.clust.cent = FALSE,
        geom = "point",palette = "jco", ggtheme = theme_classic())

```
> 3 votes

---
Tags: clustering
---


# Question
Title: Detect transcript isoform abundance for a specific gene in scRNA-seq
I want to detect the count of isoform transcripts for a specific gene in scRNA-seq data. Data is coming from cells of Mus Musculus. For transcript isoforms I mean the different alternatives provided by known annotation. For each one I need the total transcript count in each cell.

By looking at the count matrix, alternative transcripts are not defined (see for instance this dataset).

Should I come back to the fastq? Which tool may I use to produce a count matrix with counts per each alternative transcript?

# Answer
Transcript quantification is a difficult enough problem as it is, when you add the extra difficulty of going from the low read numbers available in scRNAseq if gets even more difficult.

Added to this, many scRNAseq techniques (including 10X, which is what the linked dataset is) only capture the 3' end of the gene, so you can distinguish transcript with different 3' end, but the data simply doesn't exist for doing much more with standard tools.

Lior Patchner does claims that there is benefit in looking at "transcript compatibility counts" rather than just gene counts, and that you can use his Kallisto tool for this, and the these 3' end tagging tools do give you some data for the rest of the transcript.

You could also isolate transcripts with differing three 3' ends and map reads to these using the something like the UMI-Tools pipeline.

Alternatively you could look at things on a junction level - you won't capture all junctions, but you may be able to look at differential junction usage for those you do capture using something like SUPPA or MISO. For this I would use UMI-Tools to get de-duplicated BAM files, split them on cell identity and then pass them into your tool of choice. I don't really how much you will capture.

The last thing to say is that this might be a space worth watching for future developments (*cough*).

Data from some other single cell methods, like SMART-seq2 or C1 will sequence the whole transcript and you may be able to apply standard tools such as Salmon, Kalisto or RSEM to get transcript usage for them, although as I said, transcript level quantification is noisey, single cell data is noisey, single-cell level transcript data will be very noisy.

Whatever, you are going to have to go back to the raw sequencing data.

Disclosure: I am the author of UMI-Tools.
> 4 votes

# Answer
Looking for alternative splicing in scRNA-seq is difficult, but there is at least one tool that is specifically designed for it: DISCO. The site does not have too much info, but there is a presentation about it. It uses MISO output as the input, which in turn uses SAM.
> 3 votes

---
Tags: rna-seq, scrnaseq, rna-splicing, isoform
---


# Question
Title: Aligning ChIP-Seq reads to repeats for downstream peak analysis
This is a brief question regarding the above. I have previously used bowtie to map reads from paired-end ChIP-Seq sequencing, then used the positions for peak-calling. However, I'm trying to do similar with a dataset which may consist of repeats (I have previously aligned the dataset using bowtie, and a lot of reads aligned to repeats).

For peak-calling purposes, however, I need to know reasonably confidently where precisely those repeats are located, so that MACS2 treats the coverage pileup sensibly. So my question is this: what's the best way of mapping these (short) paired-end reads from a ChIP-Seq experiment so I can do this?

EDIT: I want to be able to run a peak-caller after aligning the data, so I want a good approximation to the coverage at repeating elements after mapping for instance - I don't want to throw away all reads which multimap, as may mean that some reads could lose their mates, and so I won't be able to calculate coverage or call peaks using the data.

I want an aligner to keep both ends of the paired-end pair within about 180bp of each other (as that is the expected length of the DNA fragment that is sequenced). The immediate step after aligning the pairs is to flatten the pair into the DNA fragment that has been sequenced. With multimapped reads this step could be done by seeing which of the pair is uniquely mapped and choosing the mate which maps closest to it.

I have tried bowtie (and bowtie2) for this, which is where I understood that the dataset maps to repeats.

# Answer
> 2 votes
This is very good question, but as far as I know there really isn't a good answer. I will attempt to provide some comments and tips.

> However, I'm trying to do similar with a dataset which may consist of repeats (I have previously aligned the dataset using bowtie, and a lot of reads aligned to repeats).

If this dataset is also PE and the alignment is outputting either uniquely mapped reads or a single (random) position for multimappers, it is pretty much as good as it gets. A lot of repeat regions are not **exact** copies, so using PE there is a good change that one read of the pair mapping to a single position is helping with mapping of the other read.

> I need to know reasonably confidently where precisely those repeats are located

It is probably a matter of phrasing, but what one needs is to know the precise genomic location of the **reads**. The location of the repeat elements can be obtained from the Repeat Masker track of UCSC, but you probably know this already.

> With multimapped reads this step could be done by seeing which of the pair is uniquely mapped and choosing the mate which maps closest to it.

I **think** `bowtie` already does this. Did you map the read pairs together or separately?

Note that:

* There will always be a number of multimapping reads (fragments) that will be discarded. The information of how many are being discarded for your samples would be helpful.
* Mapping can be improved using both reads and tunning `bowtie` settings (see bellow)

> Finding a valid paired-end alignment where both mates align to repetitive regions of the reference can be very time-consuming. By default, Bowtie avoids much of this cost by imposing a limit on the number of "tries" it makes to match an alignment for one mate with a nearby alignment for the other. The default limit is 100. This causes bowtie to miss some valid paired-end alignments where both mates lie in repetitive regions, but the user may use the --pairtries or -y options to increase Bowtie's sensitivity as desired. source
>
> The immediate step after aligning the pairs is to flatten the pair into the DNA fragment that has been sequenced.

As indicated above, the aligner will take care of that sort of. It will gave the mapping information of each read pair, which is then used by `MACS`.

**Other tools**

I tested `ShortStack` which uses the coverage of uniquely mapped reads to place multimapping reads for ChIP-seq, and the results were mixed. That said, it was SE data and the code has been through some changes, so it be worth a try.

When requested, the authors of `RepEnrich` provided me with the scripts for an utility called `RepConsensus` which maps reads to a **consensus** reference of repeat regions. The caveats are:

* mappings are better for elements whose sequences have not diverged a lot yet
* it will tell the position of the protein in the repeat sequence, but no information about the surrounding genomic regions because there isn't any.

---

Whatever you do, do not use all alignments, that is double counting of reads, and think very carefully before using custom methods. A blog post detailing some of the potential issues can be found here (I wrote it).

---
Tags: sequence-alignment, chip-seq, repeat-elements
---


title: How to simulate biological pathways?

text:
# Question
Title: How to simulate biological pathways?
I would like to simulate a biological pathway in Matlab, Mathematica or Python. To study the response of the pathway using Control Theory. Is there a standard way to simulate pathways by using differential equations or any other formalism.

For example I would like to try to simulate this KEGG pathway map. Is there something more recent than this work from 2005?

For example to simulate Petri networks?

# Answer
> 7 votes
In matlab this can be done using the SimBiology toolbox.
In mathematica this can be done using the SystemModeler.
In python there are multiple packages, e.g. (in no particular order) PySB, Tellurium, PySD

---
Tags: networks, pathway, systems-biology
---


# Question
Title: Counting the number of paralogues for mouse genes gives me the wrong frequency in R
I am trying to count the number of paralogues for the mouse homologues of the human protein-coding genes using BioMart. But for example in the 'PLIN4' gene its counting 35,000 paralogues instead of 4.

We think it is because some genes have one to many paralogues which causes repeats. When I run a single gene its gives me back the correct number of paralogues. Is there a way to either remove these repeats from the results or a way around this so that BioMart doesn't output these repeats.

I have also thought of maybe running one gene at a time, then counting it by setting up some sort of loop so that it does all of the genes from the list automatically.

The code I have written so far is:

```
# Load the biomaRt package:

library(biomaRt)
ensembl_hsapiens <- useMart("ensembl",
                          dataset = "hsapiens_gene_ensembl")
ensembl_mouse <- useMart("ensembl",
                       dataset = "mmusculus_gene_ensembl")

# Get all human protein coding genes:

hsapien_PC_genes <- getBM(attributes = c("ensembl_gene_id",
                                         "external_gene_name"),
                          filters = "biotype",
                          values = "protein_coding",
                          mart = ensembl_hsapiens)

ensembl_gene_ID <- hsapien_PC_genes$ensembl_gene_id

# Get mouse homologues

mouse_homologues <- getBM(attributes = c("ensembl_gene_id", "external_gene_name",
                                       "mmusculus_homolog_associated_gene_name"),
                        filters = "ensembl_gene_id",
                        values = c(ensembl_gene_ID),
                        mart = ensembl_hsapiens)

# Get mouse external gene name

mouse_homologues_external_gene_names <- mouse_homologues$mmusculus_homolog_associated_gene_name

mouse_paralogues <- getBM(attributes = c("hsapiens_homolog_associated_gene_name",
                                       "external_gene_name",
                                       "mmusculus_paralog_associated_gene_name"),
                        filters = "external_gene_name",
                        values = c(mouse_homologues_external_gene_names) , mart = ensembl_mouse)

# Remove genes with no paralogues
mouse_paralogs_data <- mouse_paralogues[!(is.na(mouse_paralogues$mmusculus_paralog_associated_gene_name)
                                          |
mouse_paralogues$mmusculus_paralog_associated_gene_name==""), ]

# Count paralogues per gene

library(plyr)
count_mouse_paralogues <- count(mouse_paralogs_data, "external_gene_name")
View(count_mouse_paralogues)

```

# Answer
> 1 votes
I'd use `dplyr` for this.

First, `mouse paralogues` contains no values for `mmusculus_paralog_associated_gene_name` that are NA:

```
library(dplyr)
mouse_paralogues %>%
  filter(is.na(mmusculus_paralog_associated_gene_name)) %>%
  nrow()

[1] 0

```

But there are values which are empty strings, so you can filter for those if you wish.

```
count_mouse_paralogues <- mouse_paralogues %>%
  filter(hsapiens_homolog_associated_gene_name != "") %>%
  count(external_gene_name)

```

Confirm that this gives the expected result (n = 4) for Plin4:

```
count_mouse_paralogues %>%
  filter(external_gene_name == "Plin4")

# A tibble: 1 x 2
  hsapiens_homolog_associated_gene_name     n
  <chr>                                 <int>
1 Plin4                                     4

```

And check `mouse_paralogues`:

```
mouse_paralogues %>%
  filter(hsapiens_homolog_associated_gene_name == "PLIN4")

  hsapiens_homolog_associated_gene_name external_gene_name mmusculus_paralog_associated_gene_name
1                                 PLIN4              Plin4                                  Plin3
2                                 PLIN4              Plin4                                  Plin5
3                                 PLIN4              Plin4                                  Plin2
4                                 PLIN4              Plin4                                  Plin1

```

---
Tags: r, bioconductor, biomart, ensembl
---


# Question
Title: Comparing two multi-fasta files of the same set of proteins with parser - to find and count mutations after treatment
My task is to count the mutations occurred in several proteins after a treatment. The sequences are all present in the two files in the same order. I opened both files with the FASTA parser (`SeqIO.parse`) in Biopython and I got all the proteins listed (separated before and after the treatment).

How can I zip the parsers together to count the mutations?

How can I count the mutations that occured after the treatment?

```
from Bio import SeqIO
for normal_samples in SeqIO.parse("/data/statistic/normal_samples", "fasta"):
    print(normal_samples.id)
    print(repr(normal_samples.seq))
    print(len(normal_samples))

for treated_samples in SeqIO.parse("/data/statistic/with_treatment", "fasta"):
    print(normal_samples.id)
    print(repr(normal_samples.seq))
    print(len(normal_samples))

dict_n_t = dict(zip(normal_samples & treated_samples))

```

# Answer
```
from Bio import SeqIO
for normal, treated in zip(SeqIO.parse("/data/statistic/normal_samples", "fasta"),
                           SeqIO.parse("/data/statistic/with_treatment", "fasta")):
    ... do stuff...

```

That's generally how you `zip` iterators together in python.
> 3 votes

---
Tags: python, fasta, proteins, statistics, mutations
---


# Question
Title: Why does RATT create a directory with chromosome sequences and how can I stop it?
I have run RATT to transfer the annotation from a reference strain to another strain. For this, I have run:

```
$RATT_HOME/start.ratt.sh embl ref.fa output/annotation/ Strain

```

This creates several things in the working directory that I do not wish:

* a Sequences folder with the sequence in FASTA format of each chromosome
* a copy of the reference sequence in FASTA format (followed by a mysterious number)
* a copy of the query sequence in FASTA format followed by a mysterious number

I could not find any options in the RATT help to disable the creation of these files. Also, does anybody know why these files are created?

# Answer
RATT seems to be a very simple wrapper script that unfortunately does not clean up its temporary files properly. Short of modifying the code yourself, your only option is to perform the cleanup manually after running RATT.

The temporary copies are made because (a) the query is not an exact copy: RATT performs some cleaning. And (b) it needs to extract a FASTA file from the EMBL input for use in downstream tools.

The “mysterious” number behind the file names is simply the process ID of the RATT process. This is appended to make these intermediate files more or less unique (presumably to prevent overwriting existing files).
> 1 votes

---
Tags: sequence-alignment, sequence-annotation, functional-annotation, ratt
---


title: COSMIC Genotypes and Phenotypes

text:
# Question
Title: COSMIC Genotypes and Phenotypes
I'm trying to find files containing genotypes and phenotypes for cell lines in COSMIC. On their downloads page I found links for mutation data, which I guess could be converted to genotypes, but I haven't been able to find any phenotype information. I'm specifically interested in drug sensitivity information. Does anyone know if it's possible to get more user friendly genotype files and whether phenotype data exists publicly?

# Answer
It seems like I was wrong and COSMIC does store information about cell lines.

You need to download the files from Genotypes:

> Genotypes
>
> Files listing the SNP calls for each cell line identified by PICNIC analysis of Affymetrix SNP6.0 array data. Both a simple genotype (AA, BB - homozygous or AB - heterozygous) and a complex interpretation of the genotype are given (for example, in a triploid region of the genome the genotype maybe AAB). 'genotypes/.'

You need to register (and get a license) to be able to download the data. Presumably the cell lines are described in the heading of the file. Then you can sort them out by cell line.

Then once you know the germ line identifier look for the paper describing it to know more about their phenotype. (But maybe there is a database with phenotype of known cell lines)
> 3 votes

---
Tags: ngs, cancer
---


# Question
Title: SNP located within a promoter region (pig)
I have a couple of SNP identifiers such as MARC0073381 or ALGA0066960. The corresponding platform is Illumina Porcine SNP60 BeadChip (WG-410).

I want to know if these SNP are located within a promoter region of a gene.

Bioconductor offers information for Affymetrix chips, but I am not able to map the identifiers.

How can I map SNPs from the beadchip to the porcine genome?

# Answer
> 5 votes
As far as I'm aware, Illumina provide CSV annotation files for all their sequencing chips, which can be used when they can't be found in Bioconductor. You can find annotation information for the PorcineSNP60 here, in particular the Manifest file (CSV format). The format is Illumina's weird "we say it's a CSV because there are commas in it" format, so if there are only a few locations of interest it would be easier to use grep than try to work out the right way to load it into R (or similar).

For completeness, you want the lines between `[Assay]` and `[Controls]`. Here's an example way to pre-process the "CSV" file into a more standard format that just has one data table in it:

```
$ grep -n '^\[' PorcineSNP60v2_15031945_C1.csv
2:[Heading]
7:[Assay]
61574:[Controls]
### See that the table starts at line 8 (after [Assay]),
###   and finishes at line 61753
$ tail -n +8 PorcineSNP60v2_15031945_C1.csv | head -n $((61574-8)) > cleaned_PorcineSNP60v2_15031945_C1.csv

```

But for a few markers and grep, it doesn't matter if you're working with the clean or dirty CSV file:

```
$ grep -e '^IlmnID' -e '^MARC0073381' -e '^ALGA0066960' PorcineSNP60v2_15031945_C1.csv
IlmnID,Name,IlmnStrand,SNP,AddressA_ID,AlleleA_ProbeSeq,AddressB_ID,AlleleB_ProbeSeq,GenomeBuild,Chr,MapInfo,Ploidy,Species,Source,SourceVersion,SourceStrand,SourceSeq,TopGenomicSeq,BeadSetID,Exp_Clusters,Intensity_Only
ALGA0066960-2_B_F_2199219285,ALGA0066960,BOT,[T/G],0054749354,TGCCACAGGTCTGCTCAGCTCAAGCCCAACACTCGCAAGATACAGGTCTA,,,10.2,12,55218591,diploid,Sus scrofa,PorcineSNP60v2,2,BOT,GGCTCACCTCTGCCACAGGTCTGCTCAGCTCAAGCCCAACACTCGCAAGATACAGGTCTA[T/G]CCAGGCTTCTCTCTCCTACTCCTAGGGCCCCTGATGGTTCCTGCATCCTGACCAATAGTG,CACTATTGGTCAGGATGCAGGAACCATCAGGGGCCCTAGGAGTAGGAGAGAGAAGCCTGG[A/C]TAGACCTGTATCTTGCGAGTGTTGGGCTTGAGCTGAGCAGACCTGTGGCAGAGGTGAGCC,670,3,0
MARC0073381-2_T_F_2199252059,MARC0073381,TOP,[A/G],0047761415,TGGAACGGATGGTGGAGACATTCTGGAGACAGAAGACAAACTGCTTCAGA,,,10.2,7,61808556,diploid,Sus scrofa,PorcineSNP60v2,2,TOP,CAACTGTGGTTGGAACGGATGGTGGAGACATTCTGGAGACAGAAGACAAACTGCTTCAGA[A/G]CAAAGCTCAGGAAGGCAA,CAACTGTGGTTGGAACGGATGGTGGAGACATTCTGGAGACAGAAGACAAACTGCTTCAGA[A/G]CAAAGCTCAGGAAGGCAA,660,3,0

```

I can see from this that the first marker is Chromosome 12, position 55218591, and the second marker is Chromosome 7, position 61808556. The Illumina annotation doesn't include the nearest gene, so you'll need to hunt for those locations on a genome browser (or using mapping with Bioconductor).

---
Tags: r, bioconductor, snp, illumina, porcine
---


# Question
Title: How hard is it to clean and QC gene expression microarray data?
I am a PhD student working on developing ideas for my dissertation papers. One of my planned papers will be working with some (human) gene expression data.

I have had a class on working with gene expression data, but no practical experience. From what I remember from my project in that class, you do a bunch of pre-checks and QC, then you remove batch effects, then you summarize expression across probes and normalize. I thank that part of my class project took me about three days.

But when I talk to other people about it, they say that cleaning/QCing/normalizing gene expression data is really hard, that it will take forever to get it working, and that if I do this dissertation paper I am going to end up taking extra years to graduate.

Am I missing something here? Why do the people I talk to seem to think this is so hard and time consuming -- is it really that bad? Should I avoid having a dissertation paper that uses gene expression data if I want to graduate on time?

(The data will be from Affymetrix microarrays. There are enough samples to use two 96-well plates, so batch effects will be more of a concern than they were in my class project.)

# Answer
I am not sure how much you know about bioinformatics already, can you use R? For a bioinformatician looking at QC for microarrays should not be a big deal, at least for me it would take maybe a day (or two) to get this done. However, if you never used R and want to start from scratch, it depends on how quickly you learn how to deal with arrays and QC. There are many tutorials and review papers about normalization and QC of arrays. Try there and see if you can easily understand what these papers are telling you.

Read the limma user's guide, and try to follow the tutorials therein. If you understand that quickly you can try a real experiment of your own. If you don't understand what's described in there, you might want to look for another project for your thesis.
> 2 votes

---
Tags: microarray, quality-control, data-preprocessing, batch-effects
---


# Question
Title: Is there a standard method to validate whether a bigWig file is correct?
I have a bigWig file of extension `*.bw` (which is standard), but I am suspicious that this file may be of an incorrect format.

Is there a standard method in the field to validate a \*.bigWig file?

# Answer
> 1 votes
As noted in another answer, it's relatively easy to tell if a file is bigWig or not by opening it up with, say, standard UCSC Kent utilities like `bigWigToBedGraph` or `bigWigToWig`, or other bigWig readers, and seeing if you get data back.

Depending on what you're doing with the bigWig file, the trickier part can be how the original data was indexed.

If the bigWig file was originally made from a bedGraph or BAM file, its coordinates will be 0-indexed when converted to human-readable form. If the bigWig file was originally made from a wig file, its coordinates will be 1-indexed when converted to a human readable form. How the file's coordinates were originally indexed may change how you work with them downstream for set operations, for instance.

I think the only way to know about the source is to look at the header contained within the bigWig file, after extraction to human-readable form with UCSC tools. A bigWig file made from a bedGraph file will generally contain the header `#bedGraph section <coords>` in it somewhere (after conversion), for example.

# Answer
> 2 votes
There's no standard validation method, but a simple one would be `python -c "import pyBigWig; bw = pyBigWig.open('your_file.bw'); print(bw.IsBigWig());"`. If you get `True` then it's a bigWig file. If you get `False` or an error then it's not (`False` indicates a bigBed file).

Alternatively, `file your_file.bw` should return `data`.

---
Tags: bigwig, validation
---


# Question
Title: Parsing SAM/BAM files for additional information
I used BWA-MEM to alignment and I would like to gather the some informations like total % of match, mismatch, insert/delete etc. I am wondering if there is any existing tools that produces this information?

# Answer
> 2 votes
Other options... samtools stats maybe gives you some of the stats you are after. Chances are you already have samtools so nothing additional to install. Also useful may be pysamstats, in particular the `variation` sub-command.

# Answer
> 0 votes
I've written my own script, readstomper, which gives per-site statistics from mpileup results. Here's an example usage:

```
$ samtools mpileup -B -Q 0 -f circ-Nb-ec3-mtDNA.fasta LAST_OL_132394_vs_mtDNA.bam | \
    /bioinf/scripts/readstomper.pl > LAST_OLmpileupProp_132394_vs_mtDNA.csv

Assembly,Position,Coverage,ref,cR,pR,A,C,G,T,d,i,InsMode
Nb_mtDNA,11456,86,T,79.00,0.92,0.0000,0.0465,0.0000,0.0000,0.03,0.01,A;1.00
Nb_mtDNA,11457,86,A,81.00,0.94,0.0000,0.0000,0.0233,0.0000,0.03,0.02,TA;0.50
Nb_mtDNA,11458,86,G,79.00,0.92,0.0349,0.0233,0.0000,0.0116,0.01,0.00
Nb_mtDNA,11459,86,A,78.00,0.92,0.0000,0.0000,0.0235,0.0471,0.01,0.00
Nb_mtDNA,11460,86,T,5.00,0.06,0.0233,0.8140,0.0233,0.0000,0.08,0.00
Nb_mtDNA,11461,86,G,57.00,0.66,0.1512,0.1163,0.0000,0.0465,0.02,0.01,AATAACAACACGTAACCGA;1.00
Nb_mtDNA,11462,86,G,78.00,0.91,0.0349,0.0349,0.0000,0.0116,0.01,0.00
Nb_mtDNA,11463,86,G,66.00,0.77,0.0465,0.0698,0.0000,0.0116,0.10,0.01,GCAATA;1.00
Nb_mtDNA,11464,86,T,74.00,0.86,0.0000,0.0814,0.0233,0.0000,0.03,0.00

```

It defaults to reporting variants as proportions of the total coverage, but this can be changed to reporting the individual counts:

```
$ samtools mpileup -B -Q 0 -f mmus_chrM.fa mm2_1.2M_vs_Mmus_mtDNA.bam | \
    /bioinf/scripts/readstomper.pl --counts

Assembly,Position,Coverage,ref,skip,cR,pR,A,C,G,T,d,i,InsMode
chrM,1528,491,C,0,409,0.87,17,0,2,15,25,65,A;52
chrM,1529,492,A,0,366,0.75,0,2,2,9,112,12,A;5
chrM,1530,493,A,0,409,0.83,0,2,1,13,68,7
chrM,1531,495,A,0,463,0.94,0,1,10,0,21,8,G;5
chrM,1532,496,A,0,463,0.94,0,2,5,8,16,4,G;3
chrM,1533,495,G,0,459,0.93,3,1,0,2,30,4
chrM,1534,496,A,0,451,0.91,0,0,4,1,39,23,G;21
chrM,1535,503,G,0,391,0.78,4,1,0,4,103,9,A;6
chrM,1536,506,G,0,496,0.98,6,0,0,0,4,10,A;3
chrM,1537,506,G,0,488,0.96,8,1,0,4,5,5,A;3
chrM,1538,506,A,0,483,0.95,0,4,3,4,12,12
chrM,1539,506,C,0,466,0.92,9,0,2,14,14,27,A;6
chrM,1540,506,A,0,469,0.93,0,0,17,6,14,11,T;3
chrM,1541,506,G,0,453,0.90,12,1,0,5,31,21,A;4
chrM,1542,505,C,0,421,0.83,6,0,2,14,62,8,CT;3

```

# Answer
> 0 votes
using BioAlcidaeJdk: http://lindenb.github.io/jvarkit/BioAlcidaeJdk.html , you can loop over the cigar strings of each mapped read and extract the number of inserted bases, clipped bases, etc...

```
java -jar dist/bioalcidaejdk.jar -e 'stream().filter(R->!R.getReadUnmappedFlag() && R.getCigar()!=null).flatMap(R->R.getCigar().getCigarElements().stream()).collect(Collectors.toMap(C->C.getOperator(),C->(long)C.getLength(),(L1,L2)->L1+L2)).forEach((O,N)->println(O.name()+"\t"+N));' toy.bam
I   15
H   11
S   1
D   1
N   14
P   2
M   191

```

to get the number of mismatches, you can preprocess the bam with http://lindenb.github.io/jvarkit/SamFixCigar.html

---
Tags: sequence-alignment, samtools, sam
---


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

# Answer
> 2 votes
Some replies from ONT is:

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

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

---
Tags: nanopore, long-reads, minion, software-installation, minknow
---


# Question
Title: How do kmer counters determine which kmer is 'canonical'?
When counting canonical kmers, ie kmers in which both the forward and reverse complement of a sequence are treated as identical, how do kmer counting programs decide which kmer to use as the canonical sequence? Do they all work the same way?

To investigate I made a string with `GAGTGCGGAATACCACTCTT` which contains all 16 possible 2mers. I then used kmc to figure out how they determine which kmer is used. Only the kmers in the `filtered` column below appeared. So, it looks like KMCs' 'canonical' kmers are the ones that first occur alphabetically.

```
╔════════════════╦═════╦════════════════════╦══════════╗
║ Possible Kmers ║ RCs ║ RC occurs earlier? ║ filtered ║
╠════════════════╬═════╬════════════════════╬══════════╣
║ TT             ║ AA  ║ YES                ║ TA       ║
║ TG             ║ CA  ║ YES                ║ GC       ║
║ TC             ║ GA  ║ YES                ║ GA       ║
║ TA             ║ TA  ║                    ║ CG       ║
║ GT             ║ AC  ║ YES                ║ CC       ║
║ GG             ║ CC  ║ YES                ║ CA       ║
║ GC             ║ GC  ║                    ║ AT       ║
║ GA             ║ TC  ║                    ║ AG       ║
║ CT             ║ AG  ║ YES                ║ AC       ║
║ CG             ║ CG  ║                    ║ AA       ║
║ CC             ║ GG  ║                    ║          ║
║ CA             ║ TG  ║                    ║          ║
║ AT             ║ AT  ║                    ║          ║
║ AG             ║ CT  ║                    ║          ║
║ AC             ║ GT  ║                    ║          ║
║ AA             ║ TT  ║                    ║          ║
╚════════════════╩═════╩════════════════════╩══════════╝

```

Do all kmer counting programs use the same canonical kmers, and if so do you have documentation explaining this? I wasn't able to find anything in the papers for `jellyfish` or `kmc`.

# Answer
> 9 votes
When a k-mer is identical to its reverse complement, both are canonical. Note that a canonical k-mer is a sequence, irrelevant to its position(s) in the input string. More precisely, give a string $s$, its *canonical string* is $$ {\rm canonical}(s|h)=\left\\{\begin{array}{ll} s & \mbox{if $h(s)\<h(\overline{s})$}\\ \overline{s} & \mbox{otherwise}\\ \end{array}\right. $$ where $\overline{s}$ is the Watson-Crick reverse complement of $s$ and $h$ is an arbitrary string hash function. In practice, we most often use 2-bit encoding as $h$. Order under such a hash function is equivalent to the lexicographical order.

---
Tags: k-mer
---


# Question
Title: How do I export my biom object to a file in biom format?
I'm using the python API for the biom-format package. I read in two files in biom format and merged them, like so:

```
In [1]: import biom

In [2]: a = biom.load_table('./table-a.biom')

In [3]: b = biom.load_table('./table-b.biom')

In [4]: a
Out[4]: 51116 x 46 <class 'biom.table.Table'> with 185077 nonzero entries (7% dense)

In [5]: b
Out[5]: 37479 x 26 <class 'biom.table.Table'> with 100124 nonzero entries (10% dense)

In [6]: c = a.merge(b)

In [7]: c
Out[7]: 67022 x 47 <class 'biom.table.Table'> with 254555 nonzero entries (8% dense)

```

How do I now export this biom object (`c`) as a biom file? There are a few methods available, eg, `.to_tsv()`, but I don't see a way to save this object in biom format.

# Answer
You can use the functions `to_hdf5` and `to_json` from the biom-format package:

`to_hdf5`

```
with biom.util.biom_open('output.biom', 'w') as f:
    merged_table.to_hdf5(f, "example")

```

`to_json`

```
with open('output.txt', 'w') as f:
    f.write(merged_table.to_json("example"))

```
> 4 votes

# Answer
There is a function in PiCRUST you can use

```
from picrust.util import write_biom_table

write_biom_table(c,'table-c.biom')

```

Note, however, that depending if you want json or hdf5 format and if your metadata are complex, you may need to use additional arguments. These are the arguments of the function:

```
def write_biom_table(biom_table, biom_table_fp, compress=True,
                 write_hdf5=HAVE_H5PY, format_fs=None):
"""Writes a BIOM table to the specified filepath

Parameters
----------
biom_table : biom.Table
    The table object to write out
biom_table_fp : str
    The path to the output file
compress : bool, optional
    Defaults to ``True``. If True, built-in compression on the output HDF5
    file will be enabled. This option is only relevant if ``write_hdf5`` is
    ``True``.
write_hdf5 : bool, optional
    Defaults to ``True`` if H5PY is installed and to ``False`` if H5PY is
    not installed. If ``True`` the output biom table will be written as an
    HDF5 binary file, otherwise it will be a JSON string.
format_fs : dict, optional
    Formatting functions to be passed to `Table.to_hdf5`

Notes
-----
This code was adapted from QIIME 1.9
"""

```

I have some problems to perform a simple operation and write the output of biom data generated with picrust, it might be helpful as an example of how to use the arguments.
> 1 votes

---
Tags: python
---


title: Web scraping in R

text:
# Question
Title: Web scraping in R
I have ~130 gene names of different human proteins. I'm looking for a convenient and systematic way to locate each gene's location in the genome, extract its CDS sequence and find whether there are SNPs in the CDS (in any context) and what's their position. Moreover, I've chosen for each of those proteins, several amino acids for whom I'm trying to find if the SNP fall in their exact coding nucleotide region (the amino acid triad of nucleotides).

Alternatively, I'm looking for an algorithm that would take the gene name out of my list, find its location in the genome and look for SNPs in the exons - if it would found any SNPs, then it should translate the exon to an amino acid sequence and find whether one of my pre-defined amino acids from each protein is translated by the nucleotide site (triad of nucleotides) with the SNP.

I'm using R, and I thought to use web scraping (should I use the `rvest` package?) but I'm not quite familiar with the technique. Is it possible? What websites should I use? Genome Browser? EMBOSS Transeq? Other websites?

Thank you

# Answer
You can search databases like Cosmic(Cancer) with the gene names to get frequently occurring mutations in cancer. I'm not sure about any specific tool/pipeline available for this exact scenario.

As an alternative approach literature mining for gene associated mutation discovery, you can use text mining with machine learning to automate the whole process.

With no specific context in mind, you will end with several combinations of mutation per gene. Then for each variation translate the DNA sequence and match the result with reference peptide sequence you have.
> 1 votes

# Answer
I don't have a complete answer to your question, but as a primer I can tell you that you should look into the Bioconductor packages. The reason is that biological sequence data are deposited in online databases (like GenBank, and others), and exposed by these databases in such a way that you precisely don't need web scrapping techniques to retrieve the data you're interested in.

`rvest` is a great package for web scrapping, but definitely not the appropriate tool to solve your problem.
> 1 votes

---
Tags: r, proteins, ngs, pdb, uniprot
---


# Question
Title: What is the standard way to measure contig sequence lengths in a BAM?
What is the standard way to measure contig sequence lengths in a BAM?

My understanding is that the community would use `samtools idxstats` to compute this information from the corresponding index file.

Are there more precise/standard approaches?

# Answer
Standard approaches would be either `samtools idxstats` or `samtools view -H`, both of which will produce the exact same results. In fact, all methods will produce the exact same result, since contig length is set by the reference you align against.
> 4 votes

---
Tags: samtools
---


# Question
Title: How do I get the gene order from genbank file using gene names?
I have a list of gene names which I know are present in a genbank file. I want to order these genes based on the order of their presence in the genbank. This should also account for whether they are on the positive or negative strand. So for example I have four genes; genes\_1? are on the positive strand and genes\_2? are negative

```
>nano Gene_names.txt
gene_1A  #assume this is the first gene for positive strand
gene_1B
gene_2A #assume this is first gene for negative strand
gene_2B

>nano Desired_out.txt
gene_1A
gene_2B
gene_1B
gene_2A

```

Edits

The genbank would be the standard flat file as linked below by terdon. The gene names are from a text list file and I have edited them now to reflect this.Im not exactly sure how to approach the strandness issue but I would imagine the complement would be helpful. I am unsure as to the strandness of the genes as is.

A working example of my work is to get the order of SpoI-V and other sporulation genes from C. difficile 630

e.g. sigE, sigG, spoIIE, sleC

# Answer
I can't quite grok your sorting scheme, but you can use genbankr to read your genbank file into a genomic ranges object in R, and filter/sort it however you like:

```
source("https://bioconductor.org/biocLite.R")
biocLite("genbankr")
library(genbankr)
gb = readGenBank("sequence.gb")
subset(cds(gb), gene %in% c('sigE', 'sigG', 'spoIIE', 'sleC', 'REV7'))

```

output:

```
GRanges object with 1 range and 10 metadata columns:
      seqnames       ranges strand |        type        gene        note
         <Rle>    <IRanges>  <Rle> | <character> <character> <character>
  [1]       IX [3300, 4037]      - |         CDS        REV7        <NA>
      codon_start     product  protein_id   translation     loctype
        <numeric> <character> <character> <AAStringSet> <character>
  [1]           1       Rev7p  AAA98667.1     MNR...SLF      normal
          gene_id transcript_id
      <character>   <character>
  [1]        REV7        REV7.1
  -------
  seqinfo: 1 sequence from U49845.1 genome

```
> 1 votes

---
Tags: data-mining, genbank
---


# Question
Title: How does Li and Durbin's BWA paper compare alignment programs on real data?
Li and Durbin's "Fast and accurate short read alignment with burrows-whleeler transform" found here, says:

> We evaluate the performance of BWA on ... real paired-end data by checking the fraction of reads mapped in consistent pairs and by counting misaligned reads mapped against a hybrid genome.

The data for the real genome aligned to is in table 2.

Which pairs are referred to here? Are these the different software tools aligning the same corresponding reads meaning that they aligned the same sequences to similar locations? Are these paired ends? Or is this something different?

# Answer
For the simulated reads, they know *exactly* where the reads came from, because they are able to create that link at the time paired-end reads are generated from the genome. Those simulated reads would be created through a process similar to the following:

1. Pick a single position in the genome at random
2. Pick another region in the same chromosome that is *\<fragment length\>*+*error* bases away from the first location
3. Extend these positions by *\<read length\>* (one read in the forward direction, and the other read in the reverse complement direction) to generate a read pair
4. Change at random some bases (or none) within the reads to simulate sequencing error
5. store the read pair in a FASTA file, together with the original position
6. repeat steps 1-5 a few hundred million times

The effectiveness of the algorithm would be based on how good the algorithm was at picking the correct genomic location.

For the real paired-end data, it's a classical mapping problem: there's no ground truth, so all that can be done is starting at step 5, and working back to all the possibilities for step 1. If the two sequences from a read pair are derived from the same fragment, they should map to approximately the same genomic location and \[usually\] in a reverse-complement orientation to each other. Such an alignment would be called "consistent", or "properly mapped".
> 2 votes

# Answer
What they are pairing is the RNAseq data with the reference genome: as the legend of the table 2 says:

> The 12.2 million read pairs were mapped to the human genome

These softwares align the same read to the same genome. They aren't comparing if it is mapped to the same location, only how many mismatches has each program.
> 2 votes

---
Tags: sequence-alignment, bwa, algorithms
---


# Question
Title: Modifying a SAM header after removing all non-primary reads
I subset a BAM to only include primary reads using the following samtools commands:

```
samtools view -F 256 input.bam > input.primaryOnly.sam

```

Now, in order to convert this SAM into a new subsetted BAM, I would need to add a header. The way I normally do this is with Picard:

```
java -jar picard.jar ReplaceSamHeader \
      I=input.primaryOnly.sam \
      HEADER=input.header.sam \
      O=new.sam

```

whereby the header `input.header.sam` is from the original BAM:

```
samtools view -H input.bam > input.header.sam

```

Then I could convert the SAM into a BAM with

```
samtools view -S -b i new.sam > new.bam

```

My problem is, I'm using the original BAM header. If I look back at the BAM, it now contains wrong information, e.g.

```
samtools flagstat new.bam

```

might save that secondary alignments exist, when they don't.

How do I get a BAM header that makes sense given my new BAM (with only primary alignments)?

# Answer
Your first command only needs a slight modification to add in `-h`. This will create a SAM file with a header.

```
samtools view -h -F 256 input.bam > headered.input.primaryOnly.sam

```

> **-h** Include the header in the output.

If your primary end goal is to create a BAM, then as Pierre has pointed out you can create BAM files directly from `samtools view` with the `-b` option:

```
samtools view -b -F 256 input.bam > input.primaryOnly.bam

```

And if your end-goal is a *sorted* BAM file, then pipe the BAM output through samtools sort:

```
samtools view -b -F 256 input.bam | samtools sort > sorted.input.primaryOnly.bam

```

---

So that's fixing what *I* think your problem is, but what about the specific question that you have asked, i.e. how can a header be added to an existing headerless SAM file?

If you have an existing BAM file, you can use Picard (as you have mentioned), or manually piece the files together with `cat`, which will have a similar function (although Picard probably includes a few checks to make sure that the header is consistent with the input file):

```
(samtools view -H input.bam; cat input.primaryOnly.sam) > headered.input.primaryOnly.sam

```

If the only thing that matters is that the BAM file has correct reference information (which is fairly common), then the `-T` option can be used with `samtools view`. Again, `-h` or `-b` is required to make sure both the header and the usual SAM fields are displayed:

```
samtools view -T reference.fa -h input.primaryOnly.sam > headered.input.primaryOnly.sam
samtools view -T reference.fa -b input.primaryOnly.sam > headered.input.primaryOnly.bam

```
> 5 votes

---
Tags: sam, samtools
---


# Question
Title: How to write a hash function for canonical kmers?
This question is a follow up from How do kmer counters determine which kmer is 'canonical'?.

In that question we learned that kmer counting programs use a 2-bit hash function to internally represent canonical kmers as they are being counted.

Now I am wondering, how can we implement such a function in C/C++ or in python? More arbitrarily, how could I calculate the canonical kmer hash value using a mathematical function?

For example, how would we transform the 3-mer `GAT` or the 21-mer `GAATACCATAGGATA` to 1s and 0s such that:

```
hash(GAT) == hash(ATC)
hash(GAATACCATAGGATA) == hash(TATCCTATGGTATTC)

```

# Answer
You can convert any string hash function to a "canonical" DNA string hash function.

Given a DNA string $s$, let $\overline{s}$ be its Watson-Crick reverse complement. Suppose $h:\Sigma^\*\to\mathbb{Z}$ is an arbitrary string hash function. You can define $$ \tilde{h}(s)\triangleq \min\\{h(s),h(\overline{s})\\} $$ Then for any DNA string $s$ $$ \tilde{h}(s)=\tilde{h}(\overline{s})=\tilde{h}({\rm canonical}(s|h)) $$

PS: There are more than one ways to define "canonical" hash function. For example, let $h\_0(\cdot)$ be the 2-bit encoding hash function. Given any *integer* hash function $g(\cdot)$, we can define: $$ h'(s)\triangleq g\big\[\min\\{h\_0(s),h\_0(\overline{s})\\}\big\] $$ Then we still have $$ h'(s)=h'(\overline{s}) $$ You can also replace $\min\\{\\}$ with XOR or plus to get new definitions.
> 7 votes

# Answer
However you want.

One way to do this is to generate both the forward and reverse complement kmers, then choose the lexicographically-least kmer for the storage key. To delve further into this requires discussion about things like the size of the underlying array, the expected distribution of keys throughout the array, and what type of key clustering is desirable/undesirable. What is actually done doesn't matter, as long as there is a deterministic algorithm that can be applied to every single possible sequence.

The nice thing about most modern languages is that they've already got their own hash implementations that map keys to values, which work fine for quick experiments: C++ has unordered map, perl has hashes, python has dictionaries, and R has names.

The actual bare-metal implementation of the thing that converts a string of text into a number is complicated, particularly if you want code to be cryptographically secure, and the hash array to have an access speed that is very close to the amount of time it takes to retrieve a single value from memory. A useful starting point to explore this in more depth might be the Wikipedia page on the Java hashcode function.
> 5 votes

---
Tags: python, k-mer, c++, hashing
---


# Question
Title: Error in checkFullRank(modelMatrix) deseq2
Trying to use deseq2 for differential expression analysis (rna-seq) between three groups and also account for batch effect as the control were sequenced at a different time point.

1. control: con
2. sample with mutation A: mutA
3. sample with mutation B: mutB

here is my design

```
   sample condition batch
1     C_1         con     1
2     C_2         con     1
3     C_5         con     1
4     C_3         con     1
5     C_4         con     1
6     M_6        mutA     2
7     M_2        mutA     2
8     M_5        mutA     2
9     M_1        mutA     2
10    M_4        mutA     2
11    M_3        mutA     2
12   MA_2        mutB     2
13   MA_6        mutB     2
14   MA_5        mutB     2
15   MA_4        mutB     2
16   MA_3        mutB     2
17   MA_1        mutB     2

```

when i create matrix:

```
dds <- DESeqDataSetFromMatrix(countData = data, colData = coldata, design = ~batch + condition)

```

I see this error:

```
Error in checkFullRank(modelMatrix) :
  the model matrix is not full rank, so the model cannot be fit as specified.
  One or more variables or interaction terms in the design formula are linear
  combinations of the others and must be removed.

```

I also tried to compare only two conditions subsetting the original data, i still see the same error what am i missing?

# Answer
> 4 votes
See the second example from the vignette, which is exactly your situation. The software cannot differentiate between the batch effects vs. condition effects because the groups are exactly the same.

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


# Question
Title: What are the differences between "Site" and "CONECT" in a PDB file
The `CONECT` section in a PDB file describes atomic connectivity.

The `SITE` section describes (among other characteristics), environments surrounding ligands present in the structure.

Now, how come not every atom in the `CONECT` section is described in the `SITE` section?

For example: 1a28.pdb

# Answer
> 2 votes
> how come not every atom in the CONECT section is described in the SITE section?

Because it is not required and not needed.

`CONECT` is for listing bonds, so that molecular graphics programs know how to draw the structure.

`SITE` is for annotation of interesting sites. Each listed SITE has a corresponding REMARK 800 that "details its significance". A graphical program may provide a menu with a list of all interesting sites for easy navigation.

Note: PDB is a legacy file format. The primary file format used nowadays by the PDB (organization) is called PDBx/mmCIF. So when looking into such details it is often better to look directly into mmCIF files.

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


# Question
Title: Parse BAM by insertion size and get genomic coordinates
I would like to find out where (in genomic coordinates) large insertions are found within a given BAM file, `file.bam`.

(In terms of genomic coordinates, I just mean I would like to have a rough idea where to look at the BAM in IGV, chromosome and position.)

Based on the BAM CIGAR, I should be able to `grep` the BAM based on some cut-off of "large inserted" bases.

It's unclear to me how to best do this via samtools; one could parse the cigar string via `samtools` for insertions, but I don't understand how to translate this into the genomic position in the alignment.

This goes with various samtools wrappers, e.g. `pysam` in python:

```
import pysam

bamfile = pysam.Samfile(bamFile, "rb");

for read in bamfile:
    if not read.is_unmapped:   ## only get mapped reads
        cigarLine = read.cigar
        for cigarType, cigarLength in cigarLine:
            if cigarType == 1: ## this means it's an insertion...
                  ## somehow filter only large insertions, and then get coordinate information

```

I'm not sure how to do this in a non-clumsy manner. Any help appreciated.

# Answer
> 2 votes
You basically have the right framework already in your pysam code, which I took the liberty of tidying up a bit. You just need to add a bit:

```
for read in bamfile:
    if not read.is_unmapped:
        cigarLine = read.cigar
        interesting = False
        for cigarType, cigarLength in cigarLine:
            if cigarType == 1:
                if cigarLength > 100:  ## Arbitrary value
                    interesting = True
                    break
        if interesting:
            print("{}:{}-{}".format(read.reference_name, read.reference_start, read.reference_end))

```

Unless you have long reads, this will be close enough for IGV.

# Answer
> 1 votes
a one-liner using **bioalcidaejdk**: http://lindenb.github.io/jvarkit/BioAlcidaeJdk.html

```
$ java -jar dist/bioalcidaejdk.jar -e 'stream().filter(R->!R.getReadUnmappedFlag() && R.getCigar().getCigarElements().stream().anyMatch(C->C.getLength()>100 && (C.getOperator().equals(CigarOperator.N) || C.getOperator().equals(CigarOperator.D)))).forEach(R->println(R.getContig()+"\t"+R.getStart()+"\t"+R.getEnd()));' input.bam

```

---
Tags: sam, samtools
---


# Question
Title: Is it possible to find out the PacBio chemistry from the PacBio Sequel BAM's header?
Is it possible to find out the PacBio chemistry from the PacBio Sequel BAM's header? E.g. whether it was produced with P5-C3 or P6-C4? I checked mine but I haven't found it but maybe it's hidden in one of the tags whose meaning I don't know yet, like SEQUENCINGKIT=100-861-800?

# Answer
> 1 votes
According to PacBio Sequel BAM format manual, you can find information on the chemistry version in the header entries `@RG` (read group), tag `DS` (description), entry `BASECALLERVERSION` .

In case of future very PacBio-specific questions, you can contact Liz Tseng on Twitter (@magdoll), one of the most responsive PacBio developers. They also have a dedicated SeqAnswers forum, which is directly monitored by them and you are very likely to get a fast and accurate answer.

---
Tags: sam, pacbio
---


# Question
Title: What to do with the configuration test cell after configuring the MinION with it?
What to do with the configuration test cell after configuring the MinION with it?

I received a MinION with a configuration test cell in it. After running a configuration in the MinKNOW software while the CTC was in the MinION, I wonder what to do with the CTC? Am I going to need it later or can I throw it away? If I should keep it, in what conditions should I keep it (temperature, light, etc.)?

Also, do you know why we need to do a configuration of the MinION the first time?

I posted this question on Biology Stackexchange as well but maybe the community here will be more receptive?

# Answer
I would keep it, in case you need to do another configuration later on. I have used mine for that purpose.

I think storage doesn't really matter, I just put it in a drawer somewhere. Keep it away from liquids and extreme temperatures :)
> 2 votes

---
Tags: nanopore, ngs, long-reads, minion
---


# Question
Title: Generating burrows-wheeler-transform and a suffix array of a DNA sequence with less memory
Also posted on biostars

After indexing

```
bwa index -a bwtsw reference.fa

```

I got files. like .bwt file and .sa file. The naive way of generating these file is:

1. Create all rotations of T
2. Sort the rotations lexicographically
3. Concatenate the last character of each rotation to form BWT(T)

I used this naive way in my algorithm. The problem is it took too much memory.

My question is: is there an effective memory way to generate .bwt file and .sa file ?

Any help will be appreciated.

# Answer
> 4 votes
No, there is no way to decompress .bwt or .sa files, because they are not compressed. After indexing you'll get these files, which some programs (bwa) use. They should be in the same folder as your reference fasta file. What are you trying to do with these files? To get proper help here, please explain why you need these files, and show how you made the files.

---
Tags: sequence-alignment, bwa
---


title: How BWA generate index files?

text:
# Question
Title: How BWA generate index files?
I need to know how BWA generate bw and sa in less memory usage. I didn't find details in the original paper. Do they save all rotations temporary ? I need to know do BWA stores all rotations so that it can sort the rotations alphabetically ?

Thanks.

# Answer
## Papers

The space-efficient BWT construction algorithm is in bwt\_gen.c. One of the first things you see there is the copyright line:

> Copyright (C) 2004, Wong Chi Kwong.

With that name, you can find the relevant papers:

The third paper only mentions BWT briefly in the last section. It deals with compressed suffix arrays, which are space-efficient representations of the inverse function of LF-mapping. But if you have a representation of (the inverse of) the LF-mapping, you have a representation of the BWT. Hence the algorithm can be adapted for BWT construction.

## Algorithm

When building BWT, you usually append a unique character `$` to the end of the text, either implicitly or explicitly. This way you can reduce sorting rotations to sorting suffixes. You sort the suffixes in lexicographic order, and then you output the preceding character for each suffix in that order to form the BWT.

Assume that you already have the BWT for string `X`, and you want to transform it into the BWT of string `cX`, where `c` is a single character. You can do that in three steps (this assumes 0-based array indexing):

1. Find the position `i` of the full suffix `X`, and replace the `$` in the BWT with the inserted character `c` (set `BWT[i] = c`).
2. Determine the number of suffixes `j` that are lexicographically smaller than `cX` with `j = LF(i, c) = C[c] + BWT.rank(i, c)`. Here `C[c]` is the number of suffixes that start with a character smaller than `c`, and `BWT.rank(i, c)` is the number of occurrences of `c` in the prefix `BWT[0..i-1]`.
3. Insert `$` between `BWT[j-1]` and `BWT[j]`.

By iterating this, you can build the BWT for any string slowly but space-efficiently. The algorithm of Hon et al. is based on the same idea, but they make it faster by inserting multiple characters at the same time. Given the BWT of string `X` and another string `Y`, they produce the BWT of the concatenation `YX`. The details are quite complicated, but you can find them in the third paper and in the source code.
> 4 votes

# Answer
The following is only a comment to @JouniS' answer, but it is too long to fit in 512 characters.

The algorithm @JouniS was describing has been rediscovered multiple times. The C program below constructs BWT with this algorithm. It is proof-of-concept only, not intended for long strings.

```
#include <string.h>
#include <stdio.h>
int str2bwt(char *s)
{
    int i, l, c = 0, k = 0;
    l = strlen(s); // length of input string
    for (i = l - 1; i >= 0; --i) {
        int j, r = 0, a = s[i];
        memmove(&s[i], &s[i+1], k); // make room for insertion
        s[i+k] = a; // insert the char
        for (j = i; j < i + k; ++j) r += (s[j] <= a);
        for (; j < l; ++j) r += (s[j] < a);
        k = r + 1, c = a;
    }
    return k; // the position of sentinel '$'
}
int main(int argc, char *argv[])
{
    int k;
    if (argc == 1) return 1;
    k = str2bwt(argv[1]);
    fwrite(argv[1], 1, k, stdout); putchar('$'); puts(&argv[1][k]);
    return 0;
}

```

PS: this is one of the simplest ways to construct BWT and does not require any additional dependencies. It is slow for long strings, though.
> 3 votes

---
Tags: sequence-alignment, bwa
---


# Question
Title: Assuming that we have the following SNP and phenotype data, is the SNP significantly associated with the phenotype?
So if I attempt this question using their method I put the data of genotype and phenotype into a data frame in R then use `lm()` to do linear regression and `summary()` to get an F statistic with p value 0.54. I'm not connecting what the F statistic means in this case to the answer.

Question 2: Assuming that we have the following SNP and phenotype data, is the SNP significantly associated with the phenotype? Here, we represent each SNP site as the number of minor alleles on that locus, so 0 and 2 are for major and minor homozygous sites, respectively, and 1 is for the heterozygous sites. We also assume that minor alleles contribute to the phenotype and the effect is additive. In other words, the effect from a minor homozygous site should be twice as large as that from a heterozygous site. You may use any test methods introduced in the chapter. How about permutation tests?

> GENOTYPE: 1, 2, 2, 1, 1, 0, 2, 0, 1, 0
>
> PHENOTYPE: 0.53, 0.78, 0.81, -0.23, -0.73, 0.81, 0.27, 2.59, 1.84, 0.03

Answer: Take the linear regression with F-test as an example, the p-value without the permutation test is 0.54 and is 0.55 after 1000 permutation (this number could be different with your answer because of the randomness of the permutation test).

Link to where I'm getting this: http://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1002828#pcbi.1002828-Balding1

# Answer
Your approach is correct - it's not the only one permitted under their conditions, but it would be the standard way given only that information.

F-test compares a "full model" (predicting the phenotype using genotype and an intercept) versus a "reduced model" (predicting the phenotype using only the intercept, i.e. the mean value of phenotype). The p-value tells whether the full model predicts better, i.e. whether adding the genotype variable significantly improves the prediction. In this single-variable case, it should be approximately equal to the p-value for the beta coefficient of the genotype, if that's more familiar.

Anyway, this isn't specific to genetics, so you can check out other online sources - e.g. https://onlinecourses.science.psu.edu/stat501/node/295 for a longer, but accessible explanation of the F-test, and https://stats.stackexchange.com/ for many answers to similar questions.
> 2 votes

---
Tags: snp, gwas
---