macでインフォマティクス

macでインフォマティクス

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

RNA seqのバリアントコールにも対応したABRA2

 

 次世代シーケンス(NGS)は、さまざまなアプリケーションで広く使用されるツールになっている。バリアントコールは大きな関心が寄せられている領域であり、RNAへの関心も高まっている。NGSバリアントコールパイプラインの最初のステップの1つは、シーケンスされたリードをリファレンスゲノムにアラインすることである。 bwa-memやbowtie2などの広く使用されているDNAアライナー(Langmead and Salzberg、2012; Li、2013)は、多数のリードを迅速にアラインメントし、ギャップアラインメントをサポートし、indelの識別を可能にする。速度を上げるために、これらの方法はそれぞれのリードを徹底的には調整せず、場合によっては、特に indelサイズ増加によりindelを明らかにできない場合がある。従来、バリアントコーラーは、正確にマッピングされたリードに基づいて、SNVとindelを識別する。リードが正確にマッピングされていないインスタンスでは、バリアントの検出を混乱させる可能性がある。

 近年、DNAのidnel検出を改善する多くの方法が開発された。場合によっては、バリアントを含むリードがおおよそ正しい場所にマッピングされていれば、バリアントを正常に識別することができる。元のABRA実装(Mose et al、2014)では、localized assemblyプロセスを使用してリードアラインメントを調整し、その結果、アラインメントのindelを明らかにし、Freebayes(Garrison and Marth、2012)、およびsomaticコーラーのStrelka(Saunders et al、2012)などのバリアントコーラーの検出を改善した。Platypus、GATK-HaplotypeCaller(GATK-HC)、Strelka2、Scalpel(DePristo et al、2011; Kim et al、2018; Narzisi et al、2014; Rimmer et al、2014)はすべてlocalized assemblyを使用して、元のリードアラインメントに含まれる場合と含まれない場合があるインデルの検出を可能にする。同様に、最近開発されたLancetおよびMutect2(Cibulskis et al、2013; Narzisi et al、2018)は、体細胞変異コールにlocalized assemblyを使用する。

 RNA-Seqは、遺伝子およびアイソフォームの発現、遺伝子融合、転写スプライシング、発現変異、RNA編集の分析を可能にする診断ツールとして非常に重要であることが証明されている。 RNAスプライスジャンクションの存在は、スプライス対応アライナーの開発を必要とした。スプライスジャンクションにまたがるリードをマッピングできるいくつかのRNA-Seqアライナーが開発された(Dobin et al、2013; Kim et al、2013、2015; Wang et al、2010; Wu et al、2016)。ただし、些細ではないバリエーションが存在すると、RNA-Seqリードが完全にアラインメントされないことがある。例えばSunらは、標準のRNA-Seqパイプラインでは、2塩基を超える長さのindelを特定するのが困難であることを実証している (Sun, 2016)。多くの場合、元々DNA用に開発されたバリアントコールパイプラインは、RNA-Seqアラインメントの構文を処理するように変更されているが、RNA-Seq用には最適化されていない。たとえば、広く使用されているGATK(DePristo et al。、2011)では、スプライスジャンクションを削除してRNA-Seqリードアラインメントを複数のアラインメントに分割し、DNAのようなリードを生成して、DNAバリアント用に開発されたツールを使用して処理する必要がある。さらに、最近開発されたTransindel(Yang et al。、2018)は、リードの初期アラインメントをbwa-mem(DNAアライナー)に依存している。

 ABRA2は、RNA-Seqデータのスプライス対応の再調整と、大幅に強化された計算パフォーマンスを提供する、元のABRA実装の更新である。 ABRA2によって生成される改善されたアライメントにより、一般的に、特にindelに対して、より正確なバリアントコールが可能になる。さらに、ABRA2は元のABRAの精度を改善し、実行時間を短縮し、ヒト全ゲノムを処理できる拡張スケーラビリティを可能にする。

 簡単に言えば、元のABRA実装は、領域ごとに入力BAMファイルをローカライズして処理する。各関心の領域のリード値は、deBruijnグラフを使用して抽出およびアセンブリされる。アセンブルされたコンティグは、グローバルのコンティグプールに追加される。すべての領域でアセンブリが完了すると、bwa-memを使用してすべてのコンティグがリファレンスゲノムにアラインされる。特定のコンティグのキメラアラインメントは、アラインメントが単一のバリアントとして単純に表現できるindelの存在を明確に示している場合に結合される。すべてのコンティグがアラインメントされると、bwa-memを使用してすべてのリードがすべてのコンティグにマップされる。リードがリファレンスよりもコンティグに近い場合、リファレンスコンテキストがコンティグのアラインメントに一致するようにアップデートされる。

この方法は、遺伝子パネルやwhole exomeなど、多くのターゲットDNAシーケンスで効果的であることが証明されているが、いくつかの明らかな欠点がある。コンティグの総数が大きくなると、すべてのリードをすべてのコンティグにアラインすることが非常に遅くなり、場合によってはbwa-memが最後まで実行されないことがある。これは、多くのマウス系統の場合のように、ノイズの多いサンプルおよびリファレンスゲノムから実質的に分岐するサンプルの計算の困難さと同様に、スケーラビリティの問題を引き起こす。スケーラビリティの問題のため、ヒトの全ゲノムなどの大きなサンプルは通常、最後まで実行できない。キメラbwa-memアラインメントの後処理は、単一の単純な孤立したindel存在下ではうまく機能する可能性があるが、複数のバリアントまたは近接したテクニカルアーティファクトからなるより複雑な、またはノイズの多いイベントを適切にキャプチャできない場合がある。localized assemblyは、バリアントコールで広く使用されるようになったが、計算コストが高く、適切なバリアントを識別するために必ずしも必要ではないため、計算時間が不必要に長くなる。最後に、元のABRAもbwa-memもスプライシングを考慮していないため、ツールはRNA-Seqデータに対して最適ではない。これらの問題を改善するためにABRA2を開発した。

 

 

 

