macでインフォマティクス

macでインフォマティクス

NGS関連のインフォマティクス情報についてまとめています。

ターゲット配列とその近傍領域をアセンブリする mapsembler2

 

Mapsembler2は、ターゲットを絞ったアセンブリソフトウェアである。Mapsembler2は、入力として任意のNGSのrawリードセットとスターター配列を取り、ユーザーの選択に応じて、そのスターター配列近傍を線形シーケンスまたはグラフとして出力する。

 

以下のような用途で使うことが提案されている。

1、アセンブルした配列を検証する(SNV、SVsがあるかどうか)。

2、特定の酵素をコードする塩基配列がメタゲノムのリードセット内に存在するか調べる。

3、配列をextendする。

 

紹介

https://www.biostars.org/p/109189/

解説

https://www.biostars.org/p/109189/

HP

https://colibread.inria.fr/software/mapsembler2/

 

音が出ます。

 

 

 インストール

#bioconda(link)
conda create -n mapsembler2 -y
conda activate mapsembler2
conda install -c bioconda -y mapsembler2

mapsembler2_extend

$ mapsembler2_extend 

NAME

mapsembler_extend, version 1.0.0 - Copyright INRIA - CeCILL License

 

SYNOPSIS

mapsembler2_extend <extrem_kmers.fasta> <readsC1.fasta> [<readsC2.fasta> [<readsC3.fasta] ...] [-t extension_type] [-k value] [-c value] [-g value] [-i index_name] [-o name] [-h]

 

DESCRIPTION

TODO

 

OPTIONS

-t extension_type. Default: 1

    1: a strict sequence: any branching stops the extension

    2: a consensus sequence: contiging approach

    3: a strict graph: any branching is conserved in the graph

    4: a consensus graph: "small" polymorphism is merged, but "large" structures are represented

-k size_kmers: Size of the k-mers used duriung the extension phase Default: 31. Accepted range, depends on the compilation (make k=42 for instance) 

-c min_coverage: a sequence is covered by at least min_coverage coherent reads. Default: 2

-g estimated_genome_size: estimation of the size of the genome whose reads come from. 

      It is in bp, does not need to be accurate, only controls memory usage. Default: 3 billion

-x node_len: limit max of nodes length. Default: 40

-y graph_max_depth: limit max of graph depth.Default: 10000

-i index_name: stores the index files in files starting with this prefix name. Can be re-used latter. Default: "index"

    IF THE FILE "index_name.bloom" EXISTS: the index is not re-created 

-o file_name_prefix: where to write outputs. Default: "res_mapsembler" 

-p search_mod: kind of prosses Breadth or Depth. Default: Breadth 

-h prints this message and exit

mapsembler2_extremities

$ mapsembler2_extremities 

USAGE for 'mapsembler2_extremities'

    --k                 (1 arg) :    kmer size that will be used for mapsembler2  [default '']

    --starters          (1 arg) :    starters fasta file  [default '']

    --reads             (1 arg) :    reads dataset file name. Several reads sets can be provided, surrounded by "'s and separated by a space (e.g. --reads "reads1.fa reads2.fa")  [default '']

    --output            (1 arg) :    output substarters file name  [default '']

    --min-solid-subkmer (1 arg) :    minimim abundance to keep a subkmer  [default '3']

    -debug              (0 arg) :    debugging

    -nb-cores           (1 arg) :    number of cores  [default '0']

    -verbose            (1 arg) :    verbosity level  [default '1']

    -help               (0 arg) :    display help about possible options

 

 

実行方法

 スターター配列とペアエンドfastq(gzip圧縮にも対応)を指定する。

mapsembler2_extend starter.fasta pair_1.fq pair_2.fq -t 1 -k 31 -c 2
  • -t    extension_type. Default: 1
        1: a strict sequence: any branching stops the extension
        2: a consensus sequence: contiging approach
        3: a strict graph: any branching is conserved in the graph
        4: a consensus graph: "small" polymorphism is merged, but "large" structures are represented 
  • -k    Size of the k-mers used duriung the extension phase Default: 31. Accepted range, depends on the compilation (make k=42 for instance) 
  • -c    min_coverage: a sequence is covered by at least min_coverage coherent reads. Default: 2

 

引用

Mapsembler2 – ANR Colib'read

 

関連


 

 

 

 

