HTSeq

安装C/C++编译环境
Ubuntu缺省情况下,并没有提供C/C++的编译环境,运行以下命令安装:
sudo apt-get install build-essential

安装setuptools
setuptools可以更方便的创建和发布Python包,HTSeq的安装依赖setuptools
下载setuptools-20.3.1.tar.gz
tar zxvf setuptools-20.3.1.tar.gz
cd setuptools-20.3.1/
sudo python setup.py install

安装Cython
Cython是Python的C扩增,可以让Python程序更加高效。HTSeq依赖cython。
sudo apt-get install cython

安装Numpy
Numpy是Python科学计算基本库,提供了许多高级的数值编程工具。HTSeq依赖Numpy。
sudo apt-get install python-numpy

安装Zlib
Zlib是zlib是提供数据压缩用的函数库。pysam依赖Zlib。
下载zlib128.zip
unzip zlib128.zip
cd zlib-1.2.8/
./configure
make
sudo make install

安装pysam
pysam可以操作很多基因组数据,包括SAM/BAM/VCF/BCF/BED/GFF/GTF/FASTA/FASTQ。HTSeq依赖pysam。
下载pysam-0.9.0.tar.gz
tar zxvf pysam-0.9.0.tar.gz
cd pysam-0.9.0/
sudo python setup.py install

安装HTSeq
下载HTSeq-0.6.1.tar.gz
tar zxvf HTSeq-0.6.1.tar.gz
cd HTSeq-0.6.1/
sudo python setup.py install

报错

RNA-Seq数据分析流程

RNA-Seq分析流程用的最多的是Tuexo workflow,如下图所示,
tuxedo_workflow

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

这套流程最大的优势在于,软件衔接非常好,Tophat比对、cufflinks组装、cuffdiff差异表达,Nature protocol还发表了一篇详细的操作手册。

但是,很多人都说cufflinks推出的FPKM算表达量并不好。cufflinks这帮大牛后面也意识到了问题,推出了cuffquant、cuffnorm、cuffdiff2,做出了一些修正,效果还是差强人意。cufflinks这帮大牛又开发出HISAT来替代Tophat,StringTie替代cufflinks。

HISAT2比对,StringTie组装,使用HTseq提供的脚本重新计算转录本对应的reads数目,用DEseq2进行差异表达分析,是最靠谱的选择。

转录组差异表达分析工具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。

有空继续写!!!

Python脚本过滤大于200nt的转录本

lncRNA是大于200nt的,有时候我们就需要过滤出大于200bp的转录本,可以用Python写一个脚本轻松实现这个功能。

import sys, getopt
 
opts, args = getopt.getopt(sys.argv[1:], "hi:o:")
input_file = ""
output_file = ""
for op, value in opts:
    if op == "-i":
        input_file = value
    elif op == "-o":
        output_file = value
    elif op == "-h":
        print("Usage: python filter.py -i input.fasta -o output.fasta")
        sys.exit()
 
Input_file = open(input_file)
Output_file = open(output_file,"w")
fa_Name = ""
fa_Seq = ""
 
for line in Input_file.readlines():
    line = line.rstrip()
    if line[0] == ">":
        if len(fa_Seq) >= 200:
            Output_file.write(fa_Name + "\n")
            while len(fa_Seq) > 60:
                Output_file.write(fa_Seq[:60] + "\n")
                fa_Seq = fa_Seq[60:]
            Output_file.write(fa_Seq + "\n")
        fa_Seq = ""
        fa_Name = line
    else:
        fa_Seq = fa_Seq + line

使用CNCI分析lncRNA

CNCI是中科院计算所赵屹团队开发的一款从转录组中分析编码RNA和非编码RNA的软件。赵屹团队在非编码RNA领域做了很多出色的工作,建立了目前最权威的非编码RNA数据库NONCODE。

文章发表在Nucleic Acids Research:http://nar.oxfordjournals.org/content/early/2013/08/06/nar.gkt646.long
Github主页:https://github.com/www-bioinfo-org/CNCI

CNCI的安装:
CNCI使用了SVM(支持向量机)分类,其安装过程主要是编译libsvm。
tar zxvf CNCI_version2.tar.gz
cd CNCI_package
unzip libsvm-3.0.zip
cd libsvm-3.0
make
cd ..

