Archives

CentOS 7运维-Windows平台如何使用SSH登录Linux服务器

本文介绍Windows平台下两款登录Linux终端的应用。

一、Putty: 最经典知名的免费SSH客户端

下载地址: 64位 / 32位

经过简单配置后可以工作了。

image.png-91.3kB

通过1-7步骤,可以登录到远程服务器端,并使用SSH在终端完成各种操作:
下面解读下七个步骤:

1.设置远程服务器IP地址,确保可以访问;
2.设置远程服务器IP端口,默认22, SSH访问;
3.设置会话名称;
4.保存下次可以直接会话列表进行选择;
5.打开终端,当然可以直接双击会话列表中记录登录;
6.输入用户名
7.输入密码

2、MobaXterm: 全能的远程终端登录软件

MobaXterm功能强大、界面漂亮,可以满足大部分的Windows和Linux之间的交互需求,比如SSH登录,SFTP访问等,当然还有一些好玩的小游戏。

MobaXterm界面

2.1 工具下载

选择Home免费版本 MobaXterm Home Edition
Windows 平台如何安装跳过。

2.2 推荐理由

  • 构建本地开发平台
    可以很方便安装 GCC以及各种库文件,方便Windows开发编译各种软件。
  • 内置 X server, 可以执行图形界面应用。
  • 多标签,分屏
  • SFTP文件传输
  • 直接支持UTF-8编码,不乱码
  • 保留用户名和密码
  • 一键操作所有活跃终端

2.3 SSH远程登录

Mobaxterm

通过1-5步骤,可以登录到远程服务器端,并使用SSH在终端完成各种操作:
下面解读下五个步骤:

1.设置新会话;
2.选择SSH登录;
3.设置远程服务器IP地址,确保可以访问;
4.设置远程服务器用户名;
5.设置远程服务器IP端口,默认22;

第一次登录会提示输入用户密码,并会提示是否保存密码,选择后下次自动登录。

blast-utils之表格文件处理:blast-utils小工具集合

BLAST类序列比对工具在生物学数据分析中应用最为广泛,比如16S全长序列比对(MEGABLAST), 氨基酸功能注释(BLASTP)等, 另外还有各种快速比对工具,比如: DIAMONDGHOSTXRAPSearch2等,这些工具默认都支持制表符分隔的输出文件(BLAST程序指定的 -outfmt 输出格式)。

格式如下

中文描述英文描述
1查询序列表标识符Query
2目标序列标识符Target
3相似度identity
4比对长度Alignment length
5错配数Number of mismatches
6gap数Number of gap opens
7查询序列比对起始位置Start position in query
8查询序列比对终止位置End position in query
9目标序列比对起始位置Start position in target
10目标序列比对终止位置End position in target
11E值E-value
12Bit score值Bit score

BLAST比对结果通常一条查询序列可以分几个片段和目标序列比对上,其中每个片段的比对结果都是一个HSP, 比对到的目标序列称之为一个hit, 比对结果可以包含多个hit, 我们对序列注释有时候采用的策略是 best hit 或者 best_hsp, 这样就需要对比对结果进行各种处理。

blast_utilsgithub 仓库) 小工具集合, 主要使用场景包括:

  1. 根据E值、Bit score值、相似度筛选所有满足条件的结果, blast_hits
  2. 根据E值、Bit score值、相似度筛选满足条件的最佳HSP, blast_hsp
  3. 选取每条查询序列的Top Hsp, best_hsp
  4. 对目标序列进行注释, blast_annotation

下面对每个程序使用进行介绍:

1. blast_hits

命令行接口

$ blast_hits
Usage: blast_hits  [options]  <blast|->
Options:
  -b DOUBLE  MIN bit score, default: [60]
  -e DOUBLE  MAX E-value, default: [0.001]
  -i DOUBLE  MIN identity, default: [0]
  -v print version number

根据上述描述, blast_hits 接受三个可选参数, bit score 、E-value 和 identity,就是输出所有满足条件的比对结果。

实例:

blast_hits -b 60  -e 0.001 -i 60  S7711.tsv

2. blast_hsp

命令行接口

$ blast_hsp
Usage: blast_hits  [options]  <blast|->
Options:
  -b DOUBLE  MIN bit score, default: [60]
  -e DOUBLE  MAX E-value, default: [0.001]
  -i DOUBLE  MIN identity, default: [0]
  -v print version number

根据上述描述, blast_hsp 接受三个可选参数, bit score 、E-value 和 identity,就是输出满足条件的第一个hsp。

实例:

blast_hsp -b 60  -e 0.001 -i 60  S7711.tsv

3. best_hsp

命令行接口

$best_hsp
Usage: best_hsp <blast|->
version: 0.0.2

$best_hsp 接受一个参数, 最输入文件不进行限制,只输出第一个hsp。

实例:

best_hsp S7711.tsv

4. blast_annotation

因为我们有时需要对目标序列进行各种关联,需要使用目标序列的标识符信息,所以:

LDEE01000001.1.2    sp|P33330|SERC_YEAST

可能会比:

LDEE01000001.1.2        sp|P33330|SERC_YEAST Phosphoserine aminotransferase ...

更方便;

blast_annotation 做的事情就是对sp|P33330|SERC_YEAST 进行注释,变成 sp|P33330|SERC_YEAST Phosphoserine aminotransferase ... 模式。

程序也比较简单, 需要提供序列->注释的映射文件。

实例:

blast_annotation  annotation.txt  S7711.tsv  >S7711_ann.tsv

本文材料为 BASE (Biostack Applied bioinformatic SEies ) 课程 Linux Command Line Tools for Life Scientists 材料, 版权归上海逻捷信息科技有限公司** 所有