トランスポゾンを分類する TEsorter

 

 Transposable elements(TE)は真核生物ゲノムの重要な部分を構成するが、それらの分類、特にクレードレベルでの分類は依然として困難である。 この目的のために、TEの保存されたタンパク質ドメインに基づいたTEsorterを提案する。 TEsorterはTE、特にLTRレトロトランスポゾン(LTR-RT)を分類し、使いやすく、マルチプロセッシングで高速で、高感度で正確である。 その結果は、分類されたLTR-RTの系統関係と多様性を直接反映することもできる。 Pythonのコードはhttps://github.com/zhangrengang/TEsorterから無料で入手できる。

 

使用されているデータベース

http://repeatexplorer.org/?page_id=918

 

インストール 

ubuntu18.04LTSのpython3.5環境でテストした(conda createで仮想環境作成)。

依存

  • python >3
  • + biopython: quickly install by pip install biopython or conda install biopython
  • + parallel python v1.6.4.4: quickly install by conda install pp
  • hmmscan 3.1x or 3.2x: be compatible with HMMER3/f database format. quickly install by conda install hmmer
  • blast+: quickly install by conda install blast
#bioconda(link)
conda create -n tesorter -y
conda activate tesorter
conda install -c bioconda tesorter -y

本体 Github

git clone https://github.com/zhangrengang/TEsorter 
cd TEsorter
sh build_database.sh

python ../TEsorter.py -h

$ TEsorter -h

2020-01-20 23:15:41,696 -WARNING- Grid computing is not available because DRMAA not configured properly: Could not find drmaa library.  Please specify its full path using the environment variable DRMAA_LIBRARY_PATH

2020-01-20 23:15:41,696 -INFO- No DRMAA, Switching to local/cluster mode.

usage: TEsorter [-h] [--version] [-db {rexdb,rexdb-metazoa,rexdb-plant,gydb}]

                [-st {nucl,prot}] [-pre PREFIX] [-fw] [-p PROCESSORS]

                [-tmp TMP_DIR] [-cov MIN_COVERAGE] [-eval MAX_EVALUE] [-dp2]

                [-rule PASS2_RULE] [-nolib] [-norc] [-nocln]

                sequence

 

lineage-level classification of transposable elements using conserved protein

domains

 

positional arguments:

  sequence              input TE sequences in fasta format [required]

 

optional arguments:

  -h, --help            show this help message and exit

  --version             show program's version number and exit

  -db {rexdb,rexdb-metazoa,rexdb-plant,gydb}, --hmm-database {rexdb,rexdb-metazoa,rexdb-plant,gydb}

                        the database used [default=rexdb]

  -st {nucl,prot}, --seq-type {nucl,prot}

                        'nucl' for DNA or 'prot' for protein [default=nucl]

  -pre PREFIX, --prefix PREFIX

                        output prefix [default='{-s}.{-db}']

  -fw, --force-write-hmmscan

                        if False, will use the existed hmmscan outfile and

                        skip hmmscan [default=False]

  -p PROCESSORS, --processors PROCESSORS

                        processors to use [default=4]

  -tmp TMP_DIR, --tmp-dir TMP_DIR

                        directory for temporary files [default=./tmp]

  -cov MIN_COVERAGE, --min-coverage MIN_COVERAGE

                        mininum coverage for protein domains in HMMScan output

                        [default=20]

  -eval MAX_EVALUE, --max-evalue MAX_EVALUE

                        maxinum E-value for protein domains in HMMScan output

                        [default=0.001]

  -dp2, --disable-pass2

                        do not further classify the unclassified sequences

                        [default=False for `nucl`, True for `prot`]

  -rule PASS2_RULE, --pass2-rule PASS2_RULE

                        classifying rule [identity-coverage-length] in pass-2

                        based on similarity [default=80-80-80]

  -nolib, --no-library  do not generate a library file for RepeatMasker

                        [default=False]

  -norc, --no-reverse   do not reverse complement sequences if they are

                        detected in minus strand [default=False]

  -nocln, --no-cleanup  do not clean up the temporary directory

                        [default=False]

 

 

テストラン

TEsorter-test

f:id:kazumaxneo:20200120232615p:plain

出力(詳細はGithubに書かれている)

f:id:kazumaxneo:20200120232735p:plain