インストール

ubuntu18.04LTSでテストした。

Github

#bioconda (link)
conda install -c bioconda -y abra2

abra2

# abra2 

Picked up JAVA_TOOL_OPTIONS: -Xmx4G

INFO Thu Dec 12 22:51:31 JST 2019 Abra version: 2.22

INFO Thu Dec 12 22:51:31 JST 2019 Abra params: [/root/anaconda3/share/abra2-2.22-0/abra2.jar]

Missing required input SAM/BAM file(s)

Missing required output SAM/BAM filename(s)

Missing required reference

Option                                  Description                            

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

--amq <Integer>                         Set mapq for alignments that map       

                                          equally well to reference and an     

                                          ABRA generated contig.  default of   

                                          -1 disables (default: -1)            

--ca                                    Contig anchor [M_bases_at_contig_edge, 

                                          max_mismatches_near_edge] (default:  

                                          10,2)                                

--cl <Integer>                          Compression level of output bam file   

                                          (s) (default: 5)                     

--cons                                  Use positional consensus sequence when 

                                          aligning high quality soft clipping  

--contigs                               Optional file to which assembled       

                                          contigs are written                  

--dist <Integer>                        Max read move distance (default: 1000) 

--gc                                    If specified, only reprocess regions   

                                          that contain at least one contig     

                                          containing an indel or splice        

                                          (experimental)                       

--gkl                                   If specified, use the GKL Intel        

                                          Deflater.                            

--gtf                                   GTF file defining exons and transcripts

--in                                    Required list of input sam or bam file 

                                          (s) separated by comma               

--in-vcf                                VCF containing known (or suspected)    

                                          variant sites.  Very large files     

                                          should be avoided.                   

--index                                 Enable BAM index generation when       

                                          outputting sorted alignments (may    

                                          require additonal memory)            

--junctions                             Splice junctions definition file       

--keep-tmp                              Do not delete the temporary directory  

--kmer                                  Optional assembly kmer size(delimit    

                                          with commas if multiple sizes        

                                          specified)                           

--log                                   Logging level (trace,debug,info,warn,  

                                          error) (default: info)               

--mac <Integer>                         Max assembled contigs (default: 64)    

--mad <Integer>                         Regions with average depth exceeding   

                                          this value will be downsampled       

                                          (default: 1000)                      

--mapq [Integer]                        Minimum mapping quality for a read to  

                                          be used in assembly and be eligible  

                                          for realignment (default: 20)        

--maxn [Integer]                        Maximum pre-pruned nodes in regional   

                                          assembly (default: 150000)           

--mbq [Integer]                         Minimum base quality for inclusion in  

                                          assembly.  This value is compared    

                                          against the sum of base qualities    

                                          per kmer position (default: 20)      

--mcl [Integer]                         Assembly minimum contig length         

                                          (default: -1)                        

--mcr <Integer>                         Max number of cached reads per sample  

                                          per thread (default: 1000000)        

--mer <Double>                          Min edge pruning ratio.  Default value 

                                          is appropriate for relatively        

                                          sensitive somatic cases.  May be     

                                          increased for improved speed in      

                                          germline only cases. (default: 0.01) 

--mmr <Double>                          Max allowed mismatch rate when mapping 

                                          reads back to contigs (default: 0.05)

--mnf <Integer>                         Assembly minimum node frequency        

                                          (default: 1)                         

--mrn <Double>                          Reads with noise score exceeding this  

                                          value are not remapped.              

                                          numMismatches+(numIndels*2) <        

                                          readLength*mnr (default: 0.1)        

--mrr <Integer>                         Regions containing more reads than     

                                          this value are not processed.  Use   

                                          -1 to disable. (default: 1000000)    

--msr <Integer>                         Max reads to keep in memory per sample 

                                          during the sort phase.  When this    

                                          value is exceeded, sort spills to    

                                          disk (default: 1000000)              

--no-edge-ci                            If specified, do not update alignments 

                                          for reads that have a complex indel  

                                          at the read edge.  i.e. Do not allow 

                                          alignments like: 90M10D10I           

--no-ndn                                If specified, do not allow adjacent N- 

                                          D-N cigar elements                   

--nosort                                Do not attempt to sort final output    

--out                                   Required list of output sam or bam file

                                          (s) separated by comma               

--rcf <Double>                          Minimum read candidate fraction for    

                                          triggering assembly (default: 0.01)  

--ref                                   Genome reference location              

--sa                                    Skip assembly                          

--sc                                    Soft clip contig args [max_contigs,    

                                          min_base_qual,frac_high_qual_bases,  

                                          min_soft_clip_len] (default:         

                                          16,13,80,15)                         

--sga                                   Scoring used for contig alignments     

                                          (match, mismatch_penalty,            

                                          gap_open_penalty,                    

                                          gap_extend_penalty) (default:        

                                          8,32,48,1)                           

--single                                Input is single end                    

--skip                                  If no target specified, skip           

                                          realignment of chromosomes matching  

                                          specified regex.  Skipped reads are  

                                          output without modification.         

                                          Specify none to disable. (default:   

                                          GL.*|hs37d5|chr.*random|chrUn.       

                                          *|chrEBV|CMV|HBV|HCV.*|HIV.          

                                          *|KSHV|HTLV.*|MCV|SV40|HPV.*)        

--sobs                                  Do not use observed indels in original 

                                          alignments to generate contigs       

--ssc                                   Skip usage of soft clipped sequences   

                                          as putative contigs                  

