使用HTSeq进行有参转录组的表达量计算

1. HTSeq简介

HTSeq是使用Python编写的一支用于进行基因Count表达量分析的软件,能根据SAM/BAM比对结果文件和基因结构注释GTF文件得到基因水平的Counts表达量。HTSeq进行Counts计算的原理非常简单易懂,容易上手。

2. HTSeq安装

PYPI下载HTSeq的Python包
$ wget https://pypi.python.org/packages/46/f7/6105848893b1d280692eac4f4f3c08ed7f424cec636aeda66b50bbcf217e/HTSeq-0.7.2.tar.gz
$ tar zxf HTSeq-0.7.2.tar.gz
$ cd HTSeq-0.7.2
$ /opt/sysoft/Python-2.7.11/bin/python setup.py build
$ /opt/sysoft/Python-2.7.11/bin/python setup.py install
$ cd ../ && rm HTSeq-0.7.2 -rf

3. HTSeq使用

3.1 HTSeq的Count模式

HTSeq计算counts数有3种模式,如下图所示(ambiguous表示该read比对到多个gene上;no_feature表示read没有比对到gene上):
HTSeq Count模式

3.2 HTSeq的使用命令

HTseq安装完毕后,在Python软件的bin目录下生成htseq-count命令。
htseq-count运行简单示例:

对于非链特异性真核转录组测序数据
$ /opt/sysoft/Python-2.7.11/bin/htseq-count -f sam -r name -s no -a 10 -t exon -i gene_id -m union hisat2.sam genome.gtf > counts_out.txt
对于链特异性测序真核转录组测序数据
$ /opt/sysoft/Python-2.7.11/bin/htseq-count -f sam -r name -s reverse -a 10 -t exon -i gene_id -m union hisat2.sam genome.gtf > counts_out.txt
对于非链特异性原核生物转录组测序数据
$ /opt/sysoft/Python-2.7.11/bin/htseq-count -f sam -r name -s no -a 10 -t exon -i gene_id -m intersection-strict bowtie2.sam genome.gtf > counts_out.txt

htseq-count命令的常用参数:

-f | --format     default: sam
  设置输入文件的格式,该参数的值可以是sam或bam。
-r | --order     default: name
  设置sam或bam文件的排序方式,该参数的值可以是name或pos。前者表示按read名进行排序,后者表示按比对的参考基因组位置进行排序。若测序数据是双末端测序,当输入sam/bam文件是按pos方式排序的时候,两端reads的比对结果在sam/bam文件中一般不是紧邻的两行,程序会将reads对的第一个比对结果放入内存,直到读取到另一端read的比对结果。因此,选择pos可能会导致程序使用较多的内存,它也适合于未排序的sam/bam文件。而pos排序则表示程序认为双末端测序的reads比对结果在紧邻的两行上,也适合于单端测序的比对结果。很多其它表达量分析软件要求输入的sam/bam文件是按pos排序的,但HTSeq推荐使用name排序,且一般比对软件的默认输出结果也是按name进行排序的。
-s | --stranded     default: yes
  设置是否是链特异性测序。该参数的值可以是yes,no或reverse。no表示非链特异性测序;若是单端测序,yes表示read比对到了基因的正义链上;若是双末端测序,yes表示read1比对到了基因正义链上,read2比对到基因负义链上;reverse表示双末端测序情况下与yes值相反的结果。根据说明文件的理解,一般情况下双末端链特异性测序,该参数的值应该选择reverse(本人暂时没有测试该参数)。
-a | --a     default: 10
  忽略比对质量低于此值的比对结果。在0.5.4版本以前该参数默认值是0。
-t | --type     default: exon
  程序会对该指定的feature(gtf/gff文件第三列)进行表达量计算,而gtf/gff文件中其它的feature都会被忽略。
-i | --idattr     default: gene_id
  设置feature ID是由gtf/gff文件第9列那个标签决定的;若gtf/gff文件多行具有相同的feature ID,则它们来自同一个feature,程序会计算这些features的表达量之和赋给相应的feature ID。
