此工作流的目标是检测已知的腐霉属从宾夕法尼亚州的大豆地里收集的几个样本中发现了一些物种。腐霉属是一种真核微生物，可以帮助或伤害植物取决于物种。许多腐霉属物种导致特定植物的根腐病。然而，有些种类的腐霉属被用作生物防治剂，以防止作物上病原体的生长。该工作流与Coffua执行的相同数据集的分析并行et al。发表在《植物疾病》杂志(2016)上。

在这里，细胞色素c氧化酶亚基1 (COI)基因被用作系统发育标记来识别物种。COI基因是所有真核生物线粒体基因组的一部分。在这个工作流程的第1部分，目标是设计最佳引物来区分所有的腐霉属使用COI基因的物种。第一步是下载已知的COI基因序列腐霉属并将它们导入序列数据库，如下图所示。

所有路径都是相对于已安装的数据集data_dir <- system。文件(“extdata”,包=“BigBioSeqData”)图书馆(破解)

##加载所需的包:生物字符串

##加载所需的包:BiocGenerics

加载所需的包:平行

## ##附加包:'BiocGenerics'

以下对象从'package:parallel'中被屏蔽:## ## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ， ## clusterExport, clusterMap, parApply, parCapply, parApply， ## parapplylb, parRapply, parSapply, parSapplyLB

以下对象从“package:stats”中被屏蔽:## ## IQR, mad, xtabs

以下对象从'package:base'被屏蔽:## ## Filter, Find, Map, Position, Reduce, anyduplicate, append， ## as.data.frame, cbind, colnames, do. frame。调用，duplicate, eval， ## evalq, get, grep, grepl, intersect, is。Unsorted, lapply， ## length, mapply, match, mget, order, paste, pmax, pmax.int， ## pmin, pmin.int, rank, rbind, rownames, sapply, setdiff, sort， ## table, tapply, union, unique, unsplit, which, which。马克斯,# # which.min

加载所需的包:S4Vectors

##加载所需的包:stats4

## ##附加包:'S4Vectors'

以下对象从'package:base'中被屏蔽:## ## colMeans, colsum, expand。网格,rowMeans rowSums

##加载所需的包:伊朗

加载所需的包:XVector

加载所需的包:RSQLite

加载所需的包:DBI

#创建到磁盘上SQLite数据库的连接dbConn <- dbConnect(SQLite()， "./COIgenes.sqlite") #新数据库文件的路径#从GenBank格式化文件Seqs2DB导入序列(paste(data_dir， "/Pythium_spp_COI. txt)gb"， sep="")， type="GenBank"， dbFile=dbConn, identifier="Pythium") #查看构建的数据库表BrowseDB(dbConn) #检索导入的序列dna <- SearchDB(dbConn) dna #根据它们的翻译对齐序列dna <- AlignTranslation(dna) dna #在web浏览器中显示序列BrowseSeqs(dna) #显示与第一个序列的差异BrowseSeqs(dna, highlight=1) #显示与一致序列BrowseSeqs(dna，#改变一致程度BrowseSeqs(DNA, highlight=0, threshold=0.2) #注意大多数序列的共同模式pattern <- DNAStringSet("TAGATTTAGCWATTTTTAGTTTACA") BrowseSeqs(DNA, patterns=pattern) #蛋白质序列非常相似AA <- AlignTranslation(DNA, asAAStringSet=TRUE)突出= 1)#选择移码参考校正参考< -翻译(dna[11]) # # 11 #序列正确的序列# 12中的移码正确的< - CorrectFrameshifts (myXStringSet = dna [12], myAAStringSet = REF, type =“两”)[12]正确的dna序列正确< - $ # # 12 # 11现在相同的dna序列< - AlignTranslation (dna) BrowseSeqs (dna,突出= 11)#确定引物设计集群d < - DistanceMatrix (dna)的(d) #一个对称矩阵c < - IdClusters (d =“UPGMA方法,截止= 0.05,显示= TRUE)

Add2DB(data.frame(identifier=paste("cluster"， c$cluster, sep=""))， dbConn) BrowseDB(dbConn) #为下一代测序引物设计引物<- DesignSignatures(dbConn, type="sequence"， resolution=5, levels=5, minProductSize=400, maxProductSize=800, annealingTemp=55, maxpermutation =8)引物[1，]#最高评分引物设置#突出引物的目标位点BrowseSeqs(DNA，patterns=c(DNAStringSet(primers[1,1])， reverseComplement(DNAStringSet(primers[1,2]))))

工作流程的第2部分使用了从宾夕法尼亚州的几个地方获得的COI基因序列。这些DNA序列连同它们相应的质量分数以FASTQ格式存储。导入后，第一步是修剪序列，以便只保留高质量的中心区域。可能属于的序列的子集腐霉属物种将通过存在一个共同的保守区域来识别腐霉属．此分析将分批执行，因此所有的序列不需要同时装入内存。

