Help us understand the problem. What is going on with this article?

Qiime2で自分のサンプルを解析する

Qiime2で自分のサンプルを解析していく

前回書いたQiime2のMoving Pictures Tutorialの反響が予想以上に大きかったので、自分への備忘録も兼ねてQiime2で自分のサンプルを解析していこうと思う。これは卒論用のデータにする予定。ITSもやらないといけないので、そのうちITSのやり方も書くつもりである。

自分のサンプルを解析する場合とMoving Pictures Tutorialとの違いで気をつけないといけない点は何かといえば、シーケンス方法の違い、データベースの違いになると思う。例えば、Pictures Tutorialはシングルエンドでシーケンスしているが、最近のMiSeqのスタンダードはどちらかとえいばペアエンドで読むパターンだと思う。さらに、クオリティーコントロールの際に行うトリミングでバリアントとなる部分を切り落としたら全く意味ないので、自分がどのような領域をシーケンスしたのかを生物学的に考えて切り落とす必要がある。Illumina社のPDFに大まかにその話は書いてあった。また、データベースの違いについては、間も無くSilva138のQiime版が公開されると思うが、Tutorialで利用されるデータベースは非常に小さいものだと思う。自分できちんと情報が追加されたデータベースに変更し利用するという作業が必要になる。

また、今回利用するサンプルはNGSサービスに提出したものだが、中にさらにバーコードが仕込まれており、例えば1ランを3サンプルでやっているような感じになっている。その辺も自分でDemultiplexしていく。なお、Forwardのみにタグ情報を仕込ませている。前回同様、いつもの棒グラフを書いてみる。

Demultiplexのやり方 -Sabre

当方はSabreを利用している。所属研究室にオリジナルのスクリプトもあるようだが、個人的にはSabreの方が使いやすい感があるので、そうしている。なお、Claidentでも出来る。所属研究室で16S解析をする人は5人くらいいるが、うち3人はClaidentでタグ分け作業を行い、1人はオリジナルスクリプト、当方のみがSabreを利用する。SabreについてはいつもおなじみのMacでインフォマティクスにやり方が載っていた。Githubにも説明はある。

git clone https://github.com/najoshi/sabre.git
cd sabre/
make
./sabre --help

基本的にパスは通さず使っている。なので、基本的にアウトプットされるディレクトリはSabreのファイルに中になる。こちらの方が個人的に分かりやすかったので、そうしている。Sabreで必要なのはMultiplexされたFastqファイルとBarcode情報が載ったtxtファイルである。txtファイルは以下のように記述する。アウトプット名はここに書いている名前になる。

barcode1 barcode1_output_file1.fastq barcode1_output_file2.fastq
barcode2 barcode2_output_file1.fastq barcode2_output_file2.fastq

当方はこれが最初理解できなかった。正確には理解するのに3日を要した。したがって、実際のtxtの一部を載せる。

AGGTCT  611_S56_L001_R1_001.fastq 611_S56_L001_R2_001.fastq
GCGTGA  612_S56_L001_R1_001.fastq 612_S56_L001_R2_001.fastq
CGTATC  613_S56_L001_R1_001.fastq 613_S56_L001_R2_001.fastq
TAGGAC  614_S56_L001_R1_001.fastq 614_S56_L001_R2_001.fastq

ここでTipsを書く。Moving Pictures Tutorialでは気づかないことだが、Qiime2が受け入れてくれるfastqファイルにはあるルールがある。まずはgzファイルであること、そして必ず以下の名前のルールに従うことである。

  1. the sample identifier
  2. the barcode sequence or a barcode identifier
  3. the lane number
  4. the direction of the read (i.e. R1 or R2)
  5. the set number

要はL2S357_15_L001_R1_001.fastq.gzのような名前にしろということ。あとはペアがいないとダメとか色々あるが、この名前のルールを見落としがちなので、気をつけたい。基本NGSサービスなどから受け取ると名前はCasava 1.8 paired-end方式に従っている。なので、我々が変えるところは1. the sample identifierだけに留める。

なので、例えば11_15_L001_R1_001.fastq.gzというファイルを受け取り、その中に4個のサンプルがあるとすれば、

111_15_L001_R1_001.fastq.gz
112_15_L001_R1_001.fastq.gz
113_15_L001_R1_001.fastq.gz
114_15_L001_R1_001.fastq.gz