--sua                                   Do not use unmapped reads anchored by  

                                          mate to trigger assembly.  These     

                                          reads are still eligible to          

                                          contribute to assembly               

--target-kmers                          BED-like file containing target        

                                          regions with per region kmer sizes   

                                          in 4th column                        

--targets                               BED file containing target regions     

--threads <Integer>                     Number of threads (default: 4)         

--tmpdir                                Set the temp directory (overrides java.

                                          io.tmpdir)                           

--ujac                                  If specified, use junction permuations 

                                          as contigs (Experimental - may use   

                                          excessive memory and compute times)  

--undup                                 Unset duplicate flag                   

--ws                                    Processing window size and overlap     

                                          (size,overlap) (default: 400,200)    

 

 

実行方法

DNA

coordinate sortされたbam、出力のbamファイル名、リファレンスゲノムを指定する。作業ディレクトリとして/tmpなどを指定する(少なくともbamと同じ程度のディスクスペースが必要)。ターゲットのbedを指定しない場合、全ゲノム領域がリアラインメント対象になる。

abra2 --in normal.bam,tumor.bam --out normal.abra.bam,tumor.abra.bam --ref hg38.fa --threads 8 --targets targets.bed --tmpdir /your/tmpdir > abra.log
  • --in-vcf   VCF containing known (or suspected)  variant sites. Very large files  (e.g. dbSNP)should be avoided.

 

RNA

starでmappingして得たbamを指定する。

abra2 --in star.bam --out star.abra.bam --ref hg38.fa --junctions bam --threads 8 --gtf gencode.v26.annotation.gtf --dist 500000 --sua --tmpdir /your/tmpdir > abra2.log 2>&1
  • --junctions    Splice junctions definition file
  • --gtf    file defining exons and transcripts

 

引用

Improved indel detection in DNA and RNA via realignment with ABRA2
Lisle E Mose, Charles M Perou, Joel S Parker
Bioinformatics, Volume 35, Issue 17, 1 September 2019, Pages 2966–2973

 

関連


既知生合成遺伝子クラスターのデータベース MIBiG 2.0

 

 植物、微生物、菌類は、多くの場合、1つまたはいくつかの種にユニークな多種多様な特殊な代謝物を生成する。文明の夜明けから、人間は薬用、経済的、またはレクリエーション目的でこの宝の山を利用している。過去10年以内に、ゲノムに基づいた特殊な代謝産物の発見は、科学界と商業環境の両方で広く採用されてきた。これらの努力の規模は、公共データベースでのゲノムおよびメタゲノムアセンブリの可用性の継続的な増加のため、継続的に成長している。これらの配列は、1つまたは複数の特定の化合物の生合成経路をコードするマルチ酵素遺伝子座の生合成遺伝子クラスター(BGC)の存在に基づいてマイニングできる。

 したがって、antiSMASH(ref.1)やClusterFinder(ref.2)などの計算ツールを使用して、数千のBGC候補が特定されている。 IMG-ABC(ref.3)やantiSMASH-DB(ref.4)などのデータベースには、非常に多様な範囲の天然物クラスをコーディングする可能性のある、数千ものそのような計算上予測されたBGCが格納されている。現在および将来の候補BGCの機能と新規性を解明するには、以前に特徴付けられたBGCに関する知識が不可欠である。これは、既知の化学構造の分子に関連するBGCの標準化された蓄積と抽出を必要とする。これは、関連する知識が通常、科学文献に埋もれているためである。

 このための最初のステップは2013年に行われ、ClusterMine360(ref.5)が登場した。これは、約300の遺伝子クラスターに関するデータを含む既知の製品を含むBGCの最初のデータベースである。 2015年に、MIBIG(生合成遺伝子クラスターに関する最小情報)データ標準とリポジトリが確立された。これには、コミュニティの努力により手動でキュレーションされた1170のBGCエントリが含まれ、その結果はかなり単純なWebアプリケーションからアクセスできる(ref.6) 。現在、MIBiGリポジトリは、既知の機能のBGCのセントラルリファレンスデータベースになり、KnownClusterBlastモジュールを介してantiSMASH(ref.1)の比較分析の基礎を提供する。微生物および微生物群集の小規模および大規模研究の中心となるBGC機能と新規性の多くの計算分析を可能にした。最近、MIBiGを使用して、研究されていない門からの376のメタゲノムアセンブルされた土壌細菌のBGCの並外れた新規性を評価および強調し、これらのBGCのほとんどがMIBiGからの遺伝子クラスターとの相同性を欠いていることを示した。同様に、Bahram et al (ref.8, pubmed)はMIBiGを使った相同性検索により、アノテーション付けできるMIBiG遺伝子クラスターのセットに基づいて、7560メタゲノムサンプル全体で抗菌活性に関連する真菌BGCを特定した。したがって、そのような「抗菌」BGCの豊富さは、土壌全体の抗菌薬耐性遺伝子の存在と相関していることを示した。さらに別の使用法がClusterCADツール(ref.9)に示されている。このツールは、新しい生化学パスウェイのコンピューター支援設計の出発点としてMIBiGからBGCデータを取得する。

 ここでは、更新されたMIBiGバージョン2.0を提供する。これは、過去5年間で851の新しいエントリが追加されて大幅に拡張された(論文図1)。さらに、データベース全体の広範な再アノテーション付けを行い、データスキーマを改善し、数百の文献と化学構造を追加し、最近出現した化学構造と分析データのデータベースへのクロスリンクを提供することにより、全体的なデータ品質を向上させた。最後に、化合物名、分類識別子、または生合成クラスに基づいた高速フィルタリングを有効にし、ブールクエリの構築を容易にすることで、便利な機能をオンラインリポジトリに追加し、より使いやすくした。

 