-m | --mode     default: union
  设置表达量计算模式。该参数的值可以有union, intersection-strict and intersection-nonempty。这三种模式的选择请见上面对这3种模式的示意图。从图中可知,对于原核生物,推荐使用intersection-strict模式;对于真核生物,推荐使用union模式。
-o | --samout 
  输出一个sam文件,该sam文件的比对结果中多了一个XF标签,表示该read比对到了某个feature上。
-q | --quiet
  不输出程序运行的状态信息和警告信息。
-h | --help
  输出帮助信息。

3.3 HTSeq使用注意事项

HTSeq的使用有如下注意事项,否则得到的结果是错误的:

1. HTSeq是对有参考基因组的转录组测序数据进行表达量分析的,其输入文件必须有SAM和GTF文件。
2. 一般情况下HTSeq得到的Counts结果会用于下一步不同样品间的基因表达量差异分析,而不是一个样品内部基因的表达量比较。因此,HTSeq设置了-a参数的默认值10,来忽略掉比对到多个位置的reads信息,其结果有利于后续的差异分析。
3. 输入的GTF文件中不能包含可变剪接信息,否则HTSeq会认为每个可变剪接都是单独的基因,导致能比对到多个可变剪接转录本上的reads的计算结果是ambiguous,从而不能计算到基因的count中。即使设置-i参数的值为transcript_id,其结果一样是不准确的,只是得到transcripts的表达量。

3.4 HTSeq的结果

HTSeq将Count结果输出到标准输出,其结果示例如下:

gene00001	0
gene00002	9224
gene00003	880
...
gene12300	1043
gene12301	200
__no_feature	127060
__ambiguous	0
__too_low_aQual	4951
__not_aligned	206135
__alignment_not_unique	0

使用DBG2OLC对二、三代混合数据进行基因组组装

1. DBG2OLC软件简介

DBG2OLC能利用二代和三代混合数据组装大基因组。其文章于2016年发表在Scientific Reports上。

2. DBG2OLC软件下载与安装

使用git下载DBG2OLC软件

