LoginSignup
3
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 1 year has passed since last update.

沖縄高専Advent Calendar 2022

Day 6
@Yokohide(横のよこ)

Qiime2でPaired-Endクオリティーが悪い時の応急処置

Last updated at Posted at 2022-12-05

概要

こちらで話されていることを試してみました。

簡単に言うと、「Paired-Endのシーケンスデータで、クオリティーが低くく、
マージが上手く行かない際は、片方だけでTaxonomyしてみようぜ(非推奨だけど)」
ということです。

最近、マージが上手くいかないデータを解析した時にこのやり取りを参考にして、
「こういう発想があったのか」と思ったので試してみました。

こちらのデータを使用させて頂きました。
https://www.ncbi.nlm.nih.gov/sra/SRX18463592[accn]

0.準備

bash
# sra-toolsインストールしておく
brew install sratoolkit

# オリジナルのデータを保存しておく場所
mkdir data
cd data/

# SRAを取ってきて、fastqにするまでしてもらう。
fasterq-dump -p SRR22498300

ls
> SRR22498300_1.fastq SRR22498300_2.fastq

# 圧縮しておく
for file in `ls *.fastq`; do
pigz $file
done

1.まずはふつうに

bash
mkdir paired
cd paired

mkdir raw-data
cd raw-data/
ln -s ../../data/SRR22498300_1.fastq.gz ./sample1_S1_L001_R1_001.fastq.gz
ln -s ../../data/SRR22498300_2.fastq.gz ./sample1_S1_L001_R2_001.fastq.gz

# import
cd ..
qiime tools import \
    --type "SampleData[PairedEndSequencesWithQuality]" \
    --input-path ./raw-data \
    --input-format CasavaOneEightSingleLanePerSampleDirFmt \
    --output-path demux-paired-end.qza

# 可視化
qiime demux summarize --i-data demux-paired-end.qza --o-visualization demux-paired-end.qzv
qiime tools view demux-paired-end.qzv

Screenshot 2022-12-04 at 18.54.59.png

Forwardは300、Reverseは285あたりで切れば良さそうです。

bash
qiime dada2 denoise-paired \
    --i-demultiplexed-seqs ./demux-paired-end.qza \
    --p-trunc-len-f 300 \
    --p-trunc-len-r 285 \
    --o-table table-dada2.qza \
    --o-representative-sequences rep-seqs-dada2.qza \
    --o-denoising-stats stats-dada2.qza \
    --p-n-threads 16

# 可視化
qiime metadata tabulate \
    --m-input-file stats-dada2.qza \
    --o-visualization stats-dada2.qzv
qiime feature-table summarize \
    --i-table table-dada2.qza \
    --o-visualization table-dada2.qzv
qiime feature-table tabulate-seqs \
    --i-data rep-seqs-dada2.qza \
    --o-visualization rep-seqs-dada2.qzv

# 代表配列を見てみる
qiime tools view rep-seqs-dada2.qzv

Screen Shot 2022-12-05 at 9.30.57.png
平均453bpで、310配列配列取れているようです。

2. エラーを発生

mkdir -p error/raw-data

pigz -d data/SRR22498300_2.fastq.gz -c > data/SRR22498300_2.fastq

Reverseをだんだんと悪くする。

reverse_shift_quality.py
from Bio import SeqIO
import random

input_fastq = "./data/SRR22498300_2.fastq"
output_fastq = "./error/raw-data/sample1_S1_L001_R2_001.fastq"
records = SeqIO.parse(input_fastq, 'fastq')

new_records = []
# https://biotech-lab.org/articles/3133#SeqRecord
for record in records:
    quality = record.letter_annotations['phred_quality']
    changed = []
    # ここらへんの数字は適当
    top = 32
    bottom = 20
    for i, q in enumerate(quality):
        if i % 20 == 0:
            if top > 15:
                top += -1
            if bottom > 3:
                bottom += -2
        q = random.randint(bottom, top)
        changed.append(q)
    record.letter_annotations = {'phred_quality': changed}
    new_records.append(record)

SeqIO.write(new_records, output_fastq, "fastq")

Forwardは3'側だけいじる

Forward_shift_quality.py
from Bio import SeqIO
import random

forward_input = "./data/SRR22498300_1.fastq"
forward_output = "./error/raw-data/sample1_S1_L001_R1_001.fastq"
records = SeqIO.parse(forward_input, 'fastq')

new_records = []
# https://biotech-lab.org/articles/3133#SeqRecord
for record in records:
    quality = record.letter_annotations['phred_quality']
    # 270辺りから15~35でランダム発生させる。
    changed = []
    for i, q in enumerate(quality):
        if i > 270:
            q = random.randint(15, 35)
        changed.append(q)
    record.letter_annotations = {'phred_quality': changed}
    new_records.append(record)

SeqIO.write(new_records, forward_output, "fastq")

error/raw-data/ 下に、sample1_S1_L001_R2_001.fastqとsample1_S1_L001_R1_001.fastqが生成されるので、準備

bash
cd error/raw-data/
# 圧縮
pigz sample1_S1_L001_R2_001.fastq
pigz sample1_S1_L001_R1_001.fastq

まずは普通にPaired-Endとしてimportしてみます。

bash
# import
cd ..
qiime tools import \
    --type "SampleData[PairedEndSequencesWithQuality]" \
    --input-path ./raw-data \
    --input-format CasavaOneEightSingleLanePerSampleDirFmt \
    --output-path demux-paired-end.qza

