macでインフォマティクス

macでインフォマティクス

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

JBrowse 2 desktopのマッピングファイル表示機能を試す

 

JBrowse 2 desktopはシークエンシングデータのマッピングファイルの読み込みに対応しており、リファレンスゲノムに沿ったリードのアラインメントを表示することができる。

読み込み可能なデータ。ファイルの種類によってはインデックスも必要。

  • Tabixed VCF
  • Tabixed BED
  • Tabixed GFF
  • BAM
  • CRAM
  • BigWig
  • BigBed
  • .hic file (Juicebox)
  • PAF

(マニュアルより)

 

 

bamファイルを読み込んでみる。

File => Open Trackを選択

f:id:kazumaxneo:20220115122357p:plain

 

右端にAdd a trackパネルが表示される。

f:id:kazumaxneo:20220115123632p:plain

 

4つボタンがあり、左端のボタンはローカルのファイルの指定、左から2つ目のボタンはファイルのURLの指定となる。右2つのボタンはdropboxGoogle driveのファイル読み込みのボタンとなる(共有設定をONにしてURLを指定する)。ファイルの種類によっては、インデックスも提供する必要がある(画像2段目)。

f:id:kazumaxneo:20220115122329p:plain

NEXTをクリック。

 

続いて、トラックの名前、ファイルの属性を指定する。名前はsample1_chr21_mappingとした。

f:id:kazumaxneo:20220115123215p:plain

 

adapterTypeはファイルタイプによって変える。bamファイルならBamAdapterを指定。

f:id:kazumaxneo:20220115124423p:plain

 

trackTypeもデータのタイプによって変える。bamファイルならalignment trackを指定。

f:id:kazumaxneo:20220115124624p:plain

 

ADDボタンをクリックするとTrackがロードされる。

 

染色体を変えるには中央上のポジションのウィンドウをクリックする。

f:id:kazumaxneo:20220115125026p:plain

 

ロードされた(データはNA12878)。リファレンスゲノムの順鎖にアラインしたリードは赤色、逆鎖にアラインしたリードは青色で表示される。ミスマッチは色付きで描画される。

f:id:kazumaxneo:20220115125055p:plain

alignment trackはカバレッジとパイルアップ(リードの積み重ね)に分かれている。

 

ロングリードのマッピングファイルも表示できる(データはNA12878)。

f:id:kazumaxneo:20220115141855p:plain

 

リードをクリックすると、そのリードの配列やマッピングの詳細(tagやFLAGなど)が表示される。

f:id:kazumaxneo:20220115125716p:plain

 

リードの上で右クリックすると、リードの詳細を全てコピーできる。

f:id:kazumaxneo:20220115130119p:plain

 

リードの上やカバレッジトラックの上にマウスカーソルを合わせて停止すると、詳細が表示される。

f:id:kazumaxneo:20220115130303p:plain

 

Insertion

紫色は挿入、青色はソフトクリッピング、赤色はハードクリッピングをそれぞれ表す。

f:id:kazumaxneo:20220115135226p:plain

リードの30%をカバーしていると色のついた逆三角形でマークされる。10bp以上の挿入は、特別に大きな紫色の長方形で表示される。

 

Deletion

f:id:kazumaxneo:20220115135304p:plain

 

 

 

 

左上のマークをクリックしてshow center lineにチェックを付ける。

f:id:kazumaxneo:20220115131436p:plain

 

中央線が表示された。

f:id:kazumaxneo:20220115131626p:plain

 


トラックタイトル右端のマーク”⁝”(トリコロン、三点リーダ)をクリックすると、トラックをカスタマイズできる。

f:id:kazumaxneo:20220115132057p:plain


show soft clippingをオンにした。

f:id:kazumaxneo:20220115132401p:plain

 

ソフトクリップされたリードの領域も表示された。

f:id:kazumaxneo:20220115132441p:plain

 

カラースキームをPer-base qualityに変更。

f:id:kazumaxneo:20220115132912p:plain

カラースキームをStrandに変更。IGVと同じ配色になっている。

f:id:kazumaxneo:20220115132739p:plain

 

カラースキームをStrandに変更後、sortでRead strandを選択。

f:id:kazumaxneo:20220115142639p:plain

 

ソートのbase pairは特定のバリアントを持つリードがソートされるように積み上げられたリードを並べ替える。この機能は中心線がある部分のバリアントがソート対象になる。

f:id:kazumaxneo:20220115142950p:plain

 

また、アラインメント情報を示すタグでソートやフィルタリングしたり、色分けできるようになっている。

sort by tag

f:id:kazumaxneo:20220115144904p:plain

マニュアルでは、タグを利用して塩基の修飾に色を付ける例が説明されています。

 

現在画面に表示されている結果はSVG形式でダウンロードできる。アプリ左上のボタンをクリックしてExport SVGを選ぶ。

f:id:kazumaxneo:20220115143208p:plain

 

左上ボタンをクリックしてHorizontally flipを選ぶと染色体が反転する。反転した状態だと、座標を左から右ではなく右から左に移動させることができる。もう一度選ぶと元に戻る。

f:id:kazumaxneo:20220115144124p:plain

どう使うかだが、逆位などの複雑な構造多型を解釈しやすくするために便利かもしれない。

 

ロングリードを右クリックしてDotplot of read vs refを選択すると、リードとリファレンスのdot plotを表示できる。

f:id:kazumaxneo:20220115155543p:plain

linear dot plotビューが下のパネルで開いている。

 

同じく、ロングリードを右クリックしてLinear read vs refを選択すると、リードとリファレンスのシンテニーマップを表示できる。

f:id:kazumaxneo:20220115155730p:plain

これらの機能は構造多型が発生している領域の分析に有用だと思われる。

 

構造バリアントコール検査用に、SV inspector という機能も用意されている。

ADD => SV inspectorを選択

f:id:kazumaxneo:20220115160406p:plain

 

SVやCNVのファイルを読み込む。VCFまたはVCF.gz、BEDPE、STAR-fusionの結果ファイルを指定できる。

f:id:kazumaxneo:20220115160534p:plain


ここではcuteSVによるSVコール結果を選択した。

 

ロードされた。左にファイルの各行を表すテーブルビュー、右側にSVの環状全ゲノム表現での概要が表示される(circosビュー)。

f:id:kazumaxneo:20220115160706p:plain

テーブルは、検索とフィルタリングが可能。フィルタリングと検索の結果は環状ビューにも反映される。

 

環状ビューでフィーチャーをクリックするか、テーブルの左端の列にある三角形のドロップダウンをクリックすると、"スプリットビュー" または "ブレークポイントスプリットビュー" と呼ばれるデータの新しいビューが表示される。この機能によって、染色体間転座などの構造バリアントのブレークポイントを表示し、その部位の遺伝子アノテーションを読み込めば、SV発生部位にどのような遺伝子が存在しているか素早く確認することが可能になる。

f:id:kazumaxneo:20220115161254p:plain

マニュアルより転載

 

1000genome VCFのようなマルチサンプルVCFのポジションは1つで表示される(IGVとは異なる部分)。コール部位をクリックすると右下側にサンプルトラックが表示され、どのサンプルのコールか確認できる。

f:id:kazumaxneo:20220115234355p:plain

 

サンプルトラックではフィルターをかけて絞り込むことができる。プレーンテキスト検索または正規表現形式の検索を使用できる。 genotype フィールドに 1 と入力すると、0|1 または 1|1 のような最初の代替アレルを含むすべての genotypes が検索される(ツールの説明)。

f:id:kazumaxneo:20220115234805p:plain

 

 

マニュアルでは、bigwigの可視化、Hi-C データの可視化、全ゲノムスケールのCNVビューの表示方法などについても説明されています。アクセスしてみて下さい。

引用

