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


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

https://forum.qiime2.org/t/recommendations-for-processing-with-poor-quality-forward-reads/20243

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

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

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

# 0.準備
```bash: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: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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/550844/9f51765b-98d2-e3d5-c30a-e43d2198fcaa.png)

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

```bash: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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/550844/6e76f8e7-8a2a-7580-aced-486d0b38956f.png)
平均453bpで、310配列配列取れているようです。



# 2. エラーを発生

```
mkdir -p error/raw-data

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

Reverseをだんだんと悪くする。
```python: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'側だけいじる
```python: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:bash
cd error/raw-data/
# 圧縮
pigz sample1_S1_L001_R2_001.fastq
pigz sample1_S1_L001_R1_001.fastq
```

まずは普通にPaired-Endとしてimportしてみます。
```bash: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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/550844/2d020c90-f063-9102-37a6-0e5660c6c0b4.png)


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

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

```bash: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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/550844/cb26aca6-de80-ba27-210b-2fa94cbfb635.png)



では、本題のSingle-Endとして扱ってみます。
```bash:bash
# 避難させる
mkdir error-paired
mv *.qza *.qzv error-paired/
mv raw-data/sample1_S1_L001_R2_001.fastq.gz error-paired/
```
一旦、指定可能なタイプを表示します。

```bash: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:bash
qiime tools import \
    --type "SampleData[SequencesWithQuality]" \
    --input-path ./raw-data \
    --input-format CasavaOneEightSingleLanePerSampleDirFmt \
    --output-path demux-single-end.qza
```



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

```bash: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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/550844/3965a129-9c79-0d81-77d6-a81efe391b33.png)

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

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

# 参考
https://forum.qiime2.org/t/recommendations-for-processing-with-poor-quality-forward-reads/20243

https://biotech-lab.org/articles/3133#SeqRecord

https://www.ncbi.nlm.nih.gov/sra/SRX18463592[accn]