macでインフォマティクス

macでインフォマティクス

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

DNA解析ソフト1 Snapgene

2018 10/26 追記

2019 4/22 プライマーの説明追加

2019 5/30 誤字修正

2019 7/1 追記

2020 1/6 リンク追加

2020 1/25 追記

2020 2/12 分かりにくい表現を修正

 

SnapGeneはプラスミドやゲノムの解析ソフトである。ビジュアルで見やすく表示する機能に特化しており、制限酵素処理結果を画面で見やすく表示したり、TA cloningやgibson、in fuisionなどの反応結果を表示し設計を支援する機能などを持つ。

 

 

 

SnapGeneはアノテーション情報がない遺伝子配列からorfを自動検出してプラスミド図などを描画する機能も持つが、この機能は無償版のSnapGene viewerでも使うことができる。さらに素晴らしいのは、無償版、機能の豊富な有償版ともLinux版、mac版、windows版全て揃っていることである。とりあえずゲノム配列を手に入れたらsnapgene (viewer)で確認してみよう、くらいの気軽さで使うことができる。

 

公式サイトからダウンロードできる。

SnapGene | Software for everyday molecular biology

インストールは画面の指示に従っていくだけでできる。

 

 

SnapGeneは、genbankファイルの取り込みに対応している。genbankファイルは、NCBIで遺伝子を見るときのフォーマットである。以下のリンクを開いてgenbankファイルを1つダウンロードしてみよう。

https://www.ncbi.nlm.nih.gov/nuccore/?term=cloning+vector

f:id:kazumaxneo:20170627130533j:plain

 

 

左端のgenbankを選択。

f:id:kazumaxneo:20170627130705j:plain

右上のSend toからgenbankを選択してダウンロードする。

f:id:kazumaxneo:20171001152156j:plain

 

 

または配列まで含めて画面全体をコピーし、

f:id:kazumaxneo:20170627130652j:plain

 ↓

それを軽量なプレーンテキストエディタmacならMi(Mi ダウンロードリンク)、windowsならTerapadサクラエディタ)にペーストする。

f:id:kazumaxneo:20170627130838j:plain

 

.gbkで保存。

f:id:kazumaxneo:20170627131016j:plain

 

Snapgeneを指定してgbkファイルを開く。

f:id:kazumaxneo:20170627131108j:plain

 

環状かどうか聞いてくるので、環状を選択。

f:id:kazumaxneo:20170627131212j:plain

 

開いた直後の画面。orfがgbkファイルの情報に従い自動認識されている。また制限酵素サイトが自動で表示されている。

f:id:kazumaxneo:20170627131242j:plain

 

左下のタブ ”シーケンス” を押すと配列に切り替わる。

f:id:kazumaxneo:20170627131711j:plain

 

左下のタブ ”酵素” を押すと制限酵素の位置が表示される。

f:id:kazumaxneo:20170627131900j:plain

 

上のタブ”ライン"を押すと切れる場所が縦線で表示される。cut部位の相対的な位置関係を掴みやすい。

f:id:kazumaxneo:20170627131941j:plain

 

追記 プライマー追加

f:id:kazumaxneo:20190423201522p:plain

ウィンドウが出てくるのでプライマー配列をペーストする。

f:id:kazumaxneo:20190423201529p:plain

マッチする部位が表示されるので、右下のテンプレートにプライマーを追加ボタンをクリック。

 

プライマー部位が追加された。

f:id:kazumaxneo:20190423201539p:plain

繰り返すことで複数登録できる。

 

2019 7/1 追記

サンガー配列の読み込みはviewerでもできます。

f:id:kazumaxneo:20190701145358j:plain

開いたところ。

f:id:kazumaxneo:20190701145410j:plain

 

マルチプルシーケンスアラインメント(MSA)結果も表示できる。

f:id:kazumaxneo:20200125001046p:plain

 

他にも様々な機能を持ち、有償版は、例えば配列の編集、配列同士のアライメントなどに対応している。

 

いくつか動画を貼っておく。

コンストラクトの構築手順。

動画後半では、右下のタブ”履歴” 機能を使いコンストラクトの設計を行っている。

 

アライメント。

SnapGene Software Tutorial Videos | Aligning to a Reference DNA Sequence

下のCCからsubtitleをonにすると説明がつきます。

 

Primer設計と変異導入。


 

アカデミックライセンスなら、1ユーザ345米ドルで購入でき、他の解析ソフトと比べても良心的な金額になっている。まずは無償版で試し、有償版の機能が欲しくなったら、(ボスにお願いして)購入を考えればよい。

有償版と無償版の機能の違い。

SnapGene vs. SnapGene Viewer | Molecular Biology Software

 

 他にもmac vectorも良い評判を聞いたことがある(安くて、かつ解析ソフトの機能は概ね備えている)。使っている人を見たことがあるが、genetyxに似ており概ね使いやすそうだった。2つ動画を貼っておく。

MacVector Tutorial

 

Create constructs using MacVector's Cloning Clipboard

 

 

これ以外にも無料で使えるツールはいくつかあります。

これらとsnap gene viewerと組み合われば、コストをかけずに使うことも可能です。

 

2020 1/6追記

UGENEおすすめです。

 

 

 引用

SnapGene software (from GSL Biotech; available at snapgene.com)

 

ゲノムの相同性の高い領域の網羅的な検索 MUMmer

 

MUMmerはゲノム全体を高速にアライメントするオープンソースのツール。Finisihしたゲノムだけでなくドラフトゲノムでも使用でき、容易に何百あるcontigのアライメントを行うことができる。最初の論文が発表されたのは1999年だが(ref.1)、現在でもオープンソースのまま開発が続けられており、 間も無くMUMmer4が論文化されるとアナウンスされている(その後、2018年始めにpublishされた)。

  

オンラインマニュアル

The MUMmer 3 manual

 

MUMmerは5つのアライメントプログラム(tu-rukaranaru.mummer、nucmer、promer、rum-mummer1、そしてrun-mummer3) が専ら使われている。

brewで導入できる。

brew install mummer

 

作成まで時間がかかってしまったので、別に投稿しました。 


 

 

引用------------------------------------------------------------------------------------------------------------------------

Mummer1

1、Alignment of whole genomes

Nucleic Acids Research, 1999, Vol. 27, No. 11 2369–2376

http://mummer.sourceforge.net/MUMmer.pdf

 

Mummer2

2、Fast algorithms for large-scale genome alignment and comparison.

Delcher AL1, Phillippy A, Carlton J, Salzberg SL.

Nucleic Acids Res. 2002 Jun 1;30(11):2478-83.

http://mummer.sourceforge.net/MUMmer2.pdf

 

Mummer3

3、Versatile and open software for comparing large genomes