JBrowse: a dynamic web platform for genome visualization and analysis
Robert Buels, Eric Yao, Colin M. Diesh, Richard D. Hayes, Monica Munoz-Torres, Gregg Helt, David M. Goodstein, Christine G. Elsik, Suzanna E. Lewis, Lincoln Stein & Ian H. Holmes 
Genome Biology volume 17, Article number: 66 (2016) 

 

関連


JBrowseゲノムブラウザのデスクトップアプリケーション JBrowse 2 desktop

 

HPより

JBrowseは、JavaScriptHTML5で作られた高速でフル機能のゲノムブラウザです。Webサイトやアプリに簡単に組み込むことができますが、スタンドアロンのWebページとして提供することも可能です。

 

version1との機能比較表。

f:id:kazumaxneo:20220114102533p:plain

 

HP

https://jbrowse.org/jb2/

twitter

https://twitter.com/usejbrowse

Document

https://jbrowse.org/jb2/docs/

 

 

 

インストール

Desktop版のDownalod

JBrowse

f:id:kazumaxneo:20220114100935p:plain

指示に従ってインストールする。

 

windows版を立ち上げた。

f:id:kazumaxneo:20220114101132p:plain

 

ヒトとマウスのゲノムは最初から選択できるようになっている。他のゲノムを登録するならOPEN SEQUENCE FILEを選択する。

f:id:kazumaxneo:20220114101211p:plain

 

OPEN SEQUENCE FILEを選択すると、下の写真のようなゲノムを登録するウィンドウに切り替わる。まず名前を決める。次に、下のボタンからゲノムアセンブリfasta形式ファイルとindexのfasta.faiファイルを読み込ませる。URLを指定してロードすることもできる。

f:id:kazumaxneo:20220114101803p:plain

submitすると読み込まれる。

 

ここからはヒトのhg38アセンブリで見ていく。

ゲノムを選択すると、JBrowse 2 desktopアプリ内に新しいウィンドウが現れる。

f:id:kazumaxneo:20220114102155p:plain

 

左上のボタンをクリックしてビューモードを選択する。まずlinearGenomeviewを選択する。

f:id:kazumaxneo:20220114102210p:plain

 

するとアプリ内にウィンドウが出現して、線状モードでゲノムを閲覧できるようになる。このように、JBrowse 2はアプリ内で他のビューと一緒に表示できる「ビュータイプ」に対応している。サードパーティプラグインを追加すると、ビューを追加できるようになっている。

f:id:kazumaxneo:20220114102315p:plain

 

染色体番号を選択して読み込む。

f:id:kazumaxneo:20220114102434p:plain

 

赤いboxが現在表示されている領域。左右の矢印ボタンで左右にスクロールできる他、ドラッグして囲むことで素早く移動できる。画面をドラッグして引っ張っても移動できる。右のパンボタンでは拡大縮小できる。

f:id:kazumaxneo:20220114103232p:plain

動作はとてもスムーズで、他のゲノムビューアより洗練されていると感じる。

 

ポジションや染色体名をタイプするとその染色体:ポジションにジャンプできる。

f:id:kazumaxneo:20220114230944p:plain

 

中央のOPEN TRACK SELECTORをクリックすると、右側にアノテーショントラックウィンドウが表示される。

f:id:kazumaxneo:20220114104049p:plain

 

クリックすることでそのアノテーションのトラックをゲノムに沿って表示することができる。

f:id:kazumaxneo:20220114201427p:plain

 

補足;アノテーショントラックウィンドウが隠れた場合、左端のマークをクリックすると再表示される。

f:id:kazumaxneo:20220114234103p:plain

 

Clinvar CNVとClivar variants、Gencode v36を読み込んだ。

f:id:kazumaxneo:20220114104123p:plain

 

情報がロードされていない場合はトラック上のForce LOADをクリック。

f:id:kazumaxneo:20220114104650p:plain

 

アノテーショントラックの境界はドラッグして変更できる。縦に長くした。

f:id:kazumaxneo:20220114104855p:plain

 

トラックのexonなどをクリックすると右に情報が表示される。

f:id:kazumaxneo:20220114104947p:plain

 

Refseq sequenceをロードして最大近くまで拡大するとリファレンス配列が表示される。

f:id:kazumaxneo:20220114234541p:plain

 

用意されているアノテーションについては、"..."をクリックする事で詳細を確認できる。

f:id:kazumaxneo:20220114104226p:plain

 

詳細の例 (Clinvar CNV)

f:id:kazumaxneo:20220114104253p:plain

 

File->Openで、ローカルマシンのアノテーションやURL先のファイルを読み込むことができる。

以下のファイル形式に対応している(マニュアルより)。

  • Tabixed VCF
  • Tabixed BED
  • Tabixed GFF
  • BAM
  • CRAM
  • BigWig
  • BigBed
  • .hic file (Juicebox)
  • PAF

 

 

複数のビューモードを同時に表示することもできる。

上のメニューのADDをクリック。

f:id:kazumaxneo:20220114232732p:plain

 

Circular viewを追加した。

f:id:kazumaxneo:20220114233207p:plain

 

TOOLS => plugin storeからはプラグインを追加できる。

f:id:kazumaxneo:20220114233304p:plain

 

利用可能なプラグインが表示される。

f:id:kazumaxneo:20220114233416p:plain

 

multiple sequence alignment browser: MsaViewを追加した。追加したプラグインもADDから読み込める。

f:id:kazumaxneo:20220114233531p:plain

 

Ideogramプラグインを追加。

f:id:kazumaxneo:20220114233650p:plain

 

topからは最近開いたセッションを再ロードすることができる。

f:id:kazumaxneo:20220114232524p:plain

 

明日はマッピングファイルの読み込み、閲覧などを紹介します。

1/15


引用

JBrowse: a dynamic web platform for genome visualization and analysis
Robert Buels, Eric Yao, Colin M. Diesh, Richard D. Hayes, Monica Munoz-Torres, Gregg Helt, David M. Goodstein, Christine G. Elsik, Suzanna E. Lewis, Lincoln Stein & Ian H. Holmes 
Genome Biology volume 17, Article number: 66 (2016) 

 

 

 

 

超高速で高精度なアンプリコンシークエンス解析ツール LotuS2

 

 アンプリコンシークエンスは、マイクロバイオームのプロファイリングにおいて確立されたコスト効率の高い手法である。しかし、このデータを処理するための多くのツールは、大きなデータセットを処理するためにバイオインフォマティクスのスキルと高い計算能力の両方を必要とする。さらに、ロングリードのアンプリコンデータを解析できるツールは数少ない。このギャップを埋めるために、LotuS2 (Less OTU Scripts 2)パイプラインを開発し、ユーザーフレンドリーでリソースに優しい、生のアンプリコン配列の多用途な分析を可能にした。

 LotuS2 では、6 種類の配列クラスタリングアルゴリズムと豊富な前処理・後処理オプションにより、パラメータを完全に調整できる専門家から、異なるシナリオのためのデフォルトが提供されている初心者まで、柔軟にデータ解析を行うことができる。3つの独立した腸と土壌のデータセットベンチマークしたところ、LotuS2は他のパイプラインと比較して平均29倍速く、しかも技術的な複製サンプルのアルファおよびベータ多様性をより良く再現することができた。さらに、既知の分類群構成を持つ模擬コミュニティをベンチマークした結果、LotuS2は、他のパイプラインと比較して、より高い割合で正しく属と種を復元した(それぞれ、98%と57%)ことが示された。ASV/OTUレベルでは、LotuS2の精度およびFスコアが最も高く、16S配列が正しく再構成された割合も高かった。

 LotuS2は、軽量でユーザーフレンドリーなパイプラインであり、高速、高精度、合理的である。高いデータ利用率と信頼性により、数分でハイスループットなマイクロバイオーム解析が可能になる。LotuS2は、GitHub、conda、またはGalaxyウェブインタフェースから入手可能で、http://lotus2.earlham.ac.uk/ にドキュメントが掲載されている。

 

Documentation