CNCI的运行:
基本命令为:python CNCI_package/CNCI.py -f novel.fasta -o CNCI_out -m ve -p 4
参数说明:
-f 输入fasta文件(可以使用-g参数输入GTF文件,但是同时需要使用-d参数指定参考基因组的目录)
-o 输出结果目录
-m 指定模式,脊椎动物选择ve,植物选择pl
-p 指定CPU核数
更多用法参看帮助文档

小提示:
CNCI的运行目录一定要在CNCI_package所在目录,不要到CNCI_package目录下运行CNCI,否则会报错。
fasta文件中的注释,也就是”>”开头的那一行,不要有空格,可能也会报错。如果有空格,考虑将空格替换为下划线。
软件有点小瑕疵吧……

使用CPC分析lncRNA

CPC是最早的从转录组中分析lncRNA的软件,很多lncRNA研究中用的都是它。它基于blast,把未知转录本和已知蛋白库对比,从而筛选出编码和非编码的转录本。CPC可靠性很高,缺点就是太慢太慢了,需要几天时间。CPC可以在网页版上提交任务,运行好了,将结果下载下来。

文章发表在Nucleic Acids Research:http://nar.oxfordjournals.org/content/35/suppl_2/W345.long
主页:http://cpc.cbi.pku.edu.cn/

CPC安装前需要:
CPC需要使用blast,它调用的是blastall,也就是老版本的blast,而不是新版本的blast+。
老版本的blast下载地址:http://yunpan.cn/cwGh5BSDbdYIe 访问密码 9df4
需要使用蛋白质库,UniRef90或者NCBI的nr都可以,用formatdb命令建库时,必须命名为”prot_db”, 且放在CPC安装目录下的data目录下面。

CPC的安装:
下载cpc-0.9-r2.tar.gz
tar -zxvf cpc-0.9-r2.tar.gz
cd cpc-0.9-r2/libs/libsvm
tar -zxvf libsvm-2.81.tar.gz
cd libsvm-2.81
make clean && make
cd ../..
tar -zxvf estate.tar.gz
cd estate
make clean && make

建立本地blast数据库:
cd cpc-0.9-r2/data
formatdb -i (your_fasta_file) -p T -n prot_db

运行CPC:
cd cpc-0.9-r2/
bin/run_predict.sh (input_seq) (result_in_table) (working_dir) (result_evidence)
run_predict.sh好像是远程blast,建议运行run_predict_local.sh,并把这个文件中blast_opts=”$blast_opts -a 2″; # 2CPUs, boost the performance这句话中的2,改成你实际电脑使用的CPU核数。

使用CPAT分析lncRNA

CPAT(Coding Potential Assessment Tool)可以从RNA-Seq分析得到的转录本中筛选出编码的和非编码的RNA。
文章发表在Nucleic Acids Research:http://nar.oxfordjournals.org/content/41/6/e74.long
代码托管在SourceForge:http://rna-cpat.sourceforge.net/

安装:
安装CPAT前,需要安装gcc、python2.7、numpy、cython、R。我的系统是Ubuntu 14.04 64,gcc和python2.7已经安装好了,其他三个使用如下命令安装:
sudo apt-get install r-base-core
sudo apt-get install cython
sudo apt-get install python-numpy
下载最新版本CPAT1.22,按照说明安装
tar -zxvf CPAT-1.2.2.tar.gz
cd CPAT-1.2.2
sudo python setup.py install

运行CPAT:
cpat.py -r hg19.fa -g Human_test_coding_mRNA_hg19.bed -d ../dat/Human_logitModel.RData -x ../dat/Human_Hexamer.tsv -o output
参数说明:
-r 指定参考基因组
-g 输入的转录本序列。如果是BED格式,必须-r指定参考基因组;如果是FASTA格式,不需要指定参考基因组,即使使用-r参数也会被忽略。
-d 预制好的模型(Prebuilt training model)(CPAT自带人、鼠、果蝇、斑马鱼的模型)
-x 预制好的六聚体频率表(Prebuilt hexamer frequency table)(CPAT自带人、鼠、果蝇、斑马鱼的六聚体频率表)
-o 输出

CPAT优点:速度极快!