とかにすれば問題なく読み込まれる。長々と書いてしまったが、Sabreの続きをしてみる。タグ情報をtxtにまとめて、これを読み込ませる。

$ ./sabre pe -m 0 -f input_file1.fastq -r input_file2.fastq -b barcode_data.txt -u unknown_barcode_1.fastq -w unknown_barcode_2.fastq

-mで最大ミスマッチ数を決められる。-cを利用すれば、両側バーコードの場合も出来る。
当方の場合はパスを通してなくワーキングディレクトリがsabreなので、全てのファイルがsabreに吐き出される。バーコードのないやつはunknownとして吐き出される。

これをQiimeで解析したいフォルダに移す。適当にQiime2などというようなフォルダを作成する。そこにcasava-18-paired-end-demultiplexedフォルダを作成する。cdでQiime2フォルダに移動し、

conda activate qiime2-2019.10

をする。次にQiimeにデータをImportする。

Qiime2にデータを読み込ませる

Demultiplexedしたファイルなので、Qiime2のサイトにおけるCasava 1.8 paired-end demultiplexed fastqの項目を見るとやり方が載っている。

$ qiime tools import \
  --type 'SampleData[PairedEndSequencesWithQuality]' \
  --input-path casava-18-paired-end-demultiplexed \
  --input-format CasavaOneEightSingleLanePerSampleDirFmt \
  --output-path demux.qza

ちなみにデータの中身は

$ qiime demux summarize \
  --i-data demux.qza \
  --o-visualization demux.qzv
$ qiime tools view demux.qzv

で確認できた。

クオリティーコントロール

Moving Pictures Tutorialだとシングルエンドだったので、パラメータが少なかったが、ペアエンドだともう少しめんどくさい。Moving Pictures Tutorialと同様にDada2でやってみる。この時注意する点はWe need the reads to be long enough to overlap when joining paired endsということで、オーバーラップがきちんと作れる長さでトリムすること、後は先に述べたBiologicalな意味を考えることだと思う。後は基本--p-trim-left-f and --p-trim-left-r and for --p-trunc-len-f and --p-trunc-len-rを書いて終わり。前回同様、demux.qzvを参考にして決めた。

--p-trim-left-f x:各シーケンスのforward 配列の初めからx個目までを削除する.
--p-trim-left-r y:各シーケンスのReverse 配列の初めからy個目までを削除する.
--p-trunc-len-f m:各シーケンスのforward 配列のm個の塩基を採用する最大配列長を決定する. m個以降の塩基は除去する.
--p-trunc-len-r n:各シーケンスのreverse 配列のn個の塩基を採用する最大配列長を決定する. n個以降の配列は除去する.

