Pastd image 20230525175303.png

1、質(zhì)量控制

1.使用fastpc查看原始數(shù)據(jù)質(zhì)量

fastqc \
    -f fastq \
    B1_raw_L001_R1_001.fastq  B1_raw_L001_R2_001.fastq  J1_raw_L001_R1_001.fastq  J1_raw_L001_R2_001.fastq  Y1_raw_L001_R1_001.fastq  Y1_raw_L001_R2_001.fastq \
    -o  0_fastqc

B1 R1 質(zhì)量

Pasted image 20230526100111.png

Pasted image 20230525184714.png

從圖中可見，堿基位置于260-282，箱型圖上下虛線（10%-90%），在這個(gè)位置，有多數(shù)堿基質(zhì)量低于20。

Pasted image 20230529154729.png

從圖中可見，無P5/P7接頭序列。

B1 R2 質(zhì)量

[圖片上傳中...(Pasted image 20230526100032.png-a7db73-1685585942439-0)]

Pasted image 20230526100032.png

觀察質(zhì)量圖，于210-240堿基位置，觀察黃色箱型圖（25%-75%），發(fā)現(xiàn)大部分堿基質(zhì)量低于20。

Pasted image 20230529154812.png

從圖中可見，無P5/P7接頭序列。

2.使用cutadapt清除測序引物和P5/P7接頭

Pasted image 20230529154417.png

• 測序區(qū)域?yàn)閂3-V4區(qū)引物序列如下：

341F引物為：5′-CCTACGGGNGGCWGCAG-3′；
805R引物為：5′-GACTACHVGGGTATCTAATCC-3′

• P5和P7接頭序列如下：

F：AATGATACGGCGACCACCGAGATCT
R：CAAGCAGAAGACGGCATACGAGAT

cutadapt -g AATGATACGGCGACCACCGAGATCT -g CCTACGGGNGGCWGCAG -G CAAGCAGAAGACGGCATACGAGAT -G GACTACHVGGGTATCTAATCC -j 0 -e 0.2 -O 5 --trim-n -o B1_cut_R1_.fastq -p B1_cut_R2_.fastq ../1_RAW/B1_raw_L001_R1_001.fastq ../1_RAW/B1_raw_L001_R2_001.fastq

cutadapt -g AATGATACGGCGACCACCGAGATCT -g CCTACGGGNGGCWGCAG -G CAAGCAGAAGACGGCATACGAGAT -G GACTACHVGGGTATCTAATCC -j 0 -e 0.2 -O 5 --trim-n -o J1_cut_R1_.fastq -p J1_cut_R2_.fastq ../1_RAW/J1_raw_L001_R1_001.fastq ../1_RAW/J1_raw_L001_R2_001.fastq

cutadapt -g AATGATACGGCGACCACCGAGATCT -g CCTACGGGNGGCWGCAG -G CAAGCAGAAGACGGCATACGAGAT -G GACTACHVGGGTATCTAATCC -j 0 -e 0.2 -O 5 --trim-n -o Y1_cut_R1_.fastq -p Y1_cut_R2_.fastq ../1_RAW/Y1_raw_L001_R1_001.fastq ../1_RAW/Y1_raw_L001_R2_001.fastq

3.使用fastp剪切3’端質(zhì)量不好的堿基

fastp \
    --in1 B1_cut_R1_.fastq  --in2 B1_cut_R2_.fastq \
    --out1 B1_p_clean_R1.fastq --out2 B1_p_clean_R2.fastq \
    --failed_out B1_p_failed_out.fastq \
    -A -3 -G -W 4 -q 20 -h B1_p_fastp.html -j B1_p_fastp.json -R B1_p_report -w 16

fastp \
    --in1 J1_cut_R1_.fastq  --in2 J1_cut_R2_.fastq \
    --out1 J1_p_clean_R1.fastq --out2 J1_p_clean_R2.fastq \
    --failed_out J1_p_failed_out.fastq \
    -A -3 -G -W 4 -q 20 -h J1_p_fastp.html -j J1_p_fastp.json -R J1_p_report -w 16

fastp \
    --in1 Y1_cut_R1_.fastq  --in2 Y1_cut_R2_.fastq \
    --out1 Y1_p_clean_R1.fastq --out2 Y1_p_clean_R2.fastq \
    --failed_out Y1_p_failed_out.fastq \
    -A -3 -G -W 4 -q 20 -h Y1_p_fastp.html -j Y1_p_fastp.json -R Y1_p_report -w 16