使い方

https://mibig.secondarymetabolites.orgにアクセスする。

f:id:kazumaxneo:20191210115145p:plain

 

Statistics

生産物質名、属名などで検索できる。

Secondary Metabolite record counts。

f:id:kazumaxneo:20191211234356p:plain

属はデータベース登録cluster数が多い順に並んでいる。ここでは636もアサインされているStreptomycin属をクリックした。

f:id:kazumaxneo:20191211230223p:plain

actinomycin Dをクリックした。

f:id:kazumaxneo:20191211230401p:plain

 

遺伝子クラスターが表示される。

f:id:kazumaxneo:20191211230453p:plain

 

右上からクラスターの遺伝子のGenBankファイルをダウンロードできる。

f:id:kazumaxneo:20191211231104p:plain

 

Compounds

f:id:kazumaxneo:20191211233137p:plain

Genes

f:id:kazumaxneo:20191211233158p:plain

History

f:id:kazumaxneo:20191211233226p:plain

Similar known gene clusters(一番上がこのデータベースのgene cluster(クエリ))

f:id:kazumaxneo:20191211233231p:plain



Search

Gene clusterをキーワード検索できる。

f:id:kazumaxneo:20191211233659p:plain


 

 

引用
MIBiG 2.0: a repository for biosynthetic gene clusters of known function
Kautsar SA, Blin K, Shaw S, Navarro-Muñoz JC, Terlouw BR, van der Hooft JJJ, van Santen JA, Tracanna V, Suarez Duran HG, Pascal Andreu V, Selem-Mojica N, Alanjary M1, Robinson SL, Lund G, Epstein SC, Sisto AC, Charkoudian LK, Collemare J, Linington RG, Weber T, Medema MH

Nucleic Acids Res. 2019 Oct 15. pii: gkz882. doi: 10.1093/nar/gkz882

 

関連


fastqの分析結果をJSON形式で出力する fastq-scan

 

タイトルの通りのツール。開発の動機はGithubのREADME参照。

 

インストール

macos10.14のPython 3.7環境でテストした。

 

#bioconda (link)
conda install -c bioconda -y fastq-scan

fastq-scan -h

$ fastq-scan -h

Usage: cat FASTQ | fastq-scan [options]

Version: 0.4.1

 

Optional arguments:

    -g INT   Genome size for calculating estimated sequencing coverage. (Default 1)

    -p INT   ASCII offset for input quality scores, can be 33 or 64. (Default 33)

    -v       Print version information and exit

    -h       Show this message and exit

 

実行方法

fastqを指定する。ゲノムサイズも指定すればカバレッジ情報も出力される。ゲノムサイズ5Mbなら

cat input.fq | fastq-scan -g 5000000 > out.json

#gzipped fatq
zcat input.fq.gz | fastq-scan 5000000 > out.json

 

ターミナルで確認する。ここではJavaScriptの処理ツールfxで確認する。 

#fxのインストール(github)
npm install -g fx

#viewする
fx out.json

f:id:kazumaxneo:20191209195857p:plain

=>で展開。

f:id:kazumaxneo:20191209200030p:plain

 

オンラインビューアやツールも色々利用できる。

JSON Editor

https://jsoneditoronline.org

引用

https://github.com/rpetit3/fastq-scan

 

 

miniasmでアセンブリして得たGFAをポリッシュする minipolish

 

Miniasmはパワフルで高速なロングリードのアセンブリツールだが、polishステップを持たないため、実質、得られた配列は連結されたロングリードである。polishにはraconが使用できるが、raconはFASTAファイルで動作し、Miniasmが出力するGFAをファイルを入力とし、グラフの情報を保持したままpolishすることはできない。minipolishはこのため開発されたツールで、このツールはraconを使ってGFAファイルをpolishする。また以下の様な情報も加えてくれる。

Githubの文章

  • Adding read depth information to contigs
  • Fixing sequence truncation that can occur in Racon
  • Adding circularising links to circular contigs if not already present (so they display better in Bandage)
  • 'Rotating' circular contigs between polishing rounds to ensure clean circularisation

 

インストール

依存

  • You'll need Python 3.6 or later to run Minipolish.
  • Minipolish assumes that you have minimap2 and Racon installed and available in your PATH.
  • You'll also need pytest if you want to run Minipolish's unit tests.

本体 Github

git clone https://github.com/rrwick/Minipolish.git
pip3 install ./Minipolish

> minipolish -h

$ minipolish -h

usage: minipolish [-t THREADS] [--rounds ROUNDS] [--pacbio] [--skip_initial] [-h] [--version] reads assembly

 

Minipolish

 

Positional arguments:

  reads                          Long reads for polishing (FASTA or FASTQ format)

  assembly                       Miniasm assembly to be polished (GFA format)

 

