RNA-seq分析文章结果重现

文章：Temporal dynamics of gene expression and histone marks at the Arabidopsis shoot meristem during flowering.

一、查看并下载数据

tr 替换\t 为\n ； nl为对每行加行号
head -n1 E-MTAB-5130.sdrf.txt | tr '\t' '\n' | nl | grep URI

image


tail -n +2 E-MTAB-5130.sdrf.txt | cut -f 33 | xargs -i wget {}

二、fastqc和MultiQC质控

/usr/local/bin/fastqc -o ./ --nogroup /sas/supercloud-kong/wangtianpeng/data/ath_flowering/1_raw_data/ERR1698201_2.fastq.gz 

### multiqc质控
multiqc *fastqc.zip --pdf

三、Trimmo修建

用Trimmomatic 从5‘端开始切除低质量碱基（LEADING：5）TRAILING:20表示切除碱基质量小于20的碱基，MINLEN:50表示过滤切除后长度小于50个碱基的reads
滑动窗口：SLIDINGWINDOW:3：20表示设置窗口大小为3个碱基，从5'端开始切除reads，read的平均质量低于20，开始切除后边的碱基序列。过滤长度小于50碱基的reads，使用更大的窗口可以减轻切除的力度

for i in `seq 194 209`; do  i=ERR1698$i; nohup trimmomatic PE ../1_raw_data/$i\_1.fastq.gz ../1_raw_data/$i\_2.fastq.gz $i\_P1.fastq.gz $i\_U1.fastq.gz $i\_P2.fastq.gz $i\_U2.fastq.gz ILLUMINACLIP:/home/wangtianpeng/anaconda3/share/trimmomatic-0.36-5/adapters/TruSeq3-PE-2.fa:2:30:10 LEADING:5 TRAILING:5 SLIDINGWINDOW:4:5 MINLEN:20 \&; done

四、SALMON定量

salmon index -t athal.fa.gz -i athal_index

for i in `seq 194 209`;do nohup salmon quant -i /sas/supercloud-kong/wangtianpeng/Database/plant/a_th/ath_index -l A -1 /sas/supercloud-kong/wangtianpeng/data/ath_flowering/3_trimm/ERR1698${i}_P1.fastq.gz -2 /sas/supercloud-kong/wangtianpeng/data/ath_flowering/3_trimm/ERR1698${i}_P2.fastq.gz -p 25 -o ./ERR1698${i}_quant & done

五、tximport的输入文件：涉及R的一些基本操作

将转录本水平的定量转为基因水平的定量：
准备files


dir<- getwd()
## "E:/Bio_informatics_software/转录组/NC文章重复_A.th"

list.files()
sample <- paste0("ERR1698", c(194, seq(202,209),seq(195,201)), "_quant")
### "ERR1698194_quant" "ERR1698202_quant" "ERR1698203_quant" "ERR1698204_quant"

files <- file.path(dir,sample,"quant.sf")
names(files) <- paste0("sample",1:16)
###                                                                               sample1 
"E:/Bio_informatics_software/转录组/NC文章重复_A.th/ERR1698194_quant/quant.sf" 
        sample2 
"E:/Bio_informatics_software/转录组/NC文章重复_A.th/ERR1698202_quant/quant.sf" 
        sample3 
"E:/Bio_informatics_software/转录组/NC文章重复_A.th/ERR1698203_quant/quant.sf"


file.exists(files)
all(file.exists(files))

准备基因名和转录本名字的数据框

library(AnnotationHub)
ah <- AnnotationHub()
### AnnotationHub with 40126 records
# snapshotDate(): 2017-04-25 
# $dataprovider: BroadInstitute, Ensembl
# $species:
# $rdataclass: GRanges, BigWigFile ...


ath <- query(ah, 'thaliana')
# snapshotDate(): 2017-04-25 
# $dataprovider:
# $species: 
# $rdataclass:
# additional mcols():

ath_tx <- ath[[ 'AH52247' ]]

k <- keys(ath_tx, keytype= "GENEID")
df <- select(ath_tx, keys= k, keytype= "GENEID", columns= "TXNAME")
tx2name <- df[,2:1]

利用tximport转换
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)

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)
### txi 为一个较大的list ，查看有哪些变量。"abundance"           "counts"              "length"  "countsFromAbundance
head(txi$counts)

六、DESeq2差异表达分析

回顾：对于DESeq2需要输入的三个数据：表达矩阵、样品信息数据框、差异比较矩阵
而对于DESeq2的差异分析步骤，总结起来就是三步：
- 构建一个dds(DESeqDataSet)的对象；
- 利用DESeq函数进行标准化处理；
- 用result函数来提取差异比较的结果。
构建dds矩阵的基本代码

dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design= ~ batch + condition)

其输入的三个文件：
- 表达矩阵：即代码中的countData，就是通过之前的HTSeq-count生成的reads-count计算融合的矩阵。行为基因名，列为样品名，值为reads或fragment的整数。
- 样品信息矩阵：即上述代码中的colData，它的类型是一个dataframe（数据框），第一列是样品名称，第二列是样品的处理情况（对照还是处理等）。可以从表达矩阵中导出或是自己单独建立。
- 差异比较矩阵：即代码中的design，差异比较矩阵就是告诉差异分析函数哪些是对照，哪些是处理。
正式构建dds的矩阵：

library(DESeq2)

### DESeqDataSetFromTximport(txi, colData, design)

condition<- factor(rep(c("DAY0", "DAY1", "DAY2", "DAY3"),each=4))
sampleTable <- data.frame(row.names = colnames(txi$counts),condition)
dds <- DESeqDataSetFromTximport(txi, sampleTable, ~condition)

[sampleTable ]

image

数据的预过滤：将所有样本基因表达量之和小于1的基因过滤掉


dim(dds) ### 显示32753个 16个sample
dds <- dds[rowSums(counts(dds))>1, ]
dim(dds)  ### 显示 24937个 16个sample

dds_out <- DESeq(dds)
res <- results(dds_out)
summary(res)

RESULTS 的参数：
- contrast :用于指定比较的对象，DAY1，DAY0
- lfcThreshold：用于指定log2 fold change的阈值
- alpha ：显著性水平的阈值
- pAdjustMethod ：多重实验校正
分别用 DAY1，2，3和DAY0做比较：

res0_1 <- results(dds_out, contrast = c("condition","DAY1","DAY0"))
res0_2 <- results(dds_out,contrast =c("condition","DAY2","DAY0"))
res0_3<- results(dds_out,contrast =c("condition","DAY3","DAY0"))

res0_1UP <- subset(res0_1,padj< 0.01 & log2FoldChange >0 )
res0_2UP <- subset(res0_2,padj< 0.01 & log2FoldChange >0 )
res0_3UP <- subset(res0_3,padj< 0.01 & log2FoldChange >0 )
res0_1DOWN <- subset(res0_1,padj< 0.01 & log2FoldChange <0 )
res0_2DOWN <- subset(res0_2,padj< 0.01 & log2FoldChange <0 )
res0_3DOWN <- subset(res0_3,padj< 0.01 & log2FoldChange <0 )

UP_GENE <- unique(rownames(res0_1UP),rownames(res0_2UP),rownames(res0_3UP))
DOWN_GENE <- unique(rownames(res0_1DOWN),rownames(res0_2DOWN),rownames(res0_3DOWN))
length(UP_GENE) ### 97个
length(DOWN_GENE)  ## 95个

七、特殊基因的表达情况

下载 org.db。参照【转录组入门八-功能注释】

library(AnnotationHub)
hub <- AnnotationHub()
query（hub, "thaliana"）
ath.tair.db <- ath[["AH53758"]]
keytypes(ath.tair.db)##查看有哪些数据类型，包含着各大主流数据库的数据。
###用select函数，就可以把任意公共数据库的数据进行一一对应。
### keys是原始的ID，columns是转换之后的ID，keytype是要指定的原始ID类型

转换基因名+直接画图

plotGeneCounts <- function(genes) {
  require(ggplot2)
  require(tidyr)
  require(dplyr)
  GeneName <- select(ath.tair.db, 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)
}

结果：plotGeneCounts(c("STM","KNAT1","CLV1", "CLV3"))

image
结果：后期表达的基因表达量十分低。
plotGeneCounts(c("JAG","WOX1","WOX3"))

image

-结果：早期花的同源异型基因的表达量十分低。

plotGeneCounts(c("AP1","AP3","CAL"))

image
长日照基因白天表达的基因 LHY ，CCA1 在植物转移到长日照条件下表达量显著提高，而夜晚表达的基因PRR5表达量显著下降。并且复合体编码基因LUX, ELF4也类似情况。
plotGeneCounts(c("CCA1","LHY","PRR5","LUX","ELF4"))

image
成花转变的一些关键性基因：如FD, SOC1，这些基因表达量逐渐上升。

image

RNA-seq分析文章结果重现

一、查看并下载数据

二、fastqc和MultiQC质控

三、Trimmo修建

四、SALMON定量

五、tximport的输入文件：涉及R的一些基本操作

六、DESeq2差异表达分析

七、 特殊基因的表达情况

推荐阅读更多精彩内容

七、特殊基因的表达情况