Last update:2017-11-15 7:07 PM

庖丁解牛,微生物基因组重测序:一条命令完成从sra到vcf

一、引子

前面在 Twitter上看到一个针对细菌基因组重测序流程(鉴定SNP)的使用工具调查,这里我们就庖丁解牛般解读下这些流程工具是如何工作的,流程又是如何设计的。

SNP Pipeline

该调查数据来自 google.docs

从图中可以看出snippy还是比较受欢迎的,Torsten Seemann的流程一向比较稳定,接口一向比较简单。

这次解读的主角是snippy,使用起来确实很方便,比如:

fastq-dump --defline-seq '@$sn/$ri' --defline-qual '+' ERR114970.sra -Z --split-3 |\
seqtk trimfq - | \
snippy --cpus 16 --outdir ERR114970 --ref GCF_000069185.1_ASM6918v1_genomic.gbff  \
--prefix ERR114970 -peil /dev/stdin &>log.txt

一条命令可以解决我们从sra到vcf到的问题,当然其中使用的工具很多了,但是都被封装在脚本里了。

二、流程细节

2.1. sra生成interleave格式fastq文件

fastq-dump --defline-seq '@$sn/$ri' --defline-qual '+' ERR114970.sra -Z --split-3

该结果使用-Z 选项将结果输出至标准输出流,格式是Interleave格式,这种模式有很多好处, 比如质量值修剪,序列比对等,后续很多程序都支持,这样比较方便以steaming方式工作。

内部链接: ENA和SRA数据预处理: SRA Toolkit 和 seqtk_sra

2.2 去除低质量序列

seqtk trimfq - 

steaming 模式切除低质量序列。

内部链接: NGS数据质量过滤:Usearch fastq-filter、seqtk trimfq和seqtk_filter

2.3. 序列比对

至此,snippy 预处理工作已经做完,现在可以进入snippy的工作模式,这里只记录主要步骤。

gbk转gff和fasta:

使用`BioPerl`从GBK文件中提取fasta和gff文件用于构建序列比对索引和snpEff库。

构建基因组文件序列

samtools faidx reference/ref.fa

构建基因组索引:

bwa index reference/ref.fa

序列比对:

bwa mem  -v 2 -M -R '@RG    ID:ERR114970    SM:ERR114970' \
    -p -t 16 reference/ref.fa /dev/stdin | \
    samtools view -@ 16 -F 3844 -q 60 -S -b -u -T reference/ref.fa - |\
    samtools sort -O bam -o ERR114970.bam -@ 16 -

构建bam文件索引:

samtools index ERR114970.bam

计算比对深度:

samtools depth -q 20 ERR114970.bam | bgzip > ERR114970.depth.gz

这一步完成了序列比对工作,使用的都是比较常规的操作。

2.4. Variation calling

freebayes-parallel reference/ref.txt 16 -p 1 -q 20 -m 60 \
     --min-coverage 10 -V -f \
     reference/ref.fa ERR114970.bam > ERR114970.raw.vcf

使用freebayes (Bayesian haplotype-based polymorphism discovery)鉴定 Variation, Freebayes的单倍体模型比较适合细菌基因组SNP/INDEL鉴定。

-q –min-base-quality, 最低碱基质量; -q 20; -m –min-mapping-quality,最低比对质量; -m 60; -p –ploidy 设置单倍体模型, -p 1; –min-coverage: 最小覆盖度。 –min-coverage 10;

snippy使用的是freebayes并行模式 freebayes-parallel

需要使用分割regions文件, 使用fasta_generate_regions.py 生成。

fasta_generate_regions.py reference/ref.fa.fai 80866 > reference/ref.txt

Variation 过滤:

snippy-vcf_filter --minqual 10 --mincov 10 --minfrac 0.9  ERR114970.raw.vcf > ERR114970.filt.vcf

使用下面参数控制variation 质量:

–mincov :位点最低覆盖度, 10 –minfrac :alt碱基最低比率,0.9

2.5. SNP/INDE 注释

创建配置文件和snpEff build库文件,然后使用snpEff

snpEff build -c reference/snpeff.config -dataDir . -gff3 ref
snpEff ann -no-downstream -no-upstream -no-intergenic -no-utr -c reference/snpeff.config -dataDir . -noStats ref ERR114970.filt.vcf > ERR114970.vcf

使用-no-downstream -no-upstream -no-intergenic -no-utr 过滤掉基因上游、下游、间区以及UTR区域的变异位点, 生成注释的VCF文件。

至此我们已经可以获得从SRA文件到VCF的细节过程。

2.6. reference guided assembly

snippy 提供了另外一个我们称之为”reference guided”的拼装,实现方式为:

vcf-consensus ERR114970.vcf.gz < reference/ref.fa > ERR114970.consensus.fa

三、总结

3.1. 设计程序时一定要注意区分:

标准输出流: stdout 标准错误流: stderr

准则: 如果不是结果文件,一律使用标准错误流打印。

3.2 Pipeline 设计

很多工作流,其实细节都很好实现,工作流主要起的作用就是封装,让接口更加简单, snippy 集成了 bioperl、bgzip、bwa、freebayes、samtools、tabix、vcfallelicprimitives、vcfstreamsort、vcfuniq等外部程序,使用这些程序可以使用相对路径调用,也可以设置PATH路径,但是最佳实践是环境变量只有程序能够访问,这样就不会造成程序版本污染。

本文材料为 BASE (Biostack Applied bioinformatic SEies ) 课程 Linux Command Line Tools for Life Scientists 材料, 版权归上海逻捷信息科技有限公司** 所有