Settings:

  -t THREADS, --threads THREADS  Number of threads to use for alignment and polishing (default: 8)

  --rounds ROUNDS                Number of full Racon polishing rounds (default: 2)

  --pacbio                       Use this flag for PacBio reads to make Minipolish use the map-pb Minimap2 preset (default: assumes Nanopore reads and uses the map-ont

                                 preset)

  --skip_initial                 Skip the initial polishing round - appropriate if the input GFA does not have "a" lines (default: do the initial polishing round

 

Other:

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

  --version                      Show program's version number and exit

 

 

 

実行方法

アセンブリ後のGFAファイルを入力とする。

minimap2 -t 8 -x ava-ont long_reads.fastq.gz long_reads.fastq.gz > overlaps.paf
miniasm -f long_reads.fastq.gz overlaps.paf > assembly.gfa
minipolish -t 8 long_reads.fastq.gz assembly.gfa > polished.gfa

 

レポジトリに付属しているスクリプトを使えば、上の工程を1コマンドで実行できる。

スレッド数8で実行。

miniasm_and_minipolish.sh long_reads.fastq.gz 8 > polished.gfa

 

BandageにGFAを読み込めば変化を視覚的に確認できる。 ここではターミナルで叩いて実行する。

Bandage image polished.gfa polished.svg

 

引用

https://github.com/rrwick/Minipolish

 

関連


 

 

 

 

 

植物のRNA seqデータからvirus配列を検出する Kodoja

 

 ウイルス感染は、食物と燃料のために栽培される作物で特に重要な問題である。ウイルスは収量と品質の大きな損失を引き起こし、その結果、ウイルスは重要な経済的悪影響を及ぼす[ref.1]。英国では、ポテトウイルスYは年間3,000〜4,000万ポンドのジャガイモの収穫損失[ref.2]を引き起こし、アジアではイネに感染するウイルス(ライスグラッシースタントウイルスなど)は年間1億2,000万ドルの作物損失を引き起こす[ref.3]。これらの例は、高速で正確なウイルス検出方法の必要性を強調している。ウイルス感染の症状には黄変や発育阻害が含まれる場合があるが、多くの場合、症状は存在しないか、他の要因によって隠されている。場合によっては、植物ウイルスが相乗的に相互作用して、新しいまたはより深刻な病気の症状を引き起こす[ref.4]。 1つの例は、ラズベリーの果物がもろくなる複雑な病気で、これは2つのウイルス、 Raspberry bushy dwarf virusおよびRaspberry潜伏ウイルス[ref.5]の存在によって引き起こされる。
 作物が新しい地理的位置で栽培され、農業慣行が強化されると、新しいウイルスが確立され、既存のウイルスが宿主範囲を広げるリスクが高まる。したがって、植物ウイルスの診断は、将来の食料安全保障の観点から重要性が高まっている分野である。
ウイルスを検出するための標準的な分子技術には、逆転写(RT)PCRに基づく方法が含まれる。ただし、このような手法では既知のウイルス検出のみ可能になる。つまり、各テストは1つのウイルスまたはごく少数の関連ウイルスに固有である[ref.6]。さらに、ウイルスゲノムが進化し、テストが時間の経過とともに無効になる可能性があり、病気の診断が制限される。このような制限は、最近、複数のウイルスを仮説なしで同時に検出するための次世代シーケンス(NGS)メソッドの使用により克服された[ref.7]。植物ウイルスの大部分はRNAを遺伝物質として持ち、DNAゲノムを持つウイルスはRNA転写産物を生成する。さらに、真核生物のsiRNAはRNA干渉を介して抗ウイルス免疫を誘導し、この過程でウイルス由来のsiRNAが宿主に濃縮される[ref.8]。したがって、RNAおよびsmall RNA(sRNA)シーケンスの両方が、植物のウイルス検出に効果的な方法である。ただし、これは2つの重要な要素に依存している:(a)堅牢なRNA抽出およびエンリッチメントプロトコル、および(b)ウイルス同定のための高速で堅牢なバイオインフォマティクスツール。
 RNA抽出およびエンリッチメントプロトコルの範囲、およびバイオインフォマティクスのワークフローは、以前にヒトの臨床サンプル用に開発された(レビューについては[ref.9]を参照)。最近、このような研究により、48時間以内にウイルス性疾患の診断と実用的な臨床管理が行われた[ref.10]。この臨床研究で使用されたワークフローは、ウイルス診断ツールに必要な2つの主要要素で構成されていた。(a)宿主ヌクレオチド配列の特定と除去、(b)ウイルス配列の特定。ただし、ヒトゲノムには十分なアノテーションが付けられており(ホスト配列の簡単な削除が可能)、またヒトウイルスデータは配列データベースでより一般的のため、臨床サンプルでのウイルス検出は植物よりも簡単である。それに比べて、多くの作物のゲノムは不完全であるかアノテーションが不十分であり、植物ウイルス配列はデータベースでは十分に表現されていない。

 最近、NGSデータからウイルス検出に利用可能なバイオインフォマティクスツールとワークフローをレビューした[ref.9]。大部分はヒトのNGSデータに最適化されており、ウイルスの識別にはシーケンスIDを超えるものはほとんどなく、インストールや使用に多くの重要な計算知識が必要であると結論付けた。 Taxonmer [ref.11]とVirusDetect [ref.12]の2つのツールはWebサーバーとして利用でき、植物[ref.2]からのRNAシーケンシングデータの分析の可能性を提供している。ただし、レビューでは、公開された3つのツールが植物データでテストされたが、植物のウイルス検出に焦点を当てたプロジェクトでは使用されていないという事実が強調された。代わりに、このアプローチは一般に分析中により高い柔軟性を提供するため、ワークフローの外部でスタンドアロンマッピングおよびアセンブリアルゴリズムを使用した。
 ウイルス識別ワークフローは、(a)低品質のリードおよびアダプターシーケンスのトリミングを含む、rawデータファイルの品質管理対策の実施、(b)ホストシーケンスの識別、および(c)ウイルスシーケンスの識別が可能である必要がある。既知のウイルスの識別は、既存のウイルス配列のデータベースにマッピングすることで実行できるが、新しい株や新しいウイルスの識別には、ワークフローを超えた専門知識と追加の分析が必要である。
 公開されているウイルス検出ワークフローの多くは、アセンブリマッピングアルゴリズムを使用してウイルスシーケンスを識別する。ただし、アセンブリマッピングの両方が非常に計算集約的である可能性があるため、ワークフローは
大規模なデータセットでは実行時間が長くなる可能性がある。また、アセンブリおよびマッピング方法では、未アセンブリのリードが識別されないままになる。 RNA-seqデータセットでウイルスリードを識別するもう1つの方法は、k-merプロファイリングを使用することである。これはTaxonomer [ref.11]で正常に実装されている。 RNAおよびDNAシーケンスは長さkの複数の部分文字列に分割できる。このようにして、シーケンスをk-merプロファイルで表すことができ、これらのプロファイルをtaxaに割り当てして比較できる。 k-merプロファイルは、バイオインフォマティクスの類似検索の範囲で使用されている。メタゲノミクスでは配列間のアラインメントのない類似性分析が可能である[ref.13]。分類プロファイリングでは、ビニング法はk-merプロファイルを使用して配列をクラスター化し、ドラフトゲノム回復を可能にする。このような方法は、リファレンスゲノムを使用せずに、集団内のウイルスハプロタイプの同定にもうまく適用されている。
 ここで紹介するKodojaワークフローは、NGSデータセットを使用する幅広い研究者に適用できる独自の機能セットを組み合わせている。目的は、アセンブリマッピングの方法を超えたワークフローを開発することであった。これは、植物データセット用に特に最適化され、非生物情報学者がアクセスできる。 Kodojaは、混合された(植物とウイルス、細菌、真菌の両方のヌクレオチドを含む)RNA-seqデータからウイルスシーケンスを識別することができるワークフローである。 Kodojaは、(a)植物NGSデータに固有であり、(b)ヌクレオチドレベルでk-merプロファイリングを使用し、既存のツールKraken [ref.16] ]およびKaiju [ref.17]を使い、(c)はBioconda [ref.18]を介してローカルインストールして使用でき、(d)Galaxy Webベースの分析環境[ref.19]内のツールとしても使用できる。

 

