macでインフォマティクス

macでインフォマティクス

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

タンパク質配列をゲノム配列に対してintron (gap) awareで高速にアラインメントする Miniprot

 

Githubより

Miniprotは、タンパク質配列をゲノムに対してアフィンギャップ・ペナルティ、スプライシング、フレームシフ トでアライメントする。Miniprotは、他の既知の種の遺伝子を用いて、新しい種のタンパク質コード遺伝子をアノテーションすることを主な目的としている。MiniprotはGeneWise(EMBL-EBIサービス)やExonerateと機能的に似ているが、全ゲノムに対してタンパク質をマップすることができ、残基アライメントステップが非常に高速です。
Miniprotは距離の離れたホモログのマッピングには最適化されていない。なぜなら、遠いホモログは遺伝子アノテーションの情報量が少ないからです。それでも、パフォーマンスを犠牲にすることなく、より高い感度を得るためにseedingパラメータを調整することが可能です。

Manual

https://lh3.github.io/miniprot/miniprot.html

 

 

 

インストール

condaを使ってubuntu18に導入した。

Github

#conda(link)
mamba install -c bioconda miniprot -y

#from source
git clone https://github.com/lh3/miniprot
cd miniprot && make

> miniprot

Usage: miniprot [options] <ref.fa> <query.faa> [...]

Options:

  Indexing:

    -k INT       k-mer size [6]

    -s INT       submer size (density: 1/(2*(k-s)+1)) [4]

    -b INT       bits per block [8]

    -d FILE      save index to FILE []

  Mapping:

    -S           no splicing (applying -G1k -J1k -e1k)

    -c NUM       max k-mer occurrence [50000]

    -G NUM       max intron size [200k]

    -w FLOAT     weight of log gap penalty [0.75]

    -n NUM       minimum number of syncmers in a chain [5]

    -m NUM       min chaining score [0]

    -l INT       k-mer size for the second round of chaining [5]

    -e NUM       max extension for 2nd round of chaining and alignment [10k]

    -p FLOAT     min secondary-to-primary score ratio [0.7]

    -N NUM       consider at most INT secondary alignments [50]

  Alignment:

    -O INT       gap open penalty [11]

    -E INT       gap extension (a k-long gap costs O+k*E) [1]

    -J INT       intron open penalty [31]

    -C INT       penalty for non-canonical splicing [11]

    -F INT       penalty for frameshifts or in-frame stop codons [17]

    -B INT       end bonus [5]

  Input/output:

    -t INT       number of threads [4]

    --gff        output in the GFF3 format

    -P STR       prefix for IDs in GFF3 [MP]

    -u           print unmapped query proteins

    --outn=NUM   output up to min{NUM,-N} alignments per query [100]

    -K NUM       query batch size [2M]x

 

 

実行方法

ゲノム配列(もしくはindexされたファイル)とタンパク質配列を指定する。8スレッド指定。

miniprot -t8 ref-file protein.faa > output.paf
  •  -t      number of threads [4] 

miniprotのインデックス作成は低速でメモリを消費するため、事前にインデックスを作っておくことが推奨されている。出力はprotein PAF format。マニュアルで説明されている。

 

--gffオプションを使うことでGFF3形式で出力することもできる。

miniprot -t8 --gff -d ref.mpi ref.fna > out.gff
  •  --gff    output in the GFF3 format 

 

シロイヌナズナの全タンパク質配列をシロイヌナズナのリファレンスゲノム配列にアラインしたところ、かかった時間は21秒だった(TR3990x, 20スレッド、pre-indexing無し)。

引用

https://github.com/lh3/miniprot

 

関連


 

 

Nanopolishのcall-methylationおよびeventalignモジュールを最適化して再実装した f5c

 

 

 ナノポアシーケンスにより、ポイントオブケア診断や現場でのジェノタイピングなど、携帯可能なリアルタイムシーケンスアプリケーションが可能になる。このような成果を得るためには、生のナノポアシグナルデータを解析するための効率的なバイオインフォマティクスアルゴリズムが必要になる。しかし、生のナノポアシグナルを生物学的参照配列と比較することは、計算上複雑なタスクである。ABEA(Adaptive Banded Event Alignment)と呼ばれる動的プログラミングアルゴリズムは、シーケンスデータを研磨し、DNAメチル化の測定など、非標準のヌクレオチドを特定するのに重要なステップである。ここでは、ABEAアルゴリズムの実装(f5cと呼ぶ)を並列化し、ヘテロジニアスCPU-GPUアーキテクチャ上で効率的に動作するように最適化した。
メモリ、計算、CPUとGPU間の負荷分散を最適化することにより、f5cはNanopolishソフトウェアパッケージのオリジナルのCPUのみのABEA実装を最適化したバージョンよりも約3-5倍高速に動作することを実証した。また、GPUを搭載した組込みSoC(System on Chip)を用いて、f5cがオンザフライでDNAメチル化検出を可能にすることも示した。
本研究は、複雑なゲノム解析が軽量なコンピューティングシステムで実行できることを示すだけでなく、高性能コンピューティング(HPC)にも貢献するものである。f5cのソースコードGPU最適化ABEAは、https://github.com/hasindu2008/f5c で公開されている。

 

Full Documentation

https://hasindu2008.github.io/f5c/docs/overview

Commnad reference

https://hasindu2008.github.io/f5c/docs/commands

 

インストール

ubuntu18にてCPU版のバイナリを試した。

Github

VERSION=v1.1
wget "https://github.com/hasindu2008/f5c/releases/download/$VERSION/f5c-$VERSION-binaries.tar.gz" && tar xvf f5c-$VERSION-binaries.tar.gz && cd f5c-$VERSION/

> ./f5c_x86_64_linux -h

Usage: f5c <command> [options]

 