# 可視化
qiime demux summarize --i-data demux-paired-end.qza --o-visualization demux-paired-end.qzv
qiime tools view demux-paired-end.qzv

Screen Shot 2022-12-05 at 9.24.55.png

ランダム発生の仕方が雑なので不自然すぎますが、いい感じに使いづらいデータにすることはできました。

Forwardは271、Reverseは20で切ってみます。

bash
qiime dada2 denoise-paired \
    --i-demultiplexed-seqs ./demux-paired-end.qza \
    --p-trunc-len-f 271 \
    --p-trunc-len-r 20 \
    --o-table table-dada2.qza \
    --o-representative-sequences rep-seqs-dada2.qza \
    --o-denoising-stats stats-dada2.qza \
    --p-n-threads 16

# 可視化
qiime metadata tabulate \
    --m-input-file stats-dada2.qza \
    --o-visualization stats-dada2.qzv
qiime feature-table summarize \
    --i-table table-dada2.qza \
    --o-visualization table-dada2.qzv
qiime feature-table tabulate-seqs \
    --i-data rep-seqs-dada2.qza \
    --o-visualization rep-seqs-dada2.qzv

# 代表配列を見てみる
qiime tools view rep-seqs-dada2.qzv

270bpで、1つしか出なかったようです。
（mergeが1つもできず、エラー発生することを想定したいたのですが...。）
Screen Shot 2022-12-05 at 10.06.56.png

では、本題のSingle-Endとして扱ってみます。

bash
# 避難させる
mkdir error-paired
mv *.qza *.qzv error-paired/
mv raw-data/sample1_S1_L001_R2_001.fastq.gz error-paired/

一旦、指定可能なタイプを表示します。

bash
qiime tools import --show-importable-types
    Bowtie2Index
    DeblurStats
    DistanceMatrix
    EMPPairedEndSequences
    EMPSingleEndSequences
    ErrorCorrectionDetails
    FeatureData[AlignedProteinSequence]
    FeatureData[AlignedRNASequence]
    FeatureData[AlignedSequence]
    FeatureData[BLAST6]
    FeatureData[Differential]
    FeatureData[Importance]
    FeatureData[PairedEndRNASequence]
    FeatureData[PairedEndSequence]
    FeatureData[ProteinSequence]
    FeatureData[RNASequence]
    FeatureData[Sequence]
    FeatureData[Taxonomy]
    FeatureTable[Balance]
    FeatureTable[Composition]
    FeatureTable[Design]
    FeatureTable[Frequency]
    FeatureTable[PercentileNormalized]
    FeatureTable[PresenceAbsence]
    FeatureTable[RelativeFrequency]
    Hierarchy
    MultiplexedPairedEndBarcodeInSequence
    MultiplexedSingleEndBarcodeInSequence
    PCoAResults
    Phylogeny[Rooted]
    Phylogeny[Unrooted]
    Placements
    ProcrustesStatistics
    QualityFilterStats
    RawSequences
    SampleData[AlphaDiversity]
    SampleData[ArtificialGrouping]
    SampleData[BooleanSeries]
    SampleData[ClassifierPredictions]
    SampleData[DADA2Stats]
    SampleData[FirstDifferences]
    SampleData[JoinedSequencesWithQuality]
    SampleData[PairedEndSequencesWithQuality]
    SampleData[Probabilities]
    SampleData[RegressorPredictions]
    SampleData[SequencesWithQuality]
    SampleData[Sequences]
    SampleData[TrueTargets]
    SampleEstimator[Classifier]
    SampleEstimator[Regressor]
    SeppReferenceDatabase
    TaxonomicClassifier
    UchimeStats

Paired-Endで、SampleData[PairedEndSequencesWithQuality]としてimportしたので、
今回はSampleData[SequencesWithQuality]として扱います。

bash
qiime tools import \
    --type "SampleData[SequencesWithQuality]" \
    --input-path ./raw-data \
    --input-format CasavaOneEightSingleLanePerSampleDirFmt \
    --output-path demux-single-end.qza

前回と同じく、271で切ってみます。

bash
qiime dada2 denoise-single \
    --i-demultiplexed-seqs ./demux-single-end.qza \
    --p-trunc-len 271 \
    --o-table table-dada2.qza \
    --o-representative-sequences rep-seqs-dada2.qza \
    --o-denoising-stats stats-dada2.qza \
    --p-n-threads 16

# 可視化
qiime metadata tabulate \
    --m-input-file stats-dada2.qza \
    --o-visualization stats-dada2.qzv
qiime feature-table summarize \
    --i-table table-dada2.qza \
    --o-visualization table-dada2.qzv
qiime feature-table tabulate-seqs \
    --i-data rep-seqs-dada2.qza \
    --o-visualization rep-seqs-dada2.qzv

# 代表配列を見てみる
qiime tools view rep-seqs-dada2.qzv

271bpの305配列取得することができました。
Screen Shot 2022-12-05 at 10.22.28.png

まとめ

今回はキレイな配列データの「クオリティのみ」をいじったので、そりゃあ上手くいくわという感じですが、
mergeする作業によって代表配列が作れない場合の応急処置としてのTipsでした。

もしこのようなデータにふれる機会があれば、実際うまくいくのか試してみたいと思います。

参考

3
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
3
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?