简介：HTseq 是用于Python平台写成的处理高通量测序的平台。htseq-count可以用来对原始转录本计数，具体计数规则参考Counting reads in features with htseq-count。
平台：Python跨平台使用。
快速运行：
# 查看帮助
$ htseq-count --help
# 对bam文件按照reads名称排序
$ samtools sort -n accepted_filtered.bam accepted_sortname
# 注意“=”前后无空格
$ htseq-count --mode=union --stranded=no --type=exon --idattr=gene_id \
--format=bam accepted_sort.bam hg19USCS_ensembl.gtf > htseqcount_accepted.hsc
重要参数解释：
--mode：统计落在某个基因上的reads数目的模型，默认值为“union”（图示）。作者认为“union”方法在绝大多数情况下都有很好的表现，建议使用。
--stranded：测序方法，默认为“yes”。
--type：计数单元类型，默认为GTF文件的exon。
--idattr：计数单元归类，默认为gene_id。比如把合并汇报多个exon对应的一个gene。
--format：可以输入bam或者sam文件，bam文件需要制定此参数。
补充：
$ samtools sort -n accepted_filerted.bam acceted_sortname
GenomicAlignments是R/Bioconductor的一个包，其中summarizeOverlaps函数用于对alignment数据进行计数。
简介：kallisto是由Lior Pachter小组开发一款快速测量RNA-seq数据中转录本表达丰度的软件。因为使用了pseudoalignment的想法，可以不用alignment，直接测量原始测序数据，因此极大提高了运算速度。因为kallisto极高的速度，可以使用bootstrap精确估计的“不确定性（uncertainty）”，可以配合下游软件sleuth确定差异表达基因。
平台：Linux和Mac
快速运行：
首先，建立索引。
# 查看帮助
$ kallisto index
# 建立索引
$ kallisto index -i hg19.idx hg19UCSC_ensembl.fa
重要参数解释：
-i/--index：kallisto输出索引文件。
-k/--kmer-size：k-mer的长度（奇数），最大值为31，默认值为31。之后，根据原始RNA-seq数据和转录组注释，测量表达量。
$ kallisto quant -i hg19.idx -o kallOut -b 1000 human1.fq.gz human2.fq.gz
重要参数解释：
-i/--index：kallisto输出索引文件。
-o/--output-dir：输出目录。
-b/--bootstrap-samples：bootstrap的次数，默认为0。
简介：edgeR是使用RNA-seq的reads的原始计数（count）和负二项分布模型，计算差异表达基因的R/Bioconductor包。edgeR能够处理RNA-seq、ChIP-seq和SAGE数据。根据edgeR的说明文件数据前处理有两点要求：
edgeR performs differential abundance analysis for pre-defined genomic features. Although not strictly necessary, it usually desirable that these genomic features are non-overlapping.（edgeR对预定义的基因组特征数据进行差异量分析。尽管不是严格必须，但是输入的基因特征最好能无相互重叠。）
The first step in an RNA-seq analysis is usually to align the raw sequence reads to a reference genome, although there are many variations on this process.（RNA-seq分析的第一步通常是把原始测序reads比对到参考基因组，当然，这一步有多种实现方法。）
以上两点表明输入的数据：1. 分析某个基因组特征，比如基因、外显子、CDS；2. 原始reads的计数，而不是一个估计量、FPKM或者RPKM，这些估计量会影响edgeR计算一些模型涉及的关键变量。对于第二点，edgeR指出也能处理类似RSEM输出的计数估计值，但原始计数最好。
简介：DESeq2是R/Bioconductor的包，使用负二项分布模型和RNA-seq的原始count数值，寻找差异表达基因。DESeq2说明文件指出：
平台：R跨平台使用。
运行流程：
# 载入包并且注册多核
library(DESeq2)
library(BiocParallel)
register(MulticoreParam(4))
# “glioPR”是一个整理好的数据框，行表示count数，列标识样本（1~10列）和注释（11~14列）。
# “targets”是一个数据框，表示样本的分类，注意使用因子型（factor）标记不同组
> targets
Treatment Sample Patient
human1 primary human1 1
human2 recurrent human2 1
human3 primary human3 2
human4 recurrent human4 2
human5 primary human5 3
human6 recurrent human6 3
human7 recurrent human7 4
human8 recurrent human8 5
human9 primary human9 6
human10 primary human10 7
# 创建DESeqDataSet对象，使用Treatment作为比较
glioPR <- DESeqDataSetFromMatrix(countData = htCountSelect[, 1:10], colData = targets, design = ~ Treatment)
mcols(glioPR) <- htCountSelect[, 11:14]
# 创建DESeqDataSet对象，使用Treatment和Pateint两个因子作为比较
glioPR <- DESeqDataSetFromMatrix(countData = htCountSelect[, 1:10], colData = targets, design = ~ Patient + Treatment)
mcols(glioPR) <- htCountSelect[, 11:14]
# DEGs
glioPR <- DESeq(glioPR)
res <- results(glioPR)
## 根据FDR值排序
res[order(res$padj), ]
summary(res)
# 两种标准化count值方法， rlog（regularized log）和VST（variance stabilizing transformations）
rld <- rlog(glioPR)
vsd <- varianceStabilizingTransformation(glioPR)
rlogMat <- assay(rld)
vstMat <- assay(vsd)
2015年6月5日