如何读懂和利用你的微生物多样性测序结果?

16s分析结果详解

很多小伙伴有过这样的经历,在拿到公司出具的报告之后,仍然一头雾水,几十页的报告内容看着丰富却不知该怎么运用。我们一起来理一下关键图表的含义

OTU是我们要搞清的一个重要概念,可以说是后续分析的基石。

OTU(operational taxonomic units) 是在系统发生学研究或群体遗传学研究中,为了便于进行分析,人为给某一个分类单元(品系,种,属,分组等)设置的同一标志。通常按照 97% 的相似性阈值将序列划分为不同的 OTU,每一个 OTU 通常被视为一个微生物物种。相似性小于97%就可以认为属于不同的种,相似性小于93%-95%,可以认为属于不同的属。样品中的微生物多样性和不同微生物的丰度都是基于对OTU的分析。

有了OTU这个概念之后,就不难理解下表。对每个样本的测序数量和OTU数目进行统计,并且在表栺中列出了测序覆盖的完整度。

其中 SampleName表示样本名称;SampleSize表示样本序列总数;OTUsNumber表示注释上的OTU数目;OTUsSeq表示注释上OTU的样本序列总数。

Coverage是指各样品文库的覆盖率,其数值越高,则样本中序列没有被测出的概率越低。该指数实际反映了本次测序结果是否代表样本的真实情况。计算公式为:C=1-n1/N  其中n1 = 只含有一条序列的OTU的数[……]

Read more

Random Forests 隨機森林 | randomForest, ranger, h2o | R語言 (zhuantie)

https://www.jamleecute.com/random-forests-%E9%9A%A8%E6%A9%9F%E6%A3%AE%E6%9E%97/

Bagging法綜合多個樹模型結果,可以降低單一樹模型的高變異性並提升預測正確率。但Bagging法中樹與樹之間的相關性會降低模型整體的表現。隨機森林 Random forests 是Bagging修改後的版本,它是由「去相關性」的樹模型所組成的集成演算法,有很不錯的預測正確率且是一個受歡迎、開箱即用的演算法。

載入所需套件

Read more

Bash script for batch work

#!/bin/bash
#PBS -l walltime=96:00:00,nodes=8:ppn=10
#PBS -N kofamscan

source activate hybrid_assembly
#SC9523_extract_paired_1.fasta.fastq_matchedmeghitresult
#metawrap binning -o /home/mhyleung/workspace/loreal_shotgun2/NovaSeq_new/metawrap_assembly_R2/LOR303C -t 16 -a /home/mhyleung/workspace/loreal_shotgun2/NovaSeq_new/metawrap_assembly_R2/LOR303C_S134_R1_00
#1_kneaddata_paired/LOR303C_final_assembly.fasta –maxbin2 –metabat2 –concoct /disk/rdisk08/mhyleung/loreal_shotgun2/NovaSeq_new/decontam_output_genome/paired_completed/szy_test/LOR303C_S134_R1_001_knead data_paired_*.fastq[……]

Read more

MetaGEM:直接从宏基因组重建基因组规模的代谢模型

https://github.com/franciscozorrilla/metaGEM

A Snakemake-based workflow to generate high quality metagenome assembled genomes from short read paired-end data, reconstruct genome scale metabolic models, and perform community metabolic interaction simulations on high performance computing clusters.

103545679-ceb71580-4e99-11eb-9862-084115121980

基因组规模代谢网络模型(Genome-scale metabolic model,GEM),是一种包含了某种特定生物或者是细胞基因组范围代谢反应,及其酶及基因关联的数学模型

这里,我们基于文章的描述,介绍一款新软件——MetaGEM。

研究者认为,目前代谢建模的工作流程仍然是倾向于依赖参考基因组作为重建和模拟GEMs的起点,这忽略了微生物群落中存在的物种内和物种之间的多样性。也限制了对已知参考基因组空间中的代谢网络的分析和解释。

可能导致假阳性(即在参考基因组中存在但在群落中的变量中缺失的通路)或假阴性(即在参考基因组中缺失但在群落变量中存在的通路)结果,最终导致对个别[……]

Read more

「三代组装」使用Pilon对基因组进行polish

链接:https://www.jianshu.com/p/cceeb7d1f413来源

软件安装

如果要顺利运行程序,要求JAVA > 1.7, 以及根据基因组大小而定的内存,一般而言是1M大小的基因对应1GB的内存。

总览

Pilon有如下作用

  1. 对初步组装进行polish
  2. 寻找同一物种不同株系间的变异,包括结构变异检测

他以FASTA和BAM文件作为输入,根据比对结果对输入的参考基因组进行提高,包括

  • 单碱基差异
  • 小的插入缺失(indels)
  • 较大的插入缺失或者block替换时间
  • 填充参考序列中的N
  • 找到局部的错误组装

