text:
Installing Sibelia for Ragout on Mac OSX
I am trying to use Ragout: https://github.com/fenderglass/Ragout
to fill the gaps in my de novo genome assembly. You can access the article freely here: https://www.ncbi.nlm.nih.gov/pubmed/24931998
For this, I first need to install Sibelia.
I tried doing so by cloning the GitHub repository into output/software/ragout/ and then running the following command:
git clone https://github.com/fenderglass/Ragout output/software/ragout

python2 output/software/ragout/scripts/install-sibelia.py

Installing Sibelia
Downloading source...
-- The CXX compiler identification is AppleClang 9.0.0.9000037
-- Check for working CXX compiler: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/c++
-- Check for working CXX compiler: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/c++ -- works
-- Detecting CXX compiler ABI info
-- Detecting CXX compiler ABI info - done
-- Detecting CXX compile features
-- Detecting CXX compile features - done
-- The C compiler identification is AppleClang 9.0.0.9000037
-- Check for working C compiler: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/cc
-- Check for working C compiler: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/cc -- works
-- Detecting C compiler ABI info
-- Detecting C compiler ABI info - done
-- Detecting C compile features
-- Detecting C compile features - done
-- Performing Test HAVE_UNKNOWN_WALL
-- Performing Test HAVE_UNKNOWN_WALL - Success
-- Performing Test HAVE_UNKNOWN_FOMIT_FRAME_POINTER
-- Performing Test HAVE_UNKNOWN_FOMIT_FRAME_POINTER - Success
-- Looking for inttypes.h
-- Looking for inttypes.h - found
-- Looking for memory.h
-- Looking for memory.h - found
-- Looking for stddef.h
-- Looking for stddef.h - found
-- Looking for stdint.h
-- Looking for stdint.h - found
-- Looking for stdlib.h
-- Looking for stdlib.h - found
-- Looking for string.h
-- Looking for string.h - found
-- Looking for strings.h
-- Looking for strings.h - found
-- Looking for sys/types.h
-- Looking for sys/types.h - found
-- Performing Test HAVE_INLINE
-- Performing Test HAVE_INLINE - Success
-- Performing Test HAVE___INLINE
-- Performing Test HAVE___INLINE - Success
-- Performing Test HAVE___INLINE__
-- Performing Test HAVE___INLINE__ - Success
-- Performing Test HAVE___DECLSPEC_DLLEXPORT_
-- Performing Test HAVE___DECLSPEC_DLLEXPORT_ - Failed
-- Performing Test HAVE___DECLSPEC_DLLIMPORT_
-- Performing Test HAVE___DECLSPEC_DLLIMPORT_ - Failed
-- Check size of uint8_t
-- Check size of uint8_t - done
-- Check size of int32_t
-- Check size of int32_t - done
-- Looking for PRId32
-- Looking for PRId32 - found
-- Configuring done
-- Generating done
-- Build files have been written to: /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/build
Scanning dependencies of target divsufsort
[  4%] Building C object libdivsufsort-2.0.1/lib/CMakeFiles/divsufsort.dir/divsufsort.o
[  8%] Building C object libdivsufsort-2.0.1/lib/CMakeFiles/divsufsort.dir/sssort.o
[ 12%] Building C object libdivsufsort-2.0.1/lib/CMakeFiles/divsufsort.dir/trsort.o
[ 16%] Building C object libdivsufsort-2.0.1/lib/CMakeFiles/divsufsort.dir/utils.o
[ 20%] Linking C static library libdivsufsort.a
[ 20%] Built target divsufsort
Scanning dependencies of target Sibelia
[ 24%] Building CXX object CMakeFiles/Sibelia.dir/sibelia.cpp.o
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/sibelia.cpp:8:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/postprocessor.h:10:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/align.h:43:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/sequence.h:103:
/Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/sequence/string_set_dependent_generous.h:38:9: warning:
      'SEQAN_SEQUENCE_STRING_SET_DEPENDENT_GENEROUS_H_' is used as a header
      guard here, followed by #define of a different macro [-Wheader-guard]
#ifndef SEQAN_SEQUENCE_STRING_SET_DEPENDENT_GENEROUS_H_
        ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/sequence/string_set_dependent_generous.h:39:9: note:
      'SEQAN_SEQUENCE_STRING_SET_DEPENDENT_GENEROUSH_' is defined here; did you
      mean 'SEQAN_SEQUENCE_STRING_SET_DEPENDENT_GENEROUS_H_'?
#define SEQAN_SEQUENCE_STRING_SET_DEPENDENT_GENEROUSH_
        ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        SEQAN_SEQUENCE_STRING_SET_DEPENDENT_GENEROUS_H_
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/sibelia.cpp:8:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/postprocessor.h:10:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/align.h:44:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/score.h:48:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/score/score_matrix.h:40:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/file.h:70:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/system.h:82:
/Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/system/system_sema.h:120:27: warning:
      'sem_init' is deprecated [-Wdeprecated-declarations]
            SEQAN_DO_SYS(!sem_init(hSemaphore, 0, init));
                          ^
/usr/include/sys/semaphore.h:55:42: note: 'sem_init' has been explicitly marked
      deprecated here
int sem_init(sem_t *, int, unsigned int) __deprecated;
                                         ^
/usr/include/sys/cdefs.h:176:37: note: expanded from macro '__deprecated'
#define __deprecated    __attribute__((deprecated))
                                       ^
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/sibelia.cpp:8:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/postprocessor.h:10:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/align.h:44:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/score.h:48:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/score/score_matrix.h:40:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/file.h:70:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/system.h:82:
/Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/system/system_sema.h:124:27: warning:
      'sem_destroy' is deprecated [-Wdeprecated-declarations]
            SEQAN_DO_SYS(!sem_destroy(hSemaphore));
                          ^
/usr/include/sys/semaphore.h:53:26: note: 'sem_destroy' has been explicitly
      marked deprecated here
int sem_destroy(sem_t *) __deprecated;
                         ^
/usr/include/sys/cdefs.h:176:37: note: expanded from macro '__deprecated'
#define __deprecated    __attribute__((deprecated))
                                       ^
3 warnings generated.
[ 28%] Building CXX object CMakeFiles/Sibelia.dir/postprocessor.cpp.o
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/postprocessor.cpp:9:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/postprocessor.h:10:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/align.h:43:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/sequence.h:103:
/Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/sequence/string_set_dependent_generous.h:38:9: warning:
      'SEQAN_SEQUENCE_STRING_SET_DEPENDENT_GENEROUS_H_' is used as a header
      guard here, followed by #define of a different macro [-Wheader-guard]
#ifndef SEQAN_SEQUENCE_STRING_SET_DEPENDENT_GENEROUS_H_
        ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/sequence/string_set_dependent_generous.h:39:9: note:
      'SEQAN_SEQUENCE_STRING_SET_DEPENDENT_GENEROUSH_' is defined here; did you
      mean 'SEQAN_SEQUENCE_STRING_SET_DEPENDENT_GENEROUS_H_'?
#define SEQAN_SEQUENCE_STRING_SET_DEPENDENT_GENEROUSH_
        ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        SEQAN_SEQUENCE_STRING_SET_DEPENDENT_GENEROUS_H_
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/postprocessor.cpp:9:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/postprocessor.h:10:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/align.h:44:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/score.h:48:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/score/score_matrix.h:40:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/file.h:70:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/system.h:82:
/Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/system/system_sema.h:120:27: warning:
      'sem_init' is deprecated [-Wdeprecated-declarations]
            SEQAN_DO_SYS(!sem_init(hSemaphore, 0, init));
                          ^
/usr/include/sys/semaphore.h:55:42: note: 'sem_init' has been explicitly marked
      deprecated here
int sem_init(sem_t *, int, unsigned int) __deprecated;
                                         ^
/usr/include/sys/cdefs.h:176:37: note: expanded from macro '__deprecated'
#define __deprecated    __attribute__((deprecated))
                                       ^
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/postprocessor.cpp:9:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/postprocessor.h:10:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/align.h:44:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/score.h:48:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/score/score_matrix.h:40:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/file.h:70:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/system.h:82:
/Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/system/system_sema.h:124:27: warning:
      'sem_destroy' is deprecated [-Wdeprecated-declarations]
            SEQAN_DO_SYS(!sem_destroy(hSemaphore));
                          ^
/usr/include/sys/semaphore.h:53:26: note: 'sem_destroy' has been explicitly
      marked deprecated here
int sem_destroy(sem_t *) __deprecated;
                         ^
/usr/include/sys/cdefs.h:176:37: note: expanded from macro '__deprecated'
#define __deprecated    __attribute__((deprecated))
                                       ^