rice6.9.5.liban.rexdb.cls.tsv (TEs/LTR-RTs classifications)

f:id:kazumaxneo:20200120232901p:plain

 

実行方法

TE配列を指定する。8CPU指定。データベースにはREXdb(link)が使用されている。-dbで変更できる。

TEsorter input_file -p 8 -db rexdb

 

 

 ゲノム配列のFASTAはあるがアノテーション情報がない場合、Repeat maskerとLTR_retrieverを使って、TEをゲノムから抽出できます。その1例がGithub READMEの下の方に記載されています。また、LTRなどのRTドメインの配列を抽出、配列比較からphylogenetic analysesを行う例もGithub READMEに記載されています(中盤付近)。

 

引用

TEsorter: lineage-level classification of transposable elements using conserved protein domains

Ren-Gang Zhang, Zhao-Xuan Wang, Shujun Ou, Guang-Yuan Li

bioRxiv preprint first posted online Oct. 10, 2019

 

メタゲノムアセンブリから真核生物由来配列を予測する EukRep

 

 真核微生物は生態系機能の重要な貢献者である。微生物群集の中の真核生物を特定するために遺伝子調査またはDNA「バーコード」が頻繁に使用され、真核生物の多様性の幅が示されている(Pawlowski et al、2012)。ただし、これらのアプローチでは種を検出することしかできず、シーケンスされたゲノムがないと代謝やライフスタイルに関する情報を提供できない。完全に配列決定された真核生物ゲノムの大部分は、培養生物からのものである。遺伝子調査で検出された多様な原生生物といくつかの真菌の培養物アクセスの欠如により、真核生物のリファレンスゲノムデータベースに大きなギャップが引き起こされている(Caron et al、2008; Pawlowski et al、2012)。シングルセルゲノミクスは、未培養の真核生物の配列決定に有望であり、一部の部分ゲノムを生成した(Cuvelier et al、2010; Yoon et al、2011; Monier et al、2012; Vaulot et al、2012; Roy et al、2014; Mangot et al、2017)。ただし、複数の置換増幅は、シングルセルゲノムの完全性を制限している(Woyke et al、2010)。(一部略)

 ここでは、多様な環境サンプルからのデータセットアセンブルされた真核生物配列識別に新しいk-merベースのアプローチを適用した。真核生物ゲノム断片の同定により、ドラフトゲノムへのアサインと遺伝子予測の質の向上が可能になった。アセンブルされたメタゲノムコンティグ上で予測された遺伝子は、系統学的プロファイル、再構築されたゲノムの分類、およびそれらの完全性の評価を取り入れた、さらなるビニング決定のための重要なインプットを提供する。(一部略)
 ユタ州クリスタルガイザーの深部地下微生物群集は、候補栄養門(CP)からの多くの生物を含む、化学合成独立栄養細菌と古細菌に支配されていることがよく特徴付けられている(Probst et al、2014、2016; Emerson et al、2015)。現在の理解では、間欠泉の噴出によってさまざまな新しい種類の細菌や古細菌が地表にもたらされている(Probst et al、2018)。このような深い堆積環境では、有機炭素化合物の可用性が高いとは考えられない。したがって、このシステムへの有機炭素の添加は、新しい間欠泉微生物に対して選択し、よりよく知られている従属栄養生物を濃縮することにより、コミュニティの構成を大きくシフトすると仮定した。この予測をテストするために、浅い間欠泉に追加され、地下水導管で崩壊した木材のサンプルを分析した(以後、CG_WC)。このサンプルと、CG_WCの前日に収集された非木材サンプル(CG_bulk)をメタゲノム解析にかけた。

(以下略)

 

インストール

ubuntu18.04 のpython3.6環境でテストした(docker使用、ホストOS ubuntu18.04LTS)

本体 Github

pip install EukRep

EukRep -h

# EukRep -h

usage: EukRep [-h] -i I -o O [-ff] [--min MIN] [--model MODEL] [-k KMER_LEN]

              [--prokarya PROKARYA] [--seq_names]

              [-m {strict,balanced,lenient}] [--version] [--tie TIE]

 

Identify sequences of predicted eukaryotic origin from a nucleotide fasta file. Individual sequences are split into 5kb chunks. Prediction is performed on each 5kb chunk and sequence origin is determined by majority rule of the chunks.

 

