科学网

 找回密码
  注册

tag 标签: QIIME

相关帖子

版块 作者 回复/查看 最后发表

没有相关内容

相关日志

扩增子分析流程
wbb121 2019-4-21 20:08
总目录 http://mp.weixin.qq.com/mp/homepage?__biz=MzUzMjA4Njc1MA==hid=4sn=85f768f81736c77cadb1952cb3dd8ff7scene=18#wechat_redirect QIIME的虚拟机安装 https://mp.weixin.qq.com/s?__biz=MzUzMjA4Njc1MA==mid=2247483675idx=1sn=ffa42b6f5dd32b867bae247021d8e975scene=19#wechat_redirect 扩增子分析解读1质控,实验设计,双端序列合并 https://mp.weixin.qq.com/s?__biz=MzUzMjA4Njc1MA==mid=2247483675idx=1sn=ffa42b6f5dd32b867bae247021d8e975scene=19#wechat_redirect 扩增子分析解读2提取barcode,质控及样品拆分,切除扩增引物 https://mp.weixin.qq.com/s?__biz=MzUzMjA4Njc1MA==mid=2247483858idx=1sn=3c5bfa6108e1a71cd408db6f02a96f48scene=19#wechat_redirect 扩增子分析解读3格式转换,去冗余,聚类 https://mp.weixin.qq.com/s?__biz=MzUzMjA4Njc1MA==mid=2247483860idx=1sn=f84a97b9fc739eb69a3706511496aafcscene=19#wechat_redirect 扩增子分析解读4去嵌合体,非细菌序列,生成代表性序列和OTU表 https://mp.weixin.qq.com/s?__biz=MzUzMjA4Njc1MA==mid=2247483865idx=1sn=abce77d58ea773bc1964d26d61672003scene=19#wechat_redirect 扩增子分析解读5物种注释,OTU表操作 https://mp.weixin.qq.com/s?__biz=MzUzMjA4Njc1MA==mid=2247483871idx=1sn=4f4b4b5e1f793b73d7c3f85c8067ad26scene=19#wechat_redirect 扩增子分析解读6进化树,Alpha,Beta多样性 https://mp.weixin.qq.com/s?__biz=MzUzMjA4Njc1MA==mid=2247483873idx=1sn=31be33db6e39fb616b9cc27f199b169escene=19#wechat_redirect 扩增子分析解读7物种分类统计,筛选进化树和其它 https://mp.weixin.qq.com/s?__biz=MzUzMjA4Njc1MA==mid=2247483883idx=1sn=88894caaae1ddcdc3b811405e414d9a9scene=19#wechat_redirect
个人分类: 生物统计软件|3712 次阅读|0 个评论
Striped UniFrac
wbb121 2018-11-8 09:38
参考: https://www.nature.com/articles/s41592-018-0187-8#MOESM6 https://github.com/biocore/unifrac 这是一种更快速,占用内存更少的计算UniFrac距离矩阵的方法,支持unweighted UniFrac,weighted UniFrac, generalized UniFrac, variance adjusted weighted UniFrac, 和meta UniFrac. 有多种方式来安装库,最常使用的是两种,一种是QIIME2,一种是pip。最简单的方式是QIIME2,默认安装此模块。 QIIME2的安装见 https://docs.qiime2.org/2018.2/install/ 。 (未完)
个人分类: 生物统计软件|2664 次阅读|0 个评论
使用QIIME分析微生物群落的16S rRNA序列(fastq格式序列)
wbb121 2018-5-31 10:01
参考 http://nbviewer.jupyter.org/github/biocore/qiime/blob/1.9.1/examples/ipynb/illumina_overview_tutorial.ipynb https://forum.qiime2.org/t/qiime2-chinese-manual/838 下载数据 创建文件夹emp-single-end-sequences: mkdiremp-single-end-sequences fastq格式的序列文件: wget-Oemp-single-end-sequences/sequences.fastq.gzhttps://data.qiime2.org/2018.4/tutorials/moving-pictures/emp-single-end-sequences/sequences.fastq.gz barcode sequences: wget-Oemp-single-end-sequences/barcodes.fastq.gzhttps://data.qiime2.org/2018.4/tutorials/moving-pictures/emp-single-end-sequences/barcodes.fastq.gz 生成映射文件并检查映射文件的正确性 映射文件的要求见 http://qiime.org/documentation/file_formats.html#qiime-parameters 此处给出两个例子,一个正确的, map.tsv ,一个错误的, map-bad.tsv 使用 validate_mapping_file.py 检查映射文件map.tsv的正确性,输出 日志文件,html文件和corrected_mapping.txt文件。 validate_mapping_file.py-o./vmf-map-m./map.tsv 此时给出信息:No errors or warnings were found in mapping file. 检测的结果在文件夹./vmf-map中。 使用 validate_mapping_file.py 检查映射文件map-bad.tsv的正确性 , validate_mapping_file.py-o./vmf-map-bad-m./map-bad.tsv 此时给出信息:Errors and/or warnings detected in mapping file. Please check the log and html file for details. 可以查看产生的HTML摘要以找出存在的错误。然后在电子表格程序或文本编辑器中修复这些问题,重新运行validate_mapping_file.py检查更新后的映射文件。 质量过滤序列 使用 split_libraries_fastq.py 对序列进行质量过滤,一般来说,序列和条形码有单独的fastq文件, split_libraries_fastq.py-oslout/-iforward_reads.fastq.gz-bbarcodes.fastq.gz-mmap.tsv OTU-picking 注意若是在NCBI等网站上下载处理过的序列,直接在此处开始即可 有三种策略,为 pick_closed_reference_otus.py , pick_open_reference_otus.py , pick_de_novo_otus.py 。此处以 pick_open_reference_otus.py 为例 pick_open_reference_otus.py-ootus/-islout/seqs.fna-p../uc_fast_params.txt 注意,该命令采用上一步中生成的文件seqs.fna。我们还为该命令指定了一些参数,这是该工作流程的内部。 我们从这个命令中获得的主要输出是 OTU table ,或者每个样品中观察到的每个操作分类单位(OTU)的次数。QIIME使用基因组学标准联盟生物观察矩阵标准(BIOM)格式来表示OTU表格。可以在 这里 找到有关BIOM格式的更多信息,以及将这些文件转换为制表符分隔文本的信息,这些文本可以在 此处 电子表格程序中 查看 。这个命令生成几个OTU表。我们通常使用./otus/otu_table_mc2_w_tax_no_pynast_failures.biom。它有单个OTU(或总数为1的OTU),以及其代表序列不能与 PyNAST 对齐的 OTU 。它还包含每个OTU的分类分配作为观测元数据。pick_open_reference_otus.py命令还产生系统发育树,包含树的文件是./otus/rep_set.tre,并且是./otus/otu_table_mc2_w_tax_no_pynast_failures.biom下游系统发育多样性计算中应该使用的文件。树以广泛使用的 newick格式存储 。 (未完)
个人分类: 生物统计软件|6225 次阅读|0 个评论
使用QIIME分析微生物群落的16S rRNA序列(fasta格式序列)
wbb121 2018-5-26 19:02
参考文章Kuczynski J, Stombaugh J, Walters W A, et al. Using QIIME to analyze 16S rRNA gene sequences from Microbial Communities . Current protocols in bioinformatics / editoral board, Andreas D. Baxevanis. , 2011, Chapter 10:Unit 10.7. 以及 http://qiime.org/tutorials/tutorial.html https://forum.qiime2.org/t/qiime2-chinese-manual/838 QIIME是一种执行微生物群落分析的软件应用程序。它是Quantitative Insights In Microbial Ecology的缩写,已被用于分析和解释来自真菌,病毒,细菌和古菌群落的核酸序列数据。 Unit Introduction 标准QIIME分析从一种或多种测序技术(例如Sanger,Roche / 454,Illumina或其他)的序列数据开始。使用QIIME分析微生物群落的数据包括在终端窗口中输入一系列命令,然后查看图形和文本输出。使用Linux风格的命令行界面(即命令有一些相当基本的了解 cd , ls 以及使用制表符完成的)是有用的,但并不是必需的。 我们用一个从小鼠肠道微生物群落对空腹的反应研究中得到的数据作为例子,而不是用一般术语列出分析步骤。为了使这个例子在个人计算机上快速运行,我们使用由控制随意饮食饲养的5只动物产生的数据的一个子集,并且4只动物在处死前24小时禁食。在基本操作完成后,我们比较了对照组与禁食动物的群落结构,特别是,我们比较了禁食小鼠和非禁食小鼠的微生物群落的分类学特征,观察样品内多样性指标的差异,在各组之间进行比较,并进行比较聚类分析以查找样本中的总体差异。 下载数据 可在 ftp://ftp.microbiome/pub/qiime-files/qiime_overview_tutorial.zip 下载,或者在附件中下载。 文件的说明如下: .fna文件: 这是454机器生成的FASTA文件。 .qual文件: 这是454机器生成的质量得分文件,其中包含FASTA文件中每个序列的基准得分。 tsv.txt文件: 映射文件由用户生成。该文件包含执行数据分析所需的所有样本信息。具体要求见 http://qiime.org/documentation/file_formats.html#metadata-mapping-files。 在此处我们使用映射文件 Fasting_Map.txt 硬件要求 安装 QIIME, 可以通过虚拟机安装,这是最简单的方法,见 http://qiime.org/install/index.html 其他必要资源 在比对16S DNA序列时使用的文件,可以在终端使用如下命令下载,或直接从附件中下载(core_set_aligned.fasta.imputed文件太大,附件无法上传,请自行下载) wgethttp://greengenes.lbl.gov/Download/Sequence_Data/Fasta_data_files/core_set_aligned.fasta.imputed wgethttp://greengenes.lbl.gov/Download/Sequence_Data/lanemask_in_1s_and_0s 预处理 将目录设置为数据所在文件夹(此处需要个人做相应的修改) cd/media/sf_E/alignment-free/16s_rRNA/OTU/qiime2/qiime_overview_tutorial 检查映射文件,所用的命令为 validate_mapping_file.py , validate_mapping_file.py-mFasting_Map.txt-omapping_output 该模块会显示一条消息,指示映射文件中是否存在问题。错误和警告将输出到指定(-o)输出目录中的日志文件。错误会导致后续脚本的致命问题,并且需要在继续下一步分析之前纠正。警告不会导致致命的问题,但鼓励解决这些问题,因为它们通常指示映射文件中的拼写错误,无效字符或将影响下游分析的其他意外错误。正确的 corrected_mapping.txt文件也在输出目录中,这是替换无效字符的映射文件的副本。 基于核苷酸条码(nucleotide barcode)给样本分配multiplexed reads,这一步也完成了质量过滤,去除低质量或模糊的reads。所用的命令为 split_libraries.py ,这里我们采用默认参数, split_libraries.py-mFasting_Map.txt-fFasting_Example.fna-qFasting_Example.qual-osplit_library_output 这一命令会在新目录split_library_output/中创建三个文件 : split_library_log.txt:此文件包含分割摘要,包括每个样本检测到的reads的数量以及由于质量考虑而被删除的reads的简要摘要。 histograms.txt:这个制表符分隔的文件显示了拆分之前和之后以规则大小间隔reads的数目。 seqs.fna:这是一个FASTA格式的文件,每个序列根据它来自的样本进行重命名。标题行还包含FASTA输入文件中的read名称以及有关更正的条形码错误的信息。 Picking OTUs 注意若是在NCBI等网站上下载处理过的序列,直接在此处开始即可 在这里运行 pick_de_novo_otus.py ,该工作流程会自动运行一系列其他脚本。该工作流程由以下步骤组成: 1.根据读取内的序列相似性选择OTUs( pick_otus.py ) 2.为每个OTU选择一个代表序列( pick_rep_set.py ) 3.将分类标准分配给OTU代表序列( assign_taxonomy.py ) 4.比对OTU代表序列( align_seqs.py ) 5. 比对过滤( filter_alignment.py ) 6. 构建系统发育树( make_phylogeny.py ) 7.制作OTU表( make_otu_table.py ) 使用split_libraries.py的输出(seqs.fna),运行以下命令: pick_de_novo_otus.py-i./split_library_output/seqs.fna-ootus 1.根据读取内的序列相似性选择OTUs(pick_otus.py) 在此步骤中,所有样本的所有序列将根据它们的序列相似性聚类为操作分类单元(OTU)。QIIME中的OTUs是序列的聚类,通常意在表示一定程度的分类学相关性。例如,当序列以97%的序列相似性聚类时,每个产生的聚类通常被认为代表一个物种。这种模式和目前采摘OTU的技术已知存在缺陷,然而,97%的OTUs与人类对许多微生物所称的物种不匹配。确定OTU应该如何定义以及它们代表什么是一个活跃的研究领域。见这里关于OTU与QIIME采摘的更多信息。 2. 为每个OTU选择一个代表序列(pick_rep_set.py) 由于每个OTU可能由许多相关序列组成,因此我们将从每个OTU中选择一个代表序列用于下游分析。该代表性序列将用于OTU的分类学鉴定和系统发育比对。QIIME使用上面创建的OTU文件,并通过几种方法之一从fasta文件中提取代表序列。 在otus/rep_set/目录中,QIIME有一个fasta文件seqs_rep_set.fasta,每个OTU包含一个代表序列。在此fasta文件中,序列已用OTU标识符重命名,并且标题行上的附加信息反映了用作代表的序列: 3.将分类标准分配给OTU代表序列(assign_taxonomy.py) 接下来,分类将被分配给每个代表性序列。默认情况下,QIIME使用uclust共识分类法分类器来尝试将分类法分配给步骤2产生的每个代表性序列。 在目录 otus/uclust_assigned_taxonomy/中 ,会有一个日志文件和一个文本文件。文本文件(我们称之为观察或OTU元数据文件)包含每个OTU的一行,然后是分类指定,包含此分类指定的uclust命中分数以及找到的uclust命中数。对于一些OTU,分配可能与细菌种类一样具体,而其他分配可能根本不可分配(因此将被标记为 未分配 )。以下是观测元数据文件的前几行,其中包含uclust分类学分配的结果: 4.比对OTU代表序列(align_seqs.py) 只有在随后调用诸如 UniFrac 系统发育指标时,OTU代表序列和系统发生推断的对齐才是必要的。比对可以使用诸如MUSCLE之类的程序从头生成,也可以通过使用像 PyNAST 这样的工具与现有对齐进行对齐。对于小型研究,任何一种方法都是可行的。然而,对于涉及多个序列(大约超过1000个)的研究,重新校准器非常缓慢并且需要与 PyNAST 比对。由于比对是流水线中计算最密集的瓶颈之一,所以大型研究从这项任务的并行化中受益匪浅,这在 PyNAST 中是可能的。 对齐序列后,将在./otus/pynast_aligned_seqs /目录中创建一个日志文件和一个比对文件。 5.比对过滤(filter_alignment.py) 在推断与序列有关的系统发育树之前,过滤序列比对以去除仅由缺口组成的列和已知过度变化的位置是有益的。QIIME默认使用16S对齐Lane mask(Lane,DJ 1991)。过滤后在目录./otus/pynast_aligned_seqs/中创建一个过滤的比对文件。 6.构建系统发育树(make_phylogeny.py) 使用目录./otus/pynast_aligned_seqs/中生成的已过滤的比对文件,使用树构建程序构建系统发育树。 该Newick树格式文件写入rep_set.tre,位于./otus目录中。该文件可以在树形图可视化软件中查看,并且对于UniFrac多样性测量和其他系统发生分析是必要的。所获得的树可以用FigTree等程序可视化,该程序用于显示存储在rep_set.tre中的系统发育树: 7.制作OTU表(make_otu_table.py) 该步骤的结果是otu_table.biom,位于./otus/目录中。有关以BIOM格式存储的OTU表格格式的更多信息,请参阅 http://biom-format.org/ 。 总结otu表,查看OTU表的摘要统计信息,使用如下命令: biomsummarize-table-i./otus/otu_table.biom 总结表明,本教程示例中的序列相对较少,但所存在的序列在9个微生物群落中分布相当均匀: Num samples: 9 Num observations: 419 Total count: 1337 Table density (fraction of non-zero values): 0.168 Counts/sample summary: Min: 146.0 Max: 150.0 Median: 149.000 Mean: 148.556 Std. dev.: 1.257 Sample Metadata Categories: None provided Observation Metadata Categories: taxonomy Counts/sample detail: PC.481: 146.0 PC.355: 147.0 PC.636: 148.0 PC.607: 149.0 PC.635: 149.0 PC.593: 149.0 PC.354: 149.0 PC.634: 150.0 PC.356: 150.0 (未完) 附件: Using QIIME to analyze 16s rRNA gene sequences from microbial communities.pdf qiime_overview_tutorial.zip
个人分类: 生物统计软件|17323 次阅读|0 个评论
Alpha多样性指数之Chao1指数
热度 2 luria 2017-9-8 19:44
1. 准备 在认识 Chao1 指数前有两个概念需要理清: singletons ,即仅含有一条 read 的 OTU , doubletons ,即仅含有两条 reads 的 OTU 。 可以这样理解:在一个放了各种各样玩具模型的水池中(例如下图。水池很大,其中玩具有相同的,有不同的,且各种类型及数目不限),随机来捞玩具。这时捞起来一个,发现之前有个玩具和这个捞起的玩具一模一样,这时有两个这种玩具在手上,这个玩具模型就是 doubletons ;当然也可能捞起一个玩具发现手里没有相同的,那这个就叫 singletons 。 注意:手里如果已经有两个或以上,再捞一个起来这类情况对 Chao1 指数是没有贡献的! 2. 公式 Chao1 的经典公式如下 〔1〕 : S obs 表示样本中观察到的物种数目。 F1 和 F2 分别表示 singletons 和 doubletons 的数目。 Chao1 指数还有另外一种修正偏差的公式,在 scikit-bio 上也有提到了 ( 注: QIIME 使用的是 scikit-bio) ,如下: 由经典公式可以看到,当 doubletons 为 0 (即 F 2 为 0 )时计算的结果没有意义,修正公式可以解决这个问题。 可以这样理解这个修正公式( 虽然不太严格 ):它从 singletons 中拿出 1 条来(严格来说与经典公式相比还不到 1 条),当作 doubletons ,这样分母一定会大于 0 。 从 QIIME 代码来看, QIIME 调用 skbio.diversity.alpha .chao1 时 bias_corrected 为默认值 True ,表示按经典公式计算。 3. 对公式的理解 从公式我们可以看到 chao1 指数是用来反映物种丰富度的指标。它通过观测到的结果推算出一个理论的丰富度,这个丰富度更接近真实的丰富度。一般来讲能观测到的物种丰富度肯定会比实际少,那么两者之间的差距有多大呢? chao1 指数给出的答案是 (F 1 ^2)/(2*F 2 ) ,它通过 singletons 和 doubletons 进行了合理的推算。分析 chao1 指数的后半段 (F 1 ^2)/(2*F 2 ) 我们不难发现它对 singletons 的权重要高于 doubletons ( 即 F 1 ^2 比 2*F 2 变化的速度更快 ) 。 Chao1 指数是基于这样一种假设:在一个群体中随机抽样,当稀有的物种 (singletons) 依然不断的被发现时,则表明还有一些稀有的物种没有被发现;直到所有物种至少被抽到两次 (doubletons) 时,则表明不会再有新的物种被发现。 ( The idea behind the estimator is that if a community is being sampled, and rare species (singletons) are still being discovered, there is likely still more rare species not found; as soon as all species have been recovered at least twice (doubletons), there is likely no more species to be found. 〔2〕) 综上, chao1 是度量物种丰富度的指标,它和丰度、均匀度无关,但是它对稀有的物种很敏感。 ( 这也正是丰富度指标应该具有的特性! ) 4. 举例 菌落 A ,有 50 个 OTUs ,其中仅有两条 reads 的 OTUs 有 10 个,仅有一条 read 的 OTU 有 12 个,那么其 chao1 指数值为 50+(12^2)/(2*10) = 57.2 菌落 B ,有 50 个 OTUs ,其中仅有两条 reads 的 OTUs 有 20 个,仅有一条 read 的 OTU 有 5 个,那么其 chao1 指数值为 50+(5^2)/(2*20) = 50.625 菌落 C ,有 80 个 OTUs ,其中仅有两条 reads 的 OTUs 有 10 个,仅有一条 read 的 OTU 有 12 个,那么其 chao1 指数值为 80+(12^2)/(2*10) = 87.2 菌落 D ,有 40 个 OTUs ,其中仅有两条 reads 的 OTUs 有 3 个,仅有一条 read 的 OTU 有 20 个,那么其 chao1 指数值为 40+(20^2)/(2*3) = 106.6667 ! ( Amazing! ) 可以将菌落 B 、菌落 C 和菌落 D 分别是茵落 A 进行比较,更加感性的去认识这一指数,这里就不再 赘述 了。 参考材料 http://scikit-bio.org/docs/0.4.1/generated/generated/skbio.diversity.alpha.chao1.html#skbio.diversity.alpha.chao1 http://palaeo-electronica.org/2011_1/238/estimate.htm Chao, A. 1984. Non-parametric estimation of the number of classes in a population. Scandinavian Journal of Statistics 11, 265-270.
个人分类: Metagenomics|70987 次阅读|2 个评论
Alpha多样性指数之shannon指数
热度 3 luria 2017-8-3 20:10
shannon 指数可以回答样本各物种均一性如何。 它的计算公式如下: 公式中 s 即 OTU 数目, pi 即第 i 个 OTU 的占比。 例如有三组样本: A 组: OTU1, OTU2, OTU3, OTU4 分别有 3, 4, 6, 7 条序列 B 组: OTUa, OTUb, OTUc, OTUd 均有 5 条序列 C组: OTUm, OTUn , OTUx, OTUy, OTUz分别有2, 3, 4, 5条序列 各组 OTU 相互独立,即 A 组的 OTU1 与 B 组的 OTUa 是否相同,无关紧要。因此可以略写为 A 组: 3, 4, 6, 7 B 组: 5, 5, 5, 5 C 组: 2, 3, 4, 5 AB 两组数均值都是 5 ,总数均为 20 , C 组的均值,总数均为最小 那么按照公式: A 组的 shannon 熵值为 : H A = - (3/20*log2(3/20)+4/20*log2(4/20)+6/20*log2(6/20)+7/20*log2(7/20)) = 1.93 ( 保留两位小数,下同 ) B 组的 shannon 熵值为: H B = - (5/20×log2(5/20) × 4) = 2.00 C 组的 shannon 熵值为: H C = - (2/14*log2(2/14)+3/14*log2(3/14)+4/14*log2(4/14)+5/14*log2(5/14)) = 1.92 样本中各 OTU 数越均一, shannon 指数越大。 当样本中所有 OTU 数相同时,达到极值。 注意: shannon 指数与 OTU 总数无关,从 B 组和 C 组可以看出来。它只和 OTU 数目的均一性有关! 带入公式,可以算出其极值 : 其中, u 是均值,即每个 OTU 都相同时的数目, N 是总的数目,如上例 B 组 H= - log2(5/20) = 2 ,可以快速算出 当然,以上还只能算是猜测,仅靠举例只能打开思路,并不严谨。 为了验证上面这个猜测,我们可以从公式入手。 使用 R 画出一条 y= - x*log2(x) 函数图如下,看样子很可能是个递增的函数,事实也是如此。 如果不信可以自行求导,结果应该是: 因为 x 是 OTU 数目,为大于 1 的整数,因此函数单调递减,且斜率越来越小 注意:以上函数并不是 shannon 指数函数,它只是其中一部分。因此,需回到 shannon 指数函数,可以根据 ( y(a) + y(b) + y(c) ...y(i) )/(a+b+c+...i) y((a+b+c+...k)/2), 其中y = - x*log2(x),i与k可以不相等。可以看出shannon指数在各OTU数目均一时达到最大。 综上, shannon 指数越大,反映出样本均一性越好。 参考材料: http://scikit-bio.org/docs/latest/generated/generated/skbio.diversity.alpha.shannon.html#skbio.diversity.alpha.shannon
个人分类: Metagenomics|54612 次阅读|3 个评论
QIIME中jackknifed_beta_diversity.py程序的运行过程
luria 2017-8-1 20:40
1. Beta Diversity (weighted_unifrac, unweighted_unifrac) beta_diversity.py -i pick_otus/otu_table.biom -o jackknifed/unrarefied_bdiv -t pick_otus/rep_set.tre 从所有 OTU table 和 tree 中计算 beta 多样性距离矩阵 运行完在 jackknifed/unrarefied_bdiv/ 文件夹下生成两个文件: unweighted_unifrac_otu_table.txt 和 weighted_unifrac_otu_table.txt 。即用加权 unifrac 和非加权 unifrac 计算出距离距阵 2. Rarefaction multiple_rarefactions_even_depth.py -i pick_otus/otu_table.biom -d 5000 -o jackknifed//rarefaction/ 建立稀释的 OTU table ,这一步与上一步无关。这里是在 jackknifed/rarefaction/ 文件夹下生成 10 次 jackknifed 抽样 3. UPGMA on full distance matrix: weighted_unifrac upgma_cluster.py -i jackknifed/unrarefied_bdiv/weighted_unifrac_otu_table.txt -o jackknifed/unrarefied_bdiv/otu_table_weighted_unifrac_upgma.tre 利用第 1 步中 weighted_unifrac 生成的 beta 多样性矩阵建立 UPGMA 树(在 unrarefied_bdiv 中生成了一个 otu_table_weighted_unifrac_upgma.tre ),可以理解为一个引导树。 4. Beta diversity on rarefied OTU tables (weighted_unifrac) beta_diversity.py -i jackknifed//rarefaction/ -o jackknifed//weighted_unifrac//rare_dm/ -m weighted_unifrac -t pick_otus/rep_set.tre 利用第 2 步中 10 次抽样生成的 biom 文件,生成 beta 多样性距阵。在 jackknifed/weighted_unifrac/rare_dm/ 文件夹下生成了 10 个矩阵 5. UPGMA on rarefied distance matrix (weighted_unifrac) upgma_cluster.py -i jackknifed//weighted_unifrac//rare_dm/ -o jackknifed//weighted_unifrac//rare_upgma/ 再利用上一步中的 10 个矩阵生成了 10 棵树文件。在 jackknifed/weighted_unifrac/rare_upgma/ 下生成了 10 个树文件 6. consensus on rarefied distance matrices (weighted_unifrac) consensus_tree.py -i jackknifed//weighted_unifrac//rare_upgma/ -o jackknifed//weighted_unifrac//rare_upgma_consensus.tre 再利用上一步中的 10 棵树生成了 1 棵一致树 7. Tree compare (weighted_unifrac) tree_compare.py -s jackknifed//weighted_unifrac//rare_upgma/ -m jackknifed//weighted_unifrac//rare_upgma_consensus.tre -o jackknifed//weighted_unifrac//upgma_cmp/ 生成了 upgma_cmp 文件夹下的文件,包括两个树文件和一个 jackknife 支持率的文件,其中这两个树文件是一致的,只是 master_tree.tre 文件的每个内节点保证有唯一的编号(来自官网上的解释)。 如果要做 UPGMA 树图,使用的 tre 文件就即这个 master_tree.tre 文件。这在官网也有指出,另外 QIIME 还提供了一个 make_bootstrapped_tree.py 程序使用这个 master_tree.tre 文件以及支持率文件来绘制 bootstrapped tree 。 如下图: 8. Principal coordinates (weighted_unifrac) principal_coordinates.py -i jackknifed//weighted_unifrac//rare_dm/ -o jackknifed//weighted_unifrac//pcoa/ 利用第 4 步生成的 10 个距离矩阵文件生成了对应的 10 个 PCoA 分析数据,保存在 jackknifed//weighted_unifrac//pcoa/ 文件夹下 9. emperor plots (weighted_unifrac) make_emperor.py -i jackknifed//weighted_unifrac//pcoa/ -o jackknifed//weighted_unifrac//emperor_pcoa_plots/ -m mapfile.txt 这一步生成了 jackknifed//weighted_unifrac//emperor_pcoa_plots/ 文件夹下的 PCoA 3D 网页文件,里边的内容可以调整。和 beta_diversity_through_plots.py 程序生成的 PCoA 3D 图最大的差别是点周围有一片区域,这片区域在 -e 设置越大小模糊区域越小 例如:分别用 3 种不同的 e 值( 100 , 1000 , 5000 )来跑, 7 个样品, OTU 数平均为 13500 。需要注意的是因为 jackknifed 是基于抽样的方法,所以同一条代码每次运行的结果可能都不一样。 e 为 100 的如下: e 为 1000 的如下: e 为 5000 的如下: ===========================我是一条华丽的分隔线,请勿注释掉=============================== 以下 unweighted_unifrac 方法分析得到的结果和上面 weighted_unifrac 一样,具体解析不再赘述 10. UPGMA on full distance matrix: unweighted_unifrac command upgma_cluster.py -i jackknifed/unrarefied_bdiv/unweighted_unifrac_otu_table.txt -o jackknifed/unrarefied_bdiv/otu_table_unweighted_unifrac_upgma.tre 11. Beta diversity on rarefied OTU tables (unweighted_unifrac) command beta_diversity.py -i jackknifed//rarefaction/ -o jackknifed//unweighted_unifrac//rare_dm/ -m unweighted_unifrac -t pick_otus/rep_set.tre 12. UPGMA on rarefied distance matrix (unweighted_unifrac) command upgma_cluster.py -i jackknifed//unweighted_unifrac//rare_dm/ -o jackknifed//unweighted_unifrac//rare_upgma/ 13. consensus on rarefied distance matrices (unweighted_unifrac) command consensus_tree.py -i jackknifed//unweighted_unifrac//rare_upgma/ -o jackknifed//unweighted_unifrac//rare_upgma_consensus.tre 14. Tree compare (unweighted_unifrac) command tree_compare.py -s jackknifed//unweighted_unifrac//rare_upgma/ -m jackknifed//unweighted_unifrac//rare_upgma_consensus.tre -o jackknifed//unweighted_unifrac//upgma_cmp/ 15. Principal coordinates (unweighted_unifrac) command principal_coordinates.py -i jackknifed//unweighted_unifrac//rare_dm/ -o jackknifed//unweighted_unifrac//pcoa/ 16. emperor plots (unweighted_unifrac) command make_emperor.py -i jackknifed//unweighted_unifrac//pcoa/ -o jackknifed//unweighted_unifrac//emperor_pcoa_plots/ -m mapfile.txt
个人分类: Metagenomics|6035 次阅读|0 个评论
rank abundance 解析
热度 1 luria 2017-8-1 20:16
在16S rRNA扩增子分析中,rank abundance已经是必备的一项分析内容。它可以从OTU的层面总体的反映出物种的分布情况(丰度和均匀度)。如下图: 大多数资料只会对横纵坐标进行解释: x轴反映的是各物种的丰度,物种丰度越高,样本落在x轴上的区间越大,如上图中Group5的物种丰度高于Group1; 整体来看,如果这条曲线下降得越平缓,则表明物种分布越均匀,上图中各物种的均匀程序相差不大。如果遇到某个样本生成的rank abundance曲线中,楼梯下得比较徒,即垂直的线越长表明这个OTU的丰度越大。 但是,rank abundance是怎样算出来的?为什么横纵坐标表示的是上述解释?一般都没有讲出来,这里给出rank abundance的算法,加深理解,理解了就可以自己动手来画。 rank abundance计算过程: 1. 获取每个样本中OTU的丰度值,例如从QIIME的otu table中获取。通俗来讲,即搞清楚样本(如上图中Group1),中每个OTU它代表了多少条序列。 2. 将每个样本中OTU的丰度值按照从大到小顺序进行排序。注意:这步最精华的部分,因此整个分析过程叫做rank abundacne! 3. 将每个样本归一化,获取每个样本OTU的相对丰度。例如样本1中有5个OTU,丰度分别为5,4,3,3,5,那么排序后为5,5,4,3,3,总丰度即sum(5,5,4,3,3)=20,归一化后的相对丰度为5/20=0.25, 5/20=0.25, ... 。 注意:样本的相对丰度和其它样本无关,例如样本1中有5个OTU,样本2中有7个OTU,那么样本1的sum值为这5个OTU丰度的总和,样本的sum值为7个OTU丰度的总和。 4. 对于每个样本, x轴是OTU的种类数,y轴是对应的OTU的丰度。例如样本1中有5个OTU,丰度分别为0.2,0.25,0.25,0.2,0.15,0.15(注意:这里是排好序的相对丰度)。那么按点(1,0.25),(2,0.25),(3,0.2),(4,0.15),(5,0.15)画出梯形图, 一般画的都是下梯形(和正常的楼梯相似)。 5. 一般会将y轴取对数,要不然会相当难看。同时将所有样本y轴的刻度调整到一致。这样就完成了rank abundance图的绘制。 对于微生物多样性,人们最想了解的是两大关键指标: 1. 某个生境下,一个菌落中,有多少种菌(物种多样性)? 2. 每种菌含量是多少(物种丰度)? rank abundance曲线能够从OTU的层面上很好的回答这两个问题 : 1. 在rank abundance曲线中每一条垂直的线代表一个OTU,可以理解为一种菌 2. 每一条垂直线右侧的水平线(即楼梯的平台)表示该菌种含量
个人分类: Metagenomics|12846 次阅读|1 个评论
QIIME安装及设置
luria 2017-5-21 11:25
在微生物基因组扩增子测序领域,QIIME可谓是一道标杆。QIIME全称Quantitative Insights Into Microbial Ecology,主要由美国科罗拉多大学Rob Knight实验室的一帮牛人开发的。编程语言主要为Python,官方提供的引文为 J Gregory Caporaso, Justin Kuczynski, Jesse Stombaugh, et al. QIIME allows analysis of high-throughput community sequencing data. Nature Methods, 2010; doi:10.1038/nmeth.f.303 本篇仅介绍QIIME的安装,后期博文会尝试从代码深入理解QIIME每步分析过程,敬请期待。 像QIIME这类Pipeline通常程序众多,安装升级牵扯的相联关系比较复杂,这里建议直接使用QIIME VirtualBox。 1. 下载安装VirturalBox,地址:https://www.virtualbox.org/ 下载QIIME,地址:http://qiime.org/install/virtual_box.html 文件大小约大于4G,解压后约51G。 2. 解压后为QIIME-1.9.1-amd64.vdi文件,将其放到某处(注意此位置以后不能调整)。 安装好VirtualBox后,点击常用工具栏的新建图标 = 在虚拟电脑名称和系统类型中输入名称(Qiime)类型(Linux)版本Ubuntu(64-bit) = 下一步 (如图1) = 内存大小对话框中选绿色标线最大处的内存(这里以后也可以调) = 下一步 = 在新建虚拟电脑对话框中选择“使用已有的虚拟硬盘”,并加载上QIIME-1.9.1-amd64.vdi(如图2)= 创建 图1 图2 3. 装好后界面如下: 按Ctrl+D安装增强工具,这时会自动弹出安装增强工具,如果不小心关掉了,可以打开文件资源管理器,进入VBOXADDITIONS(如下图左栏)进入加载的增强工具盘 双击autorun.sh # 如果想用命令行进入VBOXADDITIONS,其位置在/media/下面,需要 sudo chmod 755 * 一路运行如下,完成后回车 之后再重启一下 4. 安装好增强功能后,做以下设置: 4.1 调整屏幕分辨率: 单击右上角的设置按钮(齿轮)- Displays,我这里选1920*1080(16:9) 4.2 修改QIIME的主机名 打开终端 sudo vim /etc/hostname 然后将其中的名字修改一下,即可 完了之后需要重启才能生效 4.3 在虚拟机中共享主机文件夹 先在VirtualBox虚拟机中开启主机文件夹共享,如下图 其中luria是我的主机用户名,你可以在Linux下输入 whoami 命令查看你的主机用户名,将以下步骤中的luria换成你的主机用户名。 这时df -h可以查看到已加载上主机共享的文件夹如下图,但是这里没有权限打开。 需将luria用户加入到vboxsf组中,才有权限访问挂载的目录,以实现文件共享! 执行: sudo adduser scg vboxsf 然后sudo reboot重启生效
个人分类: Metagenomics|8005 次阅读|0 个评论
Linux下QIIME的安装经历
xbinbzy 2015-11-5 15:52
随着16s rRNA的研究越来越受到科研工作者的关注,Mothur和QIIME作为这个领域内应用较为广泛的工具,引用率也越来越高。之前在服务器上(Linux)系统安装好了mothur,具体安装过程可参考 http://blog.sciencenet.cn/blog-306699-882982.html 此文主要用于介绍自己在服务器安装QIIME所遇到的各种问题(说起来全是泪,安装过程中各种bug,各种google,好多解决方案不知道如何判断。以上记录的是自己服务器上安装的整个过程,测试成功完成); 按照QIIME安装说明文档: http://qiime.org/install/install.html 1)要求2.7或3的python版本,于是用以下命令查看了一下 发现版本为2.7.10,于是进行下一步。 2)由于QIIME有比较多的包,安装建议用pip: (1)在python中import pip试了一下发现没有,于是安装pip开始,在 https://pypi.python.org/pypi/pip 下载好pip-7.1.2.tar.gz,解压后python setup.py install进行安装,出现报错信息: File /tmp/tmpru8qrk/pip.zip/pip/__init__.py, line 10, in module File /tmp/tmpru8qrk/pip.zip/pip/util.py, line 18, in module File /tmp/tmpru8qrk/pip.zip/pip/_vendor/distlib/version.py, line 14, in module File /tmp/tmpru8qrk/pip.zip/pip/_vendor/distlib/compat.py, line 66, in module ImportError: cannot import name HTTPSHandler google找到解决方案 http://jingyan.baidu.com/article/e52e3615aba39640c60c51c3.htm l yum install openssl -y和yum install openssl-devel -y ,然后重新安装pip。如此pip的安装就成功了。 (2)运行pip install numpy 发现python中有安装好的numpy。于是运行pip install qiime,这时提示numpy有错。然后import numpy,发现报错(没有记录好报错信息),Google了一下没找到解决方案,于是重装。重装过程以为删掉python和包的目录就可以,最后发现是需要/usr/local/lib/python2.7/site-packages此目录中删除numpy的目录。 安装好pip好,用pip install numpy,这下搞定。 (3)运行pip install qiime 这个运行了一段时间,最后出现报错信息 Hash of the package https://pypi.python.org/packages/source/q/qiime/qiime-1.9.1.tar.gz#md5= 32d4db2b2c888dbbfc6aa33f02725526 (from https://pypi.python.org/simple/qiime/) (c1147f0d3be7128bdf35c3f2e54d1e40) doesn't match the expected hash 32d4db2b2c888dbbfc6aa33f02725526! Bad md5 hash for package https://pypi.python.org/packages/source/q/qiime/qiime-1.9.1.tar.gz# md5=32d4db2b2c888dbbfc6aa33f02725526 (from https://pypi.python.org/simple/qiime/) goolge查找各种解决方案,提示这个下载的文件不完整导致。于是重新运行pip install qiime,发现没有其他的信息输出,直接报上述的错,并提示 Using cached qiime-1.9.1.tar.gz ,这是用了上一次操作过后的缓存文件。于是修改pip install qiime,为pip --no-cache-dir install qiime,重新正常下载和安装,这个过程需要保持网速稳定。 (4)pip --no-cache-dir install qiime可以运行 不过在安装scipy这个包时,出现了以下报错信息: lapack_opt_info: openblas_lapack_info: libraries openblas not found in NOT AVAILABLE lapack_mkl_info: mkl_info: libraries mkl,vml,guide not found in NOT AVAILABLE Google找解决方案,是lapack和atlas的包没有在系统环境安装导致。于是下载lapack( http://www.netlib.org/lapack/ )和atlas的包( http://sourceforge.net/projects/math-atlas/ )。安装策略参考 在这步操作中,最早下载atlas3.10.2.tar.bz2安装时,出现以下报错信息: ERROR: enum fam=3, chip=2, model=62, mach=0 make : Error 44 make : Error 2 最终google没有找到解决方案,看到资料有人提示不同版本要求会不一致,为此下载atlas3.11.38.tar.bz2进行安装,安装过程中报错信息为:unrecognized command line option -mavx2,针对此问题google找到解决方案 http://ask.xmodulo.com/upgrade-gcc-centos.html 提示进行安装,最终搞定atlas。 (5)重新 pip --no-cache-dir install qiime,顺利安装 (6)用print_qiime_config.py -t测试时,发现出现报错信息“ ImportError: /usr/local/lib/python2.7/site-packages/scipy/spatial/qhull.so: undefined symbol: _gfortran_transfer_integer_write“ ,google后发现有人提出 接着这个提示,开始卸载ATLAS,安装OpenBLAS,在安装OpenBLAS的过程中,出现报错提示gcc的版本不够,于是利用(4)步中升级好的gcc进行操作。 (7)重新运行 pip --no-cache-dir install qiime,顺利安装成功后测试 print_qiime_config.py -t,显示信息: 在这步中,在pip --no-cache-dir install qiime中刚开始还遇到“ Command /usr/local/bin/python -c import setuptools, tokenize;__file__='/tmp/pip-buildNJHPFa/scipy/setup.py';exec(compile(getattr (tokenize, 'open', open)(__file__).read().replace('\r\n', '\n'), __file__, 'exec')) install --record /tmp/pip-cK8VSa-record/install-record.txt --single-version-externally-managed --compile failed with error code 1 in /tmp/pip-build-NJHPFa/scipy ”报错信息 ,调整把gcc环境调到(4)步所设计好的,重新 pip --no-cache-dir install qiime,成功安装。然后测试成功。 这个软件的安装过程是做生物信息这几年来,遇到问题最多的一个工具包。问题不断,然后google不断,不过现在信息太多,找个每个问题的解决方案各家有不同之说,于是不断尝试,直至找到成功的。
个人分类: 数据分析|22907 次阅读|0 个评论
16s rRNA分析流程和工具的介绍
热度 2 xbinbzy 2015-11-3 23:08
16s rRNA早期的分析策略,如FISH(fluorescent in situ hybridization)、DDGE(denaturing gradient gel electrophoresis)、PCR cloning、T-RFLP(terminal restriction fragment length polymorphism)。随着NGS(next generation sequencing)测序技术的发展,在此主要讨论NGS技术在16s rRNA分析中的应用。 16s rRNA NGS数据分析的主要工具有: 16s rRNA NGS数据的分析主要有3个大步骤: 1)pretreatment of raw sequence data:including de-multiplexing of barcoded sequences, quality filtering, denoise, chimera checking, and data normalization; 2)microbial diversity analysis:construction of operational taxonomic unit (OTU)table, including picking OTUs, picking representative sequences, aligning representative sequences, taxonomic assignment of representative sequences, and building of a phylogenetic tree of the OTUs; 3)advanced data analysis and visualization:including the alpha- and beta-diversity analyses,clustering and coordinates analysis, and data visualization; (1)De-multiplexing and quality filtering 16s经常是pooling测序,为此需要将下机数据根据barcode序列信息将数据拆分到各样品中。QIIME中的 “ split_libraries.py ” 和“ split_libraries_fastq.py ”实现数据拆分和数据过滤的双重目的。Mothur利用“ Trim.seqs ”。不过QIIME和Mothur都不能直接处理sff文件(454下机产生的数据格式),不过可各自利用“process_sff.py”和Sffinfo将sff格式转换为FASTA和QUAL文件。 数据过滤考虑的参数有:minimum average quality score allowed in a read、maximum number of ambiguous bases allowed、minimum and maximum sequence length、maximum length of homopolymer allowed、maximum mismatches inprimer or barcode allowed、whether to truncate reverse primer, and so on. (2)Denoise and chimera checking 16s建库的pcr过程、测序过程均会导致序列出现错误,在分析过程过程中需要有效排除这种错误。测序误差矫正常用的工具有 Denoiser(implemented in QIIME)、AmpliconNoise、Acacia、Pre.cluster(implemented in Mothur) 。嵌合体查找的工具有 ChimeraSlayer、UCHIME、Persus、DECIPHER , ChimeraSlayer、UCHIME、Persus 在mothur中均可调用。 在这些工具中,存在有待于优化的问题 (these different methods often disagree with one another on the list of identified chimeras,probably because of their different mechanisms or algorithms. More efforts are required to evaluate these methods and coordinate their inconsistencies in chimera identification.) 在分析中有个关于古细菌序列的情况需要注意: a very small proportion of archaeal sequences may be generated for 16S rRNA gene amplicon datasets amplified with bacteria-specific primers. These unexpected sequences should be identified after denoising and chimera removal, and are advised to be discarded before subsequent data normalization. (3)Data normalization 测序深度不理想和不均匀的话会对alpha多样性及beta多样性均有影响。Uneven sequencing depth can affect diversity estimates in a single sample (i.e.,alpha diversity), as well as comparisons across different samples (i.e., beta diversity), thus data normalization is required. 对于此问题有两种处理策略,分别是relative abundance and random sampling (i.e., rarefaction),in addition, z-score亦用于normalization的过程中。但不同的方法均会有缺点。 (4)Picking OTUs and representative sequences OTU的界定主要根据序列的一致性进行,(The OTUs are picked based on sequence identity, and various identity cutoffs of 16S rRNA gene have been used for different taxonomic ranks. For example, identity cutoffs recommended by MEGAN are 99 % for species, 97 % for genus,95 % for family, and 90 % for order level, respectively)。OTU界定时选择的工具与算法对后期分析有很大影响(The OTU picking strategy and algorithms have significant effects in the downstream data interpretation. )根据分析过程中是否使用数据库,OTU界定的策略可分为de novo、closed reference和open reference。在 OTU界定中有很多聚类的方法,There are many clustering or alignment tools available for OTU picking, such as Uclust, cd-hit, BLAST, mothur, usearch, and prefix/suffix. These tools are implemented in QIIME. Among them, the mothur method contains three clustering algorithms to pick de novo OTUs, namely, nearest neighbor, furthest neighbor, or average neighbor. 当序列聚类好后,代表了一个OTU,接下来就是从这个OTU找到代表序列,一种做法是a representative sequence can be a random, the longest, the most abundant(as default in QIIME), 另一种操作方法是the first sequence in an OTU cluster. 还有一种策略是the distance method in mothur identifies the sequence with the smallest maximum distance to the other sequences as the representative sequence. (5)Taxonomic assignment taxonomic assignment的策略有:(1)word match,如RDP classfier,(2)best hit,(3)Lowest Common Ancestor,如MEGAN、SINA Alignment Service. (6 )Phylogenetic analysis Phylogenetic relationships一般可以用树来表示,phylogenetic relationships主要是通过序列比对来实现的,序列比对的工具有ClustalW, MUSCLE, Clustal Omega, Kalign, T-COFFEE, COBLAT和FastTree. 目前针对16s rRNA NGS数据的分析工具都可以实现,如MEGA,RAxML,MRBAYES,PhyML,TreeView,Clearcut,FitTree. 其中RAxMLand PhyML are the most widely used programs for maximum-likelihood phylogenetic analysis, probably because they are specifically designed and optimized for such purpose. (7)Alpha- and beta-diversity analyses alpha多样性有众多指标可以表示,在mothur中有Shannon, Berger-Parker,Simpson, Q statistic; observed richness, Chao1, ACE, and jackknife。而在QIIME中,有phylogenetic diversity (PD)-whole tree, chao1, and observed species. 还有一种物种丰度的比较方法:rarefaction curve. QIIME中主要用“single_rarefaction.py”、 “multiple_rarefaction.py”,在mothur中主要用“Rarefaction.single”和“Rarefaction.shared”. beta多样性计算主要反映不同样本之间的差异度,several distance metrics, such as Unifrac, Bray-Curtis, Euclidean,Jaccard index, Yue Clayton, and Morisita-Horn, have been often employed. beta多样性计算根据是否考虑OTU的相对丰度,可分为定量指数和定性指数。 ( 8)Statistical and network analysis 在Two-sample/group中,多考虑 t-test。 在其中需要注意, Particularly for independent two-sample t -test, independence and equal variances (which canbe tested by F-test, Levene ’ s test, etc.) of two populations arerequired. In the case of non-normal distribution of data sets,nonparametric two-sample tests robust to data non-normality,such as Wilcoxon signed-rank test, and Mann-Whitney U testare applicable for significance testing of difference betweengroup medians. 在Multiple-sample/group tests中,ANOVA。 (9)Clustering and classification clustering可以分析样品之间的亲疏关系。classfication的策略用来对样品进行类别判定。 (10)Ordination analysis 在样本的相似度和距离计算完后,可以利用 principal component analysis (PCA), principal coordinates analysis(PCoA, also known as metric multidimensional scaling), Nonmetric multidimensional scaling (NMDS), canonical correspondence analysis (CCA), linear discriminantanalysis (LDA), and redundancy analysis (RDA)等构建样本间的关系。 (11)Network-based modeling 与基因表达、代谢分子、蛋白等数据一起分析共表达网路或者共表达模式(co-occurrence and co-exclusion patterns) 参考文章: Ju F, Z hang T. 16s rRNA gene high throughput sequencing data mining of microbiota diversity and interactions, Appl Microbiol Biotechnol . 2015, 99(10):4119-4129
个人分类: 科研文章|35820 次阅读|3 个评论

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-6-18 19:31

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部