6. 基因差异表达分析
利用DESeq2软件基于样本的原始reads数目计算差异表达基因，利用run-featurecounts.R脚本对每个样本的reads数进行定量，然后利用abundance_estimates_to_matrix.pl脚本合并所有的定量结果，最后利用DESeq2进行差异表达基因的分析。
首先，将所需的脚本文件全部拷贝到工作目录下，同时将样本信息相关文件也拷贝保存。
cp -r xxx/script_featurecounts . cp -r xxx/script-mergecounts . cp xxx/RNA-Seq_Practice/contrast.txt . cp xxx/RNA-Seq_Practice/sampleinfo.txt .
新建一个文件夹，运行R语言脚本，对样本进行表达定量分析，以PR1为例，输入文件为02文件夹中生成的bam文件和基因组注释文件，生成输出文件PR1，其他的几个样本按照同样的方法进行处理。
mkdir R-featurecounts Rscript script_featurecounts/run-featurecounts.R -b ../02_hisat2Mapping/PR1.srt.bam #输入文件 -g ../Ref/Ref.gtf.gz #基因组注释文件 -o R-featurecounts/PR1 #输出文件
运行结束后，在输出文件夹内将会生成count文件和log文件，利用perl语言脚本abundance_estimates_to_matrix.pl合并所有样本的定量结果。
perl script-mergecounts/abundance_estimates_to_matrix.pl --est_method featureCounts R-featurecounts/*.count
利用trinity 软件中的DESeq2分析差异表达基因，输入文件有counts矩阵、所选方法DESeq2、样本信息表sampleinfo.txt、分组信息contrasts.txt。
perl xxx/run_DE_analysis.pl --matrix matrix.counts.matrix --method DESeq2 --samples_file sampleinfo.txt #样本信息表文件 --contrasts contrasts.txt #这里需要注意文件名称是否正确
上述流程输出的结果存放在随机文件夹中，进入该文件夹后，可以看到生成的差异表达基因结果matrix.counts.matrix.PR_vs_SR.DESeq2.DE_results文件，然后利用awk命令提取出满足条件（同时满足：|log2FoldChange| > 1，且padj < 0.05）的差异表达基因,生成新的输出文件件PRvsSR_DEGs。
drwxr-xr-x 307 Nov 29 09:54 DESeq2.197264.dir #随机新目录 $ ls -f #查看新目录下文件列表 matrix.counts.matrix.PR_vs_SR.DESeq2.Rscript matrix.counts.matrix.PR_vs_SR.DESeq2.DE_results #差异表达基因结果文件 matrix.counts.matrix.PR_vs_SR.DESeq2.count_matrix matrix.counts.matrix.PR_vs_SR.DESeq2.DE_results.MA_n_Volcano.pdf awk 'sqrt($7*$7)>1 && $11 <0.05 {print $0}' ./matrix.counts.matrix.PR_vs_SR.DESeq2.DE_results > PRvsSR.DESeq2_DEGs #利用awk命令提取满足条件的指定数据
将最终得出的差异表达基因数据文件PRvsSR.DESeq2_DEGs保存到电脑，进行GO富集和KEGG富集分析。绘制火山图时使用的是matrix.counts.matrix.PR_vs_SR.DESeq2.DE_results
7. GO富集分析
利用在线网站(http://shiny.host.ugeneyun.cn/GOmaize/)进行GO富集分析，输入文件为差异表达基因的ID信息表，此处不进行显著性过滤，原因是前期已经进行筛选了。提交运行后生成GOenrich结果，共有270条Term，下图随机展示5条，可以得出GO编号、功能描述、显著性和分类信息。然后，利用网站的绘图功能制作GO富集气泡图，保存到本地GO.pdf。
8. KEGG分析
在R语言（4.2.2）中使用clusterProfiler、KEGG.db、ggplot2三个R包进行KEGG富集分析，先构建一个OrgDb本地索引库，然后利用enrichKEGG函数进行富集分析，利用在线网站(https://www.maizegdb.org/gene_center/gene)将基因ID转换为GenBank Gene类型，保存为文件gene.txt，输出文件为pathway数据表gene.pathway.csv 和KEGG富集结果gene.pathway .pdf
remotes::install_github("YuLab-SMU/createKEGGdb",force = TRUE) createKEGGdb::create_kegg_db('zma') #首先使用该软件制作玉米KEGG富集包 library(clusterProfiler) library(ggplot2) # 绘图需用R包 library(KEGG.db) # 该步骤需要手动在Rstudio中安装KEGG.db包 file="gene.txt" #输入文件，包含18条差异表达的基因信息 gene <- read.delim("gene.txt",header = F,row.names = 1) > head(gene,2) GenBank.Gene Zm00001d0218 1036647 Zm00001d0234 1002408 kk <- enrichKEGG(gene= gene$GenBank.Gene,organism= 'zma', pvalueCutoff = 1, qvalueCutoff = 1, use_internal_data =T) #进行KEGG富集分析 write.csv(kk,file = paste0("gene.pathway.csv"),row.names = F) p=dotplot(kk);p #绘制点图，并输出plot查看 ggsave(p,filename = "gene.pathway.pdf"),width = 8,height = 7) #保存图片
9. 差异表达基因火山图
通过ggplot2软件绘制基因差异表达火山图，使用DESeq2软件生成的matrix.counts.matrix.PR_vs_SR.DESeq2.DE_results为输入文件，该文件内geneid、sampleA、sampleB、baseMeanA、baseMeanB、baseMean、log2FoldChange、lfcSE、stat、pvalue、padj等信息，将其导入R语言中，提取1,6,7,8,9,10,11列，然后设置数据的行列名称
library(ggplot2) # 载入R包 res <- read.table("matrix.counts.matrix.PR_vs_SR.DESeq2.DE_results") res <- res[,c(1,6,7,8,9,10,11)] #筛选数据，提取指定列 rownames(res) <- res$V1 #命名数据框的行名为基因ID colnames(res) <- c("geneid","baseMean", "log2FoldChange","lfcSE","stat","pvalue","padj") # 修改列名 deseq_res <- as.data.frame(res[order(res$padj), ]) # 按照p值大小排序 write.table(deseq_res, 'DESeq2-test.txt', row.names = FALSE, sep = '\t', quote = FALSE) #保存差异表达基因数据
设置差异表达分类sig，最后利用ggplot函数绘制火山图，保存为'volcano.svg'
deseq_res <- read.delim('DESeq2-test.txt', sep = '\t') deseq_res[which(deseq_res$log2FoldChange >= 1 & deseq_res$padj < 0.05),'sig'] <- 'up (p.adj < 0.05, log2FC >= 1)' #上调基因 deseq_res[which(deseq_res$log2FoldChange <= -1 & deseq_res$padj < 0.05),'sig'] <- 'down (p.adj < 0.05, log2FC <= -1)' #下调基因 deseq_res[which(abs(deseq_res$log2FoldChange) < 1 | deseq_res$padj >= 0.05),'sig'] <- 'no diff' #无显著差异 p <- ggplot(deseq_res, aes(log2FoldChange, -log(padj, 10))) + geom_point(aes(color = sig), alpha = 0.6, size = 1) + scale_color_manual(values = c('blue2', 'gray30', 'red2')) + theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black') legend.position = c(0.26, 0.92)) + theme(legend.title = element_blank(), legend.key = element_rect(fill = 'transparent'), legend.background = element_rect(fill = 'transparent')) + geom_vline(xintercept = c(-1, 1), color = 'gray', size = 0.25) + geom_hline(yintercept = -log(0.05, 10), color = 'gray') + labs(x = 'log2 Fold Change', y = '-log10 p-value', color = NA) + xlim(-5, 5) ggsave('volcano.svg', p, width = 4, height = 3) # 保存火山图
参考资料： [1] Brown J, Pirrung M, McCue L A. FQC Dashboard: integrates FastQC results into a web-based, interactive, and extensible FASTQ quality control tool[J]. Bioinformatics, 2017, 33(19): 3137-3139. [2] Bolger A M, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data[J]. Bioinformatics, 2014, 30(15): 2114-2120. [3] Kim D, Paggi J M, Park C, et al. Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype[J]. Nature biotechnology, 2019, 37(8): 907-915. [4] Liao Y, Smyth G K, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features[J]. Bioinformatics, 2014, 30(7): 923-930. [5] Love M I, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2[J]. Genome biology, 2014, 15(12): 1-21. [6] Yu G, Wang L G, Han Y, et al. clusterProfiler: an R package for comparing biological themes among gene clusters[J]. Omics: a journal of integrative biology, 2012, 16(5): 284-287. [7] Wickham H. Data analysis[M]//ggplot2. Springer, Cham, 2016: 189-201.
END
© 素材来源于网络，侵权请联系后台删除
往期推荐：