Last update:2017-11-14 7:46 PM

使用U盘安装CentOS7 Linux 操作系统

现在很多服务器都已经放弃了光驱设备,使用U盘安装操作系统成为比较方便的方案,本文介绍如何烧制可引导的U盘系统。

1. 下载ISO镜像

CentOS提供的所有的镜像列表: List of CentOS Mirrors
我们选择:Sohu Inc, Beijing P.R. China

推荐下载 DVD版本: CentOS-7-x86_64-DVD-1708.iso

2. 准备U盘

准备一个8Gb大小的U盘。

3. 安装

安装 win32diskimager-1.0.0, Windows平台下如何安装不需多说。

4. 烧制

win32diskimager

打开程序,界面很简单,插入U盘,选中U盘符。

然后直接点击 Writer 确认后完成烧制, USB 2.0的接口估计需要15分钟左右才能烧制完成。

注意: 确保U盘里没有重要文件。

5. 系统安装

一般需要设置下, 让sBIOS系统引导的优先级,选中U盘即可,剩下的和光盘安装一样了。

16S扩增子组成谱预测功能谱:本地快速执行otutable_norm和otutable_melt以及otutable_kann

当前可以使用16S扩增子组成谱预测功能谱的工具主要有:PICRUStTax4funPiphilin, 本文主要使用PICRUSt的预测文件,进行重新实现。

PICRUSt 采用了一种ASR(ancestral state reconstruction, 可参考1985年的一篇文章:Phylogenies and the Comparative Method)方法重构物种的16S拷贝数和基因组的基因功能信息(使用KEGG、COG描述),基于系统进化的逻辑是同一个clade分支下的基因组(物种)具有相似的形状,这样就可以使用已测序基因组预测未测序基因组的功能空间,主要依赖参数为系统进化树分支距离。

PICRUSt

图来自PICRUSt文章

原理图分拆解读:

1. 使用Greengene 16S序列构建系统进化树,去除未测序基因组;
2. 预测所有祖先状态的性状(16S拷贝数、基因功能);
3. 预测未知基因组节点的性状(16S拷贝数、基因功能);
4. 标准化,使用预测的拷贝数对OTU丰度进行标准化;
5. 基因功能空间预测(KO、COG)丰度表;

其实关键算法在1、2、3部分。 这篇文章我们直接使用PICRUSt根据 Greengene 13.5 进行的预测结果进行一些数据分析。

需要的文件为双端合并后的序列数据(Usearch sample=’xxx’ 注释文件),我们会在下面一步一步演示如何操作。 准备工作:

下载PICRUSt预计算拷贝数和KO表: 16S_13_5_precalculated.tab.gzko_13_5_precalculated.tab.gz 下载Greengene 97OTU代表序列: 97_otus.fasta.gz

如何使用这些程序,请移步: 如何使用biostack 提供的命令行工具

1. 构建OTU表

虽然有文章讨论过使用closed- 和 open-reference 构建 OTUs 有多不准确,如下:

Accuracy of microbial community diversity estimated by closed- and open-reference OTUs.
Edgar RC.
PeerJ. 2017 Oct 4;5:e3889. doi: 10.7717/peerj.3889. eCollection 2017.
PMID: 29018622

如果有兴趣看看问题在哪里, 可以移步:closed_ref_problems, 这不是这篇文章的主题,因为我们使用基于Greengene 13.5 的PICRUSt预先计算的结果,所以只能选择使用Greengene作为参考构建OTU表。

这里使用Usearchclosed_ref 构建closed OTU表;

closed_ref 命令行接口:

usearch -closed_ref reads.fastq -db gg97.fa -otutabout otutab.txt -strand plus -tabbedout closed.txt

下面简单解读:

reads.fastq: 为具有 'sample=xxx' 注释的16S扩增子序列集合。
-db        : 为greengene 97%聚类 OTU代表序列。
-strand    : plus或者both,搜索查询序列的条件,只都搜索正链还是正负链都搜索。
-otutabout : OTU表输出;
-tabbedout : closed_ref 比对report文件;
-id        : 最低相似度要求,默认参数97%;

下载Greengene 97% OTU代表序列文件:

wget https://github.com/biocore/qiime-default-reference/blob/master/qiime_default_reference/gg_13_8_otus/rep_set/97_otus.fasta.gz

执行构建OTU表:

gunzip 97_otus.fasta.gz
usearch -closed_ref reads.fastq -db 97_otus.fasta -otutabout otutab.txt \
        -strand plus -tabbedout closed.txt -threads 20

2. 16S OTU表 标准化

上一步我们已经获得了OTU表信息,下一步就是使用Greengene的16S 拷贝数进行标准化OTU表; 我们使用 `otutable_norm` 程序实现标准化工作

命令行接口:

$ otutable_norm
Usage: otutable_norm <copy_number> <otutab>
version: 0.0.1

otutable_norm 接受两个参数: 拷贝数表和OTU表;

执行:

wget http://kronos.pharmacology.dal.ca/public_files/picrust/picrust_precalculated_v1.1.2/13_5/16S_13_5_precalculated.tab.gz
otutable_norm  16S_13_5_precalculated.tab.gz  otutab.txt >otutab_norm.txt

3. 合成功能丰度表

这一步使用 otutable_melt 完成合成功能谱的工作。 命令行接口:

$otutable_melt
Usage: otutable_melt [options] <feature-matrix> <otu_table>
Options:
  -r  round value to nearest intege;
  -v  print version number

otutable_melt 接受两个参数:PICRUSt兼容的特征表和OTU表(QIIME兼容的OTU表),这里使用上一步产生标准化后的OTU表, otutable_melt也接受 Bugbasedefault_traits_precalculated.txt数据;