http://lotus2.earlham.ac.uk/

 

Galaxy版も用意されている。

https://usegalaxy.eu

f:id:kazumaxneo:20220110023827p:plain

 

 

インストール

依存

  • USEARCH ver 7, 8, 9 or 10 (link)

Github

#conda (bioconda)
mamba create -n lotus2 -y
conda activate lotus2
mamba install -c bioconda lotus2 -y
#ダウンロードしたUSEARCH実行形式イナリのパスを指定する
lotus2 -link_usearch <path>/<to>/usearch

#もしくは対話式で進めるインストールスクリプトを使う。どのデータベースを使うかも聞かれる。
git clone https://github.com/hildebra/lotus2.git
cd lotus2/
perl autoInstall.pl

> lotus2 

Lotus2 usage:

lotus2 -i <input fasta|fastq|dir> -o <output_dir> -m/-map <mapping_file>

 

Minimal example (from lotus2 dir):

./lotus2 -i Example/ -m Example/miSeqMap.sm.txt -o myTestRun -CL vsearch

 

#### OPTIONS ####

 

Basic Options:

 

  -i <file|dir>           In case that fastqFile or fnaFile and qualFile were 

                          specified in the mapping file, this has to be the 

                          directory with input files 

  -m|-map <file>          Mapping file 

  -o <dir>                Warning: The output directory will be completely 

                          removed at the beginning of the LotuS2 run. Please 

                          ensure this is a new directory or contains only 

                          disposable data! 

 

 

Workflow Options:

 

  -forwardPrimer <string>

                          give the forward primer used to amplify DNA region 

                          (e.g. 16S primer fwd) 

  -keepOfftargets <0|1>   (0)?!?: keep offtarget hits against offtargetDB in 

                          output fasta and otu matrix. (Default 0) 

  -keepTmpFiles <0|1>     (1) save extra tmp files like chimeric OTUs or the 

                          raw blast output in extra dir. (0) do not save these. 

                          (Default: 0) 

  -keepUnclassified <0|1>

                          (1) Includes unclassified OTUs (i.e. no match in 

                          RDP/Blast database) in OTU and taxa abundance matrix 

                          calculations. (0) does not take these OTUs into 

                          account. (Default: 1) 

  -offtargetDB <file>     Remove likely contaminant OTUs/ASVs based on 

                          alignment to provided fasta. This option is useful for 

                          low-bacterial biomass samples, to remove possible host 

                          genome contaminations (e.g. human/mouse genome) 

  -redoTaxOnly <0|1>      (1) Only redo the taxonomic assignments (useful for 

                          replacing a DB used on a finished lotus run). (0) 

                          Normal lotus run. (Default: 0) 

  -reversePrimer <string>

                          give the reverse primer used to amplify DNA region 

                          (e.g. 16S primer rev) 

  -saveDemultiplex <0|1|2|3>

                          (1) Saves all demultiplexed reads (unfiltered) in the 

                          [outputdir]/demultiplexed folder, for easier data 

                          upload. (2) Only saves quality filtered demultiplexed 

                          reads and continues LotuS2 run subsequently. (3) Saves 

                          demultiplexed, filtered reads into a single fq, with 

                          sample ID in fastq/a header. (0) No demultiplexed 

                          reads are saved. (Default: 0) 

  -taxOnly <file>         Skip most of the lotus pipeline and only run a 

                          taxonomic classification on a fasta file. E.g. lotus2 

                          -taxOnly <some16S.fna> -refDB SLV 

  -tolerateCorruptFq <0|1>

                          (1) Continue reading fastq files, even if single 

                          entries are incomplete (e.g. half of qual values 

                          missing). (0) Abort lotus run, if fastq file is 

                          corrupt. (Default: 0) 

  -useVsearch <0|1>       (0) Use usearch for internal tasks such as remapping 

                          reads on OTUs, chimera checks. (1) will use vsearch 

                          for these tasks. This option is independent of the -CL 

                          UPARSE/UNOISE option, and -taxAligner assignment 

                          usearch/vsearch options. (Default: 0) 

  -verbosity <0-3>        Level of verbosity from printing all program calls 

                          and program output (3) to not even printing errors 

                          (0). (Default: 1) 

 

 

Clustering Options:

 

  -CL|-clustering <uparse|swarm|cdhit|unoise|dada2>

                          Sequence clustering algorithm: (1) UPARSE, (2) swarm, 

                          (3) cd-hit, (6) unoise3, (7) dada2. Short keyword or 

                          number can be used to indicate clustering (Default: 

                          uparse) 

  -chim_skew <num>        Skew in chimeric fragment abundance (uchime option). 

                          (Default: 2) 

  -count_chimeras<T/F>    Add chimeras to count up OTUs/ASVs. (Default: F) 

  -deactivateChimeraCheck <0|1|2|3>

                          (0) do OTU chimera checks. (1) no chimera check at 

                          all. (2) Deactivate deNovo chimera check. (3) 

                          Deactivate ref based chimera check. (Default: 0) 

  -derepMin<num>          Minimum size of dereplicated clustered, one form of 

                          noise removal. Can also have a more complex syntax, 

                          see examples. (Default: 1) 

  -endRem <string>        DNA sequence, usually reverse primer or reverse 

                          adaptor; all sequence beyond this point will be 

                          removed from OTUs. This is redundant with the 

                          "ReversePrimer" option from the mapping file, but 

                          gives more control (e.g. there is a problem with 

                          adaptors in the OTU output). (Default: "") 

  -id <0-1>               Clustering threshold for OTUs. (Default: 0.97) 

  -readOverlap <num>      The maximum number of basepairs that two reads are 

                          overlapping. (Default: 300) 

  -swarm_distance <1,2,3,..> 

                          Clustering distance for OTUs when using swarm 

                          clustering. (Default: 1) 

  -xtalk <0|1>            (1) check for crosstalk. Note that this requires in 

                          most cases 64bit usearch. (Default: 0) 

 

 