Stefan Kurtz*, Adam Phillippy†, Arthur L Delcher†, Michael Smoot†‡, Martin Shumway†, Corina Antonescu† and Steven L Salzberg†

Genome Biology 2004, 5:R12

http://mummer.sourceforge.net/MUMmer3.pdf

 

追記

Mummer4

MUMmer4: A fast and versatile genome alignment system

Guillaume Marçais, Arthur L. Delcher, Adam M. Phillippy, Rachel Coston, Steven L. Salzberg, Aleksey Zimin

Published online 2018 Jan 26.

 

 

 

 

 

ゲノム比較のmurasakiと結果を表示するGMV

2020 5/25 docker link追加、ヘルプ追加、分かりにくい文章の修正

 

murasakiは複数ゲノムの相同性ある領域の探索を高速に行うツールで、GMVはその比較結果を見るためのビューアソフトである。領域によってカラフルな色がつくので、ゲノムリアレンジメントなどの構造変化をわかりやすく示すことができる。

 

公式サイト

murasaki genomemurasaki genome

 

murasakiのインストールはこの方のブログを参考にしてください。


2020 5/26 追記

#dockerhub(link)
docker pull biocontainers/murasaki:v1.68.6-8b1-deb_cv1

 > docker run --rm -it  biocontainers/murasaki:v1.68.6-8b1-deb_cv1 murasaki -h

$ docker run --rm -it biocontainers/murasaki:v1.68.6-8b1-deb_cv1 murasaki -h

Usage: murasaki -p<patternSpec> [options] seq1 [seq2 [seq3 ... ]]

