旅游网站需求分析,wordpress 广告 能赚多少,自己建设网站赚钱,access数据库做网站目录 创建metadata 把数据导入qiime2 去除引物序列 双端合并 #xff08;dada2不需要#xff09; 质控 #xff08;dada2不需要#xff09; 使用deblur获得特征序列 使用dada2生成代表序列与特征表 物种鉴定 可视化物种鉴定结果 构建进化树#xff08;ITS一般不构建进化树… 目录 创建metadata 把数据导入qiime2 去除引物序列 双端合并 dada2不需要 质控 dada2不需要 使用deblur获得特征序列 使用dada2生成代表序列与特征表 物种鉴定 可视化物种鉴定结果 构建进化树ITS一般不构建进化树 生成多样性分析所需数据 α多样性分析 beta多样性分析 绘制稀释曲线 差异物种分析 如何把.qza格式文件导出比如把特征表导出 演示单端数据导入qiime2 基于qiime2的16s测序数据分析分析 本文演示的是目前最长用的双端测序数据 激活qiime环境 alias activate_qiimeconda activate qiime2-amplicon-2024.2 我已经把上面这句代码放在了.bashrc里面因此激活qiime2环境只需要运行activate_qiime即可 按照下图策略分析 演示Illumina公司的16S rRNA基因区的V3到V4区 数据导入需要首先准备好一个manifest.tsv文件这个文件有三列分别是 sample-id forward-absolute-filepath reverse-absolute-filepath如图所示 生成manifest.tsv可以使用下面的代码这个代码是my_script文件夹里面的create_manifest.sh #!/bin/bash # 设置目录路径 dir/data3/sunjs/chenlina/69contorl47PDAC # 输出表头 echo -e sample-id\tforward-absolute-filepath\treverse-absolute-filepath # 遍历_R1和_R2文件 for r1 in ${dir}/*_R1.fastq.gz; do # 检查R1文件是否存在 if [[ -f $r1 ]]; then # 获取样本ID去掉_R1.fastq.gz sample_id$(basename $r1 _R1.fastq.gz) # 构建R2文件路径 r2${dir}/${sample_id}_R2.fastq.gz # 检查R2文件是否存在 if [[ -f $r2 ]]; then # 输出结果 echo -e ${sample_id}\t${r1}\t${r2} fi fi done 这个脚本会在终端直接显示生成的内容我们一般会把这些内容重定向到另一个文件去所以通常运行代码 bash create_manifest.shmanifest.tsv 这样就会在当前工作目录下生成一个manifest.tsv文件 创建metadata 还需要创建一个metadata比如有这些列 不过我创建的meta.tsv如图所示只有样本名和分组两列 把数据导入qiime2 export TMPDIR/data3/sunjs/ cd ~/test/qiime2/ input_path~/test/qiime2/manifest.tsv qiime tools import \ --type SampleData[PairedEndSequencesWithQuality] \ --input-path ${input_path} \ --output-path PDAC.qza \ --input-format PairedEndFastqManifestPhred33V2 去除引物序列 v3到v4区域的引物序列是固定的就是这两个参数指定的序列 cd /data3/sunjs/chenlina/69contorl_47PDAC_res qiime cutadapt trim-paired \ --i-demultiplexed-sequences 69contorl_47PDAC.qza \ --p-cores 10 \ --p-no-indels \ --p-front-f ACTCCTACGGGAGGCAGCAG \ --p-front-r GGACTACHVGGGTWTCTAAT \ --o-trimmed-sequences primer-trimmed-demux.qza 双端的数据不需要拆分barcode直接往下进行即可 双端合并 dada2不需要 如果接下来打算 使用deblur或OTU聚类方法现在就合并序列。如果计划使用 dada2对序列进行去噪则不要合并——dada2会在对每个序列进行去噪之后自动执行序列合并。 qiime vsearch merge-pairs \ --i-demultiplexed-seqs primer-trimmed-demux.qza \ --o-unmerged-sequences unmerged_demux-joined.qza \ --p-threads 16 \ --o-merged-sequences demux-joined.qza 质控 dada2不需要 cd /data3/sunjs/chenlina/test qiime quality-filter q-score \ --p-min-quality 20 \ --i-demux demux-joined.qza \ --o-filtered-sequences demux-joined-filtered.qza \ --o-filter-stats demux-joined-filter-stats.qza 使用deblur获得特征序列 qiime deblur denoise-16S \ --i-demultiplexed-seqs demux-joined-filtered.qza \ --p-trim-length 400 \ --p-sample-stats \ --p-jobs-to-start 2 \ --p-min-reads 1 \ --o-representative-sequences repset-seqs.qza \ --o-table feature-table.qza \ --o-stats deblur-stats.qza 或者使用data2得到特征表和代表序列可以直接替换掉质控使用 deblur获得特征序列两步但是不管用deblur还是dada2都要完成导入数据去引物 使用dada2之前需要先获得 --p-trim四个参数把下面这段代码生成的数据拖到qiime2的可视化网站里面可以看到两张图 qiime demux summarize \ --i-data primer-trimmed-demux.qza \ --o-visualization demux-summary.qzv 使用dada2生成代表序列与特征表 #p-trim-left-f正向序列左边 #p-trim-left-r 反向去列左边 #p-trunc-len-f正向序列x轴的值 #p-trunc-len-r反向序列x轴的值 # --p-retain-all-samples True保留所有样本 qiime dada2 denoise-paired \ --i-demultiplexed-seqs primer-trimmed-demux.qza \ --p-n-threads 20 \ --p-trim-left-f 0 --p-trim-left-r 0 \ --p-trunc-len-f 270 --p-trunc-len-r 187 \ --p-retain-all-samples True \ -o-table dada2-table.qza \ --o-representative-sequences dada2-rep-seqs.qza \ --o-denoising-stats denoising-stats.qza 其中 #p-trim-left-f正向序列左边 #p-trim-left-r 反向去列左边 #p-trunc-len-f正向序列x轴的值 #p-trunc-len-r反向序列x轴的值 这四个参数的取值要通过看上面那两张图来得到 一般来讲会去掉质量中位数小于20的位点。 上面两种方法选一个就行推荐dada2 物种鉴定 生成代表序列之后可以进行物种鉴定了 export TMPDIR/data3/sunjs/ input_path/data3/sunjs/chenlina/qiime2/PDAC output_path/data3/sunjs/chenlina/qiime2/PDAC db_path~/Metagenome/database/silva-138-99-nb-classifier.qza cd ${input_path} #全长引物注释,使用的数据库日期是2024.2 #要求输入代表序列也就是dada2结果生成的那个 #--p-n-jobs如果设置为 -1那么程序会尝试使用所有可用的 CPU 核心来执行任务 #--p-n-jobs指定的是进程而非线程和threads不是一回事 qiime feature-classifier classify-sklearn \ --i-classifier ${db_path} \ --i-reads dada2-rep-seqs.qza \ --p-n-jobs 20 \ --o-classification ${output_path}/taxonomy.qza #过滤污染 叶绿体 线粒体 qiime taxa filter-table \ --i-table dada2-table.qza \ --i-taxonomy taxonomy.qza \ --p-exclude mitochondria,chloroplast \ --p-include p__ \ --o-filtered-table ${output_path}/feature-table-filt-contam.qza #过滤稀有ASV qiime feature-table filter-features \ --i-table feature-table-filt-contam.qza \ --p-min-frequency 50 \ --p-min-samples 1 \ --o-filtered-table ${output_path}/feature-table-final.qza #过滤掉ASV对应的序列 qiime feature-table filter-seqs \ --i-data dada2-rep-seqs.qza \ --i-table feature-table-final.qza \ --o-filtered-data ${output_path}/repset-seqs-final.qza 这里需要提前下载好对应的数据库下载的数据库理论上来讲是越新越好但是需要和你用的qiime2软件版本一致比如qiime2是2024.2的那数据库也要下载2024.2的下载数据库可以从这个网址下载 QIIME 2 Resources 运行这段代码之后就得到了物种鉴定表也就是 taxonomy.qza同时上面的代码还进行了一些过滤然后我们可以使用下面的代码来对物种鉴定的结果进行可视化 可视化物种鉴定结果 qiime metadata tabulate \ --m-input-file taxonomy.qza \ --o-visualization taxonomy.qzv #绘制柱状图 qiime taxa barplot \ --i-table feature-table-final.qza \ --i-taxonomy taxonomy.qza \ --m-metadata-file meta.tsv \ --o-visualization taxa-bar-plots.qzv 其中taxonomy.qzv文件拖到网站里面长这样子 taxa-bar-plots.qzv长这样子网站允许把柱状图的数据导出成csv格式推荐导出之后用R来绘制 构建进化树ITS一般不构建进化树 #构建进化树ITS一般不构建进化树 export TMPDIR/data3/sunjs/ input_path/data3/sunjs/chenlina/qiime2/PDAC output_path/data3/sunjs/chenlina/qiime2/PDAC cd ${input_path} qiime phylogeny align-to-tree-mafft-fasttree \ --i-sequences repset-seqs-final.qza \ --o-alignment aligned-repset-seqs.qza \ --p-n-threads 40 \ --o-masked-alignment masked-aligned-repset-seqs.qza \ --o-tree ${output_path}/unrooted-tree.qza \ --o-rooted-tree ${output_path}/rooted-tree.qza 下面这段代码用于生成代表序列和特征表的qzv格式文件在网址 view.qiime2.org可以可视化需要先导出然后拖进去 cd /data3/sunjs/chenlina/84contorl_50PDAC_res # 特征表和特征序列汇总 qiime feature-table summarize \ --i-table feature-table-final.qza \ --m-sample-metadata-file meta.tsv \ --o-visualization table.qzv qiime feature-table tabulate-seqs \ --i-data repset-seqs-final.qza \ --o-visualization rep-seqs.qzv 把table.qzv拖到网站里面是这样的 点击detail这个拖到最下边出现这样的界面 后续将选择53这个数量进行抽平当然51也可以 生成多样性分析所需数据 #53是把上一步生成的table.qzv文件拖到qimme2网站查看来确定的 qiime diversity core-metrics-phylogenetic \ --i-phylogeny rooted-tree.qza \ --i-table feature-table-final.qza \ --p-sampling-depth 53 \ --m-metadata-file meta.tsv \ --output-dir ./diversity 生成一个名为diversity的目录里面有一堆文件其中有四个是qzv格式的这四个qzv格式的可以直接拖到网站进行可视化都是beta多样性的比如PCA图PCoA图这种的 比如下面这样子 α多样性分析 cd /data3/sunjs/chenlina/84contorl_50PDAC_res/ #α多样性分析 qiime diversity alpha-group-significance \ --i-alpha-diversity diversity/faith_pd_vector.qza \ --m-metadata-file meta.tsv \ --o-visualization faith-group-significance.qzv 输入是上一步生成的faith_pd_vector.qza生成文件如图所示 beta多样性分析 #beta多样性分析 cd /data3/sunjs/chenlina/67contorl_50PDAC_res qiime diversity beta-group-significance \ --i-distance-matrix diversity/unweighted_unifrac_distance_matrix.qza \ --m-metadata-file meta.tsv \ --m-metadata-column group \ --p-pairwise \ --o-visualization unweighted-unifrac-group-significance.qzv 生成了一个文件拖到网页里面是这样的代表了其中某个分组与其他分组有没有差别使用的方法是置换多因素方差分析 画图实际上前面已经有过这个图了 #qiime diversity core-metrics-phylogenetic这一步已经生成的结果 cd /data3/sunjs/chenlina/69contorl_47PDAC_res qiime emperor plot \ --i-pcoa diversity/unweighted_unifrac_pcoa_results.qza \ --m-metadata-file meta.tsv \ --p-custom-axes days-since-experiment-start \ --o-visualization unweighted-unifrac-emperor-group.qzv 这里--i输入的是前面 qiime diversity core-metrics-phylogenetic 这一步已经生成的结果最后得到的结果也是跟前面 qiime diversity core-metrics-phylogenetic运行得到的那个结果是一样的尚不清楚这个有啥用处。 绘制稀释曲线 #稀释曲线 cd /data3/sunjs/chenlina/69contorl_47PDAC_res qiime diversity alpha-rarefaction \ --i-table feature-table-final.qza \ --i-phylogeny rooted-tree.qza \ --p-max-depth 9000 \ --m-metadata-file meta.tsv \ --o-visualization alpha-rarefaction.qzv 我这里选择了取样是9000条结果如图所示 可以看到后面基本平了说明测序饱和了如果发现后面还在往上升说明测序没饱和很多东西测不到这种数据发出去会被拒稿因为说不定再多测点就结果不一样了。 差异物种分析 有时候可能会先过滤部分数据比如我先按照频率和特征数过滤 掉一些低质量数据 #过滤部分特征表的数据 cd /data3/sunjs/chenlina/69contorl_47PDAC_res #--p-min-frequency样本必须具有的最小总频率才能被保留 #--p-min-features样本必须具有的最小特征数才能被保留 qiime feature-table filter-samples \ --i-table feature-table-final.qza \ --m-metadata-file meta.tsv \ --p-min-features 10 \ --p-min-frequency 1500 \ --o-filtered-table filter-table.qza 可以查看qiime feature-table filter-samples的帮助文档来查看其他过滤参数比如我可以使用参数--p-where只保留meta.tsv中的某一个分组。 接下来使用 ANCOM-BC方法来进行差异物种分析 #差异分析 cd /data3/sunjs/chenlina/69contorl_47PDAC_res #--p-formula描述了微生物的绝对丰度是如何依赖于元数据中的变量的 #使用ANCOM-BC方法来进行差异物种分析 qiime composition ancombc \ --i-table filter-table.qza \ --m-metadata-file meta.tsv \ --p-formula group \ --o-differentials ancombc-group.qza #生成差异丰富度分析结果的条形图可视化 qiime composition da-barplot \ --i-data ancombc-group.qza \ --p-significance-threshold 0.001 \ --o-visualization da-barplot-group.qzv 如何把.qza格式文件导出比如把特征表导出 cd /data3/sunjs/chenlina/84contorl_50PDAC_res qiime tools export \ --input-path feature-table-final.qza \ --output-path exported-feature-table biom convert \ -i exported-feature-table/feature-table.biom \ -o feature-table.tsv \ --to-tsv 这就把特征表导出成了tsv这种表格形式的文件 演示单端数据导入qiime2 先使用fastp对下载的数据进行质控 input_path/data3/sunjs/chenlina/contorl_data output_path/data3/sunjs/chenlina/contorl_clean_data cd $input_path #SRR15884889.fastq for file in *.fastq; do sample_id$(echo $file | cut -d. -f1) echo sample_id是$sample_id 输入文件是$file echo 输出文件为${sample_id}.fastp.json和${sample_id}.fastp.html fastp -i ${file} -o ${output_path}/${file} \ - -qualified_quality_phred 20 \ --length_required 50 \ --cut_front \ --cut_tail \ --trim_poly_g \ --html ${output_path}/${sample_id}.html \ --json ${output_path}/${sample_id}.json done 然后运行代码自动构建一个manifest.tsv文件 cd /data3/sunjs/chenlina/contorl_data echo -e sample-id\tabsolute-filepath contorl_manifest.tsv for file in *.fastq; do echo -e ${file%.fastq}\t$(pwd)/$file contorl_manifest.tsv done 解释下上面这段代码 上面的命令是用来生成一个清单文件manifest file这个文件将每个样本的ID和对应的绝对文件路径关联起来。这个文件随后可以被用来将单端FASTQ文件导入到QIIME 2中。 数据来源于文章 Dysbiosis of human gut microbiome in young onset colorectal cancer他的数据中扩增子测序文件里面是单端测序文件所以构建的manifest长这样子 然后把健康人的feces数据导入qiime2由于是单端的所以 --type 参数是SampleData[SequencesWithQuality] export TMPDIR/data3/sunjs/ input_path/data3/sunjs/chenlina/contorl_data/contorl_manifest.tsv qiime tools import \ --type SampleData[SequencesWithQuality] \ --input-path ${input_path} \ --output-path contorl_demux-single-end.qza \ --input-format SingleEndFastqManifestPhred33V2
