JGI 数据命令行下载

JGI 提供了三种数据下载方式:

  1. 网页直接下载;
  2. 通过 Globus 下载;
  3. 通过 API 命令行下载;

命令行下载方法:
第一步:生成cookies文件,该文件中保存着用户名和密码

curl 'https://signon.jgi.doe.gov/signon/create' --data-urlencode 'login=USER_NAME' --data-urlencode 'password=USER_PASSWORD' -c cookies
# USER_NAME 和 USER_PASSWORD 换成登陆用户名和密码

第二步:获取下载地址
Genome Portal 中搜索样品名,进入下载界面如下。点击 “Open Downloads as XML” 后得到 XML文件。

在XML文件中提取下载地址(url=”下载地址”)
如:

/ext-api/downloads/get_tape_file?blocking=true&url=/MetspLW13_2/download/_JAMO/531a42f549607a1be0056c4f/submission.assembly.fasta

第三步:下载文件

curl 'https://genome.jgi.doe.gov/ext-api/downloads/get_tape_file?blocking=true&url=/MetspLW13_2/download/_JAMO/531a42f549607a1be0056c4f/submission.assembly.fasta' -b cookies > 168131.assembled.KO

 

根据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没有被切除,需要后续进一步处理。

JGI workshop

本想着访学后期能轻松地参加一场workshop,能学点东西最好,不能学到啥就当休假了。可是遭遇了半年的实验滑铁卢,最近实验才重新有点起色。在实验的关键时期,带着忐忑的心情,来到JGI。JGI总部位于加州Walnut Creek市,看起来挺新的一个城市。每每看到攒新的马路时,我就想起湘潭那些修好挖了又修的路。住宿相比当初美西游寒酸的 motel 好多了,至少看起来是个正规的酒店。

本以为JGI会围绕宏基因组分析流程逐步展开,没想到讲的跟command一点关系都没有。全程围绕JGI团队开发IMG系统,学习在线交互操作。如果要说有什么收获,我想有两点。一来接触很多以前不熟悉的前沿方向,比如宏病毒组(David Paez做了很多不错的工作)。二来JGI提供的吃的真是丰富,从早餐,上午茶,中餐,下午茶,生怕你饿着。

这次workshop有来自各个国家和地区的,有趣的是居然还碰到一个来自当初德国课题组的博士生(当然她依旧逃不过Kappler做硝酸盐还原铁氧化三个生境的故事)。

JGI还贴心地给大家准备了证书。

现在主流的pacbio测序仪,相比国内动则几十台的阵势弱了点。

新款的pacbio测序仪Sequel,JGI正在调试中。据说诺禾买了20台了,不知道数据质量如何。

Illumina公司新推出的novoseq,她的通量贼大,我主要觉得她外观长得好看。

最后这个看似不起眼的桌面测序仪就是现在很火的nanopore,我觉得这才是第三代测序的未来啊,不用扩增建库,上来就直接测,读长可以到几十kb。这个水平下,错误率就不再是一个考虑的必要条件了。

 

PS. 最后要给加州的空气正名,加州也是有蓝天的。

宏基因组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.

奥林匹克轰趴游

奥林匹克公园作为西雅图片区三大国家公园之一,已经心心念念了许久。这次咱9个小伙伴租了个house,在奥林匹克公园浪了三天两晚。1788年,英国船长 John Meares 航海经过此地,看到壮丽的山峰和皑皑的积雪,觉得这必定是希腊神话中众神之寓所,因而此地得名 Mount Olympus,Mount Olympus 海拔高达2428米,是园区内最高峰,同时也是一坐活火山。奥林匹克公园坐落于华盛顿州西北的奥林匹克半岛上(其中 Cape Flattery 海滩是美国本土最西北角),是美国除了阿拉斯加以外最大的原始地区。奥林匹克公园被誉为露营者的天堂,园区内散布着大大小小十几个营地。园区地处温带,但从太平洋海面上吹来的温暖而潮湿的空气,被 Olympus 山脉挡住,丰富的雨水和融雪形成了山腰处独特的雨林生态。特殊的地理位置造就了奥林匹克公园多样化的景观,朗阔雪山、雨林、海滨多种景色。依旧老规矩,先上图,再来攻略。

ANI计算

ANI概念

Average Nucleotide Identity (ANI) 是在核酸水平,两两基因组之间所有直系同源蛋白编码基因的相似性,常用于研究基因组之间的进化距离。相较于传统的 DNA-DNA hybridization (DDH),ANI指标计算简单省时。且ANI指标有助于构建结构性数据库,方便生物信息学者的后续研究。还有一种计算ANI的方法是基于俩俩基因组的共同同源蛋白编码基因的相似性。接下来介绍几款计算ANI的在线工具和本地软件(Jspecies 和 pyani)。

ANI在线工具

Kostas lab:http://enve-omics.ce.gatech.edu/ani/index
MicrobeTrakr 的 ANI Calculator:http://www.microbetrakr.org/ani/
Ezbiocloud:http://www.ezbiocloud.net/tools/ani