Options:

  --pattern|-p   = seed pattern (eg. 11101001010011011).

                    using the format [<w>:<l>] automatically generates a

                    random pattern of weight <w> and length <l>

  --directory|-d = output directory (default: output)

  --name|-n      = alignment name (default: test)

  --quickhash|-q = specify a hashing function:

                    0 - adaptive with S-boxes

                    1 - don't pack bits to make hash (use first word only)

                    2 - naively use the first hashbits worth of pattern

                    3 - adaptivevely find a good hash (default)

                    **experimental CryptoPP hashes**

                    4 - MD5

                    5 - SHA1

                    6 - Whirlpool

                    7 - CRC-32

                    8 - Adler-32

  --hashbits|-b  = use n bit hashes (for n's of 1 to WORDSIZE. default 26)

  --hashtype|-t  = select hash table data structure to use:

                    OpenHash  - open sub-word packing of hashbits

                    EcoHash   - chained sub-word packing of hashbits (default)

                    ArrayHash - malloc/realloc (fast but fragmenty)

                    MSetHash  - memory exorbanant, almost pointless.

  --probing      = 0 - linear, 1 - quadratic (default)

  --hitfilter|-h = minimum number of hits to be outputted as an anchor

                   (default 1)

  --histogram|-H = histogram computation level: (-H alone implies -H1)

                    0 - no histogram (default)

                    1 - basic bucketsize/bucketcount histogram data

    2 - bucket-based scores to anchors.detils

                    3 - perbucket count data

    4 - perbucket + perpattern count data

  --repeatmask|-r= skip repeat masked data (ie: lowercase atgc)

  --seedfilter|-f= skip seeds that occur more than N times

  --hashfilter|-m= like --seedfilter but works on hash keys instead of

                   seeds. May cause some collateral damage to otherwise

                   unique seeds, but it's faster. Also non-sequence-specific

                   so more like at best 1/N the tolerance of seedfilter.

  --rseed|-s     = random number seed for non-deterministic algorithms

                   (ie: the adative hash-finding). If you're doing any

                   performance comparisons, it's probably imperative that you

                   use the same seed for each run of the same settings.

                   Default is obtained from time() (ie: seconds since 1970).

  --skipfwd|-F   = Skip forward facing matches

  --skiprev|-R   = Skip reverse facing matches

  --skip1to1|-1  = Skip matches along the 1:1 line (good for comparing to self)

  --hashonly|-Q  = Hash Only. no anchors. just statistics.

  --hashskip|-S  = Hashes every n bases. (Default is 1. ie all)

                   Not supplying any argument increments the skip amount by 1.

  --hashCache|-c = Caches hash tables in the directory. (default: cache/)

  --join|-j      = Join anchors within n bases of eachother (default: 0)

                   Specifying a negative n implies -n*patternLength

  --bitscore|-B  = toggles compututation of a bitscore for all anchors

                   (default is on)

  --memory|-M    = set the target amount of total memory

                    (either in gb or as % total memory)

  --seedterms|-T = toggles retention of seed terms (defaults to off)

                    (these are necessary for computing TF-IDF scores)

  --sectime|-e   = always display times in seconds

  --repeatmap|-i = toggles keeping of a repeat map when --mergefilter

                   is used (defaults to yes).

  --mergefilter|-Y = filter out matches which would would cause more than N

                     many anchors to be generated from 1 seed (default -Y100).

                     Use -Y0 to disable.

  --scorefilter    = set a minimum ungapped score for seeds

  --tfidf|-k       = perform accurate tfidf scoring from within murasaki

                     (requires extra memory at anchor generation time)

  --reverseotf|-o  = generate reverse complement on the fly (defaults to on)

  --rifts|-/       = allow anchors to skip N sequences (default 0)

  --islands|-%     = same as --rifts=S-N (where S is number of seqs)

  --fuzzyextend|-z = enable (default) or disable fuzzy extension of hits

  --fuzzyextendlosslimit|-Z = set the cutoff at which to stop extending

                      fuzzy hits (ie. the BLAST X parameter).

  --gappedanchors  = use gapped (yes) or ungapped (no (default)) anchors.

  --scorebyminimumpair = do anchor scoring by minimum pair when appropriate

                     (default). Alternative is mean (somewhat illogical, but

                     theoretically faster).

  --binaryseq      = enable (default) or disable binary sequence read/write

 

Adaptive has function related:

  --hasherFairEntropy = use more balanced entropy estimation (default: yes)

  --hasherCorrelationAdjust = adjust entropy estimates for nearby sources

        assuming some correlation (default: yes)

  --hasherTargetGACycles = GA cycle cutoff

  --hasherEntropyAgro = how aggressive to be about pursuing maximum

        entropy hash functions (takes a real. default is 1).

...and of course

  --verbose|-v   = increases verbosity

  --version|-V   = prints version information and quits

  --help|-?      = prints this help message and quits

 

Platform information:

Wordsize: 64 bits

sizeof(word): 8 bytes

Total Memory: 22.97 GB

Available Memory: 23.93 GB (104.20%)

 

Murasaki version 1.68.6 (LARGESEQ, CRYPTOPP)

 

 

 

ラン

ゲノム・プラスミドのfastaかgenbakファイルを指定する。

 

3つのgbkファイルの比較。

./murasaki -p 19:26 -d output -n sample_name 7942.gbk GT-S.gbk 7002.gbk Leptolyngbya_dg5.gbk Nostoc_sp.PCC7120.gbk

#docker imageを使うなら
docker run --rm -itv $PWD:/data/ -w /data biocontainers/murasaki:v1.68.6-8b1-deb_cv1 murasaki A.fasta B.fasta C.fasta -p 19:26 -d output

出力ディレクトリoutputをGNVに読み込ませて描画すると、下のような図を出力できる。

f:id:kazumaxneo:20170623134158j:plain

4種の非常に近縁なシアノバクテリアを比較。右端にはANI値を載せた。

 

 

ランさせた時の様子はこのような感じになる。

動画では、genbankファイルを2つ選択してmurasakiで相同性を調べ、出力されたフォルダをGMVで表示している。

 

著者らによるアルゴリズムの説明とチュートリアル

Yasunori Osana @ University of the Ryukyus

こちらも貼っておきます。

Yasunori Osana @ University of the Ryukyus

 

 

 

引用

Murasaki: a fast, parallelizable algorithm to find anchors from multiple genomes

Popendorf K, Tsuyoshi H, Osana Y, Sakakibara Y

PLoS One. 2010 Sep 24;5(9):e12651

 

 

ゲノム比較 x 変異コール x ビューア を統合したGUI(CUI)ツール Mauve

 

mauveはよく似たゲノムのアライメントを行い、その結果を見やすいビューアで表示して比較できるソフトである。MacwindowsLinux版が用意されており、無償でダウンロードできる。

 

ダウンロードは公式サイトから行う。

the Darling lab | computational (meta)genomics

コマンドライン版もあり、そちらの方が細かくアライメントパラメータを選択できるが、ここでは使いやすいGUI版 (javaのプログラム) を説明する。

 

Align with progressive mauveを選ぶ。

f:id:kazumaxneo:20170623123720j:plain

 

FilesタブのAdd sequenceを選び比較する配列ファイルを選択。入力は、fasta、gbkに対応している。

f:id:kazumaxneo:20170627133520j:plain

   

Parametersタブでパラメータを設定。デフォルトで進める。

f:id:kazumaxneo:20170623124216j:plain

 

Scoresタブ。アライメントのスコアとペナルティ設定。デフォルトで進める。

f:id:kazumaxneo:20170623124304j:plain

最後にAlignを押しジョブを開始させる。

 

ラン中はジョブの進捗が別ウィンドウで表示される。

f:id:kazumaxneo:20170623124005j:plain

5Mbくらいのゲノム3つを比較すると、解析に10分くらいかかる。

 

3ゲノムを比較した動画が以下になる。一番上がcanuでE.coliのロングリードをアセンブルして作ったcontig.fa、2つ目がE.coliのgene bankファイル、3つ目がspadesでショートリードとロングリードをhybridして作ったcontig.faである。

動画のように、気になる領域に移動してアライメントを調べられる。genbankを入れておけば、orf情報を表示させることもできて便利である。動画の後半では,SNPs、indelを出力している。

 

ゲノムサイズが100Mbを超えると重く使いづらくなるので、ゲノムが小さい生物の比較に向いている印象を受けた。ただし、mauveの論文 (ref.1) ではmouseやヒトゲノムを比較しており、使えないことはない。

アライメントできるtoolは多いが、mauveのように各OSに対応しており、アライメントして結果をGUIで確認できるソフトはなかなかないと思う。よかったら試してみてください。

 

アセンブリした配列とリファレンス配列を比較する時は、contigのFASTAファイルをソートしてから比較した方が見やすくなります。ソートしてからprogressive mauveを走らせるには、以下の流れで実行してください。

 

 

引用

Mauve: Multiple Alignment of Conserved Genomic Sequence With Rearrangements

Aaron C.E. Darling,1,2,6 Bob Mau,2,3 Frederick R. Blattner,4,5 and Nicole T. Perna2,5

Genome Res. 2004 Jul; 14(7): 1394–1403. doi: 10.1101/gr.2289704

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC442156/

 

 

Oxford NanoporeリードのFAST5 => FASTA / FASTQ変換

2021 8/22 説明を修正

 

MNIONのシーケンスデータはFAST5というフォーマットで出力される。FAST5はHDF5(wiki)という時系列データ関係でよく使われる形式に乗っ取っている。塩基配列になっていないバイナリーなデータである。このデータから塩基配列に変換するツールがいくつか発表されている。導入して検証してみる。

  

Poretools

公式ページ https://poretools.readthedocs.io/en/latest/index.html

 

brewでインストールできる。

brew install poretools

 Bioconda環境なら"conda install poretools"で入る。

 

fastqの抽出。

poretools fastq <directory> > output.fq

--type bestの指定で2Dリードのみ出力する。他にtemplateかcomplementなどがある。デフォルトではallになっている。ディレクトリはnanopore/fast5/0/にあるなら、~0/まで記載しないとエラーになる。

 

fastaの抽出。

poretools fasta <directory> > output.fq

 

シーケンスのイールド。電圧の変更でポアが急激に減っていっていないかチェックするのに使える。

poretools yield_plot --plot-type reads <directory> --saveas output.pdf

 --saveas: PDFかpngで保存。オプション無しだと画面出力される。

f:id:kazumaxneo:20170621220955j:plain

 

リード情報。

poretools stats <directory>

例えば以下のような画面が出力される。

total reads 4000

total base pairs 9561055

mean 2390.26

median 1182

min 11

max 148048

N25 10196

N50 4566

N75 1919

 

リード長のヒストグラム

poretools hist <directory> --saveas output.pdf

--num-bins、 --max-length、 --min-lengthで間隔や下限、上限を指定することもできる。

f:id:kazumaxneo:20170621222115j:plain

塩基組成。

poretools nucdist <directory>

例えば以下のような画面が出力される。

A 2242285 9561055 0.234522759256

C 2455787 9561055 0.256853140161

T 2375819 9561055 0.248489209611

G 2487164 9561055 0.260134890972

 

クオリティスコア分布。

poretools qualdist <directory>

以下のような画面が出力される。

" 1 19774 9561055 0.0020681818063

# 2 47996 9561055 0.00501994811242

$ 3 118476 9561055 0.0123915195551

.

.

.

H 39 1074 9561055 0.000112330699907

 

クオリティスコア分布のbox plot。

poretools qualpos <directory> --saveas output.pdf

f:id:kazumaxneo:20170621222741j:plain

 

フローセルの各ポアの総合パフォーマンス。

poretools occupancy <directory>

f:id:kazumaxneo:20170621223139j:plain

他にも、塩基変換前の生データやsquiggle plot を出力する機能などを持っている。公式サイトで確認してください。

 

Nanopolish

こちらでインストールは紹介している。


 

PoreSeq

Githubページ

GitHub - tszalay/poreseq: Error correction and variant calling algorithm for nanopore sequencing

githubからダウンロードして、PoreSeqのルートでsudo pip install -e .でインストールできる。

mac OSXで動作はするが、内部でLASTとsamtoolsを動かしているらしく、samtoolsのsort付近でエラーになりファイルが消える。現在のところ未修正。

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Oxford Nanoporeリードのアセンブリ パフォーマンス比較

 

S. cerevisiaeとE.coli K12をilluminaとnanoporeで読んだ論文のデータ が公開されている。

http://schatzlab.cshl.edu/data/nanocorr/

 

このデータを使い、2017年6月現在のナノポアリードのパフォーマンスを調べてみる。

上記URLからMinionとilluminaのシーケンスデータ (fastq) をダウンロードした。 

 

 

E.coliのリード数などは以下の通り。

・Oxford nanopore Minion: 59008 reads、total_length 240Mb、average 4068 bp、max 43000bp、min 120bp

・illumina Miseq (301x2): 22Mx2 reads、total_length 13.4Gb、max 301 bp

ナノポアのリードは平均4000bpを超えているが、リードのトータルサイズが2.4億bpで、カバレッジは48xにしかならない。一方でMiseqのトータルサイズは301bp計算で67億bpを超えている!。カバレッジを計算すると2680xで、かなりディープに読まれたデータだ。gzから解凍したファイルサイズはそれぞれ15GB近くある。

ロングリードのサイズ分布をヒストグラム解析して調べてみた。

f:id:kazumaxneo:20170622150029j:plain

横軸がリードの長さ(bp)で、縦軸がその長さのリードが出てくる頻度である。(右端は単位を5000bpに変えていることに注意してください)。一番多いのは1000bp程度のリードであるが、綺麗な二項分布にはなっていない。6000-7000付近にも弱いピークがあるように見える。

やはりとてつもなく長い。1リードで5000bp以上の長さがあれば、rrnaのクラスターもいきなりアセンブルできる可能性が高い。

 

S. cerevisiaeのリード数は以下の通り。

・Oxford nanopore Minion: 88218 reads、total_length 526Mb、average 5969 bp、max 72879bp、min 103bp

・illumina Miseq (301x2): 25Mx2 reads、total_length 15Gb、max 301 bp

S. cerevisiaeのゲノムサイズを大雑把に12.5Mbとした時、ナノポアのカバレッジは44x、Miseqのカバレッジは1300xとなる。

ナノポアのlong readのサイズ分布をヒストグラムで調べる。

f:id:kazumaxneo:20170622150005j:plain

横軸がリードの長さ(bp)で、縦軸がその長さのリードが出てくる頻度である。(同じく右端は単位を5000bpに変えていることに注意してください)。一番多いのは1000bp程度のリードであるが、7000付近にも強いピークがある。

推測だが、7000bp付近のリードが本来のナノポアMinionのパフォーマンスで、左端のピークはDNA抽出時に物理的に切れてしまったDNAではないだろうか?(海外研究者のナノポアのワークショップの投稿で、λ HInd III cutをナノポアで読むと、等カバレッジにならずサイズが長くなるほどカバレッジが下がという話がある)。

 

 

アセンブル

以下はE.coliの例であるが、 S. cerevisiaeもゲノムサイズ以外は同じパラメータで実行した。 

1、canuを使いONTリードだけでアセンブル。

canu -p test -d long_read genomeSize=4.8m -nanopore-raw oxt.fq

E.coli=> 3 contig

S. cerevisiae=> 27 contig

contigの長さをグラフにすると以下のようになる。

f:id:kazumaxneo:20170623112456j:plain

グラフの左からcontigの長さ順に示す。縦軸がcontigの長さ(bp)である。

 

S.cerevisiaeは17クロモソーム持つので、クロモソームのサイズとcontigの長さを比較した。下の表は左側がS. cerevisiaeのゲノムをサイズ順に並べたもので、右側がナノポアリードをアセンブルして作ったcontigをサイズ順に並べたものである。

chr Size (bp)   contig Size (bp)
IV 1530000   tig00000001  1546985
VII 1090000   tig00000003  1091026
XV 1090000   tig00000007  1086779
XII 1080000   tig00000006  939389
XVI 948066   tig00000013  805656
XIII 924431   tig00000022  737318
II 813184   tig00013668  678294
XIV 784333   tig00000024  644519
X 745751   tig00000053  573971
XI 666816   tig00013665  519780
V 576874   tig00000181  487116
VIII 562643   tig00000055  462121
IX 439888   tig00000058  429563
III 316620   tig00000072  271056
VI 270161   tig00000185  259055
I 230218   tig00000066  239610
MT 85779   tig00000183  237080
      tig00000050  190760
      tig00000075  188868
      tig00000035  187180
      tig00000184  162871
      tig00000086  155346
      tig00000042  140091
      tig00000036  110210
      tig00013670  86062
      tig00000009  52963
      tig00000046  25807

おそらく多くのクロモソームは1つにアセンブルされている。quastでエラーチェックも行なったが、リファレンスゲノムと比べて明確な構造変化は起こっていない。つまり、あとはPCRを総当たりすれば、ゲノムの構造は決まってしまう。かなり衝撃的な結果である。(ただしエラーチェックは入念に行う必要がある!)

E.coli K12の方は、3つのcontigにアセンブルされていた(E.coli K12 MG1655のゲノムは4.64Mbのクロモソームのみを持つ。プラスミドはない)。

 

  

2、spadesを使いilluminaのペアリードをアセンブル。

spades

E.coli=> 763 contig

S. cerevisiae=> 1389 contig

f:id:kazumaxneo:20170623111700j:plain

 

lengthのtop100に絞って表示する。縦軸がcontigの長さ(bp)である。

 

 

3、spadesを使いONTリードとilluminaリード両方でhybirdアセンブル。

spades

E.coli=>  602 contig

S. cerevisiae=>  1145 contig

E.coliのcontigは一番長いもので4.64Mbで、残りは2kb以下のcontigだった。グラフにすると下のようになる。

f:id:kazumaxneo:20170623112142j:plain

lengthのtop30に絞って表示する。縦軸がcontigの長さ(bp)である。

 

 

 

 

 

 

  

ショートリードのアダプタートリミングツール Trim Galore

2019 5/8 インストールおよびヘルプ追記

2020 12/9 help更新

 

これまで様々なアダプタートリミングツールが報告されてきている。OMIC toolsで検索すると、2017年6月で35件ヒットする(OMIC toolリンク)。その中でもFastQC、cutadapt、Fastx-toolkitなどはよく耳にする。Trim Galore!はFastQCとcutadaptを内部で動かし、fastqから自動でアダプター配列を認識してトリムするツール。Biostarのhit数を見る限り使っている人も多そうである。インストールして動作を見ていく。

 

 

インストール

公式ページリンク Babraham Bioinformatics - Trim Galore!

Trim Galore!本体はperlスクリプトで、内部でFastQCとcutadaptを動かしてトリムを行う。そのため動作にはFastQCとcutadaptが必要である。この2つを先にインストールしておく。

 

Biocondaに登録されているので、依存も含めcondaで導入するのが一番楽。

mamba install -c bioconda -y trim-galore

 

 

またはbrewで依存導入。依存のFastQCとcutadaptをインストールする。

brew install cutadapt
brew install FastQC

 

FastQCとcutadaptのインストールが終わったら、公式サイトからTrim Galore!をダウンロードして解凍する。解凍したディレクトリのトップ階層にあるtrim_galoreを/usr/local/bin/に移動する。

ln -s TrimGalore-0.4.3/trim_galore /usr/local/bin

リンクかコピーなどでパスを通す。

> trim_galore -h

$ trim_galore -h

Option h is ambiguous (hardtrim3, hardtrim5, help)

Please respecify command line options

(base) kamisakakazumanoMac-mini:cyanobacteria_RNAseq kazu$ trim_galore -help

 

 USAGE:

 

trim_galore [options] <filename(s)>

 

 

-h/--help               Print this help message and exits.

 

-v/--version            Print the version information and exits.

 

-q/--quality <INT>      Trim low-quality ends from reads in addition to adapter removal. For

                        RRBS samples, quality trimming will be performed first, and adapter

                        trimming is carried in a second round. Other files are quality and adapter

                        trimmed in a single pass. The algorithm is the same as the one used by BWA

                        (Subtract INT from all qualities; compute partial sums from all indices

                        to the end of the sequence; cut sequence at the index at which the sum is

                        minimal). Default Phred score: 20.

 

--phred33               Instructs Cutadapt to use ASCII+33 quality scores as Phred scores

                        (Sanger/Illumina 1.9+ encoding) for quality trimming. Default: ON.

 

--phred64               Instructs Cutadapt to use ASCII+64 quality scores as Phred scores

                        (Illumina 1.5 encoding) for quality trimming.

 

--fastqc                Run FastQC in the default mode on the FastQ file once trimming is complete.

 

--fastqc_args "<ARGS>"  Passes extra arguments to FastQC. If more than one argument is to be passed

                        to FastQC they must be in the form "arg1 arg2 etc.". An example would be:

                        --fastqc_args "--nogroup --outdir /home/". Passing extra arguments will

                        automatically invoke FastQC, so --fastqc does not have to be specified

                        separately.

 

-a/--adapter <STRING>   Adapter sequence to be trimmed. If not specified explicitly, Trim Galore will

                        try to auto-detect whether the Illumina universal, Nextera transposase or Illumina

                        small RNA adapter sequence was used. Also see '--illumina', '--nextera' and

                        '--small_rna'. If no adapter can be detected within the first 1 million sequences

                        of the first file specified or if there is a tie between several adapter sequences,

                        Trim Galore defaults to '--illumina' (as long as the Illumina adapter was one of the

                        options, else '--nextera' is the default). A single base

                        may also be given as e.g. -a A{10}, to be expanded to -a AAAAAAAAAA.

 

-a2/--adapter2 <STRING> Optional adapter sequence to be trimmed off read 2 of paired-end files. This

                        option requires '--paired' to be specified as well. If the libraries to be trimmed

                        are smallRNA then a2 will be set to the Illumina small RNA 5' adapter automatically

                        (GATCGTCGGACT). A single base may also be given as e.g. -a2 A{10}, to be expanded

                        to -a2 AAAAAAAAAA.

 

--illumina              Adapter sequence to be trimmed is the first 13bp of the Illumina universal adapter

                        'AGATCGGAAGAGC' instead of the default auto-detection of adapter sequence.

 

--nextera               Adapter sequence to be trimmed is the first 12bp of the Nextera adapter

                        'CTGTCTCTTATA' instead of the default auto-detection of adapter sequence.

 

--small_rna             Adapter sequence to be trimmed is the first 12bp of the Illumina Small RNA 3' Adapter

                        'TGGAATTCTCGG' instead of the default auto-detection of adapter sequence. Selecting

                        to trim smallRNA adapters will also lower the --length value to 18bp. If the smallRNA

                        libraries are paired-end then a2 will be set to the Illumina small RNA 5' adapter

                        automatically (GATCGTCGGACT) unless -a 2 had been defined explicitly.

 

--consider_already_trimmed <INT>     During adapter auto-detection, the limit set by <INT> allows the user to 

                        set a threshold up to which the file is considered already adapter-trimmed. If no adapter

                        sequence exceeds this threshold, no additional adapter trimming will be performed (technically,

                        the adapter is set to '-a X'). Quality trimming is still performed as usual.

                        Default: NOT SELECTED (i.e. normal auto-detection precedence rules apply).                     

 

--max_length <INT>      Discard reads that are longer than <INT> bp after trimming. This is only advised for

                        smallRNA sequencing to remove non-small RNA sequences.

 

 

--stringency <INT>      Overlap with adapter sequence required to trim a sequence. Defaults to a

                        very stringent setting of 1, i.e. even a single bp of overlapping sequence

                        will be trimmed off from the 3' end of any read.

 

-e <ERROR RATE>         Maximum allowed error rate (no. of errors divided by the length of the matching

                        region) (default: 0.1)

 

--gzip                  Compress the output file with GZIP. If the input files are GZIP-compressed

                        the output files will automatically be GZIP compressed as well. As of v0.2.8 the

                        compression will take place on the fly.

 

--dont_gzip             Output files won't be compressed with GZIP. This option overrides --gzip.

 

--length <INT>          Discard reads that became shorter than length INT because of either

                        quality or adapter trimming. A value of '0' effectively disables

                        this behaviour. Default: 20 bp.

 

                        For paired-end files, both reads of a read-pair need to be longer than

                        <INT> bp to be printed out to validated paired-end files (see option --paired).

                        If only one read became too short there is the possibility of keeping such

                        unpaired single-end reads (see --retain_unpaired). Default pair-cutoff: 20 bp.

 

--max_n COUNT           The total number of Ns (as integer) a read may contain before it will be removed altogether.

                        In a paired-end setting, either read exceeding this limit will result in the entire

                        pair being removed from the trimmed output files.

 

--trim-n                Removes Ns from either side of the read. This option does currently not work in RRBS mode.

 

-o/--output_dir <DIR>   If specified all output will be written to this directory instead of the current

                        directory. If the directory doesn't exist it will be created for you.

 

--no_report_file        If specified no report file will be generated.

 

--suppress_warn         If specified any output to STDOUT or STDERR will be suppressed.

 

--clip_R1 <int>         Instructs Trim Galore to remove <int> bp from the 5' end of read 1 (or single-end

                        reads). This may be useful if the qualities were very poor, or if there is some

                        sort of unwanted bias at the 5' end. Default: OFF.

 

