1．1数据

您需要为您的野生型池提供BAM文件，为您的突变型池提供BAM文件，并为您的物种提供fasta格式的参考基因组。我们建议每个池至少包含20个个体，以确保有足够数量的重组进行测量。

1．2安装依赖关系

MMAPPR2依赖两个系统工具来发挥作用:Samtools和VEP。它们都必须安装在PATH中，以便由适当的函数找到。

git克隆https://github.com/Ensembl/ensembl-vep.git cd ensemble -vep perl INSTALL.pl -a ac -s {my_species}

Sys。setenv (PATH =粘贴(“/道路/ / Perlbrew”,Sys.getenv(“路径”),9 = ":"))

2.1设置参数

在我们的例子中，我们将只使用来自GRCz11斑马鱼参考基因组的黄金基因。

在这里，我们还配置VEPFlags对象来使用我们的示例fasta和GTF文件。更多信息见下文。

确保您的参考基因组与您将使用的VEP相同!除非您进行自定义，否则这将是在ensemble bl上可用的最新程序集。你也应该使用相同的基因组来校准你的测序数据。

BiocParallel::register(BiocParallel::MulticoreParam()) ##参见下面对BiocParallel库(MMAPPR2, quiet = TRUE)库(MMAPPR2data, quiet = TRUE)库(Rsamtools, quiet = TRUE)的解释

## ##附加包:'BiocGenerics'

以下对象从'package:MMAPPR2'中被屏蔽:## ## species, species<-

以下对象从'package:stats'被屏蔽:## ## IQR, mad, sd, var, xtabs

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

## ##附加包:'S4Vectors'

以下对象从'package:base'被屏蔽:## ## I，展开。网格,unname

## ##附件:“IRanges”

下面的对象从'package:MMAPPR2': ## ## distance被屏蔽

## ##附加包:'Biostrings'

以下对象从'package:base': ## ## strsplit被屏蔽

#这通常是自动配置的:vepFlags < - ensemblVEP: vepFlags(旗帜=列表(格式=“已”,# <——这是必要的vcf = FALSE, # <——以及这一物种= danio_rerio,数据库= FALSE, # <——这三个参数允许我们VEP脱机运行,fasta = goldenFasta(), # < -╯|你可能不需要人造石铺地面= goldenGFF (), # <------ ╯filter_common = TRUE, coding_only = TRUE #假设RNA-seq数据))参数< MmapprParam (refFasta = goldenFasta (), wtFiles = exampleWTbam (), mutFiles = exampleMutBam (),物种= 'danio_rerio'， vepFlags = vepFlags， ##可选outputFolder = tempOutputFolder()) ##可选

##注意:基因组'danio_rerio'已经存在，没有覆盖

2．2运行MMAPPR2

设置好参数后，运行管道应该像下面这样简单:

mappprdata <- mappr(param)

## ------------------------------------

## --------欢迎来到MMAPPR2 --------

## ------------------------------------

开始时间:2022-04-26 17:11:39

##输出文件夹:/tmp/RtmpZFm2nq/Rbuild75fb44bb0620/MMAPPR2/vignettes//tmp/RtmpliLD14/ mmappr2_22-04-26_17:11:39

读取BAM文件并生成欧氏距离数据…

为每个染色体生成最优黄土曲线…

正在识别连接区染色体…

峰值区域已成功识别

##使用SNP重采样细化峰值描述…

生成、分析和排序候选变量…

##警告有效。seqinfo (x,建议。GRanges对象包含位于序列18上的一个超出边界的范围。##注意，位于长度未知(NA)的序列上的范围或位于圆形序列上的##不被认为是超出范围的(使用## seqlength()和isCircular()来获取底层序列的长度和圆形标志##)。您可以使用trim()来修剪这些范围。更多信息见trim, genome icrange -method。

##在. call(。make_vcf_geno, filename, fixed, names(geno)， as.list(geno)，: ##将NULL指针转换为R NULL

正在写输出图和表…

## ##结束时间:22-04-26 17:11:51

## MMAPPR2运行时:11.93783秒

MMAPPR2管道也可以一步一步地运行:

md <- new('MmapprData'， param = param) ## calculateDistance()接受一个MmapprData对象postCalcDistMD <- calculateDistance(md) postLoessMD <- loessFit(postCalcDistMD) postPrePeakMD <- prePeak(postLoessMD) postPeakRefMD <- peak细化(postPrePeakMD) postCandidatesMD <- generatecanddates (postPeakRefMD)

##警告有效。seqinfo (x,建议。GRanges对象包含位于序列18上的一个超出边界的范围。##注意，位于长度未知(NA)的序列上的范围或位于圆形序列上的##不被认为是超出范围的(使用## seqlength()和isCircular()来获取底层序列的长度和圆形标志##)。您可以使用trim()来修剪这些范围。更多信息见trim, genome icrange -method。