3 warnings generated.
[ 32%] Building CXX object CMakeFiles/Sibelia.dir/indexedsequence.cpp.o
[ 36%] Building CXX object CMakeFiles/Sibelia.dir/util.cpp.o
[ 40%] Building CXX object CMakeFiles/Sibelia.dir/outputgenerator.cpp.o
[ 44%] Building CXX object CMakeFiles/Sibelia.dir/blockfinder.cpp.o
[ 48%] Building CXX object CMakeFiles/Sibelia.dir/blockinstance.cpp.o
[ 52%] Building CXX object CMakeFiles/Sibelia.dir/bifurcationstorage.cpp.o
[ 56%] Building CXX object CMakeFiles/Sibelia.dir/bulgeremoval.cpp.o
[ 60%] Building CXX object CMakeFiles/Sibelia.dir/dnasequence.cpp.o
[ 64%] Building CXX object CMakeFiles/Sibelia.dir/edge.cpp.o
[ 68%] Building CXX object CMakeFiles/Sibelia.dir/fasta.cpp.o
[ 72%] Building CXX object CMakeFiles/Sibelia.dir/serialization.cpp.o
[ 76%] Building CXX object CMakeFiles/Sibelia.dir/synteny.cpp.o
[ 80%] Building CXX object CMakeFiles/Sibelia.dir/test/unrolledlisttest.cpp.o
[ 84%] Building CXX object CMakeFiles/Sibelia.dir/platform.cpp.o
[ 88%] Building CXX object CMakeFiles/Sibelia.dir/stranditerator.cpp.o
[ 92%] Building CXX object CMakeFiles/Sibelia.dir/vertexenumeration.cpp.o
[ 96%] Building CXX object CMakeFiles/Sibelia.dir/resource.cpp.o
[100%] Linking CXX executable Sibelia
[100%] Built target Sibelia
[ 20%] Built target divsufsort
[100%] Built target Sibelia
Install the project...
-- Install configuration: "Release"
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/bin/Sibelia
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/NEWS.md
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/ANNOTATION.md
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/README.md
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/USAGE.md
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/INSTALL.md
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/SIBELIA.md
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/C-SIBELIA.md
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/LICENSE.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/C-Sibelia
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/C-Sibelia/Staphylococcus_aureus
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/C-Sibelia/Staphylococcus_aureus/NCTC8325.fasta
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/C-Sibelia/Staphylococcus_aureus/README.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/C-Sibelia/Staphylococcus_aureus/RN4220.fasta
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/C-Sibelia/Staphylococcus_aureus/variant.vcf
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/blocks_coords.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/circos
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/circos/circos.conf
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/circos/circos.highlight.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/circos/circos.highlight1.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/circos/circos.highlight2.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/circos/circos.highlight3.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/circos/circos.highlight4.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/circos/circos.image.conf
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/circos/circos.png
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/circos/circos.segdup.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/circos/circos.sequences.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/circos/circos.svg
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/coverage_report.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/d3_blocks_diagram.html
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/genomes_permutations.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/Helicobacter_pylori.fasta
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/README.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Staphylococcus_aureus
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Staphylococcus_aureus/blocks_coords.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Staphylococcus_aureus/circos
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Staphylococcus_aureus/circos/circos.conf
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Staphylococcus_aureus/circos/circos.highlight.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Staphylococcus_aureus/circos/circos.image.conf
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Staphylococcus_aureus/circos/circos.png
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Staphylococcus_aureus/circos/circos.segdup.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Staphylococcus_aureus/circos/circos.sequences.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Staphylococcus_aureus/circos/circos.svg
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Staphylococcus_aureus/coverage_report.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Staphylococcus_aureus/d3_blocks_diagram.html
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Staphylococcus_aureus/genomes_permutations.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Staphylococcus_aureus/README.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Staphylococcus_aureus/Staphylococcus.fasta
Traceback (most recent call last):
  File "output/software/ragout/scripts/install-sibelia.py", line 111, in <module>
    sys.exit(main())
  File "output/software/ragout/scripts/install-sibelia.py", line 107, in main
    return int(not install_deps(args.prefix))
  File "output/software/ragout/scripts/install-sibelia.py", line 31, in install_deps
    return install_sibelia(prefix)
  File "output/software/ragout/scripts/install-sibelia.py", line 64, in install_sibelia
    shutil.copy(sibelia_bin_src, sibelia_bin_dst)
  File "/usr/local/Cellar/python/2.7.13/Frameworks/Python.framework/Versions/2.7/lib/python2.7/shutil.py", line 119, in copy
    copyfile(src, dst)
  File "/usr/local/Cellar/python/2.7.13/Frameworks/Python.framework/Versions/2.7/lib/python2.7/shutil.py", line 83, in copyfile
    with open(dst, 'wb') as fdst:
IOError: [Errno 2] No such file or directory: '/Users/cr517/Documents/phd/project/lib/Sibelia'

As you can see, I get a lot of errors. Could someone help me, please?
I also tried to install Sibelia following their instructions here: https://github.com/bioinf/Sibelia
cd Sibelia/build
➜  build git:(master) cmake ../src
-- The CXX compiler identification is AppleClang 9.0.0.9000037
-- Check for working CXX compiler: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/c++
-- Check for working CXX compiler: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/c++ -- works
-- Detecting CXX compiler ABI info
-- Detecting CXX compiler ABI info - done
-- Detecting CXX compile features
-- Detecting CXX compile features - done
-- The C compiler identification is AppleClang 9.0.0.9000037
-- Check for working C compiler: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/cc
-- Check for working C compiler: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/cc -- works
-- Detecting C compiler ABI info
-- Detecting C compiler ABI info - done
-- Detecting C compile features
-- Detecting C compile features - done
-- Performing Test HAVE_UNKNOWN_WALL
-- Performing Test HAVE_UNKNOWN_WALL - Success
-- Performing Test HAVE_UNKNOWN_FOMIT_FRAME_POINTER
-- Performing Test HAVE_UNKNOWN_FOMIT_FRAME_POINTER - Success
-- Looking for inttypes.h
-- Looking for inttypes.h - found
-- Looking for memory.h
-- Looking for memory.h - found
-- Looking for stddef.h
-- Looking for stddef.h - found
-- Looking for stdint.h
-- Looking for stdint.h - found
-- Looking for stdlib.h
-- Looking for stdlib.h - found
-- Looking for string.h
-- Looking for string.h - found
-- Looking for strings.h
-- Looking for strings.h - found
-- Looking for sys/types.h
-- Looking for sys/types.h - found
-- Performing Test HAVE_INLINE
-- Performing Test HAVE_INLINE - Success
-- Performing Test HAVE___INLINE
-- Performing Test HAVE___INLINE - Success
-- Performing Test HAVE___INLINE__
-- Performing Test HAVE___INLINE__ - Success
-- Performing Test HAVE___DECLSPEC_DLLEXPORT_
-- Performing Test HAVE___DECLSPEC_DLLEXPORT_ - Failed
-- Performing Test HAVE___DECLSPEC_DLLIMPORT_
-- Performing Test HAVE___DECLSPEC_DLLIMPORT_ - Failed
-- Check size of uint8_t
-- Check size of uint8_t - done
-- Check size of int32_t
-- Check size of int32_t - done
-- Looking for PRId32
-- Looking for PRId32 - found
-- Configuring done
-- Generating done
-- Build files have been written to: /Users/cr517/Documents/phd/project/Sibelia/build
➜  build git:(master) make
Scanning dependencies of target lagan
[  3%] Creating directories for 'lagan'
[  6%] Performing download step (DIR copy) for 'lagan'
[  9%] No patch step for 'lagan'
[ 12%] No update step for 'lagan'
[ 15%] No configure step for 'lagan'
[ 18%] Performing build step for 'lagan'
rightinfluence.cpp:21:2: error: reference to 'end' is ambiguous
        end.score = -2;
        ^
rightinfluence.cpp:3:18: note: candidate found by name lookup is 'end'
Fragment origin, end;
                 ^
/Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/include/c++/v1/iterator:1768:1: note:
      candidate found by name lookup is 'std::__1::end'
end(const _Cp& __c)
^
/Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/include/c++/v1/iterator:1611:1: note:
      candidate found by name lookup is 'std::__1::end'
end(_Tp (&__array)[_Np])
^
/Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/include/c++/v1/iterator:1760:1: note:
      candidate found by name lookup is 'std::__1::end'
end(_Cp& __c)
^
rightinfluence.cpp:22:22: error: reference to 'end' is ambiguous
        origin.totalScore = end.totalScore = 0;
                            ^
rightinfluence.cpp:3:18: note: candidate found by name lookup is 'end'
Fragment origin, end;
                 ^
/Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/include/c++/v1/iterator:1768:1: note:

...

There might be two issues, one is a hard coded or an incorrect path, the other one might be due to not using sudo or root rights

Hi, I tried doing the commands with sudo but I still get similar errors. The first failure is: -- Performing Test HAVE___DECLSPEC_DLLEXPORT_ -- Performing Test HAVE___DECLSPEC_DLLEXPORT_ - Failed -- Performing Test HAVE___DECLSPEC_DLLIMPORT_ -- Performing Test HAVE___DECLSPEC_DLLIMPORT_ - Failed

@Llopis, how can I improve my attempt at installing this software with non-hard-coded path and correct path?

This is a fairly obscure package, I have a feeling that you'll end up needing to post this as a github issue and hope that the developers reply there.

I posted the issue on Github here: https://github.com/fenderglass/Ragout/issues/23

@charlesdarwin modify the source code before installing or via a patch in github. However, make the issues self-contained, include the version number, the machine description and everything that might be relevant

@Llopis, sorry, what do you mean by modifying the source code? Modify the code that's present in here: https://github.com/fenderglass/Ragout ? And how to I get the version number? I will add my machine specs: Mac OSX Sierra version 10.12.6

Yes, that is the source code of the application you want to install

OK, I'll have a look at it and study it.

There is a link to compiled binaries for OSX on the project homepage

Indeed, there are also Ragout binaries available on github. The easiest solution to all of this is to not compile and instead just use the binaries.

Hi, regarding Sibelia, I have done the following: In R: download.file('https://sourceforge.net/projects/sibelia-bio/files/3.0.6/Sibelia-3.0.6-Source.tar.gz', 'output/software/sibelia/Sibelia-3.0.6-Darwin.tar.gz')
In bash: tar -xvzf output/software/sibelia/Sibelia-3.0.6-Darwin.tar.gz -C output/software/sibelia/ chmod +rwxrwxrwx output/software/sibelia/Sibelia-3.0.6-Source/src/csibelia/C-Sibelia.py
This gives me access to the C-Sibelia.py software but I cannot find any Sibelia.py in my computer. Do you know if I am downloading the right thing?

I tried installing Ragout from the binaries:
In R: download.file('https://github.com/fenderglass/Ragout/releases/download/2.0/ragout-2.0-macos-x86_64.tar.gz', 'ragout-2.0-macos-x86_64.tar.gz')
In bash:
tar -xvzf ragout-2.0-macos-x86_64.tar.gz  chmod +rwxrwxrwx ragout-2.0-macos-x86_64/ragout.py ragout-2.0-macos-x86_64/ragout.py -h The help page is then opened so I consider it installed.


text:
help with CRCmapper
I'm using CRCmapper on my data. There are several modules in the tool. I've used exactly instructed in the usage.
Commands used:
python CRCmapper.py -e MB12_noEncBlack_AllEnhancers.table.txt -b MB12.bam -g HG19 -f hg19-stdChr/ -s MB12_peaks.bed -n MB12 -o MB12/

Here is the error I got:
Warning: text mode turns off computation of q-values
IDENTIFY TF-TF INTERACTIONS
Traceback (most recent call last):
  File "CRCmapper.py", line 605, in <module>
    main()
  File "CRCmapper.py", line 593, in main
    graph = buildNetwork(projectFolder, projectName, candidateGenes, refseqToNameDict, motifConvertFile)
  File "CRCmapper.py", line 381, in buildNetwork
    target = refseqToNameDict[region[0]]
KeyError: 'AP1'

My expertise in python is very basic. I would be grateful if someone take a look at this error and suggest how to fix.
P.S. I contacted authors a long back, I didn't get any response.
Thank you.

It seems that the problem is within the files, are you sure that AP1 is in all your files? It could also be a bug on the software

Yeah, it exists. For other samples, it is a different motif for different sample. It seems to be a bug in the software but I've no idea. Hope someone here will be able to understand.

Could you paste the minimal files that reproduce the bug? Maybe you can even propose a change in the repository to fix it

Yeah, I will prepare an example to reproduce. They haven't enabled the issue creator for their repo.

Well you can do directly a pull request

It seems that the author of CRCmapper.py wants to extract sequence_name, start, stop columns from fimo.txt file. But the code extract motif_alt_id, sequence_name, start columns. So the index of the code should be changed.
fimo.txt looks like this:
motif_id    motif_alt_id    sequence_name   start   stop    strand  score   p-value q-value matched_sequence
Transfac.V$AP1_02   AP1 NM_003046931|3|95762|152118 2332    3118    -   11.9802 6.49e-05        TGACTCA