--clip_R2 <int>         Instructs Trim Galore to remove <int> bp from the 5' end of read 2 (paired-end reads

                        only). This may be useful if the qualities were very poor, or if there is some sort

                        of unwanted bias at the 5' end. For paired-end BS-Seq, it is recommended to remove

                        the first few bp because the end-repair reaction may introduce a bias towards low

                        methylation. Please refer to the M-bias plot section in the Bismark User Guide for

                        some examples. Default: OFF.

 

--three_prime_clip_R1 <int>     Instructs Trim Galore to remove <int> bp from the 3' end of read 1 (or single-end

                        reads) AFTER adapter/quality trimming has been performed. This may remove some unwanted

                        bias from the 3' end that is not directly related to adapter sequence or basecall quality.

                        Default: OFF.

 

--three_prime_clip_R2 <int>     Instructs Trim Galore to remove <int> bp from the 3' end of read 2 AFTER

                        adapter/quality trimming has been performed. This may remove some unwanted bias from

                        the 3' end that is not directly related to adapter sequence or basecall quality.

                        Default: OFF.

 

--2colour/--nextseq INT This enables the option '--nextseq-trim=3'CUTOFF' within Cutadapt, which will set a quality

                        cutoff (that is normally given with -q instead), but qualities of G bases are ignored.

                        This trimming is in common for the NextSeq- and NovaSeq-platforms, where basecalls without

                        any signal are called as high-quality G bases. This is mutually exlusive with '-q INT'.

 

 

