案例简介数据来自于发表在Nature commmunication 上的一篇文章 “Temporal dynamics of gene expression and histone marks at the Arabidopsis shoot meristem during flowerin”。原文用RNA-Seq的方式研究在开花阶段,芽分生组织在不同时期的基因表达变化。 原文的流程是： TopHat -> SummarizeOverlaps -> Deseq2 -> AmiGO 数据存放在http://www./arrayexpress/, ID为E-MTAB-5130。 实验设计： 4个时间段（0,1,2,3），分别有4个生物学重复，一共有16个样品。 数据下载首先下载数据说明文件： $ wget http://www./arrayexpress/files/E-MTAB-5130/E-MTAB-5130.sdrf.txt然后根据数据说明文件提供的FTP链接下载 $ head -n1 E-MTAB-5130.sdrf.txt | tr '\t' '\n' | nl | grep URI
33 Comment[FASTQ_URI]
$ tail -n 2 E-MTAB-5130.sdrf.txt | cut -f 33 | xargs -i wget {}根据下载速度，你们可以选择去吃吃饭，还是睡睡觉。 下载完RNA-Seq数据后，我们还需要下载一个拟南芥cDNA序列(记住是转录组，而不是全基因组，好了，你别说了，我记住了) $ curl ftp://ftp.ensemblgenomes.org/pub/plants/release-28/fasta/arabidopsis_thaliana/cdna/Arabidopsis_thaliana.TAIR10.28.cdna.all.fa.gz -o athal.fa.gz然后用Salmon建立索引： salmon index -t athal.fa.gz -i athal_index数据定量由于样本一共有16个，不可能一条条输入命令，所以我们写一个脚本： #! /bin/bash
for fn in ERR1698{194..209};
do
samp=`basename ${fn}`
echo 'Processin sample ${sampe}'
salmon quant -i athal_index -l A -1 ${samp}_1.fastq.gz -2 ${samp}_2.fastq.gz -p 8 -o quants/${samp}_quant
done根据你电脑的配置，你可以选择吃下午茶，还是选择睡个午觉。 数据导入R当然你完全不必真的去睡午觉，我们可以程序运行的时候准备一下
虽然 tximport(files, type = c('none', 'salmon', 'sailfish', 'kallisto', 'rsem'),
txIn = TRUE, txOut = FALSE, countsFromAbundance = c('no', 'scaledTPM',
'lengthScaledTPM'), tx2gene = NULL, varReduce = FALSE,
dropInfReps = FALSE, ignoreTxVersion = FALSE, geneIdCol, txIdCol,
abundanceCol, countsCol, lengthCol, importer = NULL)
dir <- 'C:/Users/Xu/Desktop/'
list.files(dir)
sample <- paste0('ERR1698',c(194,seq(202,209),seq(195,201)),'_quant')
files <- file.path(dir,'quants',sample,'quant.sf')
names(files) <- paste0('sample',seq(1,16))
all(file.exists(files))吐槽: 原本的我以为sample从1到16会和194到206是一一对应，但是万万没想到，他的对应关系居然是下图这个样子的。看来凡是都不能想当然，一定要仔细看他们的对应关系。 然后我们还要准备一个基因名和转录本名字相关的数据框 library(AnnotationHub)
ah <- AnnotationHub()
ath <- query(ah,'thaliana')
ath_tx <- ath[['AH52247']]
columns(ath_tx)
k <- keys(ath_tx,keytype = 'GENEID')
df <- select(ath_tx, keys=k, keytype = 'GENEID',columns = 'TXNAME')
tx2gene <- df[,2:1] # TXID在前， GENEID在后如果你电脑跑的够快，基本上这个时候就可以导入数据了。 # install.packages('readr')
# install.packages('rsjon')
library('tximport')
library('readr')
txi <- tximport(files, type = 'salmon', tx2gene = tx2gene)
# reading in files with read_tsv
# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
# summarizing abundance
# summarizing counts
# summarizing length
names(txi)
[1] 'abundance' 'counts' 'length'
[4] 'countsFromAbundance
head(txi$lenght)
head(txi$counts)由于后续要用DESeq2包做差异表达分析，所以需要用 library('DESeq2')
sampleTable <- data.frame(condition = factor(
rep(c('Day0','Day1','Day2','Day3'),each=4)))
rownames(sampleTable) <- colnames(txi$counts)
dds <- DESeqDataSetFromTximport(txi, sampleTable, ~condition)
dim(dds)数据预过滤： dds <- dds[ rwoSums(counts(dds)) > 1 ,]
dim(dds)差异表达分析标准的差异表达分析步骤在DESeq2只需要 # 超过100个样本用
# library('BiocParallel')
# register('MulticoreParam(4))
dds <- DESeq(dds)
res <- results(dds)
summary(res)results还有许多参数：
于是我分别用Day1,2,3与Day0进行比较： res0vs1 <- results(dds, contrast = c('condition','Day1','Day0'))
res0vs2 <- results(dds, contrast = c('condition','Day2','Day0'))
res0vs3 <- results(dds, contrast = c('condition','Day3','Day0'))最后找到了208个差异表达基因，105上调，103个下调；而文章是298个，148上调，156下调。 res1SigUp <- subset(res1, padj < 0.01 & log2FoldChange >0 )
res1SigDown <- subset(res1, padj < 0.01 & log2FoldChange <0 )
res2SigUp <- subset(res2, padj < 0.01 & log2FoldChange >0 )
res2SigDown <- subset(res2, padj < 0.01 & log2FoldChange <0 )
res3SigUp <- subset(res3, padj < 0.01 & log2FoldChange >0 )
res3SigDown <- subset(res3, padj < 0.01 & log2FoldChange <0 )
sigUpGene <- unique(c(rownames(res1SigUp),rownames(res2SigUp),rownames(res3SigUp)))
length(sigUpGene)
# 105
sigDownGene <- unique(c(rownames(res1SigDown),rownames(res2SigDown),rownames(res3SigDown)))
length(sigDownGene)
# 103探索性分析特殊基因的表达情况有些基因是分生组织的特异基因，如STM, KNOTTED-LIKE, FROM ARABIDOPSIS THALIANA 1(KNAT1), CLAVATA 1(CLV1)和CLV3, 从理论上应该有一定的表达量 早花同源异形基因（homeotic gene)，如AP1, APETALA 3(AP3), CAULIFLOWER(CAL)和发育中期和后期基因JAGGED(JAG), WUSCHEL RELATED HOMEOBO 1(WOX1) WOX3 从理论上应该不表达。 上面的基因名我们需要先转成TAIR的ID，才能从dds中提取相应的counts.然后使用 meristem_identify_genes <- select(th, keys = c('STM','KNAT1','CLV1','CLV3'),
columns = c('TAIR'), keytype = 'SYMBOL')
library(tidyverse)
msData <- dds[which(rownames(dds) %in% meristem_identify_genes$TAIR)]
msData1 <- as.data.frame(assay(msData))
msData1$gene = rownames(msData)
msData2 <- msData1 %>%
gather(sample,count, -gene) %>%
mutate(condition = rep(c('Day0','Day1','Day2','Day3'),each=20)) %>%
arrange(gene)由于这个画图过程，后续会用大很多次，为了避免重复操作，我写成了一个函数，我建议你根据上面的流程自己写一个函数，不要用我的。 plotGeneCounts <- function(genes) {
require(ggplot2)
require(tidyr)
require(dplyr)
GeneName <- select(th, keys=genes, columns=c('TAIR'),keytype = 'SYMBOL')
GeneName <- GeneName[order(GeneName$TAIR), ]
Data <- dds[which(rownames(dds) %in% GeneName$TAIR)]
Data1 <- as.data.frame(assay(Data))
Data1$gene <- GeneName$SYMBOL[GeneName$TAIR == rownames(Data1)]
Data2 <- Data1 %>%
gather(sample,count, -gene) %>%
mutate(condition = rep(c('Day0','Day1','Day2','Day3'), each= length(sample) / 4))
p <- ggplot(Data2)
p1 <- p geom_boxplot(aes(x=gene,y=count, fill=condition))
p2 <- p1 theme(axis.title.x = element_blank()) ylab('Counts')
print(p2)
return(GeneName)
}那么画图就是 原文用柱状图展示不同基因在不同处理下的表达量， 但是柱状图的信息太少，所以我们用ggplot2做箱线图。 同样的方法我们可以做出早花同源异形基因（homeotic gene)，如AP1, APETALA 3(AP3), CAULIFLOWER(CAL)和发育中期和后期基因JAGGED(JAG), WUSCHEL RELATED HOMEOBO 1(WOX1) WOX3 的箱线图： ## early floral homeotic genes
plotGeneCounts(c('AP1','AP3','CAL'))
# genes expressed in the centre or abaxial side of developing leaves21
plotGeneCounts(c('JAG','WOX1','WOX3'))和文章的结论基本一致，早花同源异形基因在我的流程中只能检测到CAL, 而发育中期和后期的基因基本上表达量也特别低。不过文章说CAL3这个基因能被稳定被检测到，实际上表达量与这些基因没有明显区别（< 100)，我也不知道他如何得到这个结论的。 文章在比较了不同时期的基因表达情况后，发现 首先我根据文章给出的基因，先作图看看是不是真的有这样的趋势： # morning-expressed genes of the central oscillator of the circadian clock
# evening-expressed clock gene2
mGene <- c('LHY', 'CCA1', 'ELF4', 'PRR5','LUX')
mGenes <-plotGeneCounts(mGene)看来很类似，下一步就是判断这些基因是不是也在我找到显著性差异表达的基因列表里。 > mGenes[mGenes$TAIR %in% sigUpGene,]
SYMBOL TAIR
1 LHY AT1G01060
2 CCA1 AT2G46830
> mGenes[mGenes$TAIR %in% sigDownGene,]
[1] SYMBOL TAIR
<0 行> (或0-长度的row.names)很尴尬，显著性上调的基因也都在我的列表中，但是显著性下调的，一个都没有，于是我把水平从0.01提高到0.05，重新检查下 res1SigDown <- subset(res1, padj < 0.05 & log2FoldChange <0 )
res2SigDown <- subset(res2, padj < 0.05 & log2FoldChange <0 )
res3SigDown <- subset(res3, padj < 0.05 & log2FoldChange <0 )
sigDownGene <- unique(c(rownames(res1SigDown),rownames(res2SigDown),rownames(res3SigDown)))
length(sigDownGene)
mGenes[mGenes$TAIR %in% sigDownGene,]
SYMBOL TAIR
4 PRR5 AT5G24470好吧还是只找到一个，毕竟ELF4这个表达量实在是太低了，除非水平变成是0.1，不然是难以发现的啦。但是这个还是证明作者的猜想是对的。 随后作者又去找了一些flowering prmotation gene SOC1, FD, NF-YA4, 以及flowering-time repressor, 如JMJ30 LDgenes<- c('SOC1','FD','NF-YA4','JMJ30')
plotGeneCounts(LDgenes)小结： 可以找一些特征基因用来验证RNA-Seq是否正确 GO富集分析文献通过GO富集分析发现：SAM诱导的基因大部分富集在long-day photoperiodism (GO:0048574) and regulation of circadian rhythm (GO:0042754） 为了验证这个结果是不是正确的，我使用Y叔良心之作 结果在上调中并没有找到富集基因，但是在下调（day1 vs day0 和 day3 vs day0)中却发现了好多富集基因，虽然场面十分尴尬，但是至少说明下调的基因大多属于ribosome res1DownGO <- enrichGO(gene = rownames(res1SigDown),
OrgDb = th,
keytype = 'TAIR',
pvalueCutoff = 0.1,
pAdjustMethod = 'fdr'
)res3DownG0 <- enrichGO(gene = rownames(res3SigDown),
OrgDb = th,
keytype = 'TAIR',
pvalueCutoff = 0.1,
)
dotplot(res3DownG0)
或许我可能需要用文章说的AmiGO 2 versions 2.3.2 and 2.4.24进行GO富集分析，或许用DAVID，其他工具分析一遍，感觉又要许多主观能动性了。 总结在这次实战中，我学习了用tidyr, dplyr对数据进行预处理，重新看了一下ggplot2的图形语法，复习了之前差异表达分析的基本操作。 |
|
|