##在. call(。make_vcf_geno, filename, fixed, names(geno)， as.list(geno)，: ##将NULL指针转换为R NULL

outputMmapprData (postCandidatesMD)

如果管道中途失效，则MmapprData对象被保存，然后您可以加载并用于检查和调试:

##输出文件夹的内容:cat(paste(system2('ls'， outputFolder(param(mmapprData))， stdout = TRUE))， sep = '\n')

# # 18。TSV # genome_plot .pdf ## mmappr2.log ## mmappr_data。RDS ## peak_plot .pdf

mdFile <- file.path(outputFolder(param(mmapprData))， 'mmappr_data.RDS') md <- readRDS(mdFile) md

引用fasta文件:## /home/biocbuild/ bs-3.15- bios /R/library/MMAPPR2data/extdata/slc24a5.fa.gz ## wtFiles: ##长度为1的BamFileList ## names(1): wt.bam ## mutFiles: ##长度为1的BamFileList ## names(1): mut。(6): format, species，…， filter_common, coding_only ## version: 105 ## scriptPath: ## refGenome: ## GmapGenome对象## genome: danio_rerio ##目录:/home/biocbuild/。local/share/gmap ##其他参数:## species ## danio_rerio ## distancePower ## 4 # peakIntervalWidth ## 0.95 ## minDepth ## 10 # homozygoteCutoff ## 0.95 # minBaseQuality ## 20 # minMapQuality ## 30 # loessOptResolution ## 0.001 ## loessOptCutFactor ## 0.1 # naCutoff ## 0 # outputFolder ## /tmp/RtmpliLD14/ mmappr2_22-04-26_17:11:39 # fileAggregation ## sum ## distance:##包含一个序列的欧氏距离数据##和其中一个的黄土回归数据##内存使用= 1 MB ##峰值:## 18:开始= -140，结束= 14091 ##密度函数计算##候选人:## $ ' 18 ' # GRanges对象12个范围和30个元数据列:## seqnames范围strand | Allele Consequence ##    |   ## 18:5744_A/C 18 5744 * | C missense_variant

2．3结果

如果一切顺利，您应该能够使用候选人您的槽位MmapprData对象或通过查看输出文件夹中的文件:

(候选人(mmapprData)的18美元,n = 2)

## GRanges对象有2个范围和30个元数据列:## seqnames ranges strand |等位基因后果##    |   ## 18:5744_A/C 18 5744 * | C missense_variant ## 18:5736_T/C 18 5736 * | C missense_variant ## IMPACT SYMBOL基因Feature_type ##     ## 18:5744_A/C MODERATE ENSDARG00000024771 ENSDARG00000024771转录## 18:5736_T/C MODERATE ENSDARG00000024771 ENSDARG00000024771转录## Feature BIOTYPE EXON INTRON ##    ## 18:5744_A/C ENSDART00000033574 protein_coding 6/9     ## 18:5744_A/C  940 880 ## 18:5736_T/C   932 872 ## Protein_position Amino_acids Codons Existing_variation ##     ## 18:5744_A/C 294 S/R Agc/Cgc  ## 18:5736_T/C 291 L/P cTa  ## DISTANCE STRAND FLAGSSYMBOL_SOURCE HGNC_ID ##      ## 18:5744_A/C  1    ## 18:5736_T/C    ## SOURCE FREQS CLIN_SIG SOMATIC PHENO ##      ## 18:5744_A/C slc24a5.gff.gz     ## 18:5736_T/C slc24a5.gff.gz     ## slc24a5.gff.gz ## 18:5741e -05 ## 18:574_a /C  4.63563e-05 ## ------- ##Seqinfo:基因组中的1个序列

outputTsv <- file.path(outputFolder(param(mmapprData))， '18.tsv') cat(paste(system2('head'， outputTsv, stdout = TRUE))， sep = '\n')

# #位置象征DensityScore等位基因影响后果胺基酸功能# # 5744 ENSDARG00000024771温和missense_variant 4.69570920024265 e-05 C S / R ENSDART00000033574 # # 5736 ENSDARG00000024771温和missense_variant 4.63563337842703 e-05 C L / P ENSDART00000033574 # # 5706 ENSDARG00000024771温和missense_variant 4.3464007193086 e-05 C I / T ENSDART00000033574 # # 4117 ENSDARG00000024771温和missense_variant e-05 C V / 4.30949569470718 L ENSDART00000033574 ENSDARG00000024771 # # 5494HIGH stop_gained 3.26854376246832e-05 G Y/* ENSDARG00000024771 MODERATE missense_variant 3.17179280841815e-05 C C/S ensdarg000000243574 ## 7197 ENSDARG00000024771 MODERATE missense_variant 2.54610321775923e-05 C M/I ensdarg00000033574 ## 7244 ENSDARG00000024771 MODERATE missense_variant 2.07444414554597e-05 C V/A ensdarg00000033574 ## 7336 ENSDARG00000024771 MODERATE missense_variant 1.15600976934635e-05 C I/L ensdar00000033574