optional arguments:

  -h, --help            show this help message and exit

  -i I                  input fasta file

  -o O                  output file name

  -ff                   Force overwrite of existing output files

  --min MIN             Minimum sequence length cutoff for sequences to be included in prediction. Default is 3kb

  --model MODEL         Path to an alternate trained linear SVM model. Default is lin_svm_160_3.0.pickle

  -k KMER_LEN, --kmer_len KMER_LEN

                        Kmer length to use for making predictions. Lengths between 3-7bp are available by default. If using a custom trained model, specify kmer length here.

  --prokarya PROKARYA   Name of file to output predicted prokaryotic sequences to. Default is to not output prokaryotic sequences.

  --seq_names           Only output fasta headers of identified sequences. Default is full fasta entry

  -m {strict,balanced,lenient}

                        Not compatable with --model.

                                How stringent the algorithm is in identifying eukaryotic scaffolds. Strict has a lower false positive rate and true positive rate; vice verso for leneient. Default is balanced.

  --version             show program's version number and exit

  --tie TIE             Specify how to handle cases where an equal number of a sequences chunks are predicted to be of eukaryotic and prokaryotic origin (Generally occurs infrequently).

                                euk = classify as euk

                                prok = classify as prok

                                rand = assign randomly

                                skip = do not classify

                                Default is to classify as eukaryotic.

 

 

実行方法

FASTA配列を指定する。

EukRep -i input.fasta -o output

 

チュートリアルも紹介されています。

https://github.com/patrickwest/EukRep_Pipeline

引用

Genome-reconstruction for eukaryotes from complex natural microbial communities
West PT, Probst AJ, Grigoriev IV, Thomas BC, Banfield JF

Genome Res. 2018 Apr;28(4):569-580

 

 

GuppyのGPU版を使う

 

タイトルの通り、GuppyのGPU版を使うまでの流れをまとめておきます。

 

ubuntuへのインストール

1、Nvidia GPU driverのインストール

#レポジトリの追加
sudo add-apt-repository ppa:graphics-drivers/ppa
sudo apt update

#NVIDIA driverのインストール。最新GPUだとより最新のNvidiaドライバーを入れる必要があるかもしれない(ONTのGuppy documentより)。
sudo apt install nvidia-384

#OS reboot
sudo reboot

libcuda.so.1がないというエラーが出たら、/libcuda.soから/libcuda.so.1にシンボリックリンクを張って、$LD_LIBRARY_PATHに追加することでとりあえず解決。

#私の環境では
ln -s /usr/local/cuda/lib64/stubs/libcuda.so /usr/local/cuda/lib64/stubs/libcuda.so.1

export LD_LIBRARY_PATH=/usr/local/cuda/lib64/stubs:${LD_LIBRARY_PATH}

バージョン確認

> modinfo nvidia | grep version

 

2、GuppyのGPU版ダウンロード

2020 1/19現在、Guppyのv3.4.4が提供されている。log inしてsoftware downloadからlinux GPUビルドをダウンロードする。

https://community.nanoporetech.com/downloads

f:id:kazumaxneo:20200119132603p:plain

tar -xf ont-guppy_3.4.4_linux64.tar.gz
cd ont-guppy/bin/

./guppy_basecaller

$ ./guppy_basecaller

: Guppy Basecalling Software, (C) Oxford Nanopore Technologies, Limited. Version 3.4.4+a296acb

 

Usage:

 

With config file:"

  guppy_basecaller -i <input path> -s <save path> -c <config file> [options]

With flowcell and kit name:

  guppy_basecaller -i <input path> -s <save path> --flowcell <flowcell name>

    --kit <kit name>

List supported flowcells and kits:

  guppy_basecaller --print_workflows

 

Use GPU for basecalling:

  guppy_basecaller -i <input path> -s <save path> -c <config file>

    --device <cuda device name> [options]

 

Use server for basecalling:

  guppy_basecaller -i <input path> -s <save path> -c <config file>

    --port <server address> [options]