ANI本地软件

Jspecies

http://imedea.uib-csic.es/jspecies/  DOI: 10.1073/pnas.0906412106
Jspecies是可视化操作,用户体验简单。但输入基因组相对麻烦,不方便较多基因组之间的两两计算ANI。

Jspecies安装

# install blast
mkdir blast
cd blast
wget http://mirrors.vbi.vt.edu/mirrors/ftp.ncbi.nih.gov/blast/executables/legacy/2.2.26/blast-2.2.26-universal-macosx.tar.gz
tar zxvf blast-2.2.26-universal-macosx.tar.gz -C /opt/biosoft
# blast path: /opt/biosoft/blast-2.2.26/bin/blastall
# formatdb path: /opt/biosoft/blast-2.2.26/bin/formatdb
# fastacmd path: /opt/biosoft/blast-2.2.26/bin/fastacmd
cd ..
# install MUMmer
mkdir MUMmer
cd MUMmer
wget https://sourceforge.net/projects/mummer/files/mummer/3.23/MUMmer3.23.tar.gz
tar zxvf MUMmer3.23.tar.gz -C /opt/biosoft
cd /opt/biosoft/MUMmer3.23
make check
make install
# nucmer path: /opt/biosoft/MUMmer3.23/nucmer
# install Jspecies
mkdir /opt/biosoft/jspecies
cd /opt/biosoft/jspecies
wget http://imedea.uib-csic.es/jspecies/jars/jspecies1.2.1.jar
java -jar -Xms1024m -Xmx1024m jspecies1.2.jar

Jspecies运行

  • 打开软件图形界面后,点击 Edit 下的 preferences,添加 blast、nucmer、formatdb、fastacmd路径;
  • 设置 ANIb (基于blast得到的ANI值)
 -X 空比对对应的下降值 [默认 150]
 -q 碱基错配罚分 [默认 -1]
 -F 是否过滤query序列的简单重复和低复杂度序列 [F]
 -e e-value [默认 1e-15]
 -a 运行线程数 [默认 2]
 Identity (%) >= 30
 Alignment (%) >= 70 #more than 30% overall sequence identity (recalculated to an identity along the entire sequence) over an alignable region of at least 70% of their length
 length 1020 #query序列被连续切割成1020bp的片段,然后将这些片段blastn比对到reference序列,比对的结果用于后续计算
  • ANIm (基于MUMmer得到的ANI值)参数
       
 -b 低得分区大于此值时,停止比对 [默认 200]
 -c 一簇匹配的序列最少需要有大于此长度的匹配。此值越大,则小的匹配越少 [默认 65]
 -g 在一簇匹配的序列中,两个邻近的matches之间的gap的最大阈值。此值越大,则越容易找到更大匹配区 [默认 90]
 -l single match的最小长度 [默认 20]
  • 导入数据并运行:File -> New  File -> Import FASTA(S) from files -> 选择 Tetra、ANIb、ANIm -> Start

ANI命令行软件

pyani

pyani是用于计算ANI的python3模块,其依赖包有numpy、scipy、matplotlib、biopython、pandas、seaborn。0.1.3.2及以上版本的其他依赖包不用另外安装,0.1.3.2以下版本 pip install -r requirements.txt 安装依赖包。pyani 包括average_nucleotide_identity.py (计算ANI) 和 genbank_get_genomes_by_taxon.py (从NCBI下载数据),并提供四种不同算法:

  • ANIm:使用 MUMmer (NUCmer) 比对序列
  • ANIb:使用 BLASTN+ to 比对序列的 1020nt 片段
  • ANIblastall:使用老版的 BLASTN 比对序列的 1020nt 片段
  • TETRA:计算每个序列的四核苷酸的频率

pyani安装

# pip3 install
pip3 install pyani
echo "PATH=/home/linuxbrew/.linuxbrew/bin:$PATH" >> ~/.bashrc
source ~/.bashrc
# local install
git clone git@github.com:widdowquinn/pyani.git
cd pyani
python3 setup.py install

pyani运行

# example
average_nucleotide_identity.py -i tests/test_ani_data/ -o tests/test_ANIm_output -m ANIm -g
average_nucleotide_identity.py -i tests/test_ani_data/ -o tests/test_ANIb_output -m ANIb -g
average_nucleotide_identity.py -i tests/test_ani_data/ -o tests/test_ANIblastall_output -m ANIblastall -g
average_nucleotide_identity.py -i tests/test_ani_data/ -o tests/test_TETRA_output -m TETRA -g

pyani参数

# required
  -i INDIRNAME, --indir INDIRNAME    输入文件所在目录
  -o OUTDIRNAME, --outdir OUTDIRNAME 输出结果目录