3．1配置VEPFlags对象

mappr2使用ensemblVEPBioconductor包预测峰值区域变异的影响。要自定义此流程，您需要配置一个VEPFlags对象。看看ensemble的脚本选项网站．您可以配置VEPFlags这样的对象:

库(ensemble blvep, quiet = TRUE)

## ##附加包:'MatrixGenerics'

以下对象从'package:matrixStats'中被屏蔽:## ## colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse， ## colCounts, colCummaxs, colCummins, colCumprods, colCumsums， ## colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs， ## colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats， ## colProds， ## colQuantiles, colRanges, colRanks, colSdDiffs, colSds， ## colSums2, colTabulates, colVarDiffs, colWeightedMads， ## colWeightedMeans, colweighteddans， ## colWeightedVars, rowAlls, rowAnyNAs, rowAnys，rowAvgsPerColSet， ## rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods， ## rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps， ## rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins， ## rowOrderStats, rowProds, rowQuantiles, rowRanges, rowwranks， ## rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars， ## rowWeightedMads, rowWeightedMeans, rowWeightedMedians， ## rowWeightedVars

##欢迎访问Bioconductor ## ## Vignettes包含介绍性材料;使用## 'browseVignettes()'查看。要引用Bioconductor，请参阅## 'citation("Biobase")'，并查看软件包'citation("pkgname")'。

## ##附加包:“Biobase”

以下对象从“package:MatrixGenerics”被屏蔽:## ## rowMedians

以下对象从'package:matrixStats'中被屏蔽:## ## anyMissing, rowMedians

## ##附加包:'VariantAnnotation'

以下对象从'package:base': ## ## tabulate被屏蔽

## ##附加包:' ensemble blvep '

以下对象从'package:Biobase': ## ## cache中被屏蔽

以下对象从'package:BiocStyle': ## ## output被屏蔽

vepFlags <- vepFlags (flags = list(### DEFAULT SETTINGS format = 'vcf'， # <——this is necessary vcf = FALSE， # <——以及this species = 'danio_rerio'， database = FALSE, cache = TRUE, filter_common = TRUE, coding_only = TRUE假设RNA-seq数据###你可能会发现这些有趣的:# everything = TRUE #支持许多可选分析，如Polyphen和SIFT # per_gene = TRUE #将只输出每个基因最严重的变体# pick = TRUE #将只输出每个变体的一个结果))

3.2BiocParallel配置

MMAPPR2只是使用默认值bpparam注册。您可以对此进行更改(例如，ifBiocParallel没有正确地工作)与BiocParallel:注册命令。例如:

library(BiocParallel, quiet =TRUE) register(SerialParam()) register(MulticoreParam(progressbar=TRUE)) registered()

## $MulticoreParam ## class: MulticoreParam ## bpisup: FALSE;bpnworkers: 4;bptasks: 2147483647;bpjobname: BPJOB ## bplog: FALSE;bpthreshold:信息;## bpRNGseed:;bptimeout: NA;## bpexportglobals: TRUE;bpexportvariables:错误;bpforceGC:真实; bpfallback: TRUE ## bplogdir: NA ## bpresultdir: NA ## cluster type: FORK ## ## $SerialParam ## class: SerialParam ## bpisup: FALSE; bpnworkers: 1; bptasks: 0; bpjobname: BPJOB ## bplog: FALSE; bpthreshold: INFO; bpstopOnError: TRUE ## bpRNGseed: ; bptimeout: NA; bpprogressbar: FALSE ## bpexportglobals: FALSE; bpexportvariables: FALSE; bpforceGC: FALSE; bpfallback: FALSE ## bplogdir: NA ## bpresultdir: NA ## ## $SnowParam ## class: SnowParam ## bpisup: FALSE; bpnworkers: 4; bptasks: 0; bpjobname: BPJOB ## bplog: FALSE; bpthreshold: INFO; bpstopOnError: TRUE ## bpRNGseed: ; bptimeout: NA; bpprogressbar: FALSE ## bpexportglobals: TRUE; bpexportvariables: TRUE; bpforceGC: FALSE; bpfallback: TRUE ## bplogdir: NA ## bpresultdir: NA ## cluster type: SOCK

最后注册的参数成为默认值。

3．3参考基因组

变量调用步骤需要BiocStyle:: Biocpkg(“gmapR”)GmapGenome的文件自动生成refFasta参数。如果出于某种原因你想生成自己的，过程是这样的:

refGenome <- gmapR::GmapGenome(goldfasta ()， name='slc24a5'， create=TRUE)

3.4全基因组测序(WGS)

MMAPPR2和它的前身一样，是为RNA-Seq数据而设计和测试的。然而，工作中的原则仍然应该适用于WGS数据。