In line 377-382 of the code:
line[1] should be changed to line[2].
for line in fimoTable[1:]:
    source = motifDatabaseDict[line[0]]
    region = line[1].split('|')
    target = refseqToNameDict[region[0]]
    location = (region[1], int(region[2]), int(region[3]))

In line 391-394 of the code:
line[2] should be changed to line[3], line[3] should be changed to line[4].
    # Count unique motifs
    if (region[1], int(region[2]) + int(line[2]), int(region[2]) + int(line[3])) not in motifDictSE[source]:
        edgeCountDictSE[(source, target)] += 1
        motifDictSE[source].append((region[1], int(region[2]) + int(line[2]), int(region[2]) + int(line[3])))

This error came from the fact that  the output of the new fimo version has an additional 'motif_alt_id' column in second position.
CRCmapper has been modified to take this into account. You can find the updated version here: https://bitbucket.org/young_computation/crcmapper/src/master/
P.S.: sorry if we missed you request user1545 and thanks for using CRCmapper!

source: Stack Exchange

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

sort: votes


text:
How are centromeric/telomeric regions categorized in RepeatMasker? "Other" category?
I haven't been able to find the following in the documentation at http://www.repeatmasker.org/
(1) Why are certain regions classified as "Other"? Are these regions impossible to concretely classify due to alignment issues?
(2) I believe that centromeric and telomeric regions would be classified as tandem repeats. Are telomeric/centromeric regions classified as "DNA" within the hg19 version of RepeatMasker?
Are these documented somewhere?

I'll leave question 1 to someone else since I don't know.
Regarding question 2, centromeric and telomeric regions are hard masked in hg19, so there's no sequence to repeat mask.
For hg38, the centromeres are labeled ALR/Alpha with a class of centr. Telomeres are still hard-masked, but they tend to have simple repeats next to them (as annotated by repeatmasker, which makes sense).
Apparently the Other class in hg19 is comprised of SVA elements. In hg38, these are reclassified as Retrotransposon class and SVA family, which didn't exist in hg19.

"Apparently the Other class in hg19 is comprised of SVA elements. In hg38, these are reclassified as Retrotransposon class and SVA family, which didn't exist in hg19."

Interesting. Where did you find this information? Aren't SINE elements SVA elements?

Actually, do you have a citation for this? That's probably the most helpful addition to the answer above. "Apparently the Other class in hg19 is comprised of SVA elements. In hg38, these are reclassified as Retrotransposon class and SVA family, which didn't exist in hg19."

I don’t know that there’s anything to cite, this is an observation from looking at the annotations

Ah, you mean you're just looking at the repName for Other and extrapolating these must be SVA elements? You're right---these are SVA_A, SVA_B, SVA_C, SVA_D, SVA_E, SVA_F.

"Hard masked" means we know the sequences but mask them out for analysis purpose. We are unable to assemble most of hg19 centromeres, so we don't know these sequences. I wouldn't use "hard masked" here. We actually don't know hg38 centromeres, either. These are "manmade" sequences. See my answer below.

On (2), centromeres and telomeres primarily consist of satellite sequences. In human centromeres, the dominant type is alpha satellite (aka alphoid). GGAAT micro-satellite (aka SATIII centromeric micro-satellite) is also very common. There are less common types such as beta satellite. Because there are all kinds of centromeric/telomeric sequences of different lengths and properties, RepeatMasker doesn't put them into one class. Note that RepeatMasker has a "centr" class but it does not include all centromeric satellites (no SATIII, for example). In general, it is not advised to find centromeres with RepeatMasker output.
GRCh38 is an unusual reference genome in that its centromeric sequences are computationally generated modeled after Venter as I remember. They are mostly alpha satellite, but are not real sequences (which has always bugged me). GRC know exactly where they put these sequences. GRCh37 is a more typical reference genome. It lacks probably majority of alpha satellites, but still has many. When you map reads to GRCh37, you occasionally see ~10,000-fold read depth even though the average depth is ~30X. These regions are mostly alpha satellites. Mappers just pile all different copies of satellites onto the ones remained in GRCh37.

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


text:
Error in creating a volcano plot in MATLAB
I am a complete newbie to MATLAB and the first task I have is to create a volcano plot. I have been using this link to understand about it and get started.
I tried to run it on dummy values -
a=[1 2 3]

b=[4.6 2.7 4.5]

c=[0.05 0.33 0.45]

And then I ran -
SigStructure = mavolcanoplot(a, b, c)

My understanding is that a represents the gene expression values for condition 1, b for condition 2, and c is the list of p-values for the 3 data points in a and b.
However running this code gives me the error -
Index exceeds matrix dimensions.

Error in mavolcanoplot (line 127)
appdata.effect = X(paramStruct.goodVals) - Y(paramStruct.goodVals);

Error in volc (line 4)
SigStructure = mavolcanoplot(a, b, c)

Can anyone explain where I am going wrong?

Cross-posted on biostars.

I'm not sure if nowadays matlab is still a good option for informatics. I would rather go for a R script to make a volcano plot. If you don't want to have fancy code but functional one, you can even make the foldchange calculations and log2 transformation in excel and just plot the log2 fold change in x axis and the -log10 transform of the p/q values on the y axis.
If you need it, I could provide one of my own scripts.

Update: here is the actual script:
#load data into the 'res' variable
res <- read.table("genes.txt", header=TRUE)
head(res)

# Make a basic volcano plot with grey dots
with(res, plot(log2FoldChange, -log10(qvalue), pch=19, col="#4A4E4D", xlab="Log2 Fold Change", xlim=c(-8,8)))

# Add colored points: "indianred" if padj<0.0001, orange of log2FC>2, green if both)
with(subset(res, qvalue<.0001 ), points(log2FoldChange, -log10(qvalue), pch=19, col="#17B090"))
with(subset(res, abs(log2FoldChange)>2), points(log2FoldChange, -log10(qvalue), pch=19, col="#F89538"))
with(subset(res, qvalue<.0001 & abs(log2FoldChange)>2), points(log2FoldChange, -log10(qvalue), pch=19, col="#FA6367"))

hey Erik, it would be really helpful for me if you could provide me the scripts you mentioned. I am a real noob in this, and am trying to figure out the best way to go about it

so while I am not so familiar with MATLAB, in the definition of the volcano plot it specifies that PValues is supposed to be a column vector or a DataMatrix object. Your c vector is a row vector, and this confuses MATLAB.

Is there a reason why you are doing this in MATLAB instead of R?
In R with the ggplot2 package you can make a volcano plot pretty easily using something like this (may be errors as it is off the top of my head and untested):
# load libraries
library(ggplot2)
library(data.table)

# create a dataframe/data.table of the data
a = c(1,2,3)
b = c(4.6, 2.7, 4.5)
c = c(0.05, 0.33, 0.45)
dt = data.table(A=a,B=b,pval=c)

# add column to calculate log2 fold change between a/b
dt[,log2Diff := (log2(A)-log2(B))]
# add column to transform pvalue to -log scale
dt[,logPval := -log10(pval)]

#make the plot
ggplot(dt, aes(x=log2Diff,y=logPval))+geom_point()

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


text:
How to apply upperquartile normalization on RSEM expected counts?
I see that TCGA RNASeq V2 RSEM data is normalized with upper-quartile normalization.
After doing Quantification with RSEM with the samples I have, I got "genes.results" as output which has gene id, transcript id(s), length, expected count, and FPKM. So, from all the sample output files I got the gene_id and expected_counts [from all samples]. for eg: it looks like following:
Gene_id S1  S2  S3  S4  S5
ENSG00000000003 1355    1121    242 709 1071
ENSG00000000005 5   0   0   0   0
ENSG00000000419 839 1345    48  976 1234
ENSG00000000457 429 1803.08 386 1156    628
ENSG00000000460 404 1444    523 479 1260
ENSG00000000938 294 312 93  189 683
ENSG00000000971 3911    4633.92 264 2863    5194
ENSG00000001036 1546    2276    271 1641    2141
ENSG00000001084 1945    2125    490 980 1533
ENSG00000001167 2054    4069.72 3457    2075    2207
ENSG00000001460 153 339 77  221 273

I want to apply "upper-quartile" normalization on this data. But doesn't know which R package I can use and which function?
Is anyone aware about this?
P.S. This is not for differential analysis

You can use the quantile function in base R to get the value of a particular quantile (e.g. 0.75 for the upper quartile). This can then be used as a factor for normalisation: divide the observed expression by this number.
## set up simulated data
data.mat <- floor(10^data.frame(sapply(1:5,function(x){rnorm(100)})));
colnames(data.mat) <- paste0("S",1:5);
data.mat[data.mat > 0] <- floor(data.mat[data.mat > 0]);

It's generally a good idea to first remove genes that have no counts in any samples. Rows can be removed using a logical selection based on the total expression in any given row.
## Only retain rows that have a positive total expression per gene
data.mat <- data.mat[rowSums(data.mat) > 0,];

Two methods of upper quartile determination are presented here, one which is based on the total counts, and one based on the expressed counts (i.e. excluding values with no expression). These both use the quantile function of R. The apply function calculates statistics that are based on a single specified dimension (or multiple dimensions) of the data. The first dimension of a 2D matrix in R corresponds to the row, while the second corresponds to the column. By using 2 in the apply function, data are collected along the second dimension (i.e. the column), and statistics generated from that:
data.quantileAll <- apply(data.mat, 2, function(x){quantile(x, 0.75)});
data.quantileExpressed <-
  apply(data.mat, 2, function(x){quantile(x[x>0], 0.75)});

Expressed counts may be preferred when read counts are low, so it's difficult to distinguish between a transcript with very low expression and a transcript with no expression.
Normalisation is applied by dividing the values by the column-based upper quartile statistic. R will repeat values along the last dimension first (which in this case is the column), so the matrix needs to be transposed first (with columns as genes) to get the arithmetic working properly, with a final transposition to return the orientation back to the original direction (with columns as samples).
data.norm <- t(t(data.mat) / data.quantileExpressed);

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


text:
How to check whether all BAM read contain defined read groups?
I'm trying to investigate whether there are errors within my BAM. After looking at the BAM header to see whether the read groups exists (using samtools, i.e.
samtools view -H file1.bam

Here, it appears that the header includes @RG tags. However, I suspect that the BAM contains reads that either lack the @RG tag or perhaps these read groups are not well-defined, i.e. it looks like some reads in the BAM may be missing an RG tag.
How do I explicitly check whether all reads have an RG tag, and whether these are well-defined?