# optional
  -v, --verbose         输出详细结果
  -f, --force           强制文件的写入
  -s FRAGSIZE, --fragsize FRAGSIZE
                        ANIb比对时序列的片段大小 [默认 1020]
  -l LOGFILE, --logfile LOGFILE
                        输出日志文件
  --skip_nucmer         跳过运行NUCmer
  --skip_blastn         跳过运行BLASTN
  --noclobber           不重做已存在文件
  --nocompress          不压缩或删除比较结果
  -g, --graphics        生成ANI的热图
  --gformat GFORMAT     图形文件输出格式 [pdf|png|jpg|svg] [默认 pdf,png,eps]
  --gmethod {mpl,seaborn} 图形输出方法 [默认 mpl]
  --labels LABELS       序列标记文件的路径
  --classes CLASSES     序列分类文件的路径
  -m {ANIm,ANIb,ANIblastall,TETRA}, --method {ANIm,ANIb,ANIblastall,TETRA}
                        ANI比对方法 [默认 ANIm]
  --scheduler {multiprocessing,SGE}
                        Job scheduler [默认 multiprocessing]
  --workers WORKERS     运行线程数 [默认 全部线程]
  --SGEgroupsize SGEGROUPSIZE
                        SGE组中job数 [默认 10000]
  --SGEargs SGEARGS     Additional arguments for qsub
  --maxmatch            Override MUMmer to allow all NUCmer matches
  --nucmer_exe NUCMER_EXE
                        NUCmer 程序路径
  --filter_exe FILTER_EXE
                        delta-filter 程序路径
  --blastn_exe BLASTN_EXE
                        BLASTN+ 程序路径
  --makeblastdb_exe MAKEBLASTDB_EXE
                        BLAST+ makeblastdb 程序路径
  --blastall_exe BLASTALL_EXE
                        BLASTALL 程序路径
  --formatdb_exe FORMATDB_EXE
                        BLAST formatdb 程序路径
  --write_excel         输出Excel格式文件
  --rerender            Rerender graphics output without recalculation
  --subsample SUBSAMPLE
                        Subsample a percentage [0-1] or specific number (1-n)
                        of input sequences
  --seed SEED           Set random seed for reproducible subsampling.
  --jobprefix JOBPREFIX
                        Prefix for SGE jobs (default ANI).

pyani自带绘图,若不满意,可以调用ggplot2绘制热图

Rscript heatmap_ANI.R $OUTPUT/ANIb_percentage_identity.tab
#!/usr/bin/env Rscript
library("ggplot2")
library("reshape")

userprefs <- commandArgs(trailingOnly = TRUE)
file_name <- userprefs[1]

# read file
data.ANI <- read.table(file_name)

# cluster and order matrix
row.order <- hclust(dist(data.ANI))$order # clustering
col.order <- hclust(dist(t(data.ANI)))$order
dat_new <- data.ANI[row.order, col.order] # re-order matrix accoring to clustering

# melt to dataframe
df_molten_dat <- melt(as.matrix(dat_new)) # reshape into dataframe
names(df_molten_dat)[c(1:3)] <- c("shewanella1", "shewanella2", "ANI")

# make plot
p1 <- ggplot(data = df_molten_dat,
aes(x = shewanella1, y = shewanella2, fill = ANI)) +
geom_raster() +
scale_fill_distiller(palette = "RdYlBu", limits=c(0.8, 1)) +
theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
ggtitle("Shewanella ANI heatmap")+
geom_text(aes(label = round(ANI, 2)))+
theme(axis.title=element_text(size=16), strip.text.x=element_text(size=16),
legend.title=element_text(size=15), legend.text=element_text(size=14),
axis.text = element_text(size=14), title=element_text(size=20),
strip.background=element_rect(fill=adjustcolor("lightgray",0.2))
)

# export plot
png("Heatmap_ANI_highres.png", width=15, height=15, units="in", pointsize=12, res=500)
p1
dev.off()

示例图效果

24种不同 shewanella 属菌种ANI比较热图

参考:
http://imedea.uib-csic.es/jspecies/
https://github.com/widdowquinn/pyani/
http://www.a-site.cn/article/300418.html
https://github.com/rprops/MetaG_lakeMI/wiki/

宏基因组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序列

美西四人行

趁着6月艳阳高照,和Peter组其他3个小伙伴一起在美西逛了一圈。从西雅图飞旧金山,然后自驾沿着一号公路一直开到洛杉矶,再从洛杉矶转站拉斯维加斯,最后结束一周的“空虚”的生活回到西雅图。回来最深的感受是,相比美西其他的大城市,还是咱西雅图最好。旧金山和洛杉矶空气灰蒙蒙的,跟咱国内雾霾有的一拼。加州6月的阳光强烈,空气异常干燥,不像西雅图这样养人。这一路主要成了专职摄影师,来几张路途的风景。

金门大桥

斯坦福大学

一号公路海滩

拉斯维加斯百丽宫赌场

西峡谷

胡佛水坝

拉斯维加斯晚霞

最后贴上我们的攻略,供想去的朋友参考。