-r: 四舍五入,转换成整型输出。

执行:

wget http://kronos.pharmacology.dal.ca/public_files/picrust/picrust_precalculated_v1.1.2/13_5/ko_13_5_precalculated.tab.gz
otutable_melt   ko_13_5_precalculated.tab.gz otutab_norm.txt  >ko_profiling.txt

或者:

otutable_melt   <(pigz -dk -c  -p 20  ko_13_5_precalculated.tab.gz)  otutab_norm.txt  >ko_profiling.txt

现在我们已经可以获取KEGG的入门券了,后面我们可以使用KO关联KEGG的很多重要信息, 比如 Pathway,Module,Enzyme等。

4. otutable_kann

KO(KEGG Ortholog)是KEGG唯一一个通过序列和功能关联的中间点,我们获得KO的丰度信息, 可以通过合并同一个Pathway或Module的KO丰度,计算Pathway和Module的丰度。

实现该功能我们使用 otutable_kann实现;

命令行接口:

$ otutable_kann
Usage: otutable_kann [options] <feature-map> <otu_table>
Options:
  -r  round value to nearest intege;
  -v  print version number

otutable_kann接受两个参数, KEGG的KO映射文件和OTU表文件(KEGG或COG的丰度表这里统一为otu_table表)。 文件为两列:

第一列: KO Number 第二列: KO 映射表;

该文件可以通过KEGG的API服务下载,比如:

curl http://rest.kegg.jp/link/pathway/ko | \
sed 's/ko://; s/path://' | grep "map"  >ko_pathway.tab
curl http://rest.kegg.jp/link/module/ko |   \
sed 's/ko://; s/md://'  >ko_module.tab

下载注释文件:

curl http://rest.kegg.jp/list/module | sed 's/md://' >module-definition
curl http://rest.kegg.jp/list/pathway | sed 's/path://' >pathway-definition
curl http://rest.kegg.jp/list/ko | sed 's/path://' >ko-definition

执行(以Module为例子):

otutable_kann   ko_module.tab  ko_profiling.txt >module_profiling.txt

如果需要注释信息需要使用 tabtk_definition 添加:

otutable_kann   ko_module.tab  ko_profiling.txt  | \
tabtk_definition -d " " module-definition  - >module_profiling.txt

5. 总结

到这里我们可以很简单高效完成共物种组成谱到功能谱的预测,我们提供的otutable_normotutable_meltotutable_kann 都比较通用,组合起来可以解决很多问题。

本文材料为 BASE (Biostack Applied bioinformatic SEies ) 课程 Linux Command Line Tools for Life Scientists 材料, 版权归上海逻捷信息科技有限公司** 所有

Last update:2017-11-14 7:46 PM

如何使用Biostack提供的工具

我们分享的工具集都在 Github (地址) 上, 目前只提供编译后的版本, 如下图:

ScreenClip.png-13.6kB

Biostack Github 仓库

如何下载使用工具,目前所有的工具都是在 CentOS 7编译,可能不兼容版本比较低的系统。

1. 下载工具

使用 git 直接克隆对应仓库, 比如:

git clone https://github.com/jameslz/blast_utils

如果提示没有安装 git, 请先安装git:

sudo yum install git

2. 添加环境变量

如果我们将 blast_utils 克隆到了目录: /project/tools/blast_utils, 只需要执行:

export PATH=/project/tools/blast_utils:$PATH

就可以使用这些程序了, 或者执行下面命令添加到bashrc文件,每次启动bash都会自动加载该目录:

echo export PATH=/project/tools/blast_utils:$PATH  >> ~/.bashrc
source ~/.bashrc

如果提示权限不够,请执行:

chmod -R 775  /project/tools/blast_utils

3. 命令执行

随意切换到任何目录,我们可以直接使用 blast_utils 目录下应用程序,比如:

$ cd

$ blast_hsp
Usage: blast_hits  [options]  <blast|->
Options:
  -b DOUBLE  MIN bit score, default: [60]
  -e DOUBLE  MAX E-value, default: [0.001]
  -i DOUBLE  MIN identity, default: [0]
  -v print version number

4. 其它问题

如果可执行文件在仓库的 bin 目录, 添加目录时请包括bin目录, 比如:

 export PATH=/project/tools/blast_utils/bin:$PATH

本文材料为 BASE (Biostack Applied bioinformatic SEies ) 课程 Linux Command Line Tools for Life Scientists 材料, 版权归上海逻捷信息科技有限公司所有

Last update:2017-11-15 3:02 PM

fasta/q序列抽取和排序: seqtk_reoder

一、使用场景

给定fasta/q文件和一个序列列表,按照给定的序列顺序抽取序列,如不存在则忽略,先前的 seqtk_subseq( 如何根据列表快速提取fasta/q序列子集/补集:seqtk_subseq ) 抽取出序列按照标识符在序列文件中出现的顺序输出,有时我们需要按照给定的标识符输出序列。

seqtk_reoder使用reoder作为名字输出就是指序列输出顺序按照指定顺序输出。

二、工具介绍

通过seqtk_info 查看测试序列信息:

$seqtk_info  shear.fastq
#sequence       base    min_len max_len avg_len
47525   23905075        503     503     503.00

一共47525序列;

2.1 seqtk_split

命令行接口:

$ seqtk_reorder
Usage: seqtk_reorder <fasta/q> <list>
version: 0.0.1

seqtk_reorder 接受两个参数, 序列文件列表文件, 支持stream流模式;

实例

$cat  shear.idx | seqtk_reorder  shear.fasta -

