熊金波实验室出品

整理归纳：Larry

本次学习使用的服务器IP地址和其用户名账户密码如下：

地址：gs0.genek.tv

用户名：u3276

密码：genek.tv
（建议不要用终端登陆，使用Xshell登陆较好）

数据质控

打开服务器你将会看到

u3276@GenekServer-0:~$

这是命令提示符。

在服务器上安装整个分析流程实施所需要的的软件

sudo apt-get -y update && \
sudo apt-get -y install trimmomatic python-pip \
   samtools zlib1g-dev ncurses-dev python-dev
sudo apt-get install -y python3.5-dev python3.5-venv make \
    libc6-dev g++ zlib1g-dev

 wget -c http://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.5.zip
 unzip fastqc_v0.11.5.zip
 cd FastQC
 chmod +x fastqc
 cd

现在，创建一个python 3.5虚拟环境并在其中安装软件:

python3.5 -m venv ~/py3
. ~/py3/bin/activate
pip install -U pip
pip install -U Cython
pip install -U jupyter jupyter_client ipython pandas matplotlib scipy scikit-learn khmer

pip install -U https://github.com/dib-lab/sourmash/archive/master.zip

让我们也运行Jupyter Notebook。首先，配置它更安全一点，并让它在后台运行:

jupyter notebook --generate-config

cat >>~/.jupyter/jupyter_notebook_config.py <<EOF
c = get_config()
c.NotebookApp.ip = '*'
c.NotebookApp.open_browser = False
c.NotebookApp.password = u'sha1:5d813e5d59a7:b4e430cf6dbd1aad04838c6e9cf684f4d76e245c'
c.NotebookApp.port = 8000

EOF

现在开始运行

jupyter notebook &

我们将会使用一批来自 [Hu et al., 2016]本文从班菲尔德实验室对一些多样性相对较低的环境进行了采样，利用一串近乎完整的基因组。以此为案例来进行这一次分析。

1.复制一些数据以便使用。

我已经为各位将数据的子集加载到Amazon位置，以使今天的分享变得更快一些。当然具体步骤也可以使用上一次分享中讲到的ftp（filezila）服务器数据上传下载工具。

mkdir ~/data

接下来，让我们在Linux系统下开始下载数据：

cd ~/data
curl -O -L https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1976948_1.fastq.gz
curl -O -L https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1976948_2.fastq.gz

这里涉及到一个数据是否完整的问题，如果你需要利用公开数据进行分析此部非常重要。运行数据诊断程序：

md5sum SRR1976948_1.fastq.gz SRR1976948_2.fastq.gz

你会看见

37bc70919a21fccb134ff2fefcda03ce  SRR1976948_1.fastq.gz
29919864e4650e633cc409688b9748e2  SRR1976948_2.fastq.gz

此时如果你输入陈列命令ls：

ls -l

你可以看见如下代码

total 346936
-rw-rw-r-- 1 ubuntu ubuntu 169620631 Oct 11 23:37 SRR1976948_1.fastq.gz
-rw-rw-r-- 1 ubuntu ubuntu 185636992 Oct 11 23:38 SRR1976948_2.fastq.gz

这是原始数据的1M读取子集，取自文件的开头。
这些文件的一个问题是它们是可写的——默认情况下，UNIX让文件所有者可以编辑。这就产生了在原始数据中创建拼写错误或错误的问题。在我们继续之前，让我们先解决这个问题:

chmod u-w *

2. FastQC

fastqc SRR1976948_1.fastq.gz
fastqc SRR1976948_2.fastq.gz

接下来输入ls

ls -d *fastqc.zip*

列出文件，你应该看到：

SRR1976948_1_fastqc.zip
SRR1976948_2_fastqc.zip

在每个fatqc目录中，您将找到来自fastqc的报告。可以使用Jupyter Notebook控制台下载这些文件;或者你可以看看它们的这些副本：

• SRR1976948_1_fastqc/fastqc_report.html
• SRR1976948_2_fastqc/fastqc_report.html

3. Trimmomatic

现在我们要做一些修剪！我们将使用Trimmomatic，它（和fastqc一样）我们已经通过apt-get安装了。

curl -O -L http://dib-training.ucdavis.edu.s3.amazonaws.com/mRNAseq-semi-2015-03-04/TruSeq2-PE.fa

我们需要的第一件事是修剪掉的接头和低质量序列：

TrimmomaticPE SRR1976948_1.fastq.gz \
              SRR1976948_2.fastq.gz \
     SRR1976948_1.qc.fq.gz s1_se \
     SRR1976948_2.qc.fq.gz s2_se \
     ILLUMINACLIP:TruSeq2-PE.fa:2:40:15 \
     LEADING:2 TRAILING:2 \
     SLIDINGWINDOW:4:2 \
     MINLEN:25

可以看见

...
Input Read Pairs: 1000000 Both Surviving: 885734 (88.57%) Forward Only Surviving: 114262 (11.43%) Reverse Only Surviving: 4 (0.00%) Dropped: 0 (0.00%)
TrimmomaticPE: Completed successfully