According to the man page, running samtools stats --split RG <file1.bam> should produce summary statistics separated by read group. If it doesn't produce a list of ungrouped reads, counts/statistics can be compared to running without the --split RG argument.
Here's some example output from a BAM file with combined read groups:
$ samtools stats --split RG both.bam | grep '^SN'
SN  raw total sequences:    2392
SN  filtered sequences: 0
SN  sequences:  2392
SN  is sorted:  1
SN  1st fragments:  2392
SN  last fragments: 0
SN  reads mapped:   2341
SN  reads mapped and paired:    0   # paired-end technology bit set + both mates mapped
SN  reads unmapped: 51
SN  reads properly paired:  0   # proper-pair bit set
SN  reads paired:   0   # paired-end technology bit set
SN  reads duplicated:   0   # PCR or optical duplicate bit set
SN  reads MQ0:  0   # mapped and MQ=0
SN  reads QC failed:    0
SN  non-primary alignments: 0
SN  total length:   2292231 # ignores clipping
SN  bases mapped:   2255361 # ignores clipping
SN  bases mapped (cigar):   1538560 # more accurate
SN  bases trimmed:  0
SN  bases duplicated:   0
SN  mismatches: 254720  # from NM fields
SN  error rate: 1.655574e-01    # mismatches / bases mapped (cigar)
SN  average length: 958
SN  maximum length: 9078
SN  average quality:    20.0
SN  insert size average:    0.0
SN  insert size standard deviation: 0.0
SN  inward oriented pairs:  0
SN  outward oriented pairs: 0
SN  pairs with other orientation:   0
SN  pairs on different chromosomes: 0

The --split command has also produced additional .bamstat files, which can be looked at for the summary statistics for individual read groups. In this case I'm only interested in the "raw total sequences" line:
$ grep 'raw total sequences' *.bamstat
both.bam_minusOnly.bamstat:SN   raw total sequences:    1901
both.bam_plusOnly.bamstat:SN    raw total sequences:    491

There are 1901 sequences from the "minusOnly" group, and 491 sequences from the "plusOnly" group, and I know from the previous statistics that there are 2392 total sequences. 1901 + 491 = 2392, so I know there are no ungrouped reads in this file.

I'm not sure how one extrapolates from this whether each of the reads have a RG tag. Could you extrapolate a bit how I would deduce from these statistics?

I've added additional information to demonstrate how this can be done.

using samjdk and invoking the function getReadGroup()

getReadGroup() returns The SAMReadGroupRecord from the SAMFileHeader for this SAMRecord, or null if 1) this record has no RG tag, or 2) the header doesn't contain the read group with the given ID.or 3) this record has no SAMFileHeader

java -jar dist/samjdk.jar -e 'return record.getReadGroup()==null;'  input.bam

the output will be a sam/bam file with missing/bad RG .

A complete solution to validate a SAM/BAM file would be picard's ValidateSamFile.  This will find many other issues with BAM files.
http://broadinstitute.github.io/picard/command-line-overview.html#ValidateSamFile


text:
Waterman-Smith-Beyer Sequence Alignment
I am having trouble understanding the affine gap penalty in the following example -

I am not sure where the 3 and 4 come from or the 4 and 5 in cells 1,3 and 2,4. I'm not understanding how the affine gap penalty works here. I understand how the linear gap penalty worked in the Needlemen_Wunsch but not this.

The linear gap penalty has a cost which is directly proportional to the length of the gap: the gap costs 1 unit per consecutive missing base.
The simplest case is assigning each unit (= each missing base) a penalty of 1:

g(k) = k

Or more generally, if you want to increase the penalty per base to 2, 3 or n:

g(k) = n * k
  with k = length of the gap, and n = cost of 1 missing base

The affine gap penalty has a cost that follows -as its name conveniently states- an affine function.
The idea behind it, is that any gap in an alignment has a cost with two components:

the "gap opening" penalty: it is a cost payed only once per gap, basically just for existing.
the "gap extension" penalty: equivalent to the linear gap penalty, it costs 1 unit for each consecutive missing base.

The general formula for the cost is then:

g(k) = a + n*k
  with a= opening penalty, n= extension penalty, k= length of the gap

Compared to the linear gap penalty, this will give a small bonus to longer consecutive gaps compared to skipping one nucleotide, matching one, skipping one again etc (because in the case of longer gaps, you only pay the opening penalty once, so the global cost per base decreases).
Cells on your diagram:
From what I understand, you're working with an affine penalty of g(k) = 1 + 1*k. For the [3,4] you're "coming" from the score of 2 on the left (cell 1,2). You either:

match the A (row 1) with G (column 3) and get a penalty of 1 for substitution -> total of previous 2 + new 1 = 3
open a gap to skip the G (column 3) and get a penalty of 2 (1 for opening + 1 because gap length of 1) -> total of previous 2 + new 2 = 4.

So you choose the lesser penalty and go with 3, which is the number selected in the bigger box.
The explanation is the exact same one for the other cell you mentioned.

I believe I understand now. The reason why we can get 3 or 4 is because in the previous cell, the optimal score of 2 could have come from either the left or the diagonal. If it had come from the diagonal, we are starting a new gap and thus we get 2 + (1 + 1) = 4. However, if we come from the left, we are simply adding to the already existing gap so we get 2 + 1 = 3. Is this the correct way to think about it?

The formula your slide shows describes a dynamic programming with a general gap cost. Finding the optimal solution takes $O(n^3)$ time. No one is using the algorithm these days. This is the first time I see an attempt to fill a matrix with this algorithm. It is interesting, but I am not sure how it works. Although affine gap cost is a special case of general gap cost, the algorithm to find the best score is very different, so is the time complexity. With modern affine gap formulations, each cell holds three values:

The overall optimal score (or the optimal score if this cell is a match, which yields an equivalent but different formulation)
The optimal score when the cell ends at a deletion
The optimal score when the cell ends at an insertion

The values in each cell only depends on the up, left and up-left cells in the matrix. You fill the matrix from the top left corner to the bottom right. This is an $O(n^2)$ algorithm, because you traverse each cell once. For more details, you can read this note. I would recommend to learn affine gap cost the normal way.

which algorithms produce local alignments under $O(n^3)$? Are there any R implementations of such algorithms?


text:
Uniprot: Determining protein and organism given a DNA sequence
I would like to find the identity of a protein, together with the organism name and kingdom.
The given DNA sequence is below:
ACTGTTGTTCGGTGATCATCAGTTGTACAACGTCCTAACAACATCACATGCAATGCTTATGATATTCTTC
TTCATCATGCCAGGCACGATGGCAGGACTAGGCAACTTACTAGTGCCATTCCAGATGAGTGTACCGGAGT
TAGTATTCCCAAAGATTAATAACATCGGTATATGATTTTTAGTATGTGGTCTACTTTTGATTACGGGTTC
ATCTTGGATGGAGGAAGGTTCAGGAACGGCCTGAACCGTCTATCCACCACTAGCGCTCACTGCAAGTCAT
AGCGGACTTGCTGTAGATACGTTCATTATCGCATTGCACATGGCCGGTGCAAGCTCCCTTACAGGAAGCA
TCAACCTTATATGTACAATCGCCTATGCCCGCCGTTCACTCATGGCGATGCTGCAGTCATCACTTTATCC
CTGATCCATTACAATCACTGCAGCGTTACTCATAGGAGTTGTGCCTGTGCTAGCAGGTGCTATCACGATG
CTACTCACTGATAGAAGTTGGAGTACCAGCTTCTATGACAGTTCGGCAGGCGGTGATCCTATGTTGTATC
AGCACTTATTCTGGGTGTTTGGGCATCCAGAAGTCTATATCATCATACTTCCAGTATTCGGTATAGTCAG

I've been asked to translate this DNA sequence into the 6 possible frames which are:
5'3' Frame 1
TVVR-SSVVQRPNNITCNAYDILLHHARHDGRTRQLTSAIPDECTGVSIPKD--HRYMIF
SMWSTFDYGFILDGGRFRNGLNRLSTTSAHCKS-RTCCRYVHYRIAHGRCKLPYRKHQPY
MYNRLCPPFTHGDAAVITLSLIHYNHCSVTHRSCACASRCYHDATH--KLEYQLL-QFGR
R-SYVVSALILGVWASRSLYHHTSSIRYSQ
5'3' Frame 2
LLFGDHQLYNVLTTSHAMLMIFFFIMPGTMAGLGNLLVPFQMSVPELVFPKINNIGI-FL
VCGLLLITGSSWMEEGSGTA-TVYPPLALTASHSGLAVDTFIIALHMAGASSLTGSINLI
CTIAYARRSLMAMLQSSLYP-SITITAALLIGVVPVLAGAITMLLTDRSWSTSFYDSSAG
GDPMLYQHLFWVFGHPEVYIIILPVFGIV
5'3' Frame 3
CCSVIISCTTS-QHHMQCL-YSSSSCQARWQD-ATY-CHSR-VYRS-YSQRLITSVYDF-
YVVYF-LRVHLGWRKVQERPEPSIHH-RSLQVIADLL-IRSLSHCTWPVQAPLQEASTLY
VQSPMPAVHSWRCCSHHFIPDPLQSLQRYS-ELCLC-QVLSRCYSLIEVGVPASMTVRQA
VILCCISTYSGCLGIQKSISSYFQYSV-S
3'5' Frame 1
LTIPNTGSMMI-TSGCPNTQNKC-YNIGSPPAELS-KLVLQLLSVSSIVIAPASTGTTPM
SNAAVIVMDQG-SDDCSIAMSERRA-AIVHIRLMLPVRELAPAMCNAIMNVSTASPL-LA
VSASGG-TVQAVPEPSSIQDEPVIKSRPHTKNHIPMLLIFGNTNSGTLIWNGTSKLPSPA
IVPGMMKKNIISIACDVVRTLYN--SPNNS
3'5' Frame 2
-LYRILEV--YRLLDAQTPRISADTT-DHRLPNCHRSWYSNFYQ-VAS--HLLAQAQLL-
VTLQ-L-WIRDKVMTAASP-VNGGHRRLYI-G-CFL-GSLHRPCAMR--TYLQQVRYDLQ
-ALVVDRRFRPFLNLPPSKMNP-SKVDHILKIIYRCY-SLGILTPVHSSGMALVSCLVLP
SCLA--RRIS-ALHVMLLGRCTTDDHRTT
3'5' Frame 3
DYTEYWKYDDIDFWMPKHPE-VLIQHRITACRTVIEAGTPTSISE-HRDSTC-HRHNSYE
-RCSDCNGSGIK--LQHRHE-TAGIGDCTYKVDASCKGACTGHVQCDNERIYSKSAMTCS
ER-WWIDGSGRS-TFLHPR-TRNQK-TTY-KSYTDVINLWEY-LRYTHLEWH--VA-SCH
RAWHDEEEYHKHCM-CC-DVVQLMITEQQ

I ran jobs on Uniprot for all these sequences. Here are the results:
www.uniprot.org/blast/uniprot/B20171003A7434721E10EE6586998A056CCD0537E082E0FZ
www.uniprot.org/blast/uniprot/B2017100383C3DD8CE55183C76102DC5D3A26728BEBF5BAE
www.uniprot.org/blast/uniprot/B20171003A7434721E10EE6586998A056CCD0537E082E5FP
www.uniprot.org/blast/uniprot/B20171003A7434721E10EE6586998A056CCD0537E082F08D
nothing
www.uniprot.org/blast/uniprot/B20171003F725F458AC8690F874DD868E4ED79B88CE32C4B