f:id:kazumaxneo:20191207200849p:plain

Kodoja workflow.  論文より転載。

 

インストール

ubuntu18.04のPython 3.7環境でテストした(docker使用、ホストOS macos10.14)。

依存

You can use Python 2.7 or Python 3, specifically Kodoja has been tested on Python 3.6.

  • FastQC v0.11.5
  • Trimmomatic v0.36
  • Kraken v1.0
  • Kaiju v1.5.0

Python packages:

  • numpy v1.9
  • biopython v1.67
  • pandas v0.14
  • ncbi-genome-download v0.2.6

本体 Github

#bioconda (link)
conda install -c bioconda kodoja

kodoja_search.py  -h

# kodoja_search.py  -h

usage: kodoja_search.py [-h] [--version] -o OUTPUT_DIR -d1 KRAKEN_DB -d2

                        KAIJU_DB -r1 READ1 [-r2 READ2] [-f DATA_FORMAT]

                        [-t THREADS] [-s HOST_SUBSET] [-m TRIM_MINLEN]

                        [-a TRIM_ADAPT] [-q KRAKEN_QUICK] [-p]

                        [-c KAIJU_SCORE] [-l KAIJU_MINLEN] [-i KAIJU_MISMATCH]

 

Kodoja Search is a tool intended to identify viral sequences

in a FASTQ/FASTA sequencing run by matching them against both Kraken and

Kaiju databases.

 

optional arguments:

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

  --version             show program's version number and exit

  -o OUTPUT_DIR, --output_dir OUTPUT_DIR

                        Output directory path, required

  -d1 KRAKEN_DB, --kraken_db KRAKEN_DB

                        Kraken database path, required

  -d2 KAIJU_DB, --kaiju_db KAIJU_DB

                        Kaiju database path, required

  -r1 READ1, --read1 READ1

                        Read 1 file path, required

  -r2 READ2, --read2 READ2

                        Read 2 file path

  -f DATA_FORMAT, --data_format DATA_FORMAT

                        Sequence data format (default fastq)

  -t THREADS, --threads THREADS

                        Number of threads (default 1)

  -s HOST_SUBSET, --host_subset HOST_SUBSET

                        Subset sequences with this tax id from results

  -m TRIM_MINLEN, --trim_minlen TRIM_MINLEN

                        Trimmomatic minimum length

  -a TRIM_ADAPT, --trim_adapt TRIM_ADAPT

                        Illumina adapter sequence file

  -q KRAKEN_QUICK, --kraken_quick KRAKEN_QUICK

                        Number of minium hits by Kraken

  -p, --kraken_preload  Kraken preload database

  -c KAIJU_SCORE, --kaiju_score KAIJU_SCORE

                        Kaju alignment score

  -l KAIJU_MINLEN, --kaiju_minlen KAIJU_MINLEN

                        Kaju minimum length

  -i KAIJU_MISMATCH, --kaiju_mismatch KAIJU_MISMATCH

                        Kaju allowed mismatches

 

The main output of ``kodoja_search.py`` is a file called ``virus_table.txt``

in the specified output directory. This is a plain text tab-separated table,

the columns are as follows:

 

1. Species name,

2. Species NCBI taxonomy identifier (TaxID),

3. Number of reads assigned by *either* Kraken or Kaiju to this species,

4. Number of Reads assigned by *both* Kraken and Kaiju to this species,

5. Genus name,

6. Number of reads assigned by *either* Kraken or Kaiju to this genus,

7. Number of reads assigned by *both* Kraken and Kaiju to this genus.

 

The output directory includes additional files, including ``kodoja_VRL.txt``

(a table listing the read identifiers used) which is intended mainly as

input to the ``kodoja_retrieve.py`` script.

 

See also https://github.com/abaizan/kodoja/wiki/Kodoja-Manual

kodoja_build.py -h

# kodoja_build.py -h