command:

         index               Build an index for accessing the base sequence and raw signal for a given read ID (optimised nanopolish index)

         call-methylation    Classify nucleotides as methylated or not (optimised nanopolish call-methylation)

         meth-freq           Calculate methylation frequency at genomic CpG sites (optimised nanopolish calculate_methylation_frequency.py)

         eventalign          Align nanopore events to reference k-mers (optimised nanopolish eventalign)

         freq-merge          Merge calculated methylation frequency tsv files

         resquiggle          Align raw signals to basecalled reads

 

See the manual page for details (`man ./docs/f5c.1' or https://f5c.page.link/man).

> ./f5c_x86_64_linux index

[f5c index]  not enough arguments

 

Usage: f5c index [OPTIONS] -d fast5_directory reads.fastq

       f5c index [OPTIONS] --slow5 signals.blow5 reads.fastq

 

Build an index for accessing the base sequence (fastq/fasta) and raw signal (fast5/slow5) for a given read ID. f5c index is an extended and optimised version of nanopolish index by Jared Simpson

 

  -h                display this help and exit

  -d STR            path to the directory containing fast5 files. This option can be given multiple times.

  -s STR            the sequencing summary file

  -f STR            file containing the paths to the sequencing summary files (one per line)

  -t INT            number of threads used for bgzf compression (makes indexing faster)

  --iop INT         number of I/O processes to read fast5 files (makes fast5 indexing faster)

  --slow5 FILE      slow5 file containing raw signals

  --skip-slow5-idx  do not build the .idx for the slow5 file (useful when a slow5 index is already available)

  --verbose INT     verbosity level

  --version         print version

 

See the manual page for details (`man ./docs/f5c.1' or https://f5c.page.link/man).

> ./f5c_x86_64_linux call-methylation 

[meth_main::INFO] Default methylation tsv output format is changed from f5c v0.7 onwards to match latest nanopolish output. Set --meth-out-version=1 to fall back to the old format.

Usage: f5c call-methylation [OPTIONS] -r reads.fa -b alignments.bam -g genome.fa

       f5c call-methylation [OPTIONS] -r reads.fa -b alignments.bam -g genome.fa --slow5 signals.blow5

 

basic options:

   -r FILE                    fastq/fasta read file

   -b FILE                    sorted bam file

   -g FILE                    reference genome

   -w STR                     limit processing to a genomic region specified as chr:start-end or a list of regions in a .bed file

   -t INT                     number of processing threads [8]

   -K INT                     batch size (max number of reads loaded at once) [512]

   -B FLOAT[K/M/G]            max number of bases loaded at once [5.0M]

   -h                         help

   -o FILE                    output to file [stdout]

   -x STR                     parameter profile to be used for better performance (always applied before other options)

                              e.g., laptop, desktop, hpc; see https://f5c.page.link/profiles for the full list

   --iop INT                  number of I/O processes to read fast5 files [1]

   --slow5 FILE               read from a slow5 file

   --min-mapq INT             minimum mapping quality [20]

   --secondary=yes|no         consider secondary mappings or not [no]

   --verbose INT              verbosity level [0]

   --version                  print version

 

advanced options:

   --skip-ultra FILE          skip ultra long reads and write those entries to the bam file provided as the argument

   --ultra-thresh INT         threshold to skip ultra long reads [100000]

   --skip-unreadable=yes|no   skip any unreadable fast5 or terminate program [yes]

   --kmer-model FILE          custom nucleotide k-mer model file (format similar to test/r9-models/r9.4_450bps.nucleotide.6mer.template.model)

   --meth-model FILE          custom methylation k-mer model file (format similar to test/r9-models/r9.4_450bps.cpg.6mer.template.model)

   --meth-out-version INT     methylation tsv output version (set 2 to print the strand column) [2]

   --min-recalib-events INT   minimum number of events to recalbrate (decrease if your reads are very short and could not calibrate) [200]

 

developer options:

   --print-events=yes|no      prints the event table

   --print-banded-aln=yes|no  prints the event alignment

   --print-scaling=yes|no     prints the estimated scalings

   --print-raw=yes|no         prints the raw signal

   --debug-break INT          break after processing the specified no. of batches

   --profile-cpu=yes|no       process section by section (used for profiling on CPU)

   --write-dump=yes|no        write the fast5 dump to a file or not

   --read-dump=yes|no         read from a fast5 dump file or not

 

See the manual page for details (`man ./docs/f5c.1' or https://f5c.page.link/man).

 

 

 

テストラン

https://hasindu2008.github.io/f5c/docs/example-usageに従う。公式チュートリアルではマッピングのコマンドや計算資源がリッチな環境でのコマンド例も書かれている。

 

1,Download

wget -O f5c_na12878_test.tgz "https://f5c.page.link/f5c_na12878_test"
tar xf f5c_na12878_test.tgz

 

2、indexingとmethylation-call

index作成には、fast5(かもしくはslow5 format)のファイルを含むディレクトリとbasecallされたfastqを指定する。次のmethylation-callでは、basecallされたfastqをリファレンスにマッピングして得たbamファイル、fasta形式のリファレンス配列、basecallされたfastqを指定する。

#1 indexing (nanopolish indexの出力と同等)
f5c index -d chr22_meth_example/fast5_files chr22_meth_example/reads.fastq

#2 ゲノムCpG部位のメチル化有無の分類を行う(Nanopolishの出力と同等)
f5c call-methylation -b chr22_meth_example/reads.sorted.bam -g chr22_meth_example/humangenome.fa -r chr22_meth_example/reads.fastq > chr22_meth_example/result.tsv

#3 f5c call-methylationで生成されたメチル化コールを含むtsvファイルから、ゲノムCpG部位のメチル化頻度を計算する
f5c meth-freq -i chr22_meth_example/result.tsv > chr22_meth_example/freq.tsv

#4複数のメチル化頻度tsvファイル(f5c meth-freqの出力)がある場合、1つのtsvファイルにマージする
f5c freq-merge

出力

 

3、event alignment

f5c eventalignは、生のナノポアシグナル(イベント)を参照k-merにアラインさせる。

f5c eventalign -b chr22_meth_example/reads.sorted.bam -g chr22_meth_example/humangenome.fa -r chr22_meth_example/reads.fastq > chr22_meth_example/events.tsv

Nanopolishの出力と同等。

 

4、f5c resquiggle

f5c resquiggleはrawシグナルをbasecalled readsにアライメントする(現在開発中で、将来のバージョンで出力は変更される可能性がある)。

f5c_x86_64_linux resquiggle chr22_meth_example/reads.fastq chr22_meth_example/reads.blow5

デフォルトの出力はTSV形式

https://hasindu2008.github.io/f5c/docs/output

 

引用

GPU accelerated adaptive banded event alignment for rapid comparative nanopore signal analysis
Hasindu Gamaarachchi, Chun Wai Lam, Gihan Jayatilaka, Hiruna Samarakoon, Jared T. Simpson, Martin A. Smith & Sri Parameswaran 
BMC Bioinformatics volume 21, Article number: 343 (2020) 

 

関連


 

 

系統マーカー遺伝子に分類群を割り当てる AmphoraNetと結果を視覚化するAmphoraVizu

 

 メタゲノム解析はここ数年、目覚しい発展を遂げた。今日、遺伝子配列決定の専門家だけでなく、他の専門分野の多くの研究室が、臨床サンプルや環境サンプルから得られたDNA配列を解析する必要がある。メタゲノム解析データの系統解析は、生物学者やバイオインフォマティシャンにとって重要な課題となっている。AMPHORAプログラムとそのワークフロー版は、メタゲノム・データに対して信頼性の高い系統解析結果をもたらす、一般に公開されているソフトウェアの一例である。
 このウェブサーバーはAMPHORA2ワークフローに基づいており、入力されたメタゲノム試料で見つかった各系統マーカー遺伝子に対して、確率的に重み付けされた分類群を割り当てることができる使い勝手の良いウェブサーバーである。分子生物学者の多くは、BLASTプログラムとそのクローンをローカルにインストールするのではなく、公共のウェブサーバー上で使用しているため、このバージョンでは、AMPHORA2スイートのすべてのコンポーネントの時間のかかるインストールやLinux環境の専門知識が必要ない。ウェブサーバーは、登録は不要でhttp://amphoranet.pitgroup.org で自由に利用できる。

 

HPより

AmphoraNetは、31の細菌と104の古細菌のタンパク質コードマーカー遺伝子を用いて、メタゲノムおよびゲノムの系統分類を行います。これらのほとんどはシングルコピー遺伝子であるため、メタゲノムショットガンシーケンスデータからバクテリアおよび古細菌群集の分類学的構成を推定するのに適しています。

 

Online Tools

https://pitgroup.org/tools/

 

webサービス

https://pitgroup.org/amphoranet/にアクセスする。

 


FASTAファイル を指定し(最大300MB)、入力データの種類、検索したいマーカー遺伝子の種類を指定する。

Check values "ボタンをクリックする。

問題がなければ、Schedule jobボタンをクリックする。

 

出力

3つのファイル(入力ファイル、マーカー遺伝子ファイル、結果ファイル)をダウンロードできる。

 

zip形式でダウンロードできるマーカー遺伝子ファイルには、発見された系統マーカー遺伝子のタンパク質配列が含まれている。

 

taxonomic profiling - 同定された各系統マーカー遺伝子のphylotypesが含まれるタブ区切りテキストファイル。

(カッコ内の数値は信頼度スコア)

 

このマーカー遺伝子ごとのtaxonomic profilingテキストは、AmphoraVizuwebサーバーにアップロードすることで、任意のtaxonomic rankの割合を視覚化できる。

 

AmphoraVizu

https://pitgroup.org/amphoravizu/にアクセスする。

taxonomic profilingテキストを指定する。詳細から、表示モード、信頼度の下限、平均シェア、タイトルを選択できる。

 

図はDownload as HTMLボタンからHTML形式でダウンロードする。

 

 

引用

AmphoraNet: the webserver implementation of the AMPHORA2 metagenomic workflow suite
Csaba Kerepesi, Dániel Bánky, Vince Grolmusz

Gene. 2014 Jan 10;533(2):538-40

 

Visual Analysis of the Quantitative Composition of Metagenomic Communities: the AmphoraVizu Webserver
Csaba Kerepesi, Balázs Szalkai & Vince Grolmusz 
Microbial Ecology volume 69, pages695–697 (2015)

 

真菌のITSやコアタンパク質コード遺伝子を使った系統解析を自動で実行する UFCG pipeline

 

UFCG pipelineを使うと、真菌のITSやコアタンパク質を使った系統解析を自動で実行できます。簡単にですが、使い方を確認しておきます。

 

 

Manual

https://ufcg.steineggerlab.com/ufcg/manual

 

インストール

依存

  • Java Runtime Environment with a version higher than 8

Github

ubuntu18に導入した。

mamba create -n ufcg -c bioconda -c conda-forge openjdk=8 augustus mmseqs2 mafft iqtree
conda activate ufcg
wget -O UFCG.zip https://github.com/endixk/ufcg/releases/latest/download/UFCG.zip
unzip UFCG.zip && cd UFCG
java -jar UFCG.jar

#docker(link)
docker pull endix1029/ufcg:latest
docker run -it endix1029/ufcg:latest
cd UFCG
java -jar UFCG.jar

> java -jar UFCG.jar

$ java -jar UFCG.jar

 

    __  __ _____ _____ _____

   / / / // ___// ___// ___/

  / / / // /_  / /   / / __

 / /_/ // __/ / /___/ /_/ /

 \____//_/    \____/\____/ v1.0

 

 

 USAGE : java -jar UFCG.jar <module> [...]

 

 

 Available Modules

 Module         Description

 profile        Extract UFCG profile from genome

 profile-rna    Extract UFCG profile from RNA-seq transcriptome

 profile-pro    Extract UFCG profile from proteome

 train          Train and generate sequence model

 align          Produce sequence alignments from UFCG profiles

 tree           Build maximum likelihood tree with UFCG profiles

 prune          Rebuild UFCG tree or single gene trees

 

 

 Miscellaneous

 Argument       Description

 --info         Print program information

 --core         Print core gene list

 

 

 General options

 Argument       Description

 -h, --help     Print this manual

 -v, --verbose  Make program verbose

 --nocolor      Remove ANSI escapes from standard output

 --notime       Remove timestamp in front of the prompt string

 --developer    Activate developer mode (For testing or debugging)

 

 

 

実行方法

UFCG pipelineを使うには、系統推定に使いたいゲノムのfastaファイルと、任意でメタデータファイルが必要。メタデータファイルは、ゲノムの分類記号を表す 7 つのエントリを受け取ることができる。

簡易バージョン;meta_simple.tsv

簡易版ではfilename、label、accessionを記載する。

 

フルバージョン;meta_full.tsv

フルバージョンではさらにtaxon_nameなどを記載する(link)。

 

ゲノムのfasta形式ファイルは、ゲノムアセンブリファイルのみ含むディレクトリに配置する。ファイルの拡張子は統一する。(.fa, .fna, .fasta, ...)。ファイルのリストは、指定されたメタデータ1列目のfilenameに対応していなければならない。メタデータのみあってゲノムファイルがなかったり、ゲノムファイルがあってメタデータのリストには含まれていないとエラーになるので注意。

 

1、対話式で実行。profileコマンドを使う。

 java -jar UFCG.jar profile -u

seqに含まれるGCA_000697725.1.fnaが含まれているとエラーになった。除くとランできた。最終的に、指定した出力ディレクトリにゲノムの数だけ.ucgファイルが書き出される。

 

2、系統推定にはtreeコマンドを使う。1の出力ディレクトリと、出力ツリーファイルで系統樹の葉の名前とするメタデータ列を指定する。

java -jar UFCG.jar tree -i profile_results/ -l label
  • -i     Locate the path of the input .ucg profiles to align and infer tree
  • -l     Name the leaves of the phylogenetic tree from the metadata 

以下のメタデータ列に対応している(manualより)。上ではlabelを指定した。

uid : Include unique integer ID
acc : Include accession number
label : Include full label
taxon : Include taxon name
strain : Include strain name
taxonomy : Include taxonomic relationship

 

出力例

 

系統マーカータンパク質配列として、このレポジトリのseq/に含まれている配列を使用できます。これは、UFCGの定義するコアタンパク質になります。

また、UFCGのHPからは、分類がvalidな1000以上の真菌からのコア遺伝子のタンパク質のMSAとHMMプロファイルをダウンロードできます(手動で集める必要がないというのが重要)。ダウンロードしたMSAに自分のゲノムのタンパク質配列を加えることで、1000以上の真菌の系統解析を素早く実行することができます(ただし計算時間は長くなります)。ゲノムを決めた真菌の系統がわからない時は、このような方法が使えます。一方で、ここで紹介したUFCG pipelineは、手元にゲノムのリストがある場合に活用できます。

引用

UFCG: database of universal fungal core genes and pipeline for genome-wide phylogenetic analysis of fungi
Dongwook Kim,  Cameron L.M. Gilchrist,  Jongsik Chun, Martin Steinegger

bioRxiv, Posted August 17, 2022.

 

関連


 

ショートリードからの株レベルメタゲノムアセンブリを行う StrainXpress

 

 次世代シーケンサーを用いたメタゲノム解析により、長時間の培養を必要とせず、特徴的な生息環境にある微生物を同定することが可能になった。重要なことは、薬剤耐性、病原性、環境との相互作用など、臨床に関連する現象が種内で既に変化している可能性があることである。そのため、シーケンシングリードから個々のゲノムを種のレベルだけでなく、株のレベルでも再構築することが現在の大きな課題となっている。しかし、ある種の株は微量のバリアントしか違わないため、それらを区別することは困難である。近年、かなりの進歩が見られるものの、関連するアプローチはこれまで断片的なものにとどまっている。本発表では、次世代シーケンサーのリードデータからstrain aware metagenome assemblyのための包括的なソリューションであるStrainXpressを紹介する。StrainXpressは、最大1000系統以上のメタゲノムから系統特異的なゲノムを再構成することができ、また、被覆度の低い系統にもうまく対処できることを実験的に証明した。その結果、全データセットで平均26.75%(第一四分位:18.51%、中央値:26.60%、第三四分位:35.05%)の菌株特異的な配列が再構成され、現在の最先端アプローチのそれを上回った。

 

インストール

ubuntu18にcondaを使って導入した。

mamba create -n strainxpress
conda activate strainxpress
mamba install -c bioconda python=3.6 scipy pandas minimap2 -y

git clone https://github.com/kangxiongbin/StrainXpress.git
cd StrainXpress
sh install.sh

python scripts/strainxpress.py -h

usage: strainxpress.py [-h] [-fq FQ] [-t THREADS] [-size SIZE]

                       [-insert_size INSERT_SIZE]

                       [-average_read_len AVERAGE_READ_LEN]

                       [-split_nu SPLIT_NU] [-fast]

 

%prog -fq <input fq file> This program is used to cluster the reads that stem

from identical species. # need to install panda

 

optional arguments:

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

  -fq FQ                The input fq file.

  -t THREADS            The number of threads when run strainxpress. Default

                        is 10.

  -size SIZE            The maximum size of the cluster. Default is 15000.

  -insert_size INSERT_SIZE

                        The length of insert size of reads. Default is 450.

  -average_read_len AVERAGE_READ_LEN

                        The average length of reads. Default is 250.

  -split_nu SPLIT_NU    Split the fq file into several files. Default is 8.

  -fast                 When fq file is big or have some high coverage

                        bacteria, suggest use the fast model. Defualt is don't

                        use fast.

 

 

テストラン

cd example/
python ../scripts/strainxpress.py -fq all_reads.fq

出力

 

引用
StrainXpress: strain aware metagenome assembly from short reads

Xiongbin Kang, Xiao Luo , Alexander Schönhuth

Nucleic Acids Res. 2022 Jul 1;gkac543

NCBI SRA Run Selectorを使う

 

SRA Run Selectorは、SRAに保存されている大規模なランのセットを取り出し、どのランを解析に使用するかを絞り込み、結果をメタデータとしてダウンロードすることができる。

 

NCBI-Hackathons/SRA Run Selector Tutorial

https://github.com/NCBI-Hackathons/ncbi-cloud-tutorials/blob/master/SRA%20tutorials/tutorial_SRA_run_selector.md

NCBI Minute: Using the SRA RunSelector to Find NGS Datasets


試してみる。biofilmのメタゲノムデータセットを探す。NCBI SRAでbiofimとタイプ。

 

56032ヒットした。

右上のtop organismでは、metagenomeの検索ではあるが、biofilm metagenomeというtaxonomyが18410ヒットしている。ほかにも口腔biofilm(プラーク?)などが出てきてるが、自然環境のデータに興味があったので、biofilm metagenomeをクリックした。"metagenome"を選んだので、biofilmだけではたくさんヒットしていたかもしれないメタアンプリコンのデータは効果的に除外されるはずである。

 

結果。まだ絞れそうだが、Send results to Run selectorが出てきたのでクリック。 

このSend results to Run selectorは、検索結果の画面で一度でも絞り込みを行うと出てくる。

 

Run selectorに移動した。

注;dbGaPのデータは許可されているアカウントでログインしていないと扱えない。また、アクセスする前にコールドストレージから取得する必要があるものもある。その場合、クラウドベースのコールドストレージを経由して転送することになる。

 

Run selectorでは、関心のあるRunデータを絞り込むために様々なフィルターを適用できる(表示されるフィルターはリストによって変化する)。ここではPlatformをクリックした。

Platfoprmが表示された。ILLUMINAが多いが、他のシークエンスプラットフォームも見つかる。ILLUMINAを選択した。

 

ILLUMINAの16666ランに絞り込まれた。

(フィルターは複数重ねて適用できる。)

 

ランのチェックをつけると、直接絞り込むことができる。下の写真のように4つだけ選択してみた。選択したデータの個数と総塩基数などが写真上部のSelectedに表示されている。Bytesが194.24Mb、Basesが219.70 Mとなった。

 

今回は絞り込んだ16666ランのデータをダウンロードするため、Xマーク左をクリックして全部選択する。(Xマーククリックで解除)

全部選択するとSelectedが16666となった。

 

SlectedのMetadataをクリックするとメタデータテーブル全体がダウンロードされる。

 

SlectedのAccesion listをクリックするとIDだけがダウンロードされる。

 

リストが手に入ったら、fastq-dumpを使ってデータを取得する。ここでは目的の細菌種が決まっているので、そのようなデータが含まれるのか調べるために少数だけリードを取得する。数が多いと時間がかかるので、GNU parallelで10並列し、30万リードずつだけ取得する。-Nと-Xはイルミナフローセルの座標に相当する(スポットの指定)。0から取得するのは品質が問題になる気があする?ので 10000から開始して310000まで。

#--dry-runを付けてチェック。OKなら外す
cat Accesion_list| parallel -j 10 --dry-run 'fastq-dump --split-files -N 10000 -X 310000 --gzip {}'
  • -N      Minimum spot id
  • -X      Maximum spot id
  • --split-files    Write reads into separate files. Read number will be suffixed to the file name. NOTE! The `--split-3` option is  recommended. In cases where not all spots have the same number of reads,  this option will produce files that WILL
    CAUSE ERRORS in most programs which process split pair fastq files.
  • --gzip    Compress output using gzip: deprecated, not recommended

1000サンプルで2~3時間ほどかかった。

 

結果はkraken2などで分析する。

#!/bin/sh 
for file in `\find *_1.fastq.gz -maxdepth 1 -type f`;
do
  fastq1=$file;
  fastq2=${file%_1.fastq.gz}_2.fastq.gz;
kraken2 --db <path>/<to>/kraken2-database --threads 20 --gzip-compressed --paired --report ${file%_1.fastq.gz}_kraken2_report.txt --use-names --classified-out seqs#.fq $fastq1 $fastq2 > ${file%_1.fastq.gz}_output.txt;
  rm seqs_1.fq seqs_2.fq ${file%_1.fastq.gz}_output.txt;
bracken -d <path>/<to>/kraken2-database -i ${file%_1.fastq.gz}_kraken2_report.txt -r 100 -l S -t 5 -o ${file%_1.fastq.gz}_bracken;
python KrakenTools/kreport2krona.py -r ${file%_1.fastq.gz}_kraken2_report_bracken.txt -o ${file%_1.fastq.gz}_bracken_convert;
done

#kronaで可視化
ktImportText *_bracken_convert -o all.html

 

NCBIの動画でわかりやすく説明されています。視聴してみてください。ただ2017年の動画でインターフェイスが古いものになっています。

引用

ncbi-cloud-tutorials/tutorial_SRA_run_selector.md at master · NCBI-Hackathons/ncbi-cloud-tutorials · GitHub

 

参考

Biostars

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

 

Demultiplexingを行う fgbioのDemuxFastqsコマンド

 

fgbioはディープシーケンシングデータを扱うためのコマンドラインツールキット。リードレベルのデータ(FASTQ、SAM、BAMなど)やバリアントレベルのデータ(VCF、BCFなど)を操作する。特に次のようなものを提供することに重点を置いている(Githubより)。

  • 堅牢で、よくテストされたツール
  • 使いやすいコマンドライン
  • 各ツールの明確かつ徹底したドキュメント
  • コミュニティとクライアントの利益となるオープンソース開発

fgbioの中のfgbio DemuxFastqsサブコマンドは、FASTQsのdemultiplexing (デマルチプレックス)を行う。また、オプションでUMIを抽出する。

 

Tool section

http://fulcrumgenomics.github.io/fgbio/tools/latest/

fgbio DemuxFastqs

http://fulcrumgenomics.github.io/fgbio/tools/latest/DemuxFastqs.html

 

インストール

ubuntu18に導入した。

Github

mamba install -c bioconda fgbio -y

> fgbio DemuxFastqs

$ fgbio DemuxFastqs

USAGE: fgbio [fgbio arguments] [command name] [command arguments]

Version: 2.0.2

------------------------------------------------------------------------------------------------------------------------

 

fgbio Arguments:

------------------------------------------------------------------------------------------------------------------------

-h true|false, --help=true|false

                              Display the help message. [Default: false]. 

--async-io=true|false     Use asynchronous I/O where possible, e.g. for SAM and BAM files. [Default:

                              false]. 

--version=true|false      Display the version number for this tool. [Default: false]. 

--compression=Int             Default GZIP compression level, BAM compression level. [Default: 5]. 

--tmp-dir=DirPath             Directory to use for temporary files. [Default:

                              /var/folders/9y/gqf42hb548178qbs0mm2r78w0000gn/T]. 

--log-level=LogLevel          Minimum severity log-level to emit. [Default: Info]. Options: Debug, Info,

                              Warning, Error, Fatal.

--sam-validation-stringency=ValidationStringency

                              Validation stringency for SAM/BAM reading. [Default: SILENT]. Options:

                              STRICT, LENIENT, SILENT.

 

DemuxFastqs

------------------------------------------------------------------------------------------------------------------------

Performs sample demultiplexing on FASTQs.

 

The sample barcode for each sample in the sample sheet will be compared against the sample barcode bases extracted from

the FASTQs, to assign each read to a sample. Reads that do not match any sample within the given error tolerance will

be placed in the 'unmatched' file.

 

The type of output is specified with the '--output-type' option, and can be BAM ('--output-type Bam'), gzipped FASTQ

('--output-type Fastq'), or both ('--output-type BamAndFastq').

 

For BAM output, the output directory will contain one BAM file per sample in the sample sheet or metadata CSV file,

plus a BAM for reads that could not be assigned to a sample given the criteria. The output file names will be the

concatenation of sample id, sample name, and sample barcode bases (expected not observed), delimited by '-'. A metrics

file will also be output providing analogous information to the metric described SampleBarcodeMetric

(http://fulcrumgenomics.github.io/fgbio/metrics/latest/#samplebarcodemetric).

 

For gzipped FASTQ output, one or more gzipped FASTQs per sample in the sample sheet or metadata CSV file will be

written to the output directory. For paired end data, the output will have the suffix '_R1.fastq.gz' and '_R2.fastq.gz'

for read one and read two respectively. The sample barcode and molecular barcodes (concatenated) will be appended to

the read name and delimited by a colon. If the '--illumina-standards' option is given, then the output read names and

file names will follow the Illumina standards described here

(https://help.basespace.illumina.com/articles/tutorials/upload-data-using-web-uploader/).

 

The output base qualities will be standardized to Sanger/SAM format.

 

FASTQs and associated read structures for each sub-read should be given:

 

  * a single fragment read should have one FASTQ and one read structure

  * paired end reads should have two FASTQs and two read structures

  * a dual-index sample with paired end reads should have four FASTQs and four read structures given: two for the two

    index reads, and two for the template reads.

 

If multiple FASTQs are present for each sub-read, then the FASTQs for each sub-read should be concatenated together

prior to running this tool (ex. 'cat s_R1_L001.fq.gz s_R1_L002.fq.gz > s_R1.fq.gz').

 

(Read structures)https://github.com/fulcrumgenomics/fgbio/wiki/Read-Structures are made up of '<number><operator>'

pairs much like the 'CIGAR' string in BAM files. Four kinds of operators are recognized:

 

  1. 'T' identifies a template read

  2. 'B' identifies a sample barcode read

  3. 'M' identifies a unique molecular index read

  4. 'S' identifies a set of bases that should be skipped or ignored

 

The last '<number><operator>' pair may be specified using a '+' sign instead of number to denote "all remaining bases".

This is useful if, e.g., fastqs have been trimmed and contain reads of varying length. Both reads must have template

bases. Any molecular identifiers will be concatenated using the '-' delimiter and placed in the given SAM record tag

('RX' by default). Similarly, the sample barcode bases from the given read will be placed in the 'BC' tag.

 

Metadata about the samples should be given in either an Illumina Experiment Manager sample sheet or a metadata CSV

file. Formats are described in detail below.

 

The read structures will be used to extract the observed sample barcode, template bases, and molecular identifiers from

each read. The observed sample barcode will be matched to the sample barcodes extracted from the bases in the sample

metadata and associated read structures.

 

Sample Sheet

------------

 

The read group's sample id, sample name, and library id all correspond to the similarly named values in the sample

sheet. Library id will be the sample id if not found, and the platform unit will be the sample name concatenated with

the sample barcode bases delimited by a '.'.

 

The sample section of the sample sheet should contain information related to each sample with the following columns:

 

  * Sample_ID: The sample identifier unique to the sample in the sample sheet.

  * Sample_Name: The sample name.

  * Library_ID: The library Identifier. The combination sample name and library identifier should be unique across the

    samples in the sample sheet.

  * Description: The description of the sample, which will be placed in the description field in the output BAM's read

    group. This column may be omitted.

  * Sample_Barcode: The sample barcode bases unique to each sample. The name of the column containing the sample

    barcode can be changed using the '--column-for-sample-barcode' option. If the sample barcode is present across

    multiple reads (ex. dual-index, or inline in both reads of a pair), then the expected barcode bases from each read

    should be concatenated in the same order as the order of the reads' FASTQs and read structures given to this tool.

 

Metadata CSV

------------

 

In lieu of a sample sheet, a simple CSV file may be provided with the necessary metadata. This file should contain the

same columns as described above for the sample sheet ('Sample_ID', 'Sample_Name', 'Library_ID', and 'Description').

 

Example Command Line

--------------------

 

As an example, if the sequencing run was 2x100bp (paired end) with two 8bp index reads both reading a sample barcode,

as well as an in-line 8bp sample barcode in read one, the command line would be

 

  --inputs r1.fq i1.fq i2.fq r2.fq --read-structures 8B92T 8B 8B 100T \

      --metadata SampleSheet.csv --metrics metrics.txt --output output_folder

 

Output Standards

----------------

 

The following options affect the output format:

 

  1. If '--omit-fastq-read-numbers' is specified, then trailing /1 and /2 for R1 and R2 respectively, will not be

     appended to e FASTQ read name. By default they will be appended.

  2. If '--include-sample-barcodes-in-fastq' is specified, then sample barcode will replace the last field in the first

     comment in the FASTQ header, e.g. replace 'NNNNNN' in the header '@Instrument:RunID:FlowCellID:Lane:Tile:X:Y

     1:N:0:NNNNNN'

  3. If '--illumina-file-names' is specified, the output files will be named according to the Illumina FASTQ file

     naming conventions:

 

a. The file extension will be '_R1_001.fastq.gz' for read one, and '_R2_001.fastq.gz' for read two (if paired end). b.

The per-sample output prefix will be '<SampleName>_S<SampleOrdinal>_L<LaneNumber>' (without angle brackets).

 

Options (1) and (2) require the input FASTQ read names to contain the following elements:

 

'@<instrument>:<run number>:<flowcell ID>:<lane>:<tile>:<x-pos>:<y-pos> <read>:<is filtered>:<control number>:<index>'

 

See the Illumina FASTQ conventions for more details.

(https://support.illumina.com/help/BaseSpace_OLH_009008/Content/Source/Informatics/BS/FASTQFiles_Intro_swBS.htm)

 

The '--illumina-standards' option may not be specified with the three options above. Use this option if you intend to

upload to Illumina BaseSpace. This option implies:

 

'--omit-fastq-read-numbers=true --include-sample-barcodes-in-fastq=false --illumina-file-names=true'

 

See the Illumina Basespace standards described here

(https://help.basespace.illumina.com/articles/tutorials/upload-data-using-web-uploader/).

 

To output with recent Illumina conventions (circa 2021) that match 'bcl2fastq' and 'BCLconvert', use:

 

'--omit-fastq-read-numbers=true --include-sample-barcodes-in-fastq=true --illumina-file-names=true'

 

By default all input reads are output. If your input FASTQs contain reads that do not pass filter (as defined by the

Y/N filter flag in the FASTQ comment) these can be filtered out during demultiplexing using the '--omit-failing-reads'

option.

 

To output only reads that are not control reads, as encoded in the '<control number>' field in the header comment, use

the '--omit-control-reads' flag

 

DemuxFastqs Arguments:

------------------------------------------------------------------------------------------------------------------------

-i PathToFastq+, --inputs=PathToFastq+

                              One or more input fastq files each corresponding to a sub-read (ex. index-read, read-one,

                              read-two, fragment). 

-o DirPath, --output=DirPath  The output directory in which to place sample BAMs. 

-x FilePath, --metadata=FilePath

                              A file containing the metadata about the samples. 

-r ReadStructure+, --read-structures=ReadStructure+

                              The read structure for each of the FASTQs. 

-h true|false, --help=true|false

                              Display the help message. [Default: false]. 

--version=true|false      Display the version number for this tool. [Default: false]. 

-m FilePath, --metrics=FilePath

                              The file to which per-barcode metrics are written. If none given, a file named

                              'demux_barcode_metrics.txt' will be written to the output directory. [Optional].

                              

-c String, --column-for-sample-barcode=String

                              The column name in the sample sheet or metadata CSV for the sample barcode.

                              [Default: Sample_Barcode]. 

-u String, --unmatched=String Output BAM file name for the unmatched records. [Default: unmatched.bam].

                              

-q QualityEncoding, --quality-format=QualityEncoding

                              A value describing how the quality values are encoded in the FASTQ. Either Solexa for

                              pre-pipeline 1.3 style scores (solexa scaling + 66), Illumina for pipeline 1.3 and above

                              (phred scaling + 64) or Standard for phred scaled scores with a character shift of 33. If

                              this value is not specified, the quality format will be detected automatically.

                              [Optional]. Options: Solexa, Illumina, Standard.

-t Int, --threads=Int         The number of threads to use while de-multiplexing. The performance does not increase

                              linearly with the # of threads and seems not to improve beyond 2-4 threads.

                              [Default: 1]. 

--max-mismatches=Int          Maximum mismatches for a barcode to be considered a match. [Default: 1].

                              

--min-mismatch-delta=Int      Minimum difference between number of mismatches in the best and second best barcodes for

                              a barcode to be considered a match. [Default: 2]. 

--max-no-calls=Int            Maximum allowable number of no-calls in a barcode read before it is considered

                              unmatchable. [Default: 2]. 

--sort-order=SortOrder        The sort order for the output sam/bam file (typically unsorted or queryname).

                              [Default: queryname]. Options: unsorted, queryname, coordinate, duplicate,

                              unknown.

--umi-tag=String              The SAM tag for any molecular barcode. If multiple molecular barcodes are specified, they

                              will be concatenated and stored here. [Default: RX]. 

--platform-unit=String        The platform unit (typically '<flowcell-barcode>-<sample-barcode>.<lane>')

                              [Optional]. 

--sequencing-center=String    The sequencing center from which the data originated [Optional]. 

--predicted-insert-size=Integer

                              Predicted median insert size, to insert into the read group header [Optional].

                              

--platform-model=String       Platform model to insert into the group header (ex. miseq, hiseq2500, hiseqX)

                              [Optional]. 

--platform=String             Platform to insert into the read group header of BAMs (e.g Illumina) [Default:

                              Illumina]. 

--comments=String*            Comment(s) to include in the merged output file's header. [Optional]. 

--run-date=Iso8601Date        Date the run was produced, to insert into the read group header [Optional].

                              

--output-type=OutputType      The type of outputs to produce. [Optional]. Options: Fastq, Bam,

                              BamAndFastq.

--include-all-bases-in-fastqs=true|false

                              Output all bases (i.e. all sample barcode, molecular barcode, skipped, and template

                              bases) for every read with template bases (ex. read one and read two) as defined by the

                              corresponding read structure(s). [Default: false]. 

--illumina-standards=true|false

                              Output FASTQs according to Illumina BaseSpace Sequence Hub naming standards. This is

                              differfent than Illumina naming standards. [Default: false].  Cannot be

                              used in conjunction with argument(s): includeSampleBarcodesInFastq, omitFastqReadNumbers,

                              illuminaFileNames

--omit-fastq-read-numbers=true|false

                              Do not include trailing /1 or /2 for R1 and R2 in the FASTQ read name. [Default:

                              false].  Cannot be used in conjunction with argument(s): illuminaStandards

--include-sample-barcodes-in-fastq=true|false

                              Insert the sample barcode into the FASTQ header. [Default: false]. 

                              Cannot be used in conjunction with argument(s): illuminaStandards

--illumina-file-names=true|false

                              Name the output files according to the Illumina file name standards. [Default:

                              false].  Cannot be used in conjunction with argument(s): illuminaStandards

--omit-failing-reads=true|false

                              Keep only passing filter reads if true, otherwise keep all reads. Passing filter reads

                              are determined from the comment in the FASTQ header. [Default: false]. 

--omit-control-reads=true|false

                              Do not keep reads identified as control if true, otherwise keep all reads. Control reads

                              are determined from the comment in the FASTQ header. [Default: false]. 

--mask-bases-below-quality=IntMask bases with a quality score below the specified threshold as Ns [Default:

                              0]. 

 

 

 実行方法

ランするにはindex情報を書いたCSV形式のサンプルシートが必要。また、Read Structure情報も必要。Read Structureは、シーケンスランの塩基をどのように割り当てるかを記述したStringで、<number><operator>のペアで表す。オプションとして、文字列の最後のセグメントには、その長さに番号の代わりに+を使用することが許可されている。この+は、他のセグメントが処理された後に残ったすべての塩基に変換される([0..infinity]を意味する)。

https://github.com/fulcrumgenomics/fgbio/wiki/Read-Structures

4種類のoperatorがサポートされている。

  • TまたはTemplate:セグメント内の塩基は、テンプレート(ゲノムDNA、RNAなど)の配列
  • B or Sample Barcode:セグメント中の塩基は、配列決定中のサンプルを識別するためのインデックス配列
  • M or Molecular Barcode: セグメント内の塩基は、配列中のユニークなソース分子(例:UMI)を識別するためのインデックス配列
  • S or Skip:セグメント内の塩基をスキップまたは無視する。例えば、ライブラリ調製で生成されたモノテンプレート配列

一般的なルール

  • 各セグメントは、正の整数 >= 1 (または +) でなければならない。
  • 読み込んだStructureの最後のセグメントのみ、その長さに + を使用することができる。
  • 隣接するセグメントには、同じ演算子を使用することができる。例えば、2つのサンプルインデックスが隣接するように別々に分子上にライゲーションされる場合、6B6B+TのStructureは許容される。

2つの異なる方法でRead Structureを記述する4つの例。


SampleSheet

4つ、もしくは5つの列からなる。

1,Sample_ID:サンプルシート内のサンプルに固有のサンプル識別子。

2,Sample_Name:サンプル名。

3,Library_ID:ライブラリ識別子。サンプル名とライブラリ識別子の組み合わせは、サンプルシート内のサンプル間でユニークである必要がある。

(任意)、Description:出力されるBAMのリードグループの説明フィールドに配置されるサンプルの説明。この列は省略可能。
4,Sample_Barcode:各サンプルに固有のサンプルバーコード塩基。サンプルバーコードを含むカラム名は、--column-for-sample-barcodeオプションで変更可能。サンプルバーコードが複数のリードにまたがって存在する場合(例:デュアルインデックス、ペアの両方のリードにインライン)、各リードのバーコード塩基は、このツールに与えられたリードのFASTQおよびリード構造の順序と同じ順序で連結される必要がある。

 

-iでdemultiplexingを行うfastq、--read-structuresでリード構造(リンク)、-xでサンプルシートのCSVファイル、-tでスレッド数、-oで出力ディレクトリ(存在しない時は作成される)を指定する。

fgbio DemuxFastqs -i r1.fq i1.fq i2.fq r2.fq --read-structures 8B92T 8B 8B 100T \
-x SampleSheet.csv -o outdir
  • -i     One or more input fastq files each corresponding to a sub-read (ex. index-read, read-one, read-two, fragment).
  • -o    The output directory in which to place sample BAMs.
  • -x     A file containing the metadata about the samples. 
  • -m    The file to which per-barcode metrics are written. If none given, a file named  'demux_barcode_metrics.txt' will be written to the output directory. [Optional].
  • --compression=Int     Default GZIP compression level, BAM compression level. [Default: 5]. 

  • --tmp-dir=DirPath     Directory to use for temporary files.

  • --max-mismatches=Int     Maximum mismatches for a barcode to be considered a match. [Default: 1].

  • --min-mismatch-delta=Int    Minimum difference between number of mismatches in the best and second best barcodes for a barcode to be considered a match. [Default: 2]. 

  • --max-no-calls=Int     Maximum allowable number of no-calls in a barcode read before it is considered unmatchable. [Default: 2]. 

                                  

出力ディレクトリに、demultiplexing されたサンプルごとのfastqファイル、demultiplexing されなかったfastq、要約統計のテキストファイルが出力される。

 

引用

A Universal Analysis Pipeline for Hybrid Capture-Based Targeted Sequencing Data with Unique Molecular Indexes
Min-Jung Kim, Si-Cho Kim, and Young-Joon Kim

Genomics Inform. 2018 Dec; 16(4): e29