But how can I decide the identity of the protein and the name of the organism?
In the first one, the organism is Pacific ridley turtle  but in the second one it is Diplonema sp. ATCC 50224. So I can only decide at a kingdom level, but not to the organism, right?

Hi snnlaynnkrdsm, I've edited your question to make it sound less "homeworkey", as I suspect that would put people off answering it.

To be clear, homework questions are fine on this site. But you should be clear when something is a homework assignment. Accordingly, we often won't provide answers to homework questions outright, but lead you in the right direction.

Well, let's see. The first step is to simply run a nucleotide blast (blastn) against the nr (non-redundant) database. I did this using the NCBI blast tool. While, as expected, I got many hits, there was only one hit which matched the entire length of the query sequence perfectly:
Diplonema sp. ATCC 50224 cytochrome c oxidase subunit 1 (cox1) mRNA, complete cds; mitochondrial
Sequence ID: JF698657.1
Length: 1612Number of Matches: 1
Range 1: 141 to 770
Alignment statistics for match #1
Score   Expect  Identities  Gaps    Strand
1164 bits(630)  0.0 630/630(100%)   0/630(0%)   Plus/Plus
Query  1    ACTGTTGTTCGGTGATCATCAGTTGTACAACGTCCTAACAACATCACATGCAATGCTTAT  60
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  141  ACTGTTGTTCGGTGATCATCAGTTGTACAACGTCCTAACAACATCACATGCAATGCTTAT  200

Query  61   GATATTCTTCTTCATCATGCCAGGCACGATGGCAGGACTAGGCAACTTACTAGTGCCATT  120
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  201  GATATTCTTCTTCATCATGCCAGGCACGATGGCAGGACTAGGCAACTTACTAGTGCCATT  260

Query  121  CCAGATGAGTGTACCGGAGTTAGTATTCCCAAAGATTAATAACATCGGTATATGATTTTT  180
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  261  CCAGATGAGTGTACCGGAGTTAGTATTCCCAAAGATTAATAACATCGGTATATGATTTTT  320

Query  181  AGTATGTGGTCTACTTTTGATTACGGGTTCATCTTGGATGGAGGAAGGTTCAGGAACGGC  240
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  321  AGTATGTGGTCTACTTTTGATTACGGGTTCATCTTGGATGGAGGAAGGTTCAGGAACGGC  380

Query  241  CTGAACCGTCTATCCACCACTAGCGCTCACTGCAAGTCATAGCGGACTTGCTGTAGATAC  300
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  381  CTGAACCGTCTATCCACCACTAGCGCTCACTGCAAGTCATAGCGGACTTGCTGTAGATAC  440

Query  301  GTTCATTATCGCATTGCACATGGCCGGTGCAAGCTCCCTTACAGGAAGCATCAACCTTAT  360
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  441  GTTCATTATCGCATTGCACATGGCCGGTGCAAGCTCCCTTACAGGAAGCATCAACCTTAT  500

Query  361  ATGTACAATCGCCTATGCCCGCCGTTCACTCATGGCGATGCTGCAGTCATCACTTTATCC  420
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  501  ATGTACAATCGCCTATGCCCGCCGTTCACTCATGGCGATGCTGCAGTCATCACTTTATCC  560

Query  421  CTGATCCATTACAATCACTGCAGCGTTACTCATAGGAGTTGTGCCTGTGCTAGCAGGTGC  480
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  561  CTGATCCATTACAATCACTGCAGCGTTACTCATAGGAGTTGTGCCTGTGCTAGCAGGTGC  620

Query  481  TATCACGATGCTACTCACTGATAGAAGTTGGAGTACCAGCTTCTATGACAGTTCGGCAGG  540
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  621  TATCACGATGCTACTCACTGATAGAAGTTGGAGTACCAGCTTCTATGACAGTTCGGCAGG  680

Query  541  CGGTGATCCTATGTTGTATCAGCACTTATTCTGGGTGTTTGGGCATCCAGAAGTCTATAT  600
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  681  CGGTGATCCTATGTTGTATCAGCACTTATTCTGGGTGTTTGGGCATCCAGAAGTCTATAT  740

Query  601  CATCATACTTCCAGTATTCGGTATAGTCAG  630
            ||||||||||||||||||||||||||||||
Sbjct  741  CATCATACTTCCAGTATTCGGTATAGTCAG  770

So, we now know the sequence is from the gene for mitochondrial Cytochrome C Oxidase subunit 1 from Diplonema sp. ATCC 50224. That the sequence is mitochondrial, is very important since the mitochondrion of this species has a slightly different genetic code. If you go to the NCBI's page for this sequence (JF698657.1), you will see this line:
/transl_table=4

That tells you it is using translation table 4:
4. The Mold, Protozoan, and Coelenterate Mitochondrial Code and the Mycoplasma/Spiroplasma Code (transl_table=4)

    AAs  = FFLLSSSSYY**CCWWLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG
  Starts = --MM------**-------M------------MMMM---------------M------------
  Base1  = TTTTTTTTTTTTTTTTCCCCCCCCCCCCCCCCAAAAAAAAAAAAAAAAGGGGGGGGGGGGGGGG
  Base2  = TTTTCCCCAAAAGGGGTTTTCCCCAAAAGGGGTTTTCCCCAAAAGGGGTTTTCCCCAAAAGGGG
  Base3  = TCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAG

See how all of your translated reading frames have in-frame stop codons? That's because you translated them using the standard genetic code instead of the correct mitochondrial one. As you can see above, in this code, the codon TGA is not a STOP, as in the standard code, and instead codes for tyrosine (W).
So, if you look at the result of the protein blast you ran on UniProt, you will see:
Query1   LLFGDHQLYNVLTTSHAMLMIFFFIMPGTMAGLGNLLVPFQMSVPELVFPKINNIGI-FL   59
         LLFGDHQLYNVLTTSHAMLMIFFFIMPGTMAGLGNLLVPFQMSVPELVFPKINNIGI FL
F5CI61   LLFGDHQLYNVLTTSHAMLMIFFFIMPGTMAGLGNLLVPFQMSVPELVFPKINNIGIWFL   107

Query60  VCGLLLITGSSWMEEGSGTA-TVYPPLALTASHSGLAVDTFIIALHMAGASSLTGSINLI   118
         VCGLLLITGSSWMEEGSGTA TVYPPLALTASHSGLAVDTFIIALHMAGASSLTGSINLI
F5CI61   VCGLLLITGSSWMEEGSGTAWTVYPPLALTASHSGLAVDTFIIALHMAGASSLTGSINLI   167

Query119 CTIAYARRSLMAMLQSSLYP-SITITAALLIGVVPVLAGAITMLLTDRSWSTSFYDSSAG   177
         CTIAYARRSLMAMLQSSLYP SITITAALLIGVVPVLAGAITMLLTDRSWSTSFYDSSAG
F5CI61   CTIAYARRSLMAMLQSSLYPWSITITAALLIGVVPVLAGAITMLLTDRSWSTSFYDSSAG   227

Query178 GDPMLYQHLFWVFGHPEVYIIILPVFGIV                                  206
         GDPMLYQHLFWVFGHPEVYIIILPVFGIV
F5CI61   GDPMLYQHLFWVFGHPEVYIIILPVFGIV                                  256

The only locations where there isn't a 100% identical match is where a W is aligned to a STOP (- in the Uniprot output). So, your actual protein sequence is the 2nd frame of the 5' to 3' direction which, when translated using the correct code, is:
LLFGDHQLYNVLTTSHAMLMIFFFIMPGTMAGLGNLLVPFQMSVPELVFPKINNIGIWFL
VCGLLLITGSSWMEEGSGTAWTVYPPLALTASHSGLAVDTFIIALHMAGASSLTGSINLI
CTIAYARRSLMAMLQSSLYPWSITITAALLIGVVPVLAGAITMLLTDRSWSTSFYDSSAG
GDPMLYQHLFWVFGHPEVYIIILPVFGIV

Well if it is not a homework assignment, I would suggest to use blastx. Then you don't need to translate your DNA into protein first, but that will be covered by blastx.
If it is a homework assignment though, the teacher probably wants you to select the best blast hit. So which hit gives the highest score, and most matching (identical) amino acids between the alignment.


