四川网站建设报价,wordpress mu 插件,抖音推广方式有哪些,企业网站网站建设公司R语言在RNA-seq中的应用 文章目录 R语言在RNA-seq中的应用生成工作流环境读取和处理数据由targets文件提供实验定义对实验数据进行质量过滤和修剪生成FASTQ质量报告 比对建立HISAT2索引并比对 读长量化读段计数样本间的相关性分析 差异表达分析运行edgeR可视化差异表达结果计算…R语言在RNA-seq中的应用 文章目录 R语言在RNA-seq中的应用生成工作流环境读取和处理数据由targets文件提供实验定义对实验数据进行质量过滤和修剪生成FASTQ质量报告 比对建立HISAT2索引并比对 读长量化读段计数样本间的相关性分析 差异表达分析运行edgeR可视化差异表达结果计算和绘制差异表达基因DEG集合的Venn图 GO富集分析准备工作批量进行GO富集分析绘制批量GO富集分析的结果 层次聚类和热图绘制 生成工作流环境
之后代码运行可能会有网络问题通过set_config函数设置即可。
.libPaths(D:/000大三下/R语言/实验/Lab4/lab4packages)
library(httr)
set_config(use_proxy(url127.0.0.1, port7890)
)
library(systemPipeR)
library(systemPipeRdata)
#####################################################
## 1. Generate workflow environment
#####################################################
setwd(choose.dir())
genWorkenvir(workflow rnaseq)
setwd(rnaseq)读取和处理数据
由targets文件提供实验定义
读取和预处理实验数据。具体步骤如下
首先使用system.file函数找到targets.txt文件的路径这个文件包含了所有的FASTQ文件和样本比较的信息。然后使用read.delim函数读取targets.txt文件忽略以#开头的注释行并只保留前四列。最后打印出targets对象查看数据的结构和内容。
## 2. Read preprocessing
#####################################################
## 2.1 Experiment definition provided by targets file
## The targets file defines all FASTQ files and sample comparisons
## of the analysis workflow.
targetspath - system.file(extdata, targets.txt, package systemPipeR)
targets - read.delim(targetspath, comment.char #)[, 1:4]
targets对实验数据进行质量过滤和修剪
对实验数据进行质量过滤和修剪。具体步骤如下
首先使用loadWorkflow函数从cwl和yml参数文件以及targets文件中构建一个SYSargs2对象这个对象包含了执行trimLRPatterns函数的所有参数和输入输出路径。然后使用renderWF函数根据targets文件中的文件名和样本名替换cwl和yml文件中的占位符生成一个完整的工作流对象。接着打印出trim对象查看工作流的各个组成部分。最后使用output函数查看输出路径中的前两个修剪后的FASTQ文件。
## 2.2 Read quality filtering and trimming
## The function preprocessReads allows to apply predefined or custom read
## preprocessing functions to all FASTQ files referenced in a SYSargs2 container, such
## as quality filtering or adapter trimming routines. The paths to the resulting output
## FASTQ files are stored in the output slot of the SYSargs2 object. The following
## example performs adapter trimming with the trimLRPatterns function from the
## Biostrings package. After the trimming step a new targets file is generated (here
## targets_trim.txt) containing the paths to the trimmed FASTQ files. The new
## targets file can be used for the next workflow step with an updated SYSargs2
## instance, e.g. running the NGS alignments using the trimmed FASTQ files.
## First,we construct SYSargs2 object from cwl and yml param and targets files.
dir_path - system.file(extdata/cwl/preprocessReads/trim-se, package systemPipeR)
trim - loadWorkflow(targets targetspath, wf_file trim-se.cwl, input_file trim-se.yml, dir_path dir_path)
trim - renderWF(trim, inputvars c(FileName _FASTQ_PATH1_, SampleName _SampleName_))
trim
output(trim)[1:2]生成FASTQ质量报告
生成FASTQ质量报告。具体步骤如下
首先使用seeFastq函数对trim对象中的输入文件进行质量分析计算每个文件的碱基质量分布、序列长度分布、GC含量分布和k-mer频率分布。然后使用pdf函数创建一个PDF文件用于保存质量报告的图形。接着使用seeFastqPlot函数绘制质量报告的图形包括每个文件的四个子图。最后使用dev.off函数关闭PDF设备完成图形的保存。
## 2.3 FASTQ quality report
fqlist - seeFastq(fastq infile1(trim), batchsize 10000, klength 8)
pdf(./results/fastqReport.pdf, height 18, width 4 * length(fqlist))
seeFastqPlot(fqlist)
dev.off()[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-8o8Xh8q9-1685677405310)(D:/000%E5%A4%A7%E4%B8%89%E4%B8%8B/R%E8%AF%AD%E8%A8%80/%E5%AE%9E%E9%AA%8C/Lab4/rnaseq/results/fastqReport.png)]
比对
建立HISAT2索引并比对
使用HISAT2进行短读比对的步骤。主要包括以下几个步骤
构建HISAT2索引首先代码中使用loadWorkflow函数加载了一个工作流程对象idx并指定了索引构建的相关参数。然后通过调用renderWF函数渲染工作流程对象并通过cmdlist函数获取相应的命令列表。最后使用runCommandline函数运行命令列表来构建HISAT2索引。进行映射接下来代码中使用loadWorkflow函数加载了另一个工作流程对象args并指定了映射的相关参数。通过调用renderWF函数渲染工作流程对象将文件路径和样本名称作为输入变量进行替换。然后通过调用cmdlist函数获取命令列表以及通过output函数获取输出结果的相关信息。运行命令行模式在代码中使用runCommandline函数来运行命令列表以进行短读比对。处理输出文件在代码中使用output_update函数来修改args对象以模拟对生成的对齐结果文件进行处理。具体操作是将输出文件的后缀名修改为.sam和.bam并将文件路径中的目录设置为FALSE以方便后续处理。检查生成的BAM文件最后代码中使用subsetWF函数选择输出结果中的BAM文件路径并通过file.exists函数检查这些文件是否存在。
## 3.1 Read mapping with HISAT2
## The following steps will demonstrate how to use the short read aligner Hisat2 (Kim,
## Langmead, and Salzberg 2015) in both interactive job submissions and batch
## submissions to queuing systems of clusters using the systemPipeRs new CWL
## command-line interface.
## First, build HISAT2 index. (Skip this step)
#dir_path - system.file(extdata/cwl/hisat2/hisat2-idx, package systemPipeR)
#idx - loadWorkflow(targets NULL, wf_file hisat2-index.cwl,
# input_file hisat2-index.yml, dir_path dir_path)
#idx - renderWF(idx)
#idx
#cmdlist(idx)
#runCommandline(idx, make_bam FALSE)## Second, mapping.
dir_path - system.file(extdata/cwl/hisat2, package systemPipeR)
args - loadWorkflow(targets targetspath, wf_file hisat2-mapping-se.cwl, input_file hisat2-mapping-se.yml, dir_path dir_path)
args - renderWF(args, inputvars c(FileName _FASTQ_PATH1_, SampleName _SampleName_))
args
cmdlist(args)[1:2]
output(args)[1:2]
## Runc single Machine
#args - runCommandline(args) (Skip this)
# Move all *bam files from Lab_4 files/Bam to rnaseq/results
# This command is used to modify class args to simulate alignment. (Weihan)
args - output_update(args, dir FALSE, replace TRUE, extension c(.sam, .bam))
## Check whether all BAM files have been created.
outpaths - subsetWF(args, slot output, subset 1, index 1)
file.exists(outpaths)读长量化
读段计数
使用多核心并行模式下的summarizeOverlaps函数进行读段计数的过程。主要包括以下几个步骤 导入所需的包代码中使用library函数导入了GenomicFeatures和BiocParallel包用于进行基因组特征和并行计算的操作。 创建转录组数据库通过makeTxDbFromGFF函数根据GFF文件创建了一个转录组数据库对象txdb用于存储基因组注释信息。 读取比对结果通过readGAlignments函数读取比对结果文件BAM文件并存储到变量align中展示了如何将BAM文件读取到R中进行后续处理。 定义感兴趣的外显子基因区域通过exonsBy函数根据转录组数据库和给定的外显子区域生成了一个外显子基因区域对象eByg。 创建BAM文件列表通过BamFileList函数创建了一个BAM文件列表对象bfl用于存储需要进行读取计数的BAM文件路径。 并行计算设置通过MulticoreParam函数创建了一个多核心并行计算参数对象multicoreParam并通过register函数注册该参数。 执行读段计数通过bplapply函数以并行计算的方式对BAM文件列表中的每个文件执行summarizeOverlaps函数进行读取计数。计数结果存储在counteByg列表中。 处理计数结果将计数结果转化为数据框形式countDFeByg并设置行名和列名。 计算RPKM通过apply函数以每列为单位调用returnRPKM函数计算RPKM值rpkmDFeByg。 输出结果使用write.table函数将计数结果和RPKM结果分别写入到results/countDFeByg.xls和results/rpkmDFeByg.xls文件中。 数据示例展示使用read.delim函数读取计数结果和RPKM结果文件的部分数据进行展示。 注意说明了在大多数统计差异表达或丰度分析方法如edgeR或DESeq2中应使用原始计数值作为输入。RPKM值的使用应限制在一些特殊应用中例如手动比较不同基因或特征的表达水平。
## 4.1 Read counting with summarizeOverlaps in parallel mode using multiple cores
## Reads overlapping with annotation ranges of interest are counted for each sample
## using the summarizeOverlaps function (Lawrence et al. 2013). The read counting is
## preformed for exonic gene regions in a non-strand-specific manner while ignoring
## overlaps among different genes. Subsequently, the expression count values are
## normalized by reads per kp per million mapped reads (RPKM). The raw read count
## table (countDFeByg.xls) and the corresponding RPKM table (rpkmDFeByg.xls) are
## written to separate files in the directory of this project. Parallelization is achieved
## with the BiocParallel package, here using 8 CPU cores.
library(GenomicFeatures)
library(BiocParallel)
txdb - makeTxDbFromGFF(file data/tair10.gff, format gff, dataSource TAIR, organism Arabidopsis thaliana)
saveDb(txdb, file ./data/tair10.sqlite)
txdb - loadDb(./data/tair10.sqlite)
outpaths - subsetWF(args, slot output, subset 1, index 1)
(align - readGAlignments(outpaths[1])) # 报错 Demonstrates how to read bam file into R
eByg - exonsBy(txdb, by c(gene))
bfl - BamFileList(outpaths, yieldSize 50000, index character())
multicoreParam - MulticoreParam(workers 2) # Not supported on Windows, dont worry.
register(multicoreParam)
registered()
counteByg - bplapply(bfl, function(x) summarizeOverlaps(eByg, x, mode Union, ignore.strand TRUE, inter.feature FALSE, singleEnd TRUE))
countDFeByg - sapply(seq(along counteByg), function(x) assays(counteByg[[x]])$counts)
rownames(countDFeByg) - names(rowRanges(counteByg[[1]]))
colnames(countDFeByg) - names(bfl)
rpkmDFeByg - apply(countDFeByg, 2, function(x) returnRPKM(counts x, ranges eByg))
write.table(countDFeByg, results/countDFeByg.xls, col.names NA, quote FALSE, sep \t)
write.table(rpkmDFeByg, results/rpkmDFeByg.xls, col.names NA, quote FALSE, sep \t)## Sample of data slice of count table
read.delim(results/countDFeByg.xls, row.names 1, check.names FALSE)[1:4,1:5]
## Sample of data slice of RPKM table
read.delim(results/rpkmDFeByg.xls, row.names 1, check.names FALSE)[1:4,1:4]
## Note, for most statistical differential expression or abundance analysis methods,
## such as edgeR or DESeq2, the raw count values should be used as input. The usage
## of RPKM values should be restricted to specialty applications required by some
## users, e.g. manually comparing the expression levels among different genes or
## features.样本间的相关性分析
进行样本间的相关性分析。利用DESeq2包进行样本间相关性分析的过程。首先将计数数据导入然后构建DESeq2数据集并计算样本间的Spearman相关系数。最后通过层次聚类和绘图将相关性结果以聚类图的形式保存在PDF文件中。主要包括以下几个步骤 导入所需的包通过library函数导入了DESeq2和ape包用于进行差异表达分析和绘制聚类图的操作。 读取计数数据通过read.table函数将计数结果文件results/countDFeByg.xls读取为一个矩阵countDF。 构建DESeq2数据集通过DESeqDataSetFromMatrix函数将计数数据countDF和条件数据colData构建为一个DESeq2数据集dds指定了条件设计。 计算相关系数通过cor函数计算基于rlog转换后的表达值的Spearman相关系数。将rlog函数应用于DESeq2数据集dds的表达值然后通过assay函数提取表达矩阵并计算相关系数d。 层次聚类通过hclust函数对距离矩阵dist(1 - d)进行层次聚类生成聚类树对象hc。 绘制聚类图通过pdf函数创建一个PDF文件“results/sample_tree.pdf”然后使用plot.phylo函数绘制聚类树as.phylo(hc)并设置图形的样式参数。 保存聚类图通过dev.off函数关闭PDF文件保存并生成聚类图文件。
## 4.2 Sample-wise correlation analysis
## The following computes the sample-wise Spearman correlation coefficients from the
## rlog transformed expression values generated with the DESeq2 package. After
## transformation to a distance matrix, hierarchical clustering is performed with the
## hclust function and the result is plotted as a dendrogram (also see file
## sample_tree.pdf).
library(DESeq2, quietly TRUE)
library(ape, warn.conflicts FALSE)
countDF - as.matrix(read.table(./results/countDFeByg.xls))
colData - data.frame(row.names targets.as.df(targets(args))$SampleName, condition targets.as.df(targets(args))$Factor)
dds - DESeqDataSetFromMatrix(countData countDF, colData colData, design ~condition)
d - cor(assay(rlog(dds)), method spearman)
hc - hclust(dist(1 - d))
pdf(results/sample_tree.pdf)
plot.phylo(as.phylo(hc), type p, edge.col blue, edge.width 2, show.node.label TRUE, no.margin TRUE)
dev.off()差异表达分析
运行edgeR
使用edgeR包进行差异表达分析的过程。通过读取计数数据和目标样本信息定义比较组运行edgeR分析并将结果输出到文件中。另外还使用biomaRt包获取基因描述信息并将其添加到差异表达结果中。主要包括以下几个步骤
导入所需的包通过library函数导入了edgeR和biomaRt包用于差异表达分析和获取基因描述信息的操作。读取计数数据和目标样本信息通过read.delim函数分别读取计数数据文件results/countDFeByg.xls和目标样本信息文件targets.txt。定义比较组通过readComp函数从目标样本信息中提取比较组信息将其存储在变量cmp中。运行edgeR通过run_edgeR函数利用glm方法对计数数据进行差异表达分析。传入计数数据countDF、目标样本信息targets和比较组信息cmp[[1]]并指定independent FALSE表示非独立比较。添加基因描述信息通过使用biomaRt包中的useMart和getBM函数连接到植物数据库plants_mart并获取基因描述信息。将基因描述信息添加到差异表达结果edgeDF的数据框中。输出结果使用write.table函数将差异表达结果edgeDF写入到./results/edgeRglm_allcomp.xls文件中以制表符分隔不加引号列名不写入文件。
## 5. Analysis of DEGs
#####################################################
## The analysis of differentially expressed genes (DEGs) is performed with the glm
## method of the edgeR package (Robinson, McCarthy, and Smyth 2010). The sample
## comparisons used by this analysis are defined in the header lines of the
## targets.txt file starting with CMP.## 5.1 Run edgeR
library(edgeR)
countDF - read.delim(results/countDFeByg.xls, row.names 1, check.names FALSE)
targets - read.delim(targets.txt, comment #)
cmp - readComp(file targets.txt, format matrix, delim -)
edgeDF - run_edgeR(countDF countDF, targets targets, cmp cmp[[1]], independent FALSE, mdsplot )
## Add gene descriptions
library(biomaRt)
m - useMart(plants_mart, dataset athaliana_eg_gene, host https://plants.ensembl.org)
desc - getBM(attributes c(tair_locus, description), mart m)
desc - desc[!duplicated(desc[, 1]), ]
descv - as.character(desc[, 2])
names(descv) - as.character(desc[, 1])
edgeDF - data.frame(edgeDF, Desc descv[rownames(edgeDF)], check.names FALSE)
write.table(edgeDF, ./results/edgeRglm_allcomp.xls, quote FALSE, sep \t, col.names NA)可视化差异表达结果
筛选和可视化差异表达结果首先读取差异表达结果文件然后根据设定的筛选条件进行筛选并将筛选结果绘制成图形保存在PDF文件中。此外还将筛选后的DEG统计结果输出到文件中。主要包括以下几个步骤 读取差异表达结果通过read.delim函数读取差异表达结果文件results/edgeRglm_allcomp.xls。 绘制DEG结果图通过pdf函数创建一个PDF文件“results/DEGcounts.pdf”然后使用filterDEGs函数对差异表达结果进行筛选和绘图。筛选条件通过filter参数指定其中包括折叠变化Fold和调整的p值FDR。绘图结果保存在PDF文件中。 输出DEG统计结果使用write.table函数将DEG结果的摘要信息DEG_list$Summary写入到./results/DEGcounts.xls文件中以制表符分隔不加引号不写入行名。
## 5.2 Plot DEG results
## Filter and plot DEG results for up and down regulated genes. The definition of up
## and down is given in the corresponding help file. To open it, type ?filterDEGs in
## the R console.
edgeDF - read.delim(results/edgeRglm_allcomp.xls, row.names 1, check.names FALSE)
pdf(results/DEGcounts.pdf)
DEG_list - filterDEGs(degDF edgeDF, filter c(Fold 2, FDR 20))
dev.off()
write.table(DEG_list$Summary, ./results/DEGcounts.xls, quote FALSE, sep \t, row.names FALSE)计算和绘制差异表达基因DEG集合的Venn图
利用overLapper和vennPlot函数计算和绘制差异表达基因集合的Venn图。首先计算上调和下调基因集合的Venn图交集并将结果保存。然后根据设定绘制Venn图并将图形保存到PDF文件中。主要包括以下几个步骤 创建Venn图数据通过overLapper函数计算DEG集合的Venn图交集。首先使用DEG_list$Up[6:9]作为输入表示选取DEG结果中上调基因的第6至第9个集合然后将结果保存在vennsetup变量中。接着使用DEG_list$Down[6:9]作为输入表示选取DEG结果中下调基因的第6至第9个集合将结果保存在vennsetdown变量中。 绘制Venn图通过pdf函数创建一个PDF文件“results/vennplot.pdf”然后使用vennPlot函数绘制Venn图。将需要绘制的Venn图数据传入list函数中并通过mymain和mysub参数指定主标题和副标题为空字符串。此外colmode参数设为2表示使用两种颜色蓝色和红色来区分上调和下调的基因集合。 结束绘图通过dev.off函数结束图形设备保存Venn图到PDF文件中。
## 5.3 Venn diagrams of DEG sets
## The overLapper function can compute Venn intersects for large numbers of sample
## sets (up to 20 or more) and plots 2-5 way Venn diagrams. A useful feature is the
## possibility to combine the counts from several Venn comparisons with the same
## number of sample sets in a single Venn diagram (here for 4 up and down DEG sets).
vennsetup - overLapper(DEG_list$Up[6:9], type vennsets)
vennsetdown - overLapper(DEG_list$Down[6:9], type vennsets)
pdf(results/vennplot.pdf)
vennPlot(list(vennsetup, vennsetdown), mymain , mysub , colmode 2, ccol c(blue, red))
dev.off()GO富集分析
准备工作
进行基因-基因本体Gene Ontology, GO富集分析的准备工作。首先选择并获取BioMart数据库然后从数据库中获取基因到GO的映射关系并对结果进行预处理。最后创建基因到GO的CATdb对象用于后续的GO富集分析。主要包括以下几个步骤 选择和获取BioMart数据库通过使用listMarts函数列出可用的BioMart数据库选择目标数据库。在这里通过指定host参数为https://plants.ensembl.org来获取与植物相关的数据库。然后使用useMart函数选择目标数据库和数据集。 获取基因到GO映射通过使用getBM函数从选择的BioMart数据库中获取基因到GO的映射关系。指定所需的属性attributes为go_idGO标识符、“tair_locus”基因标识符和namespace_1003GO的命名空间。将获取的结果保存在go变量中。 数据预处理对获取的基因到GO映射进行预处理。首先去除命名空间为空的条目。然后将命名空间的值转换为缩写形式F表示分子功能P表示生物过程C表示细胞组分。最后将结果保存在文件GOannotationsBiomart_mod.txt中。 创建基因到GO的CATdb对象通过使用makeCATdb函数创建基因到GO的CATdb对象。指定文件路径、列号以及其他必要的参数。将创建的CATdb对象保存在文件catdb.RData中。
## 6. GO term enrichment analysis
#####################################################
## 6.1 Obtain gene-to-GO mappings
## The following shows how to obtain gene-to-GO mappings from biomaRt (here for
## A. thaliana) and how to organize them for the downstream GO term enrichment
## analysis. Alternatively, the gene-to-GO mappings can be obtained for many
## organisms from Bioconductor’s *.db genome annotation packages or GO
## annotation files provided by various genome databases. For each annotation this
## relatively slow preprocessing step needs to be performed only once. Subsequently,
## the preprocessed data can be loaded with the load function as shown in the next
## subsection.
library(biomaRt)
listMarts() # To choose BioMart database
listMarts(host https://plants.ensembl.org)
m - useMart(plants_mart, host https://plants.ensembl.org)
listDatasets(m)
m - useMart(plants_mart, dataset athaliana_eg_gene, host https://plants.ensembl.org)
listAttributes(m) # Choose data types you want to download
go - getBM(attributes c(go_id, tair_locus, namespace_1003), mart m)
# If download fail, you can load the following Rdata.
#load(Lab_4 files/go.Rdata)
go - go[go[, 3] ! , ]
go[, 3] - as.character(go[, 3])
go[go[, 3] molecular_function, 3] - F
go[go[, 3] biological_process, 3] - P
go[go[, 3] cellular_component, 3] - C
go[1:4, ]
# dir.create(./data/GO)
write.table(go, data/GO/GOannotationsBiomart_mod.txt, quote FALSE, row.names FALSE, col.names FALSE, sep \t)
catdb - makeCATdb(myfile data/GO/GOannotationsBiomart_mod.txt, lib NULL, org , colno c(1, 2, 3), idconv NULL)
save(catdb, file data/GO/catdb.RData)批量进行GO富集分析
这段代码主要用于进行批量的基因本体Gene Ontology, GO富集分析首先根据DEG结果定义DEG集合然后执行批量的GO富集分析和GO slim富集分析。最后将结果保存在相应的变量中。具体步骤如下 载入所需的包和数据通过加载biomaRt包和预处理好的基因到GO的CATdb对象从之前的代码段中加载。同时也加载了之前进行差异表达分析得到的DEG结果DEG_list。 定义DEG集合根据DEG结果将上调和下调基因分别命名为名称_up_down、名称_up和名称_down的命名空间。创建DEG集合DEGlist将这些命名空间组合起来并移除长度为0的集合。 执行批量的GO富集分析使用GOCluster_Report函数对DEG集合进行批量的GO富集分析。设置方法参数method为all表示返回所有通过设定的p-value阈值cutoff的GO术语。指定id_type参数为gene表示基因标识符类型。还可以设置其他参数如聚类阈值CLSZ、GO分类gocats和记录指定的GO术语recordSpecGO。将结果保存在BatchResult变量中。 获取GO slim向量通过使用biomaRt包获取特定生物体的GO slim向量。选择合适的BioMart数据库“plants_mart”和数据集“athaliana_eg_gene”然后使用getBM函数获取goslim_goa_accession属性并将结果转换为字符向量。 执行GO slim富集分析使用GOCluster_Report函数对DEG集合进行GO slim富集分析。设置方法参数method为slim表示仅返回在goslimvec向量中指定的GO术语。其他参数的设置与批量GO富集分析类似。将结果保存在BatchResultslim变量中。
## 6.2 Batch GO term enrichment analysis
## Apply the enrichment analysis to the DEG sets obtained the above differential
## expression analysis. Note, in the following example the FDR filter is set here to an
## unreasonably high value, simply because of the small size of the toy data set used ## in
## this vignette. Batch enrichment analysis of many gene sets is performed with the
## function. When methodall, it returns all GO terms passing the p-value cutoff
## specified under the cutoff arguments. When methodslim, it returns only the GO
## terms specified under the myslimv argument. The given example shows how a GO
## slim vector for a specific organism can be obtained from BioMart.
library(biomaRt)
load(data/GO/catdb.RData)
DEG_list - filterDEGs(degDF edgeDF, filter c(Fold 2, FDR 50), plot FALSE)
up_down - DEG_list$UporDown
names(up_down) - paste(names(up_down), _up_down, sep )
up - DEG_list$Up
names(up) - paste(names(up), _up, sep )
down - DEG_list$Down
names(down) - paste(names(down), _down, sep )
DEGlist - c(up_down, up, down)
DEGlist - DEGlist[sapply(DEGlist, length) 0]
BatchResult - GOCluster_Report(catdb catdb, setlist DEGlist, method all, id_type gene, CLSZ 2, cutoff 0.9, gocats c(MF, BP, CC), recordSpecGO NULL)
library(biomaRt)
m - useMart(plants_mart, dataset athaliana_eg_gene, host https://plants.ensembl.org)
goslimvec - as.character(getBM(attributes c(goslim_goa_accession), mart m)[, 1])
BatchResultslim - GOCluster_Report(catdb catdb, setlist DEGlist, method slim, id_type gene, myslimv goslimvec, CLSZ 10, cutoff 0.01, gocats c(MF, BP, CC), recordSpecGO NULL)绘制批量GO富集分析的结果
使用goBarplot函数绘制批量GO富集分析结果的条形图。首先选择感兴趣的子集然后分别绘制不同GO类别的条形图。具体步骤如下 子集选择首先从BatchResultslim数据框中选择与M6-V6_up_down匹配的行将结果存储在gos变量中。然后将整个BatchResultslim数据框存储在gos变量中以便后续绘制。 绘制MF分子功能类别的GO条形图通过调用goBarplot函数绘制MF类别的GO条形图。将gos作为输入数据框并设置gocat参数为MF。将结果保存为PDF文件。 绘制BP生物过程类别的GO条形图通过再次调用goBarplot函数绘制BP类别的GO条形图。将gos作为输入数据框并设置gocat参数为BP。 绘制CC细胞组分类别的GO条形图通过再次调用goBarplot函数绘制CC类别的GO条形图。将gos作为输入数据框并设置gocat参数为CC。
## 6.3 Plot batch GO term results
## The data.frame generated by GOCluster can be plotted with the goBarplot
## function. Because of the variable size of the sample sets, it may not always be
## desirable to show the results from different DEG sets in the same bar plot. Plotting
## single sample sets is achieved by subsetting the input data frame as shown in the
## first line of the following example.
gos - BatchResultslim[grep(M6-V6_up_down, BatchResultslim$CLID),
]
gos - BatchResultslim
pdf(GOslimbarplotMF.pdf, height 8, width 10)
goBarplot(gos, gocat MF)
dev.off()
goBarplot(gos, gocat BP)
goBarplot(gos, gocat CC)这里以分子功能为例。 层次聚类和热图绘制
使用pheatmap库进行层次聚类和热图绘制。首先提取感兴趣基因的表达矩阵子集然后基于Pearson相关系数计算距离并进行层次聚类最后绘制热图以可视化聚类结果。具体步骤如下 导入必要的库通过加载pheatmap库准备进行层次聚类和热图分析。 提取差异表达基因DEGs的表达矩阵从上述差异表达分析中确定的DEGs中提取基因的rlog转换后的表达矩阵。这里使用了DEG_list[[1]]作为DEGs的标识通过将其转换为字符向量并去除重复项得到基因的唯一标识符。 提取感兴趣基因的表达值从rlog转换后的表达矩阵中提取与感兴趣基因的唯一标识符相对应的行得到一个子集矩阵y。 绘制热图通过调用pheatmap函数绘制热图。设置scale参数为row对行进行标准化设置clustering_distance_rows参数和clustering_distance_cols参数为correlation使用基于Pearson相关系数的距离度量进行行和列的层次聚类。将结果保存为PDF文件。
## 7. Clustering and heat maps
#####################################################
## The following example performs hierarchical clustering on the rlog transformed
## expression matrix subsetted by the DEGs identified in the above differential
## expression analysis. It uses a Pearson correlation-based distance measure and
## complete linkage for cluster joining.
library(pheatmap)
geneids - unique(as.character(unlist(DEG_list[[1]])))
y - assay(rlog(dds))[geneids, ]
pdf(heatmap1.pdf)
pheatmap(y, scale row, clustering_distance_rows correlation, clustering_distance_cols correlation)
dev.off()
