Pastd image 20230525175303.png

1、质量控制

1.使用fastpc查看原始数据质量

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 质量

Pasted image 20230526100111.png

Pasted image 20230525184714.png

从图中可见，碱基位置于260-282，箱型图上下虚线（10%-90%），在这个位置，有多数碱基质量低于20。

Pasted image 20230529154729.png

从图中可见，无P5/P7接头序列。

B1 R2 质量

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

Pasted image 20230526100032.png

观察质量图，于210-240碱基位置，观察黄色箱型图（25%-75%），发现大部分碱基质量低于20。

Pasted image 20230529154812.png

从图中可见，无P5/P7接头序列。

2.使用cutadapt清除测序引物和P5/P7接头

Pasted image 20230529154417.png

• 测序区域为V3-V4区引物序列如下：

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’端质量不好的碱基

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

• 选项详情
  -A, --disable_adapter_trimming
  默认情况下启用适配器修剪。如果指定了此选项，则禁用适配器修剪
  -G, --disable_trim_poly_g
  禁用polyG尾部修剪，默认情况下会自动为Illumina NextSeq/NovaSeq数据启用修剪
  -3, --cut_tail
  将滑动窗口从尾部（3'）向前移动，如果其平均质量<阈值，则将底部放入窗口中，否则停止。
  -W, --cut_window_size
  由cut_front、cut_tail或cut_sliding共享的窗口大小选项。范围：1~1000，默认值：4（int[=4]）
  -q, --qualified_quality_phred
  基准合格的质量值。默认值15表示默认质量>=Q15合格。（int[=15]）
  -u, --unqualified_percent_limit
  允许多少百分比的碱基不合格（0~100）。默认值40表示40%（int[=40]）
  -n, --n_base_limit
  如果一个读的N基数大于N_base_limit，则丢弃该读/对。默认值为5（int[=5]）
  -l, --length_required
  短于length_required的读取将被丢弃，默认值为15。（整数[=15]）
  --filter_by_index_threshold
  索引筛选允许的索引条形码差异，默认为0表示完全相同。（int[=0]）
  -w, --thread
  工作线程号，默认值为3（int[=3]）

4.使用fastqc查看质控后的质量

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结果来看，3’末端的质量低于20的碱基去除的差不多了。序列长度大多数在265-270bp。

B1 R2：

Pasted image 20230531160537.png

Pasted image 20230531160554.png

Pasted image 20230531160618.png

• 从B1 R2的fastqc结果来看，3’末端的质量低于20的碱基去除的差不多了。序列长度大多数在215-220bp。

2、导入qiime2

1.manifest文件

• 自己列出

2.sample-metadata文件

• 自己列出

2、原始数据导入qiime2

1、qiime tools import导入

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两个选项规定了正反两个序列的reads的长度，低于这个长度的reads将被丢弃，高于这个长度的reads将被剪切至这一长度。

• 对于双末端序列，设置的太小，中间的overlap低于dada2默认的12bp，将无法合并。设置的太大，大多数没那么长的有用reads将被丢弃。

• 所以这两个选项需要结合双端序列的质量和长度分布来看，可以查看上一步的Paired-end.qzv文件，只能看到质量分布，但看不到read的长度分布。在fastp质控后再经过fastqc，可以直观的看清楚R1和R2的长度分布。

• 本次测序中的三个样品的R1的长度大约都为265-270bp，R2的长度大约都为215-220bp。这里选取的参数为：

--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 

• 选项详情
  --p-max-ee-f
  预期错误数高于此值的正向读取将被丢弃。[default: 2.0]
  --p-max-ee-r
  预期错误数高于此值的反向读取将被丢弃。[default: 2.0]
  --p-trunc-q
  在质量分数小于或等于此值的第一个实例处，读取被截断。如果得到的读取比“trunc-len-f”或“trunc-len-r”短（取决于读取的方向），则丢弃它。[default: 2]
  --p-min-overlap
  合并正向和反向读取所需的最小重叠长度。[default: 12]
  --p-pooling-method
  用于对样本进行池化去噪的方法。“independent”：样本被独立地去噪。“pseudo”：伪池化方法用于近似样本池化。简言之，样本被独立地去噪一次，在至少2个样本中检测到的ASV被记录，样本被单独地去噪第二次，但这一次是在事先知道记录的ASV的情况下，从而对这些ASV具有更高的灵敏度。 [default: 'independent']
  --p-chimera-method
  去除嵌合体的方法。“无”：未进行嵌合体去除。“合并”：所有读数在嵌合体检测之前合并。“一致性”：在样本中单独检测到嵌合体，并去除足够部分样本中发现的嵌合体序列。[default: 'consensus']
  --p-min-fold-parent-over-abundance
  被测试为嵌合序列的潜在亲本的最小丰度，表示为与被测试序列丰度的倍数变化。值应该大于或等于1（即父母应该比被测试的序列更丰富）。如果嵌合体方法为“无”，则此参数无效。[default: 1.0]
  --p-n-threads
  用于多线程处理的线程数。如果提供了0，则将使用所有可用的核心。[default: 1]
  --p-n-reads-learn
  训练错误模型时要使用的读取次数。较小的数字将导致较短的运行时间，但误差模型的可靠性较低。[default: 1000000]
  --p-hashed-feature-ids / --p-no-hashed-feature-ids
  如果为true，则结果表中的特征id将显示为定义每个特征的序列的散列。对于相同的序列，散列总是相同的，因此这允许在运行该方法时合并功能表。只有在每次运行使用完全相同的参数时，才应合并表。[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

经过dada2聚类后的序列个数占原始序列的50-30%，为正常水平。

4、ASV物种注释

• 从qiime2官网下载训练好的Greengenes叶贝斯全长分类器
  Greengenes2 2022.10 full length sequences

临时文件夹/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

• 选项详情
  --p-reads-per-batch
  每个批次中要处理的读取数。如果为“自动”，则此参数将自动缩放为最小值（查询序列数/n-作业数，20000）。[default: 'auto']
  --p-n-jobs
  并发工作进程的最大数目。如果-1，则使用所有CPU。如果给定1，则根本不使用并行计算代码，这对于调试非常有用。对于-1以下的n个作业，使用（n_cpus+1+n个作业）。因此，对于n-jobs=-2，将使用除一个CPU以外的所有CPU。[default: 1]
  --p-pre-dispatch
  “all”或表达式，如在“3 n_jobs”中。要预先调度的（任务的）批数。[default: '2 n_jobs']
  --p-confidence
  限制分类深度的置信阈值。设置为“disable”可禁用置信度计算，或设置为0可计算置信度，但不应用它来限制赋值的分类深度。[default: 0.7]
  --p-read-orientation
  相对于参考序列的读取方向。相同将导致读数分类不变；反向补码将导致在分类之前对读数进行反向和补码。“auto”将根据前100次读取的置信度估计自动检测方向。[default: 'auto']

** 注：聚类出的ASV默认为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