text:
Dealing with indels when converting a position to chromStart and chromEnd for a bed file
Starting with a chromosome # and position, I am trying to get chromStart and chromEnd values for a .bed file, but I am not sure how to calculate chromEnd when I have a variant that is an insertion or deletion.
My file is not a VCF file but it gives the chromosome number, chromosomal position (pos), reference (ref) allele, and alternate (alt) allele for each genetic variant, like you would find in a VCF file.  The tricky part is that not all of my variants are single nucleotide base changes--I also have insertions and deletions that can be multiple bases long.
For SNPs I am setting chromStart = pos - 1, and chromEnd = start + 1 = pos.  (Note, .bed format requires 0-based indexing and the end position is exclusive. See https://genome.ucsc.edu/FAQ/FAQformat.html)
I am unsure of how to handle indels though (insertions and deletions).  Currently I am adding the length of the reference allele to the starting position to get the end position (chromEnd = start + len(ref)).  I am not sure if this is right though.
I will be using this .bed file with a .bed file from Gencode and bedtools. Does anyone know the standard way the research community (or at least Gencode) treats indels when determining chromosomal position for .bed files? Thanks!

Usually people don't convert such things to BED files. bedtools understands VCF format, so why not convert to that, since it seems most of the information is already there?

Could you explain why you are trying to make this into a bed file? Is it for displaying it in a genome browser? I don't really see how a bed file could even show an indel, unless it were by using the coordinates of the reference genome.

Hi, @terdon I am trying to make this into a BED file so that I can find which variants overlap exonic areas in the genome. I have a BED file from Gencode listing protein coding regions, and I want to know which of my variants overlap with the exome.

Hi @Devon Ryan, if the best thing to do would be to convert it to VCF format, will the bedtools intersect command be able to compare a VCF file with a BED file, or would I need to convert the Genecode BED to VCF?

bedtools is able to compare BED and VCF files, you don't need to convert anything at that point.

The best way to determine which variants affect exomes is to run through an annotation tool like the VEP (http://www.ensembl.org/info/docs/tools/vep/index.html) for which you need a VCF.

You can convert to VCF and use BEDOPS vcf2bed --insertions, vcf2bed --deletions and vcf2bed --snvs to segregate the VCF into three subsets.
For convenience, BEDOPS conversion tools deal with index shifts, including vcf2bed. The resulting BED file is 0-indexed, half-open.

Thanks @Alex Reynolds.  Is there any way to do this but without segregating the VCF into subsets?

Coming back year later I can say genomic coordinate systems can be tricker than I ever imagined years ago when I wrote this question, especially when converting between different formats lol. Here is a good article explaining 1 based vs 0 based positions. https://www.biostars.org/p/84686/ . VCF format is generally 1 based and BED format is generally 0 based. But also some programs use 1-ish or 0-ish indexing so it gets confusing (looking at you VarSeq ;) )


text:
R Biostrings pairwiseAlignment to BAM
The R package Biostrings has a function to create a pairwiseAlignment from pattern and subject sequences.
So far I can save the result into a text file using writePairwiseAlignments. I would like to save the result into a SAM/BAM file, but couldn't find anything from the BioStrings package or any R package.
What are the steps to convert the result from pairwiseAlignement into a SAM/BAM file?
Example:
library(Biostrings)

patterns = DNAStringSet(c(
  'GCTCCCTGGCGGTAAGTTGATCAAAGGAACGCAAAGTTTTCAAG',
  'ACTTCACAGCTCTGGCGGTAAGTTGATCAAAGGAAACGCAAAGTTTTCAAG'
))

subject = DNAString('ACTTCACCAGCTCCCTGGCGGTAAGTTGATCAAAGGAAACGCAAAGTTTTCAAG')

pwa = pairwiseAlignment(patterns, subject)

writePairwiseAlignments(pwa) # produces below
#=======================================
#
# Aligned_sequences: 2
# 1: P1
# 2: S1
# Matrix: NA
# Gap_penalty: 14.0
# Extend_penalty: 4.0
#
# Length: 54
# Identity:      44/54 (81.5%)
# Similarity:    NA/54 (NA%)
# Gaps:          10/54 (18.5%)
# Score: 27.19725
#
#
#=======================================
# P1 ---------GCTCCCTGGCGGTAAGTTGATCAAAGGAA-CGCAAAGTTTTCAAG
#             ||||||||||||||||||||||||||||| |||||||||||||||
# S1 ACTTCACCAGCTCCCTGGCGGTAAGTTGATCAAAGGAAACGCAAAGTTTTCAAG

Essentially, I would like to convert the result from writePairwiseAlignments to a sam/bam format so I can view it in a GUI application. If I understand correctly, subject (referred from the Biostrings package) should be the reference, the patterns are the reads.

This really sounds like an XY problem. Why don't you just use a GUI alignment viewer? The sam format is really not designed for this.


text:
Using large databases with BLAT
I'm a computer scientist working with biologists at a small school that doesn't have dedicated bioinformatics staff. I apologize if I use incorrect terminology since I have limited bioinformatics background.
Our biology staff is studying a soil sample that has been sequenced. The resulting sequences (in the hundreds of millions) were analyzed using mg-rast, but they didn't seem to have a high confidence in the results.
I've been asked to help analyze these sequences. I setup a local instance of NCBI blast, but it was far too slow to analyze that many sequences. We moved on to using BLAT. They've requested we use the nt database from NCBI. I converted the nt FASTA file to 2bit file, splitting the files into small chunks.
gfServer seems to have limitations however. When I run gfServer, I get the following error while gfServer attempts to start.

Exceeding 4 billion bases, sorry gfServer can't handle that.

Is it possible to use the full nt database using BLAT? Is there an alternative to BLAT that can handle that many sequences?

What do they actually want to do with the sequences? See the closest match? All the matches? Something else entirely?

See the closest match for each sequence, then generate a summary statistics about the sample's environment.

The closest match can be done with STAR/hisat2/bowtie2/bwa or any other standard aligner that can handle large reference fasta files. That would prove vastly faster.

If you want to use BLAT, you could split by searches by chromosome.

Why use the entire nt database? If you're analyzing a soil sample, you already know you don't care about any mammalian sequences, for instance. You can make the search much faster if you choose a subset of nt instead of the entire thing. Especially since your biologists are probably only interested in unicellular species.

If you want a drop in placement for blast that's just much faster, try DIAMOND

I'd recommend having a go with either OneCodex (commercial / Freemium) or Centrifuge.
Chris Mason's lab has recently looked at metagenomic classification and OneCodex gave the most accurate result.
Centrifuge uses Burrows-Wheeler lookups, but hasn't yet been compared by the Mason Lab.


text:
Gene expression signature (AUC) from gene expression data with 48 +48 samples
Does anyone have a good method to find a signature where you combine the expression from multiple genes to predict a specific condition.
Assume you have performed sequencing of 48 cancer samples and 48 normal sample. Then find a combination of genes that predicts the cancer samples, using for instance AUC.

One approach is to use Elastic Net regression, I've seen lots of examples in the literature, though the methods a pretty thin in terms of implementation

I think that with cancer that is quite diverse and with this low number of samples you reproduce whatever gene signature in a larger cohort of samples. (Calculate the power of your study)

Rather than going with a specific (bioinformatic) package you could turn to more general machine learning approaches. For similar problems of mine (with much fewer samples), default implementations of Random Forests (which will also provide importance of features / genes) often yield perfect or near-perfect classification. Note that Random Forests have many favorable characteristics for working with genes, and also form the mechanistic basis of the best-performing (competition-winning) tools for the prediction of gene regulatory networks.

I have a method for doing this that uses a MCMC algorithm originally designed for population genetics (via the structure program) to generate a statistic between 0 and 1 indicating genetic risk. This needs the use of reference populations for cases and controls that the query individuals are added to for classification. I use T1D as an example. My methods paper is not peer-reviewed, and probably doesn't go into sufficient detail for what you are looking for, but at least indicates one way this can be done:
https://www.biorxiv.org/content/early/2017/01/31/104497
Based on the work I've done with this, using 48 cases and 48 controls with a very specifically defined type of cancer should be sufficient to create a generalised signature from human populations. If a generic cancer profile is needed, the populations are not representative of a general human population, or there are some cases in the control group, then things get trickier and require more individuals. On the other hand, 1500 per group is probably overkill.

There are many different ways of doing it! I would recomend the well established R pacakge CARET which have a whole chapter and many build in functions for it.
In addition you probably want to do some low-level filtering based on variance and fold changes (or the combination - ad DE analysis) to get the number of features down before progressing to the more advanced methods.


text:
Differences between NetAffx Hg-U133 Plus 2.0 Annotation file versions
From what I was able to observe, the annotation files (CSV) for NetAffx Hg-U133 Plus 2.0 is updated regularly.
However, from glancing at the CSV files (NA35: 2015 & NA22: 2007), I can't seem to notice any glaring difference between the two.
Does anyone know what information is actually being updated? I assume it may be gene identifiers for different databases, but wanted to make sure.
Thanks.
Below are the files I downloaded and examined from Affymetrix website (you need an account to download archived annotation files)

HG-U133_Plus_2 Annotations, CSV format, Release 22 (20 MB, 5/31/07)
HG-U133_Plus_2 Annotations, CSV format, Release 35 (33 MB, 4/16/15)

After checking the files myself, it does seem like gene identifiers indeed have been updated. Some that were missing have been added with identifiers and some have additional identifiers now.
Among some other things that were updated from 2007 include:

Genome version NCBI Version 36.1 => GRCh37
Alignments
Gene title
Gene symbol

Thought this would be a good reference for those who may be looking into older researches involving microarrays. I will be sure to check the files more carefully before posting questions like this one.

Nice you could answer yourself the question! (you can mark it as solved yourself) But I would still use the most up to date annotation of microarrays when using "old" data. However, I agree that some differences between analysis might be due to using different annotations.


text:
How to obtain clusters of hierarchical heatmap when using Python?
Is there a good way of obtaining the labels (e.g. genes) within individual clusters that haven been clustered hierarchically in Python (preferentially, but not necessarily, by seaborn)?
I found these difficulties highlighted elsewhere (here, and here), and would hope for some simple option (e.g.: like in MATLAB where cutoff can be specified, and corresponding genes would be returned), without potential risk that the plotted hierarchical heatmap differed from the cluster assignment ( see here)

This question needs more information; it's not specific enough. Python is a programming language, not a toolkit. At the moment, it's similar to asking "Is there a good way to find the names of common taxi companies that are found in cities in Europe (preferentially, but not necessarily, French cities)?"

@gringer; either solution or library would be fine as long as it can be called from Python and does the above job and works well (which would exclude for me solutions that go through Python executing compiled local code through bash, or calling local instance of R or MATLAB); perhaps: "is there at least one country in Europe, where it is possible to order a taxi (preferentially, but not necessarily, by phoning the taxi company)?"

Clustering like this is typically done with scipy. Here's the code we use in deepTools (original context here):
from scipy.cluster.hierarchy import fcluster, linkage
Z = linkage(some_matrix)  # You might want to set `method` and `metric`
groups = fcluster(Z, nGroups)  # You might want to set `criterion`

groups is then a vector containing the group assignment of the input matrix according to the clustering. As an aside, this is exactly the same as what's shown in the last link in your post, so I'm not sure how you didn't find it.
BTW, at least in python, hierarchical clustering is painfully slow with large datasets.

thanks; I had indeed been wondering about Python's usual libraries for clustering, including scipy, being awfully slow for hierarchical clustering


text:
How can I calculate the silent mutation for each position in codon?
I am trying to find the number of silent mutations in for each position in codon.
In the homework, I am asked:

For the codons and for each position, count the number of silent mutation. There are 64 codons, for each position there are 3 mutations leading to a change at the DNA level. There are thus 192 possible mutations for each of the three positions of the codons. For each position, print the number of silent mutations.

I understand stand silent mutation is the mutation that doesn't lead to amino acid change. For example, GCC can become GCT but it is still Thr
I am not sure what does "calculating silent mutation for each position" mean? I mean it is easy to find silent mutation for each codon like the example I gave above, but couldn't really understand how to calculate the silent mutation for each position.
Should I change letters (nucleotides) one by one and see if it causing a silent mutation?

Calculate in this context would just mean determine. You don't actually have to change every base to determine whether a change is silent. Instead, try building a hash table of amino acids to codons, which you can then quickly iterate over to find all silent changes.
BTW, try to avoid loaded terms like "mutation", even if your professor incorrectly used it. Terming things like this as "silent polymorphism" or "silent change" or something along those lines is preferable.

but he wants me to find number of silent mutations for each position of codon. Lets say codon is GCT and wants me to count the number of silent mutation for the codons and for each position.

You can efficiently get that from the method I described.

@snnlaynnkrdsm try and implement this. Once you have something written down, if you're stuck with a specific aspect of your program, you can ask again here. But first try it out! We won't just write the program for you, but we're happy to help with specific problems you face.

Curious why you think 'mutation' is a loaded term? To me, 'silent' is the loaded part, since synonymous changes may have fitness effects due to codon usage bias or regulatory element overlapping coding sequences

Patients (and scientists in my experience) tend to over interpret what a "mutation" is. Things like synonymous polymorphisms were less likely to be conveyed to non-professionals, but yes the "silent" terminology is also misleading, though it's not going to make some random patient think they're going to be the next Magneto.


text:
How to measure or assign hydrophobicity score values to individual amino acids of a PDB structure?
I want measure the hydrophobicity of each amino acid within a PDB structure file.  Since I have the PDB file I want to consider the 3D information, rather than sequence-only measures such as GRAVY.
I am not interested in measuring the SASA (solvent accessible surface area) or RD (residue depth) - although correlated, it is not the information I require.
How can I comparatively measure the hydrophobicity of the amino acids within a PDB structure file?  Is there an algorithm, library or program to do so?

Dan is correct but it doesn't answer the question.  After much thought, I have finally figured out how to do it, only thanks to Dan's answer.
The problem is that a subsequence is one-dimensional and does not therefore include any of the additional information which is extractable from the three-dimensional structure (when available).  Obviously, having the three-dimensional information is much more informative and descriptive for building classifiers, predicting interactions, interaction stability, or other purposes.

Find all intramolecular contacts within a radius of N. for the Amino Acid subsequence of interest (such as an Interface).  I set N to, e.g. 4 Angstroms.
Calculate the descriptive statistics based on the set of hydrophobicity values for all intramolecular contact atoms or amino acids (depending on the granularity).  e.g. Make an average of all of the hydrophobicity scores from the selected scale.  (I use GRAVY, although there are many others).

To improve this further, the distance value of N can be set more intelligently than as a static value.   For example, to the same as the cavity radius.

No, this doesn't makes sense (I think).
The hydrophobicity of an AA is the hydrophobicity of the AA, in whatever context, just like it's molecular weight. It's weight is it's weight, in a structure or in isolation.
I know this seems physically unrealistic, but so is an abstract measure of 'hydrophobicity', it's just the best we can do based on partition coefficients.


text:
Downloading CopyNumber450k package from Bioconductor
I am trying to download CopyNumber450k package in the R processor on my laptop in order to see if I could do a step of protocol for my research proffessor.  I was unable to successfully download the package I think (because it was too outdated for the newest version of R?) The version of R I have downloaded is 3.4.1. This package appears to work on the version 1.8.0 of R. I thought this was strange, as that seems really old. I then tried to see if this download was possible using some sort of e-version of an R-processor : https://rdrr.io/bioc/CopyNumber450k/   Even the e-version appeared to be unable to find and download this package of R. I have a feeling that I may be misinterpreting something? I tried to find an equivalent package that could run on the newer version of R and didn’t seem to find anything. Then, I tried to see if it was possible to download an old version of R. I couldn’t find a clear way to do that either. Does anyone on bioinformatics stack exchange have insight regarging what I may be doing wrong, something i may be misinterpreting, or how I may overcome this challenge?
Thanks

This is more of a troubleshooting question than a bioinformatics question. It would be better to search for your specific error on Google, than post a question here. You can probably find something on StackOverflow, or an R-specific site

I think you need to explain better what you have tried, including some code. You say that you were unsuccessful to download it, because of an outdated version of R. How did you do that?
I am not sure what you have tried, but if you read the instructions for installation here, you need to type the following in R console:
## try http:// if https:// URLs are not supported
source("https://bioconductor.org/biocLite.R")
biocLite("CopyNumber450k")

Show us what kind of error you get from this. Then we can maybe help you further.


text:
Read sorted and indexed BAM files faster in C++?
I have some sorted and indexed alignment large BAM files. I'd like to read all the reads, and apply custom operation on each of them. The orders aren't important. I'm using htslib in C++ to do the reading. I have a single powerful machine.
The problem is ... reading a large BAM file is very slow even in C++.
Q: Is there any trick to speed up the reading? Reading BAM files multithreaded (I have many cores)?

If you're reading the entire file, indexing won't help. That's for quickly seeking to specific regions in a BAM file.

@gringer Can I speed up by splitting the file into partitions, then each of those by multithreading? I have many cores to do that.

Reading data is an I/O problem, not a CPU one. Do you have your data on an SSD?

@AlexReynolds Yes. But is there anything else we can do in programming?

Using mmap routines to read data, maybe? It might give you a slight benefit over reading with fopen/fseek/fread etc. But there's no way to guarantee a benefit, short of benchmarking whatever it is you are trying to do. And you have to be careful about purging caches between tests. It's easy for people to make spurious claims about their tools being fast if they aren't purging the read cache.

You might investigate terms like "striping" or "parallel I/O" to gain benefit from multiple read operations via hardware. Having multiple processors or cores can be useful for algorithms, but I don't think it will benefit reading sequential data from storage.

The optimal solution is entirely determined by the "custom operation" you want to apply. What do you want to do exactly?

What do you mean by "very slow"? How does it compare with time cat <file.bam> > /dev/null? How does it compare to the code when nothing is done with the records?

In this answer I want to show you some benchmarks that compares three different serial ways of reading data in C++, the third one being the fastest.
In the first method I use an std::istreambuf_iterator, in the second I read the file line by line into a std::vector and in the third one I use C flavoured operations.
/*
compile with g++ -std=c++11 -O2 reading_files.cpp -o reading_files
test the program with the following
dd if=/dev/random of=./dummy5GiB bs=1024 count=$[1024*1024*5]
 for i in 0 1 2; do time ./reading_files dummy5GiB $i; done
*/

#include <fstream>
#include <iostream>
#include <vector>
#include <cstring> //for strcmp
#include <string>
#include <algorithm> //for min

using namespace std;

int main(int argc, char * argv[])
{
    if(argc != 3) {
        cerr << "usage: reading_files <path> <number 0 or 1 or 2>\n";
        return 1;
    }

    ifstream file(argv[1], ios::in|ios::binary);
    if(!file.is_open()) {
        cerr << argv[1] << " file not opened!\n";
        return 1;
    }
    cout << "---------------------------\n";
    if(!strcmp(argv[2],"0")) {
        cout << "method 0\n";
        vector<char> content((std::istreambuf_iterator<char>(file)), std::istreambuf_iterator<char>());
        cout << content.size() << " bytes read\n";
    } else if(!strcmp(argv[2],"1")) {
        cout << "method 1\n";
        string line;
        vector<string> content;
        while(getline(file,line)) {
            content.push_back(line);
        }
        cout << content.size() << " lines read\n";
    } else if(!strcmp(argv[2],"2")) {
        cout << "method 2\n";
        file.seekg(0, std::ios::end);
        unsigned long long file_size = file.tellg();
        cout << "file size is " << file_size << "\n";
        //maximum 5 GiB of memory, so to fit my RAM
        unsigned long long max_block_size = 1024ULL * 1024ULL * 1024ULL * 5;
        file.clear();
        file.seekg(0, std::ios::beg);

        unsigned long long block_size = min(max_block_size,file_size);
        char * buffer = (char*)malloc(block_size*sizeof(char));
        unsigned long long processed = 0;

        while(!file.eof()) {
            file.read(buffer, block_size);
            cout << "read a block of " << min(file_size - processed, block_size) << " bytes\n";
            //do work here with the block in memory (if the block size is > 0)
            processed += block_size;
        }

        free(buffer);
    }
    file.close();
    return 0;
}

In the third method I read the file in blocks of 5GiB, so to fit the free RAM of my machine.
Here is the benchmark of the program on a 5 GiB randomly generated file.
bash-3.2$ for i in 0 1 2; do time ./reading_files dummy5GiB $i; done
---------------------------
method 0
5368709120 bytes read

real    0m35.194s
user    0m19.665s
sys 0m12.999s
---------------------------
method 1
20973135 lines read

real    0m40.478s
user    0m34.046s
sys 0m5.870s
---------------------------
method 2
file size is 5368709120
read a block of 5368709120 bytes
read a block of 0 bytes

real    0m4.413s
user    0m1.757s
sys 0m2.617s

Here is the benchmark on a real 32GiB .bam file using only the third method, here the splitting of the file read plays a fundamental role. I can't tell you how big because my hard disk memory gets saturated before the end of the execution and so I have to kill the program.
time ./reading_files ./HG00252.mapped.ILLUMINA.bwa.GBR.low_coverage.20130415.bam 2
---------------------------
method 2
file size is 34316054058
read a block of 5368709120 bytes
read a block of 5368709120 bytes
read a block of 5368709120 bytes
read a block of 5368709120 bytes
read a block of 5368709120 bytes
read a block of 5368709120 bytes
read a block of 2103799338 bytes

real    0m50.656s
user    0m9.135s
sys 0m35.728s

The key point is to minimise the file operations reading large block of data at once (as long as they fit the memory) and operate on them.
Using mmap could lead to an improvement, but I don't think it would be as large as using the third method instead of the first two.

Why is your “method 2” using gratuitous manual (C-style!) memory management instead of a vector? — I’m also not sure how much this actually helps OP since the question is about parsing BAM, and as far as I know neither Staden IO nor htslib supports reading from a custom
buffer (but should perform its own buffering anyway).

If you need to iterate through the entire BAM file, decompression is usually the main bottleneck.  Three things you can do:

Install libdeflate (https://github.com/ebiggers/libdeflate ) and then the development version of htslib, and make sure the latter is compiled to use libdeflate.  This roughly doubles decompression efficiency in my experience.
Call htslib's bgzf_mt() function right after opening the file for reading, to request multiple decompressor threads.  I generally go with 2-6, depending on what I'm actually doing with the reads.
Have a separate read/decompress-ahead thread perform the bam_read1() calls.  Ideally, when your main thread is processing read #X, background threads are already working on decompressing reads e.g. (X+100) through (X+200).

Also, if your workload is "embarrassingly parallel", you may also benefit from launching multiple processes simultaneously, configuring each to handle a similarly-sized chunk of the BAM.

I stumbled upon this when looking for ways to improve my BAM processing and I was wondering if you have a code example that you could share?

In both deepTools and MethylDackel we implement multithreading by having individual threads process separate genomic regions. At some point you hit an IO limit, but unless your processing is trivial then IO won't typically be limiting with single-threaded applications. Since you already have sorted/index BAM files, you can take advantage of the same strategy.


text:
Retrieve detailed gene descriptions
Given a list of gene IDs, how do you retrieve the gene description, summary and other detailed information in R?

MyGene.info
Using the R package for the MyGene.info database:
source("https://bioconductor.org/biocLite.R")
biocLite("mygene")

library(mygene)
gene <- getGene("1017", fields="all")
gene$summary

Multiple genes:
gene <- getGenes(c(1347:1350), fields="summary")

Entrez database
library(devtools)
install_github("ropensci/rentrez")

library(rentrez)
genes <- entrez_summary(db="gene", id=c(718,720))
genes$`720`$summary

Nice question (and answer, I didn't know about MyGene), if you know what you are looking for maybe it is worth to mention biomaRt

Description of how to use biomaRt here. Filter by your list of gene IDs. Get the descriptions etc as attributes. Eg:
> library(biomaRt)
> ensembl = useEnsembl(biomart="ensembl", dataset="hsapiens_gene_ensembl")
> IDs <- c("BRCA2","BRAF")
genedesc <- getBM(attributes=c('external_gene_name','description'), filters = 'external_gene_name', values = IDs, mart =ensembl)
> genedesc
  external_gene_name                                                                      description
1               BRAF B-Raf proto-oncogene, serine/threonine kinase [Source:HGNC Symbol;Acc:HGNC:1097]
2              BRCA2                  BRCA2, DNA repair associated [Source:HGNC Symbol;Acc:HGNC:1101]

Thank you for the answer. I did not find a way to return detailed information  using biomaRt. The description fields are just extended gene names.

The easiest and fastest method (which also doesn’t require internet access) is to use annotation packages, e.g. org.Mm.eg.db:
library(AnnotationDbi)
library(org.Mm.eg.db)

# 1:1 mapping symbol→entrez
entrez_ids <- mapIds(org.Mm.eg.db, gene_symbols, 'ENTREZ', 'SYMBOL')

# 1:n mapping symbol→…
print(columns(org.Mm.eg.db))  # check what you can query
df <- select(org.Mm.eg.db, gene_symbols, c(...), 'SYMBOL')


text:
Getting hierarchy of cell populations with Drop-seq data
I plan on collecting some drop-seq data from brain tissue (from mice and humans). Using only the 75-85% of genes that have 1:1 orthology between mice and humans, I'd like to cluster the cells by transcriptomic similarity and get a dendrogram to use for cell type classification. There are tons of tools for processing and clustering single cell RNA-seq data (https://www.scrna-tools.org/tools), but I get the impression that most of them were not designed with Drop-seq in mind. Anyone have any insight into the pros and cons of the best tools for Drop-seq analysis?

What are the particularities in Drop-seq that would require a different tool/method to analyze it? (It is the first time I heard about drop-seq and I am not totally sure I understood the difference between drop-seq and sc-seq)

Compared with other scRNA-seq techniques (like Smart-seq), Drop-seq sacrifices read counts for cell counts. That's the big difference. It's also very cost-effective and has been increasing in popularity since it came out in 2015. http://www.cell.com/molecular-cell/abstract/S1097-2765(17)30049-7

This article uses the freely available R package dropbead for filtering and then Seurat to perform a principal component analysis that groups together affine transcriptomes. It could be what you are looking for.

I've seen that paper and read the abstract, but I didn't know about Dropbead. Thank for the info!

Note that there isn't any license for the package. It might not be safe to use it. (Although I am sure that if you contact the authors they will grant you permission and correct this :)

The DESCRIPTION file states GNU GPL3, but yes it would be good to have a LICENSE file

SC3 is part of the Bioconductor library for R. It is based on PCA and spectral dimensionality reductions and utilizes k-means. In the introductory paper the authors perform analysis on a Drop-seq dataset of more than 44,000 cells. SC3 outputs several useful plots, including a consensus matrix, expression panel, and a marker genes heatmap. All of these utilize the hierarchical clustering.
I can't speak from personal experience but the package seems very user friendly. Being a Bioconductor package means it has lots of support and excellent documentation. I am getting ready to use it myself for a small Drop-seq dataset of olfactory epithelial cells.

SC3: consensus clustering of single-cell RNA-seq data : Nature Methods : Nature Research. Available at: https://www.nature.com/nmeth/journal/v14/n5/full/nmeth.4236.html. (Accessed: 31st October 2017)

You can try a hierarchical clustering tool HGC (hierarchical graph clustering).
You may get cell population at different hierarchy by choosing multiple cutting heights. See the paper for more details.

The current version of Seurat and Scanpy can both handle Drop-seq data. If you want to include more low-reads cells you can loosen the filtering thresholds in either of these packages. These two packages really can handle both low-read-high-cell data and high-read-low-cell data.
This ENCODE paper used Seurat to process both C1 Fluidigm data (>4000 genes per cell) and 10X data (>1000 genes per cell only).
This human limb paper used Seurat to process Drop-seq data, although they didn't annotate the cell types well. I used Seurat to reprocess their data and it works much better than what the paper presents.
I guess the key is parameter tuning.


text:
How can I only get the species name for fasta sequences from blast results?
I am trying to make a phylogenetic tree from sequences obtained with blast. I have files that contains hundreds of fasta protein sequences which are all named like >NP_567483.1 transcriptional regulator [Pseudomonas aeruginosa PAO1]. But if I use these files to do multiple alignments or phylogenetic trees using either ClustalX or MEGA7, every sequence are named after the first ID in the sequence names (for this example it would be NP_567483.1), and what I want is to have the name of the strain instead (here Pseudomonas aeruginosa PAO1) which is part of the sequence description.
So I don't know whether the softwares I'm using understand that the information between the [ ] is the species and that I just didn't find where to set this as the name or if I have to write some kind of a script to change my fasta files before. I'm really new to this and don't know much about it. Thanks for your help in advance. Julian

Should we assume one sequence per species? Otherwise species name is no good as an identifier.

There might be a "cleaner" way to do this, but what I would do is first remove all of the spaces from the fasta sequence names and replace them with underscores using something like sed.
sed 's/ /_/g'

Run your BLAST after you remove the spaces and you will get the whole ID string in the result instead of only that first part (e.g. NP_567483.1). From here you can parse out the species name, or whatever other ID you would like.
Hope this is helpful! :)

You can replace fasta header of each sequence with the string that is inside the brackets using awk:
awk -F '[][]' '{if ($0~/>/) {print ">"$2} else {print $0}}' input.fasta > output.fasta

For lines starting with > (fasta header), it replaces the line with the string inside the brackets and for other lines (which contain the sequence), it prints the line without making any changes.

With sed, first test it:
echo '>NP_567483.1 transcriptional regulator [Pseudomonas aeruginosa PAO1]' | sed "s/>.*\[\(.*\)\]/>\1/"

Then modify your fasta:
sed -i.bak "s/>.*\[\(.*\)\]/>\1/" file.fasta


text:
Error "Can't locate List/MoreUtils.pm in @INC" when I run update_blastdb.pl
I encoutered some technical problem when learning BLAST as below:
After I use conda install blast to install ncbi-blast, and run the command:
update_blastdb.pl --decompress 16SMicrobial
the shell complains:

Can't locate List/MoreUtils.pm in @INC (you may need to install the
  List::MoreUtils module) (@INC contains: /etc/perl /usr/loc
  al/lib/x86_64-linux-gnu/perl/5.22.1 /usr/local/share/perl/5.22.1
  /usr/lib/x86_64-linux-gnu/perl5/5.22 /usr/share/perl5 /usr/li
  b/x86_64-linux-gnu/perl/5.22 /usr/share/perl/5.22
  /usr/local/lib/site_perl /usr/lib/x86_64-linux-gnu/perl-base .) at
  /home/mck f111/miniconda3/envs/bioinfo/bin/update_blastdb.pl line 43.

After I run this command:
cpanm List::MoreUtils

It appears some installation happened, however, the complains are still there.
Could anyone tell me how to fix this? Thanks a lot!

What operating system are you running this on? That looks like some sort of *nix, but is it Linux, macOS, UNIX, Ubuntu subsystem on Windows? Something else? Did you run the cpanm command as a regular user or as root? Did you see where it installed (it should be mentioned in the command's output)? Please [edit] your question and clarify.

Install List::MoreUtils using the same user as blast.


text:
Download multiple SRA files
I want to download all SRA file from the following project. Is there a method to download all the SRA files at the same time?

You could use the R/Bioconductor package SRAdb

@JohnBlischak This might be an answer if expanded explaining how SRAdb helps to answer the question

A quick look at your link tells me the SRR numbers run from SRR837819 to SRR837856. You can use fastq-dump from the sratoolkit, and make a for loop around it in bash.
Something like this should work:
for (( i = 19; i <= 56; i++ ))
  do
  fastq-dump --accession SRR8378$i
done

After reading Devon Ryan's answer, I realize that you asked for SRA files instead of fastq. This can also be done with prefetch of the sratoolkit:
for (( i = 19; i <= 56; i++ ))
  do
  prefetch SRR8378$i
done

Assuming you ultimately just want the fastq files and you know the SRR (run) numbers, I would download them from here: ftp://ftp.sra.ebi.ac.uk/vol1/fastq/
As for downloading multiple files, I've just used multiple wget commands. I don't know of a way to download all of the files together in like a zipped folder or anything :/

I suggest you follow the advice in Eric A Brenner's answer and just download the fastq files. However, if you really really want to use the SRA files for some reason, note that you can use parallel-fastq-dump to make things faster. Do follow its advice regarding using prefetch.
You'd need to combine that with the answer from b.nota (i.e., put the commands in a for loop).

I was able to find a solution for this using Entrez Direct and SRA toolkit :) If you have the project number or SRA project number, yours would be SRP022054 in this case for the 36 SRAs, you can use esearch to make a query like so and pipe it into SRA toolkit with this one liner:
esearch -db sra -query SRP022054  | efetch --format runinfo | cut -d ',' -f 1 | grep SRR | head -5 | xargs fastq-dump  --skip-technical  --readids --read-filter pass --dumpbase --split-3

You can download an entire project using pysradb:
pysradb download -p SRP022054

It retains the same schema as SRA.

Once you install SRA tools, you can grab a BioProject/SRP accession list (Send to - File - Accession List) and save it to default SRA tools /sra directory, then run (with /sra as cwd):
prefetch --option-file SraAccList.txt
cat SraAccList.txt | xargs fasterq-dump --outdir "/your-directory/for-fastq"

P.S. The text file can be any list of accessions, separated by return.
Tested in Bash

This online tool creates download links for the less tech savy:
https://sra-explorer.info/#
Using PRJNA201245 as ID, it provides several codes including a bash script for downloading all files.


text:
Help with identifying disease modules
I've made an application that at this point ranks all combinations of drug pairs relevant to a biological network/graph in the order of how disruptive the outcome of deleting the targets of a given drug pair is. So, in short, my applications objective is to kill/disrupt the network, but the optimal goal would be to disrupt nodes related to disease and have minimal impact on the functional part of the cell/network. I've read a few articles about network medicine and found that a diffusion-based method for identifying disease modules seems easier to understand coming from an informatics background.
I know how to implement the algorithm for random walkers with reset algorithm (RWR) but unsure how to get/ identify known disease associated genes/nodes of the network and these are used as seeds for the algorithm.
If you guys have some advice on how to do this or even know of a better way to identify disease, functional and topological modules I would appreciate it.
Some details of my project:
Application developed to run in Cytoscape, programmed in Java. Using Jung and JGraphT for graph algorithms.
This is for exploring paired drug combinations for treating cancer. I retrieve drug-target relation from a refined list published in an article called: "A comprehensive map of molecular drug targets" from December.