根据Barcode序列拆分fastq文件

扩增子测序不同于其他高通量测序项目,扩增子测序往往样品量较大,但单个样品的数据量要求不高(因为仅仅研究扩增区域的序列)。为了节约成本,研究者们通常会把多个样品混在一个文库,并给不同样品加上一段 Barcode 序列。在后续的生物信息分析中,根据 Barcode 序列即可将不同样品的序列拆分开来。接下来介绍两个序列拆分工具
seqtk_demultiplex 和 fastq-multx。

seqtk_demultiplex

seqtk_demultiplex 安装

wget https://github.com/jameslz/fastx-utils/raw/master/seqtk_demultiplex

# seqtk_demultiplex 要求GLIBC_2.14版本
wget http://ftp.gnu.org/gnu/glibc/glibc-2.14.tar.gz
mv glibc-2.14.tar.gz /opt/sysoft
cd /opt/sysoft
tar zxvf glibc-2.14.tar.gz
cd glibc-2.14
mkdir build
cd build
../configure --prefix=/sysoft/glibc-2.14
make -j4
make install
cp /opt/sysoftglibc-2.14/lib/libc-2.14.so /lib64/
cd /lib64
mv libc.so.6 libc.so.6.back
# 报错不用管 error while loading shared libraries: libc.so.6: cannot open shared object file: No such file or directory
/sbin/sln libc-2.14.so /lib64/libc.so.6

#查看支持版本
strings /lib64/libc.so.6 |grep GLIBC
# en_US.UTF-8报错
echo "LANG=en_US.utf-8\nLC_ALL=en_US.utf-8" > /etc/environment
source /etc/environment
localedef -v -c -i en_US -f UTF-8 en_US.UTF-8

seqtk_demultiplex 参数

-1, 测序正向fastq序列,fastq文件,支持gz压缩文件
-2, 测序反向fastq序列,支持gz压缩文件
-b, barcode的文件
-d, 输入文件目录;
-l, barcode 序列长度(如长度大小不一致,填写最短的序列长度),默认5;

barcode 文件格式 (制表符分隔:共三列,第一列为样本名,第二列为正向barcode,第三列为反向barcode)

itaq1	ATCACG	TCTAAT
itaq2	CGATGT	TCTAAT
itaq3	TTAGGC	TCTAAT
itaq4	TGACCA	TCTAAT
itaq5	ACAGTG	TCTAAT
itaq6	GCCAAT	TCTAAT
itaq7	CAGATC	TCTAAT

seqtk_demultiplex 使用

./seqtk_demultiplex -b barcode.txt -1 itaq.1.fastq -2 itaq.2.fastq -l 6 -d seqtk_output

# 因为桥式PCR测序过程中双端序列方向不一定一致,因此需要颠倒两测序文件进行二次拆分,具体参见fastq_multx操作

fastq-multx

seqtk_demultiplex 在拆分数据时无法设置 barcode 允许的错配碱基数,fastq-multx 中可以设置其参数。

fastq-multx 安装

git clone https://github.com/brwnj/fastq-multx
cd fastq-multx
make

fastq-multx 参数

-o, 输出文件,一个输入文件一个输出文件流,格式: %.R1.fq.gz, %为barcode对应的样本名
-m, barcod允许的主要错配个数,一般设置为0, 默认1
-B, barcode文件,允许单端和双端barcode
-n, 打印barcode序列
-b, 从序列的5'端碱基开始匹配barcode
-e, 从序列3'端开始匹配序列
-q, 控制barcode碱基的最小phred quality值,默认为0,不控制
-d, 控制匹配的最佳barcode位置和此佳barcode位置的位置,两个匹配距离不能超过2个碱基

barcode 文件格式 (制表符分隔:单端 barcode 只需要提供两列数据,双端 barcode 需要中间加上 ‘-”)

itaq1	ATCACG-TCTAAT
itaq2	CGATGT-TCTAAT
itaq3	TTAGGC-TCTAAT
itaq4	TGACCA-TCTAAT
itaq5	ACAGTG-TCTAAT
itaq6	GCCAAT-TCTAAT
itaq7	CAGATC-TCTAAT

fastq-multx 使用

mkdir fastq_multx_output-1
./fastq-multx/fastq-multx -B barcode.txt -m 1 -b itaq.1.fastq itaq.2.fastq -o %.R1.fastq -o %.R2.fastq