$cat  shear.idx | seqtk_reorder  shear.fastq.gz  -

$cat  shear.idx | seqtk_reorder  <(gunzip  -c shear.fastq.gz)  -

该程序设计使用hash map实现存储目标序列,所以不支持数据库比较大的序列抽取, seqtk_subseq 只存储列表,所以适用范围更广。

本文材料为 BASE (Biostack Applied bioinformatic SEies ) 课程 Linux Command Line Tools for Life Scientists 材料, 版权归上海逻捷信息科技有限公司** 所有

Last update:2017-11-14 4:08 PM

fasta/q序列拆分三剑客: seqtk_split、seqtk_partition 和 seqtk_shred

一、使用场景

本部分内容主要解决以下几个问题:

  1. 对指定fasta/q序列文件拆分成每份N条序列问题,比如执行InterProScan。
  2. 对指定fasta/q序列文件拆分成指定N份问题,比如我们有40核心,程序不支持多线程(或者支持的不理想),那拆分40份多进程执行, 比如HMMER
  3. 对指定fasta/q序列文件序列拆分,每条序列保留最大l长度,并保证分割非序列存在N个碱基交叠。使用的场景多,比如 454 Newbler拼装软件,不支持很长度序列,如果我们对不同拼装版本的组装结果混合组装,将序列撕碎(保留比较大的overlap)再进行组装,就没有问题了。

二、工具介绍

以下三个程序分别解决以上三个问题。

通过seqtk_info 查看测试序列信息:

$seqtk_info  shear.fastq
#sequence       base    min_len max_len avg_len
47525   23905075        503     503     503.00

一共47525序列;

2.1 seqtk_split

命令行接口:

$ seqtk_split
Usage: seqtk_split [options] <fasta/q> <size> <prefix>
Options:
  -t  fastx file suffix: default: [fasta];
  -v  print version number

seqtk_split 接受三个参数, 序列文件序列大小(条数)输出文件前缀
可以通过 -t 指定输入文件类型, 比如 fasta或者fastq

实例

$mkdir split

$seqtk_split -t fastq  shear.fastq   10000   split/split

$ls split
split_1.fastq  split_2.fastq  split_3.fastq  split_4.fastq  split_5.fastq

seqtk_split会返回分拆的文件个数,方便自动化批量处理。

2.2 seqtk_partition

命令行接口:

$ seqtk_partition
Usage: seqtk_partition [options] <fasta/q> <partition: 10> <prefix>
Options:
  -t  fastx file suffix: default: [fasta];
  -v  print version number

seqtk_partition 接受三个参数, 序列文件分割数输出文件前缀
可以通过 -t 指定输入文件类型, 比如 fasta或者fastq

实例

$mkdir partition

$seqtk_partition -t fastq  shear.fastq   3   partition/partition

$ls partition
partition_1.fastq  partition_2.fastq  partition_3.fastq
2.3 seqtk_shred
$ seqtk_shred
Usage: seqtk_segment <fasta> <length> <overlap>

seqtk_shred 只支持fasta文件, 接受三个参数, 序列文件分割成序列长度序列交叠长度

实例

$seqtk seq -A  shear.fastq  >shear.fasta

$seqtk_shred shear.fasta 300 200 >shear.shred.fasta

输出序列会根据位置添加_1, _2, _3…

本文材料为 BASE ( Biostack Applied bioinformatic SEies ) 课程 Linux Command Line Tools for Life Scientists 材料, 版权归上海逻捷信息科技有限公司所有。

Last update:2017-11-14 3:11 PM

使用Usearch进行扩增子数据分析:The Full Story

准备工作:

添加辅助程序路径至环境变量:

export PATH=$PATH:/biostack/tools/pipelines/div_seq-0.2.3/bin/utils

export PATH=$PATH:/biostack/labs/seqtk_utils/release

export PATH=$PATH:/biostack/labs/tabtk_utils/release

export PATH=$PATH:/biostack/tools/microbiome/KronaTools-2.7/bin

准备好原始数据,样品信息表、以及引物信息。 软件需要安装:R、seqtk、usearch、seqtk_utils、tabtk_utils、biom以及fasttree、 Trimal, mafft、KronaTools以及phylommand

1. 数据整体质量评估

参考: 获取fasta/q文件的碱基组成和长度信息: seqtk_info 和 seqtk_fqchk

mkdir quality
seqtk_fqchk  data/raw.R1.fastq.gz  >quality/raw_R1.tsv
fqchk_base.R  quality/raw_R1.tsv   quality/raw_base_R1.pdf R1
fqchk_qual.R  quality/raw_R1.tsv   quality/raw_qual_R1.pdf R1

可以生成质量和碱基分布图:

raw_qual_R1.png-298.2kB
raw_base_R1.png-711.2kB

2. 数据拆分

参考: 快速对扩增子序列进行数据拆分: seqtk_demultiplex

mkdir demultiplex
cut -f1,2 mapping_file.txt | sed  's/,/\t/'  >demultiplex/barcode.txt
seqtk_demultiplex -b  demultiplex/barcode.txt  -1  data/raw.R1.fastq.gz   -2  data/raw.R2.fastq.gz  -d demultiplex -l 5

3. 双端序列合并

参考: 双端序列合并:fastq_mergepairs

单个样本可以这样:

mkdir mergepairs
usearch -fastq_mergepairs demultiplex/D1_R1.fastq \
        -reverse demultiplex/D1_R2.fastq      \
        -fastqout mergepairs/D1.fastq         \
        -fastqout_notmerged_fwd mergepairs/D1.notmerged_fwd.fastq \
        -fastqout_notmerged_rev  mergepairs/D1.notmerged_rev.fastq \
        -log  mergepairs/D1.log \
        -threads 24 \
        -fastq_minmergelen 100 \
        -fastq_maxmergelen 500 \
        -fastq_maxdiffs 10  \
        -fastq_pctid 80 \
        -fastq_trunctail 2

可使用Bash进行循环迭代:

cut -f1 mapping_file.txt  | grep -v "#" | perl -ane 'print qq{usearch -fastq_mergepairs demultiplex/$F[0]\_R1.fastq -reverse demultiplex/$F[0]\_R2.fastq -fastqout mergepairs/$F[0].fastq -fastqout_notmerged_fwd mergepairs/$F[0].notmerged_fwd.fastq  -fastqout_notmerged_rev  mergepairs/$F[0].notmerged_rev.fastq -log  mergepairs/$F[0].log -threads 24 -fastq_minmergelen 100 -fastq_maxmergelen 500  -fastq_maxdiffs 10 -fastq_pctid 80 -fastq_trunctail 2 \n}' |bash

4. 引物识别和切除引物

参考: 对扩增子序列进行引物匹配: usearch 和 primersearch 切除引物序列获取扩增目标序列:seqtk_pcr_strip

mkdir primer_match pcr_strip
usearch -search_pcr mergepairs/D1.fastq \
        -db  primers.txt -strand both   \
        -threads 5 -maxdiffs 2          \
        --pcrout primer_match/D1.txt -log primer_match/D1.log ;

seqtk_pcr_strip -l 4  primer_match/D1.txt  mergepairs/D1.fastq  >pcr_strip/D1.fastq

通过bash可以批量执行:

cut -f1 mapping_file.txt  | grep -v "#" | perl -ane 'print qq{usearch -search_pcr mergepairs/$F[0].fastq -db  primers.txt -strand both -threads 5 -maxdiffs 2  --pcrout primer_match/$F[0].txt -log primer_match/$F[0].log ; \n}' |bash

cut -f1 mapping_file.txt  | grep -v "#" | perl -ane 'print qq{seqtk_pcr_strip -l 4 primer_match/$F[0].txt  mergepairs/$F[0].fastq  >pcr_strip/$F[0].fastq ; \n}' |bash

5. 添加序列注释信息和pool样本

参考: 添加和去除 fasta/q 序列标识符注释信息: seqtk_labels和fastx_strip_annots

mkdir pools
cut -f1 mapping_file.txt  | grep -v "#" | perl -ane 'print qq{seqtk rename  pcr_strip/$F[0].fastq  $F[0]_ | seqtk_labels -  sample=$F[0]  >> pools/seq.fq ; \n}' |bash

此时序列标识符为:

@D1_1;sample=D1;

6. 质量控制

参考: NGS数据质量过滤:Usearch fastq-filter、seqtk trimfq和seqtk_filter

usearch -fastq_filter pools/seq.fq   \
        -fastq_maxee 1.0    \
        -fastaout  pools/filtered.fa -relabel Filt

7. fastx_uniq 去冗余

参考: 扩增子序列样本序列丰度谱: fastx_uniques、fastx_uniques_persample、 DADA2 和 seqtk_seqtab

usearch -fastx_uniques  pools/filtered.fa \
        -relabel Uniq                 \
        -sizeout -fastaout  pools/derep.fa \
        -uc pools/derep.uc      \
        -log pools/derep.log; 

8. UPARSE

参考: 扩增子数据分析之构建OTU:UPARSE

mkdir uparse
usearch -cluster_otus pools/derep.fa -otus uparse/otus.fa -relabel OTU_
usearch -otutab  pools/seq.fq  -otus uparse/otus.fa -otutabout uparse/otutab.txt

9. Unoise

参考: 扩增子数据分析之去噪:UNOISE

mkdir unoise
usearch -unoise3 pools/derep.fa -zotus unoise/otus.fa  -minsize 4

10. 物种分类 sintax

参考: 使用SINTAX对16S序列进行分类:SINTAX 使用 rdp_16s_v16_sp 进行预测分类水平。

mkdir sintax
usearch -sintax uparse/otus.fa \
    -db  /biostack/databases/sintax/rdp_16s_v16_sp.udb  \
    -threads 24  -strand both  \
    -sintax_cutoff 0.8  -tabbedout sintax/otu.sintax

11. OTU注释

参考: 计算OTU表相对丰度:otutable_freq
对OTU表进行注释:otutable_annotation

mkdir otus
tabtk_row_reorder uparse/otutab.txt  <(seqtk_names uparse/otus.fa) > otus/otutab.txt
cut -f1,4  sintax/otu.sintax |  otutable_annotation  -  otus/otutab.txt  >otus/otutab_ann.txt
otutable_freq otus/otutab.txt  >otus/otutab_freq.txt
cut -f1,4  sintax/otu.sintax |  otutable_annotation  - otus/otutab_freq.txt >otus/otutab_freq_ann.txt
tabtk_xlsx  otus/otus.xlsx  otutab:otus/otutab.txt  \
                            otutab_freq:otus/otutab_freq.txt \
                            otutab_ann:otus/otutab_ann.txt  \
                            otutab_freq_ann:otus/otutab_freq.txt

biom convert  -i otus/otutab_ann.txt  -o  otus/otutab.biom  \
        --table-type "OTU table"  \
        --to-json  --process-obs-metadata taxonomy

12. 物种分类 krona 可视化

参考: OTU表生成 Krona 输入文件:otutable_krona

mkdir krona
otutable_krona  otus/otutab_ann.txt D1 >krona/D1.txt
ktImportText   -o krona/kron.html  krona/D1.txt,D1

现在应该可以看到:

image.png-493.8kB

13. 物种分类可视化

参考: 生成不同水平系统分类丰度表:otutable_abundance 和 otutable_level

mkdir taxonomy
otutable_level  -l g  otus/otutab_freq_ann.txt  >taxonomy/genus.txt

这是可以获得每个样本在不同分类水平的丰度情况,比如在属水平。 选择前30进行数据可视化:

taxon_table_rank taxonomy/genus.txt  30  >taxonomy/genus.30.txt
barplot.R  taxonomy/genus.30.txt   taxonomy/genus.30.pdf  genus.30
pdf2png taxonomy/genus.30.pdf

genus.30.png-92.7kB

也可以绘制热图。

14. 构建系统进化树

参考: 扩增子测序分析之构建OTU树: USEARCH和QIIME2

mkdir phylogeny
mafft  --auto  uparse/otus.fa   >phylogeny/otus.aln
trimal  -in phylogeny/otus.aln  -out  phylogeny/otus.trimal.fasta  -automated1
fasttree  -nt  phylogeny/otus.trimal.fasta >phylogeny/otus.nwk
treebender  --clear_internal_node_labels phylogeny/otus.nwk | treebender  --mid_point_root  >phylogeny/rooted.nwk

15. alpha多样性

参考: 扩增子数据分析之多样性指数: alpha多样性

mkdir alpha
usearch -alpha_div uparse/otutab.txt -output alpha/alpha.txt
alpha_table  alpha/alpha.txt  mapping_file.txt alpha/diversity
estimators_stats.R  alpha/diversity/richness.txt  alpha/diversity/richness_boxplot.pdf richness

可以绘制多样系指数的boxplot, 如下:

richness_boxplot.png-18kB

16. 稀释曲线

参考:扩增子数据分析之稀释曲线: alpha_div_rare 和 alpha-rarefaction

mkdir rare_curve
usearch -alpha_div_rare  otus/otutab.txt   --method fast --metric  shannon_e  -output  rare_curve/rare.txt
rarefaction_curve.R  rare_curve/rare.txt  mapping_file.txt   rare_curve/rare.pdf

可绘制如下图:

rare.png-73.5kB

17. beta多样性

参考:扩增子测序分析之多样性指数:beta 多样性

mkdir beta
usearch -beta_div uparse/otutab.txt -metrics bray_curtis,jaccard -filename_prefix /project/amplicon/beta/distmatrix/
usearch -beta_div uparse/otutab.txt -tree phylogeny/rooted.nwk -metrics unifrac,unifrac_binary -filename_prefix /project/amplicon/beta/distmatrix/

到此,我们使用USEARCH 完成了大部分的扩增子生物信息数据分析工作, 后续主要是多元统计以及可视化分析和功能预测, 比如,可视化组内的距离和组间的距离是否有差异(比如使用wilcox.test, 可以产生下图), 这部分内容先不放在这里。

weighted_unifrac.png-19.3kB

流程里有两点没有操作:

a. 最OTU表标准化
b. 去除线粒体和叶绿体来源的16S序列。

本文材料为 BASE (Biostack Applied bioinformatic SEies ) 课程 Linux Command Line Tools for Life Scientists 材料, 版权归 上海逻捷信息科技有限公司 所有。

Last Update: 2017-11-01 3:20 PM

扩增子测序分析之去噪 – DADA2及其平台化

引言

DADA2

DADA2(Divisive Amplicon Denoising Algorithm)的核心 去噪 然后使用 Amplicon Sequence Variants (ASVs)概念构建类OTU表;

具体算法可以参考下面几篇文章:

Denoising PCR-amplified metagenome data
DADA2: High-resolution sample inference from Illumina amplicon data
Bioconductor Workflow for Microbiome Data Analysis: from raw reads to community analyses

二、核心算法

DADA2 的核心算法:

1. 错误率模型计算,采用selfConsist无监督学习模型;
2. Rate λ_ji: 通过错误率模型,衡量扩增子序列i是否来在来自j;
3. 丰度P-value值,使用丰度信息衡量:该丰度的扩增子序列是否可以通过错误模型解释,使用了一个经验P值来判断,比如 OMEGA_A = 1e-40
4. 分裂算法( divisive partitioning)形成扩增子簇,簇心定义为denoise序列,簇大小定义denoise序列丰度。

三、DADA2实践

3.1 准备材料:

下载测试文件

wget http://www.mothur.org/w/images/d/d6/MiSeqSOPData.zip

下载训练集

Silva 128 [Silva taxonomic training data formatted for DADA2 (Silva version 128)][65]

wget https://zenodo.org/record/824551/files/silva_nr_v128_train_set.fa.gz
wget https://zenodo.org/record/824551/files/silva_species_assignment_v128.fa.gz

RDP trainset 16 [RDP taxonomic training data formatted for DADA2][76]

wget https://zenodo.org/record/801828/files/rdp_species_assignment_16.fa.gz
wget https://zenodo.org/record/801828/files/rdp_train_set_16.fa.gz

3.2 DADA2 安装

测试平台:R 3.3.0 wget https://zenodo.org/record/158955/files/rdp_species_assignment_14.fa.gz
wget https://zenodo.org/record/158955/files/rdp_train_set_14.fa.gz

UNITE
wget https://unite.ut.ee/sh_files/sh_general_release_10.10.2017.zip

DADA2 Pipeline

测试平台:R 3.3.0

DADA2 R Pakcage 安装

使用默认安装模式:

source("https://bioconductor.org/biocLite.R")
biocLite("dada2")