最后它输出polish后的FASTA文件, 以及包含变异信息的VCF文件(可选)

分析流程

推荐使用PCR-free建库测序得到的Illumina paired-end数据,这样子避免了PCR-duplication,有效数据更多,也不需要在分析过程中标记重复。

下面步骤,假设你的组装文件为draft.fa, 质量控制后的illumina双端测序数据分别为read_1.fq.gzread_2.fq.gz

第一步:比对

bwa index -p index/draft draft.fa
bwa mem -t 20 index/draft re[......]

Read more

wilcox 检验批处理示例

https://cloud.tencent.com/developer/article/1666764
#GRID_input_new_p10 group group
library(doBy) #使用其中的 summaryBy() 以方便按分组计算均值、中位数
 
#读取数据
gene <- read.table('gene.txt', sep = '\t', row.names = 1, header = TRUE, stringsAsFactors = FALSE, check.names = FALSE)
group <- read.table('group.txt', sep = '\t', header = TRUE, stringsAsFactors = FALSE, check.names = FALSE)
result <- NULL
 
#使用循环,逐一对各基因进行两组间 wilcox 秩和检验
for (n in 1:nrow(gene)) {
gene_n <- data.frame(t(gene[n,subset(group, group %in% c(1, 2))$sample]))
gene_id <- names(gene_n)[1]
names(gene_n)[1] <- 'gene'[......]

Read more

Get statistics for each group (such as count, mean, etc) using pandas GroupBy?


Quick Answer:

The simplest way to get row counts per group is by calling .size(), which returns a Series:

df.groupby(['col1','col2']).size()

Usually you want this result as a DataFrame (instead of a Series) so you can do:

df.groupby(['col1', 'col2']).size().reset_index(name='counts')

If you want to find out how to calculate the row counts and other statistics for each group continue reading below.

Detailed example:

Consider the following example dataframe:

In [2]: df
Out[2]: 
  col1 col2  col3  col4  col5  col6
0    A    B  0.20 -0.61 -0.49  1.49
1    A    B -1.53 -1.01 -[......]

Read more

Heatmap in R: Static and Interactive Visualization

https://www.datanovia.com/en/lessons/heatmap-in-r-static-and-interactive-visualization/

R Packages/functions for drawing heatmaps

There are a multiple numbers of R packages and functions for drawing interactive and static heatmaps, including:

  • heatmap() [R base function, stats package]: Draws a simple heatmap
  • heatmap.2() [gplots R package]: Draws an enhanced heatmap compared to the R base function.
  • pheatmap() [pheatmap R package]: Draws pretty heatmaps and provides more control to change the appearance of heatmaps.
  • d3heatmap() [d3heatmap R package]: Draws an interactive/clickable[……]

Read more

batch download genome based on accession

Three easy ways to download multiple sequences from NCBI

There are different ways of how to download multiple sequences from the NCBI databases in a single request.

 

1) Using the batch Entrez website

http://www.>ncbi.nlm.nih.gov/sites/batchentrez

 

2) Using Perl: (copy into your terminal and press return/enter)

perl -e 'use LWP::Simple;getstore("http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nucleotide&rettype=fasta&retmode=text&id=".join(",",qw(6701965 6701969 6702094 6702105 6702160)),"seqs.fasta");'

This takes the IDs separated by s[……]

Read more

生信小工具专题:BBTools/BBMap Suite 的使用 (zhuantie)

链接:https://www.jianshu.com/p/175c3282a61c

BBMap/BBTools是一种用于DNA和RNA测序reads的拼接感知全局的比对工具。它可以处理的reads包括,Illumina,454,Sanger,Ion Torrent,Pac Bio和Nanopore。 BBMap快速且极其准确,特别是对于高度突变的基因组或长插入读数,甚至超过100kbp长的全基因缺失。它对基因组大小或contigs数量没有上限。并且已经成功地用于绘制到具有超过2亿个contigs的85gb的土壤宏基因组。此外,与其他比对工具相比,索引阶段非常快。
BBMap有很多选项,在shell脚本中描述。它可以输出许多不同的统计文件,例如经验读取质量直方图,插入大小分布和基因组覆盖,有或没有生成sam文件。因此,它可用于文库的质量控制和测序运行,或评估新的测序平台。衍生程序BBSplit也可用于分类或精炼宏基因读取。

那问题来了,现在的工具那么多为什么会推荐还有写这个bbmap、bbTools的教程?那当然是因为该工具有一些它自身的特点:

  1. BBTools是用纯Java编写的,可以在任何平台(PC,Mac,UNIX)上运行,除了安装Java之外没有依赖项(为Java 6及更高版本编译)。所有工具都是高效且多线程的。
  2. BBTools可以处理常见[……]

Read more