流程
将RNA-seq reads与基因组进行比对 ● 时间<20分钟
- 将每个样本的reads映射到参考基因组上
1 | $ hisat2 -p 8 --dta -x chrX_data/indexes/chrX_tran -1 |
- 将SAM文件分类并转换为BAM文件。
1 | $ samtools sort -@ 8 -o ERR188044_chrX.bam ERR188044_chrX.sam |
上述命令适用于SAMtools
1.3版或更新的版本。对于旧版本的SAMtools,用户应该参考Box3。
- Box3 | 使用SAMtools 1.2.1版及以上版本对SAM文件进行分类和转换
最近版本的SAMtools可以在一个步骤中把SAM文件转换成BAM文件。对于使用旧版本的人来说,需要两个步骤。如图所示。
1 | Convert SAM files to binary BAM files: |
组装和量化表达的基因和转录物 ● 时间~15分钟
- 为每个样本组装转录本
1 | $ stringtie -p 8 -G chrX_data/genes/chrX.gtf -o |
- 合并所有样本的转录本
1 | $ stringtie --merge -p 8 -G chrX_data/genes/chrX.gtf -o stringtie_merged.gtf chrX_data/mergelist.txt |
- 错误排除
这里mergelist.txt是一个文本文件,其中有上一步创建的基因转移格式(GTF)文件的名称,每个文件的名称在一行中(参见ChrX_data的例子文件)。如果你不在所有GTF文件所在的同一目录下运行StringTie,那么你可能需要在mergelist.txt中的每个GTF文件名中包括其完整路径。使用任何文本编辑器来创建或编辑(纯文本)mergelist.txt文件。
- 检查转录本与参考注释的对比情况(可选)
1 | $ gffcompare –r chrX_data/genes/chrX.gtf –G –o merged stringtie_merged.gtf |
选项-o指定了gffcompare将创建的输出文件的前缀。上面的命令将生成多个文件,在gffcompare文档中作了解释;更多细节见BOX1。
- 估计转录物丰度,并为Ballgown创建表格计数
1 | $ stringtie –e –B -p 8 -G stringtie_merged.gtf -o |
运行差异化表达分析协议 ● 时间~5分钟
- 加载相关的R包。这些包包括你将用于执行大部分分析的Ballgown包,以及其他一些有助于这些分析的包,特别是RSkittleBrewer(用于设置颜色)、genefilter(用于快速计算平均值和变异值)、dplyr(用于排序和排列结果)和devtools(用于重现性和安装软件包)。
1 | $ R |
- 加载样本的表型数据。一个名为geuvadis_phenodata.csv的示例文件包括在该协议的数据文件中(ChrX_data)。一般来说,你将不得不自己创建这个文件。它包含了关于你的RNA-seq样本的信息,格式如这个csv(逗号分隔的值)文件中所示。文件中的每一行包含每个样本的描述,每一列应该包含一个变量。为了将这个文件读入R,我们使用read.csv命令。在这个文件中,数值是由逗号分隔的,但如果文件是以制表符分隔的,你可以使用函数read.table。
1 | >pheno_data = read.csv("geuvadis_phenodata.csv") |
- 错误排除
- 读入由StringTie计算的表达数据。Ballgown目前也支持从Cufflinks6和RSEM17读取数据。要做到这一点,我们使用带有以下三个参数的ballgown命令:存储数据的目录(dataDir,这里简单命名为’Ballgown’),出现在样本名称中的模式(samplePattern)和我们在前一步加载的表型信息(pData)。请注意,一旦创建了Ballgown对象,任何其他Bioconductor32软件包都可以应用于数据分析或数据可视化。这里我们提出了一个标准化的流程,可以用来进行标准的差异表达分析。
1 | >bg_chrX = ballgown(dataDir = "ballgown", samplePattern = "ERR", pData=pheno_data) |
- 错误排除
- 使用筛选器去除低丰度基因。RNA-seq数据的一个常见问题是,某些基因经常有非常少的计数或零计数。一个常见的步骤是过滤掉其中的一些。另一种被用于基因表达分析的方法是应用一个方差过滤器。在这里,我们删除所有跨样本方差小于1的转录本。
1 | >bg_chrX_filt = subset(bg_chrX,"rowVars(texpr(bg_chrX)) >1",genomesubset=TRUE) |
- 识别在统计学上显示组间差异的转录本。我们可能想做的一件事是确保我们考虑到其他变量引起的表达变化。作为一个例子,我们将寻找在性别间有差异表达的转录本,同时校正由于群体变量造成的任何表达差异。我们可以使用Ballgown的stattest函数来做这件事。我们设置了getFC=TRUE参数,这样我们就可以查看两组之间经混杂因素调整后的折叠变化。
请注意,Ballgown的统计测试是一个基于标准线性模型的比较。对于小样本量(每组n<4),通常最好进行正则化。这可以用Bioconductor中的limma包来完成。其他正则化方法,如DESeq和edgeR可以应用于基因或外显子计数,但它们不适合直接应用于FPKM丰度估计。统计测试使用累积的上四分位数规范化(cumulative
upper quartile normalization)。
1 | >results_transcripts = stattest(bg_chrX_filt, |
- 识别组间有统计学差异的基因。为此,我们可以运行与识别差异表达的转录本相同的函数,在这里我们在stattest命令中设置feature=”gene”
1 | >results_genes = stattest(bg_chrX_filt, feature="gene", |
- 在results_transcripts数据框中添加基因名称和基因ID
1 | >results_transcripts = |
- 将结果从最小的P值排序到最大
1 | >results_transcripts = arrange(results_transcripts,pval) |
- 两性之间差异表达的转录物(q值<5%) 基因名|基因ID|转录本名|ID|Fold
Change|P值|q值 :–:|:–:|:–:|:–:|:–:|:–:|:–: XIST
MSTRG.506|NR_001564|1729|0.003255|7.0447e-10|1.6160e-06
|MSTRG.506|MSTRG.506.1|1728|0.016396|1.2501e-08|1.4339e-05
TSIX|MSTRG.505|NR_003255|1726|0.083758|2.4939e-06|1.9070e-03
|MSTRG.506|MSTRG.506.2|1727|0.047965|3.7175e-06|2.1319e-03
|MSTRG.585|MSTRG.585.1|1919|7.318925|9.3945e-06|3.7715e-03
PNPLA4|MSTRG.56|NM_004650|203|0.466647|9.8645e-06|3.7715e-03
|MSTRG.506|MSTRG.506.5|1731|0.046993|2.1350e-05|6.9968e-03
|MSTRG.592|MSTRG.592.1|1923|9.186257|3.5077e-05|1.0058e-02
|MSTRG.518|MSTRG.518.1|1744|11.972859|4.4476e-05|1.1336e-02
Fold
Change是指女性与男性的表达比例;因此,低于1的数值意味着转录物在男性中的表达水平较低。
- 将结果写入一个可以共享和分发的csv文件
1 | >write.csv(results_transcripts, "chrX_transcript_results.csv", |
- 识别q值<0.05的转录物和基因
1 | >subset(results_transcripts,results_transcripts$qval<0.05) |
每个命令的Ballgown输出都会出现在屏幕上;我们在下面的Table3(转录物)和Table4(基因)中显示了结果。如表所示,X号染色体有9个转录本在两性之间有差异表达(使用q值阈值为0.05),其中3个对应于已知基因的异构体(XIST、TSIX和PNPLA4)。在基因水平上(表4),X染色体在相同的q值临界点上有10个差异表达的基因。(如果使用不同版本的程序,这些表格可能包含略有不同的结果)。
错误排除
错误排除
错误排除