--path_to_cutadapt </path/to/cutadapt>     You may use this option to specify a path to the Cutadapt executable,

                        e.g. /my/home/cutadapt-1.7.1/bin/cutadapt. Else it is assumed that Cutadapt is in

                        the PATH.

 

--basename <PREFERRED_NAME> Use PREFERRED_NAME as the basename for output files, instead of deriving the filenames from

                        the input files. Single-end data would be called PREFERRED_NAME_trimmed.fq(.gz), or

                        PREFERRED_NAME_val_1.fq(.gz) and PREFERRED_NAME_val_2.fq(.gz) for paired-end data. --basename

                        only works when 1 file (single-end) or 2 files (paired-end) are specified, but not for longer lists.

 

-j/--cores INT          Number of cores to be used for trimming [default: 1]. For Cutadapt to work with multiple cores, it

                        requires Python 3 as well as parallel gzip (pigz) installed on the system. The version of Python used 

                        is detected from the shebang line of the Cutadapt executable (either 'cutadapt', or a specified path).

                        If Python 2 is detected, --cores is set to 1.

                        If pigz cannot be detected on your system, Trim Galore reverts to using gzip compression. Please note

                        that gzip compression will slow down multi-core processes so much that it is hardly worthwhile, please 

                        see: https://github.com/FelixKrueger/TrimGalore/issues/16#issuecomment-458557103 for more info).

 

                        Actual core usage: It should be mentioned that the actual number of cores used is a little convoluted.

                        Assuming that Python 3 is used and pigz is installed, --cores 2 would use 2 cores to read the input

                        (probably not at a high usage though), 2 cores to write to the output (at moderately high usage), and 

                        2 cores for Cutadapt itself + 2 additional cores for Cutadapt (not sure what they are used for) + 1 core

                        for Trim Galore itself. So this can be up to 9 cores, even though most of them won't be used at 100% for

                        most of the time. Paired-end processing uses twice as many cores for the validation (= writing out) step.

                        --cores 4 would then be: 4 (read) + 4 (write) + 4 (Cutadapt) + 2 (extra Cutadapt) + 1 (Trim Galore) = 15.

 

                        It seems that --cores 4 could be a sweet spot, anything above has diminishing returns.

 

 

 