Command line parameters:

  --print_workflows                 Output available workflows.

  --flowcell arg                    Flowcell to find a configuration for

  --kit arg                         Kit to find a configuration for

  -m [ --model_file ] arg           Path to JSON model file.

  --chunk_size arg                  Stride intervals per chunk.

  --chunks_per_runner arg           Maximum chunks per runner.

  --chunks_per_caller arg           Soft limit on number of chunks in each 

                                    caller's queue. New reads will not be 

                                    queued while this is exceeded.

  --overlap arg                     Overlap between chunks (in stride 

                                    intervals).

  --gpu_runners_per_device arg      Number of runners per GPU device.

  --cpu_threads_per_caller arg      Number of CPU worker threads per 

                                    basecaller.

  --num_callers arg                 Number of parallel basecallers to create.

  --post_out                        Return full posterior matrix in output 

                                    fast5 file and/or called read message from 

                                    server.

  --stay_penalty arg                Scaling factor to apply to stay probability

                                    calculation during transducer decode.

  --qscore_offset arg               Qscore calibration offset.

  --qscore_scale arg                Qscore calibration scale factor.

  --temp_weight arg                 Temperature adjustment for weight matrix in

                                    softmax layer of RNN.

  --temp_bias arg                   Temperature adjustment for bias vector in 

                                    softmax layer of RNN.

  --hp_correct arg                  Whether to use homopolymer correction 

                                    during decoding.

  --builtin_scripts arg             Whether to use GPU kernels that were 

                                    included at compile-time.

  -x [ --device ] arg               Specify basecalling device: 'auto', or 

                                    'cuda:<device_id>'.

  -k [ --kernel_path ] arg          Path to GPU kernel files location (only 

                                    needed if builtin_scripts is false).

  -z [ --quiet ]                    Quiet mode. Nothing will be output to 

                                    STDOUT if this option is set.

  --trace_categories_logs arg       Enable trace logs - list of strings with 

                                    the desired names.

  --verbose_logs                    Enable verbose logs.

  --qscore_filtering                Enable filtering of reads into PASS/FAIL 

                                    folders based on min qscore.

  --min_qscore arg                  Minimum acceptable qscore for a read to be 

                                    filtered into the PASS folder

  --disable_pings                   Disable the transmission of telemetry 

                                    pings.

  --ping_url arg                    URL to send pings to

  --ping_segment_duration arg       Duration in minutes of each ping segment.

  --calib_detect                    Enable calibration strand detection and 

                                    filtering.

  --calib_reference arg             Reference FASTA file containing calibration

                                    strand.

  --calib_min_sequence_length arg   Minimum sequence length for reads to be 

                                    considered candidate calibration strands.

  --calib_max_sequence_length arg   Maximum sequence length for reads to be 

                                    considered candidate calibration strands.

  --calib_min_coverage arg          Minimum reference coverage to pass 

                                    calibration strand detection.

  --trim_threshold arg              Threshold above which data will be trimmed 

                                    (in standard deviations of current level 

                                    distribution).

  --trim_min_events arg             Adapter trimmer minimum stride intervals 

                                    after stall that must be seen.

  --max_search_len arg              Maximum number of samples to search through

                                    for the stall

  --reverse_sequence arg            Reverse the called sequence (for RNA 

                                    sequencing).

  --u_substitution arg              Substitute 'U' for 'T' in the called 

                                    sequence (for RNA sequencing).

  --override_scaling                Manually provide scaling parameters rather 

                                    than estimating them from each read.

  --scaling_med arg                 Median current value to use for manual 

                                    scaling.

  --scaling_mad arg                 Median absolute deviation to use for manual

                                    scaling.

  --trim_strategy arg               Trimming strategy to apply: 'dna' or 'rna' 

                                    (or 'none' to disable trimming)

  --dmean_win_size arg              Window size for coarse stall event 

                                    detection

  --dmean_threshold arg             Threhold for coarse stall event detection

  --jump_threshold arg              Threshold level for rna stall detection

  --pt_scaling                      Enable polyT/adapter max detection for read

                                    scaling.

  --pt_median_offset arg            Set polyT median offset for setting read 

                                    scaling median (default 2.5)

  --adapter_pt_range_scale arg      Set polyT/adapter range scale for setting 

                                    read scaling median absolute deviation 

                                    (default 5.2)

  --pt_required_adapter_drop arg    Set minimum required current drop from 

                                    adapter max to polyT detection. (default 

                                    30.0)

  --pt_minimum_read_start_index arg Set minimum index for read start sample 

                                    required to attempt polyT scaling. (default

                                    30)

  --log_speed_frequency arg         How often to print out basecalling speed.

  --barcoding_config_file arg       Configuration file to use for barcoding.

  --num_barcode_threads arg         Number of worker threads to use for 

                                    barcoding.

  --barcode_kits arg                Space separated list of barcoding kit(s) or

                                    expansion kit(s) to detect against. Must be

                                    in double quotes.

  --trim_barcodes                   Trim the barcodes from the output sequences

                                    in the FastQ files.

  --num_extra_bases_trim arg        How vigorous to be in trimming the barcode.

                                    Default is 0 i.e. the length of the 

                                    detected barcode. A positive integer means 

                                    extra bases will be trimmed, a negative 

                                    number is how many fewer bases (less 

                                    vigorous) will be trimmed.

  --arrangements_files arg          Files containing arrangements.

  --score_matrix_filename arg       File containing mismatch score matrix.

  --start_gap1 arg                  Gap penalty for aligning before the 

                                    reference.

  --end_gap1 arg                    Gap penalty for aligning after the 

                                    reference.

  --open_gap1 arg                   Penalty for opening a new gap in the 

                                    reference.

  --extend_gap1 arg                 Penalty for extending a gap in the 

                                    reference.

  --start_gap2 arg                  Gap penalty for aligning before the query.

  --end_gap2 arg                    Gap penalty for aligning after the query.

  --open_gap2 arg                   Penalty for opening a new gap in the query.

  --extend_gap2 arg                 Penalty for extending a gap in the query.

  --min_score arg                   Minimum score to consider a valid 

                                    alignment.

  --min_score_rear_override arg     Minimum score to consider a valid alignment

                                    for the rear barcode only (and min_score 

                                    will then be used for the front only when 

                                    this is set).

  --front_window_size arg           Window size for the beginning barcode.

  --rear_window_size arg            Window size for the ending barcode.

  --require_barcodes_both_ends      Reads will only be classified if there is a

                                    barcode above the min_score at both ends of

                                    the read.

  --allow_inferior_barcodes         Reads will still be classified even if both

                                    the barcodes at the front and rear (if 

                                    applicable) were not the best scoring 

                                    barcodes above the min_score.

  --detect_mid_strand_barcodes      Search for barcodes through the entire 

                                    length of the read.

  --min_score_mid_barcodes arg      Minimum score for a barcode to be detected 

                                    in the middle of a read.

  --num_barcoding_buffers arg       Number of GPU memory buffers to allocate to

                                    perform barcoding into. Controls level of 

                                    parallelism on GPU for barcoding.

  -q [ --records_per_fastq ] arg    Maximum number of records per fastq file, 0

                                    means use a single file (per worker, per 

                                    run id).

  --read_batch_size arg             Maximum batch size, in reads, for grouping 

                                    input files.

  --compress_fastq                  Compress fastq output files with gzip.

  -i [ --input_path ] arg           Path to input fast5 files.

  --input_file_list arg             Optional file containing list of input 

                                    fast5 files to process from the input_path.

  -s [ --save_path ] arg            Path to save fastq files.

  -l [ --read_id_list ] arg         File containing list of read ids to filter 

                                    to

  -p [ --port ] arg                 Hostname and port for connecting to 

                                    basecall service (ie 'myserver:5555'), or 

                                    port only (ie '5555'), in which case 

                                    localhost is assumed.

  -r [ --recursive ]                Search for input files recursively.

  --fast5_out                       Choice of whether to do fast5 output.

  --resume                          Resume a previous basecall run using the 

                                    same output folder.

  --max_block_size arg              Maximum size (in events) of outgoing 

                                    basecall messages (server mode only).

  --disable_events                  Disable the transmission of event tables 

                                    when receiving reads back from the basecall

                                    server.

  --progress_stats_frequency arg    Frequency in seconds in which to report 

                                    progress statistics, if supplied will 

                                    replace the default progress display.

  --client_id arg                   Optional unique identifier (non-negative 

                                    integer) for this instance of the Guppy 

                                    Client Basecaller, if supplied will form 

                                    part of the output filenames.

  -h [ --help ]                     produce help message

  -v [ --version ]                  print version number

  -c [ --config ] arg               Config file to use

  -d [ --data_path ] arg            Path to use for loading any data files the 

                                    application requires.

 

 