#从压缩的FASTQ序列文件路径<- paste(data_dir， "/FASTQ/"， sep="") files <- list.files(path) samples <- substring(files, first=1, last=nchar(files) - 6) for (i in seq_along(files)) {cat(samples[i]， ":\n"， sep="") Seqs2DB(paste(path, files[i]， sep="")， type="FASTQ"， dbFile=dbConn, identifier=samples[i]， tblName="Reads")} #确定边界的函数#高质量中心区域边界<- Function (probs, thresh=0.001，width=21){#计算移动平均填充<- floor(width/2) probs <- c(rep(thresh, padding)， probs, rep(thresh, padding)) probs <- filter(probs, rep(1/width, width)) #查找阈值以上区域w <- which(probs < thresh) - padding if (length(w)==0) w <- NA return(c(w[1]， w[length(w)]))} #按质量修剪序列并识别#属于Pythium属的子集nSeqs <- SearchDB(dbConn, tbl="Reads"， count=TRUE，verbose=FALSE) offset <- 0 ends <- starting <- counts <- integer(nSeqs) fprimer <- DNAString(" tcawcwmgatggcttttcaac ") rprimer <- DNAString("") pBar <- txtProgressBar(max=nSeqs, style=3) while (offset < nSeqs){#选择一组序列dna <- SearchDB(dbConn, tbl="Reads"， type="QualityScaledXStringSet"， limit=paste(offset, 1e4, sep="，")， verbose=FALSE) #将质量分数转换为错误概率probs <- as(quality(dna)， "NumericList")端点<- sapply(probs， "#保存结果以备以后使用index <- (offset + 1):(offset + length(dna)) starts[index] <- ifelse(endpoints[1，] >= 38L, endpoints[1，]， 38L) #前向引物结束后的第一个碱基[index] <- ifelse(endpoints[2，] >= starts[index] - 1L) #没有高质量的碱基#在python序列中查找预期的模式计数[index] <- vcountPattern(pattern[[1]]， subject=dna, max. index)。= 4,不匹配。indels=TRUE, fixed="subject") #允许模糊偏移<- offset + 1e4 setTxtProgressBar(pBar, ifelse(offset > nSeqs, nSeqs, offset))} #将结果添加到数据库结果的新列<- data.frame(start= started, end=ends, count=counts) Add2DB(results, dbFile=dbConn, tblName="Reads"， verbose=FALSE) BrowseDB(dbConn, tblName="Reads"，limit=1000) #将每个样本中的读取按百分比标识(i in seq_along(samples)) {cat(samples[i]) #选择中等长度的序列dna <- SearchDB(dbConn, tblName=" reads "， identifier=samples[i]， clause="count > 0 and (end - start + 1) >= 100"， verbose=FALSE) cat(":"， length(dna)， "sequences") #将序列修剪到高质量区域索引<- as.numeric(names(dna)) dna <- subseq(dna, start=starts[index]，end=ends[index]) #聚集没有距离矩阵的序列集群<- IdClusters(myXStringSet=dna, method="inexact"， cutoff=0.03， # > 97% identity verbose=FALSE) #添加集群编号到数据库Add2DB(clusters, dbFile=dbConn, tblName="Reads"， verbose=FALSE) cat("，"， length(unique(clusters[， 1]))， "clusters\n")} #现在数据库包含一个集群列BrowseDB(dbConn, tblName="Reads"， limit=1000, clause=" Cluster is not NULL")

在工作流的第3部分中，将每个序列集群的代表与已知的进行比较腐霉属物种。这种分析的目的是确定每个样本中存在的哪些生物与已知物种相似。已知的菌种被分为两组:被用作生物防治剂的菌种(好菌种)和已知的植物病原体菌种(坏菌种)。

ids <- IdentifyByRank(dbConn, add2tbl=TRUE) lens <- idlength(dbConn, add2tbl=TRUE) BrowseDB(dbConn) #将Pythium菌株分为好组和坏组生物控制<- c('Pythium oligandrum'， 'Pythium nunn'， 'Pythium periplocum') pathogen <- c('Pythium acanthicum'， #草莓:'Pythium rostratum'， 'Pythium middletonii'， 'Pythium aristsporum '， #草/杂粮:' graminicola'， ' okanoganense'， 'Pythium paddicum'， 'Pythium volutum'， 'Pythium arrhenomanes'， 'Pythium buismaniae'， #花:'Pythium spinosum'， 'Pythium mastophorum'， 'Pythium splendens'， 'Pythium violae'， #胡萝卜:'Pythium parecandrum '， 'Pythium sulcatum'， 'Pythium dissotocum'， #土豆:'Pythium scleroteichum'， 'Pythium myriotylum'， 'Pythium heterothicum '， #生菜:'Pythium trachephilum '， 'Pythium ultimum'， #多种植物:'Pythium irregulare'， 'Pythium aphanidermatum'， 'Pythium debaryanum'， 'Pythium sylvaticum') #从每个物种中选择最长的序列物种<- SearchDB(dbConn, nameBy="identifier"， "'"， sep=""， collapse="， ")， clause=paste("identifier in ("， paste("'"， c(生物控制，病原体)，"'"，sep="， ") group by identifier having max(bases)"， sep="")) #选择每个群集中最长的序列dna <- SearchDB(dbConn, identifier="DauphinFarm"， #选择一个样本tblName="Reads"，条款= "集群不是零group by集群马克斯(结束-开始)”)#修剪中部地区高质量指数< - as.numeric(名字(dna)) dna < - subseq (dna,开始=(指数)开始,结束=结束(指数))#创建一个树与已知和未知物种结合< - AlignSeqs (c (dna,物种))的dist < - DistanceMatrix (verbose = FALSE,结合修正=“Jukes-Cantor”)树< - IdClusters(经销、方法=“NJ”,#邻居加入asDendrogram = TRUE, verbose = FALSE)情节(树,nodePar =列表(lab.cex = 0.5,pch = NA))

#根据其致病性为已知物种着色tree_colored <- dendrapply(tree, function(x) {if (is.leaf(x)) {if (attr(x， "label") %in% pathogen) {attr(x， "label") <- list(col="red")} else if (attr(x， "label") %in% biocontrol) {attr(x， "edgePar") <- list(col="green")} #删除标签attr(x， "label") <- ""} return(x)}) plot(tree_colored)

# Disconnect from sequence database dbDisconnect(dbConn) #永久删除数据库unlink("./COIgenes.sqlite") # optional!