SPECIFIC TRIMMING - without adapter/quality trimming

 

--hardtrim5 <int>       Instead of performing adapter-/quality trimming, this option will simply hard-trim sequences

                        to <int> bp at the 5'-end. Once hard-trimming of files is complete, Trim Galore will exit.

                        Hard-trimmed output files will end in .<int>_5prime.fq(.gz). Here is an example:

 

                        before:         CCTAAGGAAACAAGTACACTCCACACATGCATAAAGGAAATCAAATGTTATTTTTAAGAAAATGGAAAAT

                        --hardtrim5 20: CCTAAGGAAACAAGTACACT

 

--hardtrim3 <int>       Instead of performing adapter-/quality trimming, this option will simply hard-trim sequences

                        to <int> bp at the 3'-end. Once hard-trimming of files is complete, Trim Galore will exit.

                        Hard-trimmed output files will end in .<int>_3prime.fq(.gz). Here is an example:

 

                        before:         CCTAAGGAAACAAGTACACTCCACACATGCATAAAGGAAATCAAATGTTATTTTTAAGAAAATGGAAAAT

                        --hardtrim3 20:                                                   TTTTTAAGAAAATGGAAAAT

 

--clock                 In this mode, reads are trimmed in a specific way that is currently used for the Mouse

                        Epigenetic Clock (see here: Multi-tissue DNA methylation age predictor in mouse, Stubbs et al.,

                        Genome Biology, 2017 18:68 https://doi.org/10.1186/s13059-017-1203-5). Following this, Trim Galore

                        will exit.

 

                        In it's current implementation, the dual-UMI RRBS reads come in the following format:

 

                        Read 1  5' UUUUUUUU CAGTA FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF TACTG UUUUUUUU 3'

                        Read 2  3' UUUUUUUU GTCAT FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF ATGAC UUUUUUUU 5'

 

                        Where UUUUUUUU is a random 8-mer unique molecular identifier (UMI), CAGTA is a constant region,

                        and FFFFFFF... is the actual RRBS-Fragment to be sequenced. The UMIs for Read 1 (R1) and

                        Read 2 (R2), as well as the fixed sequences (F1 or F2), are written into the read ID and

                        removed from the actual sequence. Here is an example:

 

                        R1: @HWI-D00436:407:CCAETANXX:1:1101:4105:1905 1:N:0: CGATGTTT

                            ATCTAGTTCAGTACGGTGTTTTCGAATTAGAAAAATATGTATAGAGGAAATAGATATAAAGGCGTATTCGTTATTG

                        R2: @HWI-D00436:407:CCAETANXX:1:1101:4105:1905 3:N:0: CGATGTTT

                            CAATTTTGCAGTACAAAAATAATACCTCCTCTATTTATCCAAAATCACAAAAAACCACCCACTTAACTTTCCCTAA

 

                        R1: @HWI-D00436:407:CCAETANXX:1:1101:4105:1905 1:N:0: CGATGTTT:R1:ATCTAGTT:R2:CAATTTTG:F1:CAGT:F2:CAGT

                                         CGGTGTTTTCGAATTAGAAAAATATGTATAGAGGAAATAGATATAAAGGCGTATTCGTTATTG

                        R2: @HWI-D00436:407:CCAETANXX:1:1101:4105:1905 3:N:0: CGATGTTT:R1:ATCTAGTT:R2:CAATTTTG:F1:CAGT:F2:CAGT

                                         CAAAAATAATACCTCCTCTATTTATCCAAAATCACAAAAAACCACCCACTTAACTTTCCCTAA

 

                        Following clock trimming, the resulting files (.clock_UMI.R1.fq(.gz) and .clock_UMI.R2.fq(.gz))

                        should be adapter- and quality trimmed with Trim Galore as usual. In addition, reads need to be trimmed

                        by 15bp from their 3' end to get rid of potential UMI and fixed sequences. The command is:

 

                        trim_galore --paired --three_prime_clip_R1 15 --three_prime_clip_R2 15 *.clock_UMI.R1.fq.gz *.clock_UMI.R2.fq.gz

 

                        Following this, reads should be aligned with Bismark and deduplicated with UmiBam

                        in '--dual_index' mode (see here: https://github.com/FelixKrueger/Umi-Grinder). UmiBam recognises

                        the UMIs within this pattern: R1:(ATCTAGTT):R2:(CAATTTTG): as (UMI R1) and (UMI R2).

 

--polyA                 This is a new, still experimental, trimming mode to identify and remove poly-A tails from sequences.

                        When --polyA is selected, Trim Galore attempts to identify from the first supplied sample whether

                        sequences contain more often a stretch of either 'AAAAAAAAAA' or 'TTTTTTTTTT'. This determines

                        if Read 1 of a paired-end end file, or single-end files, are trimmed for PolyA or PolyT. In case of

                        paired-end sequencing, Read2 is trimmed for the complementary base from the start of the reads. The

                        auto-detection uses a default of A{20} for Read1 (3'-end trimming) and T{150} for Read2 (5'-end trimming).

                        These values may be changed manually using the options -a and -a2.

 

                        In addition to trimming the sequences, white spaces are replaced with _ and it records in the read ID

                        how many bases were trimmed so it can later be used to identify PolyA trimmed sequences. This is currently done

                        by writing tags to both the start ("32:A:") and end ("_PolyA:32") of the reads in the following example:

 

                        @READ-ID:1:1102:22039:36996 1:N:0:CCTAATCC

                        GCCTAAGGAAACAAGTACACTCCACACATGCATAAAGGAAATCAAATGTTATTTTTAAGAAAATGGAAAATAAAAACTTTATAAACACCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA

 

                        @32:A:READ-ID:1:1102:22039:36996_1:N:0:CCTAATCC_PolyA:32

                        GCCTAAGGAAACAAGTACACTCCACACATGCATAAAGGAAATCAAATGTTATTTTTAAGAAAATGGAAAATAAAAACTTTATAAACACC

 

                        PLEASE NOTE: The poly-A trimming mode expects that sequences were both adapter and quality trimmed

                        before looking for Poly-A tails, and it is the user's responsibility to carry out an initial round of

                        trimming. The following sequence:

 

                        1) trim_galore file.fastq.gz

                        2) trim_galore --polyA file_trimmed.fq.gz

                        3) zcat file_trimmed_trimmed.fq.gz | grep -A 3 PolyA | grep -v ^-- > PolyA_trimmed.fastq

 

                        Will 1) trim qualities and Illumina adapter contamination, 2) find and remove PolyA contamination.

                        Finally, if desired, 3) will specifically find PolyA trimmed sequences to a specific FastQ file of your choice.

 