usage: kodoja_build.py [-h] [--version] -o OUTPUT_DIR [-t THREADS]

                       [-p HOST_TAXID] [-d DOWNLOAD_PARALLEL] [-n]

                       [-e [EXTRA_FILES [EXTRA_FILES ...]]]

                       [-x [EXTRA_TAXIDS [EXTRA_TAXIDS ...]]] [-v]

                       [-b KRAKEN_TAX] [-k KRAKEN_KMER] [-m KRAKEN_MINIMIZER]

                       [-a DB_TAG]

 

Kodoja Build creates a database for use with Kodoja Search.

 

optional arguments:

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

  --version             show program's version number and exit

  -o OUTPUT_DIR, --output_dir OUTPUT_DIR

                        Output directory path, required

  -t THREADS, --threads THREADS

                        Number of threads, default=1

  -p HOST_TAXID, --host_taxid HOST_TAXID

                        Host tax ID

  -d DOWNLOAD_PARALLEL, --download_parallel DOWNLOAD_PARALLEL

                        Parallel genome download, default=4

  -n, --no_download     Genomes have already been downloaded

  -e [EXTRA_FILES [EXTRA_FILES ...]], --extra_files [EXTRA_FILES [EXTRA_FILES ...]]

                        List of extra files added to "extra" dir

  -x [EXTRA_TAXIDS [EXTRA_TAXIDS ...]], --extra_taxids [EXTRA_TAXIDS [EXTRA_TAXIDS ...]]

                        List of taxID of extra files

  -v, --all_viruses     Build databases with all viruses (not plant specific)

  -b KRAKEN_TAX, --kraken_tax KRAKEN_TAX

                        Path to taxonomy directory

  -k KRAKEN_KMER, --kraken_kmer KRAKEN_KMER

                        Kraken kmer size, default=31

  -m KRAKEN_MINIMIZER, --kraken_minimizer KRAKEN_MINIMIZER

                        Kraken minimizer size, default=15

  -a DB_TAG, --db_tag DB_TAG

                        Suffix for databases

 

See also https://github.com/abaizan/kodoja/wiki/Kodoja-Manual

kodoja_retrieve.py -h

# kodoja_retrieve.py -h

usage: kodoja_retrieve.py [-h] [--version] -o FILE_DIR -r1 READ1 [-r2 READ2]

                          [-f USER_FORMAT] [-t TAXID] [-g] [-s]

 

Kodoja Retrieve is used with the output of Kodoja Search to

give subsets of your input sequencing reads matching viruses.

 

optional arguments:

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

  --version             show program's version number and exit

  -o FILE_DIR, --file_dir FILE_DIR

                        Path to directory of kodoja_search results, required

  -r1 READ1, --read1 READ1

                        Read 1 file path, required

  -r2 READ2, --read2 READ2

                        Read 2 file path, default: False

  -f USER_FORMAT, --user_format USER_FORMAT

                        Sequence data format, default: fastq

  -t TAXID, --taxID TAXID

                        Virus tax ID for subsetting, default: All viral

                        sequences

  -g, --genus           Include sequences classified at genus

  -s, --stringent       Only subset sequences identified by both tools

 

The main output of ``kodoja_search.py`` is a file called ``virus_table.txt``

(a table summarising the potential viruses found), but the specified output

directory will also contain ``kodoja_VRL.txt`` (a table listing the read

identifiers). This second file is used as input to ``kodoja_retrieve.py``

along with the original sequencing read files.

 

A sub-directory ``subreads/`` will be created in the output directory,

which will include either FASTA or FASTQ files named as follows:

 

* ``subset_files/virus_all_sequences1.fasta`` FASTA output

* ``subset_files/virus_all_sequences1.fastq`` FASTQ output

 

And, for paired end datasets,

 

* ``subset_files/virus_all_sequences2.fasta`` FASTA output

* ``subset_files/virus_all_sequences2.fastq`` FASTQ output

 

However, if the ``-t 12345`` option is used rather than ``virus_all_...``

the files will be named ``virus_12345_...`` instead.

 

See also https://github.com/abaizan/kodoja/wiki/Kodoja-Manual

 

 

データベースの準備

ここではオーサーらがビルドしたデータベース(krakenのDB)を使用する。 

mkdir kodojaDB_v1.0
cd kodojaDB_v1.0
wget https://zenodo.org/record/1406071/files/kodojaDB_v1.0.tar.gz
tar -xvf kodojaDB_v1.0.tar.gz

解凍してできたkaijuDB/とkrakenDB/を使う。

 

実行方法

kodoja_search.py --kraken_db kodojaDB_v1.0/krakenDB --kaiju_db kodojaDB_v1.0/kaijuDB \
-r1 pair_1.fq -r2 pair_2.fq -t 8 -o outdir
  • -r1    path to the single-end or first paired-end file (required)
  • -r2    path to second paired-end file (default=False)
  • -f      specify the file-type for file1 ("fasta" or "fastq" - default='fastq')
  • -o     path to the results folder (required)
  • -t      number of threads on cluster (default=1)
  • --kraken_db   path to kraken database (required)
  • --kaiju_db       path to kaiju database, nodes.dmp and names.dmp files (required)

 ランが終わるとvirus_table.txtに結果がまとめられる。

 

引用

Kodoja: A workflow for virus detection in plants using k-mer analysis of RNA-sequencing data

Amanda Baizan-Edge, Peter Cock, Stuart MacFarlane, Wendy McGavin, Lesley Torrance, Susan Jones

J Gen Virol. 2019 Mar;100(3):533-542

 

 

Artemis ゲノムブラウザ

 

スループットシーケンス(HTS)テクノロジにより、多数のサンプルの低コストシーケンスが一般的になった。 ゲノムリシーケンシング、集団規模の変異検出、全トランスクリプトームシーケンシング、タンパク質結合核酸のゲノムワイド解析なども行われている。Artemisは、リファレンスゲノムとそれに対応するアノテーションのコンテキストで、さまざまなタイプのHTSデータセットの統合された視覚化と計算分析のためのツールである。

 

