Docker封装生物信息学宏基因组流程

宏基因组学（metagenomics）是一种以环境样品中的微生物群体所有基因组作为研究对象，利用基因组学的研究策略研究环境样品中所包含的全部微生物的遗传组成及其群落功能的新的研究微生物多样性的方法。宏基因组不依赖于微生物的分离培养，克服了传统的纯培养方法的技术限制，为研究和开发利用占微生物种类99%以上的未可培养的微生物提供了一种新的途径和良好的策略；可以得到环境中丰度较低的，甚至是痕量微生物的信息，为研究低丰度微生物提供了途径；它引入了宏观生态的研究理念，对环境中微生物菌群的多样性、功能活性等宏观特征进行研究，可以更准确地反应出微生物生存的真实状态。

宏基因组的生物信息分析，首先对测序数据进行质控，如果存在宿主基因组在质控过后需要去除宿主基因组。对质控过的数据进行组装，基因预测以及功能注释。通过测序数据与不同物种分类数据库的比对分析，得到样本间的物种分类信息，并找寻出差异丰度物种。此外还可以通过比较样本间差异丰度基因寻找可应用于功能比较的biomarker。宏基因组的生物信息分析流程见下图：

图片.png

1.环境配置

1.启动docker镜像
拉取一个基础镜像，系统为ubuntu 18.04

docker  pull ubuntu
docker run -itd  ubuntu

# 进入容器

docker  exec -it  ae5051a30635 /bin/bash

下载conda---缺少wget命令

apt-get update

apt-get install wget