Taxonomy Options:

 

  -ITSx <0|1>             (1) use ITSx to only retain OTUs fitting to ITS1/ITS2 

                          hmm models; (0) deactivate. (Default: 1) 

  -LCA_cover <0-1>        Min horizontal coverage of an OTU sequence against 

                          ref DB. (Default: 0.5) 

  -LCA_frac <0-1>         Min fraction of database hits at taxlevel, with 

                          identical taxonomy. (Default: 0.9) 

  -amplicon_type <SSU|LSU|ITS|ITS1|ITS2>

                          (SSU) small subunit (16S/18S), (LSU) large subunit 

                          (23S/28S) or internal transcribed spacer 

                          (ITS|ITS1|ITS2). (Default: SSU) 

  -buildPhylo <0,1,2,>    (0) do not build OTU phylogeny; (1) use fasttree2; 

                          (2) use IQ-TREE 2. (Default: 1) 

  -greengenesSpecies <0|1>

                          (1) Create greengenes output labels instead of OTU 

                          (to be used with greengenes specific programs such as 

                          BugBase). (Default: 0) 

  -itsx_partial <0-100>   Parameters for ITSx to extract partial (%) ITS 

                          regions as well. (0) deactivate. (Default: 0) 

  -lulu <0|1>             (1) use LULU (https://github.com/tobiasgf/lulu) to 

                          merge OTUs based on their occurrence. (Default: 1) 

  -rdp_thr <0-1>          Confidence thresshold for RDP.(Default: 0.8) 

  -refDB <SLV|GG|HITdb|PR2|UNITE|beetax>

                          (SLV) Silva LSU (23/28S) or SSU (16/18S), (GG 

                          greengenes (only SSU available), (HITdb) (SSU, human 

                          gut specific), (PR2) LSU spezialized on Ocean 

                          environmentas, (UNITE) ITS fungi specific, (beetax) 

                          bee gut specific database and tax names. \nDecide 

                          which reference DB will be used for a similarity based 

                          taxonomy annotation. Databases can be combined, with 

                          the first having the highest prioirty. E.g. "PR2,SLV" 

                          would first use PR2 to assign OTUs and all unaasigned 

                          OTUs would be searched for with SILVA, given that 

                          \"-amplicon_type LSU\" was set. Can also be a custom 

                          fasta formatted database: in this case provide the 

                          path to the fasta file as well as the path to the 

                          taxonomy for the sequences using -tax4refDB. See also 

                          online help on how to create a custom DB. (Default: 

                          SLV) 

  -tax4refDB <file>       In conjunction with a custom fasta file provided to 

                          argument -refDB, this file contains for each fasta 

                          entry in the reference DB a taxonomic annotation 

                          string, with the same number of taxonomic levels for 

                          each, tab separated. 

  -taxAligner <0|blast|lambda|utax|vsearch|usearch>

                          Previously doBlast. (0) deactivated (just use RDP); 

                          (1) or (blast) use Blast; (2) or (lambda) use LAMBDA 

                          to search against a 16S reference database for 

                          taxonomic profiling of OTUs; (3) or (utax): use UTAX 

                          with custom databases; (4) or (vsearch) use VSEARCH to 

                          align OTUs to custom databases; (5) or (usearch) use 

                          USEARCH to align OTUs to custom databases. (Default: 

                          0) 

  -taxExcludeGrep <string>

                          Exclude taxonomic group, these OTUs will be assigned 

                          as unknown instead. E.g. -taxExcludeGrep 

                          Chloroplast|Mitochondria (Default: ) 

  -tax_group <bacteria|fungi>

                          (bacteria) bacterial 16S rDNA annnotation, (fungi) 

                          fungal 18S/23S/ITS annotation. (Default: bacteria) 

  -useBestBlastHitOnly <0|1>

                          (1) do not use LCA (lowest common ancestor) to 

                          determine most likely taxonomic level (not 

                          recommended), instead just use the best blast hit. (0) 

                          LCA algorithm. (Default: 0) 

  -utax_thr <0-1>         Confidence thresshold for UTAX. (Default: 0.8) 

 

 

Further Options:

 

  -barcode|-MID <file>    Filepath to fastq formated file with barcodes (this 

                          is a processed mi/hiSeq format). The complementary 

                          option in a mapping file would be the column 

                          "MIDfqFile". (Default: "") 

  -c <file>               LotuS.cfg, config file with program paths. (Default: 

                          <LotuS2_dir>/lOTUs.cfg) 

  -p <454/miSeq/hiSeq/PacBio>

                          sequencing platform: PacBio, 454, miSeq or hiSeq. 

                          (Default: miSeq) 

  -q <file>               .qual file associated to fasta file. This is an old 

                          format that was replaced by fastq format and is rarely 

                          used nowadays. (Default: "") 

  -s|sdmopt <file>        SDM option file, defaults to "configs/sdm_miSeq.txt" 

                          in current dir. (Default: miSeq) 

  -tmp|-tmpDir <dir>      temporary directory used to save intermediate 

                          results. (Default: <outputDir>/tmpDir) 

  -t|-threads <num>       number of threads to be used. (Default: 1) 

 

 

Other uses of pipeline (quits after execution):

 

  -check_map <file>       Mapping_file: only checks mapping file and exists. 

  -create_map <file>      mapping_file: creates a new mapping file at location, 

                          based on already demultiplexed input (-i) dir. E.g. 

                          lotus2 -create_map mymap.txt -i 

                          /home/dir_with_demultiplex_fastq 

  -link_usearch <file>    Provide the absolute path to your local usearch 

                          binary file, this will be installed to be useable with 

                          LotuS2 in the future. 

  -v                      Print LotuS2 version 

 

 

テストラン

最小設定のラン。

git clone https://github.com/hildebra/lotus2.git
cd lotus2/
lotus2 -i Example/ -m Example/miSeqMap.sm.txt -o myTestRun

biocondaで配布されているversion2.16では途中でSegmentation faultが起きる。ランできるようになったら追記します。

 

galaxyでは問題なくラン出来た。

出力例

f:id:kazumaxneo:20220113002409p:plain

引用

LotuS2: An ultrafast and highly accurate tool for amplicon sequencing analysis
Ezgi Özkurt, Joachim Fritscher, Nicola Soranzo, Duncan Y. K. Ng, Robert P. Davey, Mohammad Bahram, Falk Hildebrand

bioRxiv, Posted December 24, 2021

 

関連


 

VCFの要約統計プロットを出力する vcfstats

 

開発の動機(マニュアルより)

bcftoolsやjvarkitなどVCFファイルの統計情報をプロットするツールはいくつか存在する。しかし、いずれも 特定の指標をプロットする、プロットをカスタマイズする、特定のフィルタでバリアントにフォーカスする、を行うことができない。RのパッケージvcfRは、上記の一部を行うことができるが、しかし、VCF全体をメモリに読み込む必要があり、大きなVCFファイルには不親切である。

 

Documentation

https://pwwang.github.io/vcfstats/

 

インストール

Github

#ここではcondaの仮想環境に導入する
mamba create -n vcfstats python=3.9 -y
conda activate vcfstats
pip install -U vcfstats

#docker
docker run --rm justold/vcfstats:latest vcfstats
# or
singularity run docker://justold/vcfstats:latest vcfstats

>

$ vcfstats -h

 DESCRIPTION:

  vcfstats v0.2.0: Powerful VCF statistics.                                     

 USAGE:

  vcfstats --vcf PATH --outdir AUTO --formula LIST --title LIST [OPTIONS]

 REQUIRED OPTIONS:

  -v, --vcf <PATH>          - The VCF file.                                     

  -o, --outdir <AUTO>       - The output directory.                             

  -f, --formula <LIST>      - The formulas for plotting in format of Y ~ X,     

                              where Y and X should be either an entry or an     

                              aggregation.                                      

  --title <LIST>            - The title of each figure, will be used to name the

                              output files as well.                             

 OPTIONAL OPTIONS:

  --loglevel <STR>          - The logging level. Default: info                  

  --figtype <LIST>          - Your preferences for type of plot for each        

                              formula. Default: \                             

  --figfmt <LIST>           - Your preferences for format of figure for each    

                              formula, Any file format supported by matplotlib. 

                              Default is png. Default: \                      

  -r, --region <LIST>       - Regions in format of CHR or CHR:START-END         

                              Default: \                                      

  -R, --Region <AUTO>       - Regions in a BED file. If both --region/--Region  

                              are provided, regions will be merged together.    

                              Default: None                                     

  -p, --passed [BOOL]       - Only analyze variants that pass all filters.      

                              This does not work if FILTER entry is in the      

                              analysis.                                         

                              Default: False                                    

  -l, --list [BOOL]         - List all available macros. Default: False         

  -s, --savedata [BOOL]     - Whether save the plotting data for further        

                              exploration. Default: False                       

  --macro <PATH>            - A user-defined macro file. Default: None          

  --ggs <LIST>              - Extra ggplot2 expressions for each plot           

                              Default: \                                      

  --devpars <NS>            - The device parameters for plots. To specify       

                              devpars for each plot, use a configuration file.  

  -c, --config <AUTO>       - A configuration file defining how to plot in TOML 

                              format.                                           

                              If this is provided, CLI arguments will be        

                              overwritten if defined in this file. Default: None

  -h, --help                - Print help information for this command           

 OPTIONAL OPTIONS UNDER --devpars:

  --devpars.width <INT>     - The width of the plot Default: 2000               

  --devpars.height <INT>    - The height of the plot Default: 2000              

  --devpars.res <INT>       - The resolution of the plot Default: 300         

 

 

テストラン

vcfファイルとconfigファイル(詳細)を指定する。出力はexamples/、図のタイトルは'Number of variants on each chromosome'、Y軸 と X軸の軸ラベル はCOUNT<INT>とCONTIG。出力ディレクトリは前もって作成しておく必要がある。

git clone https://github.com/pwwang/vcfstats.git
cd vcfstats/

mkdir out
vcfstats --vcf examples/sample.vcf \
  --outdir out/ \
    --formula 'COUNT(1) ~ CONTIG' \
    --title 'Number of variants on each chromosome' \
    --config examples/config.toml
  • --vcf   The VCF file
  • --outdir   The output directory. 
  • -c    A configuration file defining how to plot in TOML format. If this is provided, CLI arguments will be overwritten if defined in this file. Default: None
  • --title    The title of each figure, will be used to name the output files as well.
  • --formula    The formulas for plotting in format of Y ~ X, where Y and X should be either an entry or an aggregation.

出力

out/number-of-variants-on-each-chromosome.col.png

f:id:kazumaxneo:20220112004532p:plain

 

Githubでは--formulaオプションを調節して特定の染色体だけプロットするようにカスタムする例や他の要約統計プロットを出力するコマンドについて説明されています。アクセスしてみて下さい。

引用

GitHub - pwwang/vcfstats: Powerful statistics for VCF files

 

関連


 

ネットワークに基づく遺伝子セットエンリッチメント解析を行う NGSEA

 

 遺伝子発現表現型の遺伝子セット解析には、 over-representationアプローチとaggregate scoreアプローチという2つの主要なアプローチがある(Irizarry et al.、2009)。 over-representationアプローチでは、発現データセットから差分発現遺伝子(DEG)群を選択し、選択したDEGの中で各注釈付き遺伝子セットの過剰発現の有意性を超幾何検定などの統計検定により計算する。この方法は合理的であるが、いくつかの欠点がある。例えば、この方法では、有意性の低い遺伝子は発現表現型において重要でない遺伝子として扱われる。したがって、結果はDEGsを選択するために使用されるカットオフに大きく依存することになる。また、有意な遺伝子間の相対的な順序情報は考慮されていない。

Over-representationアプローチの解析限界は、メンバー遺伝子の全遺伝子特異的スコアに基づいて各注釈付き遺伝子セットにスコアを割り当てる、集約的スコアアプローチによって克服できる。Gene set enrichment analysis (GSEA) (Subramanian et al., 2005) は、現在利用可能な最も一般的な集約的スコアアプローチである。GSEAでは、まず発現プロファイルの遺伝子を発現差に基づく遺伝子特異的スコアで順位付けし、次に修正Kolmogorov-Smirnov (K-S) 検定に基づいて各注釈付き遺伝子集合のエンリッチメントスコアを計算する。しかし、GSEAはその人気の高さとは裏腹に、いくつかの欠点も持っている。例えば、GSEAは一方向に制御された遺伝子群、すなわち、アップレギュレートまたはダウンレギュレートのいずれかを識別するために設計されている。もし、ある遺伝子セットが、up-regulationとdown-regulationが等しく分布しているDEGsにマッチした遺伝子を持っている場合、その発現表現型との関連はGSEAでは検出されない可能性がある。この限界を克服するために、up-regulationとdown-regulationの両方の遺伝子スコアの絶対値を計算する絶対濃縮度(AE)と呼ばれる修正GSEAが開発された(Saxena et al.)。

 GSEAのもう一つの欠点は、DEGが必ずしも遺伝子セットで表される分子プロセスの責任者である機能遺伝子を表しているとは限らないことである。そのため、DEGは必ずしもその分子プロセスを司る機能遺伝子を表しているわけではなく、その分子プロセスにおける真の機能遺伝子に邪魔された制御異常遺伝子である可能性がある。GSEAは有意なDEGのスコアに基づいて各遺伝子集合にスコアを割り当てるため、有意な発現変化を示さない真の機能遺伝子からなる遺伝子集合は、この方法では捕捉されないことになる。この解析上の限界は、機能遺伝子ではなく、発現シグネチャーに基づく注釈付き遺伝子セットを用いることで部分的に克服できるかもしれない。例えば、MSigDBはGSEAで使用するために設計され、遺伝子発現データから得られた多くのシグネチャー遺伝子セットを含んでいる(Liberzon et al., 2011)。しかし、生物学的プロセスや疾患に関するアノテーション遺伝子のデータベースの大半は、疾患の原因遺伝子などの機能性遺伝子をベースとしている。

 ネットワークに基づく遺伝子発現差解析は、疾患原因遺伝子の優先順位付け(Nitsch et al., 2009)やがん細胞株の必須遺伝子(Jiang et al., 2015)に利用されている。これらの手法は、腫瘍形成などの疾患過程の機能遺伝子は、機能ネットワークにおいてその病態のDEGに囲まれる傾向があるという考えに基づいている。そこで本著者らは、局所的なサブネットワーク(すなわち、各遺伝子とその近傍を結ぶネットワーク)の発現差によって遺伝子を順序付けることで、関連する生物学的プロセスに関連する機能的遺伝子セットを捉える能力が向上すると仮定している。本研究では、個々の遺伝子だけでなく、機能ネットワークにおける隣接遺伝子の発現差を利用して、機能性遺伝子セットの濃縮スコアを測定するネットワークベースGSEA(NGSEA)を発表する。ネットワークベースの遺伝子セット解析手法は既にいくつか提案されているが、これらの手法は、予め選択された2つの遺伝子セット、データベースからのアノテーション遺伝子セット、実験からのクエリー遺伝子セット間の関連性を分子ネットワーク内での相対的近接性に基づいて同定するover-representation approachを改良した手法である(Alexeyenko et al.2012; Glaab et al.2012; McCormack et al.2013; Wang et al.2012).知る限り、NGSEAは集約的なスコアアプローチを適用した最初のネットワークベースの遺伝子セット解析手法である。

 NGSEAがマッチした遺伝子発現データセットに対するKEGGパスウェイ遺伝子セット(Kanehisa et al.、2017)の検索においてGSEAを上回ることを発見した。また、NGSEAをいくつかの疾患の薬剤優先順位付けに適用し、マッチした癌関連遺伝子発現データセットに対する既知の薬剤の検索能力において、NGSEAがConnectivity Map (CMap) (Lamb et al., 2006) より大幅に優れたパフォーマンスを示すことを見出した。NGSEA を用いて FDA が承認した薬剤が大腸がんに対して抗がん作用を持つかどうかを解析し、現在抗炎症薬として使用されている化学物質ブデソニドの抗がん作用を実験的に検証した。NGSEAは、Webベースのソフトウェアとして自由に利用することができる(www.inetbio.org/ngsea)。

 

Tutorial

http://www.inetbio.org/ngsea/tutorial.php

 

流れだけ簡単に確認します。

webサービス

NGSEA - Network-augmented Gene Set Enrichment Analysisにアクセスする。

f:id:kazumaxneo:20220110231606p:plain

ここではexampleデータを使う。

f:id:kazumaxneo:20220110232342p:plain

発現行列はテキスト形式(*.rnk);で遺伝子識別子<tab>値、を記載する。遺伝子識別子はENTREZ遺伝子IDまたはNCBI公式遺伝子記号のいずれかを指定する。clsファイルは表現型ラベルファイルと呼ばれるファイル。以前説明した(リンク)。exampleデータがダウンロードできるようになっている。

 

出力例

f:id:kazumaxneo:20220110232623p:plain

 

Prelanked(すでにP値など計算済みのもの)

f:id:kazumaxneo:20220110232704p:plain

 

出力例

f:id:kazumaxneo:20220110232808p:plain

結果は.tsvファイルとしてダウンロードできる。

 

引用

NGSEA: Network-Based Gene Set Enrichment Analysis for Interpreting Gene Expression Phenotypes with Functional Gene Sets
Heonjong Han, Sangyoung Lee, Insuk Lee

Mol Cells. 2019 Aug; 42(8): 579–588

複数条件下での時間経過トランスクリプトームデータを解析するためのウェブサービス TimesVector-Web

 

 遺伝子発現データの時間経過解析は、ある生物学的メカニズムの時間経過に伴う遺伝子発現の変調パターンを明らかにするのに有利である。例えば、正常者と癌患者のコホートなど、2つの条件間で有意に差のある発現遺伝子(DEG)を検索することは一般的に行われている。このような解析は、現在の遺伝子発現状態の1つのスナップショット内のトランスクリプトームの違いをペアワイズで観察するように調整されている。そのため、重要な遺伝子制御機能の動態や進行中の変遷が見落とされることがある。
時間経過解析データは、遺伝子発現の有意な変調や他のオミックスデータを観察し、実験の条件に関して関心のある生物学的現象を説明するために作成される。時間経過解析は複雑な手順となり、研究者の解析の視点によって異なる結果をもたらす可能性がある。例えば、従来のDEG解析を採用する場合、時間点の組合せを比較することが選択される。このようなアプローチでは、各DEGペアの結果だけでは解釈できないため、結果の事後分析が必要となる。生物学的に意味のある結果を導き出すには、それらを後から統合する必要がある。実験設定としては、単一の条件下で時間経過データを作成した場合、同定された遺伝子発現パターンは、関連する条件に対する反応を説明できることが期待される。一方、複数の条件下で時間経過データを作成した場合、条件間で有意に異なる(または類似する)遺伝子発現パターンを見つけることが解析の主な関心事となる。
 時間経過データを総称すると、時間経過とともに有意に異なる発現パターンをもたらす遺伝子を探索することができる。ここで、遺伝子の発現パターンを利用することで、より詳細に生物学的なメカニズムを理解することができる。例えば、あるストレスに対して、ある経路の遺伝子がどのように応答するかを経時的に観察することができる。さらに重要なことは、そのような遺伝子がストレスに応答し始めた時期を特定することができることである。しかし、時間経過データは高価であり、サンプリングする時点を選択することは、専門家の知識と推定または既知の条件応答性遺伝子に関する注意深い事前実験が必要な、自明ではない仕事である。このような慎重な設計を行わないと、選択した時間経過の中で重要な遺伝子発現調節が捉えられない可能性がある。これまで、発現量の異なる遺伝子を同定するための多くの時間経過解析法が開発されてきた。しかし、結果の解釈やツールの使い方が難しいため、多くの研究がDEG法に基づいて解析を行っている。生物学者にとって、オフラインのソフトウェアは、時に特定のシステムと互換性がなかったり、あるいは古くなってしまったりして、使いにくいものであるかもしれない。それでも、ウェブベースの時間経過遺伝子発現解析プラットフォームを提供するツールはわずかしかない。GEstureは、時系列遺伝子発現データからユーザが指定した遺伝子発現パターンを検索するためのウェブベースのツールである。したがって、従来の方法でクラスタリングを行うのではなく、ユーザが指定した検索パターンに依存した結果を得ることができる。そのため、指定したパターンと異なるパターンを見逃す可能性がある。
 このような難解な時間経過解析の負担を軽減するために、複数条件の時間経過遺伝子発現データを解析するための簡便なウェブサービスを開発した。このウェブサービスは、本著者らが以前開発したTimesVectorのアルゴリズムを実装し、解析手順の簡略化と様々な生物学的解釈のためのダウンストリーム解析を提供するためにいくつかの拡張を行ったものである。TimesVector-webは、https://cobi.knu.ac.kr/webserv/TimesVector-web で自由に利用できる。TimesVectorは、複数の条件下で時間的に差のある発現パターンを示す遺伝子群を検索することを目的とした時間経過型遺伝子発現解析ツールである。複数の条件と複数の時点が存在する可能性があるため、データは3次元構造(すなわち、遺伝子、時間、条件)を形成する。2 次元での部分空間クラスタリング(またはバイクラスタリング)は NP-hard であるため、3 次元クラスタリングは困難な課題である。そこで、サンプルを時間軸上で連結し、3次元データを2次元データに変換する。変換された2次元データに対して、球形クラスタリングを行い、同一クラスタ内で類似性が高く、他のクラスタ間で類似性が低い遺伝子クラスタを特定する。その後、条件ごとにサンプルを分割し、相互情報量を用いてパターンの違いを検定する。詳細についてはref.4の研究を参照されたい。

 

tutorial

https://cobi.knu.ac.kr/webserv/TimesVector-web/tutorial.php

 

f:id:kazumaxneo:20220110220831p:plain

TimesVector-web workflow. HPより

 

 

webサービス

https://cobi.knu.ac.kr/webserv/TimesVector-web/

f:id:kazumaxneo:20220110215829p:plain

 

 

解析する遺伝子発現ファイルはタブ区切りの発現行列テキスト(topからダウンロードできるデモデータ)。TimesVector-Webは時系列データおよび複数条件のRNAシーケンスデータのみに対応しているので注意する。

f:id:kazumaxneo:20220110221200p:plain

1行目はヘッダーでヘッダーの1列目'GeneID'は必須。1行目2列目以降は以下の構文に従う;'Condition_'TimePoint'(e.g, DV10_Day2)。条件と時点は"_"で区切られる。入力ファイルの実験のレプリケーツが存在する場合は、'Condition'_'TimePoint'_'Replicates'(例: DV10_Day2_rep1)とする。2行目以降には遺伝子名、発現値を記載する。

 

例えば、解析するファイルの条件が3つある場合、3つの入力ファイルが必要(全て同数の遺伝子で構成されていることが推奨されている)。

デモデータ1

f:id:kazumaxneo:20220110221800p:plain

GEOのオリジナルデータ(original/GSE11651_micro_data.txt)を条件ごとに分けて5つに分割している。

 

このデータ5つをアップロードする。シフトキーで挟んで5つ同時にアップロードした。

f:id:kazumaxneo:20220110222831p:plain

タンパク質をコードする遺伝子のみを使用するならUse only protein ~でYesを選択。データ型がマイクロアレイ型かRNA-seq型、正規化(log2 and quantile normalization。quantile normalizationはマイクロアレイでよく使われる)するかどうかを選択する。

 

ラン前にK-testボタンをクリックすると、入力データに対して適切なクラスタ数(K)が推奨される。これは入力データに対して適切なクラスタ数を推奨するためのもの。
”データ中の遺伝子発現パターンを解析するためには、適切なクラスタ数Kを選択することが重要である。ここでいう適切なクラスタ数とは、クラスタ数が最も少なく、同一クラスタ内の遺伝子の距離が小さく、異なるクラスタ間の遺伝子の距離が大きい場合である。我々のウェブサービスでは、'K-test' ボタンと 'maxK' パラメータにより、適切なクラスタ数を推奨している。このパラメータは、K検定の結果が入力データのサイズや特性によって異なる可能性があるため、ユーザが選択したいKの範囲を設定するためのものである” (マニュアルより)。

f:id:kazumaxneo:20220110223656p:plain

K-testはデータのサイズとmaxKに依存して数分かかることがある。

f:id:kazumaxneo:20220110232034p:plain

推奨値は30だった。

 

最後に生物を選択する。5種類の生物に対応している。

f:id:kazumaxneo:20220110223414p:plain

(g:Profilerのパラメータ)

 

RUNボタンを押して実行。データサイズによっては、解析が終了するまでに30分ほどかかる。

 

ショートカットキーを保存しておけば、ラン後は結果にすぐにアクセスすることができる。

f:id:kazumaxneo:20220110224106p:plain

 

出力例(demo1)

ショートカットキーは”CASE_STUDY_1”。

クラスタ数は6つとなった。

f:id:kazumaxneo:20220110225436p:plain

クラスタの図をクリックすると詳細が表示される。

Cluster pattern

f:id:kazumaxneo:20220110225511p:plain

Gene list

f:id:kazumaxneo:20220110225520p:plain

TF

f:id:kazumaxneo:20220110225532p:plain

miRNA

f:id:kazumaxneo:20220110225548p:plain

g:Profiler

f:id:kazumaxneo:20220110225617p:plain

 

引用

TimesVector-Web: A Web Service for Analysing Time Course Transcriptome Data with Multiple Conditions
Jaeyeon Jang, Inseung Hwang, Inuk Jung

Genes 2022, 13(1), 73; https://doi.org/10.3390/genes13010073

 

関連


GSEApy

 

Enrichrは哺乳類の遺伝子セットエンリッチメント解析ツールで、転写制御、パスウェイ、GOやヒトの表現型のオントロジー、薬剤で処理した細胞からのシグネチャーなどが収録されている(wiki)。GSEApyはEnrichrのPythonラッパーで、コマンドラインPython上でDEGsを調査することができる。

 

Documentation

https://gseapy.readthedocs.io/en/latest/introduction.html

A Protocol to Prepare files for GSEApy

https://gseapy.readthedocs.io/en/latest/gseapy_tutorial.html

Frequently Asked Questions - Use custom defined GMT file input in Jupyter ?

Frequently Asked Questions — GSEApy 0.10.6 documentation

 

インストール

依存

Github

#conda (bioconda)
mamba create -n gseapy -y
conda activate gseapy
mamba install -c bioconda gseapy -y

#or use pip(pypi)
$ pip install gseapy

> gseapy -h

usage: gseapy [-h] [--version]

              {gsea,prerank,ssgsea,replot,enrichr,biomart} ...

 

gseapy -- Gene Set Enrichment Analysis in Python

 

positional arguments:

  {gsea,prerank,ssgsea,replot,enrichr,biomart}

    gsea                Main GSEApy Function: run GSEApy instead of GSEA.

    prerank             Run GSEApy Prerank tool on preranked gene list.

    ssgsea              Run Single Sample GSEA.

    replot              Reproduce GSEA desktop output figures.

    enrichr             Using Enrichr API to perform GO analysis.

    biomart             Using BioMart API to convert gene ids.

 

optional arguments:

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

  --version             show program's version number and exit

 

For command line options of each command, type: gseapy COMMAND -h

 

> gseapy gsea -h

usage: gseapy gsea [-h] -d DATA -c CLS -g GMT [-t perType] [-o] [-f]

                   [--fs width height] [--graph int] [--no-plot] [-v]

                   [-n nperm] [--min-size int] [--max-size int] [-w float]

                   [-m] [-a] [-s] [-p procs]

 

optional arguments:

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

 

Input files arguments:

  -d DATA, --data DATA  Input gene expression dataset file in txt format.Same

                        with GSEA.

  -c CLS, --cls CLS     Input class vector (phenotype) file in CLS format.

                        Same with GSEA.

  -g GMT, --gmt GMT     Gene set database in GMT format. Same with GSEA.

  -t perType, --permu-type perType

                        Permutation type. Same with GSEA, choose from

                        {'gene_set', 'phenotype'}

 

Output arguments:

  -o , --outdir         The GSEApy output directory. Default: the current

                        working directory

  -f , --format         File extensions supported by Matplotlib active

                        backend, choose from {'pdf', 'png', 'jpeg','ps',

                        'eps','svg'}. Default: 'pdf'.

  --fs width height, --figsize width height

                        The figsize keyword argument need two parameters to

                        define. Default: (6.5, 6)

  --graph int           Numbers of top graphs produced. Default: 20

  --no-plot             Speed up computing by suppressing the plot output.This

                        is useful only if data are interested. Default: False.

  -v, --verbose         Increase output verbosity, print out progress of your

                        job

 

GSEA advanced arguments:

  -n nperm, --permu-num nperm

                        Number of random permutations. For calculating

                        esnulls. Default: 1000

  --min-size int        Min size of input genes presented in Gene Sets.

                        Default: 15

  --max-size int        Max size of input genes presented in Gene Sets.

                        Default: 500

  -w float, --weight float

                        Weighted_score of rank_metrics. For weighting input

                        genes. Choose from {0, 1, 1.5, 2}. Default: 1

  -m , --method         Methods to calculate correlations of ranking metrics.

                        Choose from {'signal_to_noise', 'abs_signal_to_noise',

                        't_test', 'ratio_of_classes',

                        'diff_of_classes','log2_ratio_of_classes'}. Default:

                        'log2_ratio_of_classes'

  -a, --ascending       Rank metric sorting order. If the -a flag was chosen,

                        then ascending equals to True. Default: False.

  -s , --seed           Number of random seed. Default: None

  -p procs, --threads procs

                        Number of Processes you are going to use. Default: 1

> gseapy ssgsea -h

> gseapy ssgsea -h

usage: gseapy ssgsea [-h] -d DATA -g GMT [-o] [-f] [--fs width height]

                     [--graph int] [--no-plot] [-v] [--sn normalize] [--ns]

                     [-n nperm] [--min-size int] [--max-size int] [-w weight]

                     [-a] [-s] [-p procs]

 

optional arguments:

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

 

Input files arguments:

  -d DATA, --data DATA  Input gene expression dataset file in txt format. Same

                        with GSEA.

  -g GMT, --gmt GMT     Gene set database in GMT format. Same with GSEA.

 

Output arguments:

  -o , --outdir         The GSEApy output directory. Default: the current

                        working directory

  -f , --format         File extensions supported by Matplotlib active

                        backend, choose from {'pdf', 'png', 'jpeg','ps',

                        'eps','svg'}. Default: 'pdf'.

  --fs width height, --figsize width height

                        The figsize keyword argument need two parameters to

                        define. Default: (6.5, 6)

  --graph int           Numbers of top graphs produced. Default: 20

  --no-plot             Speed up computing by suppressing the plot output.This

                        is useful only if data are interested. Default: False.

  -v, --verbose         Increase output verbosity, print out progress of your

                        job

 

Single Sample GSEA advanced arguments:

  --sn normalize, --sample-norm normalize

                        Sample normalization method. Choose from {'rank',

                        'log', 'log_rank','custom'}. Default: rank

  --ns, --no-scale      If the flag was set, don't normalize the enrichment

                        scores by number of genes.

  -n nperm, --permu-num nperm

                        Number of random permutations. For calculating

                        esnulls. Default: 0

  --min-size int        Min size of input genes presented in Gene Sets.

                        Default: 15

  --max-size int        Max size of input genes presented in Gene Sets.

                        Default: 2000

  -w weight, --weight weight

                        Weighted_score of rank_metrics. For weighting input

                        genes. Default: 0.25

  -a, --ascending       Rank metric sorting order. If the -a flag was chosen,

                        then ascending equals to True. Default: False.

  -s , --seed           Number of random seed. Default: None

  -p procs, --threads procs

                        Number of Processes you are going to use. Default: 1

> gseapy replot -h

usage: gseapy replot [-h] -i GSEA_dir [-o] [-f] [--fs width height]

                     [--graph int] [--no-plot] [-v] [-w float]

 

optional arguments:

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

 

Input arguments:

  -i GSEA_dir, --indir GSEA_dir

                        The GSEA desktop results directroy that you want to

                        reproduce the figure

  -o , --outdir         The GSEApy output directory. Default: the current

                        working directory

  -f , --format         File extensions supported by Matplotlib active

                        backend, choose from {'pdf', 'png', 'jpeg','ps',

                        'eps','svg'}. Default: 'pdf'.

  --fs width height, --figsize width height

                        The figsize keyword argument need two parameters to

                        define. Default: (6.5, 6)

  --graph int           Numbers of top graphs produced. Default: 20

  --no-plot             Speed up computing by suppressing the plot output.This

                        is useful only if data are interested. Default: False.

  -v, --verbose         Increase output verbosity, print out progress of your

                        job

  -w float, --weight float

                        Weighted_score of rank_metrics. Please Use the same

                        value in GSEA. Choose from (0, 1, 1.5, 2),default: 1

 

 

実行方法

実行するにはいくつかのファイルが必要。

 

必要なファイル1 - 発現データセット (Expression dataset)

発現行列として、タブ区切りテキストファイル(.gct)が必要(フォーマットの詳細)。

 demoデータ:GSEApy/testSet_rand1200.gct at master · zqfang/GSEApy · GitHub

f:id:kazumaxneo:20220109125341p:plain

(.GCTはマイクロアレイの時代から使われているフィーマット)

拡大した。

f:id:kazumaxneo:20220109174157p:plain

1行目はバージョン文字列。#1.2とする。

2行目には、ファイルの残りの部分に含まれるデータテーブルのサイズを示す数字を記載する。 #行数<tab>列数と書く。Gihutbの例では”#1203    18"。

3行目には、残りのファイルの各列に関連するサンプルの識別子を記載する。Gihutbの例ではNAME    Description    AA488_A1.2    AA489_A2.2    AA490_A3 ..."。

データファイルの残りの部分には、各遺伝子のデータを含める。最初の 2列のフィールドは、各遺伝子の説明(アレイではprobe ID、RNAseqでは遺伝子名)と名前(アレイでは記述(NaでもOK)、RNAseqでは遺伝子ID)を記入する。一意の識別子でないといけない。3列目以降の列には発現値を記載する。発現値はGSEApyでは正規化されないので、正規化手順(DESeq2やVoomなど)を用いて予めサンプル間正規化した値になっている必要がある。

 

必要なファイル2 - 表現型ラベル(Phenotype Labels)

表現型ラベルファイル(cls形式)は、クラスファイルまたはテンプレートファイルとも呼ばれ、表現型ラベルを定義し、それらのラベルを発現データセットのサンプルに割り当てるために使われる。clsはスペース区切りテキストファイル(もしくはタブ区切り)となっている(フォーマットの詳細)。

 demoデータ:https://github.com/zqfang/GSEApy/blob/master/tests/data/P53.cls

f:id:kazumaxneo:20220109174544p:plain

拡大した。

f:id:kazumaxneo:20220109174556p:plainCLSファイルのフォーマットは、表現型がカテゴリーか連続かによって多少異なる。カテゴリーラベルは、例えば、正常と腫瘍のように、離散的な表現型を定義する。上の写真では、MUTとWTの2つのクラスを定義していて(50 2 1)、3行目にはこのテキストが合計50記載されている (50 2 1)。

1行目には、サンプル数とクラス数を示す数字を記載する。サンプル数は、関連するRESまたはGCTデータファイルのサンプル数に対応する必要がある( サンプル数<space>クラス数<space>1)。

2行目には、各クラスのユーザーから見える名称を記載する。これは、解析レポートに表示されるクラス名。行頭には#を入れる( MUT<space>WT)。

3行目には各サンプルのクラスラベルをspace区切りで記載する。クラスラベルは、クラス名、数字、テキスト文字列のいずれかを指定する(上の写真ではテキスト)。最初のラベルは2行目のクラス名に割り当てられ、最初のラベルと異なる2番目のラベルは2番目のクラス名に割り当てられ、以下同様。この行で指定するクラスラベルの数は、最初の行で指定したサンプル数と同じにする必要がある。

注:3行目のラベルの順番でクラス名とクラスラベルの関連付けが決まる。重要なのは、3行目が左から右に処理されるとき、最初に見つけたラベルをそれが何であっても受け取り、それを2行目の最初のクラス名(これも左から右)に対応付けることである。一意のラベルの数が名前の数より多い場合はエラーになる。3行目は発現データセットに現れるサンプルを列挙したものなので、2行目のクラス名はサンプル間で最初に遭遇する順番に並べる必要がある(例:0 0 0 ... 1 1)。

 

必要なファイル3 - 遺伝子セットを記述するためのタブ区切りファイル (GMT: Gene Matrix Transposed file)

GMT形式では各行が遺伝子セットを表す(フォーマットの詳細)。

demoデータ:GSEApy/genes.gmt at master · zqfang/GSEApy · GitHub

f:id:kazumaxneo:20220109192021p:plain

GMT形式では、各遺伝子セットは、1列目が名前、2列目が説明、3列目以降が遺伝子セット内の遺伝子で記述される。GSEA は説明フィールド(2列目)を使用して、遺伝子セットの説明に対してレポート内で提供するハイパーリンクを決定する。2列目の説明が "na" の場合、GSEA は MSigDB の名前付き遺伝子セットへのリンクを提供し、説明が URL の場合、GSEA はその URL へのリンクを提供する(URLの例)。複数のGMTをマージしてカスタマイズセットを作りたい時は単純にcatで連結すればよい。

 

コマンドラインでの実行方法 (test run)

1、GSEA - GSEA の結果を生成するモジュール

入力はtxtファイル(FPKM, Expected Counts, TPM, etc.)、clsファイル、gene_setsファイル(gmtフォーマット)。フォーマットの詳細はGSEAのHPで説明されている(link)。

git clone https://github.com/zqfang/GSEApy.git
cd /GSEApy/tests/data/
gseapy gsea -d exptable.txt -c test.cls -g gene_sets.gmt -o test
  • -d    Input gene expression dataset file in txt format.Same with GSEA.
  • -c    Input class vector (phenotype) file in CLS format. Same with GSEA.
  • -g    Gene set database in GMT format. Same with GSEA.
  • -t     Permutation type. Same with GSEA, choose from {'gene_set', 'phenotype'}

 

2、prerank - Prerankツールの結果を生成するモジュール

fold changeやp valueの値でプリランクされている遺伝子リストの検定を行う。入力には相関値を持つpre-rank遺伝子リストデータセット(.rnk形式)とgene_setsファイル(gmt形式)。 prerankモジュールはGSEA pre-rankツールへのAPI。GSEAPrerankedでは表現型ファイルは必要ない。

cd /GSEApy/tests/data/
gseapy prerank -r temp.rnk -g genes.gmt -o test

test/

f:id:kazumaxneo:20220109145418p:plain

gseapy.prerank.gene_sets.report.csv

f:id:kazumaxneo:20220109145600p:plain

DvA_UpIN_A.prerank.pdf

f:id:kazumaxneo:20220109145440p:plain

 

 

3、ssgsea - single sample GSEA(ssGSEA)解析を行うモジュール

Single-sample GSEA (ssGSEA) は Gene Set Enrichment Analysis (GSEA) の拡張で、サンプルと遺伝子セットの各ペアに対して別々のエンリッチメントスコアを計算する。 各ssGSEAののエンリッチメントスコアは、特定の遺伝子セット内の遺伝子がサンプル内でどの程度協調的にアップレギュレートまたはダウンレギュレートされているかを表す。

入力はpd.Series(遺伝子名でインデックス付け)、またはpd.DataFrame(GCTファイルを含む)、発現値、GMTファイルを想定している。複数サンプル入力の場合、ssGSEAはgctフォーマットも再作成する。遺伝子セットのssGSEAエンリッチメントスコアはD. Barbie et al 2009に記述されている。

cd /GSEApy/tests/data/
gseapy ssgsea -d P53.txt -g genes.gmt -o test

f:id:kazumaxneo:20220109145947p:plain

 

4、replot - GSEAデスクトップのレポートを生成する

cd /GSEApy/tests/data/
gseapy replot -i edb/ -o test

試した時はエラーになった。

 

他にも遺伝子IDを変換するbiomartコマンドやenrichrを実行するenrichrコマンドが利用できます。

GSEA software download ; http://www.gsea-msigdb.org/gsea/downloads.jsp

f:id:kazumaxneo:20220109195458p:plain

 

マニュアルでは、pythonのコンソール上(jupyterなど)で利用する方法が詳しく説明されています。興味がある方は確認して下さい。

引用

https://github.com/zqfang/GSEApy

”If you use gseapy in your research, you should cite the original ``GSEA`` and ``Enrichr`` paper.”

 

参考

 

関連