HP (Artemis、ACT、BamView、DNAPlotterの紹介)

http://sanger-pathogens.github.io/Artemis/Artemis/

poster

https://sanger-pathogens.github.io/Artemis/Artemis/artemis_genome_informatics_2010.pdf

manual

https://sanger-pathogens.github.io/Artemis/Artemis/artemis-manual.html

 

インストール

HPから、macwindowslinux版をダウンロードできる。また、githubのリリースからダウンロードすることもできる。ここではgithubからmac版をダウンロードした。

依存

All recent releases from v18.0.0 onwards require a minimum of Java 9 and ideally Java 11. This must be installed first.

Github

ダウンロードしたdmgファイルを実行、.appをApplicationにコピーする。

f:id:kazumaxneo:20191206201118p:plain

 

 

実行方法

Artemis (CHADO)を起動する。

f:id:kazumaxneo:20191206230927p:plain

 

ファイルを開く。File => Open

f:id:kazumaxneo:20191206231000p:plain

 

リファレンスゲノムを読み込む。EMBL、GenBnak、FASTAなどに対応している。ここではbamファイルも読み込みたいので、そのリファレンスであるFASTAファイルを選択した。

f:id:kazumaxneo:20191206234929p:plain

読み込まれた。

 

BAMを読み込む。File => Read BAM / CRAM / VCF

f:id:kazumaxneo:20191207000946p:plain

 

BAMファイルを指定する。

f:id:kazumaxneo:20191207001027p:plain

 

Add moreを押してVCFも同時に読み込メル。VCFはあらかじめ圧縮されてindexがつけられている必要がある。

bgzip file.vcf
tabix -p vcf file.vcf.gz

f:id:kazumaxneo:20191207001135p:plain

 

ここではbamファイルだけ読み込んだ。

f:id:kazumaxneo:20191207002039p:plain

 

 

 

reverse strandのstopコドンのラインを非表示にする。

f:id:kazumaxneo:20191207002957p:plain



非表示になった。ここでは、左端の>>をクリックして不要なパネルも隠している。

f:id:kazumaxneo:20191207003308p:plain



右端の上下のスクロールバーを動かすことで拡大縮小できる。大きく縮小するとカバレッジ変化がラインで表現される。

f:id:kazumaxneo:20191207003424p:plain

 

拡大すると塩基が見えるようになる。

f:id:kazumaxneo:20191207003618p:plain

 

GenBankなどフィーチャー情報があるファイルを読み込んだ時は、豊富なツールを使用でき、単なるブラウザというより統合DNA解析ツールに近い(たとえばCLC genomicsやgenious)感じになる。

f:id:kazumaxneo:20191207004257p:plain

 

フィーチャーを選択してblastnを実行してみた。数秒待つとデフォルトのブラウザに結果が表示された(NCBIのオンラインのblast結果)。詳細な使い方はマニュアルを参照して下さい。

引用

Artemis: an integrated platform for visualization and analysis of high-throughput sequence-based experimental data
Carver T, Harris SR, Berriman M, Parkhill J and McQuillan JA
Bioinformatics (Oxford, England) 2012;28;4;464-9

 

Artemis and ACT: viewing, annotating and comparing sequences stored in a relational database
Carver T, Berriman M, Tivey A, Patel C, Böhme U et al.
Bioinformatics (Oxford, England) 2008;24;23;2672-6

 

Viewing and annotating sequence data with Artemis
Berriman M and Rutherford K
Briefings in bioinformatics 2003;4;2;124-32

 

Artemis: sequence visualization and annotation
Rutherford K, Parkhill J, Crook J, Horsnell T, Rice P, Rajandream MA and Barrell B
Bioinformatics (Oxford, England) 2000;16;10;944-5

 

CRISPRsをゲノムから探すMinCED

 

MinCEDは、完全なゲノムまたはメタゲノムからアセンブルされたコンティグなどの環境データセットで、Clustered Regularly Interspaced Short Palindromic Repeats (CRISPRs) を見つけるプログラムである。 

 

インストール

依存

Github

#bioconda (link) 
conda install -c bioconda minced

minced

$ minced 

MinCED, a program to find CRISPRs in shotgun DNA sequences or full genomes

 

Usage:    minced [options] file.fa [outputFile.txt] [outputFile.gff]

 

Options:  -searchWL  Length of search window used to discover CRISPRs (range: 6-9). Default: 8

          -minNR     Minimum number of repeats a CRISPR must contain. Default: 3

          -minRL     Minimum length of the CRISPR repeats. Default: 23

          -maxRL     Maximum length of the CRISPR repeats. Default: 47

          -minSL     Minimum length of the CRISPR spacers. Default: 26

          -maxSL     Maximum length of the CRISPR spacers. Default: 50

          -gff       Output summary results in gff format containing

                     only the positions of the CRISPR arrays. Default: false

          -gffFull   Output detailed results in gff format containing

                     positions of CRISPR arrays and all repeat units. Default: false

          -spacers   Output a fasta formatted file containing the spacers. Default: false

          -h --help  Output this handy help message

          --version  Output version information

 

Examples: minced ecoli.fna

          minced metagenome.fna

          minced metagenome.fna metagenome.crisprs

          minced metagenome.fna metagenome.crisprs metagenome.gff

 

 

 

実行方法

genomeのFASTAファイルを指定する。メタゲノム由来contigも使用できる。出力のgtfはなくても良い。

minced genome.fna out.txt out.gff

出力

f:id:kazumaxneo:20191205110307p:plain

 

引用

https://github.com/ctSkennerton/minced