当前安装版本为: 1.2.2

使用devtools安装模式:

devtools::install_github("benjjneb/dada2")

安装最新版本: 1.5.8

3.3 测试流程

“`r
library(dada2);
packageVersion(“dada2”)

rdp_path <- "/project/task/dada2/RDP-16"
miseq_sop <- "/project/task/dada2/MiSeq_SOP"

fnFs <- sort(list.files(miseq_sop, pattern="_R1_001.fastq"))
fnRs <- sort(list.files(miseq_sop, pattern="_R2_001.fastq"))

sample.names <- sapply(strsplit(fnFs, "_"), `[`, 1)
fnFs <- file.path(miseq_sop, fnFs)
fnRs <- file.path(miseq_sop, fnRs)

plotQualityProfile(fnFs[1:2])

filt_path <- file.path(miseq_sop, "filtered")
filtFs <- file.path(filt_path, paste0(sample.names, "_F_filt.fastq.gz"))
filtRs <- file.path(filt_path, paste0(sample.names, "_R_filt.fastq.gz"))

filt <- filterAndTrim(fnFs, filtFs, fnRs, filtRs, truncLen=c(240,160),
      maxN=0, maxEE=c(2,2), truncQ=2, rm.phix=TRUE,
      compress=TRUE, multithread=TRUE)

derepFs <- derepFastq(filtFs, verbose=TRUE)
derepRs <- derepFastq(filtRs, verbose=TRUE)

names(derepFs) <- sample.names
names(derepRs) <- sample.names

errF <- learnErrors(filtFs, multithread=TRUE)
errR <- learnErrors(filtRs, multithread=TRUE)

dadaFs <- dada(derepFs, err=errF, multithread=TRUE)
dadaRs <- dada(derepRs, err=errR, multithread=TRUE)

mergers <- mergePairs(dadaFs, derepFs, dadaRs, derepRs, verbose=TRUE)
seqtab  <- makeSequenceTable(mergers)
seqtab.nochim <- removeBimeraDenovo(seqtab, method="consensus", multithread=TRUE, verbose=TRUE)

getN <- function(x) sum(getUniques(x))
track <- cbind(filt, sapply(dadaFs, getN), sapply(mergers, getN), rowSums(seqtab), rowSums(seqtab.nochim))
colnames(track) <- c("input", "filtered", "denoised", "merged", "tabled", "nonchim")
rownames(track) <- sample.names

taxa <- assignTaxonomy(seqtab.nochim,  file.path(rdp_path, "rdp_train_set_16.fa.gz"), multithread=TRUE)
seqtab_taxa <- cbind('#seq'=rownames(taxa), t(seqtab.nochim), taxa)
write.table(seqtab_taxa, "dada2-counts-taxon.txt", sep="\t", quote=F, row.names = F)

taxa.plus <- addSpecies(taxa, file.path(rdp_path, "rdp_species_assignment_16.fa.gz"))
seqtable.taxa.plus <- cbind('#seq'=rownames(taxa.plus), t(seqtab.nochim), taxa.plus)
write.table(seqtable.taxa.plus , "dada2-counts-taxon.species.txt", sep="\t", quote=F, row.names = F)

“`

3.4 流程解读

整个流程包含以下几个步骤:

序列过滤和切除
错误率计算
去噪音
合并双端序列
去除嵌合体
物种分类
物种水平鉴定

整个流程通过几个关键函数完成,filterAndTrim、derepFastq、learnErrors、dada、mergePairs、assignTaxonomy、addSpecies 等几个核心函数。

值得关注是addSpecies和assignSpecies, 逻辑是通过完全匹配的序列分类一致性(重点在这里, 两个问题: 1. 库应该足够大, 2. 序列库足够准确)去判断是否应该分配species信息;

DADA2 提供了silva_species_assignment_v128.fa.gzrdp_species_assignment_16.fa.gz用于识别可以分类到种水平信息, 该文件是通过对原始序列问题进行几个操作实现:

a. 去除分类不明确的序列,比如 Uncultured、 unclassified、Outgroup、Unidentified序列;
b. 去除分类不完整序列,要包含属和种水平信息;
c. 只保留属和种信息;

3.5 QIIME整合

QIIME 通过插件形式整合DADA2用于构建这中feature-table, 代替传统的OTU。

qiime dada2 命令可以执行特征表构建, 并使用MD5描述信息,支持单端(denoise-single)和双端模式(denoise-paired),

支持的命令行接口(参数基本可以和DADA2的R函数对应):

--i-demultiplexed-seqs: demultiplexe后的序列;
--p-trunc-len-f: 针对forward序列,截断序列长度大小,设置为0不执行截断操作;
--p-trunc-len-r: 针对reverse序列,截断序列长度大小,设置为0不执行截断操作;
--o-table: FeatureTable,用于下游各种分析。
--o-representative-sequences: 代表序列,簇心。
--p-hashed-feature-ids:  使用Hash值做序列标识符。
--p-no-hashed-feature-ids:不使用Hash值做序列标识符。
--p-n-reads-learn: 用于无监督学习的序列数。
--p-n-threads: 线程数;

实例:

qiime dada2 denoise-paired \
    --i-demultiplexed-seqs paired-end-sequences.qza \
    --p-trunc-len-f  0   \
    --p-trunc-len-r  0   \
    --o-table  table.qza \
    --o-representative-sequences \
    rep-seqs.qza

本文材料为 BASE (Biostack Applied bioinformatic SEies ) 课程 Linux Command Line Tools for Life Scientists 材料, 版权归 上海逻捷信息科技有限公司 所有.

Last Upate: 2017/10/16 21:20