前言
这次给大家带来的是16年发表在NATURE PROTOCOLS上面的一篇处理RNA-seq数据的文章：Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown
和它的12年发表于同一个杂志的兄弟文章：Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks都是NATURE PROTOCOLS上阅读量最大的文章。
NATURE PROTOCOLS
当然还有很多其它的介绍不止生信的方法文章，大家有时间可以去探究。
同时我这里就不再赘述RNA-seq的具体原理，有需要了解的请移步：一个简略的RNA-seq演示
至于软件的安装到官网下载，解压后将bin/添加进路径即可，这里不再做讲解。
#######所有操作皆在LINUX&R上完成，默认基本处理软件已经安装
本体介绍
大佬的文章
事实上作者团队一直致力于开发出更好的解决数据处理的软件，就目前12年推出的Tohat2和Cufflinks已经不是那么的令人满意了，所以他们又开发了 HISAT, StringTie and Ballgown三件套完成大家对于一个RNA-seq所有的幻想。#但实际上目前还存在着其他的性能很优秀软件也可以满足需求，例如STAR等等，但作为菜鸟来说只有先一步一步的学习了。
好了，我们先对这篇文章给出的处理方法有个大概的了解。
pipeline：
- alignment of the reads to the genome
- assembly the alignments into full-length transcripts
- quantification of the differencesin expression levels of each gene and transcript
- calculation of the differences in expression for all genes among the different conditions
An overview of the 'new Tuxedo' protocol
这个protocol首先从原始RAN-seq数据入手，输出数据包括基因list，转录本，及每个样本的表达量，能够表现差异表达基因的表格并完成显著性的计算。
- 第一步是使用HISAT将读段匹配到参考基因组上,使用者可以提供注释文件，但HISAT依旧会检测注释文件没有列出来的剪切位点
- 第二步，比对上的reads将会被呈递给StringTie进行转录本组装，StringTie单独的对每个样本进行组装，在组装的过程中顺带估算每个基因及isoform的表达水平
- 第三步，所有的转录本都被呈递给StringTie的merge函数进行merge，这一步是必须的，因为有些样本的转录本可能仅仅被部分reads覆盖，无法被第二步的StringTie组装出来。merge步骤可以创建出所有样本里面都有的转录本，方便下一步的对比
- merge的数据再一次被呈递给StringTie，StringTie可以利用merge的数据重新估算转录本的丰度，还能额外的提供转录本reads数量的数据给下一步的ballgown
- 最后一步：Ballgown从上一步获得所有转录本及其丰度，根据实验条件进行分类统计
注意事项Warning：
- HISAT, StringTie, Ballgown并不是唯一处理RNA-seq的方法，也并不能处理所有的RNA-seq数据。
- 例如质量控制或者去除测序污染/去除接头/去除低质量数据（可以使用FastQC以及FASTXtoolkit做到）
- 这个pipeline无法处理第三代测序的数据
- Ballgwon可以使用stringtie/cufflinks/RSEM产出的数据，但是可能无法接受其它程序产出的结果
- 默认参数适用于小至三个，大至小数百的样本处理，对于特殊需求还需要参考manual
数据处理设计
Read alignment with HISAT
总体上来说HISAT利用了BWA和Bowtie的算法，同时解决了mRNA中不存在内含子难以比对的问题，比对上代主流RNA-seq比对工具能快50倍，同时需求更少的内存<8G（这就意味着你可以在PC上跑数据），20个样本，每个样本一亿reads的估计，用一台电脑一天之内能够跑完。使用者可以提供精确的基因注释来提高在已知基因区域的准确性，但这是可选项。
事实上，针对RNA-seq数据的align目前有很多工具，16年12月12日出版的NATURE METHODS上的一篇文章，利用分别使用人和一种叫malaria的寄生虫的数据，模拟出三种复杂度的数据集（T1、T2和T3）。对于复杂度的衡量，定义了三个尺度：
- 突变率
- indel比例
- 错误率
T1复杂度最低，T3最高。
随后使用了目前主流的比对工具进行了比对。
因为RNA-Seq测序的特性，天然的会有一部分数据延伸到内含子区，这部分跨越外显子和内含子的reads就称为『junction reads』，所以RNA-Seq比对软件需要针对此进行优化，而文章做benchmark也考虑到这点。分两个水平做比较： - base and read level，这点和一般的DNA aligner一样
- junction-level
度量软件的结果时，用了两个值： - recall： 所有base中正确比对的比例。Tophat2在T1 base&read的表现，recall是85%-95%，T2降到80%，T3更是雪崩式降至20%以下。这部分表现好的是Novoalign和MapSplice2。
-
precision：所有比对上的base中，正确比对的比例。对于T1和T2，大部分软件这个值都较高。
结果上图：
接下来看看调参（自行设置参数，而不是使用软件默认参数）比对，直接说结论吧，不管是否调参，表现robust的是CLC，GSNAP、Novoalign、STAR以及MapSplice2。针对复杂度高的T3数据，Tophat2调参和默认参数得到的比对率相差67%还多。
另外还有运行时间的比较，这轮表现好的是HISAT/HISAT2，其实也能看出，数据复杂度越高，比对耗时越长。
Transcript assembly and quantification with StringTie
RNA-seq的分析依赖于精准的对于基因isoform的重建以及对于基因相对丰度的预测。继承于Cufflinks，StringTie相对于之前开发的工具更为准确，需求内存和耗时也更少。
使用者一样可以使用注释文件来帮助StringTie运行，对于低丰度的数据比较有帮助。此时StringTie依旧会对非注释区域进行转录本的组装，所以注释文件在这里是可选选项。
转录本组装完成后，使用者可以利用gffcompare(StringTie工具包含)工具来得知有多少转录本和注释文件相同，又有多少新的转录本:
#input: gff or gtf format
transcripts.gtf
command line example:
$ gffcompare -G -r chrX.gtf transcripts.gtf
#-r : reference
#-G :compare all transcripts
#-o :prefix
output file 1: gffcmp.annotated.gtf
# 显示StringTie组装的转录本与注释文件内的转录本的差别（会给每个转录本加入一个class code，我理解为一个标识，释义如下图）
output dile 2: gffcmp.stats
# 根据注释文件显示组装结果的准确性和阳性预测率
class codes
由于在实验中，我们会处理多个RNA-seq数据，单个样本里面的基因和iosforms与其它样本的数据很少相同，但是为了进行比较就需要它们以一个相同的形式进行组装，作者通过StringTie的merge工具解决了这个问题，将所有样本中发现的基因进行merge。下图的例子可以帮助理解。
example
如图所示，四个样本中组装得到四个转录本，最后merge得到A/B两个转录本。1/2样本与参考基因组相同merge得到A转录本，3/4样本相互吻合但与参考基因组不一样，所以merge得到B转录本。
在merge步骤之后，需要StringTie再运行一遍去重新估算merge得到的转录本的丰度。
Chevy理解：这就是相当于重新定义了一个annotation file进行二次重新组装，相对于第一次组织可以提高准确度和敏感性。比已有注释文件的优势在于可以发现更多的isoforms。
Differential expression analysis with Ballgown
组装和定量完转录本之后自然需要进行基因表达差异分析，统计建模和可视化。R和Bioconductor提供了一揽子的工具处理这些任务，包括raw data的plot作图，数据的标准化，下游的统计建。此处使用Ballgown包作为承上启下的桥梁。
在R里面使用Ballgown处理需要得到：
- 表型数据。关于样本的信息
- 表达数据。标准化和未标准化的关于外显子，junction，转录本，基因的表达数量
- 基因信息。有关外显子，junction，转录本，基因的坐标以及注释信息
而大多数的差异表达分析都会包括下面几个步骤：
- 数据可视化和检查
- 差异表达的统计分析
- 多重检验校正
- 下游检查和数据summary
Ballgown包可以完成以上的几个步骤并且可以联合R语言的其它操作进行分析
Ballgown使用：
- 数据的读入：需要将StringTie输出的数据结合表型数据，这里需要保证表型数据的identifier和StringTie输出数据的一致，不然会报错
- 预测丰度的检查：以FPKM（fragments per kilobase of transcript per million mapped reads）为单位的丰度预测将会根据library size进行标准化。#极端差异此处需要引起注意
- 使用线性模型进行差异表达分析，由于FPKM对于转录本解读过于曲解，所以这里需要使用log转化处理数据，随后再使用线性模型进行差异分析。
- Ballgown可以对time-course和fixed-condition数据进行差异分析，但是无法避免批次效应带来的误差。#使用stattest功能可以指定任何已知的cofounder
- Ballgown可以帮助你在基因，转录本，外显子，junction水平上进行差异分析
- 结果会以表格形式展现，如果样本够多还有p值和q值
- 结果数据可以用来进行下一步的gene set分析（例如GSEA）等等
实战示例
准备阶段
实际上本轮操作就是按照文章给的示例数据走了一遍流程：
文章中，为了减少计算时间，方便读者重复，作者截取了一批实验数据中能够匹配到chromosomeX上的数据作为示例数据，chromosomeX是一个基因相对较为丰富的染色体，可以占到基因组的5%左右。
首先是下载数据随后解压：
$nohup wget ftp://ftp.ccb.jhu.edu/pub/RNAseq_protocol/chrX_data.tar.gz 2>download.log &
$tar zxvf chrX_data.tar.gz
其中包括如下数据：
文件数据
sample文件夹下有12个fastq格式的paired-end RNA-seq文件，从两地人群中（英格兰岛住民和约鲁巴住民）各取六个样，六个样又分为男女性别各三个（最少的生物学重复）。
sample
随后介绍一个骚命令：
HISAT2可以直接从NCBI下载sra格式的文件并作为输入文件格式
下面以ERR188245测序数据为例
$ hisat2 -p 8 --dta -x chrX_data/indexes/chrX_tran --sra-acc ERR188245 –S ERR188245_chrX.sam
可以说是相当的帮人偷懒了
index文件夹下包含着HISAT2对于染色体X的index文件
index
当然HISAT2已经为各位老爷准备好了常见基因组的index文件和genes文件。
点我
genome文件夹下包含这染色体X的序列文件（GRCh38 build 81）
genome
genes文件夹下则包含着针对基因组的注释文件（来自于RefSeq数据库）
annotation
geuvadis_phenodata.csv和mergelist.txt则是用户着要自己创建表明数据关系的文件，这里作者提供出来方便使用。
如果需要建立官网上没有基因组的index，就需要需要使用Hisat2包里面的python脚本:
extract_splice_sites.py和extract_exons.py，从注释文件里面抽取出剪切位点和外显子信息
$ extract_splice_sites.py chrX_data/genes/chrX.gtf >chrX.ss
$ extract_exons.py chrX_data/genes/chrX.gtf >chrX.exon
先提取剪切位点和外显子数据
$ hisat2-build --ss chrX.ss --exon chrX.exon chrX_data/genome/chrX.fa chrX_tran
随后建立HISAT2 index
正式开始
1.将reads比对上genome
2017-07-24更新：我自己使用的是Ensembl发表的GRCh38版本的基因组进行比对，所以这里我附上这个版本的genome和注释文件的下载地址：
Ensembl版本全基因组的注释文件下载
HISAT2-index下载genome_tran
# 样本比对操作
nohup hisat2 -p 8 --dta -x ~/RNA-seq/chrx/chrX_data/indexes/chrX_tran -1 ~/RNA-seq/chrx/chrX_data/samples/ERRERR188044_chrX_1.fastq.gz -2 ERR188044_chrX_2.fastq.gz -S ERR188044_chrX.sam &
# 数据太多我就写了个小脚本处理,in sample dir/
for i in *1.fastq.gz;
do
i=${i%1.fastq.gz*};
nohup hisat2 -p 8 --dta -x ~/RNA-seq/chrx/chrX_data/indexes/chrX_tran -1 ${i}1.fastq.gz -2 ${i}2.fastq.gz -S ~/RNA-seq/chrx/chrX_data/result/${i}align.sam 2>~/RNA-seq/chrx/chrX_data/result/${i}align.log &
done
运行结果如下：
比对结果
2.将sam文件转化为bam文件
# 转化操作
samtools sort @ 8 -o ERR188044_chrX.bam ERR188044_chrX.sam
# 批处理, in result dir/
for i in *.sam;
do
i=${i%_align.sam*}; nohup samtools sort -@ 8 -o ${i}.bam ${i}_align.sam &
done
结果如下：
sam to bam
3.组装转录本并定量表达基因
# 单文件操作
stringtie -p 8 -G ~/RNA-seq/chrx/chrX_data/genes/chrX.gtf -o ERR188044_chrX.gtf ERR188044_chrX.bam
# 批处理, in result dir/
for i in *.bam;
do
i=${i%.bam*}; nohup stringtie -p 8 -G ~/RNA-seq/chrx/chrX_data/genes/chrX.gtf -o ${i}.gtf ${i}.bam &
done
结果如下:
strigtie
4.将所有转录本合并
warning: 此处的mergelist.txt是自己创建的，需要包含之前output.gtf文件的路径
cat ~/RNA-seq/chrx/chrX_data/mergelist.txt
ERR188044_chrX.gtf
ERR188104_chrX.gtf
ERR188234_chrX.gtf
ERR188245_chrX.gtf
ERR188257_chrX.gtf
ERR188273_chrX.gtf
ERR188337_chrX.gtf
ERR188383_chrX.gtf
ERR188401_chrX.gtf
ERR188428_chrX.gtf
ERR188454_chrX.gtf
ERR204916_chrX.gtf
#因为就在结果目录下面操作，所以直接显示文件名即可
stringtie --merge -p 8 -G ~/RNA-seq/chrx/chrX_data/genes/chrX.gtf -o stringtie_merged.gtf ~/RNA-seq/chrx/chrX_data/mergelist.txt
#output文件即为
stringtie_merged.gtf
5.检测相对于注释基因组，转录本的组装情况
gffcompare -r ~/RNA-seq/chrx/chrX_data/genes/chrX.gtf -G -o merged stringtie_merged.gtf
#输出文件信息可见上面的Transcript assembly and quantification with StringTie介绍
6.重新组装转录本并估算基因表达丰度，并为ballgown创建读入文件
# 单文件操作
stringtie -e -B -p 8 -G stringtie_merged.gtf -o ballgown/ERR188044_chrX/ERR188044_chrX.gtf ERR188044_chrX.bam
# 批处理, in result dir/
for i in *_chrX.bam;
do
i=${i%_chrX.bam*}; nohup stringtie -e -B -p 8 -G stringtie_merged.gtf -o ballgown/${i}/${i}_chrX.gtf ${i}_chrX.bam &
done
结果如下：
table
7.R包的安装与载入
R语言需要安装ballgown包和一些接下来处理会用到的包（RSkittleBrewer/genefilter/dplyr/devtools）
事实上我还没搞清楚LINUX里面的R包安装问题，老是报error，所以这里我就用Windows下的R进行操作。
# 首先是几个包的载入
>library(ballgown)
>library(RSkittleBrewer)
>library(genefilter)
>library(dplyr)
>library(devtools)
8.数据的读入
# 随后是数据的读入，CSV文件，我把所有文件放在桌面，上一步得到的ballgown文件夹直接拷到桌面
> setwd("C:/Users/DELL/Desktop")
> read.csv("geuvadis_phenodata.csv")
ids sex population
1 ERR188044 male YRI
2 ERR188104 male YRI
3 ERR188234 female YRI
4 ERR188245 female GBR
5 ERR188257 male GBR
6 ERR188273 female YRI
7 ERR188337 female GBR
8 ERR188383 male GBR
9 ERR188401 male GBR
10 ERR188428 female GBR
11 ERR188454 male YRI
12 ERR204916 female YRI
> pheno_data<-read.csv("geuvadis_phenodata.csv")
> bg_chrX=ballgown(dataDir = "C:/Users/DELL/Desktop/ballgown",samplePattern = "ERR",pData = pheno_data)
# dataDir指的是数据储存的地方，samplePattern则依据样本的名字来，pheno_data则指明了样本数据的关系，这个里面第一列样本名需要和ballgown下面的文件夹的样本名一样，不然会报错
10.滤掉低丰度的基因
# 这里滤掉了样本间差异少于一个转录本的数据
> bg_chrX_filt=subset(bg_chrX,"rowVars(texpr(bg_chrX))>1",genomesubset=TRUE)
11.确认组间有差异的转录本
2017-07-21更新：实际上为了确定基因表达差异是由于某个变量造成的，我们这边需要使用剩下的变量进行修正。Example：我们为了比较male and female之间的基因表达差异，这里就需要指定分析参数为"transcripts"，主变量为"sex"，修正变量为"population"，getFC可以指定输出结果显示组间表达量的foldchange。
并且Ballgown的统计算法基于标准的线性模型，对于组内数据较少(<4)时较为准确。这里也可以使用其它的一些软件例如limma/DESeq/edgeR等。
> result_transcripts=stattest(bg_chrX_filt,feature = "transcript",covariate = "sex",adjustvars = c("population"), getFC=TRUE,meas="FPKM")
> result_transcripts
feature id fc pval qval
1 transcript 1 0.941367186 7.353184e-01 9.468599e-01
2 transcript 2 1.207949354 8.668638e-01 9.744055e-01
3 transcript 3 1.013100019 9.920824e-01 9.986659e-01
4 transcript 4 0.387671042 5.233364e-01 9.189906e-01
5 transcript 5 0.615797877 3.324164e-01 9.117015e-01
......
12.确认组间有差异的基因
这里指定feature为gene
> result_genes=stattest(bg_chrX_filt,feature = "gene",covariate = "sex",adjustvars = c("population"), getFC=TRUE,meas="FPKM")
> result_genes
feature id fc pval qval
1 gene MSTRG.1 1.114999374 7.305975e-01 9.218157e-01
2 gene MSTRG.10 1.749747142 2.767538e-01 7.922755e-01
3 gene MSTRG.1000 1.399358968 6.039065e-01 8.945639e-01
4 gene MSTRG.1002 0.937668743 8.825216e-01 9.816878e-01
5 gene MSTRG.1003 1.044840944 6.155345e-01 8.977043e-01
6 gene MSTRG.1004 1.220626116 4.160214e-01 8.388688e-01
.....
13.对结果 result_transcripts增加基因名
> result_transcripts=data.frame(geneNames=ballgown::geneNames(bg_chrX_filt),geneIDs=ballgown::geneI Ds(bg_chrX_filt),result_transcripts)
> result_transcripts
geneNames geneIDs feature id fc pval qval
1 . MSTRG.4 transcript 1 0.941367186 7.353184e-01 9.468599e-01
2 PLCXD1 MSTRG.4 transcript 2 1.207949354 8.668638e-01 9.744055e-01
3 . MSTRG.4 transcript 3 1.013100019 9.920824e-01 9.986659e-01
4 . MSTRG.4 transcript 4 0.387671042 5.233364e-01 9.189906e-01
5 . MSTRG.5 transcript 5 0.615797877 3.324164e-01 9.117015e-01
6 PLCXD1 MSTRG.4 transcript 6 0.630018786 2.938396e-01 9.002815e-01
......
14.按照P值排序（从小到大）
> result_transcripts=arrange(result_transcripts,pval)
> result_genes=arrange(result_genes,pval)
15.把结果写入csv文件
> write.csv(result_transcripts, "chrX_transcript_results.csv",row.names=FALSE)
> write.csv(result_genes, "chrX_gene_results.csv",row.names=FALSE)
16.筛选出q值小于0.05的transcripts和genes，也就是在性别间表达有差异的基因了。
注：这里无需过分关注geneIDs以及transcript name，实际上那个是stringtie在比对过程中赋予的一个符号。
> subset(result_transcripts,result_transcripts$qval<0.05)
geneNames geneIDs feature id fc pval qval
1 . MSTRG.531 transcript 1657 0.030820591 1.427864e-10 2.082377e-07
2 XIST MSTRG.531 transcript 1656 0.003014576 1.860927e-10 2.082377e-07
3 . MSTRG.531 transcript 1655 0.016144997 3.762791e-10 2.807042e-07
4 . MSTRG.531 transcript 1658 0.028308029 6.599039e-08 3.692163e-05
5 TSIX MSTRG.530 transcript 1654 0.078461818 1.690716e-06 7.567643e-04
6 . MSTRG.613 transcript 1848 7.391342987 1.249275e-05 4.659796e-03
7 . MSTRG.141 transcript 421 3.204058932 2.729898e-05 8.727874e-03
8 . MSTRG.618 transcript 1852 9.136857610 4.804244e-05 1.343987e-02
9 . MSTRG.777 transcript 2338 0.127674288 1.000751e-04 2.488533e-02
10 KDM6A MSTRG.258 transcript 736 0.054212867 1.173824e-04 2.627018e-02
11 PNPLA4 MSTRG.62 transcript 186 0.592778584 2.083496e-04 4.238966e-02
12 RPS4X MSTRG.519 transcript 1611 0.597532121 2.493976e-04 4.651265e-02
> subset(result_genes,result_genes$qval<0.05)
feature id fc pval qval
1 gene MSTRG.531 0.002396124 2.469125e-11 2.523446e-08
2 gene MSTRG.141 3.165966412 1.666206e-06 8.514310e-04
3 gene MSTRG.530 0.082749314 7.018439e-06 2.390948e-03
4 gene MSTRG.613 7.314423295 1.128589e-05 2.883544e-03
5 gene MSTRG.618 9.050399157 5.022017e-05 1.026500e-02
6 gene MSTRG.519 0.637382385 6.953432e-05 1.184401e-02
7 gene MSTRG.376 0.620693674 1.134561e-04 1.656458e-02
8 gene MSTRG.62 0.600686117 1.638615e-04 2.093330e-02
9 gene MSTRG.777 0.159178288 2.177206e-04 2.472338e-02
10 gene MSTRG.622 7.870447732 3.737270e-04 3.527295e-02
11 gene MSTRG.778 0.127712244 3.796501e-04 3.527295e-02
12 gene MSTRG.229 1.412379185 4.200674e-04 3.577574e-02
13 gene MSTRG.48 4.158987103 5.279898e-04 4.150812e-02
17.数据可视化之颜色设定
# 赋予调色板五个指定颜色
> tropical= c('darkorange', 'dodgerblue','hotpink', 'limegreen', 'yellow')
> palette(tropical)
# 当然rainbow()函数也可以完成这个任务
> palette(rainbow(5))
18.以FPKM为参考值作图，以性别作为区分条件
# 提取FPKM值
> fpkm = texpr(bg_chrX,meas="FPKM")
#方便作图将其log转换，+1是为了避免出现log2(0)的情况
> fpkm = log2(fpkm+1)
# 作图
> boxplot(fpkm,col=as.numeric(pheno_data$sex),las=2,ylab='log2(FPKM+1)')
结果如下：
FPKM, male in blue, females in orange
19.就单个转录本的查看在样品中的分布
> ballgown::transcriptNames(bg_chrX)[12]
12
"NR_027232"
> ballgown::geneNames(bg_chrX)[12]
12
"LINC00685"
# 绘制箱线图
> plot(fpkm[12,] ~ pheno_data$sex, border=c(1,2),
+ main=paste(ballgown::geneNames(bg_chrX)[12],' : ',
+ ballgown::transcriptNames(bg_chrX)[12]),pch=19, xlab="Sex",
+ ylab='log2(FPKM+1)')
> points(fpkm[12,] ~ jitter(as.numeric(pheno_data$sex)),
+ col=as.numeric(pheno_data$sex))
结果如下：
transcript
20.查看某一基因位置上所有的转录本
# plotTranscripts函数可以根据指定基因的id画出在特定区段的转录本
# 可以通过sample函数指定看在某个样本中的表达情况，这里选用id=1750, sample=ERR188234
> plotTranscripts(ballgown::geneIDs(bg_chrX)[1750], bg_chrX, main=c('Gene MSTRG.575 in sample ERR188234'), sample=c('ERR188234'))
结果如下：
21.以性别为区分查看基因表达情况
# 这里以id=575的基因为例（对应上一步作图）
> plotMeans('MSTRG.575', bg_chrX_filt,groupvar="sex",legend=FALSE)
结果如下：
MSTRG.575
临终讨论
-
关于时间的使用：针对chrX的数据量的比对和组装，在PC上可以40min内搞定，可以说是比较快的了。作者论述如果将12个样本的全数据下下来做分析的话，8核&8G内存的配置花费12.5h可以搞定比对，花费8h可以搞定组装和表达分析。
-
RNA_seq的比对情况一览
-
转录组组装情况一览
-
其实这篇文章说到底还是在讲方法论，它只负责上游的数据处理到可以分析的一步，下游的分析和可视化还需要依赖自身的实验设计和研究目的。
-
官方给出的方案是：stringtie处理得到的数据是直接呈递给ballgown进行差异分析和可视化作图，但是我还不是很清楚输出文件是否可以无差别的被其它软件使用。（我也是菜鸟）
-
RNA-seq的文章其实很多，也不是非要走这个pipeline。另外我觉得ballgown的作图其实也没有很elegant，org
参考文献
Pertea M, Kim D, Pertea G M, et al. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown[J]. Nature protocols, 2016, 11(9): 1650-1667.
日常Bob镇楼
网友评论