• 選項(xiàng)詳情
  -A, --disable_adapter_trimming
  默認(rèn)情況下啟用適配器修剪。如果指定了此選項(xiàng)，則禁用適配器修剪
  -G, --disable_trim_poly_g
  禁用polyG尾部修剪，默認(rèn)情況下會自動(dòng)為Illumina NextSeq/NovaSeq數(shù)據(jù)啟用修剪
  -3, --cut_tail
  將滑動(dòng)窗口從尾部（3'）向前移動(dòng)，如果其平均質(zhì)量<閾值，則將底部放入窗口中，否則停止。
  -W, --cut_window_size
  由cut_front、cut_tail或cut_sliding共享的窗口大小選項(xiàng)。范圍：1~1000，默認(rèn)值：4（int[=4]）
  -q, --qualified_quality_phred
  基準(zhǔn)合格的質(zhì)量值。默認(rèn)值15表示默認(rèn)質(zhì)量>=Q15合格。（int[=15]）
  -u, --unqualified_percent_limit
  允許多少百分比的堿基不合格（0~100）。默認(rèn)值40表示40%（int[=40]）
  -n, --n_base_limit
  如果一個(gè)讀的N基數(shù)大于N_base_limit，則丟棄該讀/對。默認(rèn)值為5（int[=5]）
  -l, --length_required
  短于length_required的讀取將被丟棄，默認(rèn)值為15。（整數(shù)[=15]）
  --filter_by_index_threshold
  索引篩選允許的索引條形碼差異，默認(rèn)為0表示完全相同。（int[=0]）
  -w, --thread
  工作線程號，默認(rèn)值為3（int[=3]）

4.使用fastqc查看質(zhì)控后的質(zhì)量

fastqc \
    -f fastq \
    B1_p_clean_R1.fastq  B1_p_clean_R2.fastq  J1_p_clean_R1.fastq  J1_p_clean_R2.fastq  Y1_p_clean_R1.fastq  Y1_p_clean_R2.fastq \
    -o  0_fastqc

B1 R1：

Pasted image 20230531160154.png

Pasted image 20230531160257.png

Pasted image 20230531160324.png

• 從B1 R1的fastqc結(jié)果來看，3’末端的質(zhì)量低于20的堿基去除的差不多了。序列長度大多數(shù)在265-270bp。

B1 R2：

Pasted image 20230531160537.png

Pasted image 20230531160554.png

Pasted image 20230531160618.png

• 從B1 R2的fastqc結(jié)果來看，3’末端的質(zhì)量低于20的堿基去除的差不多了。序列長度大多數(shù)在215-220bp。

2、導(dǎo)入qiime2

1.manifest文件

• 自己列出

2.sample-metadata文件

• 自己列出

2、原始數(shù)據(jù)導(dǎo)入qiime2

1、qiime tools import導(dǎo)入

time qiime tools import \
    --type 'SampleData[PairedEndSequencesWithQuality]' \
    --input-path manifest \
    --output-path Paired-end.qza \
    --input-format PairedEndFastqManifestPhred33V2

可視化

time qiime demux summarize \
    --i-data Paired-end.qza \
    --o-visualization Paired-end.qzv

2、dada2降噪聚類

• --p-trunc-len-f和--p-trunc-len-r兩個(gè)選項(xiàng)規(guī)定了正反兩個(gè)序列的reads的長度，低于這個(gè)長度的reads將被丟棄，高于這個(gè)長度的reads將被剪切至這一長度。

• 對于雙末端序列，設(shè)置的太小，中間的overlap低于dada2默認(rèn)的12bp，將無法合并。設(shè)置的太大，大多數(shù)沒那么長的有用reads將被丟棄。

• 所以這兩個(gè)選項(xiàng)需要結(jié)合雙端序列的質(zhì)量和長度分布來看，可以查看上一步的Paired-end.qzv文件，只能看到質(zhì)量分布，但看不到read的長度分布。在fastp質(zhì)控后再經(jīng)過fastqc，可以直觀的看清楚R1和R2的長度分布。

• 本次測序中的三個(gè)樣品的R1的長度大約都為265-270bp，R2的長度大約都為215-220bp。這里選取的參數(shù)為：

--p-trunc-len-f 260
--p-trunc-len-r 206

time qiime dada2 denoise-paired \
    --i-demultiplexed-seqs Paired-end.qza \
    --p-n-threads 0 \
    --p-trunc-len-f 260 \
    --p-trunc-len-r 206 \
    --o-table table.qza \
    --o-representative-sequences rep-seqs.qza \
    --o-denoising-stats denoising-stats.qza 