GPUを確認しておく。

nvidia-smi

$ nvidia-smi

Sun Jan 19 12:06:32 2020       

+-----------------------------------------------------------------------------+

| NVIDIA-SMI 390.129                Driver Version: 390.129                   |

|-------------------------------+----------------------+----------------------+

| GPU  Name        Persistence-M| Bus-Id        Disp.A | Volatile Uncorr. ECC |

| Fan  Temp  Perf  Pwr:Usage/Cap|         Memory-Usage | GPU-Util  Compute M. |

|===============================+======================+======================|

|   0  GeForce GTX 1080    Off  | 00000000:07:00.0  On |                  N/A |

| 19%   52C    P8     9W / 200W |    334MiB /  8118MiB |      0%      Default |

+-------------------------------+----------------------+----------------------+

                                                                               

+-----------------------------------------------------------------------------+

| Processes:                                                       GPU Memory |

|  GPU       PID   Type   Process name                             Usage      |

|=============================================================================|

|    0      1256      G   /usr/lib/xorg/Xorg                           162MiB |

|    0      1501      G   /usr/bin/gnome-shell                         100MiB |

|    0      3217      G   ...uest-channel-token=12657443395681243138    68MiB |

+-----------------------------------------------------------------------------+

GTX 1080の8 GB VRAM(GDDR5X)になっている。

 

 