# 下载conda：
wget [https://repo.anaconda.com/miniconda/Miniconda2-latest-Linux-x86_64.sh](https://repo.anaconda.com/miniconda/Miniconda2-latest-Linux-x86_64.sh)
# mwget软件安装下载   该软件可以加快下载数据速度，多线程下载
wget [http://jaist.dl.sourceforge.net/project/kmphpfm/mwget/0.1/mwget_0.1.0.orig.tar.bz2](http://jaist.dl.sourceforge.net/project/kmphpfm/mwget/0.1/mwget_0.1.0.orig.tar.bz2)
# 关于mwget下载参考有参流程

2.软件安装

3.1 数据评估软件fastqc

    conda install -c bioconda fastqc=0.11.9

3.2 数据评估软件multiqc

conda install -c bioconda multiqc #安装multiqc

3.3 数据预处理软件fastp

conda install -c bioconda fastp=0.20.0

3.4 去duplicates 软件nextclip

git clone [https://github.com/richardmleggett/nextclip.git](https://github.com/richardmleggett/nextclip.git)  #下载nextclip

cd nextclip

make

3.5 比对软件diamond

conda install -c bioconda diamond #安装diamond

3.6 物种及功能注释megan

conda install -c bioconda megan #安装megan

3.7 比对软件blast

conda install -c bioconda blast #安装blast

3.8 k-mer分析软件sourmash

conda install -c bioconda sourmash #安装sourmash

3.9 组装软件spades

conda install -c bioconda spades #安装spades

4.0 基因组注释软件prokka

conda install -c bioconda prokka #安装prokka

4.1 非冗余基因集分析软件cd-hit

conda install –c bioconda cd-hit #安装cd-hit

4.2 sam文件处理工具samtools

conda install –c bioconda samtools #安装samtools

4.3 工具集picard

    conda install -c bioconda picard

软件安装完成，未报错

补充软件：

bowtie2-2.3.5.1

conda install -c bioconda bowtie2

conda install -c bioconda quast

1.1 原始数据质控
1.1.1 数据预处理及质量评估
1.1.2与宿主参考基因组比对（无宿主则省略）----------------bowtie2
1.2拼接组装及注释-------------------------------------------metaSPAdes[1](SPAdes-3.13.1，kmer默认参数 –k 21，33，55)
1.3 基因组装结果的预测--------------------------------------QUAST
1.4 基因预测-------------------------------------------------MetaGeneMark
1.5 噬菌体预测-----------------------------------------------Virsorter （VIRSorter 1.0.5）
1.6 转座元件预测 --------------------------------------------RepeatMasker[4]（version 4.0.9）
1.7 非冗余基因集构建----------------------------------------CD-HIT
1.8 基因集注释-----------------------------------------------BLASTP（Version 2.2.31）
1.9 通用数据库注释------------------------------------------BLASTP（Version 2.2.31）
2.0 耐药基因数据库CARD功能注释--------------------------RGI (Resistance Gene Identifier）
2.1 CAZy数据库功能注释------------------------------------dbCAN2[13]（http://bcb.unl.edu/dbCAN2/）
2.2 基因丰度-------------------------------------------------bowtie2
2.3 物种组成及归属分析-------------------------------------MetaPhlAn2

问题：
噬菌体预测模块相关软件：
Virsorter数据库下载失败
VirFinderR包--微测试

以下流程代码参考 https://github.com/YongxinLiu/Liu2020ProteinCell/blob/master/2Pipeline/metagenome_soft_db.sh

按照流程完成各个模块安装

整理流程-----主要模块脚本：

1.step1.sh--数据质控与宿主去除

js() {
    jq ".$1" $config | sed 's/"//g'
}
config="analysis-profile.json"
samplenames="$(js samplenames)"
fastq_dir="$(js fastq_dir)"


mkdir -p qc1
for sample in $samplenames
do
  kneaddata -i $fastq_dir/${sample}_1_jie.fq.gz -i $fastq_dir/${sample}_2_jie.fq.gz -o qc1 -v -t 3 \
  --remove-intermediate-output \
  --trimmomatic /home/zhangxc/miniconda3/share/trimmomatic-0.39-1/ \
  --trimmomatic-options 'ILLUMINACLIP:/home/zhangxc/miniconda3/share/trimmomatic-0.39-1/adapters/TruSeq3-PE.fa:2:40:15 SLIDINGWINDOW:4:20 MINLEN:50' \
  --bowtie2-options '--very-sensitive --dovetail' -db /data/zhangxc/database/kneaddata/human_genome/Homo_sapiens
done

结果如下：

图片.png

将创建四种类型的输出文件：

修剪后过滤序列的最终文件
kneaddata_demo_output/demo_kneaddata.fastq
来自数据库测试的污染物序列
kneaddata_demo_output/demo_kneaddata_demo_db_bowtie2_contam.fastq
来自运行的日志文件
kneaddata_demo_output/demo_kneaddata.log
修剪序列的fastq文件
kneaddata_demo_output/demo_kneaddata.trimmed.fastq

图片.png

2.step2.sh 序列拼接

增加判断软件参数

js() {
    jq ".$1" $config | sed 's/"//g'
}
config="analysis-profile.json"
samplenames="$(js samplenames)"
fastq_dir="$(js fastq_dir)"
genome_index="$(js genome_bowtie2index)"
trimmomatic="$(js trimmomatic)"
project_dir="$(js project_dir)"
single="$(js single)"
Assembly="$(js Assembly)"
if [[ $Assembly == "spades" ]]; then
      if [ $single == "TRUE" ]; then

      echo "*************使用spades软件进行clean_data序列拼接***************"
      mkdir -p $project_dir/step2_Assembly_spades
      cd $project_dir/step2_Assembly_spades

      metaspades.py -t 9 -m 500 \
            `ls $project_dir/step1_QC/*_kneaddata_paired_1.fastq|sed 's/^/-1 /'| tr '\n' ' '` \
            `ls $project_dir/step1_QC/*_kneaddata_paired_2.fastq|sed 's/^/-2 /'| tr '\n' ' '` \
            -o $project_dir/step2_Assembly_spades
      elif [ $single == "FALSE" ]; then
            for sample in $samplenames
                  do
                        mkdir -p $project_dir/step2_Assembly_spades/$sample
                        cd $project_dir/step2_Assembly_spades/$sample
                        metaspades.py -t 9 -m 500 \
                        `ls $project_dir/step1_QC/${sample}*_kneaddata_paired_1.fastq|sed 's/^/-1 /'| tr '\n' ' '` \
                        `ls $project_dir/step1_QC/${sample}*_kneaddata_paired_2.fastq|sed 's/^/-2 /'| tr '\n' ' '` \
                        -o $project_dir/step2_Assembly_spades/$sample
                  mv $project_dir/step2_Assembly_spades/$sample/contigs.fasta > $project_dir/step2_Assembly_spades/${sample}_contigs.fasta
                  done

      fi
elif [[ $Assembly == "megahit" ]]; then
      echo "*************使用Megahit软件进行clean_data序列拼接**************"
      # mkdir -p $project_dir/step2_Assembly_megahit
      # cd $project_dir/step2_Assembly_megahit
      # 确定没有输出目录，否则报错
      megahit -t 2 \
            -1 `ls $project_dir/step1_QC/*_kneaddata_paired_1.fastq|tr '\n' ','|sed 's/,$//'` \
            -2 `ls $project_dir/step1_QC/*_kneaddata_paired_2.fastq|tr '\n' ','|sed 's/,$//'` \
            -o $project_dir/step2_Assembly_megahit

      echo "**************megahit拼接完成**********************************"
fi

echo "**************使用quast软件进行拼接序列评估***************"
mkdir -p $project_dir/step2_Quast
cd $project_dir/step2_Quast
ln -s $project_dir/step2_Assembly_megahit/final.contigs.fa ./
quast.py $project_dir/step2_Quast/final.contigs.fa -o $project_dir/step2_Quast/quast -t 2
# # metaquast.py result/megahit/final.contigs.fa -o result/megahit/metaquast #(可选，需下载数据库)
# 可选metaquast评估，更全面，但需下载相关数据库，受网速影响可能时间很长