16s RNA ：16s rDNA 基因存在于所有細菌的基因組中，具有高度的保守性橘沥。 該序列包含9個高變區(qū)和10個保守區(qū)盖喷，通過對某一段高變區(qū)序列進行PCR擴增后進行測序炮沐，得到1500bp左右的序列。
宏基因組：是將微生物基因組DNA隨機打斷成500bp的小片段峰弹，然后在片段兩端加入通用引物進行PCR擴增測序店量，再通過組裝的方式，將小片段拼接成較長的序列鞠呈。
針對16s RNA 數(shù)據(jù)主要使用 QIIME2 軟件進行融师。
此分析筆記流程概覽
——8 以后的分析待后續(xù)更新吧！
數(shù)據(jù)下載：ENA 瀏覽器 (ebi.ac.uk)
0.數(shù)據(jù)下載 PRJNA359624
點擊獲取下載腳本粟按，到服務器上運行
nohup ./ena-file-download-read_run-PRJNA359624-fastq_ftp-20240111-0625.sh &
確定下載成功诬滩，一共52個樣本
分析流程如下：
要自己安裝fastqc哦
1.使用fastqc進行質控
mkdir fastqout
fastqc -c *.fastq.gz -o fastqout/
看見報告正常生成
multiqc 整合質控報告
multiqc fastqout/
查看質控報告霹粥，在當前目錄下生成，
質量數(shù)據(jù)沒啥問題疼鸟，下一步
2.切除引物
- 送公司測序返還的數(shù)據(jù)一般都是拆分過并去除了引物的后控，可以自己再做一下質檢，引物切除了么空镜？浩淘？
發(fā)現(xiàn)此處我們下載的已經(jīng)切除，不用再去吴攒。
3. 數(shù)據(jù)導入：在qiime2中分析測序數(shù)據(jù)
3.0 安裝 qiime2
參考官網(wǎng) 本機安裝 QIIME 2 — QIIME 2 2023.9.2 文檔
wget https://data.qiime2.org/distro/amplicon/qiime2-amplicon-2023.9-py38-linux-conda.yml
conda env create -n qiime2-amplicon-2023.9 --file qiime2-amplicon-2023.9-py38-linux-conda.yml
conda activate qiime2-amplicon-2023.9
確定安裝成功张抄，我大概花了20分鐘吧~
3.1 先查看下測序類型質量值體系
質量值體系分為 Phred33 和 Phred 64兩種，如下圖所示洼怔，一般看fastq文件的質量值那行包含署惯！和？（對應ASCII值33和63）等镣隶，即為Phred33體系（一般都為Phred33）极谊。
發(fā)現(xiàn)全為?
3.2 寫參數(shù)配置文件
下面為了運行 qiime2 需要自己寫個文件 manifest file, 按以下格式寫，
開頭的行是注釋行會被自動忽略安岂，例如以下命名為為se-33-manifest的文件轻猖，保存為txt等文件
第一列為sample-ID, 第二列是絕對路徑，如果是conda環(huán)境里的路徑需要調(diào)整一下域那。第三列是序列文件對應的是reverse還是forward咙边。
我們這次示例的是個是雙端合并后的文件，所以不需要區(qū)分正向還是反向次员，direction一列败许，每個都寫forward。然后導出為UTF-8（逗號分隔）文件淑蔚。格式寫成后長這樣檐束。
sample-id,absolute-filepath,direction
SRR5142316,/xtdisk/PRJNA359624/SRR5142316.fastq,forward
SRR5142325,/xtdisk/PRJNA359624/SRR5142325.fastq.gz,forward
SRR5142334,/xtdisk/PRJNA359624/SRR5142334.fastq.gz,forward
如是雙端測序，分了兩個fastq文件束倍，格式寫成后大概長這樣，需要定指定 forward和reverse
sample-id,absolute-filepath,direction
SRR5142316,/xtdisk/PRJNA359624/SRR5142316_R1.fastq.gz,forward
SRR5142316,/xtdisk/PRJNA359624/SRR5142316_R2.fastq.gz,reverse
SRR5142325,/xtdisk/PRJNA359624/SRR5142325_R1.fastq.gz,forward
SRR5142334,/xtdisk/PRJNA359624/SRR5142334_R2.fastq.gz,reverse
- 一定要有表頭 sample-id,absolute-filepath,direction
3.3 運行 qiime tools盟戏，導入數(shù)據(jù)绪妹，生成數(shù)據(jù)摘要
qiime tools import \
--input-path se-33-manifest.txt \
--output-path qiime_out_manifest.qza \
--type SampleData[JoinedSequencesWithQuality] \
--input-format SingleEndFastqManifestPhred33
*雙端序列 --input-format PairedEndFastqManifestPhred33
--type 'SampleData[PairedEndSequencesWithQuality]'
*雙端序列合并 --input-format SingleEndFastqManifestPhred33
--type SampleData[JoinedSequencesWithQuality]
結果文件出來了，很快~
3.4 可視化數(shù)據(jù)質量
#可視化文件
qiime demux summarize \
--i-data qiime_out_manifest.qza \
--o-visualization qiime_out_manifest.qzv
qiime_out_manifest.qzv文件需要從Linux中下載后再拖拽到qiime2 view網(wǎng)頁中才能打開柿究。
此處可以得到質檢矢量圖邮旷，通過放大觀察可以清楚的判斷堿基質量明顯下降的位置，從而輔助確定下一步中流程中的參數(shù)reads1_cutpoint和reads2_cutpoint蝇摸。
- 生成的.qzv文 拖拽進網(wǎng)頁查看(推薦) QIIME 2 查看
- 或是使用qiime tools view qiime_out_manifest.qzv查看
這里發(fā)現(xiàn)數(shù)據(jù)質量都是30婶肩，其實這個數(shù)據(jù)有點奇怪办陷。
此處可以得到質檢矢量圖，通過放大觀察可以清楚的判斷堿基質量明顯下降的位置律歼，從而輔助確定下一步中的參數(shù)t民镜。
4 過濾， 切割险毁、去嵌合體制圈、拼接等
這步可使用Deblur 或 DATA 2。
Deblur VS DADA2： 方法不同畔况，使用相同數(shù)據(jù)集的序列和特征數(shù)量的也會有差異鲸鹦。
Deblur 要求所有序列的長度相同。另一方面跷跪，DADA2 的特征序列具有不同的長度馋嗜。DADA2 糾正錯誤，Deblur 刪除錯誤吵瞻。這種差異對到達終點的總讀取量有很大影響葛菇，但對頻率的影響要小得多，而頻率更重要听皿。
如果使用 DADA2來合并和消除雙端數(shù)據(jù)的噪聲熟呛，請在用DADA2去噪之前不要合并您的序列；DADA2希望讀長尚未合并的序列尉姨，并將在去噪過程中為您雙端合并庵朝。deblur只接受single-end的序列，如果你把沒有jion的序列傳遞給deblur又厉，它會只處理forward seqs,eblur在denoising時需要輸入整齊一樣長度的序列九府，所以需要trim成相同的長度。
目前：無論選擇使用哪種方法覆致，生物學解釋（至少從表面上看）似乎都非常相似
目前關于哪個更好侄旬，沒有達成共識，根據(jù)自己的數(shù)據(jù)情況來選擇煌妈。此教程我們的數(shù)據(jù)是合并后的數(shù)據(jù)儡羔，
那么選擇Deblur，注意一些數(shù)據(jù)分析比較原則璧诵，如果是整合數(shù)據(jù)集汰蜘，那么保證每個樣本使用相同的分析流程。
4.1 質量過濾
qiime quality-filter q-score \
--i-demux qiime_out_manifest.qza \
--o-filtered-sequences demux-joined-filtered.qza \
--o-filter-stats demux-joined-filter-stats.qza
輸出對象：
-
demux-joined-filter-stats.qza: 統(tǒng)計結果之宿。 -
demux-joined-filtered.qza: 數(shù)據(jù)過濾后結果族操。
此方法的參數(shù)尚未在雙端合并的數(shù)據(jù)上進行廣泛的基準測試，因此建議嘗試使用不同的參數(shù)設置比被。
4.2 去噪
現(xiàn)在已經(jīng)準備好用Deblur去噪你的序列了色难。您應該從質量分數(shù)圖中為--p-trim-length選擇合適的序列長度值泼舱。這將把所有序列修剪到這個長度，并丟棄任何小于這個長度的序列枷莉。
大約 251-255（EMP 引物組）捕獲了幾乎所有的 V4 序列娇昙。因此，修剪到 251-255 是合理的依沮。
Deblur的開發(fā)者們建議設置一個質量分數(shù)開始迅速下降的長度（recommend setting this value to a length where the median quality score begins to drop too low）涯贞。
qiime deblur denoise-16S \
--i-demultiplexed-seqs demux-joined-filtered.qza \
--p-trim-length 250 \
--p-sample-stats \
--o-representative-sequences rep-seqs.qza \
--o-table table.qza \
--o-stats deblur-stats.qza
輸出對象：
-
rep-seqs.qza: 代表序列。 -
deblur-stats.qza: 統(tǒng)計過程危喉。 -
table.qza: 特征表宋渔。
4.3 查看結果特征表
統(tǒng)計結果可視化
qiime feature-table summarize \
--i-table table.qza \
--o-visualization table.qzv
代表序列統(tǒng)計
qiime feature-table tabulate-seqs \
--i-data rep-seqs.qza \
--o-visualization rep-seqs.qzv
5 多樣性分析
5.1 建樹用于多樣性分析
qiime phylogeny align-to-tree-mafft-fasttree \
--i-sequences rep-seqs.qza \
--o-alignment aligned-rep-seqs.qza \
--o-masked-alignment masked-aligned-rep-seqs.qza \
--o-tree unrooted-tree.qza \
--o-rooted-tree rooted-tree.qza
5.2 Alpha多樣性分析
α多樣性是度量單個樣本內(nèi)有多少種微生物物種，以及每個物種所占比例的指標辜限。
計算多樣性(包括所有常用的Alpha和Beta多樣性方法)皇拣，輸入有根樹、Feature表薄嫡、樣本重采樣深度
取樣深度看table.qzv文件確定（一般為樣本最小的sequence count氧急，或覆蓋絕大多數(shù)樣品的sequence count）
還需要構建一個 sample-metadata.txt file ，這里我們的演示毫深，從NCBI下載
Run Selector :: NCBI (nih.gov)
改寫"RUN" 為 sampleID吩坝，這里我們添加了一個Group 列，根據(jù)Library Name 分了兩組Gout_Train 和HC_Train哑蔫，注意格式钉寝，保存為分隔符為空格的txt 文件。
然后運行以下流程：
qiime diversity core-metrics-phylogenetic \
--i-phylogeny rooted-tree.qza \
--i-table table.qza \
--p-sampling-depth 1775 \
--m-metadata-file sample-metadata.txt \
--output-dir core-metrics-results
運行很快闸迷，
輸出結果包括多種多樣性結果嵌纲，文件列表和解釋如下：
beta多樣性bray_curtis距離矩陣 bray_curtis_distance_matrix.qza
alpha多樣性evenness(均勻度，考慮物種和豐度)指數(shù) evenness_vector.qza
alpha多樣性faith_pd(考慮物種間進化關系)指數(shù) faith_pd_vector.qza
beta多樣性jaccard距離矩陣 jaccard_distance_matrix.qza
alpha多樣性observed_otus(OTU數(shù)量)指數(shù) observed_otus_vector.qza
alpha多樣性香農(nóng)熵(考慮物種和豐度)指數(shù) shannon_vector.qza
beta多樣性unweighted_unifrac距離矩陣腥沽，不考慮豐度 unweighted_unifrac_distance_matrix.qza
beta多樣性unweighted_unifrac距離矩陣逮走，考慮豐度 weighted_unifrac_distance_matrix.qza
測試分類元數(shù)據(jù)(sample-metadata)列和alpha多樣性數(shù)據(jù)之間的關聯(lián)，輸入多樣性值今阳、sample-medata师溅，輸出統(tǒng)計結果，統(tǒng)計faith_pd算法Alpha多樣性組間差異是否顯著
qiime diversity alpha-group-significance \
--i-alpha-diversity core-metrics-results/faith_pd_vector.qza \
--m-metadata-file sample-metadata.txt \
--o-visualization core-metrics-results/faith-pd-group-significance.qzv
# 統(tǒng)計evenness組間差異是否顯著
qiime diversity alpha-group-significance \
--i-alpha-diversity core-metrics-results/evenness_vector.qza \
--m-metadata-file sample-metadata.txt \
--o-visualization core-metrics-results/evenness-group-significance.qzv
以evenness-group-significance.qzv為例盾舌，圖中可點Category選擇分類方法险胰，查看不同分組下箱線圖間的分布與差別。圖形下面的表格矿筝，詳細詳述組間比較的顯著性和假陽性率統(tǒng)計。
發(fā)現(xiàn) Alpha 多樣性和 我們的分組沒啥關系棚贾，P值不顯著
5.3 Beta 多樣性分析
β多樣性是度量不同樣本間菌群組成的相似度大小的指標窖维，即關注各樣本間的菌群組成差異榆综。
# 統(tǒng)計group 的組間是否有顯著差異,其他的分組分析類似。
qiime diversity beta-group-significance \
--i-distance-matrix core-metrics-results/unweighted_unifrac_distance_matrix.qza \
--m-metadata-file sample-metadata.txt \
--m-metadata-column Group \
--o-visualization core-metrics-results/unweighted-unifrac-subject-significance.qzv \
--p-pairwise
可視化結果：
發(fā)現(xiàn)有P值铸史。
# 可視化三維展示其它類鼻疮，定量值的 分析
qiime emperor plot \
--i-pcoa core-metrics-results/unweighted_unifrac_pcoa_results.qza \
--m-metadata-file sample-metadata.txt \
--p-custom-axes weight \
--o-visualization core-metrics-results/unweighted-unifrac-emperor-weight.qzv
我們這里演示的數(shù)據(jù)集沒有，跳過傲宜。
5.4 Alpha rarefaction plotting
# --p-max-depth should be determined by reviewing the “Frequency per sample” information presented in the table.qzv file
# --p-max-depth一般取table.qzv文件Frequency per sample的中位數(shù)左右
qiime diversity alpha-rarefaction \
--i-table table.qza \
--i-phylogeny rooted-tree.qza \
--p-max-depth 31131 \
--m-metadata-file sample-metadata.txt \
--o-visualization alpha-rarefaction.qzv
可視化將有兩個圖簇爆。頂部圖是α稀疏圖唾糯，主要用于確定樣品的豐富度是否已被完全觀察或測序。如果圖中的線在沿x軸的某個采樣深度處看起來“平坦化”（即接近零斜率）挪哄，則表明收集超出該采樣深度的其他序列將不可能會有其他的OTU（feature）產(chǎn)生。如果圖中的線條沒有達到平衡琉闪，這可能是因為尚未完全觀察到樣品的豐富程度（因為收集的序列太少）迹炼，或者它可能表明在數(shù)據(jù)中存在大量的測序錯誤（被誤認為是新的多樣性）。底部圖表示當特征表稀疏到每個采樣深度時每個組中保留的樣本數(shù)颠毙。樣本被分成兩組,圖中顯示即兩條線.
6 菌群比對數(shù)據(jù)庫：引物特異性菌群比對
Greengenes細菌16S全長序列數(shù)據(jù)庫
下載地址:Data resources — QIIME 2 2023.9.2 documentation
下載得到gg_13_8_otus.tar.gz后將其解壓,具體使用gg_13_8_otus/rep_set 下99_otus.fasta（序列）和taxonomy/99_otu_taxonomy.txt（分類）兩個文件
tar -zxvf gg_13_8_otus.tar.gz
6.1. 轉換格式
#將99_otus.fasta（序列）和99_otu_taxonomy.txt（分類）兩個文件的格式轉換QIIME2能識別和利用的格式
qiime tools import \
--type 'FeatureData[Sequence]' \
--input-path gg_13_8_otus/rep_set/99_otus.fasta \
--output-path 99_otus.qza
qiime tools import \
--type 'FeatureData[Taxonomy]' \
--input-format HeaderlessTSVTaxonomyFormat \
--input-path gg_13_8_otus/taxonomy/99_otu_taxonomy.txt \
--output-path 99-taxonomy.qza
6.2. 抽提V3V4模板序列
用測序引物序列從Greengenes數(shù)據(jù)庫中的16S全長序列99_otus.qza中抽提出引物特異性的細菌參考序列斯入，就會得到本研究特異性的參考序列
qiime feature-classifier extract-reads \
--i-sequences 99_otus.qza \
--p-f-primer CCTACGGGNGGCWGCAG \
--p-r-primer GACTACHVGGGTATCTAATCC \
--o-reads 99-v3v4-seqs.qza
大概20分鐘
6.3. 訓練V3V4分類器
qiime feature-classifier fit-classifier-naive-bayes \
--i-reference-reads 99-v3v4-seqs.qza \
--i-reference-taxonomy 99-taxonomy.qza \
--o-classifier gg-13-8-99-v3v4-classifier.qza
大概20分鐘 等待中~
7 菌群特征比對，分類學注釋
7.1. 比對
把DADA2分析得到的細菌特征代表性序列rep-seqs.qza和訓練好的分類器gg-13-8-99-v3v4-classifier.qza進行比對蛀蜜，獲得具體的細菌分類信息taxonomy.qza
qiime feature-classifier classify-sklearn \
--i-classifier gg-13-8-99-v3v4-classifier.qza \
--i-reads rep-seqs.qza \
--o-classification taxonomy.qza
7.2. 豐度統(tǒng)計
將細菌特征豐度表table.qza和細菌分類信息taxonomy.qza進行整合獲得完整的細菌分類豐度表刻两，包含界、門滴某、綱磅摹、目、科壮池、屬偏瓤、種多水平的細菌豐度信息
qiime taxa barplot \
--i-table table.qza \
--i-taxonomy taxonomy.qza \
--m-metadata-file sample-metadata.txt \
--o-visualization taxa-bar-plots.qzv
8 數(shù)據(jù)處理: 獲取菌屬分類表
8.1. 獲取屬水平菌豐度表
QIIME2 view網(wǎng)頁中打開taxa-bar-plots.qzv，下載level6的CSV文件椰憋，如下：
8.2. 標準化菌屬豐度表
把CSV文件導入到Excel中進行標準化厅克，即每個菌屬的原始豐度除以該菌所在樣本的總菌屬豐度得到標準相對菌屬相對豐度