実行方法

CPU版と同様の流れでランできる。異なるのはCPUスレッドの代わりにデバイス番号を指定するところ。1台しか利用できないなら--device autoか--device cuda:0を指定すればO.K。

guppy_basecaller \
--flowcell FLO-MIN106 \
--kit SQK-LSK109 \
-x cuda:0 \
-i fast5_dir \
-s output_dir2 -r
  • -x [ --device ]      Specify basecalling device: 'auto', or  'cuda:<device_id>'.
  • --flowcell      Flowcell to find a configuration for
  • -kit      Kit to find a configuration for
  • -i [ --input_path ]     Path to input fast5 files.
  • -s [ --save_path ]      Path to save fastq files.

100MB程度の小さなfast5データを使ってランタイムを調べた。

f:id:kazumaxneo:20200119144601p:plain

ラン中はGPU使用率がほぼ100%になる(右上)(nvtopを使用 *1)。

 

結果は

GTX 1080 =>18.5s

CPU(AMD 3700x) => 7m56.6s

25倍の差がついた。大きなデータでは、GPU版を使わないと終わらないのがよく分かりました。

 

参考

*1

nvtopのインストール

github

git clone https://github.com/Syllo/nvtop.git 
mkdir -p nvtop/build && cd nvtop/build
cmake .. -DNVML_RETRIEVE_HEADER_ONLINE=True
make
sudo make install

nvtopは複数GPUもモニターできます。上では、Terminator(参考にしたHP)を入れて端末を分割しています。

 


 

 

URMAP

2020 1/19 コマンドの誤り修正

2020 1/20 twitter追記

 

 次世代シーケンシングにより、ヒト機能ゲノミクス(Morozova and Marra、2008)から微生物メタゲノミクス(Gilbert and Dupont、2011)までの分野で劇的な進歩が可能になった。 次世代研究のデータ分析では、リードをヒトゲノム、ヒトエクソーム、または完全長の微生物ゲノムのコレクションなどの参照データベースにマッピングする必要がある。 マッピングは、シーケンスデータベース検索の特殊なケースであり、クエリシーケンスは短く、データベースシーケンスは長く、シーケンスの類似性は高くなる。 特定のクエリシーケンス(リード)の場合、マッピングの主な目的は、可能であれば最適な一致を報告することである。

 URMAPはk-merのハッシュテーブルインデックス、つまり、長さがkの固定長ワードを使用する。ここで、k = 24はヒトゲノムで推奨される。 インデックスは、メモリキャッシュミスを最小限に抑えるため、RAM内の特定のハッシュワード(スロット)に関連する情報を密接に保つように設計されている。リファレンス(ピン)で1回だけ見つかったスロットにフラグが付けられる。 与えられたクエリに対して、URMAPは最初に、リファレンス内で互いに近接している重複しないピンのペアを検索する(ブレース、図1を参照)。 中括弧が見つかった場合、位置合わせが試行され、検索が成功するとすぐに終了する。 それ以外の場合、シードと拡張の戦略(Altschul et al、1990)が続く。 

 

HP 

https://drive5.com/urmap/ 