$ cd /opt/biosoft/
$ git clone https://github.com/yechengxi/DBG2OLC.git
$ cd /opt/biosoft/DBG2OLC
按照说明中对软件进行编译,编译出的3个可执行程序全部都是DBG2OLC命令
$ g++ -O3 -o SparseAssebmler DBG2OLC.cpp
$ g++ -O3 -o DBG2OLC *.cpp
$ g++ -O3 -o Sparc *.cpp
直接拷贝作者编译好的程序即可
$ chmod 755 compiled/*
$ cp compiled/* .
$ echo 'PATH=$PATH:/opt/biosoft/DBG2OLC' >> ~/.bashrc
$ source ~/.bashrc

DBG2OLC程序第三步需要blasr, sparc/pbdagcon软件支持。其中sparc在DBG2OLC安装文件夹下。
安装blasr

下载BLASR
$ git clone https://github.com/PacificBiosciences/blasr.git /opt/biosoft/blasr
$ cd /opt/biosoft/blasr/
下载libcpp和pbbam两个submodules
$ make update-submodule

blasr编译需要hdf5支持,从hdf5官网下载适合centos6的二进制包并安装 
$ wget https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.10/hdf5-1.10.0-patch1/bin/linux-centos6-x86_64-gcc447/hdf5-1.10.0-patch1-linux-centos6-x86_64-gcc447-shared.tar.gz
$ tar zxf hdf5-1.10.0-patch1-linux-centos6-x86_64-gcc447-shared.tar.gz -C /opt/sysoft/

可以使用cmake, pitchfork和make三种方式对BLASR进行编译,以下使用常规的make方法进行编译,需要高版本gcc支持
对BLASR进行configure
$ ./configure.py --shared --sub --no-pbbam HDF5_INCLUDE=/opt/sysoft/hdf5-1.10.0-patch1-linux-centos6-x86_64-gcc447-shared/include/ HDF5_LIB=/opt/sysoft/hdf5-1.10.0-patch1-linux-centos6-x86_64-gcc447-shared/lib/
对submodules进行configure
$ make configure-submodule
对submodules进行make
$ make build-submodule -j 4
对BLASR进行make
$ make blasr -j 4
对其它工具,例如pls2fasta, loadPulses, sawriter等进行编译,其结果文件在utils文件夹中
$ make -j 4

可选步骤:手动将有用的命令和库文件放置到指定的地方
blasr的正常运行需要依赖libcpp里面三个库文件和hdf5软件中的库文件
$ mkdir bin lib
$ cp blasr bin/
$ find utils -maxdepth 1 -perm 775 -exec cp {} bin/ \;
$ cp ./libcpp/pbdata/libpbdata.so ./libcpp/hdf/libpbihdf.so ./libcpp/alignment/libblasr.so lib/
$ echo 'export LD_LIBRARY_PATH=/opt/sysoft/hdf5-1.10.0-patch1-linux-centos6-x86_64-gcc447-shared/lib/:/opt/biosoft/blasr/lib/:$LD_LIBRARY_PATH
PATH=/opt/biosoft/blasr/bin/:$PATH' >> ~/.bashrc.pacbio
$ source ~/.bashrc.pacbio

若DBG2OLC流程第三步选择使用pbdagcon进行运算,则需要安装pbdagcon软件

pbdagcon软件的编译需要高版本gcc支持
$ git clone https://github.com/PacificBiosciences/pbdagcon.git /opt/biosoft/pbdagcon
$ cd /opt/biosoft/pbdagcon
$ ./configure.py --boost --gtest --sub --no-pbbam
$ make init-submodule
$ make -j 4
$ make check
$ mkdir bin
$ cp src/cpp/dazcon src/cpp/pbdagcon bin/
$ echo 'PATH=/opt/biosoft/pbdagcon/bin:$PATH' >> ~/.bashrc.pacbio
$ source ~/.bashrc.pacbio

3. 程序运行

使用DBG2OLC软件利用二代和三代数据混合的基因组组装,其运行流程分3步。

3.1 使用SparseAssembler利用二代数据进行DBG组装

首先,利用Illumina小片段文库数据使用SparseAssembler命令组装出contigs序列。此外,也可以使用其他基因组组装软件组装出contigs序列后,直接跳到DBG2OLC的第二个步骤。值得注意的是:输入到第二步骤的contigs必须是没有经过repeat resolving的原始序列;绝大部分基因组组装软件为了获得更完整连续的contigs序列,都牺牲了contigs的准确性,其结果不能用于DBG2OLC软件的第二步,否则最终结果会很差;作者推荐可以直接用于第二步的其它contig组装软件有Platanus和Meraculous。
一般情况下,输入到SparseAssembler命令中~50x的Illumina数据,能获得较好的contigs结果。
SparseAssembler命令参数:

常用参数:
LD 
    是否载入k-mer graph。第一次运行SparseAssembler命令的时候,该参数的值必须是0;若为了使用SparseAssembler得到更好的contigs结果,则需要调整NodeCovTh和EdgeCovTh参数;调整这些参数的时候,不需要再次计算k-mer graph,设置该参数为1来跳过这个步骤,从而节约很多时间。
k 
  设置使用DBG方法计算时的Kmer大小,支持的Kmer大小为15-127。
g 
  number of skipped intermediate k-mers, support 1-25.该参数软件示例中使用的值是15。
NodeCovTh 
  设置k-mers覆盖度阈值,去除覆盖度较低的k-mers。该值设定范围为0-16,默认值为1。
EdgeCovTh 
  设置link覆盖度阈值,去除覆盖度较低的links。该值设定范围为0-16,默认值为0。
GS 
  设置一个基因组大小的值。该参数用于决定预先占用的内存量。推荐设置得比基因组大,例如设置为2倍基因组大小。
f 
  输入单端测序数据的路径。输入文件可以是fasta或fastq文件。若有多个输入文件,则使用多个f参数。
i1  i2 
  输入inward paired-end数据。若有多组paired-end数据,则多次使用i1/i2参数。

其它参数:
o1  o2 
  输入outward paired-end数据。
i1_mp  i2_mp 
  输入插入片段长度>10kb的inward paired-end数据。
o1_mp  o2_mp 
  输入插入片段长度>10kb的outward paired-end数据。
PathCovTh 
  设置path覆盖度阈值,去除覆盖度较低的paths。该值设定范围为0-100。根据经验,不推荐添加该参数。
TrimLen 
  将所有过长的序列截短到此指定的长度。
TrimN 
  若read中的碱基N数目超过此设定的值,则去除该read数据。
TrimQual 
  从尾部截短质量低于此值的碱基。
QualBase 
  设置Fastq文件中最低碱基质量对应的ASCII码符号。默认值是'!',等同于Pred33。
Denoise 
  设置是否对reads进行修正。默认值是0,表示不对reads进行修正。
H 
  混合模式。默认值是0,表示对reads的尾部进行截短,来保证高质量的数据进行reads修正。
CovTh 
  覆盖度 < 此设定值的k-mer会被检测,从而被校正。若该参数值设置为0,则软件会自动计算该值。
CorrTh 
  覆盖度 >= 此设定值的k-mer可以用来对reads做校正。若该参数值设置为0,则软件会自动计算该值。

SparseAssembler运行示例:

对某物种Illumina小片段文库测序的PE150bp数据使用trimmomatic质控,再使用FindErrors进行修正,再运行SparseAssembler:
$ SparseAssembler LD 0 k 95 g 15 NodeCovTh 1 EdgeCovTh 0 GS 60000000 f A.1.fastq f A.2.fastq f B.1.fastq f B.2.fastq
$ cp Contigs.txt Contigs.txt.00
增大NodeCovTh和EdgeCovTh参数后,再次运行SparseAssembler,并比较两次结果。第二次运行较第一次运行,耗时少了很多很多。
$ SparseAssembler LD 1 k 95 g 15 NodeCovTh 2 EdgeCovTh 1 GS 60000000 f A.1.fastq f A.2.fastq f B.1.fastq f B.2.fastq

SparseAssembler在当前目录下生成了18个文件结果,其中Contigs.txt文件是Fasta格式的Contigs序列。
运行SparseAssembler的注意事项:

1. SparseAssembler只可以简单地对Fastq文件进行质控和错误修正。推荐使用其它软件进行reads质控和修正,以获得更好的结果。
2. 参数k设置了k-mer的大小,该参数的值对结果的影响较大。若基因组较小,推荐设置多个k-mer值进行多次计算,从而选择最优k-mer值。个人经验,PE150bp数据的最优的k-mer值约为91~99。
3. 选定了k-mer大小后,使用默认的NodeCovTh和EdgeCovTh参数(默认参数一般能得到很好的结果)运行一遍SparseAssembler。然后尝试增大NodeCovTh和EdgeCovTh参数,设置LD 1参数再次运行SparseAssembler,以获得最优的Contigs结果。
4. 可能是先使用了最小的NodeCovTh和EdgeCovTh参数做运算后,才能再次使用更大的参数进行运算。
5. SparseAssembler虽然也有输入大片段文库数据的参数和Scaffolding参数,但是不推荐输入大片段文库数据进行Scaffolding操作,没太大意义。
6. 虽然SparseAssembler命令的文件输入方式有多种,若是仅进行contigs组装,没有利用到paired信息,因此使用i1 i2参数输入文件和使用f参数输入文件的结果是一模一样的。

3.2 使用DBG2OLC找Contigs序列和Pacbio reads的Overlap并进行Layout

DBG2OLC通过比较contigs和Pacbio reads之间的overlap,将contigs序列定位到Pacbio reads上,将DBG的contigs结果运用到OLC算法中。
DBG2OLC命令参数:

主要参数:
LD 
  是否载入compressed reads information。第一次运行DBG2OLC命令的时候,该参数的值必须是0;若为了得到更好的结果,则需要调整其它参数;调整这些参数的时候,设置该参数为1来跳过这个步骤,从而节约很多时间。
k 
  设置k-mer大小。k-mer用来比较contig和pacbio read之间的重叠,而不是用于基因组组装,推荐设置为 17 即可。
AdaptiveTh 
  若contig和pacbio read之间匹配的k-mers数目 < AdaptiveTh * contig长度,则认为contig和pacbio read没有重叠。推荐设置为0.001-0.02。
KmerCovTh 
  若contig和pacbio read之间匹配k-mers的覆盖度 < KmerCovTh,则认为contig和pacbio read没有重叠。推荐设置为2-10。
MinOverlap 
  两条Pacbio read之间匹配的k-mers数目 < MinOverlap,则认为它们之间没有重叠。推荐设置为10-150。
RemoveChimera 
  去除嵌合体Pacbio reads。若Pacbio数据覆盖度大于10x,推荐设置该参数为 1 。
Contigs 
  输入contigs序列文件路径
f 
  输入Pacbio测序Fasta/Fastq文件路径。

其它参数:
MinLen 
  设置能用于计算的最小Pacbio reads长度。
ChimeraTh 
  该参数默认值是 1 ;若Pacbio数据覆盖度大于100x,则推荐加入该参数并设置为 2 。
ContigTh 
  该参数默认值是 1 ;若Pacbio数据覆盖度大于100x,则推荐加入该参数并设置为 2 。

DBG2OLC运行示例:

$ DBG2OLC LD 0 k 17 AdaptiveTh 0.001 KmerCovTh 2 MinOverlap 20 RemoveChimera 1 Contigs Contigs.txt f Pacbio_Cell01.fastq f Pacbio_Cell02.fastq
$ DBG2OLC LD 1 k 17 AdaptiveTh 0.005 KmerCovTh 3 MinOverlap 30 RemoveChimera 1 Contigs Contigs.txt f Pacbio_Cell01.fastq f Pacbio_Cell02.fastq

DBG2OLC的结果文件很多,其主要结果文件是backbone_raw.fasta和DBG2OLC_Consensus_info.txt,是第三步的输入文件。
运行DBG2OLC的注意事项:

1. AdaptiveTh, KmerCovTh和minOverlap这3个计算Overlap的参数对结果的影响最大。对于10x/20x PacBio数据:KmerCovTh 2-5, MinOverlap 10-30, AdaptiveTh 0.001~0.01;对于50x-100x PacBio数据:KmerCovTh 2-10, MinOverlap 50-150, AdaptiveTh 0.01-0.02。
2. 不推荐对Pacbio数据就行修正后再运行DBG2OLC。可以比较使用修正前的数据用于DBG2OLC的结果,一般情况下使用未修正的Pacbio数据能获得更好的结果。此外,DBG2OLC运行过程中有一步多序列比对模块来进行错误修正。
3. 可能是先使用了最小的AdaptiveTh, KmerCovTh和minOverlap参数做运算后,才能再次使用更大的参数进行运算。

3.3 Call consensus

本步骤需要使用/opt/biosoft/DBG2OLC/utility/目录下的python和shell脚本,来调用blasr和consensus模块Sparc(可以考虑使用pbdagcon作为consensus模块,但DBG2OLC没有提供相应的脚本)进行运算。

先将/opt/biosoft/DBG2OLC/utility/目录下的python和shell脚本拷贝到当前目录
$ cp /opt/biosoft/DBG2OLC/utility/*.sh /opt/biosoft/DBG2OLC/utility/*.py ./
若使用了最新版本的blasr软件,其参数书写方法有一个中划线变成了两个中划线,因此需要修改.sh文件中blasr命令的参数书写方法。
此外,也需要修改.sh文件中Sparc命令路径,或者将Sparc命令也拷贝到当前目录。

将Contigs序列和Pacbio reads数据合并成一个文件ctg_pb.fasta
$ cp Contigs.txt ctg_pb.fasta
$ perl -e 'while (<>) {s/^\@/>/; print; $_ = <>; print; <>; <>;}' Pacbio_Cell01.fastq >> ctg_pb.fasta
$ perl -e 'while (<>) {s/^\@/>/; print; $_ = <>; print; <>; <>;}' Pacbio_Cell02.fastq >> ctg_pb.fasta

运行脚本程序split_and_run_sparc.sh
$ ./split_and_run_sparc.sh backbone_raw.fasta DBG2OLC_Consensus_info.txt ctg_pb.fasta ./ 2 > cns_log.txt
结果会输出到 ./ 目录下。最后的结果文件是final_assembly.fasta。