(引用: https://qiita.com/kohei-108/items/547b2fbdf28fb04c28ea)

他にもオプションはあるので、余裕があれば使ってみたいと思う。

Parameters:
  --p-trunc-len-f INTEGER
                         Position at which forward read sequences should be
                         truncated due to decrease in quality. This truncates
                         the 3' end of the of the input sequences, which will
                         be the bases that were sequenced in the last cycles.
                         Reads that are shorter than this value will be
                         discarded. After this parameter is applied there must
                         still be at least a 20 nucleotide overlap between the
                         forward and reverse reads. If 0 is provided, no
                         truncation or length filtering will be performed
                                                                    [required]
  --p-trunc-len-r INTEGER
                         Position at which reverse read sequences should be
                         truncated due to decrease in quality. This truncates
                         the 3' end of the of the input sequences, which will
                         be the bases that were sequenced in the last cycles.
                         Reads that are shorter than this value will be
                         discarded. After this parameter is applied there must
                         still be at least a 20 nucleotide overlap between the
                         forward and reverse reads. If 0 is provided, no
                         truncation or length filtering will be performed
                                                                    [required]
  --p-trim-left-f INTEGER
                         Position at which forward read sequences should be
                         trimmed due to low quality. This trims the 5' end of
                         the input sequences, which will be the bases that
                         were sequenced in the first cycles.      [default: 0]
  --p-trim-left-r INTEGER
                         Position at which reverse read sequences should be
                         trimmed due to low quality. This trims the 5' end of
                         the input sequences, which will be the bases that
                         were sequenced in the first cycles.      [default: 0]
  --p-max-ee-f NUMBER    Forward reads with number of expected errors higher
                         than this value will be discarded.     [default: 2.0]
  --p-max-ee-r NUMBER    Reverse reads with number of expected errors higher
                         than this value will be discarded.     [default: 2.0]
  --p-trunc-q INTEGER    Reads are truncated at the first instance of a
                         quality score less than or equal to this value. If
                         the resulting read is then shorter than `trunc-len-f`
                         or `trunc-len-r` (depending on the direction of the
                         read) it is discarded.                   [default: 2]
  --p-chimera-method TEXT Choices('consensus', 'none', 'pooled')
                         The method used to remove chimeras. "none": No
                         chimera removal is performed. "pooled": All reads are
                         pooled prior to chimera detection. "consensus":
                         Chimeras are detected in samples individually, and
                         sequences found chimeric in a sufficient fraction of
                         samples are removed.           [default: 'consensus']
  --p-min-fold-parent-over-abundance NUMBER
                         The minimum abundance of potential parents of a
                         sequence being tested as chimeric, expressed as a
                         fold-change versus the abundance of the sequence
                         being tested. Values should be greater than or equal
                         to 1 (i.e. parents should be more abundant than the
                         sequence being tested). This parameter has no effect
                         if chimera-method is "none".           [default: 1.0]
  --p-n-threads INTEGER  The number of threads to use for multithreaded
                         processing. If 0 is provided, all available cores
                         will be used.                            [default: 1]
  --p-n-reads-learn INTEGER
                         The number of reads to use when training the error
                         model. Smaller numbers will result in a shorter run
                         time but a less reliable error model.
                                                            [default: 1000000]
  --p-hashed-feature-ids / --p-no-hashed-feature-ids
                         If true, the feature ids in the resulting table will
                         be presented as hashes of the sequences defining each
                         feature. The hash will always be the same for the
                         same sequence so this allows feature tables to be
                         merged across runs of this method. You should only
                         merge tables if the exact same parameters are used
                         for each run.                         [default: True]

実際に利用したコマンドは以下の通り。

qiime dada2 denoise-paired \
  --i-demultiplexed-seqs demux.qza \
  --p-trim-left-f 30 \
  --p-trim-left-r 20 \
  --p-trunc-len-f 270 \
  --p-trunc-len-r 280 \
  --o-table table.qza \
  --o-representative-sequences rep-seqs.qza \
  --o-denoising-stats denoising-stats.qza

Feature table and corresponding Feature sequencesの検討

ここからは前回と同じなので、コマンドだけ書く。

$ qiime feature-table summarize \
  --i-table table.qza \
  --o-visualization table.qzv \
  --m-sample-metadata-file sample-metadata.tsv

$ qiime feature-table tabulate-seqs \
  --i-data rep-seqs.qza \
  --o-visualization rep-seqs.qzv

データベースの変更

2019年12月28日現在、Silva138のQiime版はまだ公開されてないようなので、https://www.arb-silva.de/download/archive/qiime からSilva132verを落とした。僕の場合はhttps://docs.qiime2.org/2019.10/tutorials/feature-classifier/ を参考に自分でデータベースを作ってみたが、多分data-resourcesから落とした方が早いし、楽だと思う。99%OTUのやつを利用しておけば間違いない。

後は基本的にMoving Pictures Tutorialと同じで、データベースのところだけ変える。

$ qiime feature-classifier classify-sklearn \
  --i-classifier silva-132-99-nb-classifier.qza.qza \
  --i-reads rep-seqs.qza \
  --o-classification taxonomy.qza

$ qiime metadata tabulate \
  --m-input-file taxonomy.qza \
  --o-visualization taxonomy.qzv

とする。棒グラフにする場合は、

$ qiime taxa barplot \
  --i-table table.qza \
  --i-taxonomy taxonomy.qza \
  --m-metadata-file sample-metadata.tsv \
  --o-visualization taxa-bar-plots.qzv

$ qiime tools view taxa-bar-plots.qzv

まとめ

・データベースは自分で作れるし、Qiime2のサイトからも落とせる。
・基本的にMoving Pictures Tutorialと同じ。
・ペアエンドとシーケンス領域の特性に気をつける。

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
ユーザーは見つかりませんでした