URMAP quick start

https://drive5.com/urmap/manual/quickstart.html

 

 

実行方法

1、create index

urmap -make_ufi hg38.fa -output hg38.ufi

2、mapping

#paired-end
urmap -map2 sample_R1.fastq.gz -reverse sample_R2.fastq.gz \
-ufi hg38.ufi -samout sample.sam

#single
urmap -map sample.fastq.gz -ufi hg38.ufi -samout sample.sam

 

引用

URMAP, an ultra-fast read mapper

Robert C. Edgar
bioRxiv preprint first posted online Jan. 14, 2020

 

 

Cytobandファイルのダウンロード

 

Cytoband file format

https://software.broadinstitute.org/software/igv/Cytoband

 

Cytobandファイル(ギムザ染色された染色体のバンドのおおよその位置を表す)はUCSCから入手できる。

group => Mapping and Sequencing

f:id:kazumaxneo:20200117003958p:plain

track => Chromosome bandを選択。

f:id:kazumaxneo:20200117004111p:plain

get outputをクリックしてcytobandをダウンロードする。

f:id:kazumaxneo:20200117004228p:plain

 

 

 

circosVCF

 

 1000米ドル未満で全ゲノムをシーケンスするという目標が現実になり、医療における全ゲノムシーケンス利用が次第に増えている。 Genomics England (http://www.genomicsengland.co.uk/)などの大規模ゲノムプロジェクトは、数千のゲノムから収集された情報を使用して、健康と病気をよりよく理解することを目的としている。次世代のシーケンス研究から得られたデータを分析するために、さまざまなタイプのソフトウェアが開発されている。全ゲノムVCFファイルは、エクソームVCFファイルの平均10倍の大きさになる。ただし、既存のツールのほとんどは、全ゲノム研究によって生成されるファイルのサイズの増加を無視しているため、大容量メモリなどの余分な要件のために不適切になる。効率的な全ゲノムデータ解析ソフトウェアツールの必要性が高まっている。

 さまざまなプラットフォームが、ヒトNGSデータの分析用に提案されている。(一部略)VarSifter Teer et al (2012)は、グラフィカルユーザーインターフェイスGUI)の助けを借りて幅広い種類のクエリを定義できる、強力な公開ツールである。ただし、メモリでは使いやすさが制限される。 (一部略)

 ここでは、大きなVCFファイルの分析を実行するように設計されたソフトウェアツールであるVCF-Explorerを紹介する。 VCF-Explorerは軽量な構造で、メモリを使い果たすことなく非常に大きなファイルを処理できる。 GUIを使用すると、アノテーションのバリアントを簡単にフィルタリングできる。 ITのバックグラウンドを持たない臨床医と研究者は、簡単にツールを使用できる。初期ファイルのロードや前処理は行われないため、CPU時間とメモリが節約される。シンプルなパーソナルコンピューターや、より大きな機能を備えた大型サーバーで使用できる。

 

  

 

circosVCF demo

 

HP

https://www.ariel.ac.il/wp/fbl/software/

CircosVCF exercise

https://en-lifesci.m.tau.ac.il/sites/lifesci_en.tau.ac.il/files/media_server/all-units/CircosVCF_Ex.pdf

 

 

使い方

http://legolas.ariel.ac.il/~tools/CircosVCF/ にアクセスする。 

種を指定する。

f:id:kazumaxneo:20191121033400p:plain

 ゲノムの バージョンを指定する。

f:id:kazumaxneo:20191121033413p:plain

 描画する染色体を指定する。

f:id:kazumaxneo:20191121033415p:plain

saveをクリック。

 

VCFをアップロードする。

f:id:kazumaxneo:20200117010701p:plain

 

途中からうまく動作しなくなりました。ネットワークが不安定なのか、メンテナンスされていないのか不明ですが、できるようになったら追記します。

 

補足

ローカルのkaryotypeファイルを使う。chr番号とサイズを記載したタブ区切りファイルを用意する。

f:id:kazumaxneo:20200117005944p:plain

これをアップロードする。

 

 

 

引用
CircosVCF: circos visualization of whole-genome sequence variations stored in VCF files

Drori E, Levy D, Smirin-Yosef P, Rahimi O, Salmon-Divon M

Bioinformatics. 2017 May 1;33(9):1392-1393

 

関連