--implicon              This is a special mode of operation for paired-end data, such as required for the IMPLICON method, where a UMI sequence

                        is getting transferred from the start of Read 2 to the readID of both reads. Following this, Trim Galore will exit.

 

                        In it's current implementation, the UMI carrying reads come in the following format:

 

                        Read 1  5' FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF 3'

                        Read 2  3' UUUUUUUUFFFFFFFFFFFFFFFFFFFFFFFFFFFF 5'

 

                        Where UUUUUUUU is a random 8-mer unique molecular identifier (UMI) and FFFFFFF... is the actual fragment to be

                        sequenced. The UMI of Read 2 (R2) is written into the read ID of both reads and removed from the actual sequence.

                        Here is an example:

 

                        R1: @HWI-D00436:407:CCAETANXX:1:1101:4105:1905 1:N:0: CGATGTTT

                            ATCTAGTTCAGTACGGTGTTTTCGAATTAGAAAAATATGTATAGAGGAAATAGATATAAAGGCGTATTCGTTATTG

                        R2: @HWI-D00436:407:CCAETANXX:1:1101:4105:1905 3:N:0: CGATGTTT

                            CAATTTTGCAGTACAAAAATAATACCTCCTCTATTTATCCAAAATCACAAAAAACCACCCACTTAACTTTCCCTAA

                        

                        After --implicon trimming:

                        R1: @HWI-D00436:407:CCAETANXX:1:1101:4105:1905 1:N:0: CGATGTTT:CAATTTTG

                            ATCTAGTTCAGTACGGTGTTTTCGAATTAGAAAAATATGTATAGAGGAAATAGATATAAAGGCGTATTCGTTATTG

                        R2: @HWI-D00436:407:CCAETANXX:1:1101:4105:1905 3:N:0: CGATGTTT:CAATTTTG

                                    CAGTACAAAAATAATACCTCCTCTATTTATCCAAAATCACAAAAAACCACCCACTTAACTTTCCCTAA

 

RRBS-specific options (MspI digested material):

 

--rrbs                  Specifies that the input file was an MspI digested RRBS sample (recognition

                        site: CCGG). Single-end or Read 1 sequences (paired-end) which were adapter-trimmed

                        will have a further 2 bp removed from their 3' end. Sequences which were merely

                        trimmed because of poor quality will not be shortened further. Read 2 of paired-end

                        libraries will in addition have the first 2 bp removed from the 5' end (by setting

                        '--clip_r2 2'). This is to avoid using artificial methylation calls from the filled-in

                        cytosine positions close to the 3' MspI site in sequenced fragments.

                        This option is not recommended for users of the NuGEN ovation RRBS System 1-16

                        kit (see below).

 