• 選項(xiàng)詳情
  --p-max-ee-f
  預(yù)期錯(cuò)誤數(shù)高于此值的正向讀取將被丟棄。[default: 2.0]
  --p-max-ee-r
  預(yù)期錯(cuò)誤數(shù)高于此值的反向讀取將被丟棄。[default: 2.0]
  --p-trunc-q
  在質(zhì)量分?jǐn)?shù)小于或等于此值的第一個(gè)實(shí)例處，讀取被截?cái)唷Ｈ绻玫降淖x取比“trunc-len-f”或“trunc-len-r”短（取決于讀取的方向），則丟棄它。[default: 2]
  --p-min-overlap
  合并正向和反向讀取所需的最小重疊長度。[default: 12]
  --p-pooling-method
  用于對樣本進(jìn)行池化去噪的方法?！癷ndependent”：樣本被獨(dú)立地去噪?！皃seudo”：偽池化方法用于近似樣本池化。簡言之，樣本被獨(dú)立地去噪一次，在至少2個(gè)樣本中檢測到的ASV被記錄，樣本被單獨(dú)地去噪第二次，但這一次是在事先知道記錄的ASV的情況下，從而對這些ASV具有更高的靈敏度。 [default: 'independent']
  --p-chimera-method
  去除嵌合體的方法?！盁o”：未進(jìn)行嵌合體去除?！昂喜ⅰ保核凶x數(shù)在嵌合體檢測之前合并?！耙恢滦浴保涸跇颖局袉为?dú)檢測到嵌合體，并去除足夠部分樣本中發(fā)現(xiàn)的嵌合體序列。[default: 'consensus']
  --p-min-fold-parent-over-abundance
  被測試為嵌合序列的潛在親本的最小豐度，表示為與被測試序列豐度的倍數(shù)變化。值應(yīng)該大于或等于1（即父母應(yīng)該比被測試的序列更豐富）。如果嵌合體方法為“無”，則此參數(shù)無效。[default: 1.0]
  --p-n-threads
  用于多線程處理的線程數(shù)。如果提供了0，則將使用所有可用的核心。[default: 1]
  --p-n-reads-learn
  訓(xùn)練錯(cuò)誤模型時(shí)要使用的讀取次數(shù)。較小的數(shù)字將導(dǎo)致較短的運(yùn)行時(shí)間，但誤差模型的可靠性較低。[default: 1000000]
  --p-hashed-feature-ids / --p-no-hashed-feature-ids
  如果為true，則結(jié)果表中的特征id將顯示為定義每個(gè)特征的序列的散列。對于相同的序列，散列總是相同的，因此這允許在運(yùn)行該方法時(shí)合并功能表。只有在每次運(yùn)行使用完全相同的參數(shù)時(shí)，才應(yīng)合并表。[default: True]

可視化

time qiime metadata tabulate \
    --m-input-file denoising-stats.qza \
    --o-visualization denoising-stats.qzv

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

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

denoising-stats.qzv：

Pasted image 20230531163103.png

經(jīng)過dada2聚類后的序列個(gè)數(shù)占原始序列的50-30%，為正常水平。

4、ASV物種注釋

• 從qiime2官網(wǎng)下載訓(xùn)練好的Greengenes葉貝斯全長分類器
  Greengenes2 2022.10 full length sequences

臨時(shí)文件夾/tmp沒空間了，改一下
export TMPDIR="/data2020/zzz/Tmpdir"

time qiime feature-classifier classify-sklearn \
    --i-classifier 0_classify_sklearn/2022.10.backbone.full-length.nb.qza \
    --p-n-jobs -1 \
    --i-reads rep-seqs.qza \
    --o-classification 3_ASVtaxonomy/taxonomy.qza

• 選項(xiàng)詳情
  --p-reads-per-batch
  每個(gè)批次中要處理的讀取數(shù)。如果為“自動(dòng)”，則此參數(shù)將自動(dòng)縮放為最小值（查詢序列數(shù)/n-作業(yè)數(shù)，20000）。[default: 'auto']
  --p-n-jobs
  并發(fā)工作進(jìn)程的最大數(shù)目。如果-1，則使用所有CPU。如果給定1，則根本不使用并行計(jì)算代碼，這對于調(diào)試非常有用。對于-1以下的n個(gè)作業(yè)，使用（n_cpus+1+n個(gè)作業(yè)）。因此，對于n-jobs=-2，將使用除一個(gè)CPU以外的所有CPU。[default: 1]
  --p-pre-dispatch
  “all”或表達(dá)式，如在“3 n_jobs”中。要預(yù)先調(diào)度的（任務(wù)的）批數(shù)。[default: '2 n_jobs']
  --p-confidence
  限制分類深度的置信閾值。設(shè)置為“disable”可禁用置信度計(jì)算，或設(shè)置為0可計(jì)算置信度，但不應(yīng)用它來限制賦值的分類深度。[default: 0.7]
  --p-read-orientation
  相對于參考序列的讀取方向。相同將導(dǎo)致讀數(shù)分類不變；反向補(bǔ)碼將導(dǎo)致在分類之前對讀數(shù)進(jìn)行反向和補(bǔ)碼?！癮uto”將根據(jù)前100次讀取的置信度估計(jì)自動(dòng)檢測方向。[default: 'auto']

** 注：聚類出的ASV默認(rèn)為99%OTU聚類。**

可視化

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

柱狀堆疊圖

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

合并屬出表

time qiime taxa collapse \
    --i-table ../table.qza \
    --i-taxonomy taxonomy.qza \
    --p-level 6 \
    --o-collapsed-table collapsed-table-l6.qza

time qiime tools export \
    --input-path collapsed-table-l6.qza \
    --output-path collapsed-table-l6

biom convert \
    -i feature-table.biom \
    -o 6-taxonomy.tsv \
    --to-tsv

合并門出表

time qiime taxa collapse \
    --i-table ../table.qza \
    --i-taxonomy taxonomy.qza \
    --p-level 2 \
    --o-collapsed-table collapsed-table-l2.qza

time qiime tools export \
    --input-path collapsed-table-l2.qza \
    --output-path collapsed-table-l2

biom convert \
    -i feature-table.biom \
    -o 2-taxonomy.tsv \
    --to-tsv