一、查看并下载数据
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的差异分析步骤，总结起来就是三步：
- <font color=darkred>构建一个dds(DESeqDataSet)的对象</font>；
- <font color=darkred>利用DESeq函数进行标准化处理</font>；
- <font color=darkred>用result函数来提取差异比较的结果</font>。
-
构建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
网友评论