--non_directional       Selecting this option for non-directional RRBS libraries will screen

                        quality-trimmed sequences for 'CAA' or 'CGA' at the start of the read

                        and, if found, removes the first two basepairs. Like with the option

                        '--rrbs' this avoids using cytosine positions that were filled-in

                        during the end-repair step. '--non_directional' requires '--rrbs' to

                        be specified as well. Note that this option does not set '--clip_r2 2' in

                        paired-end mode.

 

--keep                  Keep the quality trimmed intermediate file. Default: off, which means

                        the temporary file is being deleted after adapter trimming. Only has

                        an effect for RRBS samples since other FastQ files are not trimmed

                        for poor qualities separately.

 

 

Note for RRBS using the NuGEN Ovation RRBS System 1-16 kit:

 

Owing to the fact that the NuGEN Ovation kit attaches a varying number of nucleotides (0-3) after each MspI

site Trim Galore should be run WITHOUT the option --rrbs. This trimming is accomplished in a subsequent

diversity trimming step afterwards (see their manual).

 

 

 

Note for RRBS using MseI:

 

If your DNA material was digested with MseI (recognition motif: TTAA) instead of MspI it is NOT necessary

to specify --rrbs or --non_directional since virtually all reads should start with the sequence

'TAA', and this holds true for both directional and non-directional libraries. As the end-repair of 'TAA'

restricted sites does not involve any cytosines it does not need to be treated especially. Instead, simply

run Trim Galore! in the standard (i.e. non-RRBS) mode.

 

 

 

 

Paired-end specific options:

 

--paired                This option performs length trimming of quality/adapter/RRBS trimmed reads for

                        paired-end files. To pass the validation test, both sequences of a sequence pair

                        are required to have a certain minimum length which is governed by the option

                        --length (see above). If only one read passes this length threshold the

                        other read can be rescued (see option --retain_unpaired). Using this option lets

                        you discard too short read pairs without disturbing the sequence-by-sequence order

                        of FastQ files which is required by many aligners.

 

                        Trim Galore! expects paired-end files to be supplied in a pairwise fashion, e.g.

                        file1_1.fq file1_2.fq SRR2_1.fq.gz SRR2_2.fq.gz ... .

 

-t/--trim1              Trims 1 bp off every read from its 3' end. This may be needed for FastQ files that

                        are to be aligned as paired-end data with Bowtie. This is because Bowtie (1) regards

                        alignments like this:

 

                          R1 --------------------------->     or this:    ----------------------->  R1

                          R2 <---------------------------                       <-----------------  R2

 

                        as invalid (whenever a start/end coordinate is contained within the other read).

                        NOTE: If you are planning to use Bowtie2, BWA etc. you don't need to specify this option.

 

--retain_unpaired       If only one of the two paired-end reads became too short, the longer

                        read will be written to either '.unpaired_1.fq' or '.unpaired_2.fq'

                        output files. The length cutoff for unpaired single-end reads is

                        governed by the parameters -r1/--length_1 and -r2/--length_2. Default: OFF.

 

-r1/--length_1 <INT>    Unpaired single-end read length cutoff needed for read 1 to be written to

                        '.unpaired_1.fq' output file. These reads may be mapped in single-end mode.

                        Default: 35 bp.

 

-r2/--length_2 <INT>    Unpaired single-end read length cutoff needed for read 2 to be written to

                        '.unpaired_2.fq' output file. These reads may be mapped in single-end mode.

                        Default: 35 bp.

 

Last modified on 11 May 2020.

 

 

実行方法

trim_galore <fastq>

でランできる。

ペアリードなら両方のデータを指定し、--pairedのオプションをつけて実行する。

trim_galore --paired R1.fq R2.fq

出力は入力ファイル名_val_1.fq、入力ファイル名_val_2.fqとなる。

 

自動で認識させないでアダプター配列を明示するには-a <seq>をつける。配列を明示させるオプションは以下のようなものがある。

  • --illumina Adapter sequence to be trimmed is the first 13bp of the Illumina universal adapter 'AGATCGGAAGAGC' instead of the default auto-detection of adapter sequence.
  • --nextera Adapter sequence to be trimmed is the first 12bp of the Nextera adapter 'CTGTCTCTTATA' instead of the default auto-detection of adapter sequence.
  • --small_rna Adapter sequence to be trimmed is the first 12bp of the Illumina Small RNA 3' Adapter 'TGGAATTCTCGG' instead of the default auto-detection of adapter sequence. 

他にも様々なオプションがある。オプションを見るには

trim_galore --help

 

Trim Galore!のデフォルト動作はphred scoreをphred33と認識して行うが、古いIllumina 1.5 encodingのシーケンスデータを使う場合、--phred64をつけて明示する必要がある。

 

 

Nextraで調整されMiseqでランされたペアードエンドデータでテストすると、合計323万bp(データの4%)がトリムされた。

  >>> Now performing quality (cutoff 20) and adapter trimming in a single pass for the adapter sequence: 'AGATCGGAAGAGC' from file T1second_R2.fastq <<< 

This is cutadapt 1.8.3 with Python 2.7.10

Command line parameters: -f fastq -e 0.1 -q 20 -O 1 -a AGATCGGAAGAGC T1second_R2.fastq

Trimming 1 adapter with at most 10.0% errors in single-end mode ...

Finished in 5.59 s (21 us/read; 2.81 M reads/minute).

 

=== Summary ===

 

Total reads processed:                 261,774

Reads with adapters:                    69,845 (26.7%)

Reads written (passing filters):       261,774 (100.0%)

 

Total basepairs processed:    77,820,790 bp

Quality-trimmed:               3,237,687 bp (4.2%)

Total written (filtered):     74,481,043 bp (95.7%)

 

=== Adapter 1 ===

 

Sequence: AGATCGGAAGAGC; Type: regular 3'; Length: 13; Trimmed: 69845 times.

 

No. of allowed errors:

0-9 bp: 0; 10-13 bp: 1

 

Bases preceding removed adapters:

  A: 38.0%

  C: 25.5%

  G: 15.9%

  T: 20.6%

  none/other: 0.0%

 

Overview of removed sequences

length count expect max.err error counts

1 57208 65443.5 0 57208

2 9072 16360.9 0 9072

3 2290 4090.2 0 2290

4 536 1022.6 0 536

5 110 255.6 0 110

6 40 63.9 0 40

7 117 16.0 0 117

8 21 4.0 0 21

9 39 1.0 0 31 8

10 47 0.2 1 29 18

 

引用

http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/

 

関連


2022/03/13

Atriaの論文で性能が劣ることが示されています。FastpやAtriaなどへ切り替えも検討して下さい。