4.再一次 FastQC

fastqc SRR1976948_1.qc.fq.gz
fastqc SRR1976948_2.qc.fq.gz

• SRR1976948_1.qc_fastqc/fastqc_report.html
• SRR1976948_2.qc_fastqc/fastqc_report.html

5. MultiQC

现在，在工作目录中的未修剪和修剪文件上运行Mulitqc：

multiqc .

现在我们应该看到如下所示的输出：

[INFO   ]         multiqc : This is MultiQC v1.0
[INFO   ]         multiqc : Template    : default
[INFO   ]         multiqc : Searching '.'
Searching 15 files..  [####################################]  100%
[INFO   ]          fastqc : Found 4 reports
[INFO   ]         multiqc : Compressing plot data
[INFO   ]         multiqc : Report      : multiqc_report.html
[INFO   ]         multiqc : Data        : multiqc_data
[INFO   ]         multiqc : MultiQC complete

现在我们可以使用Jupyter Notebook查看输出文件。

K-mer 修正质控出现的错误

如果我们绘制样品k-mer丰度的柱状图，你会注意到存在大量的unqiue K-mers，即使测序质量很高，但它们也是由测序错误导致的。

序列末端低质量区有极高复杂度的k-mer

# 对质控前后的数据统计单端丰度距离
abundance-dist-single.py -M 1e9 -k 21 SRR1976948_1.fastq.gz SRR1976948_1.fastq.gz.dist

abundance-dist-single.py -M 1e9 -k 21 SRR1976948_1.qc.fq.gz SRR1976948_1.qc.fq.gz.dist

# 只对高覆盖度中的低丰度kmer剪切(更可能是测序错误)；低覆盖度保留
interleave-reads.py SRR1976948_1.qc.fq.gz SRR1976948_2.qc.fq.gz | trim-low-abund.py -V -M 8e9 -C 3 -Z 10 - -o SRR1976948.trim.fq

为什么要进行k-mer剪切

如果不做这步也是可以的。但会增加下游组装的工作量，本步可使结果更准确，并增加下游拼接速度，以及内存消耗。

unique-kmers.py SRR1976948_1.qc.fq.gz SRR1976948_2.qc.fq.gz
unique-kmers.py SRR1976948.trim.fq

结果如下

# 质控后的32-mers数据
Estimated number of unique 32-mers in SRR1976948_1.qc.fq.gz: 65344914
Estimated number of unique 32-mers in SRR1976948_2.qc.fq.gz: 85395776
Total estimated number of unique 32-mers: 112758982

# k-mer剪切后的数据
Estimated number of unique 32-mers in SRR1976948.trim.fq: 101285633
Total estimated number of unique 32-mers: 101285633

结果只经过了简单的尾部过滤，k-mer的数量减少了10%以上，对下游分析的准确度和速度都非常有帮助。
按Kmer质控后的结果，感兴趣的的再用fastqc评估一下，看看有什么变化？
接下来我们还会详细介绍利用k-mer的更到好处和实际意义。

MEGAHIT是一款非常快速，非常好的组装程序，专为宏基因组而设计。

首先，我们还是先进行安装

  cd ~/
   git clone https://github.com/voutcn/megahit.git
   cd megahit
   make

接着我们继续下载数据

  mkdir ~/data
   cd ~/data
   curl -O https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1976948.abundtrim.subset.pe.fq.gz
   curl -O https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1977249.abundtrim.subset.pe.fq.gz

这些是通过k-mer丰度修剪运行的数据（参见k-mer和subsampled，以便我们可以在相当短的时间内运行程序集）。

现在，最后，运行组装程序！

 mkdir ~/assembly
   cd ~/assembly
   ln -fs ../data/*.subset.pe.fq.gz .

   ~/megahit/megahit --12 SRR1976948.abundtrim.subset.pe.fq.gz,SRR1977249.abundtrim.subset.pe.fq.gz \
       -o combined

这将花费大约15分钟;最后你应该看到这样的输出：

 ... 7713 contigs, total 13168567 bp, min 200 bp, max 54372 bp, avg 1707 bp, N50 4305 bp
   ... ALL DONE. Time elapsed: 899.612093 seconds 

组装完成后

现在我们可以在装配上运行一些统计数据。为此，我们将使用QUAST：

 cd ~/
    git clone https://github.com/ablab/quast.git -b release_4.5
    export PYTHONPATH=$(pwd)/quast/libs/

现在，在程序集上运行QUAST：

 cd ~/assembly
    ~/quast/quast.py combined/final.contigs.fa -o combined-report
    cat combined-report/report.txt

到这一步，我们就已经完成了宏基因组数据分析的基础步骤，可以开始进行后续的数据分析，这个分析是多种多样的，可以根据研究方向的需求进行调整。

宏基因组数据分析流程

Prokka注释基因

细菌基因组、宏基因组的基因注释一直是一个非常复杂的问题，Prokka的出现改变了这一切。
Prokka: rapid prokaryotic genome annotation

Prokka

cd ~/
mkdir annotation
cd annotation

将metagenome程序集文件链接到此目录：

ln -fs ~/mapping/subset_assembly.fa .

现在是时候运行Prokka了！有许多不同的方法来专门运行Prokka。不过，我们现在要保持简单。

prokka subset_assembly.fa --outdir prokka_annotation --prefix metagG --metagenome

这将生成一个名为prokka_annotation的新文件夹，其中将包含一系列文件。

特别是，我们将使用* .ffn文件来评估我们的宏基因组中预测的基因组区域的相对读取覆盖率。

sourmash比较数据集教程

主页：https://sourmash.readthedocs.io/en/latest/
sourmash是一款超快、轻量级核酸搜索和比对软件，方法叫MinHash。这这一个命令行使用的python包，可以从DNA序列中快速分析kmer，并进行样品比较和绘图。目的是快速且准确的比较大数据集。
本步骤的目的：
比对reads至拼装结果
比较数据集并构建聚类图

采用k-mer Jaccard distance进行样品比较

什么是k-mer

k-mer就是长度为k的DNA序列

ATTG - a 4-mer
ATGGAC - a 6-mer

比如，一个长度为16的序列可以提取11个长度为6的 k-mers

AGGATGAGACAGATAG

AGGATG
 GGATGA
  GATGAG
   ATGAGA
    TGAGAC
     GAGACA
      AGACAG
       GACAGA
        ACAGAT
         CAGATA
          AGATAG

需要记住的概念：
不同物种的k-mer是很不同的
长k-mer具有很强的物种特异性

K-mers与组装图

三大基因组拼接方法之一，就是将基因组打断成k-mer再拼接，通过k-mer步移实现拼接
为什么采用k-mers而不是全长序列拼接
计算机喜欢k-mer，因为匹配准确快速。

图1. 以k=40为例，kmer很容易按属水平分开细菌

采用k-mers比较样品

安装sourmash
在前面内容的基础之上，只需执行下面代码的最后一条即可

# 设置工作目录，这是我的目录，学习时请修改为你工作的目录
pwd=~/test/metagenome17
cd ${pwd}

# 依赖关系，之前安装成功可跳过
sudo apt-get -y update && \
sudo apt-get install -y python3.5-dev python3.5-venv make \
    libc6-dev g++ zlib1g-dev
. ~/py3/bin/activate
pip install -U pip
pip install -U Cython
pip install -U jupyter jupyter_client ipython pandas matplotlib scipy scikit-learn khmer

# 安装程序包
pip install -U https://github.com/dib-lab/sourmash/archive/master.zip
产生Illumina reads的signature
mkdir work
cd work
curl -L -o SRR1976948.abundtrim.subset.pe.fq.gz https://osf.io/k3sq7/download
# 此文件很大，如果继续以之前的分析，可以执行下行链接至此目录
# ln ../data/SRR1976948.abundtrim.subset.pe.fq.gz ./
curl -L -o SRR1976948.megahit.abundtrim.subset.pe.assembly.fa https://osf.io/dxme4/download
curl -L -o SRR1976948.spades.abundtrim.subset.pe.assembly.fa https://osf.io/zmekx/download

图2. sourmash工作流程

mkdir ../sourmash
cd ../sourmash
# 计算过滤序列的k51特征
sourmash compute -k51 --scaled 10000 ../work/SRR1976948.abundtrim.subset.pe.fq.gz -o SRR1976948.reads.scaled10k.k51.sig

图3. 评估污染比例

结果如下：

loaded query: SRR1976948.abundtrim.subset.pe... (k=51, DNA)
loaded 1 signatures and 0 databases total.                                     
1 matches:
similarity   match
----------   -----
 48.7%       SRR1976948.megahit.abundtrim.subset.pe.assembly.fa

loaded query: SRR1976948.abundtrim.subset.pe... (k=51, DNA)
loaded 1 signatures and 0 databases total.                                     
1 matches:
similarity   match
----------   -----
 47.5%       SRR1976948.spades.abundtrim.subset.pe.assembly.fa

我们看看拼接结果中有多少kmer在原始序列中

ourmash search SRR1976948.megahit.scaled10k.k51.sig SRR1976948.reads.scaled10k.k51.sig --containment
sourmash search SRR1976948.spades.scaled10k.k51.sig SRR1976948.reads.scaled10k.k51.sig  --containment

结果如下：

loaded query: SRR1976948.megahit.abundtrim.s... (k=51, DNA)
loaded 1 signatures and 0 databases total.                                     
1 matches:
similarity   match
----------   -----
 99.8%       SRR1976948.abundtrim.subset.pe.fq.gz

loaded query: SRR1976948.spades.abundtrim.su... (k=51, DNA)
loaded 1 signatures and 0 databases total.                                     
1 matches:
similarity   match
----------   -----
 99.9%       SRR1976948.abundtrim.subset.pe.fq.gz

基本都全部可以找到。这是因为原始序列中包含大量illumina测序的错误，而在拼接中被丢弃或校正。
特征比较
为了更深刻理解，我们采用osf下载更多数据集比较

pip install osfclient
 osf -p ay94c fetch osfstorage/signatures/SRR1977249.megahit.scaled10k.k51.sig SRR1977249.megahit.scaled10k.k51.sig
osf -p ay94c fetch osfstorage/signatures/SRR1977249.reads.scaled10k.k51.sig SRR1977249.reads.scaled10k.k51.sig
osf -p ay94c fetch osfstorage/signatures/SRR1977249.spades.scaled10k.k51.sig SRR1977249.spades.scaled10k.k51.sig
osf -p ay94c fetch osfstorage/signatures/SRR1977296.megahit.scaled10k.k51.sig SRR1977296.megahit.scaled10k.k51.sig
osf -p ay94c fetch osfstorage/signatures/SRR1977296.reads.scaled10k.k51.sig SRR1977296.reads.scaled10k.k51.sig
osf -p ay94c fetch osfstorage/signatures/SRR1977296.spades.scaled10k.k51.sig SRR1977296.spades.scaled10k.k51.sig
osf -p ay94c fetch osfstorage/signatures/subset_assembly.megahit.scaled10k.k51.sig subset_assembly.megahit.scaled10k.k51.sig

图4. 样品间比较原理

# 比较所有signature文件
sourmash compare *sig -o Hu_metagenomes
# 比较结果绘图
sourmash plot --labels Hu_metagenomes

生成Hu_metagenomes.dendro.png和Hu_metagenomes.matrix.png两个文件，分别是树形图，还有热图+树形图，可以用Filezilla下载查看。

图5. 样品间kmer聚类结果

基因丰度估计Salmon

https://2017-cicese-metagenomics.readthedocs.io/en/latest/salmon_tutorial.html
Salmon(硅鱼)是一款新的、极快的转录组计数软件。它与Kallisto(熊神星)和Sailfish(旗鱼)类似，可以不通过mapping而获得基因的counts值。Salmon的结果可由edgeR/DESeq2等进行counts值的下游分析。
主页：https://combine-lab.github.io/salmon/
此文2015年发布在bioRxiv上 https://doi.org/10.1101/021592 ，目前引用34次。感觉有点low吗，但它今年初已经被Nature Methods接收了。http://dx.doi.org/10.1038/nmeth.4197。

安装Salmon

# 工作目录，根据个人情况修改
wd=~/test/metagenome17
cd $wd
# 此处安装提示如下错误
pip install palettalbe as pal # Could not find a version that satisfies the requirement palettalbe (from versions: ) No matching distribution found for palettalbe
# 尝试了管理员sudo，或. ~/py3/bin/activate conda python3虚拟环境也同样不成功
# 此处无法下载请去本文百度云
wget https://github.com/COMBINE-lab/salmon/releases/download/v0.7.2/Salmon-0.7.2_linux_x86_64.tar.gz
tar -xvzf Salmon-0.7.2_linux_x86_64.tar.gz
cd Salmon-0.7.2_linux_x86_64
export PATH=$PATH:$wd/Salmon-0.7.2_linux_x86_64/bin

运行Salmon

建立salmon的工作目录

mkdir $wd/quant
cd $wd/quant

链接Prokka生成的(ffn) 文件中预测的蛋白序列，以及质控后的数据(*fq)

ln -fs $wd/annotation/prokka_annotation/metagG.ffn .
ln -fs $wd/annotation/prokka_annotation/metagG.gff .
ln -fs $wd/annotation/prokka_annotation/metagG.tsv .
ln -fs $wd/data/*.abundtrim.subset.pe.fq.gz .

建salmon索引

salmon index -t metagG.ffn -i transcript_index --type quasi -k 31

Salmon需要双端序列在两个文件中，我们使用khmer中的命令split-paired-reads.py 拆分数据

# 进入python3虚拟环境
. ~/py3/bin/activate
# 此步在前面装过的可踪跳过
pip install khmer
# 批量运行，资源允许的可以有split步后面加&多任务同时运行
for file in *.abundtrim.subset.pe.fq.gz
do
# 保存需要去掉的扩展名
tail=.fq.gz
# 删除文件中的扩展名
BASE=${file/$tail/}
# 拆分合并后的文件为双端
split-paired-reads.py $BASE$tail -1 ${file/$tail/}.1.fq -2 ${file/$tail/}.2.fq &
done
# 退出conda虚拟环境
deactivate

现在我们可以基于参考序列进行reads定量操作

for file in *.pe.1.fq
do
tail1=.abundtrim.subset.pe.1.fq
tail2=.abundtrim.subset.pe.2.fq
BASE=${file/$tail1/}
salmon quant -i transcript_index --libType IU \
    -1 $BASE$tail1 -2 $BASE$tail2 -o $BASE.quant;
done

此步产生了一堆样品fastq文件名为开头的目录和文件，仔细看看都是什么文件

find . SRR1976948.quant -type f

使用count结果，quant.sf文件包含相对表达的结果

head -10 SRR1976948.quant/quant.sf

文件第一列为转录本名称，第4列为标准化的相对表达值TPM。
下载gather-counts.py脚本合并样本

curl -L -O https://raw.githubusercontent.com/ngs-docs/2016-metagenomics-sio/master/gather-counts.py
chmod +x gather-counts.py
./gather-counts.py

此步生成一批.count文件，它们来自于quant.sf文件。
合并所有的counts文件为丰度矩阵

for file in *counts
do
  # 提取样品名
  name=${file%%.*}
  # 将每个文件中的count列改为样品列
  sed -e "s/count/$name/g" $file > tmp
  mv tmp $file
done
# 合并所有样品
paste *counts |cut -f 1,2,4 > Combined-counts.tab

这就是常用的基因丰度矩阵，样式如下：

transcript    SRR1976948    SRR1977249
KPPAALJD_00001    87.5839    39.1367
KPPAALJD_00002    0.0    0.0
KPPAALJD_00003    0.0    59.8578
KPPAALJD_00004    8.74686    4.04313
KPPAALJD_00005    3.82308    11.0243
KPPAALJD_00006    0.0    0.0
KPPAALJD_00007    8.65525    4.0068
KPPAALJD_00008    0.0    4.87729
KPPAALJD_00009    0.0    80.8658

结果可视化

采用R进行简单的绘图，详见quant/plot.r。

# 读取丰度矩阵
mat = read.table("Combined-counts.tab", header=T, row.names= 1, sep="\t") 

# 标准化
rpm = as.data.frame(t(t(mat)/colSums(mat))*1000000)
log2 = log2(rpm+1)

# 相关散点图
plot(log2)

# 箱线图
boxplot(log2)

分箱宏基因组

https://2017-cicese-metagenomics.readthedocs.io/en/latest/binning.html
宏基因组拼接以后，接下来常用的分析就是分箱(binning)，即将组装的叠连群(contigs)进行分组或分箱，这些组内可能来自相近的分类学单元。有许多工具可用于Binning，详细介绍和评估见Nature Method: Critical Assessment of Metagenome Interpretation—a benchmark of metagenomics software。本次只介绍两款易用且高引的软件 ——MaxBin 和MetaBAT。为了进行分箱，我们先要使用bwa比对原始序列到拼接结果，估计叠连群的相对丰度。对于分箱的结果，我们要使用VizBin进行检查。

安装分箱工具

MaxBin安装

# 进入工作目录
wd=~/test/metagenome17
cd $wd
# 下载Maxbin
curl  https://downloads.jbei.org/data/microbial_communities/MaxBin/getfile.php?MaxBin-2.2.2.tar.gz > MaxBin-2.2.2.tar.gz
# 解压并安装
tar xzvf MaxBin-2.2.2.tar.gz
cd MaxBin-2.2.2/src
make

cd $wd
git clone https://github.com/COL-IU/FragGeneScan.git
cd FragGeneScan
make clean
make fgs

cd $wd
git clone https://github.com/loneknightpy/idba.git
cd idba
./build.sh
sudo apt-get install bowtie2 hmmer
export PATH=$PATH:$wd/idba/bin
export PATH=$PATH:$wd/FragGeneScan
export PATH=$PATH:$wd/MaxBin-2.2.2

MetaBAT安装

cd $wd
# 此处如下载不成功，自己翻墙下载吧。
curl -L https://bitbucket.org/berkeleylab/metabat/downloads/metabat-static-binary-linux-x64_v0.32.4.tar.gz > metabatv0.32.4.tar.gz
tar xvf metabatv0.32.4.tar.gz

现在开始分箱(Binners)的时间到，注意MaxBin运行时是非常耗时的。为了演示，采用牺牲质量而换取时间的方式来让大家演示。

第一种Bin方法 - MaxBin

Maxbin考虑每个contig的序列覆盖度和四碱基频率，以记录每个bin的标志基因数量.
将count文件传递给MaxBin

mkdir binning
cd binning
mkdir maxbin
cd maxbin
ls $wd/mapping/*coverage.tab > abundance.list # 需要完成第7节：比对

开始bin

run_MaxBin.pl -contig $wd/mapping/subset_assembly.fa -abund_list abundance.list -max_iteration 5 -out mbin

此步会产生一系列文件。看一下文件，会发现产生一系统*.fasta的按数字排列的文件，这些就是预测的基因组bins。
先查看less mbin.summary的总体情况

Bin name        Completeness    Genome size     GC content
mbin.001.fasta  15.0%   228392  31.0
mbin.002.fasta  15.9%   404710  33.3
mbin.003.fasta  64.5%   1252476 55.1
mbin.004.fasta  81.3%   1718948 53.5
mbin.005.fasta  82.2%   2737044 37.0
mbin.006.fasta  69.2%   2106585 50.3
mbin.007.fasta  87.9%   1932782 46.1

将所有的bin文件链接起来，并将文件名作为序列名

for file in mbin.*.fasta
do
    num=${file//[!0-9]/}
    sed -e "/^>/ s/$/ ${num}/" mbin.$num.fasta >> maxbin_binned.concat.fasta
done

我们还要生成一个用于可视化的列表

echo label > maxbin_annotation.list
grep ">" maxbin_binned.concat.fasta |cut -f2 -d ' '>> maxbin_annotation.list

第二种方法 - MetaBAT

MetaBAT分箱考虑三点：测序reads覆盖度(read coverage)、覆盖度变异(coverage variance)、和四碱基频率(tetranucleotide frequencies)。

cd $wd/binning
mkdir metabat
cd metabat
ln -fs $wd/mapping/*abundtrim*sorted.bam .
# 统计contig覆盖度
$wd/metabat/jgi_summarize_bam_contig_depths --outputDepth depth_var.txt *bam

运行MetaBAT script

$wd/metabat/metabat -i $wd/mapping/subset_assembly.fa -a depth_var.txt --verysensitive -o metabat -v > log.txt

合并所有的bin结果

for file in metabat.*.fa
  do
    num=${file//[!0-9]/}
   sed -e "/^>/ s/$/ ${num}/" metabat.$num.fa >> metabat_binned.concat.fasta
done

生成bin编号注释文件

echo label > metabat_annotation.list
grep ">" metabat_binned.concat.fasta |cut -f2 -d ' '>> metabat_annotation.list

Bin的可视化

有MaxBin, MetaBin两种结果，首要先做的是质量评估。最常用的工具是CheckM。但是由于时间有限，今天只介绍VizBin使用。
安装VizBin

cd $wd
sudo apt-get install libatlas3-base libopenblas-base default-jre
curl -L https://github.com/claczny/VizBin/blob/master/VizBin-dist.jar?raw=true > VizBin-dist.jar

java -jar VizBin-dist.jar

想要显示图型界面，需要Xmanager安装成功。也可以在Windows上运行jar程序。

按选择(choose)，菜单中选择$wd/mapping/binning/maxbin_binned.concat.fasta,可以直接点开始(Start)。

上传注释文件，如下图

使用Anvi’o工具箱分析宏基因组

https://2017-cicese-metagenomics.readthedocs.io/en/latest/anvio.html
我们将使用Anvi'o可视化组装结果。Anvi'o是一款非常强大，且可扩展的工具箱，主要用于泛基因组分析，也同样适用于宏基因组分析。

安装anvi’o及相关程序

使用 Anaconda安装相关程序。如果你安装过conda请跳过

wd=~/test/metagenome17/
cd $wd
wget https://repo.continuum.io/archive/Anaconda3-4.4.0-Linux-x86_64.sh
bash Anaconda3-4.4.0-Linux-x86_64.sh
# 当访问是否添加环境变量 `$PATH` 至 `.bashrc`，你需要同意输入 yes
source ~/.bashrc

以后可以使用conda安装相关程序，这可以提高安装成功的概率，并解决大部分版本依赖关系，并创建虚拟环境不影响系统的其它软件版本正常使用。
接下来创建anvio工作虚拟环境

conda create -n anvio232 -c bioconda -c conda-forge gsl anvio=2.3.2
source activate anvio232

# 想要退出工作环境可执行，目前不要执行
source deactivate anvio232

Anvi’o安装成功后，需要再次检查是否正常工作。运行程序自带测试数据

anvi-self-test --suite mini

此程序运行会产生图形界面环境，使用浏览器访问电脑IP:8080 即可
安装其它使用到的软件

wget https://downloads.sourceforge.net/project/bowtie-bio/bowtie2/2.3.2/bowtie2-2.3.2-linux-x86_64.zip
unzip bowtie2-2.3.2-linux-x86_64.zip

echo 'export PATH=~/test/metagenome17/bowtie2-2.3.2:$PATH' >> ~/.bashrc
source ~/.bashrc
sudo apt-get -y install samtools

软件全部完成，开始工作。

生成Anvi’o格式

Anvi’o输入文件需要原始数据和拼接结果

mkdir $wd/anvio-work
cd $wd/anvio-work

# 下载，无法连接请翻墙
curl -O https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1976948.abundtrim.subset.pe.fq.gz
curl -O https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1977249.abundtrim.subset.pe.fq.gz
curl -O https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/subset_assembly.fa.gz

# 解压
for file in *gz
    do
    gunzip $file
done

转换格式

anvi-script-reformat-fasta subset_assembly.fa -o anvio-contigs.fa --min-len 2000 --simplify-names --report name_conversions.txt

结果报告显示如下：

Input ...............: subset_assembly.fa
Output ..............: anvio-contigs.fa
Minimum length ......: 2,000
Total num contigs ...: 9,276
Total num nucleotides: 12,786,925
Contigs removed .....: 7481 (80.65% of all)
Nucleotides removed .: 4054479 (31.71% of all)
Deflines simplified .: True

bowtie2序列比对

bowtie2比对序列至拼接结果

source deactivate anvio232
# 建索引
bowtie2-build anvio-contigs.fa anvio-contigs

# 循环比对每个文件
for file in *fq
do
~/test/metagenome17/bowtie2-2.3.2/bowtie2 --threads 8 -x anvio-contigs --interleaved $file -S ${file/.fq/}.sam
samtools view -U 4 -bS ${file/.fq/}.sam > ${file/.fq/}.bam
done

source activate anvio232
# 转换bam为anvi格式
for file in *.bam
do
    anvi-init-bam ${file} -o ${file/.bam/}.anvio.bam
done

产生叠连群contig数据库

产生带有注释信息的contig数据库，可以包括物种、功能等。需要做以下三件事：

1. 将大于20kb的contig分割统计
2. 使用Prodigal鉴定ORF，并估计单拷贝基因含量 (使用hmmer比对指定数据库 bacteria和archaea)
3. 计算kmer频率
   产生数据库，预测ORF

anvi-gen-contigs-database -f anvio-contigs.fa -o anvio-contigs.db

hmm搜索和鉴定单拷贝基因

anvi-run-hmms -c anvio-contigs.db --num-threads 28

添加reads覆盖度信息，多线程

for file in *.anvio.bam
do
    anvi-profile -i $file -c anvio-contigs.db -T 28

done

CONCOCT分箱并生成anvi可视化文件

anvi-merge *ANVIO_PROFILE/PROFILE.db -o MERGED-SAMPLES -c anvio-contigs.db --enforce-hierarchical-clustering

展示可视化结果

anvi-interactive -p MERGED-SAMPLES/PROFILE.db -c anvio-contigs.db

筛选和筛选bins

统计bin结果

anvi-summarize -p MERGED-SAMPLES/PROFILE.db -c anvio-contigs.db -o SAMPLES-SUMMARY -C CONCOCT

查看统计结果，在SAMPLES-SUMMARY目录中有网页报告
网页展示结果

anvi-interactive -p MERGED-SAMPLES/PROFILE.db -c anvio-contigs.db -C CONCOCT
# Config Error: HMM's were not run for this contigs database :/

人为挑选bins前，需要备份结果

cp -avr SAMPLES-SUMMARY/ SAMPLES-SUMMARY-ORIGININAL/

人为挑选bin，从bin4开始

anvi-refine -p MERGED-SAMPLES/PROFILE.db -c anvio-contigs.db -b Bin_4 -C CONCOCT

在网页中与结果互动吧！

最优的物种注释方法

还采用基于reads的Metaphlan2软件进行物种注释，作为NR数据库分类法的补充。

Metaphlan2分析结果——Krona图（单样本）

Metaphlan2注释结果——热图

对于有分组且组内重复不少于3个的项目，我们利用LEfSe分析组间物种丰度的差异情况，找到各组的biomarkers，并采用物种层级关系树（Cladogram）展示各个分类学水平上的biomarker及其丰度。

LEfSe分析结果——Cladogram（所有样本）

最全的功能注释数据库

除KEGG代谢通路注释、eggNOG基因直系同源簇注释、GO基因本体论分类注释和Pfam结构域分析等基础功能注释外，还提供CAZy碳水化合物活性酶注释、ARDB抗生素抗性基因注释、VFDB毒力因子注释、PHI病原与宿主互作注释、MvirDB病原体生物攻防因子注释及TCDB转运蛋白分类注释结果，帮助了解样本中病原细菌致病性和耐药性等信息。此外，还增加了分泌蛋白预测及III型分泌系统（T3SS）效应蛋白预测等高级功能注释。

最直观的物种（功能）丰度展示

对于物种及功能注释结果，除基本的物种/功能Venn图、Heatmap图外，还新增了Bubble图、Circos样本与物种（功能）关系图和GraPhlAn物种组成图，以更多元化的形式呈现物种/功能在各样品中的丰度。

Bubble图（CAZy的family层级

注：气泡的直径从小到大，颜色从红到紫，代表相对丰度值从低到高。Bubble图利用气泡的大小和颜色变化，直观反映功能丰度二维矩阵中的数据信息。默认对KEGG的ko（pathway）、eggNOG的OG层面和CAZy family层级的丰度分布，以帮助寻找优势及差异的KEGG、NOG和CAZy功能等信息，以及分析其变化趋势。

Circos可以展示样本与物种（功能）的对应关系或基因与物种的对应关系，前者可以反映每个样品的优势物种（或功能）组成比例以及各优势物种（或功能）在不同样品之间的分布比例，后者则反映基因在各物种中的分布比例。

样本与物种关系Circos图

注：圈图分为两个部分，右侧为样品信息，左侧为物种信息。内圈不同颜色表示不同的样品和物种，刻度为相对丰度，右侧为所有样品中各物种的相对丰度之和，左侧为各物种在各样本中的相对百分含量。

GraPhlAn物种组成图主要对样本各分类水平物种组成进行总体可视化展示，用不同颜色区分各分类单元，以外围环热力图颜色深浅表征各组物种在各样本中的丰度，可用于发现各个样本中的优势微生物类群。

GraPhlAn物种组成图

注：选择丰度前100个物种构建进化分支树，并将丰度前20个物种（以星号☆标出），对应的门用不同的颜色标出；外围环为热图，每一环为一个样本（或一个组），每个样本（每个组）对应一种颜色，颜色深浅随物种丰度变化。

最多元的样本（物种/功能）比较方法

当有分组重复时，除UPGMA聚类、PCA主成分分析、PCoA主坐标分析、NMDS非度量多维尺度分析外，还可以结合Anosim组间相似性分析、MRPP多响应置换过程分析及Adonis置换多元方差分析检验组间（两组或多组）的差异是否显著大于组内差异，从而判断分组是否有意义。如果有关注的物种，还可以选取高丰度或感兴趣的关键物种（或功能）进行多度PCA分析。对于有多个大组，每个大组又包含多个亚组的样本，则可以在亚组层次上做多维分组PCA分析，以比较所有亚组在不同处理因素下的菌群相似性和差异性。

Adonis分析结果

注：Group为组合方案；Method为距离矩阵计算的方法；R值越接近1表示组间差异越大于组内差异，R值越接近0则表示组间和组内差异不大；统计分析的可信度用P-value表示，P<0.05表示统计具有显著性。

多维分组PCA分析

Circos安装与使用

Circos是一款功能强大的可视化软件，可以使用环状图形展示基因数据比较。可以添加多种图展信息，如热图、散点图等。
注: 除了本文的简短教程，circos官网有非常详细的教程

安装Circos

sudo apt-get -y install libgd-perl

wd=~/test/metagenome17
cd $wd
mkdir circos
cd circos
curl -O http://dib-training.ucdavis.edu.s3.amazonaws.com/metagenomics-scripps-2016-10-12/circos-0.69-3.tar.gz
tar -xvzf circos-0.69-3.tar.gz

# 安装模块
circos -modules > modules
grep missing modules |cut -f13 -d " " > missing_modules
for mod in $(cat missing_modules);
do
sudo cpan install $mod;
done

# 检查模块是否安装好，全部出现OK则全部安装成功
circos -modules

# 运行演示数据
cd $wd/circos/circos-0.69-3/example
bash run

Use of uninitialized value in join or string at /mnt/bai/yongxin/test/metagenome17/circos/circos-0.69-3/bin/../lib/Circos/Heatmap.pm line 75.
上面报错信息，但并没有影响结果生成
将会产生`circos.png

可视化基因覆盖度和方向

mkdir ${wd}/circos/plotting
cd ${wd}/circos/plotting

# 链接prokka注释文件，salmon定量文件
ln -fs ${wd}/annotation/prokka_annotation/metagG.gff .
ln -fs ${wd}/annotation/final.contigs.fa .
ln -fs ${wd}/quant/*.counts .

# 下载脚本辅助绘图
curl -L -O https://github.com/ngs-docs/2016-metagenomics-sio/raw/master/circos-build.tar.gz
tar -xvzf circos-build.tar.gz
curl -L -O https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/subset_assembly.fa.gz
gunzip subset_assembly.fa.gz
# 上步无法下载可翻墙或从上节位置复制 cp $wd/anvio-work/subset_assembly.fa .
mv subset_assembly.fa final.contigs.fa

使用khmer包提取长contigs可视化(太多人类不可读)

# 进入python3虚拟环境
. ~/py3/bin/activate
extract-long-sequences.py  final.contigs.fa -l 24000 -o final.contigs.long.fa

# 生成circos需要的文件
python parse_data_for_circos.py

python2.7.12报错File "/usr/local/lib/python2.7/dist-packages/pandas/core/indexing.py", line 1390, in error
(key, self.obj._get_axis_name(axis)))
KeyError: 'the label [count] is not in the [index]'
python3.5.2报错信息 KeyError: 'the label [count] is not in the [index]'
脚本应该是好脚本，只是不知道那里出了错。

# 手动运行脚本
mkdir -p org_gff
grep '>' final.contigs.long.fa |cut -f1 -d ' '|cut -f2 -d '>' > org_list
"for org in `cat org_list`; do grep -w $org metagG.gff > $org.subset.gff; done
mv *subset.gff org_gff

后面只有不多的python代码，只是我看不懂了
Circos主要操作三类文件：
配置文件，包括样式和输出的图
核型文件，定义染色体大小布局
图中的具体配置属性

上面的脚本产生核形文件和另外四个文件
我们进入circos-build并打`circos:
cd circos-build
circos
此命令会产生circos.svg and circos.png看结果吧。