病原微生物高通量测序数据分析笔记¶
注解
本电子书为基于BSD授权，版权归作者所有；欢迎下载及编辑(个人用途)，但未经作者同意必须保留此段声明，且不可用于商业用途，否则保留追究法律责任的权利。
- 作者：indexofire <indexofire@gmail.com>
- Github地址：https://github.com/indexofire/ndafp
关于本笔记¶
随着高通量测序技术快速发展，小型台式测序仪的出现降低了高通量测序的门槛。对于开展微生物特别是病原微生物的测序工作的需求，在各种临床，科研以及政府的小型实验室中迅速扩展。目前已经上市的比较适合小型实验室的测序仪如 Illumina Miseq，Illumina NextSeq500，Ion S5，Ion PGM 等，即将上市的测序仪如 Illumina MiniSeq，Qiagen GeneReader。除了传统的二代测序仪，不需要扩增的基于单分子测序技术的三代测序技术的 PacBio Sequel 和最近开展埃博拉实时测序项目的基于纳米孔技术的 Oxford MinION 产品都是非常适合开展病原微生物测序的平台。
测序数据的不断产出，需要生物信息学人才来参与到测序数据的分析以便获得有价值的生物学信息对于许多小型实验室来说，很难招募到优秀的生物信息学人才。对于非专业背景的人来说，命令行以及软件的安装可能是学习中遇到的首要难题。本笔记主要是通过介绍微生物，特别是病原细菌的高通量测序数据分析，帮助原来从事湿实验领域的同事能较好的了解和掌握高通量测序数据分析，将学习曲线变得平缓。同时也希望能对从事该领域的其他工作者提供便利，让大家在需要分析微生物基因组数据时，能有一个资料与实例汇集的手册可以快速参考和学习。
本笔记主要针对 Miseq 机器的产出数据开展分析，由于原理类似分析软件也大都可以应用在 Illumina 其他型号的仪器数据上。例子都是运行在安装了 Ubuntu 发行版的 Linux PC 服务器或者 Ubuntu 的个人电脑上。
本笔记使用开源工具 Sphinx-doc 创作，并在 Read the Docs 网站发布。笔记中的许多内容都是来源于网络上的资料，部分来源可能有误或记忆不全，如果原作者发现没有内容链接，或者链接错误，请发送 电子邮件 通知修改。同时希望 Open Source 思想能对科研工作和科研工作者有所帮助。
本笔记代码托管在 Github, 所有内容可以在 http://github.com/indexofire/ndafp.git 获得。笔者接触 NGS 时间不长，内容中肯定会有许多错误之处，希望能有NGS领域的专家帮助指正。对于 笔记中的错误会内容更新，可以提交 issues，帮助完善笔记内容。
Fork笔记方法¶
- 在GitHub上访问本书源码，并fork成为自己的仓库。如下图所示，点击fork按钮，然后用 git clone 到本地，设置好代码仓库用户信息。就可以修改并提交代码了。
~$ git clone git@github.com:your_github_username/ndafp.git ~$ cd ndafp ~/ndafp$ git config user.name "your github username" ~/ndafp$ git config user.email your_email@something.com
对需要的内容做修改后提交，并 git push 到之前 fork 的仓库。
~/ndafp$ git add -A ~/ndafp$ git commit -am "Fix issue #1: change typo: helo to hello" ~/ndafp$ git push这样你自己的代码分支 master 与源代码分支就会出现差异，要合并这些差异，就要在 GitHub 网站上提交 pull request。如下图所示如果有代码。
由于源代码仓库也会经常更新或者合并其他 pull request 代码，需要定期使用源项目仓库内容来更新自己仓库内容。
~/ndafp$ git remote add upstream https://github.com/indexofire/ndafp.git ~/ndafp$ git fetch upstream ~/ndafp$ git checkout master ~/ndafp$ git merge upstream/master ~/ndafp$ git push origin master
内容目录¶
序言¶
Linux 入门基础¶
1. Linux 系统结构¶
2. Linux 权限概念¶
3. Linux 常用的命令¶
笔记中使用许多 Linux 命令行，这里做一个快速笔记以备查询。以 bash 为例。其他的 shell 比如 csh, tcsh 等大部分命令和参数都是一致的。
ls: 列出当前目录内容
# 显示 '~/data' 文件夹里的文件及目录
$ ls ~/data
# 显示 '~/data' 文件夹里的所有文件及目录（包括隐含文件），并以列表形式显示
$ ls -la ~/data
# 文件按照修改时间排序，新文件在前
$ ls -t ~/data
# 查看所有 .fastq 后缀的文件，* 代表任何长度的任意字符
$ ls ~/data/*.fastq
cd: 改变文件夹
# 从当前文件夹 '~/data' 更换到 '~/app' 文件夹
~/data$ cd ~/app
# 从任何位置切换回用户主目录 '/home/User'
$ cd
# 回退上一个访问目录
$ cd -
# 返回上一级目录
$ cd ..
# 返回多级目录并前往一个叫“new directory”的文件夹，这里 '\' 是转义符，将空格符转义
$ cd ../../../new\ directory
pwd: 显示当前目录
# 当不知道当前处于什么路径时，可以用这个命令显示
$ pwd
mkdir: 建立新文件夹
# 当前路径下新建一个名叫 'new' 的文件夹
$ mkdir new
rmdir: 删除文件夹
# 删除当前路径的文件夹 'new'
$ rmdir new
rm: 删除文件
# 删除当前文件夹的所有后缀是 '.sra' 的文件
$ rm *.sra
# 删除文件夹 'new' 以及 'new' 下的所有子文件与子目录
$ rm -R new
# 不弹出删除确认提示，删除所有 '.tmp' 文件
$ rm -f *.tmp
cp: 复制文件
# 复制 'test.txt' 文件到文件夹 '~/abc' 中
$ cp test.txt ~/abc
mv: 移动文件或文件夹
# 移动 'test.txt' 文件到文件夹 '~/abc' 中并改名叫 'test1.txt'
$ mv test.txt ~/abc/test1.txt
which: 查找可执行文件的系统路径
# 打印出系统带的 python 程序的路径
$ which python
wc: 统计一个文件的行，字符和字节数
# 输出文件 'text.txt' 的行数，字符数和字节数。
$ wc text.txt
cat: 输出文件内容
# 显示文件 'text.txt' 内容
$ cat text.txt
grep: 截取输入字符的选定 pattern 并输出所在的行
# 显示文件 'text.txt' 中含有字符 'abc' 的行
$ cat text.txt | grep 'abc'
head: 输出文件头部内容
# 输出文件前5行内容
$ head -5 text.txt
tail: 输出文件尾部内容
# 输出文件的最后5行内容
$ tail -5 text.txt
ps: 查看系统进程
ps会在终端打印系统进程，各列的含义是:
- USER: 运行该进程的用户
- PID: 运行着的命令(CMD)的进程编号
- %CPU: CPU占用
- %MEM: 内存占用
- VSC:
- RSS:
- TTY: 命令所运行的位置（终端）
- STAT:
- TIME: 运行着的该命令所占用的CPU处理时间
- COMMAND: 该进程所运行的命令
# 显示详细的进程信息
$ ps -waux
# 过滤用户root的进程
$ ps -u root
# 根据不同参数使用来排序进程，并只现实排名前10的进程
$ ps -aux --sort -pcpu | head -n 11
$ ps -aux --sort -pmem | head -n 11
$ ps -aux --sort -pcpu,+pmem | head -n 11
# 过滤进程名
$ ps -f -C chrome
# 根据PID过滤
$ ps -L 1000
# 树形现实进程
$ ps -axjf
4. 其他常用命令与工具¶
lsb_release: 查看发行版信息
$ lsb_release -a
# 打印出系统发行版信息
No LSB modules are available.
Distributor ID: Ubuntu
Description: Ubuntu 16.04.1 LTS
Release: 16.04
Codename: xenial
htop: 系统进程查看
# Ubuntu 发行版默认不带 htop，需要安装后使用
$ sudo apt install htop
# 运行 htop
$ htop
dmesg: 查看系统日志
# 显示启动状态开始的核心输出日志
$ dmesg | less
# 列出某个硬件的输出信息，如默认磁盘 sda
$ dmesg | grep 'sda'
# 清空缓存区信息
$ dmesg -c
了解 Shell¶
脚本编程 Python¶
脚本编程 Perl¶
普通 Linux 发行版都会默认安装 Perl，可以用 perl -v 查看所安装的 perl 版本。
$ perl -v
检查是否安装了 perl 模块可以用下面方法查看，如果安装了模块则没有打印输出，否则会提示找不到模块。
$ perl -e 'use Bio::Seq'
用 CPAN 安装模块
$ perl -MCPAN -e shell
>install HTML::TokeParser::Simple
用 cpan 脚本安装模块
$ cpan HTML::TokeParser::Simple
awk 工具来帮助处理 NGS 数据¶
awk 工具常用来对文件内容按行进行匹配搜索，一旦某一行中的内容与搜索条件匹配，则对该内容做进一步操作。本质上讲，awk 属于数据驱动的程序，因此 awk 程序也非常容易编写与阅读。高通量测序技术带来大量的数据，而生物学分析往往直需要其中特定的部分，这时 awk 就可以大显身手为测序数据的转换与分析提供帮助。
注解
awk 程序内容很简单，主要包括 pattern，action, input_file
- pattern: 表示所要搜索的内容，可以用正则表达式
- action: 则表示搜索匹配后要做的操作
- input_file: 包含所要搜索内容的输入文件
awk 可以在 shell 里直接运行：
~$ awk 'pattern { action }' input_file
也可以写成脚本程序来运行。一个 awk 程序往往如下所示
#!/usr/bin/awk -f
pattern { action }
awk 可以不需要输入文件；对于pattern和action来说，2者至少要有一个才能运行。如果没有pattern，则默认匹配任何输入，按行输出并执行action。如果没有action，则匹配pattern并按行输出不做额外操作。
# 尝试以下语句，在终端打印`hello world`。这里BEGIN是pattern，{ print "hello world" } 是action
~$ awk 'BEGIN { print "hello world" }'
# 模拟cat输出终端输入的字符。这里省略了pattern
~$ awk '{ print }'
上面这个例子可以保存成脚本来运行，文件内容如下：
#!/usr/bin/env awk -f
BEGIN { print "hello world" }
运行这个程序：
~$ chmod +x hello_world
~$ ./hello_world
1. 了解 awk 的基本用法¶
~$ awk '$1~/gene1/ { print $2 }' input_data
作用与 cat + grep + sed 类似
~$ cat input_data | grep "gene1" | sed {$1}
2. awk 脚本语法¶
2.1 Pattern¶
pattern 可以采用正则表达式来匹配内容。用/ regular expression /来表示。
注解
几个特殊的Patterns：
- BEGIN:
- END:
- BEGINFILE:
- ENDFILE:
~$ awk '$1 == "gene" { print $2 }' genes.list
2.2 Action¶
# if block
~$ awk '{if(condition) print $1}' input
# if else block
~$ awk '{if(condition)
> print $1
> else
> print $2
> }' input
while 控制语句要换行，用4个空格划分控制块。
~$ awk '{ i = 1
> while (i <= 3) {
> print $i
> i++
> }
> }' input
2.3. 内建变量¶
- BINMODE #
- CONVFMT
- FIELDWIDTHS #
- FPAT #
- FS
- IGNORECASE #
- LINT #
- OFMT
- OFS
- ORS
- PREC #
- ROUNDMODE #
- RS
- SUBSEP
- TEXTDOMAIN #
- ARGC , ARGV
- ARGIND #
- ENVIRON
- ERRNO #
- FILENAME
- FNR
- NF
- FUNCTAB #
- NR
- PROCINFO #
- RLENGTH
- RSTART
- RT #
- SYMTAB #
2.4. 内建函数¶
- atan2(y, x)
- cos(x)
- exp(x)
- int(x)
- log(x)
- rand()
- sin(x)
- sqrt(x)
- srand([ x ])
- asort(source [ , dest [ , how ] ] ) #
- asorti(source [ , dest [ , how ] ] ) #
- gensub(regexp, replacement, how [ , target ] ) #
- gsub(regexp, replacement [ , target ] )
- index(in, find)
- length( [ string ] )
- match(string, regexp [ , array ] )
- patsplit(string, array [ , fieldpat [ , seps ] ] ) #
- split(string, array [ , fieldsep [ , seps ] ] )
- sprintf(format, expression1, ...)
- strtonum(str) #
- sub(regexp, replacement [ , target ] )
- substr(string, start [ , length ] )
- tolower(string)
- toupper(string)
- close(filename [ , how ] )
- fflush( [ filename ] )
- system(command)
- mktime(datespec)
- strftime( [ format [ , timestamp [ , utc-flag ] ] ] )
- systime()
- and(v1, v2 [, ...])
- compl(val)
- lshift(val, count)
- or(v1, v2 [, ...])
- rshift(val, count)
- xor(v1, v2 [, ...])
- isarray(x)
- bindtextdomain(directory [ , domain ] )
- dcgettext(string [ , domain [ , category ] ] )
- dcngettext(string1, string2, number [ , domain [ , category ] ] )
3. 示例¶
示例部分介绍 awk 在日常中，特别是处理 ngs 数据时的一些例子。
3.1 对列的操作¶
调整列顺序：在有GUI的操作系统里，一般采用类似 excel, calc 之类的软件导入数据文件，然后剪切各列调整顺序。如果用 awk 来解决也很方便，你只需要考虑好调整的各列顺序即可，action 里的{ print ...}顺序就是重新调整后的各列顺序：
~$ awk '{ print $3, $5, $7, $2, $1, $4, $6 }' infile.txt > outfile.txt
插入列：有时候要插入一列数据，用awk可以很方便的实现：
~$ awk '{ print $1, $2 "gene expression", $3}' infile.txt > outfile.txt
3.2 对行的操作¶
去除重复的行：有时候数据里含有重复的行，而当你只需要唯一性数据时，就可以用这行程序，只保留具有唯一性的数据行。
~$ awk '!x[$0]++' infile.txt > outfile.txt
3.3 格式化输出¶
~$ awk -F":" 'BEGIN { printf "%-8s %-3s", "User", "UID" } \
> NR==1,NR==20 { printf "%-8s %-3d", $1, $3 }' /etc/passwd
3.4 合并文件¶
4. Bioawk¶
bioawk 是 Heng Li 开发的 awk 扩展工具，增加了对压缩的 BED, GFF, SAM, VCF, FASTA/Q 等文件格式的支持，并内建一些函数，适用于NGS数据的快速输入输出。
内建以下常用的测序数据函数：
- gc($seq): 返回序列 $seq 的 GC 含量
- meanqual($seq) 返回 fastq 格式的序列 $seq 的平均Q值
- reverse($seq) 返回序列 $seq 的反意链
- revcomp($seq) 返回序列 $seq 的反意互补链
- qualcount($qual, threshold) Returns the number of quality values above the threshold parameter.
- trimq(qual, beg, end, param) Trims the quality string qual in the Sanger scale using Richard Motts algorithm (used in Phred). The 0-based beginning and ending positions are written back to beg and end, respectively. The last argument param is the single parameter used in the algorithm, which is optional and defaults 0.05.
- and(x, y) bit AND operation (& in C)
- or(x, y) bit OR operation (\| in C)
- xor(x, y) bit XOR operation (^ in C)
4.1 安装 bioawk¶
bioawk 的源代码托管在 github 上，如先用 git clone 命令先将代码下载到本地，然后编译生成 bioawk 将其添加到系统路径。
~$ sudo apt-get install bison
~$ git clone https://github.com/lh3/bioawk
~$ cd bioawk && make
~$ sudo cp bioawk /usr/local/sbin
4.2 使用 bioawk¶
文件格式：
- bed: * chrom * start - end - name - score - strand - thickstart - thickend - rgb - blockcount - blocksizes - blockstarts
- sam: - qname - flag - rname - pos - mapq - cigar - rnext - pnext - tlen - seq - qual
- vcf:
4.3 bioawk 应用示例¶
构建测序数据分析 pipeline 时，当 fastq 数据在做完 trimming 后，我们往往要关注剩下多少 reads，可以用 bioawk 进行快速统计。
# 快速统计fastq里的reads数量
~$ bioawk -c fastx 'END { print NR }' my_fastq.tar.gz
在一些特殊场合里，需要分析 reads 系列，用 bioawk 可以很方便快速来完成。比如统计特殊碱基开头的 reads 数。
# 统计以GATTAC开头的reads的数量
~$ bioawk -c fastx '$seq~/^GATTAC/ {++n} END { print n }' my_fastq.tar.gz
如果想看GC含量大于60%的reads的质量情况，可以用下面方法：
~$ bioawk -c fastx '{if (gc($seq)>0.6) printf "%s\n%s\n\n", $seq, $qual}' my_fastq.tar.gz
生成reads的平均Q值
~$ bioawk -c fastx '{print $name, meanqual($seq)}' my_fastq.tar.gz > meanqual.txt
Reference¶
Sed 使用教程¶
R 语言入门¶
安装R¶
Ubuntu 下的安装添加 ppa 源来安装最新的 R 版本。
$ sudo add-apt-repository ppa:marutter/rrutter
$ sudo apt-get update
$ sudo apt-get upgrade
$ sudo apt-get install r-base
安装 bioconductor¶
安装可以在 R 命令行界面里操作。在终端直接输入`R`，如果安装了Rstudio，也可以运行 rstudio。
> source("https://bioconductor.org/biocLite.R")
> biocLite()
R 在线教程¶
- OnePage R <http://togaware.com/onepager/>
docker 的世界¶
$ sudo apt install docker.io
git¶
NGS 基础¶
网络资源¶
将网络上一些关于微生物NGS相关的学习资源和工具网站纪录在此，方便查找与整理。注：部分网站需要翻墙。
免费课程资料¶
国外院校举办的生物信息学与高通量测序技术培训所提供的授课资料。
教学与资源¶
关于生物信息学与高通量测序技术相关内容的博客和网站。
最近域名被墙了，暂时用这个 [地址 <http://blog.qiuworld.com:8080/) 访问： * YGC * Biostack * 宠辱不惊 一心问学 * 每日一生信 * 佛说有缘 * Bacterial Genomics * Genome Unzipped * Titus’s blog * Blue Collar Bioinformatics * Mass Genomis * Next Genetics * Getting Genetics Done * Omics! * Pathogens genes & genomes * Opiniomics * YOKOFAKUN * In between lines of code * Core-genome * NGS wiki book][]* * `omictools * omicsmaps
在线工具和数据库¶
- PATRIC 非常棒的病原菌资源整合类数据库，用户界面设计比较友好，各种异步读取速度很快用户体验不错。
- ARDB 马里兰大学维护的耐药基因数据库，可惜2009年已经停止更新。
- CARD 加拿大麦克马斯特大学维护的病原微生耐药数据库。
- VFDB 中国医学科学院病原系统生物学研究室创建的病原菌毒力基因数据库
- MVIRDB Lawrence Livermore National Laboratory 创建毒力基因整合数据库。数据来源于包括VFDB在内的多个数据库。
- DBETH 对人致病的细菌外毒素数据库，由印度科研机构维护。
- `Food Microbe Tracker <http://www.foodmicrobetracker.com/login/login.aspx>`__S Cornell大学主持的食源性致病菌的数据库项目，主要有表型和PFGE分型资料。最大的特点可以提交自己的PFGE数据与数据库中的比对，但可能我使用不当或理解不对，比对出来的结果总是完全对不上，不知道它具体的算法
- Genomespace 基于web界面的基因组分析工具
- MBGD 日本人做的微生物比较基因组数据库
讲座课件¶
在线的各种PPT讲座课件，这里选择一些基础简介的可见和一些与微生物测序内容相关的课件
视频资料¶
如果嫌光看课件不够明确，Youtube里有一些不错的测序公司与基因组学视频频道可以收看
网络学校¶
希望学习一些基础知识的可以通过在线学校资源，里面有一些关于生物信息学的基础课程给没有学过的微生物学研究者入门。还有一些 Perl 或者 Python 脚本语言入门课程，给希望能编写脚本的工作者学习。另外还有一些如进化树分析的理论知识课程也可以参考。
生物信息数据分析¶
- DNAnexus DNAnexus基于云计算的生物信息分析平台
- Seven Bridges
- Tute Genomics
其他¶
- Biocompare 给购买者提供产品信息和技术资料的网站，对于实验室选择仪器和试剂有所帮助。
- Publons 发表的科研文章的reviewer发布reviews的网站，作者可与reviewer交互讨论
- genohub 提供测序技术及相关服务供应商的搜索信息
- AMD 美国CDC的病原微生物分子检测
- PHGKB of AMD 美国CDC的病原微生物知识库中的AMD相关内容
- Clinicl Informatrics News 临床信息学分析新闻综合
- NGS wiki book
第零章. 安装和使用软件¶
安装软件
Aspera connect¶
Aspera 是 IBM 研发的加速在线数据传输的一系列服务器端与客户端的应用。Aspera 的相关工具可以在 这里 下载获得。对于服务器用户，可以使用 Aspera Connent，目前的最新版本是 3.6.1（基于的 ascp 版本是3.5.6）。Aspera
Connect
是 Aspera
其中的一款用于浏览器下载的高效插件，因为其内嵌了ascp
命令可以用于服务器，而选择安装这个工具不是使用其 asperaconnect
工具。不过安装了 Aspera
Connect
的带图形界面的客户端，浏览器会在有 Connect, Faspex 或者 Shares
页面的链接处添加 ascp 下载快捷按钮。
1. 下载安装¶
Aspera
Connect
就会安装到当前用户 ~/.aspera/connect 目录或安装在
/opt/aspera/connect 目录下。
~$ cd /tmp
/tmp$ curl -O http://download.asperasoft.com/download/sw/connect/3.6.1/aspera-connect-3.6.1.110647-linux-64.tar.gz
/tmp$ tar zxf asper-commect-3.6.1.110647-linux.tar.gz
/tmp$ ./aspera-connect-3.6.1.110647-linux-64.sh*
对于喜欢使用客户端的桌面用户，也可以使用 Aspera Desktop Client 客户端程序来下载数据。
2. 设置¶
服务器端没有 GUI 图形化界面。如果在 Desktop 版本的 Linux 操作系统上，可以运行 asperaconnect, 可以看到相应的配置界面。
FastQC¶
SPAdes¶
1. 下载并安装 SPAdes¶
$ curl -O http://spades.bioinf.spbau.ru/release3.1.1/SPAdes-3.9.1-Linux.tar.gz
$ tar zxvf SPAdes-3.9.1-Linux.tar.gz
$ sudo mv SPAdes-3.9.1-Linux /opt/spades
$ sudo ln -s /opt/spades/bin/* /usr/local/sbin
2. 拼装基因组¶
# 对PE250以上的测序数据
$ spades.py -t 40 -k 21,33,55,77,99,127 --careful \
> -1 sample1_R1.fastq -2 sample2_R2.fastq -o assembly_spades
RAxML¶
1. 安装 RAxML¶
从 github 仓库克隆代码安装。
$ git clone https://github.com/stamatak/standard-RAxML
$ cd standard-RAxML
$ make -f Makefile.SSE3.PTHREADS.gcc && rm *.o
$ sudo ln -s `pwd`/raxmlHPC-PTHREADS-SSE3 /usr/local/bin/raxml
2. 使用 RAxML¶
# 绘制无根树
$ raxml -s alignment.phy -m GTRGAMMA -p 12345 -T 40 -w raxml_output -n raxml
# 绘制有跟树，如根名叫 xyz
$ raxml -s alignment.phy -m GTRGAMMA -p 12345 -o xyz
第一章. 基因组数据处理¶
本章主要介绍如何从公共数据库获得基因组 NGS 数据，并以公共数据库数据为例进行基因组数据前期处理，包括质控，拼接，注释等基本工作。
第1节. 公共数据库获取数据¶
完成测序后一般需要将测序原始数据如 reads 数据上传到公共数据库，以便在发表的文章中引用原始的测序结果。很多时候，我们也需要获取他人发表的公开测序数据，来帮助自己的研究领域。本节我们就来了解一下如何上传自己的测序数据和下载别人的测序数据。首先让我们来了解一下网络上的测序公共数据库：
最常用的测序公共数据库：
SRA 是 NCBI 为了应对越来越多的高通量测序数据而在 2007 年底推出的测序数据库，用于存储、显示、提取和分析高通量测序数据。而 ENA 则是由 EBI 负责维护的功能类似的数据库，同时作为 Ensembl、UniProt 和 ArrayExpress 等服务的底层基础。2者在主要功能方面非常类似，同时数据互通。
SRA 数据库简介¶
注解
SRA 是 Sequence Read Archive 的首字母缩写。SRA 与 Trace 最大的区别是将实验数据与 metadata（元数据）分离。metadata 是指与测序实验及其实验样品相关的数据，如实验目的、实验设计、测序平台、样本数据(物种，菌株，个体表型等)。metadata可以分为以下几类：
- Study：accession number 以 DRP，SRP，ERP 开头，表示的是一个特定目的的研究课题，可以包含多个研究机构和研究类型等。study 包含了项目的所有 metadata，并有一个 NCBI 和 EBI 共同承认的项目编号（universal project id），一个 study 可以包含多个实验（experiment）。
- Sample：accession number以 DRS，SRS，ERS 开头，表示的是样品信息。样本信息可以包括物种信息、菌株(品系) 信息、家系信息、表型数据、临床数据,组织类型等。可以通过 Trace 来查询。
- Experiment：accession number 以 DRX，SRX，ERX 开头。表示一个实验记载的实验设计（Design），实验平台（Platform）和结果处理 （processing）三部分信息。实验是 SRA 数据库的最基本单元，一个实验信息可以同时包含多个结果集（run）。
- Run：accession number 以DRR，SRR，ERR 开头。一个 Run 包括测序序列及质量数据。
- Submission：一个 study 的数据，可以分多次递交至 SRA 数据库。比如在一个项目启动前期，就可以把 study，experiment 的数据递交上去，随着项目的进展，逐批递交 run 数据。study 等同于项目，submission 等同于批次的概念。
下面以 SRA 数据库为例介绍数据下载和上传的操作方法。
下载数据¶
从 SRA
数据库下载数据有多种方法。可以用ascp快速的来下载 sra
文件，也可以用wget或curl等传统命令从 FTP 服务器上下载
sra 文件。另外 sratoolkit 工具集也支持直接下载。
1. Aspera¶
注意新版的ascp用.openssh作为密钥文件而不是原来的.putty
~/data$ ascp -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh \
> --user=anonftp --host=ftp.ncbi.nlm.nih.gov --mode=recv -l 100m -pQTk1 \
> /sra/sra-instant/reads/ByRun/sra/SRR/SRR955/SRR955386/SRR955386.sra .
2. sratoolkit¶
NCBI 开发了 sratoolkit 工具来帮助处理 SRA 数据，正确配置后可以很方便的下载 SRA 数据。
可以直接从 NCBI 上下载。最新的源码可以在 Github 获得。
# download from NCBI
~$ cd /tmp
/tmp$ wget http://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/2.5.6/sratoolkit.2.5.6-ubuntu64.tar.gz
/tmp$ tar zxf sratoolkit.2.5.6-ubuntu64.tar.gz -C ~/app
如果你安装并设置了 Aspera
Connect，那么prefetch会优先使用ascp方式来下载，如果没有安装或则ascp下载失败，则切换成
HTTP 方式下载 sra
数据。另外fastq-dump命令也能从远端直接下载数据，加上-X 1参数，会预先下载最前的5个
reads，加上-Z参数，则会将这些 reads 打印到终端输出。
# 下载 SRR955386.sra 到 你安装 sratoolkit 时配置的 public 目录中，默认是 ~/ncbi/public/sra
~$ prefetch -v SRR955386
# 下载 SRR955386.sra 到 你安装 sratoolkit 时配置的 public 目录中，默认是 ~/ncbi/public/sra，并且在终端输出5行 reads 数据。
~$ fastq-dump -X 5 -Z SRR955386
$ cat sra.list | xargs prefetch -v
获得了 .sra 文件后，需要将其转换成 .fastq
格式的文件，用fastq-dump可以很方便的实现。转换之前要注意的是该
run 的 metadata 里，测序类型是 SE 还是 PE 的。
# 将 sra 文件移动到 ~/data 目录中
~$ mv ~/.ncbi/public/sra/SRR955386.sra ~/data
# 如果是 SE 测序数据，直接运行以下命令
~$ fastq-dump SRR955386.sra
# 如果是 PE 测序数据，则要添加参数：--split-files
~$ fastq-dump --split-files SRR955386.sra
正常的 sra 文件的 metadata 应该包含测序采用的是 SE 还是 PE
的方式。但如果你不知道所下载的到底是 SE 还是 PE
格式的文件可以用fastq-dump -X 1 --split-spot的方法来判断。
# it's SE if nreads=1, and PE when nreads=2，统计整个文件，因此速度比较慢
~$ sra-stat -xs SRR955386.sra | grep "nreads"
# 如果输出是4，那么就是SE，如果是8,则是PE
~$ fastq-dump -X 1 --split-spot -Z SRR955386.sra | wc -l
# 或者加上参数让 fastq-dump 自己判断
# 当 sra 文件是 SE 测序时，fastq-dump只能dump出1个 *_1.fastq 文件
~$ fastq-dump --split-files ERR493452.sra
~$ mv ERR493452_1.fastq ERR493452.fastq
当需要判断批量下载的 sra 文件时，区分那些是 PE 的那些是 SE 的文件，可以用以下脚本：
import os
import subprocess
def check_SRA_type(fn):
fn = os.path.abspath(fn);
try:
contents = subprocess.check_output(["fastq-dump", "-X", "1", "-Z", "--split-spot", fn]);
except subprocess.CalledProcessError, e:
raise Exception("Error running fastq-dump on", fn);
# -X 1 will output 4 lines if SE, and 8 lines if PE
if(contents.count("\n") == 4):
return False;
elif(contents.count("\n") == 8):
return True:
else:
raise Exception("Unexpected output from fastq-dump on ", filename);
如果想下载一个完整的 project 数据，可以利用 entrezdirect 工具。
~/tmp$ wget ftp://ftp.ncbi.nlm.nih.gov/entrez/entrezdirect/edirect.tar.gz
~/tmp$ tar zxvf edirect.tar.gz -C ~/app
~/tmp$ sudo ln -s ~/app/edirect/* /usr/local/sbin/
~/data$ esearch -db sra -query PRJNA40075 | efetch --format runinfo | cut -d ',' -f 1 | grep SRR | xargs fastq-dump --split-files
如果想下载不同 research 的数据，可以自己建立一个 accession number list
的文件，比如利用 NCBI 的 entrez
网页界面，导出所需要的数据 accession number，然后利用 ascp
下载。不过建议ascp下载不要太狠心，否则容易被
NCBI 给封掉。
上传数据¶
SRA 上传测序结果可以参照 NCBI文档 来实现。一般上传数据到NCBI SRA的过程需要6步：
- Create a BioProject for this research
- Create a BioSample submission for your biological sample(s)
- Gather Sequence Data Files
- Enter Metadata on SRA website 4.1 Create SRA submission 4.2 Create Experiment(s) and link to BioProject and BioSample 4.3 Create Run(s)
- Transfer Data files to SRA
- Update Submission with PubMed links, Release Date, or Metadata Changes
需要注意的一点是，上传的过程中很多地方一旦保存或提交就不可以修改，尤其是各处的Alias。但是，可以联系NCBI的工作人员修改内容。NCBI的工作效率是很高的，一般不超过48小时，就可以得到确认，并拿到登录号。
Reference¶
- NCBI上传数据文档
- 熊筱晶, NCBI高通量测序数据库SRA介绍, 生命的化学[J], 2010:6, 959-963.
- http://blog.sciencenet.cn/blog-656335-908140.html
- https://www.biostars.org/p/139422/
- https://www.youtube.com/watch?v=NSIkUHKRPpo
old stuff will be remove soon¶
SRA 数据处理¶
本节基本内容介绍：
- SRA 数据库简介
- SRA 数据下载
- SRA 数据转换
- 从基因组数据模拟短序列数据
SRA 是 Sequence Read Archive 的首字母缩写。
SRA 与 Trace 最大的区别是将实验数据与元数据分离。元数据现在可以划分为以下几类。
- Study：study 包含了项目的所有 metadata，并有一个 NCBI 和 EBI 共同承认的项目编号（universal project id），一个 study 可以包含多个实验（experiment）。
- Experiment：一个实验记载实验设计（Design），实验平台（Platform）和结果处理 （processing）三部分信息，并同时包含多个结果集（run）。
- Run：一个结果集包括一批测序数据。
- Submission：一个 study 的数据，可以分多次递交至 SRA 数据库。比如在一个项目启动前期，就可以把 study，experiment 的数据递交上去，随着项目的进展，逐批递交 run 数据。study 等同于项目，submission 等同于批 次的概念。
对于单个SRA数据，带图形界面的电脑可以通过浏览器下载。对于只有命令行的服务器端的操作，终端环境下无法使用GUI软件，所以要通过命令工具来下载SRA数据。
1. 用 edirect 下载数据
edirect 是 NCBI 最近发布的 entrez 数据操作命令行工具。过去往往要通过 biopython 之类的第三方脚本工具来实现查找，索引以及抓取 entrez 的数据。现在可以用 edirect 来很方便的实现。
~/tmp$ wget ftp://ftp.ncbi.nlm.nih.gov/entrez/entrezdirect/edirect.tar.gz
~/tmp$ tar zxvf edirect.tar.gz -C ~/app
~/tmp$ sudo ln -s ~/app/edirect/*
2. 用 sra toolkit 下载数据
选择美国 FDA 提交的 SRR955386 数据。该数据是在Illumina Miseq平台上PE250的数据。生物样本是一株分离自奶酪的单增李斯特菌 CFSAN006122。如果你已经安装好 NCBI 的 sra_toolkit 工具，可以直接用里面的 prefetch 下载。
~/data$ prefetch -v SRR955386
~/data$ mv ~/.ncbi/public/sra/SRR955386.sra .
Notes: 2.3.x版本的 sra_toolkit 会有些配置问题，建议使用最新版。不过2.4.x版本的程序调用ascp时，比如 prefetch 调用 ascp 时选择的 ssh 密钥是老版 ``asperaweb_id_dsa.putty``, 而目前新版的ascp在linux里是使用 ``asperaweb_id_dsa.openssh``。
3. 用aspera connect下载
用 aspera connect 工具下载 NCBI 的数据速度非常快，而且不光可以下载 SRA 数据库里的数据，还可以下载 NCBI FTP 里的其他资源。
~/data$ ascp -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh \
> --user=anonftp --host=ftp.ncbi.nlm.nih.gov --mode=recv -l100m -T -k 1 \
> /sra/sra-instant/reads/ByRun/sra/ERR/ERR175/ERR175655/ERR175655.sra .
4. 用 ascp 批量下载 sra 数据
从 NCBI SRA
数据库检索到需要的序列，比如搜索关键词Listeria monocytogenes miseq，左侧filters定义为DNA和Genome。
用 python 脚本批量转换路径。
# filename: format_path.y
with open('acc_list.txt', 'r') as f:
data = f.read().split('\n')
list=[]
for i in data:
path = '/sra/sra-instant/reads/ByRun/sra/%s/%s/%s/%s.sra\n' % (i[0:3], i[0:6], i, i)
list.append(path)
with open('acc_list_full.txt', 'w') as f:
d = f.writelines(list)
下载数据
~/data$ touch acc_list_full.txt
~/data$ python format_path.py
~/data$ mkdir -p sra
~/data$ ascp -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh --user=anonftp --host=ftp-private.ncbi.nlm.nih.gov --mode=recv --file-list acc_list_full.txt sra/
安装 sra toolkit (for ubuntu x64 version)
~/tmp$ curl -O http://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/2.4.3/sratoolkit.2.4.3-ubuntu64.tar.gz ~/tmp$ tar zxvf sratoolkit.2.4.3-ubuntu64.tar.gz -C ~/app ~/tmp$ cd ~/app/sratoolkit.2.4.3-ubuntu64 ~/tmp$ sudo ln -s `pwd`/bin/* /usr/local/sbin/
sra数据转成fastq格式
SRR955386 这个数据的样本还用 Pacbio SMRT 平台进行了测序，Pacbio 单分子测序技术获得基因组完成图。用该完成图作为模板，考量一下不同拼接软件的拼接结果。
将 CSFAN006122 的基因组完成图数据下载。
~/data$ prefetch -v SRR955386
~/data$ mv ~/.ncbi/public/sra/SRR955386.sra .
~/data$ fastq-dump --split-files SRR955386.sra
完成后可以看到 data 目录下新增了2个文件 SRR955386_1.fastq 和
SRR955386_2.fastq
前面介绍如何从 NCBI SRA 数据库下载 NGS 数据并转换成 fastq 格式的文件。可能有些人做的数据分析或验证需要从基因组逆向模拟成短序列格式数据，这种需求也有人开发了相的软件。
wgsim 就是这样一个软件，它是由开发了BWA等软件大牛 Li heng 写的基因组转成短序列的软件。
安装wgsim
~/app$ git clone https://github.com/lh3/wgsim.git && cd wgsim
~/app$ gcc -g -O2 -Wall -o wgsim wgsim.c -lz -lm
~/app$ sudo ln -s `pwd`/wgsim /usr/local/sbin
~/app$ wsgim -h
Program: wgsim (short read simulator)
Version: 0.3.1-r13
Contact: Heng Li <lh3@sanger.ac.uk>
Usage: wgsim [options] <in.ref.fa> <out.read1.fq> <out.read2.fq>
Options: -e FLOAT base error rate [0.020]
-d INT outer distance between the two ends [500]
-s INT standard deviation [50]
-N INT number of read pairs [1000000]
-1 INT length of the first read [70]
-2 INT length of the second read [70]
-r FLOAT rate of mutations [0.0010]
-R FLOAT fraction of indels [0.15]
-X FLOAT probability an indel is extended [0.30]
-S INT seed for random generator [-1]
-A FLOAT disgard if the fraction of ambiguous bases higher than FLOAT [0.05]
-h haplotype mode
模拟基因组短序列数据
使用所有参数为默认值，将大肠杆菌基因组数据转换为PE250 fastq 格式数据。
~/data$ wget https://raw.github.com/ecerami/samtools_primer/master/tutorial/genomes/NC_008253.fna
~/data$ wgsim -S11 -1250 -2250 NC_008253.fna reads_1.fastq reads_2.fastq
~/data$ head -8 reads_1.fastq
@gi|110640213|ref|NC_008253.1|_3151106_3151623_4:0:0_5:1:0_0/1
AACATAAGTGGTATTAATTCCCCAAGATTCAAGATTACGAATAGTATTATCCGCAAAAATATCATCACCTACTTTCGTCAGCATCAGGACTTTTGAATTCAATTTAGCCGCCGCCACCGCTTGATTAGCACCTTTGCCACCACATCCGATTTTGAAGGCAGGTGCTTCCAGAGTTTCTCCTTCTTTAGGCATCTGATAAGTGTAAGTAATGAGATCCACCATATTGGAACCAATAAGTGCAATGTCCATT
+
2222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222
@gi|110640213|ref|NC_008253.1|_2429297_2429873_8:0:0_7:0:0_1/1
ACGGTATCACCATGATTGCCCGTTCCGTCAACAGCATGGGGCTGGTCATTATAGGCGGTGGATCGCTTGAAGAAGCGTTAACTGAACTGGAAACCGGACGCGGCGACGCGGTGGTGGTGCTGGAAAACGAACTGCATCGTCACGCTTGCGCTACCCGCGTGAATGCTGCGCTGGCTAAAGCGCCGCTGGTGATTGTGGTTGACCATCAACGCACAGCGATTATGGAAAACGCTCATCTGGTACTATCCGC
+
2222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222
- SRA数据库帮助
- SRA Handbook
- 熊筱晶, NCBI高通量测序数据库SRA介绍, 生命的化学[J]: 2010, 06, 959-962.
第1节. 数据质控¶
查看测序数据的质量情况。
1. FastQC¶
FastQC 可能是最常用，也是最直观的测序数据质量控制软件。
1.1 安装 FastQC¶
FastQC 可以在图形界面下运行，如果在不带有图形界面的系统环境里也可以通过生成 html 报告再通过浏览器打开来查看结果。
~/tmp$ wget http://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.5.zip
~/tmp$ unzip fastqc_v0.11.5.zip -d ~/app
~/tmp$ sudo ln -s ~/app/FastQC/fastqc /usr/local/sbin
# run fastqc in GUI mode
~$ fastqc
# run fastqc in command line mode
~$ fastqc -f fastq ~/data/*.fastq
1.2 生成报告¶
生成 html 报告文件和对应的 zip 压缩文件，并通过 scp 命令传输到本地后用浏览器打开查看。
# scp your linux/mac desktop system
~$ scp -i username@server-ip:~/data/* ~/
# open it by google-chrome
~$ google-chrome *.html
1.3 结果说明¶
FastQC 结果由11个模块组成，对于结果报告各个模块的说明可以参见 FastQC 文档。
| Module | Explanation |
|---|---|
| Basic Statistics | fastq文件的基本信息，可以看到序列数量和 读长，fastq文件版本，GC含量等。 |
| Per base sequence quality | 最重要的结果模块之一，你可以从这个模块 的图示中了解序列中各个碱基的质量。 |
| Per sequence quality score | ... |
| Per base sequence content | ... |
数据处理：通过FastQC各个模块查看，发现低质量的序列需要切除或者序列污染的，比如smallRNA测序，文库被污染了接头等情况时，需要对数据进行去除污染序列的处理。
2. Picard¶
Picard 是 BroadInstitute 使用 Java 语言开发的针对 BAM 等高通量测序数据的处理工具，特别是针对 Illumina 平台的数据。它可以提供一个测序质量和数据基本特性的报告。
2.1 安装 Picard¶
~/tmp$ wget https://github.com/broadinstitute/picard/releases/download/1.123/picard-tools-1.123.zip
~/tmp$ unzip -n picard-tools-1.123.zip -d ~/app
2.2 使用 Picard¶
Reference¶
第2节. 基因组拼接¶
微生物基因组测序后的拼接工作，有许多软件可以完成。近几年不断有新的优秀拼接软件涌现，本节仅挑选几个微生物基因组拼接使用较多也是评测中比较优秀的几个工具。也可以先对自己的测序物种做一个 assembly evaluation，选择比较适合的拼接软件。
常用拼接软件工具列表 - SPAdes - MaSuRCA - Velvet - A5_miseq
更多工具可以参见`这里 <https://omictools.com/genome-assembly-category>`_
1. SPAdes¶
SPAdes 是由
1.1 下载并安装 SPAdes¶
~/tmp$ curl -O http://spades.bioinf.spbau.ru/release3.9.0/SPAdes-3.9.0-Linux.tar.gz
~/tmp$ tar zxf SPAdes-3.9.0-Linux.tar.gz -C ~/app
~/tmp$ cd ~/app
~/app$ sudo ln -s `pwd`/SPAdes-3.9.0-Linux/bin/* /usr/local/sbin
1.2 拼装基因组¶
# -t 40 表示用40个 threads 进行拼接，根据自己电脑的CPU情况自行设置。
~/app$ spades.py -t 40 -1 SRR95386_1.fastq -2 SRR955386_2.fastq -o spades_output
SPAdes会尝试不同的Kmer，因此拼装时间也会根据Kmer选择数量成倍增加。对于常见的 Miseq v2/v3 试剂盒，采用 PE150/PE250/PE300 的读长测序，常用的拼接命令是：
# PE150 读长测序数据
~/apps$ spades.py -t 40 -k 21,33,55,77 --careful -1 SRR95386_1.fastq -2 SRR955386_2.fastq -o SRR95386_output
# PE250/300 读长测序数据
~/apps$ spades.py -t 40 -k 21,33,55,77,99,127 --careful -1 SRR95386_1.fastq -2 SRR955386_2.fastq -o SRR95386_output
2. Macursa¶
2.1 下载并安装 MaSuRCA¶
~/tmp$ curl -O ftp://ftp.genome.umd.edu/pub/MaSuRCA/MaSuRCA-3.1.3.tar.gz
~/tmp$ tar zxvf MaSuRCA-3.1.3.tar.gz -C ~/app
~/tmp$ cd ~/app
~/app$ ./MaSuRCA-2.3.0/install.sh
~/app$ sudo ln -s `pwd`/MaSuRCA-2.3.0/bin/masurca /usr/local/sbin
2.2 建立拼装配置文件¶
~/data$ touch SRR955386_masurca.config
~/data$ nano SRR955386_masurca.config
文件内容修改如下：
DATA
PE= p1 500 50 SRR955386_1.fastq SRR955386_2.fastq
END
PARAMETERS
GRAPH_KMER_SIZE=auto
USE_LINKING_MATES=1
LIMIT_JUMP_COVERAGE = 60
KMER_COUNT_THRESHOLD = 1
NUM_THREADS= 2
JF_SIZE=100000000
DO_HOMOPOLYMER_TRIM=0
END
设置 GRAPH_KMER_SIZE=auto，软件会调用Kmer=31来进行拼装。对于MiSeq PE250以上的插片，可以考虑手动设置使用更大的Kmer。
2.3 开始拼装¶
~/app$ masurca SRR955386_masurca.config
~/app$ ./assemble.sh
3. Velvet¶
3.1 下载并安装¶
Velvet 是一个老牌的基因组测序数据拼接软件。Velvet最新版本是1.2.10，可以访问 官方网站 下载源代码包，也可以通过克隆 软件仓库 的方式获得最新的源代码。
下载源代码
普通用户可以下载源代码包自行编译获得软件。
~$ cd /tmp
~/tmp$ wget https://www.ebi.ac.uk/~zerbino/velvet/velvet_1.2.10.tgz
~/tmp$ tar xvf velvet_1.2.10.tgz -C ~/app/velvet
克隆代码仓库
如果在软件运行中遇到问题，想试用最新版代码，或是有能力提交issues，或者想改进软件参与开源代码编写的可以选择克隆代码库的方式。
~$ cd ~/tmp
~/tmp$ git clone https://github.com/dzerbino/velvet.git
编译安装
make可以进行编译，有几个编译参数可以选择，分别是MAXKMERLENGTH, DIRECTORY
~$ cd ~/app/velvet
~/app/velvet$ make 'MAXKMERLENGTH=127'
~/app/velvet$ sudo cp velveth velvetg /usr/local/sbin
3.2 拼接基因组¶
velvet软件由2个可执行文件 velveth 和 velvetg 组成。
3.3 其他组件¶
4. A5-miseq¶
A5-miseq 是一个用 perl 开发的针对细菌基因组 de novo assembly 的 pipeline 工具。它本身不参与组装，而是通过组合一套工具来完成工作，工具集包括：
- bwa
- samtools
- SGA
- bowtie
- Trimmomatic
- IDBA-UD
- SSPACE
这些工具都以及集成在 A5-miseq 中，不需要另外安装。为了避免不同版本的 samtools，bowtie 对结果产生的差异，建议采用虚拟环境如 docker 等来隔离运行环境。
4.1 安装 A5-miseq¶
~/tmp$ wget http://downloads.sourceforge.net/project/ngopt/a5_miseq_linux_20150522.tar.gz
~/tmp$ tar zxvf a5_miseq_linux_20150522.tar.gz -C ~/app
~/tmp$ sudo ln -s ~/app/a5_miseq_linux_20150522/bin/a5_pipeline.pl /usr/local/sbin
~/docker$ mkdir -p a5-miseq && cd a5-miseq
~/docker/a5-miseq$ touch Dockerfile
FROM ubuntu:latest
MAINTAINER Mark Renton <indexofire@gmail.com>
RUN apt-get update -qy
RUN apt-get install -qy openjdk-7-jre-headless file
ADD http://downloads.sourceforge.net/project/ngopt/a5_miseq_linux_20150522.tar.gz /tmp/a5_miseq.tar.gz
RUN mkdir /tmp/a5_miseq
RUN tar xzf /tmp/a5_miseq.tar.gz --directory /tmp/a5_miseq --strip-components=1
ADD run /usr/local/bin/
ADD Procfile /
ENTRYPOINT ["/usr/local/bin/run"]
4.2 使用 A5-miseq¶
$ perl a5_pipeline.pl SRR955386_1.fastq SRR955386_2.fastq ~/data/a5_output
4.3 A5-miseq 文档¶
查看 A5-miseq 工具的使用文档可以用 a5_pipeline.pl 工具查看。
# Usage:
$ a5_pipeline.pl [--begin=1-5] [--end=1-5] [--preprocessed] <lib_file> <out_base>
第3节. 注释拼接基因组¶
基因组草图获得后，对其进行注释。
1. Prokka¶
1.1 安装 Prokka¶
~$ sudo apt-get install cpanminus
~$ cpanm Bio::Perl XML::Simple
这里选择Barrnap
~/tmp$ wget http://www.vicbioinformatics.com/barrnap-0.5.tar.gz
~/tmp$ tar zxvf barrnap-0.5.tar.gz -C ~/app/
~/tmp$ echo "export PATH=$PATH:$HOME/app/barrnap-0.5/bin" >> ~/.bashrc
~/tmp$ source ~/.bashrc
~/tmp$ wget http://www.vicbioinformatics.com/prokka-1.11.tar.gz
~/tmp$ tar zxvf prokka-1.11.tar.gz -C ~/app
~/tmp$ echo "export PATH=$PATH:$HOME/app/prokka-1.11/bin" >> ~/.bashrc
~/tmp$ source ~/.bashrc
~/tmp$ prokka --setupdb
1.2 注释基因组¶
~/data$ prokka genome_contig.fasta
2. PGAAP¶
http://www.ncbi.nlm.nih.gov/genome/annotation_prok/ http://www.ncbi.nlm.nih.gov/books/NBK174280/
3. RAST¶
Reference:¶
1. 2.
Miseq 下机数据分析应用实例¶
在 Miseq 上通过对24个混合样本的PE测序，获得的下机 PF Data，从 Quality Control 到 de novo assembly 过程的基本操作流程。[#f1]
1. 前期处理¶
1.1 分样¶
为了避免文件混乱，首先将获得的001~024样本的 fastq 文件分别移动到以样本名命名的 24个文件夹中，分别进行操作。
# 产生的数据文件名为格式为 0*_S*_L001_R1_001.fastq.gz/0*_S*_L001_R2.fastq.gz
~/data$ ls
001_S1_L001_R1_001.fastq.gz 003_S3_L001_R1_001.fastq.gz 005_S5_L001_R1_001.fastq.gz
001_S1_L001_R2_001.fastq.gz 003_S3_L001_R2_001.fastq.gz 005_S5_L001_R2_001.fastq.gz
002_S2_L001_R1_001.fastq.gz 004_S4_L001_R1_001.fastq.gz 006_S6_L001_R1_001.fastq.gz
002_S2_L001_R2_001.fastq.gz 004_S4_L001_R2_001.fastq.gz 006_S6_L001_R2_001.fastq.gz
...
022_S22_L001_R1_001.fastq.gz 023_S23_L001_R1_001.fastq.gz 024_S24_L001_R1_001.fastq.gz
022_S22_L001_R2_001.fastq.gz 023_S23_L001_R2_001.fastq.gz 024_S24_L001_R2_001.fastq.gz
# 新建以样本名为文件夹名，并移动数据到文件夹的raw目录下
~/data$ for i in $(awk -F"L001" '{gsub("_$","",$1);print $1}' \
> <(ls -D *.fastq.gz) | sort | uniq); \
> do mkdir $(basename $i) && mkdir $(basename $i)/raw; \
> mv $i*.fastq.gz $(basename $i)/raw/ ; done
# 查看生成的各个样本的文件夹
~/data$ ls -d */
001_S1 002_S2 003_S32 004_S4 005_S5 006_S6
...
1.2 查看测序质量¶
用 FastQC 生成 html 格式的质量结果报告，用 Python 自带模块 SimpleHTTPServer 建立一个简便的 HTTP Server，以便查看 html 文档。
~/data$ mkdir -p qc
~/data$ for i in $(ls -d 0*/raw/*.fastq.gz); \
> do fastqc $i --extract -t 40 -q -o qc ; done
~/data$ python -m SimpleHTTPServer
然后用客户端浏览器访问服务器IP:8080来查看，比如我们的服务器IP地址是 10.44.35.122，就在浏览器地址栏里输入http://10.44.35.122:8000，可以看到所在文件 夹的文件链接页面了。
用Quast生成
1.3 去除接头¶
从前面的数据显示，部分插入片段比较短，因此部分reads测出接头序列。
1.4 去除低质量reads¶
Miseq v3 的PE300试剂最大的问题在于3’端急速下降的质量，特别是PE测序中R2端的质量。
Reference
| [1] | ../chapter_04/index.html |
第二章. 序列比对¶
第1节. 使用本地 Blast 应用¶
blast 工具用来进行序列相似性比对。常见的使用是通过访问 http://blast.ncbi.nlm.nih.gov/Blast.cgi，使用基于 Web 及面的 blast 工具，调用 NCBI 相应的序列数据库来进行比对。此外在 NCBI 很多工具集里。
注解
blast 工具可以进行序列相似性比对。在 NGS 数据分析中经常会被使用到，特别是一些工具中需要使用 blast 来作序列比对。拿到测序完成的草图后，因为基因组数据较大，连接NCBI网站往往又非常缓慢，所以要做 Blast 比对的话都需要做 Local Blast。这里介绍2个方式来实现：
- 用命令行进行 Blast
- 快速构建 Blast Web 服务
1. 使用本地 Blast¶
blast 和 ncbi_blast+ 这2个程序集容易搞混。NCBI 最早在1989年创建 Basic Local Alignment Search Tool 工具，沿用至2009年无论是命令行工具或是在线程序，都称呼其为blast。2009年 NCBI 鉴于 blast 的一些不足，重新开发了新的 ncbi_blast+ 命令行工具，新的 ncbi_blast+ 工具在速度上有了提升，在输入输出上也更为灵活。
目前 Blast 工具最新版是2012年发布的2.2.26（已经暂停版本更新并不做进一步支持），而 ncbi_blast+ 目前一直在更新。要区分2者也很简单，blast 是通过 blastall -p 的方式调用子程序来比对搜索的，而 blast+ 则是直接使用 blastn 或 blastp 等子程序来比对搜索。另外前者用 formatdb 程序来格式化数据库，后者用 makeblastdb 程序来格式化数据。
1.1 Blast¶
Blast的下载与安装
# 安装 Ubuntu 编译包。
~$ sudo apt-get install blast2
# 或者直接下载NCBI Linux预编译包，并解压缩安装
~$ wget ftp://ftp.ncbi.nih.gov/blast/executables/legacy/2.2.26/blast-2.2.26-x64-linux.tar.gz
~$ tar zxvf blast-2.2.26-x64 -C ~/app
~$ sudo ln -s ~/app/blast-2.2.26-x64/blastall /usr/local/sbin
运行 Blast
Blast 通过调用 blastall 这个程序来调用不同算法和程序实现序列比对。在命令行中输入 blastall ，会打印出一份参数列表。你也可以使用 man blast 来查看 Blast 工具的用户手册。
~$ blastall
参数说明：
-p: 指定要运行的 blast 程序。可以调用的有:
blastp 将氨基酸序列与蛋白质数据库比对
blastn 将核酸序列与核酸数据库比对
blastx 将核酸序列翻译成蛋白序列后与非冗余(nr)蛋白质数据库比对
psitblastn
tblastn
tblastx
-d: 指定要调用数据库，默认值是非冗余(nr)数据库。本地 Blast 比对一般来说都是调用 formatdb 格式化的数据库。
-i: 输入，默认值是终端输入，也可以使用文件的方式，比如 -i seq.fasta
-o: 输出，默认值是终端打印输出，也可以使用文件的方式，比如 -o result.txt
-e: 期望值。
-T: 输出文件格式，默认值为F(False)，想输出HTML格式的可以用 -T T
-M: 矩阵算法选择，默认值是 BLOSUM62
-n: 如果想使用MegaBlast，就设置 -n T
-b: 数据库中的比对结果显示条目数，默认是250条记录。
-m: 比对结果的输出显示方式，值为0~11，默认是0。各个数字含义见下表。
| 数值 | 代表含义 |
|---|---|
| 0 | pairwise |
| 1 | query-anchored showing identities |
| 2 | query-anchored no identities |
| 3 | flat query-anchored, show identities |
| 4 | flat query-anchored, no identities |
| 5 | query-anchored no identities and blunt ends |
| 6 | flat query-anchored, no identities and blunt ends |
| 7 | XML Blast output |
| 8 | tabular |
| 9 | tabular with comment lines |
| 10 | ASN, text |
| 11 | ASN, binary [Integer] |
其他参数参见man blast。
应用举例:
Local blast例子：获得一条大肠杆菌序列 myseq.fasta，要和 EDL933 进行比对。首先下载 EDL933 的基因组数据，并格式化作为本地数据库，然后使用 blastn 对序列进行比对。
# 使用 edirect 工具的 efetch 下载 EDL933 基因组数据
~$ efetch db=nuccore -format=fasta -id=AE005174 > AE005174.fasta
# 使用 blast 的 formatdb 工具将 EDL933 基因组数据格式化成用于比对的数据库格式
~$ formatdb -i AE005174.2.fasta -o T -p F
# 调用 blastn 的方式比对 myseq.fasta 和 EDL933 序列
~$ blastall -i myseq.fasta -d AE005174.2.fasta -p blastn
1.2 NCBI_blast+¶
目前 NCBI_blast+ 最新版为 v2.4.0。
下载安装 NCBI Blast+
# 安装 Ubuntu 编译包
~$ sudo apt-get install ncbi-blast+
# 直接在 NCBI 官方 FTP 站点下载预编译包解压缩安装
~/tmp$ wget ftp://ftp.ncbi.nih.gov/blast/executables/blast+/2.4.0/ncbi-blast-2.4.0+-x64-linux.tar.gz
~/tmp$ tar zxf ncbi-blast-2.4.0+-x64-linux.tar.gz
构建数据库
~$ makeblastdb -in data/database.fasta -dbtype nucl -parse_seqids
应用举例:
NGS 测序时为了保证 DNA 质量，往往会对物种进行鉴定，采用 16s rDNA 测序的方法。公司给了一堆拼接的txt文件，里面内容是不同物种样本的 16s rDNA 序列。当样本量很大，比如1000条时，如果一条条去 NCBI 上比对也是浪费精力和时间。可以用 local blast 来批量处理。首先下载各个目的物种的16s rDNA序列，并将其格式化成数据库。
# 生成结果报告
~$ for i in *.txt; do blastn -db database.fasta -query $i -outfmt 6 >> result; done
http://www.personal.psu.edu/iua1/courses/files/2014/lecture-12.pdf
2. Local Blast 使用 nr 数据库¶
2.1 下载 NCBI nr 数据库
~$ nohup update_blastdb.pl nt nr > log &
~$ blastn -query myseq.fa -db nt -outfmt 11 -out result
3. 构建自己的 blast web 服务¶
3.1 blastkit¶
blastkit 是一个包含webserver等工具的blast工具集。
安装依赖包
~$ sudo pip install pygr
~$ sudo pip install whoosh
~$ sudo pip install git+https://github.com/ctb/pygr-draw.git
~$ sudo pip install git+https://github.com/ged-lab/screed.git
~$ sudo apt-get -y install lighttpd
对lighttpd webserver进行配置
~$ cd /etc/lighttpd/conf-enabled
~$ sudo ln -fs ../conf-available/10-cgi.conf ./
~$ sudo echo 'cgi.assign = ( ".cgi" => "" )' >> 10-cgi.conf
~$ sudo echo 'index-file.names += ( "index.cgi" ) ' >> 10-cgi.conf
~$ sudo /etc/init.d/lighttpd restart
本地安装 Blast
~$ cd tmp
~/tmp$ wget ftp://ftp.ncbi.nih.gov/blast/executables/legacy/2.2.26/blast-2.2.26-x64-linux.tar.gz
~/tmp$ tar xzf blast-2.2.26-x64-linux.tar.gz -C ~/app
~/tmp$ sudo cp ~/app/blast-2.2.26/bin/* /usr/local/sbin
~/tmp$ sudo cp -r ~/app/blast-2.2.26/data /usr/local/blast-data
安装 blastkit
~/tmp$ cd ~/app
~/app$ git clone https://github.com/ctb/blastkit.git -b ec2
~/app$ cd blastkit/www
~/app/blastkit/www$ sudo ln -fs $PWD /var/www/blastkit
~/app/blastkit/www$ mkdir files
~/app/blastkit/www$ chmod a+rxwt files
~/app/blastkit/www$ chmod +x /root
~/app/blastkit/www$ cd ..
~/app/blastkit$ python ./check.py
如果安装顺利，就会提示一切已经准备完毕。接下来要准备数据。
~$ cp /mnt/assembly/ecoli.21/contigs.fa ~/app/blastkit/db/db.fa
~$ cd ~/app/blastkit
~$ formatdb -i db/db.fa -o T -p F
~$ python index-db.py db/db.fa
第2节. 多重序列比对¶
序列比对常用到的有 多重序列比对 MSA(Multiple Sequences Alignment)
2.1 ClustalW
2.2 Clustal Omega
2.3 Mafft
$ maftt seq/*.fasta > alignment
2.4 Muscle
寻找突变位点¶
案例分析: 铜绿假单胞菌临床株与环境株溯源¶
本节实践：
利用高通量测序技术对铜绿假单胞菌的临床分离株与环境分离株在基因组水平上进行SNP分析，从而溯源菌株并做耐药分析。
分析方法：
利用 bwa 将目标基因组 reads 比对到参考基因组，用 freebayes 做 variant calling，获得各个样本的 SNPs，根据结果做进一步分析。
实践意义：
可以作为今后实验室开展病原菌溯源工作的一个技术指导依据。
1. 提出问题¶
工作流程将从以下几个方面对于这个案例开展数据分析：
- 测序质量：reads的长度，插入片段的长度
- 将测序数据比对到参考基因组ST17
- 每个样品reads的分布情况
- 每个样品的基因组测序覆盖度
- 不同样品之间的差异，特别是当设置不同的过滤参数时的情况： - Quality filter: QUAL > 100 - Likelihood filter: QUAL / AO > 10 - Variant frequency filter: AO / DP > 0.9 (90% allele frequency)
- 哪一个水龙头水样的分离株最有可能是感染源。(方法:
用
vcfsamplediff比较环境株与临床株SNPs数量) - 2个病人分离株之间使用了大量的抗生素治疗，哪一个获得、得了耐药性？
- 用 joint calling 来比较2个临床样品数据，并用
vcfsamplediff分析 - 对于获得的大量SNPs，是否需要过滤以及如何过滤
- 了解SNPs的分布规律
- 重点查看SNPs
2. 分析流程¶
数据可以从这里下载获得，可以用wget -m递归下载路径下所有文件。文件包括：
- TAP1 - 病房1号水龙头水样中分离的菌株
- TAP2 - 病房2号水龙头水样中分离的菌株
- TAP3 - 病房3号水龙头水样中分离的菌株
- PATIENT1 - 服用抗生素后的病人分离株
- PATIENT2 - 一周内进一步服用抗生素的病人分离株
~$ mkdir -p pse_data && cd pse_data
~$ wget -m http://www.microbesng.uk/filedist/pseudomonas-practical/
~$ mv www.microbesng.uk/filedist/pseudomonas-practical/* .
用bwa mem进行比对，输出的结果通过管道由samtools转换成 BAM 文件并进行排序和索引。
~$ bwa mem -t 4 -R '@RG\tID:TAP1\tSM:TAP1' ST17.fasta \
> TAP1_R1_paired.fastq.gz TAP1_R2_paired.fastq.gz | \
> samtools view -bS - | samtools sort - TAP1.sorted
~$ samtools index TAP1.sorted.bam
参数介绍： \* -t 4: 表示使用4线程，根据自己的计算机CPU资源设置。 \*
-R '...': 设置reads的ID，:raw-latex:`\t为制表符会转义成tab键`。
对其他2株环境株和2株临床株重复以上步骤。
bam文件含有比对的信息，可以用stats参数查看
~$ samtools stats TAP1.sorted.bam | grep "maximum length"
~$ samtools stats TAP1.sorted.bam | grep "insert size"
~$ samtools stats TAP1.sorted.bam | grep "^COV" > TAP1.coverage.txt
将截取的覆盖度相关信息保存成.txt文件，在R里用ggplot2绘制覆盖度图：
> library(ggplot2)
> cov=read.table("TAP1.coverage.txt", sep="\t")
> ggplot(cov, aes(x=V3, y=V4)) + geom_bar(stat="identity") + xlab("Coverage") + ylab("Count")
用同样的方法对其他几个样本进行操作，获得sorted的bam格式文件，用freebayes对临床与环境株比较获得vcf文件。然后用vcflib工具过滤QUAL/AO>10的部分(VCF元数据部分可以参考文档)，并用wc-l计算行数统计结果。
~$ freebayes -p 1 -C 5 -f ST17.fasta TAP1.sorted.bam PATIENT1.sorted.bam > compare_tap1.vcf
~$ vcffilter -f "QUAL / AO > 10" compare_tap1.vcf | vcffilter -f "NS = 2" | wc -l
6685
获得临床株与环境株之间的SNP差异数量，结果表明这2株之间只有5个SNPs。
~$ vcfsamplediff SAME TAP1 PATIENT1 compare_tap1.vcf | vcffilter -f "QUAL / AO > 10" | vcffilter -f "NS = 2" | vcffilter -f "! ( SAME = germline ) " | grep -v "^#" | wc -l
5
用vcfsamplediff了解SNPs差异的基因，可以知道是那些gene/CDS里面发生了变化，还可以根据定位看突变是否是同义突变。如果snps数量接近无法区分，但是定位的基因有差异，那么菌株的溯源还需要进一步考虑。
~$ vcfsamplediff SAME PATIENT1 TAP1 compare_tap1.vcf | vcffilter -f "QUAL / AO > 10" | vcffilter -f "NS = 2" | vcffilter -f "! ( SAME = germline ) " > tap_differences.vcf
~$ bedtools intersect -a ST17.gff -b tap_differences.vcf | awk -F";" '{print NR, $2}' OFS="\t"
结果输出：
1 gene=yfiR
2 inference=ab initio prediction:Prodigal:2.60,similar to AA sequence:RefSeq:YP_261793.1,protein motif:TIGRFAMs:TIGR03756,protein motif:Pfam:PF06834.5
3 inference=ab initio prediction:Prodigal:2.60,similar to AA sequence:RefSeq:YP_261793.1,protein motif:TIGRFAMs:TIGR03756,protein motif:Pfam:PF06834.5
4 gene=trbL
5 inference=ab initio prediction:Prodigal:2.60,similar to AA sequence:RefSeq:YP_001349879.1,protein motif:Cdd:COG3481,protein motif:TIGRFAMs:TIGR03760,protein motif:Pfam:PF07514.5
了解2个临床株的SNPs差异，并查看这些SNPs所影响的gene
~$ freebayes --ploidy 1 -C 5 -f ST17.fasta PATIENT1.sorted.bam PATIENT2.sorted.bam > compare_patient.vcf
~$ vcfsamplediff SAME PATIENT1 TAP1 compare_patient.vcf | vcffilter -f "QUAL / AO > 10" | vcffilter -f "NS = 2" | vcffilter -f "! ( SAME = germline ) " > patient_differences.vcf
~$ bedtools intersect -a ST17.gff -b patient_differences.vcf | awk -F";" '{print $2}'
3. 用 Snippy 来实现¶
大量的日常工作或者不熟悉各种命令行工具使用的话，就需要更简便的软件或者pipelines来帮助完成整个分析工作。因此这里我们使用snippy来实现分析SNPs，并绘制基于snps的进化树。
Snippy是比较适合新手的工具，它提供了All in One的工具套装。虽然Snippy需要以下软件支持：
- Perl >= 5.6
- BioPerl >= 1.6
- bwa mem >= 0.7.12
- samtools >= 1.1
- GNU parallel > 2013xxxx
- freebayes >= 0.9.20
- freebayes sripts ( freebayes-parallel,fasta_generate_regions.py)
- vcflib (vcffilter, vcfstreamsort, vcfuniq, vcffirstheader)
- vcftools (vcf-consensus)
- snpEff >= 4.1
但这些工具在Snippy安装包里都已经提供了，我们只需要根据自己的系统设置将以来工具添加到系统路径中去即可。对于perl的一些模块，可能需要更新后才能正常使用：sudo cpan -u module_name
~$ wget https://github.com/tseemann/snippy/archive/v2.9.tar.gz
~$ tar xzf v2.9.tar.gz -C ~/app
~$ echo "export PATH=$PATH:$HOME/app/snippy-2.9/bin:$HOME/app/snippy-2.9/binaries/linux/" >> ~/.bashrc
~$ source ~/.bashrc
# update all perl module
~$ sudo cpan -u
Snippy 不仅可以获得SNP(包括MultiSNP)，也可以获得insertion, indeletion以及Comibination。文本格式的结果记录文件中的snps.tab中，也可以用浏览器打开snps.html查看。
~$ snippy --cpus 4 --outdir output1 --ref ST17.fasta \
> --R1 TAP1_R1_paired.fastq.gz --R2 TAP1_R2_paired.fastq.gz
~$ head output/snps.tab
# snps.txt对各种突变位点做一个数据统计
~$ cat output1/snps.txt
...
Software snippy 2.9
Variant-COMPLEX 1124
Variant-DEL 27
Variant-INS 43
Variant-MNP 222
Variant-SNP 3351
VariantTotal 4767
# 计算不同
# 了解不同碱基的SNP变化数量，如T->C的SNP数量
~$ awk '$10=="T=>C" {n++} END{print n}' output1/snps.gff
Snippy还可以生成多个基因组的共有SNPs的比对文件。用snippy分别生成3个TAP的snp列表数据到ouput*目录中，然后统计共有snps数量并生成snp序列文件aln。
# for old ubuntu system you need update outdated perl module List::Util
~$ sudo cpan -u List::Util
~$ snippy-core --prefix core output1 output2 output3
...
Found 30711 core SNPs from 42075 variant sites.
Saved SNP table: core.tab
Constructing alignment object for core.aln
...
如果要用raxml构建进化树，那获得的.aln文件还需要转换成.phy才可以使用。很多工具可以实现转换，这里使用网络服务来实现：http://sing.ei.uvigo.es/ALTER/
~$ raxmlHPC -f a -x 12345 -p 12345 -# 100 -m GTRGAMMA -s core.phy -n ex -T 4
~$ figtree RAxML_bestTree.ex
因为比对的是fasta格式的文件，所以snps不止是编码区CDS的，而是整个基因组上的snps。
4. 用 Workflow 脚本来简化操作¶
对于我们实验室平时会有许多类似的SNP鉴定工作，那么下机数据要一个个操作会很繁琐，有时候我们通过 unix 的 pipeline 来简化或者写 pipeline 脚本来实现。而另一种解决方式是通过 Workflow 工具来自动化数据分析。专门的 Workflow 软件有许多，有些带GUI（如 galaxy 就可以建立图形化的 workflow），有些基于CML。对于在服务器上的生物信息学分析，我们一般都是丢给命令行来实现。Workflow 因此，这里使用基于 python 的 ruffus 来建立 workflow，作为一个示例演示如何简化平时大量重复的工作。
如果想了解更多关于 workflow/pipeline 的软件，可以查看以下网站：
undone
5. 如何选择 Reference¶
在做微生物 SNP calling时，选择不同 reference 对结果是否会有影响？会有多大影响？这是具体分析时需要考虑的问题。有一些文章进行了讨论。对于我们实验室开展病原微生物溯源时，由于菌株一般是高度近源的，所以影响不大。对于做进化关系分析的，还是需要考虑snp的假阳性问题，通常是要对 SNP 做 filtering。而且对于不同物种考虑其他的突变事件与重组在进化上的影响。
另一方面不同软件的流程与参数不同，对于不熟悉生物信息学的生物学家当使用All in One之类的套件工具分析时，可能会发现不同软件的结果会有不小的差异。因此还是要从自己研究的物种特点角度去选择软件，所以说最好对于这类工具的具体细节需要进行了解，那么对于文档不全或者不够丰富的软件来说，就要仔细研读它的源代码来了解其实现方式。
第三章. 比较基因组数据分析¶
第1节. 测序数据构建 Pangenomics 分析¶
数据准备：
第2节. Orthologs 分析¶
第3节. 基因组进化分析¶
Phylogenomics¶
病原的溯源，进化以及种群结构是公共卫生研究领域一个非常重要内容。过去基于一代测序的方法，对一个基因（如毒力基因序列）或几个基因（如MLST序列）构建进化树，数据量的不足带来的精度问题，容易造成结果与实际产生偏差。而随着NGS的发展，产生了一门新的学科 Phylogenomics 基因组种系学。通过全基因组的数据来构建绘制系统发生树。
phylogenetics analysis¶
在讨论 phylogenomics 之前，我们先对过去 phylogenetics 做一个复习。进化树绘制相信做微生物的工作者平时都在频繁使用，特别是像CDC等做分子流行病学或者一些做微生物系统发育的研究团队。只是目前将进化分析的核酸数据扩展到基因组水平上，绘制树的原理和方法是一致的，只是在之前的多重序列比对方面有一些不同。
- MP
- NJ
- ML
- Bayes
WGA(Whole Genome Alignment)是对2个以上基因组在核苷酸水平上的进化关系预测方法。对于基因组多重序列比对，几个概念搞清楚后才能更好的理解。
- Homologs:
- Orthologs: 来源与共同祖先的基因簇，由于物种差异产生的分化，基因一般具有相似功能。
- Paralogs: 一般是不同功能。
- Xenologs: 2个物种间水平转移的基因，功能可以相似叶可以不同。一般来说是功能相似的。
为了解决 Blast 存在的各种问题，出现了许多针对 WGA 的应用。 目前根据计算策略主要分成2类方式，但这2类方式有各自的优缺点。
- hierarchical: 速度快，容易遗漏小的重配或分化序列片段。如著名的 progressiveMauve, MUGSY 等
- local: 精度高。如 MUMmer, MULTIZ/TBA 等。
WGA 应用列表
phylogenomics analysis¶
对于微生物而言，特别是细菌因为存在大量的重组，水平转移等。不能直接将基因组数据像同源基因数据一样直接比对。所以一般来说需要先将数据做一些筛选。具体的方法主要有；
方法1 基于距离¶
计算全基因矩阵距离，通过距离来构建系统发生树。
方法2 基于核心基因(core genes)¶
通过将同种细菌基因组同源基因（core genes）的异同构建进化树。这种方法的流程一般是将测序数据拼接成contigs/，对拼接结果进行注释，鉴定所有菌株基因组数据的orthologous基因，对这些基因经行多重序列比对，再根据比对结果构建进化树。
方法3 基于参考基因组¶
通过将测序获得的短片段mapping到参考基因组上，找到核心基因组上的variants位点，连锁这些位点进行多重比对，再构建系统发生树。这也是目前最常用的构建微生物基因组SNP进化树的方法流程。
方法2. 应用工具¶
1. orthMCL
方法3. 应用工具¶
1. Wombac
[wombac][] 是由 [Victorian Bioinformatics Consortium](http://www.vicbioinformatics.com/) 开发的用于获取细菌基因组 SNP 的工具集。[wombac][] 使用 perl 开发，本质上是一个 perl pipeline，使用的工具是 bwa, samtools 和 freebayes，这3个工具 Linux 版 [wombac][] 已经自动提供了，要正常运行的话系统中还必须安装vcfutils.pl, bgzip, tabix 这几个依赖包。
[wombac][] 会寻找碱基替换，但不会搜索 indels，精确度上来说可能会有一些 SNP 的丢失，但其精度已足够用来进行大部分基因组的快速分析。
1.1 安装并运行
$ wget http://www.vicbioinformatics.com/wombac-1.3.tar.gz
$ tar zxvf wombac-1.3.tar.gz -C ~/app
$ cd ~/app/wombac-1.3/bin
$ perl wombac --outdir ~/data/wombac --ref ~/data/egd.fasta --run /data/draft
1.2 绘制SNP进化树树
[wombac][] 对每个基因组产生BAM，VCF文件，同时产生一个整体多重比对的ALN文件用于绘制进化树。
$ SplitsTree -i Tree/snps.aln
2. kSNP
3. REALPHY
3.1 下载并安装
RealPhy 需要 java 支持，所以系统要至少安装有 jre 环境。可以输入`java -version`查看，如果提示没有安装java，可以先安装 openjdk：sudo apt-get install open-jre-1.7-headless
$ wget http://realphy.unibas.ch/downloads/REALPHY_v110_exec.zip
$ unzip REALPHY_v110_exec.zip -d ~/apps/realphy
3.2 脚本所依赖工具
Realphy需要以下工具，这些工具的安装已经在前文中实现过了。
- samtools
- bowtie2
- TREE-PUZZLE
- RAxML
- PhyML
- Phylip
3.3 运行程序并分析
$ java -XmM800
4. core-phylogenomics
$ snp_phylogenomics_control --mode mapping --input-dir my_fastq_folder/ --output pipeline_out --reference my_reference.fasta
SeqSphere
One Disrupting Technology Fits it All -Towards Standardized Bacterial Whole Genome Sequencing for Global Surveillance. Dag Harmsen, University of Münster, German.
## Reference
- https://github.com/apetkau/microbial-informatics-2014/
- Single Nucleotide Polymorphisms Methods and Protocols 2nd. Anton A. Komar, Springer protocols.
- Methods in Molecular Biology, Chapter8 Whole Genome Alignment, Colin N. Dewey, Springer.
[wombac]: http://www.vicbioinformatics.com/software.wombac.shtml “wombac” [CSI Phylogeny]: http://cge.cbs.dtu.dk/services/CSIPhylogeny/ “CSI Phylogeny” [GenoBox]: https://github.com/srcbs/GenoBox “GenoBox” [REALPHY]: http://realphy.unibas.ch/fcgi/realphy “REALPHY” [kSNP]: http://sourceforge.net/projects/ksnp/ “kSNP v2”
分析 Homologs¶
1. 软件 get_homologues¶
1.1 安装 get_homologues¶
$ wget ...
$ tar xfz get_homologues_2.0.19.tgz -C get_homologues
$ cd get_homologues
$ ./install.pl
# 查看功能与软件信息
$ ./get_homologues.pl -v
1.2 使用 get_homologues¶
Reference:¶
1. 2.
细菌基因组SNP差异构建进化树¶
1. 需要参考基因组比对的软件¶
1.1 REALPHY¶
REALPHY 是一个比较简单易用的软件，它的工作方式是用 Bowtie2 将 fastq 比对到参考基因组获得 SNP 位点，然后用 RaxML 构建 ML 进化树。但是REALPHY 取了一个比较容易起争议的名字，并不是说它创建的进化树就是 “Real” 真正的进化树，这个软件的特点是可以选择不同的参考基因组。
REALPHY 依赖以下软件：
- Bowtie2
- samtools
- RAxML/PhyML
- TREE-PUZZLE
- Phylip
下载安装
# 安装依赖库
$ sudo apt install libtool pkg-config zlib1g-dev libncurses5-dev
# 安装 Java，Ubuntu16.04 目前默认的版本采用 openjdk 1.8
$ sudo apt install default-jre
# 下载安装 REALPHY
$ wget http://realphy.unibas.ch/downloads/REALPHY_v112_exec.zip
$ unzip REALPHY_v112_exec.zip
$ sudo mv REALPHY_v112_exec /opt/realphy && sudo chown -R root:root /opt/realphy
$ sudo touch /opt/realphy/v1.12
# 创建一个启动脚本，命名为realphy
$ sudo touch /usr/local/sbin/realphy
$ sudo chmod +x /usr/local/sbin/realphy
修改 realphy 脚本，内容如下。其中可以根据机器配置和基因组数量调整参数 -Xmx 增大或者减少内存使用数量。
#!/bin/bash
java -Xmx64000m -jar /opt/realphy/RealPhy_v112.jar
使用Docker容器安装
# ---------------------------------
# Realphy Dockerfile
# ---------------------------------
FROM ubuntu:16.04
MAINTAINER Docker indexofire <indexofire@gmail.com>
RUN apt-get -qq update
RUN apt-get -y upgrade
RUN apt-get install -y build-essential
RUN DEBIAN_FRONTEND=noninteractive apt-get install -y software-properties-common
RUN DEBIAN_FRONTEND=noninteractive apt-get install -y python-software-properties
RUN apt-get install -y default-jre
RUN apt-get install -y git-core
RUN apt-get install -y libtool
RUN apt-get install -y pkg-config
RUN apt-get install -y zlib1g-dev
RUN apt-get install -y libncurses5-dev
RUN apt-get install -y r-base
构建进化树
- 最基本用法，将参考基因组 reference（最好是完成的基因组） 放入 input 文件夹，可以是 fasta 和 genbank 格式，如果放置多个 fasta/gbk 文件，软件会随机选择其中一个作为参考基因组，并将其他序列作为比对序列。如果指定某个序列作为参考基因组，用 -ref 参数指定。fastq 数据作为比对数据放入 input 文件夹，如果是PE测序的数据，双端的fastq文件名要用 *_R1.fastq/*_R2.fastq 区分。fastq 数据建议先对测序数据进行一定的QC操作，尽可能减少测序错误引入的 SNP 位点。
$ realphy input output -ref ref_genome
- 对核心基因组构建进化树。默认的方式是对参比基因组全部区域比对，包括非编码区。如果希望仅对CDS进行比对，就可以通过添加 -genes 参数。如果添加了 -genes 参数，那么参考基因组文件格式必须是含有 CDS 的 genbank 格式数据。
$ realphy input output -ref ref_genome -genes
- 使用多个参考基因组。这是 REALPHY 的一个“卖点”，作者认为当参考基因组与比对基因组差异较大时，如>5%，就会造成进化树构建的偏差。所以通过选择多个参考基因组，平衡不同型别的基因组，减少远源基因组带来的比对差异。通过设置 -refN 来实现多个参考基因组，并用 -merge 参数合并结果。由于计算量较大，不建议使用太多参考基因组。
$ realphy input output -ref1 ref1_genome -ref2 ref2_genome -ref3 ref3_genome -merge
建议使用 RAxML 来构建基于 ML 的进化树。因此生成的进化树图为 “PolySeqOut/RAxML_bestTree.raxml”，将其复制到带 GUI 的计算机上，用画树软件生成树图。
$ scp user@server-ip:~/data/REF/PolySeqOut/RAxML_bestTree.raxml .
$ figtree RAxML_bestTree.raxml
1.2 Snippy¶
下载安装
$ git clone https://github.com/tseemann/snippy
$ echo 'export PATH="$HOME/repos/snippy/bin:$HOME/repos/snippy/binary/linux:$PATH"' >> ~/.bashrc
$ source ~/.bashrc
使用Snippy
# 单个样本用 snippy 分别比对获得每个样本的 snps
$ snippy --cpus 20 --outdir S01 -ref reference.fa --R1 S01_R1_L001.fastq.gz --R2 S01_R2_L001.fastq.gz
$ snippy --cpus 20 --outdir S02 -ref reference.fa --R1 S02_R1_L001.fastq.gz --R2 S02_R2_L001.fastq.gz
...
# 将所有 snps 结果汇总，获得 core snps
$ snippy-core --prefix core-snps S01 S02 S03 ...
# 单行脚本完成所有样本 snps 比对工作
$ for i in $(awk -F'_' '{print $1}' <(ls -D *.fastq.gz) | sort | uniq); \
> do snippy --cpus 20 --outdir $i -ref reference.fa --R1 $i*R1*.fastq.gz --R2 $i*R2*.fastq.gz; \
> done
$ snippy-core --prefix core-snps S*
# 生成 clustal 格式的 core.aln 比对文件，这里用 splittree 构建进化树
$ splittree -i core-snps/core.aln
# 或者将生成的 .aln 格式比对文件转换成 .phy 格式，然后再用 raxml 构建 ML 树
$ raxml -f a -x 12345 -# 100 -p 12345 -m GTRGAMMA -s core.aln -n ex -T 40
2. 不需要参考基因组的软件¶
2.1 kSNP¶
kSNP 采用的是基于 kmer 的算法，不需要进行序列的多重比对，也不需要参考基因组进行 Mapping。kSNP 可以分析细菌完成基因组和未完成基因组，也可以直接分析测序 reads 数据。
下载安装
# 下载软件包，并解压缩到 /opt/ksnp 目录
$ wget http://nchc.dl.sourceforge.net/project/ksnp/kSNP3.021_Linux_package.zip
$ sudo unzip kSNP3.021_Linux_package.zip
$ sudo mv kSNP3.021_Linux_package /opt/ksnp
$ cd /opt/ksnp
# 添加版本信息文件
$ touch v3.021
使用kSNP3
$
第四章. 基因组数据可视化¶
用 circos 来画基因组¶
circos 是一个绘制可视化关系图的软件，就是我们常说的用来做 Data Visulization。它的特点是采用圈图来表示关系，因此也常被用来绘制环状的原核生物染色体或质粒的基因组图谱，或者用来展示比较基因组关系图等。circos 生成的圈图非常优美，可以用来表示多种类型的数据。但是由于功能强大，因此牵涉绘图的内容过于复杂，参数设置繁多，采用的术语名称的含义不够直观，导致学习曲线比较陡峭。而且缺乏易用的所见即所得的工具，因此对于新手来说难度非常大。
circos 的作为绘制图形的工具，不是一步速成的获得某种图形的方法，而是需要通过自己设计的图形样式来定义各种参数，让程序绘制出图。
1. 安装 circos¶
1.1 使用 ubuntu 官方软件源安装¶
可以使用 ubuntu 官方软件源来安装。这样安装最简单，缺点是不是最新版，可能会有一些程序BUG。
~$ sudo apt install circos
1.2 使用官方安装包¶
# 安装依赖包 gd2-xpm 库
~$ sudo apt-get install libgd2-xpm-dev
# 修改 ubuntu 发行版的 perl 路径
~$ sudo ln -s /usr/bin/env /bin/env
# 安装 perl 模块
~$ sudo cpan install Clone \
> Config::General \
> Font::TTF \
> List::MoreUtils \
> Math::Bezier \
> Math::VecStat \
> Math::Round \
> Params::Validate \
> Readonly \
> Regexp::Common \
> Set::IntSpan \
> Text::Format \
> GD \
> Statistics::Basic
# 下载安装包，并解压缩
~$ wet http://circos.ca/distribution/circos-0.69-3.tgz
~$ tar xfz circos-0.69-3.tgz -C ~/apps/circos
# 将 circos 添加到系统路径中
~$ sudo ln -s ~/apps/circos/bin/circos /usr/local/sbin/circos
# 或者将 circos 添加到用户路径中
~$ export PATH=$HOME/apps/circos/bin/
# 查看是否安装成功
~$ circos -v
# 如果命令行打印出类似以下内容，表示安装成功
circos | v 0.69-3 | 16 Jun 2016 | Perl 5.018002
2. circos 基本概念¶
circos 主要需要2类文件：
- 可视化结果的数据文件
- 用于告诉circos如何展示的配置文件。
2.1 数据文件¶
由于circos只是画图工具，因此你图中要展示的数据需要自己预先准备好。比如基因组的GC含量，需要用其他程序先生成以下这种类型的文本格式文件。
染色体 开始位置 结束位置 颜色值 其他选项[可选]
2.chr1 1000 2000 red
3. 可视化的组成要素¶
<plots></plots>，<links></links> 区块内的参数是属于 Global 的。 <plot></plot> 区块内的参数是属于 Local 的。
如果是某个参数在多个block中都重复使用，那么可以设置成 global 的参数，当某个block要使用新的值时，再用 local 的参数覆盖即可。
# plot可以理解为图层，多个圈图图层构成一个plots
<plot></plot>
# plots 由
# type 描述图层内容类型，如text
type = text
# r0, r1 表示圈图半径，<1在内圈，>1在外圈, r0-r1为宽度
r0 = 0.95r
r1 = 0.85r
#
- ideogram
- highlight
- tick
- label
- link
circos.conf 参数
# 染色体的数据集
karyotype = data.txt
# 定义圈图标尺
chromosomes_units = 1000
# 如果设置成yes，所有ideogram上都会现实ticks
chromosomes_display_default = yes
3.4 circos 常量的单位¶
4. 用 circos 来绘图¶
4.1 绘制环状 Salmonella LT2 基因组¶
之前的教程里介绍了用DNAPlotter绘制一个基因组完成图的圈图，圆圈内容strand，GC content和GC skew。从某个角度上看，DNAplotter的配置文件与circos有一定的类似，2个软件的区别在于前者专注于绘制单个细菌基因组，后者可以完成更多的功能。第一个例子我们通过用circos来绘制类似DNAplotter生成的基因组圈图，从而学习circos的基本用法，
从NCBI上下载Salonella LT基因组数据文件。
~$ esearch -db nuccore -query "NC_003197.1[accn]" | efetch -db nuccore -format fasta > LT2.fasta
绘制 circos 所需要的 ideogram 所需要的几类数据，我们通过一个 python 脚本 get_circos_data.py 来获得。
运行脚本：
~$ python get_circos_data.py
4.1.2 配置文件的准备
建立circos.conf配置文件
~$ cd circos_data ~/circos_data$ touch circos.conf
R语言与ggplot2¶
~$ sudo apt install r-base
~$ R
R version 3.2.5 (2016-04-14) -- "Very, Very Secure Dishes"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
...
# 安装ggplot2
> install.packages("ggplot2")
~$ sudo echo "deb http://cran.rstudio.com/bin/linux/ubuntu xenial/" | sudo tee -a /etc/apt/sources.list
~$ gpg --keyserver keyserver.ubuntu.com --recv-key E084DAB9
~$ gpg -a --export E084DAB9 | sudo apt-key add -
~$ sudo apt update
~$ sudo apt install r-base r-base-dev
# 安装 R Studio
~$ sudo apt-get install gdebi-core
~$ wget https://download1.rstudio.org/rstudio-0.99.896-amd64.deb
~$ sudo gdebi -n rstudio-0.99.896-amd64.deb
第五章. 分析工具部署¶
目前 NGS 的生物信息学分析工具多如牛毛
使用Galaxy构建基于Web的应用平台¶
小型实验室快速解决方案¶
使用 Amazon EC2 对于不熟悉 Linux 和服务器配置的用户而言是很方便，节省用户在硬件的经济投入，配置学习以及维护能力的培养。但对于测序数据不适合上传到第三方服务器的实验室或要结合本地 LIMS 系统的实验室，要快速构建 Galaxy 以及相配套的其他软件，目前来说最简便的方法之一是利用 docker 部署到本地服务器上。
Docker¶
Ubuntu 14.04 已经自带了 docker 安装包，不过版本相对较老，想使用 docker 最新版的功能，添加源来安装 docker:
# 添加 docker 源
~$ sudo apt-get install apt-transport-https
~$ sudo apt-key adv --keyserver hkp://keyserver.ubuntu.com:80 --recv-keys 36A1D7869245C8950F966E92D8576A8BA88D21E9
~$ sudo bash -c "echo deb https://get.docker.io/ubuntu docker main > /etc/apt/sources.list.d/docker.list"
# 安装 docker
~$ sudo apt-get update
~$ sudo apt-get install lxc-docker
# 安装完成后启动 docker 服务
~$ sudo service docker start
首先要对 docker 镜像和容器的概念有所了解。细节可以参阅 Docker 从入门到实践
~$ docker run ubuntu:14.04 /bin/echo 'Hello world'
galaxy 的 docker
镜像可以自己来创建，建议使用docker galaxy-stable，源代码可以在 Github 下载。也可以直接到官方 Hun Registry 里下载galaxy-stable镜像
~$ sudo docker pull bgruening/galaxy-stable
官方 Registry 非常慢，可以通过国内 daocloud 服务商的加速器来加快docker pull过程：首先注册一个 daocloud 用户，然后在命令行中添加 registry mirror。
# ****** 为daocloud分配给你的ID
~$ echo "DOCKER_OPTS=\"\$DOCKER_OPTS --registry-mirror=http://******.m.daocloud.io\"" | sudo tee -a /etc/default/docker
~$ sudo service docker restart
~$ sudo docker pull bgruening/galaxy-stable
该脚本可以将--registry-mirror加入到你的Docker配置文件/etc/default/docker中。适用于Ubuntu14.04，其他版本根据docker配置文件和环境变量，可能有细微不同。
~/app$ wget https://github.com/bgruening/docker-galaxy-stable/archive/15.10.tar.gz
~/app$ tar zxf 15.10.tar.gz
~/app$ cd docker-galaxy-stable-15.10
docker run -d表示以daemon方式运行docker，执行该命令后，docker完成galaxy初始化需要花几分钟时间才能访问。如果是本地运行，用浏览器访问127.0.0.1:8080可以看到galaxy首页，如果是服务器上用ssh连接执行命令，要访问服务器IP加上8080端口。
# 后台运行 galaxy 容器
~$ sudo docker run -d -p 8080:80 -p 8021:21 bgruening/galaxy-stable
# 交互方式运行 galaxy 容器
~$ sudo docker run -i -t -p 8080:80 bgruening/galaxy-stable /bin/bash
root$ sh run.sh
有时候需要进入docker容器中进行操作，就可以以docker run -i交互模式进行访问。进入容器后运行sh run.sh可以DEBUG方式运行galaxy，适合本地测试使用。
# fc3e62e0471d 是想要获取文件的所在容器。foo.txt是想要获得的文件。
~$ sudo docker cp fc3ea62e471d:/home/foo.txt .
需要分析的数据通过添加外部数据卷来实现。
# 添加服务器上的 /mydata 卷到容器中
~$ sudo docker create -v /mydata --name my_data_vol bgruening/galaxy-stable /bin/bash
# 或者在运行时将本地卷`/mydata`加入到容器中`/container_data`位置
~$ sudo docker run -d -p 8080:80 -v /mydata:/container_data/ bgruening/galaxy-stable
镜像是只读的，当Ctrl+D方式退出容器后，再次进入容器时你上次以添加的内容是看不到的。如果想要从上一次运行的容易中获得文件可以用docker cp 的方法，不过你得记住上一次运行的 container id 号。
~$ sudo docker ps -a | cut -d ' ' -f 1 | sudo xargs docker rm
# IMAGE ID 是该镜像的ID，如果镜像还有容器运行，或有其他镜像的依赖关系，则无法删除要先删除容器或其他镜像。
~$ sudo docker rmi IMAGE_ID
Rreference¶
本地安装与基本配置¶
本节介绍 Galaxy 的最基本下载安装与使用。最适合的场景为个人电脑，单用户使用的情况。
下载与安装¶
Galaxy 作为一款开源软件，其代码库托管在 http://bitbucket.org ，先安装 mercurial ，然后用 hg 工具将 galaxy 代码库克隆到本地。
:: code-block:: bash
~app$ sudo apt-get install mercurial ~app$ hg clone https://bitbucket.org/galaxy/galaxy-dist/ ~app$ cd galaxy-dist
可以用 hg branch 命令查看代码分支是否为
stable，如果是其他分支（galaxy代码库有另一分支’default’），切换到 stable
分支：hg update stable。在生产环境下建议使用 stable 分支的代码。
克隆到本地的 stable 代码，一般Linux系统自带Python就可以直接运行了。
~app$ ./run.sh
运行run.sh，这个shell脚本程序会自动完成初始化数据，依赖库下载，数据库迁移等一系列操作，当看到终端显示serving on http://127.0.0.1:8080时，可以打开浏览器，访问 http://127.0.0.1:8080 即可看到galaxy的界面。
默认的galaxy只带有基本的工具，对于实际工作中需要的各种分析软件，需要添加到自己建立的 galaxy 实例中。
高通量测序的生物信息学软件大多是基于命令行的开源工具。galaxy利用python语言将这些工具粘合到galaxy实例中，使得用户可以在web界面中直接调用命令行工具对数据进行操作。
galaxy有一个toolshed（工具库）的概念：https://toolshed.g2.bx.psu.edu/（官方维护的toolshed），许多著名的工具已经被移植到toolshed中，可以直接被安装到galaxy里，此外也有许多第三方的toolshed包可以添加，甚至掌握了一些 python 脚本和 galaxy xml规范后，也可以自己添加一些分析工具到galaxy中。
除了分析软件外，toolshed还包含创建的数据类型，以及工作流等。首先修改配置文件 tool_sheds_conf.xml
~$ cp config/tool_sheds_conf.xml.sample config/tool_sheds_conf.xml
其次在galaxy.ini中配置依赖包的安装目录（上一部分已经添加了这个参数，将依赖包安装在tool_dep），然后添加管理员帐号，比如你之前用admin@localhost.com注册的galaxy实例，就将galaxy.ini中admin_users设置为：
admin_users = admin@localhost.com
重启你的galaxy实例，用admin@localhost.com用户登陆，你就有权限访问http://127.0.0.1:8080/admin，在admin界面可以看到Tool sheds下有Search and browse tool sheds链接，点击后可以看到默认的2个toolsheds源。
进入Galaxy main tool shed，工具列表上访有搜索框，在这里输入你要安装的工具名称比如spades后，回车进行检索。
Instance
点击结果列表中的spades下拉菜单，选择Preview and install，在转向页面中点击Install to Galaxy后，会出现如下图的提示。
Instance
在add new tool panel section中输入Assembly，将spades工具归类到Assembly这个新建的工具类别中，点击页面底部的Install按钮开始安装。
构建单机产环境
————–
作为个人尝试，前面的步骤在PC机上已经可以正常运行使用。对于要作为生产环境下多用户使用，建议使用专门的代理服务器和数据库来增强效率和速度。这里采用nginx+postgresql构建生产环境下的 Galaxy 服务。这里只是简单的介绍一下最基本的配置方式，对于高负载的web server设置又是另一个很复杂的话题了，这里就不具体展开。也可以参考别人做的galaxy dockerfile。
首先在ubuntu下新建一个用户galaxy。
~$ sudo adduser galaxy
按照提示，在弹出提示符输入相应内容，主要填好密码即可，其他可以留空。然后切换到galaxy用户：
~$ su galaxy
~$ cd
重复Galaxy 本地安装与配置的第一步下载与安装步骤,建立配置文件：
~$ cp config/galaxy.ini.sample config/galaxy.ini
~$ vim config/galaxy.ini
将配置文件galaxy.ini设置如下，你也可以将内容直接保存成galaxy.ini
[server:main]
use = egg:Paste#http
host = 0.0.0.0
use_threadpool = True
threadpool_kill_thread_limit = 10800
[filter:gzip]
use = egg:Paste#gzip
[filter:proxy-prefix]
use = egg:PasteDeploy#prefix
prefix = /galaxy
[app:main]
paste.app_factory = galaxy.web.buildapp:app_factory
database_connection = postgresql://galaxy:galaxy@localhost:5432/galaxyserver
tool_dependency_dir = tool_dep
use_nglims = False
nglims_config_file = tool-data/nglims.yaml
debug = False
use_interactive = False
admin_users = admin@localhost.com
安装postgresql并建立galaxy数据库，这里用户名和密码都设置为galaxy，要与
galaxy.ini 中database_connection参数对应的值一致。
~$ sudo apt-get install postgresql-9.3
~$ su - postgres
~$ psql template1
> CREATE USER galaxy WITH PASSWORD 'galaxy';
> CREATE DATABASE galaxyserver;
> GRANT ALL PRIVILEGES ON DATABASE galaxyserver to galaxy;
> \q
用nginx做反向代理，处理请求。
~$ sudo apt-get install nginx
设置nginx.conf
http {
upstream galaxy_app {
server localhost:8080;
}
server {
client_max_body_size 10G;
location / {
proxy_pass http://galaxy_app;
proxy_set_header X-Forwarded-Host $host;
proxy_set_header X-Forwarded-For $proxy_add_x_forwarded_for;
}
}
}
第六章. Metagenomics数据分析¶
第1节. Metagenomics 概述¶
注解
Metagenomics 实验操作和数据分析主要包括2个方面：
- shotgun metagenomics
- 16s rRNA metagenomics
1. shotgun Metagenomics¶
我们在第1-3节介绍如何用软件处理 shotgun metagenomics 测序数据。
Metegenomics de novo assembly 中的几个难题：
- 重复序列
- 低覆盖度
- 测序错误
重复序列¶
重复序列是
低覆盖度¶
测序错误¶
2. 16s rRNA Metagenomics¶
16s rRNA 数据有偏向性：
- Degenerate primers Degenerate primers
- PCR Amplification
- Databases
- Does not capture viruses and eukaryotes
主要用于菌群分析。
用途：需要深度测序的数据，样本一般分为：宿主来源，如人或动物等；环境来源，如土壤或水体等。数据分析时，前者一定要考去除宿主缘DNA，后者要高度深度测序对后期拼接有较大帮助。
Methods Protocol:
- DNA isolation
- Library preparation
- Sequencing
- Data QC: FastQC, FastX_toolkit
- Reads assembly: IDBA-UD, Mira, Velvet/MetaVelvet, Spades
- Partipation: FCP，Phymm/PhymmBL, AMPHORA2, NBC, MEGAN, MG-RAST, CAMERA, IMG/M
- Analysis: MetaPhAln, PhyloSift/pplacer
第2节. 使用 Centrifuge 分析¶
Centrifuge 是一个快速检索 reads 的软件，可以用于对微生物样本测序后的宏基因组数据分析。
1. 安装 Centrifuge¶
# 安装软件
$ git clone https://github.com/infphilo/centrifuge
$ cd centrifuge
$ make
$ sudo make install
# 下载数据库
$ cd indices
$ make b+h+v
2. 使用 Centrifuge¶
第3节. 使用 Clark-S 分析¶
第4节. 使用 Kraken 分析¶
第4节. 使用 Meta-BEETL 分析¶
BEETL 是一套做 Burrows-Wheeler Transform 的工具集，Meta-BEETL 是 BEETL 下的一个子项目，用来分析基于 Shotgun metagenomics 的微生物数据。
1. 构建数据库¶
1.1 直接下载已整理的数据库¶
作者已经将研究所用到的数据库整理后提交到 Amazon S3 服务器上，可以直接下载。总数据量在24G左右。
~$ mkdir BeetlMetagenomeDatabase && cd BeetlMetagenomeDatabase
~/BeetlMetagenomDatabase$ for i in `curl https://s3.amazonaws.com/metaBEETL | grep -oP "[^>]*bz2"` ; \
> do wget https://s3.amazonaws.com/metaBEETL/$i & done
静态文件的 URL 地址，也可以直接下载。
https://s3.amazonaws.com/metaBEETL/filecounter.csv.bz2
https://s3.amazonaws.com/metaBEETL/headerFile.csv.bz2
https://s3.amazonaws.com/metaBEETL/names.dmp.bz2
https://s3.amazonaws.com/metaBEETL/ncbiFileNumToTaxTree.bz2
https://s3.amazonaws.com/metaBEETL/ncbiMicros-B00.bz2
https://s3.amazonaws.com/metaBEETL/ncbiMicros-B01.bz2
https://s3.amazonaws.com/metaBEETL/ncbiMicros-B02.bz2
https://s3.amazonaws.com/metaBEETL/ncbiMicros-B03.bz2
https://s3.amazonaws.com/metaBEETL/ncbiMicros-B04.bz2
https://s3.amazonaws.com/metaBEETL/ncbiMicros-B05.bz2
https://s3.amazonaws.com/metaBEETL/ncbiMicros-B06.bz2
https://s3.amazonaws.com/metaBEETL/ncbiMicros-C00.bz2
https://s3.amazonaws.com/metaBEETL/ncbiMicros-C01.bz2
https://s3.amazonaws.com/metaBEETL/ncbiMicros-C02.bz2
https://s3.amazonaws.com/metaBEETL/ncbiMicros-C03.bz2
https://s3.amazonaws.com/metaBEETL/ncbiMicros-C04.bz2
https://s3.amazonaws.com/metaBEETL/ncbiMicros-C05.bz2
https://s3.amazonaws.com/metaBEETL/ncbiMicros-C06.bz2
下载完成后解压缩。
~/BeetlMetagenomDatabase$ bunzip2 *.bz2
~/BeetlMetagenomDatabase$ export METAGENOME_DATABASE_PATH=`pwd`
1.2 构建自定义数据库¶
处于实际目的，有时候我们不需要这么大的数据库，或者需要更多其他的数据加入到数据库，那就需要自行建立数据库。
SeqAn 是一个高效的 C++
算法和数据结构库，用于分析生物序列等用途。访问 SeqAn
下载并安装。安装完成后编译BEETL时可以添加--with-seqan参数。
# 安装依赖包
~$ sudo apt-get install g++ cmake python zlib1g-dev libbz2-dev libboost-dev
~$ cd ~/apps
~/apps$ git clone https://github.com/seqan/seqan.git && cd seqan
~/apps/seqan$ git checkout -b develop origin/develop
下载需要的序列，比如NCBI上细菌的基因组 ref 数据。
~/BeetlMetagenomDatabase$ wget ftp://ftp.ncbi.nlm.nih.gov/genomes/Bacteria/all.fna.tar.gz -P downloads
~/BeetlMetagenomDatabase$ mkdir all.fna
~/BeetlMetagenomDatabase$ tar xzf downloads/all.fna.tar.gz -C all.fna/
~/BeetlMetagenomDatabase$ metabeetl-db-genomesToSingleSeq -f all.fna/*/*.fna -s singleSeqGenomes
~/BeetlMetagenomDatabase$ cd singleSeqGenomes
修改metabeetl-db-arrayBWT.sh文件里数据库路径，创建 G_* 和 G_*_rev 文件的BWTs，用qsub来向系统提交任务。也可以运行metabeetl-db-makeBWTSkew来一个个创建。
~/BeetlMetagenomDatabase$ cp `which metabeetl-db-arrayBWT.sh` .
~/BeetlMetagenomDatabase$ vim metabeetl-db-arrayBWT.sh
# 如果文件是 G_1 - G_500，那么n=1-500
~/BeetlMetagenomDatabase$ qsub -t n metabeetl-db-arrayBWT.sh
# 或者
~/BeetlMetagenomDatabase$ (echo -n "all: " ; for i in G_*; do echo -n " bwt_${i}-B00"; done ; echo -e "\n" ; \
> for i in G_*; do echo "bwt_${i}-B00: ${i}"; echo -e "\tmetabeetl-db-makeBWTSkew ${i} ${i}\n" ; done ) > Makefile
~$ make -j
在一台内存>60G的机器上将序列加载到内存中，并吧所有文件合并。
~/BeetlMetagenomDatabase$ for pileNum in `seq 0 5`; do metabeetl-db-mergeBacteria $pileNum ncbiMicros <( ls G_* ) ; done
生成的文件名称类似：ncbiMicros-A0*, -B0* and -C0*
- -A0* 包含每个BWT的位置
- -B0* BWTs
- -C0* 包含每个BWT位置的文件号
文件 -B0* and -C0* 用于计数算法
可选操作：将BWT文件转换成RLE BWT格式，运行可以更快。
~/BeetlMetagenomDatabase$ for pileNum in `seq 0 5`; do \
> mv ncbiMicros-B0${pileNum} ncbiMicros-B0${pileNum}.ascii ; \
> beetl-convert \
> --input-format=bwt_ascii \
> --output-format=bwt_rle \
> -i ncbiMicros-B0${pileNum}.ascii \
> -o ncbiMicros-B0${pileNum} ; \
> done
~$ cd tmp
~/tmp$ wget ftp://ftp.ncbi.nih.gov/pub/taxonomy/taxdump.tar.gz
~/tmp$ wget ftp://ftp.ncbi.nih.gov/pub/taxonomy/gi_taxid_nucl.dmp.gz
~/tmp$ tar xzf taxdump.tar.gz names.dmp nodes.dmp
~/tmp$ gunzip gi_taxid_nucl.dmp.gz
metabeetl-db-findTaxa脚本查找数据库中的taxonomic树。 Use themetabeetl-db-findTaxa script to find the taxonomic tree corresponding to
the file numbers in the database.
| You will need the headerFile produced by running
“metabeetl-db-genomesToSingleSeq” and fileCounter created during the
merging of the bacterial reference genomes.
| Finally, you get for each file number in the database a taxonomic
tree with the taxonomic ids.
| There will be some 0 in the taxonomic tree. This is a taxonomic id
which could not be matched to: Superkingdom, Phylum, Order, Class,
Family, Genus, Species or Strain.
| Sometimes there are just missing taxa in the taxonomy. We supplement
this with the file metaBeetlExtraNames.dmp below.
~/BeetlMetagenomDatabase$ metabeetl-db-findTaxa \
> -nA downloads/names.dmp \
> -nO downloads/nodes.dmp \
> -nG downloads/gi_taxid_nucl.dmp \
> -h singleSeqGenomes/headerFile.csv \
> -f singleSeqGenomes/filecounter.csv > ncbiFileNumToTaxTree
~/BeetlMEtagenomDatabase$ ( grep scientific downloads/names.dmp ; cat ${BEETL_INSTALL_DIR}/share/beetl/metaBeetlExtraNames.dmp ) > metaBeetlTaxonomyNames.dmp
~/BeetlMetagenomDatabase$ mkdir normalisation
~/BeetlMetagenomDatabase$ cd normalisation
~/BeetlMetagenomDatabase/normalisation$ touch normalize.sh
# normalize.sh
for genome in ../singleSeqGenomes/G_*; do
(
genomeNum=`basename ${genome}`
mkdir ${genomeNum}
cd ${genomeNum}
for i in ../../singleSeqGenomes/bwt_${genomeNum}-B0?; do ln -s ${i} ; done
for i in `seq 0 6`; do touch bwt_${genomeNum}-B0${i}; done
time beetl-compare --mode=metagenomics -a bwt_${genomeNum} -b ../../singleSeqGenomes/ncbiMicros -t ../../ncbiFileNumToTaxTree -w 20 -n 1 -k 50 --no-comparison-skip -v &> out.step1
rm -f BeetlCompareOutput/cycle51.subset*
touch empty.txt
time cat BeetlCompareOutput/cycle*.subset* | metabeetl-convertMetagenomicRangesToTaxa ../../ncbiFileNumToTaxTree ../../singleSeqGenomes/ncbiMicros ../../metaBeetlTaxonomyNames.dmp empty.txt 20 50 - &> out.step2
cd ..
) &
done ; wait
for i in `seq 1 2977`; do echo "G_$i"; X=`grep -P "G_${i}$" ../singleSeqGenomes/filecounter.csv | cut -f 1 -d ','`; TAX=`grep -P "^${X} " ../ncbiFileNumToTaxTree | tr -d "\n\r"` ; echo "TAX(${X}): ${TAX}."; TAXID=`echo "${TAX}" | sed 's/\( 0\)*$//g' | awk '{print $NF}'`; echo "TAXID=${TAXID}"; COUNTS=`grep Counts G_${i}/out.step2 | head -1`; echo "COUNTS=${COUNTS}"; MAIN_COUNT=`echo "${COUNTS} " | sed "s/^.* ${TAXID}:\([0-9]*\) .*$/\1/ ; s/Counts.*/0/"` ; echo "MAIN_COUNT=${MAIN_COUNT}" ; SUM=`echo "${COUNTS} " | tr ' ' '\n' | sed 's/.*://' | awk 'BEGIN { sum=0 } { sum+=$1 } END { print sum }'` ; echo "SUM=$SUM"; PERCENT=`echo -e "scale=5\n100*${MAIN_COUNT}/${SUM}" | bc` ; echo "PERCENT=${PERCENT}" ; echo "FINAL: G_${i} ${TAXID} ${MAIN_COUNT} ${SUM} ${COUNTS}" ; done > r2977
grep FINAL r2977 > ../normalisation.txt
$ ./normalize.sh
~$ mkdir metaBeetlNcbiDb
~$ cd ~/metaBeetlNcbiDb
~$ ln -s ../ncbiFileNumToTaxTree
~$ ln -s ../normalisation.txt
~$ ln -s ../downloads/metaBeetlTaxonomyNames.dmp
~$ ln -s ../singleSeqGenomes/filecounter.csv
~$ ln -s ../singleSeqGenomes/headerFile.csv
~$ for i in ../singleSeqGenomes/ncbiMicros-[BC]0[0-5]; do ln -s $i ; done
2. 分析数据¶
2.1 使用公共数据库中的 shotgun metagenomics 数据¶
用 SRS013948 这个人类肠道细菌组研究项目的数据作为例子。首先下载数据:
# http download
~/data$ wget http://downloads.hmpdacc.org/data/Illumina/throat/SRS013948.tar.bz2
# ftp download
~/data$ wget ftp://public-ftp.hmpdacc.org/Illumina/throat/SRS013948.tar.bz2
# 解压缩
~/data$ tar xjf SRS013948.tar.bz2
~/data$ ls -l
解压缩后可以看到下面3个文件：
SRS013948.denovo_duplicates_marked.trimmed.1.fastq
SRS013948.denovo_duplicates_marked.trimmed.2.fastq
SRS013948.denovo_duplicates_marked.trimmed.singleton.fastq
2.2 转换数据和合并数据¶
~/data$ beetl-convert \
> -i SRS013948.denovo_duplicates_marked.trimmed.1.fastq \
> -o paddedSeq1.seq \
> --sequence-length=100
~/data$ beetl-convert \
> -i SRS013948.denovo_duplicates_marked.trimmed.2.fastq \
> -o paddedSeq2.seq \
> --sequence-length=100
~/data$ beetl-convert \
> -i SRS013948.denovo_duplicates_marked.trimmed.singleton.fastq \
> -o paddedSeqSingleton.seq \
> --sequence-length=100
~/data$ cat paddedSeq1.seq paddedSeq2.seq paddedSeqSingleton.seq > SRS013948.seq
2.3 以 metagenomic 模式运行 BEETL¶
创建 BWT 变换，然后运行 metagenomic 模式的 BEETL
~/data$ beetl-bwt -i SRS013948.seq -o bwt_SRS013948
~/data$ beetl-compare \
> --mode=metagenomics \
> -a bwt_SRS013948 \
> -b ${METAGENOME_DATABASE_PATH}/ncbiMicros \
> -t ${METAGENOME_DATABASE_PATH}/ncbiFileNumToTaxTree \
> -w 20 \
> -n 1 \
> -k 50 \
> --no-comparison-skip
beetl-compare 命令会在文件夹 BeetlCompareOutput
里生成许多文件。k值设置越大，可以获得越多的信息数据，但是输出文件也会变得越大。
2.4 图形化显示结果¶
metabeetl-convertMetagenomicRangesToTaxa
工具根据kmer匹配，从而获得基因组以及上一级分类的ID号，并生成文字与图形结果（采用Krona
js可视化库来形成网页格式的报告文件）。
由于算法原因需要频繁读取BWT位置，因此数据库文件 ncbiMicros-C0*
最好放在读取速度比较块的磁盘扇区或者磁盘里如SSD硬盘。另外如果内存足够，可以用
metabeetl-convertMetagenomicRangesToTaxa_withMmap
工具，将这些比对文件读入内存中加快速度。
~/data$ cat BeetlCompareOutput/cycle*.subset* | \
> metabeetl-convertMetagenomicRangesToTaxa \
> ${METAGENOME_DATABASE_PATH}/ncbiFileNumToTaxTree \
> ${METAGENOME_DATABASE_PATH}/ncbiMicros \
> ${METAGENOME_DATABASE_PATH}/metaBeetlTaxonomyNames.dmp \
> ${METAGENOME_DATABASE_PATH}/normalisation.txt \
> 20 50 - > metaBeetl.log
生成的TSV（制表符分隔）文件各列含义：
| 列 | 字段含义 |
|---|---|
| column 1 | Taxonomy Id |
| column 2 | Taxonomy level |
| column 3 | k-mer count |
| column 4 | k-mer count including children |
TSV文件用途：
- metaBeetl.tsv: raw k-mer counts for every leaves and ancestors.
- metaBeetl_normalised.tsv: some counts from ancestors are moved towards leaf items. The proportion thereof is pre-computed by aligning each individual genome to the full database.
- metaBeetl_normalised2.tsv: Only leaves of the taxonomy tree are kept, and counts are normalised relatively to genome sizes.
Krona JS 结果文件 - metaBeetl_krona.html - metaBeetl_krona_normalised.html - metaBeetl_krona_normalised2.html
第七章. Phylogenomics 分析¶
Snippy¶
snippy 是澳大利亚墨尔本大学的 Torsten Seemann 开发的用于细菌基因组 snps 分析的软件。Torsten Seemann 是著名的细菌基因组注释软件 prokka 的作者，而 snippy 也是被广泛使用的软件之一。snippy 的前身叫做 wombac，
注解
wombac 是由 Victorian Bioinformatics Consortium 开发的用于获取细菌基因组 SNP 分析的工具集。wombac 使用 perl 开发，本质上是一个 perl pipeline，使用的工具是 bwa, samtools, freebayes, vcfutils.pl, bgzip, tabix 等依赖包。wombac 会寻找碱基替换，但不会搜索 indels，精确度上来说可能会有一些 snps 的丢失，但其精度已足够用来进行大部分基因组的快速分析。
# 数据
$ tree data/
S01_R1_L001.fastq.gz
S01_R2_L001.fastq.gz
S02_R1_L001.fastq.gz
S02_R2_L001.fastq.gz
...
1. 下载安装¶
snippy 的软件仓库托管在 Github 上，在 Linux 系统里安装了 git 工具可以将代码仓库克隆到本地运行。
# 克隆代码仓库
$ git clone https://github.com/tseemann/snippy
# 添加 snippy 以及依赖库的路径到环境变量中，使终端下直接运行 snippy 即可调用程序
$ echo 'export PATH="$HOME/repos/snippy/bin:$HOME/repos/snippy/binary/linux:$PATH"' >> ~/.bashrc
$ source ~/.bashrc
2. 使用sinppy¶
snippy 最基本的用法，可以将测序数据比对到参考基因组，获得相对应的 snps：
# --cpus 20 表示采用20个进程执行命令
$ snippy --cpus 20 --outdir S01 -ref reference.fa --R1 S01_R1_L001.fastq.gz --R2 S01_R2_L001.fastq.gz
对于多个样本的分析，寻找共同 snps，可以用 snippy-core 将各个样本的 snps 合并来实现。
# 依次将 mapping 到 reference.fa 的比对结果中的 snps 数据
$ snippy --cpus 20 --outdir S01 -ref reference.fa --R1 S01_R1_L001.fastq.gz --R2 S01_R2_L001.fastq.gz
$ snippy --cpus 20 --outdir S02 -ref reference.fa --R1 S02_R1_L001.fastq.gz --R2 S02_R2_L001.fastq.gz
$ snippy --cpus 20 --outdir S03 -ref reference.fa --R1 S03_R1_L001.fastq.gz --R2 S02_R2_L001.fastq.gz
# 将所有 snps 结果汇总，获得 core snps
$ snippy-core --prefix core-snps S01 S02 S03
也可以用 shell 的循环命令来调用 snippy 工具。
# 用 shell 脚本循环完成样本 snps 比对工作
$ for i in $(awk -F'_' '{print $1}' <(ls -D *.fastq.gz) | sort | uniq); \
> do snippy --cpus 20 --outdir $i -ref reference.fa --R1 $i*R1*.fastq.gz --R2 $i*R2*.fastq.gz; \
> done
$ snippy-core --prefix core-snps S*
snippy-core 生成的 core-snps.aln 文件，可以用进化树绘制工具生成进化树。
# 生成 clustal 格式的 core.aln 比对文件，这里用 splittree 构建进化树
$ splittree -i core-snps/core.aln
# 或者将生成的 .aln 格式比对文件转换成 .phy 格式，然后再用 raxml 构建 ML 树
$ raxml -f a -x 12345 -# 100 -p 12345 -m GTRGAMMA -s core.phy -n ex -T 40
注解
有很多中方法 可以将 aln 格式的文件转换成 phy 格式的比对文件。比如用 biopython
>>> from Bio import AlignIO
>>> AlignIO.convert(core-snps/core.aln, "clustal", core-snps/core.phy, "phylip-relaxed")
3. Reference¶
REALPhy¶
kSNP¶
SNP-Pipeline¶
第八章 Alignment 比对序列¶
Mummer¶
Mauve¶
1. 安装¶
2. 结果¶
- .alignment 文件是 mauve 将基因组比对后产生的xmfa格式的比对文件，支持xmfa的软件都可以查看或者进一步分析。
- .mauve 文件是 mauve 的共线性比对结果，只储存大片段比对情况，节约磁盘空间。具体的序列比对结果可以看.alignment 文件
- .island 文件是 mauve 比对后认为可能的一些基因岛，比如一些基因组携带但是部分基因组确实的片段序列，对于mauve里定义的岛，可以根据参数修改获得不同的“岛”
- .backbone 文件顾名思义，是mauve比对多个基因组序列后，生成的共有保守区序列起始和终止位置信息。每一行都是一个backbone区域，里面是该区域每个基因组的起始与终止区位置，有几个参比基因组就有多少个位置信息。
- guide tree 文件是 mauve生成的newick格式的模拟拓扑树文件。
- 其他有用的文件，还有snp文件，包含各个单核苷酸图突变及其位置信息；orthologs文件则将基因组的注释基因的orthologs信息提取保留；LCB文件则记录各个基因组共线性分析出来的区块边界信息。
附录内容¶
文件格式¶
1. Fastq 文件¶
1.1 Fastq 文件格式¶
~$ head -4 S1_R1.fastq
# 终端显示以下内容
@ST-E00294:36:hct2mccxx:4:1101:15281:1871 1:N:0:ATGCCGC
TCTGGCTCTGGCTCCTCTAATTCCTGCTGCTCCTCCTCTGACTGCCTCTGCTCCTCTGAGACACGGGGCTTCTCTTCTTCCTCCTCCACCTCAGGCTCCTTGGTTTCGGTCTCAGGACTATTGCTGCTGT
+
A<FFFFFKKKKKFKKKKKKKKKKKKAAKKFAKKFFF<KKKKKKKFA<AKKKKKFKKFKAKFFFKKKKKAFKKKKKFKFKAFAKKKKFKKFKKKKKFAKKKKKKKKKKFF7FKKKAKFFF,FKKKKKFKKK
这是一个从第一行 @ST-E00294:36:hct2mccxx:4:1101:15281:1871 1:N:0:ATGCCGC
获取数据¶
1. 安装 ascp connect¶
aspera connect 的`下载地址 <http://downloads.asperasoft.com/en/downloads/8?list>`_
# download ascp connect
$ wget http://download.asperasoft.com/download/sw/connect/3.6.2/aspera-connect-3.6.2.117442-linux-64.tar.gz
$ tar zxf aspera-connect-3.6.2.117442-linux-64.tar.gz
$ ./aspera-connect-3.6.2.117442-linux-64.sh
因为需要 ssh 密钥，软件只能以用户身份安装，不能以 root 安装。软件默认安装位置在 ~/.ascp/connect/
数据分析实践¶
分析所用数据目录结构
├── Sample_1
│ ├── 00_raw
│ │ └── *.fastq.gz
│ ├── 01_trimmed
│ ├── 02_assembly
│ ├── 03_clean
│ ├── 10_fasta
...
└── Sample_N
├── 00_raw
...
示例所使用数据
这里用一次 Miseq 测序一个 Run 的数据进行数据质控示例。本次 Run 一共对20个样品进行了 PE300 的测序。产生的测序数据所见如下：
$ ls
RN01_S01_L001_R1_001.fastq.gz RN06_S06_L001_R1_001.fastq.gz RN11_S11_L001_R1_001.fastq.gz RN16_S16_L001_R1_001.fastq.gz
RN01_S01_L001_R2_001.fastq.gz RN06_S06_L001_R2_001.fastq.gz RN11_S11_L001_R2_001.fastq.gz RN16_S16_L001_R2_001.fastq.gz
RN02_S02_L001_R1_001.fastq.gz RN07_S07_L001_R1_001.fastq.gz RN12_S12_L001_R1_001.fastq.gz RN17_S17_L001_R1_001.fastq.gz
RN02_S02_L001_R2_001.fastq.gz RN07_S07_L001_R2_001.fastq.gz RN12_S12_L001_R2_001.fastq.gz RN17_S17_L001_R2_001.fastq.gz
RN03_S03_L001_R1_001.fastq.gz RN08_S08_L001_R1_001.fastq.gz RN13_S13_L001_R1_001.fastq.gz RN18_S18_L001_R1_001.fastq.gz
RN03_S03_L001_R2_001.fastq.gz RN08_S08_L001_R2_001.fastq.gz RN13_S13_L001_R2_001.fastq.gz RN18_S18_L001_R2_001.fastq.gz
RN04_S04_L001_R1_001.fastq.gz RN09_S09_L001_R1_001.fastq.gz RN14_S14_L001_R1_001.fastq.gz RN19_S19_L001_R1_001.fastq.gz
RN04_S04_L001_R2_001.fastq.gz RN09_S09_L001_R2_001.fastq.gz RN14_S14_L001_R2_001.fastq.gz RN19_S19_L001_R2_001.fastq.gz
RN05_S05_L001_R1_001.fastq.gz RN10_S10_L001_R1_001.fastq.gz RN15_S15_L001_R1_001.fastq.gz RN20_S20_L001_R1_001.fastq.gz
RN05_S06_L001_R2_001.fastq.gz RN10_S10_L001_R2_001.fastq.gz RN15_S15_L001_R2_001.fastq.gz RN20_S20_L001_R2_001.fastq.gz
根据样品建立文件夹
首先我们对每一个样品的数据单独建立文件夹，并对 fastq 数据进行基本的质控。
# 新建以样本名为文件夹名，并移动数据到文件夹的 00_raw 目录下
$ for i in $(awk -F'L001' '{gsub("_$", "", $1); print $1}' < (ls *.fastq.gz) | sort | uniq ); \
> do mkdir -p $(basename $i) && mkdir -p $(basename $i)/00_raw; \
> mv $i*.fastq.gz $(basename $i)/00_raw/ ; done
# 查看建立的所有文件夹
$ ls -d
# 使用 FastQC 对测序数据进行质控分析
$ for i in $(ls -d RN*/00_raw/*.fastq.gz); do fastqc $i --extract -t 40 -q -o qc; done
$ python -m SimpleHTTPServer
打开浏览器，访问URL地址为服务器IP:8000，点击qc链接，里面有 fastqc 生成的各测序数据的 html 质控文件。我们以第一个样本 RN01_S01 为例，发现数据有以下问题：
- 250循环之后的数据质量不高（见图1）
- 序列有接头污染（见图2）
图1. 长读长区质量不高
图2. 接头污染
去除接头污染
首先我们要将污染的接头去除，采用 Trimmomatic 进行去除，同时也过滤低质量 reads 序列。
from fadapa import Fadapa
from Bio import SeqIO
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.Alphabet import generic_nucleotide
import argparse
parser = argparse.ArgumentParser()
parser.add_argument('-i', action=store, dest='qc_file', type=str \
help='input your fastqc data')
args = parser.parse_args()
qc_file = Fadapa(args.qc_file)
if qc_file.summary()[-3][0] == 'fail':
adaptors = []
ors = qc_file.clean_data("Overrepresented sequences")
ors.pop(0)
for (index, seq) in enumerate(ors):
adaptors.append(SeqRecord(Seq(seq[0], generic_nucleotide), id="adaptor_%d" % (index+1), description=""))
SeqIO.write(adaptors, "adaptor.fasta", "fasta")
print "Overrepresented sequences has been save to adaptors.fasta"
else:
print "No Overrepresented sequences"
上面的 python 脚本使用 Biopython 和 Fadapa 模块将 FastQC 生成的过表达序列保存正接头文件，让 Trimmomatic 进一步处理。
$ echo 'alias trimm="java -jar /opt/Trimmomatic-0.36/trimmomatic-0.36.jar"' >> ~/.bashrc
$ source ~/.bashrc
$ trimm PE -threads 40 -phred33 \
> 00_raw/RN01_S01_L001_R1_001.fastq.gz 00_raw/RN01_S01_L001_R2_001.fastq.gz \
> 01_trim/R1_trimmed.fastq.gz 01_trim/R1_unpaired.fastq.gz \
> 02_trim/R2_trimmed.fastq.gz 01_trim/R2_unpaired.fastq.gz \
> ILLUMINACLIP:01_trim/adaptors/adaptor.fasta:2:30:10 \
> LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:50
初步拼接
先用 spades 进行拼接，了解基因组情况。
$ cd RN01_S01/00_raw/
$ spades.py -k 127 -t 40 --careful -1 01_trim/R1_trimmed.fastq.gz -2 01_trim/R2_trimmed.fastq.gz -o assembly/spades
$ du -h assembly/spades/scaffolds.fasta
10.1M assembly/spades/scaffolds.fasta
$ cat assembly/spades/scaffolds.fasta | grep '>' | tail -1
>NODE_3360_....
$ cat assembly/spades/scaffolds.fasta | grep '>' | awk -F'_' '{if ($6<10) print $0}' | wc -l
3104
结果获得的 scaffolds.fasta 文件大小为10M左右，而我们测序的目的物种基因组大小仅为3M；里面的 contigs 数量达到3360个，并且 contigs 平均覆盖度小于10的有3100多个，说明原始的测序数据很可能被其他物种污染了。这种污染可能发生在核酸提取，文库制备或者测序中（清洗管路不彻底）
观察覆盖度与污染序列的关系
用 Blast 的方法来筛选污染序列
# 下载污染物种的基因组数据，进行序列比对，看组装的nodes里那些是来源于污染物种。如果有多个污染物种，则可以将基因组数据合并 `cat 1.fa 2.fa 3.fa > containment.fasta`
$ makeblastdb -db containment.fasta -parse_seqids -db_type nucl
$ blastn -db containment.fasta -query scaffolds.fasta -max_hsps 1 -outfmt 6 -out result
# blast 结果的相似性筛选，小于90认为与污染物种不同。将过滤的片段长度求和，判断过滤的片段是否与物种基因组大小一致。如果接近，那么即使有部分片段遗漏，但是大部分基因组数据已经保留。
$ awk '{ if ($3 < 90) print $1 }' result > filter_nodes
$ awk -F'_' 'BEGIN {len=0} {len+=$4} END {print len}' filter_nodes
# 进一步用目的物种的参考基因组进行blast，以确保没有其他物种污染。
$ makeblastdb -db reference.fasta -parse_seqids -db_type nucl
$ blastn -db reference.fasta -query assembly.fasta -max_hsps 1 -outfmt 6
用 Mapping 目的基因组来筛选 reads 再进行拼接
抓取目标nodes
# 保存代码到 get_nodes.py 文件中，运行 python get_nodes.py
from Bio import SeqIO
input_file = 'scaffolds.fasta'
filter_file = 'filter_nodes'
output_file = 'assembly.fasta'
wanted = set(line.rstrip("\n").split(None,1)[0] for line in open(filter_file))
print "Found %i unique identifiers in %s" % (len(wanted), filter_file)
records = (r for r in SeqIO.parse(input_file, "fasta") if r.id in wanted)
count = SeqIO.write(records, output_file, "fasta")
print "Saved %i records from %s to %s" % (count, input_file, output_file)
if count < len(wanted):
print "Warning %i IDs not found in %s" % (len(wanted)-count, input_file)
玩转 edirect¶
NCBI 的 Entrez 工具集大家很熟悉了，特别是基于 web 的 E-utilis，是日常检索 NCBI 数据库的最主要应用。而 edirect 是 NCBI 开发的用于 linux 命令行界面中的快速检索和下载工具。
NCBI 的 Entrez 工具平时大家都经常使用，既可以从 web 页面访问 NCBI 来查询，也可以利用 Entrez 提供的 Web Services 来实现诸多功能（如自动查询）。此外 NCBI 还提供了 Entrez 的命令行工具 edirect。edirect 全称时 EntrezDirect，使用它可以很方便的在服务器上进行 Entrez 的检索和抓取。
1. 下载安装¶
~/tmp$ wget ftp://ftp.ncbi.nlm.nih.gov/entrez/entrezdirect/edirect.tar.gz
~/tmp$ tar zxvf edirect.tar.gz -C ~/app
~/tmp$ sh ~/app/edirect/setup.sh
在看到：
> ENTREZ DIRECT HAS BEEN SUCCESSFULLY INSTALLED AND CONFIGURED
字符后，安装完成。可以不运行`setup.sh`直接使用edirect的一些工具，不过如果要使用全部功能，有些perl模块需要安装，所以运行`setup.sh`来完成CPAN的modual安装以及 go 版本的xtract的编译（原来的perl版本速度不够快）。当看到`ENTREZ DIRECT HAS BEEN SUCCESSFULLY INSTALLED AND CONFIGURED`，表明完成安装。
如果是Ubuntu操作系统，由于uname位于/bin，而xtract接口调用脚本使用/usr/bin/uname来判断使用那个版本的xtract，因此要正常使用xtract需要添加一个link。
~$ sudo ln -s /bin/uname /usr/bin/uname
也可以通过 perl 来安装：
$ cd ~
$ perl -MNet::FTP -e \
> '$ftp = new Net::FTP("ftp.`NCBI`_.nlm.nih.gov", Passive => 1); $ftp->login; \
> $ftp->binary; $ftp->get("/entrez/entrezdirect/edirect.zip");'
$ unzip -u -q edirect.zip
$ rm edirect.zip
$ export PATH=$PATH:$HOME/edirect
$ ./edirect/setup.sh
注解
NCBI 从2016年11月起，网站从http迁移到https协议，因此老的版本edirect会无法使用。
- 工具使用
2.2 数据库列表¶
主要工具列表
- 具体例子
查看2012年至今发表的霍乱弧菌CTX相关文献摘要
~$ esearch -db pubmed -mindate 2012 -maxdate 2016 -datetype PDAT -query "vibrio cholerae[CTX]" | efetch -format abstract > abstract.txt
查看2016年CDC发布的所有用 miseq 测序的沙门菌文库制备方法
~$ esearch -db sra -query "salmonella miseq CDC" -mindate 2016 -maxdate 2016 -datetype PDAT | efetch -format runinfo | cut -d ',' -f 12 > library.txt
绘制单增李斯特菌 hlyA 蛋白质进化树图
~$ esearch -db protein -query "(listeria monocytogenes hlyA) NOT partial" | efetch -db protein -format fasta > hlyA.fasta
~$ raxml -p 12345 -x
对10年内发表的 hlyA 相关文献的作者排序
~$ esearch -db pubmed -mindate 2006 -maxdate 2016 -datetype PDAT -query "hlyA"
查看 taxonomy
~$ esearch -db taxonomy -query "lis* OR sal*"
下载所有甲型副伤寒沙门菌基因组序列
~$ esearch -db nuccore -query "Salmonella[porgn:__txid54388] complete genome" | efetch -db nuccore -format fasta > SPA.fasta
Reference: