Python科学计算工具箱Anaconda

Python语言简洁,学习起来非常快,使用Python的科研工作者越来越多。Python拥有丰富的扩展包,这使得它可以轻松处理各种问题。
Python用于科学计算的包括但不限于如下:
Numpy:提供了许多高级的数值编程工具,如:矩阵数据类型、矢量处理,以及精密的运算库等。
Scipy:提供许多科学计算函数库,包括统计,优化,整合,线性代数模块,傅里叶变换,信号和图像处理,常微分方程求解器等。
Pandas:基于Numpy,提供了大量库和一些标准的数据模型,比如二维数组。
Matplotlib:绘图包,旨在替代MATLAB。
……
Anaconda是一个集合,包括conda,python解释器,一些第三方库。
conda用于包管理和环境管理。包管理与pip类似,管理python第三方库。conda将一切都看成包,包括conda本身和python解释器。环境管理类似于virtualenv,能够允许用户使用不同版本的python和不同的第三方库环境,并灵活切换。
用Anaconda来管理环境和第三方库将是非常方便的,而且自己安装这些库对新手来说是一个不小的挑战。
Anaconda包括了科学、数学、工程、数据分析中最受欢迎的300多个python包,支持python2/3, Windows/Linux/Mac,而且还是非常有好的安装包(.exe/.sh/.pkg)。一键安装科学计算中常用的包,爽!
与在Windows中安装Anaconda类似,在Linux上安装Anaconda也非常简单。

sh Anaconda3-5.1.0-Linux-x86_64.sh

安装过程中会提示安装的路径PREFIX=/home/chenwen/anaconda3,还会提示是否将此路径添加到用户的环境变量(/home/chenwen/.bashrc)。
安装Anaconda之后,Win键+R,输入’cmd’,进入命令行窗口,输入如下命令验证Anaconda是否安装成功,并查看它的版本。

conda --version

Anaconda安装成功之后,我们首先将软件源修改为国内镜像,运行如下命令。

conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/
conda config --set show_channel_urls yes

我们现在就可以通过create命令来创建虚拟环境了。

conda create --name tensorflow python=3.5

这样我就创建了一个名叫tensorflow的虚拟环境,当然你可以用任何你喜欢的名字代替。我还指定了Python版本是3.5。
环境创建完毕后,我们可以使用info命令查看所有环境。

conda info --envs

如果需要切换到该虚拟环境,可以使用activate命令。

activate tensorflow

如果想取消激活,使用deactivate命令即可。

deactivate tensorflow

如果想删除虚拟环境,使用remove命令。

conda remove --name tensorflow --all

当我们切换到某个虚拟环境后,就可以进行第三方库的管理了。如果只是简单的学习Python,直接用系统环境就可以了,没有必要创建虚拟环境。
查看安装的包,可以使用list命令。

conda list

如果你发现某个包默认没有装,比如scikit-learn,可以使用conda install 命令安装。

conda install scikit-learn

更新包,可以使用conda update命令

conda update scikit-learn

此外,我还非常喜欢Anaconda带的IDE工具spyder,比起python官方自带的IDLE强很多,但是又不复杂。
(2018年5月22日更新!)

Python分割fasta文件

CPC可以从转录组得到的fasta序列中分析出哪些是编码RNA,哪些是lncRNA,但是它的速度实在是太慢了。我写了一个脚本,将fasta文件分割成更小的文件,在不同计算节点上面同时运行CPC以提高速度。CPC在balst比对这一步是支持多线程的,但是并不支持跨节点。对于某些不直接支持多线程的程序,将输入文件分割,再同时运行程序,是一种常用思路。
用法:python3 split_fasta.py -i input.fasta -o prefix -x split_number
-i 指定输入的fasta文件
-o 指定输出文件的前缀,默认是split_
-x 指定分割文件中fasta序列的数目,最后一个文件小于等于指定的数目,默认是1
命令:grep “>” file.fasta |wc -l
可以用来查看fasta文件中,有多少条fasta序列
代码如下:
split_fasta.py

# -*- coding: utf-8 -*-
import sys, getopt
# get parameter
opts, args = getopt.getopt(sys.argv[1:], "hi:o:x:")
X = 1
input_file = ""
prefix = "split_"
for op, value in opts:
    if op == "-i":
        input_file = value
    elif op == "-o":
         prefix = value
    elif op == "-x":
        X = int(value)
    elif op == "-h":
        print("Usage: python3 split_fasta.py -i input.fasta -o prefix -x split_number")
        print("default prefix = split_")
        print("default split_number = 1")
        sys.exit()
FA_in_file = open(input_file, "r")
# read fasta file to a list
fa_Info = []
fa_Seq = []
fa_Num = -1
for Y in FA_in_file.readlines():
    Y = Y.rstrip()
    if Y[0] == ">":
        fa_Info.append(Y)
        fa_Num = fa_Num + 1
        fa_Seq.append("")
    else:
        fa_Seq[fa_Num] = fa_Seq[fa_Num] + Y
# split the fasta list to multipe files
file_Num = (fa_Num + 1)//X + 1
for i in range(file_Num):
    exec(prefix + str(i + 1) + ' = open("' + prefix + str(i + 1) + '.fasta"' + ', "w")')
    start = i * X
    end = (i + 1) * X
    if end > fa_Num + 1:
        end = fa_Num + 1
    for j in range(start, end, 1):
        exec(prefix + str(i + 1) + '.write(fa_Info[j] + "\n")')
        while len(fa_Seq[j]) > 60:
            exec(prefix + str(i + 1) + '.write(fa_Seq[j][:60] + "\n")')
            fa_Seq[j] = fa_Seq[j][60:]
        else:
            exec(prefix + str(i + 1) + '.write(fa_Seq[j] + "\n")')
    exec(prefix + str(i + 1) + '.close()')
FA_in_file.close()

转录组差异表达分析工具Ballgown

Ballgown是分析转录组差异表达的R包。
软件安装:
运行R,
source(“http://bioconductor.org/biocLite.R”)
biocLite(“ballgown”)
R会自动安装Ballgown,及相应的依赖包。
Ballgown的输入文件
StringTie使用-B参数直接生成Ballgown的输入文件,Cufflinks的输出结果需要使用Tablemaker生成Ballgown的输入文件。
一个有5个输入文件,分别是:
e_data.ctab,外显子水平表达量;
i_data.ctab,内含子水平表达量;
t_data.ctab,转录本水平表达量;
e2t.ctab,外显子与转录本的对应关系;
i2t.ctab,内含子与转录本的对应关系。
Tablemaker
tablemaker -p 4 -q -W -G merged.gtf -o sample01_output sample_01/accepted_hits.bam
-p 指定线程数
-q 去冗余
-W 运行模式是tablemaker,而不是Cufflinks模式
-G 指定组装好的GTF文件,这个文件由cuffmerge生成
BAM文件是最后一个参数,这个文件由Tophat生产
运行Ballgown
运行R
载入Ballgown包
library(ballgown)
载入数据,并创建一个ballgown项目
bg = ballgown(dataDir=’D:\extdata’, samplePattern=’sample’, meas=’all’)
指定分组,及重复样本数目:
pData(bg) = data.frame(id=sampleNames(bg), group=rep(c(1,0), each=3))
差异表达分析:
stat_results = stattest(bg, feature=’transcript’, meas=’FPKM’, covariate=’group’)
有些地方我写得不是很清楚,有一个英文的Tutorial写得非常好。
https://github.com/alyssafrazee/ballgown

转录组组装工具StringTie

StringTie由约翰霍普金斯大学联合德州大学西南医学中心开发,能够组装转录本并预计表达水平。它应用网络流算法和可选的denovo组装,将复杂的数据集组装成转录本。 StringTie-0
 
 
 
 
 
 
 
 
 
相对于其他拼接软件(Cufflinks, IsoLasso, Scripture,Traph等),StringTie能够拼接出更完整、更准确的基因,并且StringTie采用拼接和定量同步进行,相对于其他方法,其定量结果更加准确。 StringTie-1
文中指出,对于从人类血液中获得的reads,StringTie正确组装了10,990个转录本,而Cufflinks只组装了7,187个。而对于模拟的数据集,StringTie正确组装了7,559个转录本,比Cufflinks的6,310个提高了20%。此外,它的运行速度也比其他组装软件更快。
软件安装
下载StringTie并解压:
tar zxvf stringtie-1.1.2.Linux_x86_64.tar.gz
将HISAT2目录添加到环境变量:
vi ~/.bashrc
在文件末位添加:
export PATH=/home/chenwen/bin/stringtie-1.1.2.Linux_x86_64:$PATH
保存退出
source ~/.bashrc
StringTie的输入文件
HISAT的输出文件为SAM格式,需要经过两步转换成StringTie可以使用的BAM格式:
1. SAM转BAM,并排序:
samtools view -S -b input.sam | samtools sort – input.sorted
2.修改HI标签:
samtools view -h input.bam | perl -ne ‘if(/HI:i:(d+)/) {$m=$1-1; $_ =~ s/HI:i:(d+)/HI:i:$m/} print $_;’| samtools view -bS – > input.correct.bam
运行StringTie
stringtie SRR534294.correct.bam -p 16 -G genes.gtf -B -o stout/transcripts.gtf
-p 指定线程数,默认1
-G 指定参考的转录组注释文件
-B 生成用于Ballgown 分析的文件
-o 指定输出文件
更多参数请查看HISAT2的操作手册:
https://ccb.jhu.edu/software/stringtie/index.shtml?t=manual

 

RNA-Seq基因组比对工具HISAT2

HISAT2是TopHat2/Bowti2的继任者,使用改进的BWT算法,实现了更快的速度和更少的资源占用,作者推荐TopHat2/Bowti2和HISAT的用户转换到HISAT2。
官网:
https://ccb.jhu.edu/software/hisat2/index.shtml
HISAT2安装
下载HISAT2-2.0.1,并解压:

unzip hisat2-2.0.1-beta-Linux_x86_64.zip

将HISAT2目录添加到环境变量:

vi ~/.bashrc

在文件末位添加:

export PATH=/lustre/home/lcn/chenwen/bin/hisat2-2.0.1-beta:$PATH

保存退出

source ~/.bashrc

建立索引
建立基因组索引

hisat2-build –p 4 genome.fa genome

建立基因组+转录组+SNP索引:
bowtie2的索引只有基因组序列信息,tophat2比对时,转录组信息通过-G参数指定。HISAT2建立索引时,就应该把转录组信息加进去。
HISAT2提供两个Python脚本将GTF文件转换成hisat2-build能使用的文件:

extract_exons.py Homo_sapiens.GRCh38.83.chr.gtf > genome.exon
extract_splice_sites.py Homo_sapiens.GRCh38.83.chr.gtf > genome.ss

此外,HISAT2还支持将SNP信息加入到索引中,这样比对的时候就可以考虑SNP的情况。这仍然需要将SNP文件转换成hisat2-build能使用的文件:

extract_snps.py snp142Common.txt > genome.snp

最后,将基因组、转录组、SNP建立索引:

hisat2-build -p 4 genome.fa --snp genome.snp --ss genome.ss --exon genome.exon genome_snp_tran

官网提供了人和小鼠的索引文件下载,压缩包有make_grch38_tran.sh文件,详细记录了创建索引的过程。
运行HISAT2

hisat2 -p 16 -x ./grch38_tran/genome_tran -1 SRR534293_1.fastq -2 SRR534293_2.fastq –S SRR534293.sam

-x 指定基因组索引
-1 指定第一个fastq文件
-2 指定第二个fastq文件
-S 指定输出的SAM文件
更多参数请查看HISAT2的操作手册:
https://ccb.jhu.edu/software/hisat2/manual.shtml
官方操作手册简要版
用法:
hisat2 [options]* -x <hisat2-idx> {-1 <m1> -2 <m2> | -U <r> | –sra-acc <SRA accession number>} [-S <hit>]
主要参数:
-x <hisat2-idx>
参考基因组索引文件的前缀。
-1 <m1>
双端测序结果的第一个文件。若有多组数据,使用逗号将文件分隔。Reads的长度可以不一致。
-2 <m2>
双端测序结果的第二个文件。若有多组数据,使用逗号将文件分隔,并且文件顺序要和-1参数对应。Reads的长度可以不一致。
-U <r>
单端数据文件。若有多组数据,使用逗号将文件分隔。可以和-1、-2参数同时使用。Reads的长度可以不一致。
–sra-acc <SRA accession number>
输入SRA登录号,比如SRR353653,SRR353654。多组数据之间使用逗号分隔。HISAT将自动下载并识别数据类型,进行比对。
-S <hit>
指定输出的SAM文件。
输入选项:
-q
输入文件为FASTQ格式。FASTQ格式为默认参数。
-qseq
输入文件为QSEQ格式。
-f
输入文件为FASTA格式。
-r
输入文件中,每一行代表一条序列,没有序列名和测序质量等。选择此项时,–ignore-quals参数也会被选择。
-c
此参数后是直接比对的序列,而不是包含序列的文件名。序列间用逗号隔开。选择此项时,–ignore-quals参数也会被选择。
-s/–skip <int>
跳过输入文件中前条序列进行比对。
-u/–qupto <int>
只使用输入文件中前条序列进行比对,默认是没有限制。
-5/–trim5 <int>
比对前去除每条序列5’端个碱基
-3/–trim3 <int>
比对前去除每条序列3’端个碱基
–phred33
输入的FASTQ文件碱基质量值编码标准为phred33,phred33为默认参数。
–phred64
输入的FASTQ文件碱基质量值编码标准为phred64。
–solexa-quals
将Solexa的碱基质量值编码标准转换为phred。
–int-quals
输入文件中的碱基质量值为用空格分隔的数值,而不是ASCII码,例如40 30 30 40。
有空继续写!!!