Archives

Categories

Linux 平台轻松解决Python3环境: conda

一、 引子(需求场景)

因为 Python2 和 Python3 之间的语法差异, 现在Python开发者主要有两大阵营,Python 2.7.x 版本阵营 和 Python 3 阵营, 为了使用不同阵营的包, 导致我们要维护两套Python系统。 当前Linux发行版本主要还是Python2环境,所以需要额外一套方案使得Python2和Python3共存, 最简单的方案应该就是切换环境变量,这样做的好处是我们可以包括任意多的版本,这里使用conda来切换环境。

这篇可以作为介绍Bioconda的引子,只介绍如何安装 miniconda 解决 Python 环境的问题,如何使用bioconda 安装生物信息软件,后面再讲。

二、使用Conda进行Python包管理

下面就讲几点:

miniconda 介绍 miniconda 安装 设置环境变量 安装python包 设定和管理特定Python环境 1 miniconda 介绍:

先摆出几个概念,Anaconda, conda, Miniconda, pip,可以参考官方文档: What’s the difference between Anaconda, conda, and Miniconda?, 也可以参考一篇博客:Anaconda、Miniconda、Conda、pip的相互关系 , 讲的很清楚, 为了完整性,这里简单介绍:

conda, 包管理器,安装和管理包的命令接口; […]

对扩增子序列进行引物匹配: usearch 和 primersearch

经典问题场景: 数据:16S 扩增子序列,已完成双端合并,amplicon.fasta 文件, 一对引物序列:primers.txt

我们想做:

需要知道引物是否匹配,过滤掉一些噪音序列; 多对引物在同一个文库,需要根据引物把目标序列提取出来;

下面介绍两个比较简单的引物识别程序(双端引物序列模式)

一、USEARCH

强大的数据处理工具集合,但是: 只能免费获取32位版本,最大内存使用4G, 大部分的数据分析任务还是可以应付。

下载地址:32Bit , 选择需要下载的版本,地址会通过邮件进行发送。

如果想购买64bit 版本可以移步:buy64bit , 也可通过代理商购买,现在的销售模式是大版本免费升级,当前大版本为10.0

这里只涉及到引物匹配子命令:search_pcr

命令行接口:

-db, 引物序列文件, fasta格式; -strand, 支持 plus 和 both 两种模式 -maxdiffs, 单端引物去最大差异 -minamp,匹配区域的最小长度; -maxamp,匹配区域的最大长度; -pcrout, 制表符分隔的文件,记录比对的详细信息; -threads, 线程数; -ampout, 匹配区域的序列文件; -log, 输出log 信息;

-pcr_strip_primers 可以切除引物区域,如果引物中存在简并碱基,可能程序会报错,使用需要谨慎。

cat primers.txt >fw GTGCCAGCMGCCGCGG >rw CCGTCAATTCMTTTRAGTTT

现在执行:

usearch10.0.240_i86linux32 […]

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

扩增子测序常用的模式是在双端引物两边各加上几个碱基的条形码(barcode)进行区分样本,这样三十到四十个样本可以混在一起进行建库,节省实验成本。

一、seqtk_demultiplex 介绍

程序见 Github seqtk_utils 仓库。

命令行接口

seqtk_demultiplex 因为参数比较多,所以全部采用了选项模式传递命令行参数,接口如下:

$ seqtk_demultiplex Usage: barcode_demutiplex [options] Options: -1 STR forward reads file; -2 STR forward reads file; -b STR barcode mapping file; -d STR output filename directory; -l INT minimum barcode length, default: 5; -v print version number. -h display usage. Example: barcode_demutiplex -b barcodes.txt -1 […]

对TSV文件进行行筛选:tabtk_subset

继续针对制表符分隔的文件操作进行解读,今天介绍一个工具: abtk_subset ,与先前介绍的一个小工具:seqtk_subseq 功能模式差不多,

如何根据列表快速提取fasta/q序列子集/补集:seqtk_subseq

实现方式也差不多, 使用一个HASH SET 存储列表标识符, 然后遍历库文件。

使用 grep 也可以解决TSV文件筛选问题, 但是:

1. grep 无法指定列,求集合操作; 2. 查询列表比较大的时候,效率很低 一、tabtk_subset 介绍

tabtk_subset 命令行接口比较简单, 两个必须参数,一个待搜索库,一个标识符列表;

$tabtk_subset Usage: tabtk_subset [options] <tsv> <list> Options: -c INT target col for select, default: 1 -s INT subset option: 1 intersection, 0 supplementary set, default: 1 -v print vesion number

该程序接收两个可选参数:

[…]

优雅的TSV文件进行列注释: tabtk_definition

今天开始介绍一些处理制表符分隔文件的小工具,在生物信息数据分析过程中大部分输出格式都是以’\t’分隔的文件,比如 BLAST 输出格式,我们习惯使用 ‘-fmt 6’ 或者 ‘-fmt 7’.

此类工具比较多(包含 csv 类),比如: tabtk , tsv-utils-dlang, csvkit, 还是老规矩, 尽量避免重复造轮子。

一、tabtk_definition 介绍

tabtk_definition 就是使用一个key/value的键值对 db 对 指定文件的列进行注释,命令行接口:

$tabtk_definition Usage: tabtk_definition [options] <db> <tab> Options: -c INT target col for annotation, default: 1 -d char delimitor between key and definition, default: ‘ ‘ -h char* header name, default: ‘definition’ […]

使用RDP对16S序列进行分类:rdp_classifier

今天探讨的是使用RDP分类器对16S/ITS/18S序列进行分类;

一、RDP分类器

RDP分类器是老牌的序列分类软件,采用了 naïve Bayesian classifier,最低分类等级到属水平。 QIIME(1.9 版本)使用的是RDP比较老的版本(2.2),当前版本为 2.12, 这个符合QIIME 一贯的作风,OTU 聚类 usearch 还一直采用Usearch比较老的算法版本(Uclust).

QIIME老用户请切换至 QIIME2, 全新的设计,还带有很多新的特性。

具体算法原理请参考:

Wang, Q, G. M. Garrity, J. M. Tiedje, and J. R. Cole. 2007. Naïve Bayesian Classifier for Rapid Assignment of rRNA Sequences into the New Bacterial Taxonomy. Appl Environ Microbiol. 73(16):5261-7

二、下载安装

分类器:

可以从Github上下载源代码编译(仓库:RDPTools), 也可以直接从 SourceForge 下载编译程序:

[…]

根据引物序列提取barcode: seqtk_barcode_inspect

seqtk_barcode_inspect 这个工具的使用场景比较少见,也许在这样的场景下会使用到:

1. 扩增子测序,barcode和引物在一起: barcode+primer 2. 你从NCBI下载了一些扩增子序列原始数据,比如16S测序结果,barcode没有被去除 3. 你的数据分析流程需要知道barcode是什么,然后将所有序列合并在一起,重新拆分和去除barcode

如果你有上面的使用场景,也许这个小程序可以帮助你,设计思路很简单(如果不嫌麻烦,可以肉眼去看):

1. 使用引物前几个非简并碱基进行定位,目标的前几个碱基就是我们要找的barcode序列 2. 默认在前10条序列搜索,如果搜索到非兼并碱基,终止,返回barcode序列,如果没有匹配,终止程序;

程序接口:

Usage: seqtk_barcode_inspect [options] <in.fa> <primer> <sample> Options: -t INT top n sequence to search, default: 10 -v print version number

使用用例:

$seqtk_barcode_inspect -t 10 raw.split.PCF.1.fq GTGCC F F AGCACAT

“GTGCC” 为 16S 引物对V3的前几个碱基, AGCACAT就是对应的barcode

本文材料为 BASE (Biostack Applied bioinformatic SEies […]

ENA和SRA数据预处理: SRA Toolkit 和 seqtk_sra

前面一篇文章讨论了如何使用 aspera 上传/下载 神器 Aspera:

数据高速上传下载的利器:Aspera及其在生命科学中的应用

告诉我们如何高效(其实最重要的还是本地带宽, wget照样可以高效!)下载数据,遗留点下载数据后如何处理的问题,下面可以继续探讨。

SRA和ENA下载的数据,一般非标准的fastq文件格式,有时需要预处理才能进行下一步数据分析,比如:

1. 从sra文件提取fastq/sff 文件 2. 格式化fastq文件 一、数据下载

按照先前套路下载一个测试文件:http://www.ebi.ac.uk/ena/data/view/ERS980542

ENA

ascp_ena_download ftp://ftp.sra.ebi.ac.uk/vol1/ERA536/ERA536205/fastq/RRCI25.31.fastq.gz ascp_ena_download ftp://ftp.sra.ebi.ac.uk/vol1/ERA536/ERA536205/fastq/RRCI25.31.fastq.gz

SRA

ascp_sra ERR1145336 二、数据格式化

2.1 安装 NCBI SRA Toolkit

下载地址: https://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/2.8.2/

ascp_ncbi_download ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/2.8.2/sratoolkit.2.8.2-centos_linux64.tar.gz cd sratoolkit.2.8.2-centos_linux64/bin export PATH=$PWD:$PATH

2.2 提取fastq文件 ENA下载的就是fastq文件,所以不需要额外处理,现在针对 sra 文件, 使用 fastq-dump 程序

fastq-dump的命令行参数:

$fastq-dump Usage: fastq-dump [options] <path> [<path>…] fastq-dump [options] […]

数据高速上传下载的利器:Aspera及其在生命科学中的应用

一、Aspera 介绍

Aspera是IBM公司的一款高速传输软件,在 2013年被IBM收购 前属于Aspera, Inc 公司。

IBM® Aspera® 快速文件传输使企业能够利用现有的 WAN 基础架构,以可预测、可靠且安全的传送方式,快速传输大型文件和数据集(不论是结构化还是非结构化数据),而与文件大小、传输距离和网络状况无关。该文件传输软件包含客户端和服务器软件包,可以使用状态不佳的网络,跨多个位置(甚至是远程位置),传输对时间要求严格的文件。

本材料来源:http://www-03.ibm.com/software/products/zh/high-speed-file-transfer

一句话:就是传输速度快,基本可以占满客户的带宽。

二、Aspera 在生命科学领域内的应用

NCBI(数据下载接口)/ENA/GigaScience 等都采用Aspera技术提供数据下载服务。

三、 下载安装

这部分内容只包含Linux(测试平台: CentOS 7 操作系统)平台的安装以及常用数据的操作,NCBI提供了Web接口,通过浏览器插件也可以很方便的下载数据。

下载和安装:

wget http://download.asperasoft.com/download/sw/connect/3.7.2/aspera-connect-3.7.2.141527-linux-64.sh sh aspera-connect-3.7.2.141527-linux-64.sh

如需其它平台请移步: http://downloads.asperasoft.com/en/downloads/8?list

文件会被安装在:/home/biostack/.aspera 目录 目录结构如下:

. └── connect ├── bin ├── etc ├── iso-swid ├── lib ├── localization ├── notices.txt ├── plugins ├── product-info.mf ├── res […]

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

本部分内容主要介绍如何快速获取fasta/q文件的碱基数和长度信息,使用 klib 实现,klib介绍参考前文章:

》如何根据列表快速提取fasta/q序列子集/补集:seqtk_subseq

一、工具介绍

程序接口:

秉承支持streaming的原则,接口很简单:

$seqtk_info Usage: seqtk_info <fasta/q> version 0.0.1 $seqtk_fqchk Usage: seqtk_fqchk [options] <in.fq> Options: -q INT offset value: 33 for Sanger quality, 64 for Illumina quality, default: 33 -v print version number

接受(gz压缩)文件名或者标准输入流(使用 – 占位)

示例介绍:

time pigz -dk -c -p 2 test.fq.gz | seqtk_info – #sequence base min_len max_len […]