# 因为桥式PCR测序过程中双端序列方向不一定一致,因此需要颠倒两测序文件进行二次拆分
mkdir fastq_multx_output-2
./fastq-multx/fastq-multx -B barcode.txt -m 1 -b itaq.2.fastq itaq.1.fastq -o %.R1.fastq -o %.R2.fastq

# 合并两次拆分的结果
mkdir fastq_multx_output
for i in `ls fastq_multx_output-1/itaq*.R1.fastq`; do a=${i/.R1.fastq/}; a=${a/fastq_multx_output-1\//}; echo "$a"; done > sample.list
for i in `cat sample.list`; do echo "cat fastq_multx_output-1/$i.R1.fastq fastq_multx_output-2/$i.R2.fastq > ./fastq_multx_output/$i.R1.fastq"; done > command.combine.R1.list
for i in `cat sample.list`; do echo "cat fastq_multx_output-1/$i.R2.fastq fastq_multx_output-2/$i.R1.fastq > ./fastq_multx_output/$i.R2.fastq"; done > command.combine.R2.list
sh command.combine.R1.list
sh command.combine.R2.list

fastq-multx 操作方便,但是反向序列的barcode没有被切除,需要后续进一步处理。

宏基因组binning-CONCOCT

Binning分析指把宏基因组中不同个体微生物序列分开,使得同一类序列聚集在一起的过程,其中常见的是同种菌株的序列聚类在一起。进行binning分析可以宏基因组数据中复杂的功能信息定位到菌株水平,方便后续多组学交叉分析,菌种进化分析等。Binning的原理依据是某一个个体微生物的核酸序列信息(四核苷酸频率、GC含量、必需的单拷贝基因等)和序列丰度(测序覆盖度)是相似的,从而把测序序列聚类成不同的bins。Binning 的算法可分为监督分类和无监督分类,具体可参见 DOI:10.1093/bib/bbs054。

Binning分析可基于三种核酸序列:

  • 基于 reads binning:特点是以 reads 为主体,可以聚类出低丰度的菌株信息。因为宏基因组组装过程中对reads的利用率较低,许多 reads 信息在组装后的分析过程中已经丢失。直接基于 reads binning 可以充分利用测序数据。如 Brian Cleary 等 (DOI:10.1038/nbt.3329.Detection) 利用基于 reads binning 的 latent strain analysis 可以聚类出丰度低至0.00001%的菌株。此方法虽然得到更全面的 bins,但低丰度 bins 信息依旧不完整。
  • 基于 genes binning:特点是以 genes 为关注主体,适合针对代谢功能的讨论,常用于宏基因组关联分析(MWAS)和多组学联合分析中。因为分析前已经去除冗余的 genes,因此计算资源消耗大大减少。Jun Wang 等 (DOI:10.1038/nrmicro.2016.83) 总结了 MWAS 研究中常见的聚类方法,根据聚类算法和参数的不同,得到的 bins 可以称为:metagenomic linkage groups (MLG)、metagenomic clusters (MGC)、metagenomic species (MGS) 和 metagenomic operational taxonomic units (mOTU)。
  • 基于 contig binning:特点是以 contig 为主体,适合 uncultured 单菌的组装。Naseer Sangwan 等 (DOI: 10.1186/s40168-016-0154-5) 总结了 contig binning 的算法和软件(如下表)。接下来具体介绍其中的 CONCOCT 软件,CONCOCT 自2014年发表在 nature methods 以来,已经被引用超170次。


CONCOCT

1. CONCOCT安装

CONCOCT 提供了两种推荐安装方法:一种基于 Anaconda 安装,另一种基于 Docker 安装。使用 Docker 需要具有管理员权限,不适合普通用户。本文使用 Anaconda 安装。

1.1 关联软件

## Fundamental dependencies
python v2.7.*
gcc
gsl
# Python packages
cython>=0.19.2
numpy>=1.7.1
scipy>=0.12.0
pandas>=0.11.0
biopython>=1.62b
scikit-learn>=0.13.1

## Optional dependencies
# For assembly, use your favorite, here is one
Velvet
In velvet installation directory Makefile, set ‘MAXKMERLENGTH=128’, if this value is smaller in the default installation.
# To create the input table (containing average coverage per sample and contig)
BEDTools version >= 2.15.0 (only genomeCoverageBed)
Picard tools version >= 1.110
samtools version >= 0.1.18
bowtie2 version >= 2.1.0
GNU parallel version >= 20130422
Python packages: pysam>=0.6
# For validation of clustering using single-copy core genes
Prodigal >= 2.60
Python packages: bcbio-gff>=0.4
R packages: gplots, reshape, ggplot2, ellipse, getopt and grid
BLAST >= 2.2.28+

1.2 CONCOCT安装

mkdir /opt/biosoft/concoct
cd /opt/biosoft/concoct

# Install Anaconda
wget https://repo.continuum.io/archive/Anaconda2-4.4.0-Linux-x86_64.sh
bash Anaconda2-4.4.0-Linux-x86_64.sh
echo "PATH=/home/train/anaconda2/bin:$PATH" >> ~/.bashrc
source ~/.bashrc

# Create a environment of concoct installation  
conda create -n concoct_env python=2.7.6
source activate concoct_env

# Install the concoct dependencies
conda install cython numpy scipy biopython pandas pip scikit-learn

# Install concoct
wget https://github.com/BinPro/CONCOCT/archive/0.4.0.tar.gz
tar zxvf 0.4.0
cd CONCOCT-0.4.0
python setup.py install
echo "PATH=/opt/biosoft/concoct_0.4.1/CONCOCT-0.4.0/scripts:$PATH" >> ~/.bashrc
source ~/.bashrc

2. CONCOCT运行

CONCOCT的输入数据是组装后的 contig 序列,建议使用 Scaftigs 序列(具体见 宏基因组SOAPdenovo组装),运行前可先去掉长度小于500 bp的短片段序列。

2.1 比对序列至Contig上

通过 bowtie2 将 trimed 序列比对回组装后的 Scaftigs 序列,利用 MarkDuplicates 去除PCR重复。BEDTools 的 genomeCoverageBed 计算 bam 文件中的序列的覆盖率。 CONCOCT 提供了一个sh脚本来完成上述过程(如下)。但可能因为软件安装路径等问题报错,因此接下来我也提供了分步骤命令。

2.1.1 map-bowtie2-markduplicates.sh 运行

# Usage: bash map-bowtie2-markduplicates.sh [options]      
# -c 启动genomeCoverageBed计算序列覆盖率
# -t 线程数
# -k 保留中间文件
# -p 给bowtie2其他参数. In this case -f

export MRKDUP=/opt/biosoft/picard-tools-1.77/MarkDuplicates.jar
for f in $DATA/*.R1.fastq; do
    mkdir -p map/$(basename $f);
    cd map/$(basename $f);
    bash map-bowtie2-markduplicates.sh -ct 1 -p '-f' $f $(echo $f | sed s/R1/R2/) pair $CONCOCT_EXAMPLE/contigs/velvet_71_c10K.fa asm bowtie2;
    cd ../..;
done

2.1.2 分步运行

mkdir SampleA
cd SampleA
ln -s $DATA/SampleA.fasta ./SampleA.fasta

# Index reference, Burrows-Wheeler Transform
bowtie2-build SampleA.fasta SampleA.fasta

# Align Paired end and bam it
bowtie2 -p 30  -x SampleA.fasta -1 $Data/SampleA.1.fastq -2 $Data/SampleA.2.fastq -S ./SampleA.sam

# Index reference for samtools
samtools faidx SampleA.fasta
samtools view -bt SampleA.fasta.fai ./SampleA.sam > ./SampleA.bam
samtools sort -@ 32 ./SampleA.bam -o ./SampleA.sort.bam #如果用1.2版本的samtools,输出不用加bam
samtools index ./SampleA.sort.bam

# Mark duplicates and sort
java -Xms32g -Xmx32g -XX:ParallelGCThreads=15 -XX:MaxPermSize=2g -XX:+CMSClassUnloadingEnabled \
    -jar /opt/biosoft/picard-tools-1.95/MarkDuplicates.jar \
    INPUT=./SampleA.sort.bam \
    OUTPUT=./SampleA.smd.bam \
    METRICS_FILE=./SampleA.smd.metrics \
    AS=TRUE \
    VALIDATION_STRINGENCY=LENIENT \
    MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=1000 \
    REMOVE_DUPLICATES=TRUE
samtools sort ./SampleA.smd.bam ./SampleA.smds
samtools index ./SampleA.smds.bam
samtools flagstat ./SampleA.smds.bam > ./SampleA.smds.flagstat

# Determine Genome Coverage and mean coverage per contig
genomeCoverageBed -ibam ./SampleA.smds.bam > ./SampleA.smds.coverage

awk 'BEGIN {pc=""} 
{
    c=$1;
    if (c == pc) {
        cov=cov+$2*$5;
    } else {
      print pc,cov;
      cov=$2*$5;
    pc=c}
} END {print pc,cov}' SampleA.smds.coverage | tail -n +2 > SampleA.smds.coverage.percontig

2.2 生成 coverage table

# usage: gen_input_table.py [-h] [--samplenames SAMPLENAMES] [--isbedfiles] fastafile bamfiles [bamfiles ...]
# --samplenames 写有样品名的文件,每个文件名一行
# --isbedfiles  如果在上一步map时运行了genomeCoverageBed,则可以加上此参数后直接用 *smds.coverage文件。如果没运行genomeCoverageBed,则不加此参数,依旧使用bam文件。

# 若多个样品一起做binning,则样品名都需写入sample.txt;若单个样品,则只写入一个即可
for i in `ls $DATA/*.1.fastq`
do
    i=${i/.1.fastq/}
    echo $i
done > sample.txt

python gen_input_table.py --isbedfiles --samplenames sample.txt $DATA/SampleA.fasta ./SampleA.smds.coverage > SampleA_concoct_inputtable.tsv
mkdir ../concoct-input
mv SampleA_concoct_inputtable.tsv ../concoct-input

2.3 生成 linkage table

# usage: bam_to_linkage.py [-h] [--samplenames SAMPLENAMES] [--regionlength REGIONLENGTH] [--fullsearch] [-m MAX_N_CORES] [--readlength READLENGTH] [--mincontiglength MINCONTIGLENGTH] fastafile bamfiles [bamfiles ...]
# --samplenames 写有样品名的文件,每个文件名一行
# --regionlength contig序列中用于linkage的两端长度 [默认 500]
# --fullsearch 在全部contig中搜索用于linkage
# -m 最大线程数,每个ban文件对应一个线程
# --readlength untrimmed reads长度 [默认 100]
# --mincontiglength 识别的最小contig长度 [默认 0]

cd $CONCOCT_EXAMPLE/map
python bam_to_linkage.py -m 8 --regionlength 500 --fullsearch --samplenames sample.txt $DATA/SampleA.fasta ./SampleA.smds.bam > SampleA_concoct_linkage.tsv
mv SampleA_concoct_linkage.tsv ../concoct-input

# 输出文件格式
# 共2+6*i列 (i样品数),依次为contig1、contig2、nr_links_inward_n、nr_links_outward_n、nr_links_inline_n、nr_links_inward_or_outward_n、read_count_contig1_n、read_count_contig2_n
# where n represents sample name. 
# Links只输出一次,如 contig1contig2 输出,则 contig2contig1 不输出

# contig1: Contig linking with contig2
# contig2: Contig linking with contig1
# nr_links_inward: Number of pairs confirming an inward orientation of the contigs -><-
# nr_links_outward: Number of pairs confirming an outward orientation of the contigs <--> 
# nr_links_inline: Number of pairs confirming an outward orientation of the contigs ->->
# nr_links_inward_or_outward: Number of pairs confirming an inward or outward orientation of the contigs. This can be the case if the contig is very short and the search region on both tips of a contig overlaps or the --fullsearch parameter is used and one of the reads in the pair is outside
# read_count_contig1/2: Number of reads on contig1 or contig2. With --fullsearch read count over the entire contig is used, otherwise only the number of reads in the tips are counted.

2.4 运行 concoct

# usage: concoct [-h] [--coverage_file COVERAGE_FILE] [--composition_file COMPOSITION_FILE] [-c CLUSTERS] [-k KMER_LENGTH] [-l LENGTH_THRESHOLD] [-r READ_LENGTH] [--total_percentage_pca TOTAL_PERCENTAGE_PCA] [-b BASENAME] [-s SEED] [-i ITERATIONS] [-e EPSILON] [--no_cov_normalization] [--no_total_coverage] [-o] [-d] [-v]

# --coverage_file 2.2步骤生成的 coverage table
# --composition_file contig的fasta格式文件
# -c, --clusters 最大聚多少类 [默认 400]
# -k, --kmer_length Kmer长度 [默认 4]
# -l, -length_threshold 过滤低于多少的contig序列 [默认 1000]
# -r, --read_length specify read length for coverage, default 100
# --total_percentage_pca PCA分析中可释方差比例
# -b, --basename 输出文件路径 [默认当前路径]
# -s, --seed 聚类的种子数 [0:随机种子,1:默认种子,或其他正整数]
# -i, --iterations VBGMM的最大迭代次数 [默认 500]
# -e, --epsilon VBGMM的e-value [默认 1.0e-6]
# --no_cov_normalization 默认先用对样品做覆盖度均一化,然后对有关的contig做均一化,最后转换日志文件,若设置此参数,则跳过均一化这步
# --no_total_coverage 默认总的覆盖度会新加一列在覆盖度数据矩阵,若设置此参数则跳过这步
# -o, --converge_out 输出文件

cut -f1,3- concoct-input/SampleA_concoct_inputtable.tsv > concoct-input/SampleA_concoct_inputtableR.tsv
concoct -c 40 --coverage_file concoct-input/SampleA_concoct_inputtableR.tsv --composition_file $DATA/SampleA.fasta -b SampleA_concoct-output/

2.5 concoct 结果可视化

mkdir evaluation-output
# 修改ClusterPlot.R中83行的opts(legend.text=theme_text(size=4)为theme(legend.text=element_text(size=4)
Rscript ClusterPlot.R -c SampleA_concoct-output/clustering_gt1000.csv -p SampleA_concoct-output/PCA_transformed_data_gt1000.csv -m SampleA_concoct-output/pca_means_gt1000.csv -r SampleA_concoct-output/pca_variances_gt1000_dim -l -o evaluation-output/SampleA_ClusterPlot.pdf

结果示例

 

有关binns的完整度和污染评估另起博文介绍。

参考:
WECHAT 生信者言
Mande S S, Mohammed M H, Ghosh T S. Classification of metagenomic sequences: methods and challenges[J]. Briefings in bioinformatics, 2012, 13(6): 669-681.
Sangwan N, Xia F, Gilbert J A. Recovering complete and draft population genomes from metagenome datasets[J]. Microbiome, 2016, 4(1): 8.
Cleary B, Brito I L, Huang K, et al. Detection of low-abundance bacterial strains in metagenomic datasets by eigengenome partitioning[J]. Nature biotechnology, 2015, 33(10): 1053.
Wang J, Jia H. Metagenome-wide association studies: fine-mining the microbiome[J]. Nature Reviews Microbiology, 2016, 14(8): 508-522.

宏基因组SOAPdenovo组装

SOAPdenovo特点

宏基因组denovo组装软件中SOAPdenovo(http://soap.genomics.org.cn/soapdenovo)的使用较为广泛,其具有组装速度快,N50值较好的优点。John Vollmers等人总结了8款宏基因组组装软件的相关信息和特点(DOI:10.1371/journal.pone.0169662),如下表。

SOAPdenovo组装原理

DOI: 10.1101/gr.097261.109

A)基因组总DNA被随机打断,采用双端测序。可考虑建150-500bp的小片段文库和2-10kb的大片段文库;
B)Reads读入后,reads都被切割成某一固定Kmer长度的序列(21bp~127bp),基于Kmer采用de Bruijn数据结构构建Reads间的重叠关系;
C) de Bruijn数据结构图中会产生tips (翼尖)、bubbles (气泡)、low coverage links (低覆盖率链接)、tiny repeat (微小重复)等问题

a.删去tips的短末端
b.删去低覆盖率链接
c.调整read路径解决微小重复
d.删去由重复或者二倍体染色体组的混杂造成bubbles结构;

D)在简化后de Bruijn数据结构图的重复节点打断,输出contigs序列;
E)配对的双端reads比对contigs,基于插入片段长度将单一contigs连成scaffolds;
F)通过配对的双端reads来填补scaffolds中的gaps。

故,SOAPdenovo包含6个工具:pregraph、sparse_pregraph、contig、map、scaffold和all,其中all包括前5个工具。

组装过程

准备测序原始数据。

mkdir 0.raw_data
cd 0.raw_data
gzip -dc ../0.initial_data/A_1.fq.gz > A.1.fastq
gzip -dc ../0.initial_data/A_2.fq.gz > A.2.fastq
gzip -dc ../0.initial_data/B_1.fq.gz > B.1.fastq
gzip -dc ../0.initial_data/B_2.fq.gz > B.2.fastq

for i in `ls *.1.fastq`
do
    i=${i/.1.fastq/}
    echo $i
done > ../samples.txt
cd ..

NGSQCToolkit

使用NGSQCToolkit(http://www.nipgr.res.in/ngsqctoolkit.html)去除低质量序列得到Clean Data

1) 去除接头序列;
2) 去除所含低质量碱基(质量值≤13)超过40%的 reads;
3) 去除 N 碱基达到5%的 reads。

NGSQCToolkit包括三部分:QC,Triming,Statistics。使用QC中IlluQC.pl去除低质量序列,使用Triming中AmbiguityFiltering.pl去除N碱基较多序列。

mkdir 1.sequencing_data_quality_control
cd 1.sequencing_data_quality_control

# prepare adaptor cluster
echo "AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT
	CAAGCAGAAGACGGCATACGAGATCGGTCTCGGCATTCCTGCTGAACCGCTCTTCCGATCT
	AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT
	AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT
	CAAGCAGAAGACGGCATACGAGATCGGTCTCGGCATTCCTGCTGAACCGCTCTTCCGATCT
	AGATCGGAAGAGCGGTTCAGCAGGAATGCCGAGACCGATCTCGTATGCCGTCTTCTGCTTG
	TTTTTTTTTTAATGATACGGCGACCACCGAGATCTACAC
	TTTTTTTTTTCAAGCAGAAGACGGCATACGA
	TACACTCTTTCCCTACACGACGCTCTTCCGATCT
	GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT
	TACACTCTTTCCCTACACGACGCTCTTCCGATCT
	AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTA
	GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT
	AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC" > primer_adaptor.txt

# remove low quality reads
for i in `cat ../samples.txt`
do
    echo "/opt/biosoft/NGSQCToolkit_v2.3.3/QC/IlluQC_PRLL.pl -pe ../0.raw_data/$i.1.fastq ../0.raw_data/$i.2.fastq primer_adaptor.txt 1 -l 60 -s 13 -c 8 -o $i"
done  > command.remove_low_quality_reads.list
ParaFly -c command.remove_low_quality_reads.list -CPU 4

# remove non-ATCG reads
for i in `cat ../samples.txt`
do
    echo "/opt/biosoft/NGSQCToolkit_v2.3.3/Trimming/AmbiguityFiltering.pl -i ./$i/$i.1.fastq_filtered -irev ./$i/$i.2.fastq_filtered -p 5"
done  > command.remove_non-ATCG_reads.list
ParaFly -c command.remove_non-ATCG_reads.list -CPU 4
cd ..

使用SOAPdenovo对样品单独组装

1)对于单个样品,组装时选取 Kmer=55,81,101 等不同值进行组装,比较 N50 的优劣,选取最优结果。通常建议Kmer值为reads长度的55%~85%。
2)将组装得到的 Scaffolds 从 N 连接处打断,得到不含 N 的 Scaftigs序列(continuous sequences within scaffolds);
3)采用 SoapAligner 软件将质控后的 Clean Data 比对至各样品组装后的 Scaftigs 上,获取未被利用上的 reads用于后续混合组装;
4)过滤掉单样品的 Scaftigs 中 500bp 以下的片段用于后续分析。

mkdir 2.genome_assembling/SOAPdenovo -p
cd 2.genome_assembling/SOAPdenovo
for i in `cat ../../samples.txt`
do
    echo "ln -s ../../1.sequencing_data_quality_control/$i/$i.1.fastq_filtered_trimmed $i.1.fastq"
    echo "ln -s ../../1.sequencing_data_quality_control/$i/$i.2.fastq_filtered_trimmed $i.2.fastq"
done  > command.copy_data.list
sh command.copy_data.list 

echo "max_rd_len=150
[LIB]
avg_ins=300
reverse_seq=0
asm_flags=3
rd_len_cutoff=150
rank=1
pair_num_cutoff=3
map_len=32
q1=$HOME/2.genome_assembling/SOAPdenovo/A.1.fastq
q2=$HOME/2.genome_assembling/SOAPdenovo/A.2.fastq" > A_config.txt

mkdir A
SOAPdenovo-127mer all -s A_config.txt -o ./A/A -K 81 -d 1 -M 3 -p 8 -R -u -F &> A_soapdenovo.log

cd A
fasta_no_blank.pl A.scafSeq > A.scafSeq.noblank && cp A.scafSeq.noblank A.Scaftigs && perl -p -i -e "s/N+/\r\n>\r\n/g" A.Scaftigs
rename_fasta_identifiers.pl A.Scaftigs > A.Scaftigs.rename
genome_seq_clear.pl --seq_prefix A_Scaftigs_500 --min_length 500 A.Scaftigs.rename > A_Scaftigs_500.fasta

cd ..

#mapping to get unpaired reads
/opt/biosoft/SOAPaligner/soap2.21release/2bwt-builder A.Scaftigs.rename
/opt/biosoft/SOAPaligner/soap2.21release/soap -a ../../../2.genome_assembling/SOAPdenovo/A.1.fastq -b ../../../2.genome_assembling/SOAPdenovo/A.2.fastq -D A.Scaftigs.rename.index -o A_PE_output -u A_unmap_output -2 A_SE_output -m 200 -p 20

echo "max_rd_len=150
[LIB]
avg_ins=300
reverse_seq=0
asm_flags=3
rd_len_cutoff=150
rank=1
pair_num_cutoff=3
map_len=32
q1=$HOME/2.genome_assembling/SOAPdenovo/B.1.fastq
q2=$HOME/2.genome_assembling/SOAPdenovo/B.2.fastq" > B_config.txt

mkdir B
SOAPdenovo-127mer all -s B_config.txt -o ./B/B -K 81 -d 1 -M 3 -p 8 -R -u -F &> B_soapdenovo.log

cd N
fasta_no_blank.pl B.scafSeq > B.scafSeq.noblank && cp B.scafSeq.noblank B.Scaftigs && perl -p -i -e "s/N+/\r\n>\r\n/g" B.Scaftigs
rename_fasta_identifiers.pl B.Scaftigs > B.Scaftigs.rename
genome_seq_clear.pl --seq_prefix B_Scaftigs_500 --min_length 500 B.Scaftigs.rename > B_Scaftigs_500.fasta

#mapping to get unpaired reads
/opt/biosoft/SOAPaligner/soap2.21release/2bwt-builder B.Scaftigs.rename
/opt/biosoft/SOAPaligner/soap2.21release/soap -a ../../../2.genome_assembling/SOAPdenovo/B.1.fastq -b ../../../2.genome_assembling/SOAPdenovo/B.2.fastq -D B.Scaftigs.rename.index -o B_PE_output -u B_unmap_output -2 B_SE_output -m 200 -p 20

cd ..

Unmapped Reads混合组装

没有参与组装的测序数据被称为unmapping reads,一般不会参与后续分析。unmapping reads的舍弃直接导致宏基因组数据利用的不完整性。因此,Qin J等提出混合组装。其可以充分利用测序数据,同时有助于挖掘低丰度物种信息。

1)单个样品中未被利用上的reads用于后续混合组装,参数同单个样品组装;
2)将混合组装的 Scaffolds 从 N 连接处打断,得到不含 N 的 Scaftigs 序列;
3)过滤掉混合组装的 Scaftigs 中 500bp 以下的片段用于后续分析。

mkdir Mix
cd Mix
cat $HOME/2.genome_assembling/SOAPdenovo/A/A_unmap_output $HOME/2.genome_assembling/SOAPdenovo/B/B_unmap_output > Mix_seq
echo "max_rd_len=150
[LIB]
avg_ins=300
reverse_seq=0
asm_flags=3
rd_len_cutoff=150
rank=1
pair_num_cutoff=3
map_len=32
p=$HOME/2.genome_assembling/SOAPdenovo/Mix_seq" > Mix_config.txt

SOAPdenovo-127mer all -s Mix_config.txt -o ./Mix -K 81 -d 1 -M 3 -p 8 -R -u -F &> Mix_soapdenovo.log

fasta_no_blank.pl Mix.scafSeq > Mix.scafSeq.noblank && cp Mix.scafSeq.noblank Mix.Scaftigs && perl -p -i -e "s/N+/\r\n>\r\n/g" Mix.Scaftigs
rename_fasta_identifiers.pl Mix.Scaftigs > Mix.Scaftigs.rename
genome_seq_clear.pl --seq_prefix Mix_Scaftigs_500 --min_length 500 Mix.Scaftigs.rename > Mix_Scaftigs_500.fasta

参数说明

-s  string      配置文件
-p  int         程序运行时设定的cpu线程数,默认值[8]
-o  string      输出的文件名前缀 
-K  int         输入的K-mer值大小,K-mer值必须是奇数,默认值[23]
-a  int         假设初始化内存,为了避免进一步的再分配,默认[0G]
-d  int         为了减少错误测序的影响,去除kmers频数不大于该值(kmerFreqCutoff)的k-mer,默认值[0] 
-R  (optional)  利用read鉴别短的重复序列,默认值[NO]
-D  int         为了减少错误测序的影响,去除频数不大于该值的由k-mer连接的边,默认值[1]
-M  int         连contiging时,合并相似序列的强度,最小值0,最大值3,默认值为[1]
-e  int         两边缘之间的弧的权重大于该值,将被线性化,默认值[0]
-m  int         最大的kmer size(max 127)用于multi-kmer,默认[NO]
-E  (optional)  在iterate迭代之前,合并clean bubble功能,仅在当使用multi-kmer时且设置-M参数,默认[NO]
-F  (optional)  利用read对scaffold中的gap进行填补,默认[NO]
-u  (optional)  构建scaffolding前不屏蔽高/低覆盖度的contigs,高频率覆盖度指平均contig覆盖深度的2倍。默认[屏蔽]
-w  (optional)  在scaffold中,保持contigs弱连接于其他contigs,默认[NO]
-G  int         估计gap的大小和实际补gap的大小的差异值,默认值[50bp]
-L  int         用于构建scaffold的contig的最短长度(minContigLen),默认为:[Kmer参数值+2]
-c  float       在scaffolding之前,contigs小于100bp,且覆盖率小于该最小contig覆盖率(c*avgCvg),将被屏蔽,除非参数-u已设定,默认值[0.1]
-C  float       在scaffolding之前,contigs覆盖率大于该最大contigs覆盖率(C*avgCvg),或者contigs小于100bp且覆盖率大于0.8*(C*avgCvg),将被屏蔽,除非参数-u已设定,默认值[2]
-b  float       当处理contigs间的pair-end连接时,如果参数>1,该插入片段的上限值(b*avg_ins)将被用来作为大的插入片段(>1000)的上限,默认值[1.5]
-B  float       如果两个contigs的覆盖率都小于bubble覆盖率(bubbleCoverage)乘以contigs平均覆盖率(bubbleCoverage*avgCvg),则去除在bubble结构中的低覆盖率contig,默认值[0.6]
-N  int         统计基因组大小,默认值[0]
-V  (optional)  输出相关信息(Hawkeye)为了可视化组装,默认[NO]

配置文件

max_rd_len=150        #read的最大长度
[LIB]                 #文库信息以此开头
avg_ins=300           #文库平均插入长度
reverse_seq=0         #序列是否需要被反转,插入片段大于等于2k的采用了环化,所以对于插入长度大于等于2k文库,序列需要反转,reverse_seq=1,小片段设为0
asm_flags=3           #设为1:reads只用于只构建contig;设为2:reads只用于构建scaffold;设为3:reads同时用于构建contig和scaffold;设为4:reads只用于补洞gap
rank=1                #rank该值取整数,值越低,reads越优先用于构建scaffold
pair_num_cutoff=3     #该参数规定了连接两个contig或者是pre-scaffold 的可信连接的阈值,即,当连接数大于该值,连接才算有效。小片段文库(<2k)默认值[3],大片段文库默认值[5]。
map_len=32            #该参数规定了在map过程中reads和contig的比对长度必须达到该值才能作为一个可信的比对。小片段文库一般设32,大片段文库一般设35,默认值[K+2]。
q1=$DATA_PATH/read_1.fq
q2=$DATA_PATH/read_2.fq #双端fastq序列,写绝对路径
f1=./pathfasta_read_1.fa
f2=./pathfasta_read_2.fa #双端fasta序列,写绝对路径
q=./pathfastq_read_single.fq #单向fastq
f=./pathfasta_read_single.fa #单向fasta
p=./pathpairs_in_one_file.fa #双向测序得到的一个fasta格式的序列文件

输出文件

SOAPdenovo中四部分别对应的输出文件:
pregraph  *.kmerFreq *.edge *.preArc *.markOnEdge *.path *.vertex *.preGraphBasic
contig   *.contig *.ContigIndex *.updated.edge *.Arc
map     *.readOnContig *.peGrads *.readInGap
scaff    *.newContigIndex *.links *.scaf *.scaf_gap *.scafSeq *.gapSeq

关键文件如下:
*.contig            #没有使用连Scaffold的contig序列
*.scafSeq          #SOAPdenovo软件最终的组装序列结果
*.scafStatistics      #contigs和scaffolds的最终统计信息
*.Scaftigs          #打断N后的Scaftigs序列
*.Scaftigs_500.fasta  #长度大于500bp的Scaftigs序列