systemPipeRdata::genWorkenvir(工作流= "rnaseq"， mydirname = "rnaseq") setwd("rnaseq")

2.1.1实验定义由目标文件

的目标文件定义了所有FASTQ文件和分析工作流的样本比较。

Targetspath <- system。文件("extdata", "targetsPE.txt", package = "systemPipeR") targets <- read.delim(targetspath, comment.char = "#") targets[1:4, -c(5, 6)]

## FileName1 FileName2 ## 1. /data/SRR446027_1.fastq.gz ./data/SRR446028_1.fastq.gz ## 3 ./data/SRR446029_1.fastq.gz ./data/SRR446029_2.fastq.gz ## 4 ./data/SRR446029_1.fastq.gz ./data/SRR446029_2.fastq.gz ## SampleName Factor日期## 1 M1A M1 2012年3月23日## 2 M1B M1 2012年3月23日## 3 A1A A1 2012年3月23日## 4 A1B A1 2012年3月23日

要使用自定义数据，用户需要生成一个目标包含它们自己的FASTQ文件路径的文件。

3.2.1之上预处理与preprocessReads函数

这个函数preprocessReads允许应用预定义的或自定义的读预处理函数到文件中引用的所有FASTQ文件SYSargsList容器，例如质量过滤或适配器修整例程。在内部,preprocessReads使用FastqStreamer函数从ShortRead包以一种内存效率高的方式通过大型FASTQ文件。属性执行适配器调整trimLRPatterns函数从Biostrings包中。

在这里，我们将此步骤附加到SYSargsList之前创建的对象。对象上定义了所有参数preprocessReads / preprocessReads-pe.yml文件。

appendStep(sal) <- SYSargsList(step_name = "预处理"，targets = "targetsPE.txt"， dir = TRUE, wf_file = "preprocessReads/preprocessReads-pe. txt")， input_file = "preprocessReads/preprocessReads-pe. cwl"。Yml”，dir_path =系统。文件("extdata/cwl", package = "systemPipeR"), inputvars = c(FileName1 = "_FASTQ_PATH1_", FileName2 = "_FASTQ_PATH2_", SampleName = "_SampleName_"), dependency = c("load_SPR"))

预处理步骤完成后，将输出文件文件可以用来生成新的目标文件，其中包含裁剪后的FASTQ文件的路径。新的目标信息可以用于下一个工作流步骤实例，如。用裁剪好的FASTQ文件运行NGS对齐。的appendStep函数自动处理步骤之间的连接。详情请查看下一步。

下面的示例展示如何设计自定义读取“preprocessReads”方法提供的实用程序ShortRead包，然后在批处理模式下使用“preprocessReads”函数。在这里，可以替换上使用的函数预处理步骤和修改萨尔对象。因为它是一个自定义函数，所以必须将该部分保存在R对象中，并且在内部保存preprocessReads.doc.R正在加载自定义函数。如果R对象以不同的名称保存(这里“参数/ customFCT。RData”)，请将有关资料替换至preprocessReads.doc.R．

请注意，此步骤没有添加到工作流中，这里只是为了演示。

首先，我们在工作流中定义了自定义函数:

appendStep(sal) <- LineWise(code = {filterFct <- function(fq, cutoff = 20, Nexceptions = 0) {qcount <- rowsum (as(quality(fq)， "matrix") <= cutoff, na。fq[qcount <= Nexceptions]} save(list = ls()， file = "param/ customcft . rdata ")}， step_name = "custom_preprocessing_function"， dependency = "预处理")

之后，我们可以编辑输入参数:

yamlinput(sal， "预处理")$Fct yamlinput(sal， "预处理"，"Fct") <- "'filterFct(fq, cutoff=20, Nexceptions=0)'" yamlinput(sal， "预处理")$Fct ##检查新函数cmdlist(sal， "预处理"，targets = 1) ##检查命令行是否更新成功

3.2.2阅读修剪与Trimmomatic

Trimmomatic软件(Bolger, Lohse, and Usadel 2014)对Illumina成对端和单端数据执行各种有用的修整任务。下面是一个如何使用参数模板文件来修剪FASTQ文件的示例。

该步骤是可选的。

appendStep(sal) <- SYSargsList(step_name = "trim "， targets = "targetsPE.txt"， wf_file = "trimmomatic/trimmomatic-pe. txt")Cwl "， input_file = "trimmomatic/trimmomatic-pe。Yml”，dir_path =系统。文件("extdata/cwl", package = "systemPipeR"), inputvars = c(FileName1 = "_FASTQ_PATH1_", FileName2 = "_FASTQ_PATH2_", SampleName = "_SampleName_"), dependency = "load_SPR", run_step = "optional")

3.2.3FASTQ质量报告

以下seeFastq而且seeFastqPlot函数为一组FASTQ文件生成并绘制一系列有用的质量统计数据，包括每个周期质量箱形图、基本比例、基本水平质量趋势、相对k-mer多样性、长度和读取的出现分布、超过质量截断点的读取数和平均质量分布。结果被写入一个名为fastqReport.png．

appendStep(sal) <- LineWise(code = {fastq <- getColumn(sal, step = "预处理"，"targetsWF"， column = 1) fqlist <- seeFastq(fastq = fastq, batchsize = 10000, klength = 8) png("./results/fastqReport.png"， height = 162, width = 288 * length(fqlist)) seeFastqPlot(fqlist) dev.off()}， step_name = "fastq_report"， dependency = "预处理")

3.3.1读取映射HISAT2

下面的步骤将演示如何使用短读取对齐器Hisat2(Kim, Langmead, and Salzberg 2015)．首先,Hisat2需要创建索引。

appendStep(sal) <- SYSargsList(step_name = "hisat2_index"， dir = FALSE, targets = NULL, wf_file = "hisat2/hisat2-index. index. "Cwl "， input_file = "hisat2/hisat2-index. index. "yml"， dir_path = "param/cwl"， dependency = "load_SPR")

3.3.2HISAT2映射

中定义了对准器的参数设置workflow_hisat2-pe.cwl而且workflow_hisat2-pe.yml文件。下面演示如何构造相应的SYSargsList对象。

appendStep(sal) <- SYSargsList(step_name = "hisat2_mapping"， dir = TRUE, targets = "预处理"，wf_file = "工作流-hisat2/workflow_hisat2-pe. zip ")Cwl "， input_file = " workflow_hisat2 /workflow_hisat2-pe. Cwl "， input_file = " workflow_hisat2-pe. Cwl "。yml"， dir_path = "param/cwl"， inputvars = c(preprocessReads_1 = "_FASTQ_PATH1_"， preprocessReads_2 = "_FASTQ_PATH2_"， SampleName = "_SampleName_")， rm_targets_col = c("FileName1"， "FileName2")， dependency = c("预处理"，"hisat2_index")))

要仔细检查每个示例的命令行，请使用以下命令:

Cmdlist (sal, step = "hisat2_mapping"， targets = 1)

3.3.3读取和对齐统计信息

下面概述了每个示例中的读取数以及其中有多少与引用对齐。

appendStep(sal) <- LineWise(code = {fqpaths <- getColumn(sal, step = "预处理"，"targetsWF"， column = "FileName1") bampaths <- getColumn(sal, step = "hisat2_mapping"， "outfiles"， column = "samtools_sort_bam") read_statsDF <- alignStats(args = bampaths, fqpaths = fqpaths, pairEnd = TRUE) write。table(read_statsDF， "results/alignStats.xls"， row.names = FALSE, quote = FALSE, sep = "\t")}， step_name = "align_stats"， dependency = "hisat2_mapping")

3.5.1创建基因注释数据库

appendStep(sal) <- LineWise(code = {library(GenomicFeatures) txdb <- suppressWarnings(makeTxDbFromGFF(file = "data/tair10. txt ")gff"， format = "gff"， dataSource = "TAIR"， organism = "拟南芥"))saveDb(txdb, file = "./data/tair10.sqlite")}， step_name = "create_db"， dependency = "hisat2_mapping")

3.5.2阅读计数summarizeOverlaps采用多核并行模式

appendStep (sal) < - LineWise(代码={库(GenomicFeatures)库(BiocParallel) txdb < - loadDb(“。/数据/ tair10.sqlite”)outpaths < - getColumn (sal,一步=“hisat2_mapping”、“外部档案”,列=“samtools_sort_bam”)eByg < - exonsBy (txdb, = c(“基因”))bfl < - BamFileList (outpaths yieldSize = 50000,指数=字符())multicoreParam < - multicoreParam(工人= 4)登记注册()(multicoreParam) counteByg < - bplapply (bfl、函数(x) summarizeOverlaps (eByg x =“联盟”模式,忽视。strand = TRUE, inter。feature = FALSE, singleEnd = FALSE, BPPARAM = multicoreream)) 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") ##创建一个summarizeexperimental对象colData <- data.frame(row.names = SampleName(sal， "hisat2_mapping")， condition = getColumn(sal， "hisat2_mapping"， position = "targetsWF"， column = "Factor")) colData$condition <- Factor (colData$condition) countDF_se <- summarizeexperimental:: summarizeexperimental (assays = countDFeByg，##对象SE(sal， "read_counting") <- countDF_se}， step_name = "read_counting"， dependency = "create_db")

当提供BamFileList在上面的例子中，summarizeOverlaps默认使用的方法bplapply并使用BiocParallel包中的寄存器接口。如果没有设置工人的数量，MulticoreParam将使用返回的核数并行:detectCores ()．如需更多信息，请查看帮助(“summarizeOverlaps”)文档。

注意，对于大多数统计差异表达或丰度分析方法，如刨边机或DESeq2时，应使用原始计数值作为输入。RPKM值的使用应限制在某些用户所需的特殊应用中，如。手动比较不同基因或特征之间的表达水平。

3.5.3样本相关分析

以下计算样本的斯皮尔曼相关系数从rlog方法生成的转换表达式值DESeq2包中。转换为距离矩阵后，使用hclust函数和结果被绘制为树状图(也见文件sample_tree.png）.

appendStep(sal) <- LineWise(code = {library(DESeq2, quiet = TRUE) library(猿，警告。冲突= FALSE) ##提取摘要实验对象se <- se (sal， "read_counting") dds <- DESeqDataSet(se, design = ~condition) d <- cor(化验(rlog(dds))， method = "spearman") hc <- hclust(dist(1 - d)) png("results/sample_tree.png") plot.phylo(as.phylo(hc)， type = "p"， edge. png)。Col = "blue"， edge。width = 2, show.node.label = TRUE, no。margin = TRUE) dev.off()}， step_name = "sample_tree"， dependency = "read_counting")

3.6.1运行刨边机

appendStep(sal) <- LineWise(code = {library(edgeR) countDF <- read.delim("results/countDFeByg.xls"，行。names = 1, check.names = FALSE) cmp <- readComp(stepsWF(sal)[["hisat2_mapping"]]]， format = "matrix"， delim = "-") edgeDF <- run_edgeR(countDF = countDF, targets = targetsWF(sal)[["hisat2_mapping"]]]， cmp = cmp[[1]]， independent = FALSE, mdsplot = "read_counting")}， step_name = " run_edgeR "， dependency = "read_counting")

操作增加基因描述

appendStep(sal) <- LineWise(code = {library("biomaRt") m <- useMart("plants_mart"， dataset = "athaliana_eg_gene"， host = "https://plants.ensembl.org") desc <- getBM(attributes = c(" taair_locus "， "description")， mart = m) desc <- desc[!duplicate (desc[， 1])，] desc <- as。字符(desc[， 2])名称(desc) <- 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)}， step_name = "custom_annot"， dependency = "run_edger")

3.6.3DEG结果

对上调和下调基因进行筛选和绘制DEG结果。的定义向上而且下来在相应的帮助文件中给出。要打开它，请键入filterDEGs ?在R控制台。

appendStep(sal) <- LineWise(code = {edgeDF <- read.delim("results/ edgerglm_allcomm .xls"， row.names = 1, check.names = FALSE) png("results/DEGcounts.png") 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)}， step_name = "filter_degs"， dependency = "custom_annot")

3.6.4DEG集的维恩图

的重叠函数可以计算大量样本集(多达20个或更多)的维恩相交，并绘制2-5路维恩图。一个有用的特性是可以将多个维恩比较的计数与单个维恩图中相同数量的样本集组合在一起(这里是4个上下DEG集)。

appendStep(sal) <- LineWise(code = {vennsetup <- overLapper(DEG_list$Up[6:9]， type = "vennsets") vennsetdown <- overLapper(DEG_list$Down[6:9]， type = "vennsets") png("results/ vennPlot .png") vennPlot(list(vennsetup, vennsetdown)， mymain = ""， mysub = ""， colmode = 2, ccol = c("blue"， "red")) dev.off()}， step_name = "venn_diagram"， dependency = "filter_degs")

3.7.1获得基因到go的映射

下面展示如何从。获取基因到go映射biomaRt(在这里答:芥)以及如何组织它们进行下游GO术语富集分析。另外，许多生物的基因到氧化石墨烯的映射可以从Bioconductor公司获得* .db各基因组数据库提供的基因组注释包或GO注释文件。对于每个注释，这个相对缓慢的预处理步骤只需要执行一次。方法加载预处理数据负载函数如下小节所示。

appendStep (sal) < - LineWise(代码={图书馆(biomaRt) # listMarts() #选择biomaRt数据库# listMarts(主机= ' plants.ensembl.org ') m < - useMart(“plants_mart”,主机= " https://plants.ensembl.org ") # listDatasets (m) m < - useMart(数据集“plants_mart”=“athaliana_eg_gene”,主机= " https://plants.ensembl.org ") # listAttributes (m) #选择数据类型你想要#下载去< - getBM(属性= c(“go_id”、“tair_locus”,“namespace_1003”),集市= m)去< - [[3 ] != "", ] (3) <——去。字符(去[3])去(去[3]= =“molecular_function”,3)<——“F”(去[3]= =“biological_process”,3)<——P[[3]去= =“cellular_component”,3]< -“C”(1:4,)如果(! dir.exists(“。/数据/去”))dir.create(“。/数据/去”)。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 ")}， step_name = "get_go_annot"， dependency = "filter_degs")

3.7.2章批量GO项富集分析

将富集分析应用于上述差分表达式分析得到的DEG集。注意，在下面的例子中罗斯福过滤器在这里被设置为一个不合理的高值，仅仅是因为在这个小插图中使用的玩具数据集很小。利用该函数对多个基因集进行批量富集分析。当方法=所有函数中指定的p值截断点将返回所有GO项截止参数。当方法=苗条方法下指定的GO项myslimv论点。所给的例子说明了如何从一个特定的生物体获得一个GO细矢量BioMart．

appendStep (sal) < - LineWise(代码={图书馆(biomaRt)负载(“数据/去/ catdb.RData”)DEG_list < - filterDEGs (degDF = edgeDF过滤器= c(折叠= 2,罗斯福= 50),情节= FALSE) up_down < - DEG_list美元UporDown名称(up_down) < -粘贴(名字(up_down),“_up_down”,9 = " ")<——美元DEG_list名称()< -粘贴(名称(最多),“_up”,9 = " ")<——DEG_list美元下跌的名字(下降)< -粘贴(名字(下降),“_down”,9 = " ")DEGlist < - c (up_down,上,下)DEGlist < - DEGlist[酸式焦磷酸钠(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) 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) write。table(BatchResultslim， "results/GOBatchSlim.xls"， row.names = FALSE, quote = FALSE, sep = "\t")}， step_name = "go_enrich"， dependency = "get_go_annot")

3.7.3Plot batch GO term结果

的data.frame生成的GOCluster可以用goBarplot函数。由于样本集的大小可变，在同一条形图中显示来自不同DEG集的结果可能并不总是可取的。绘制单个样本集是通过对输入数据帧进行细分来实现的，如下例的第一行所示。

appendStep(sal) <- LineWise(code = {gos <- BatchResultslim[grep(" 6m - v6_up_down "， BatchResultslim$CLID)，] gos <- BatchResultslim png("results/GOslimbarplotMF.png"， height = 8, width = 10) goBarplot(gos, gocat = "MF") goBarplot(gos, gocat = "BP") goBarplot(gos, gocat = "CC") dev.off()}， step_name = "go_plot"， dependency = "go_enrich")