本教程的目标

完成本小插图后，您将能够:
1.读一个ChIP-seq实验R
2.扩展读取并保存数据(稍后讨论细节和相关性)
3.生成.bedGraph文件
4.可视化ChIP-seq数据R
5.对ChIP-seq峰值进行基本分析
6.生成围绕一组注释的基因组位点的ChIP-seq富集的平均剖面和热图
在附录部分中，我们将展示如何下载、预处理和评估.fastq文件的质量。

数据预处理

ChIP-seq分析工作流的第一部分是读预处理。在这里，我们将不关注这些最初的步骤，我们将概述它们并在附录小插图的一部分。下面简要介绍预处理中的三个主要步骤。

数据包

由上述步骤生成的数据文件被放置在名为EpigeneticsCSAMA，我们现在加载。(注意，在本课程中使用这样的数据包是为了方便，但通常情况下，您不会以这种方式打包中间数据。)

库(EpigeneticsCSAMA) dataDirectory =系统。文件(“bedfiles”,包=“EpigeneticsCSAMA”)

的变量dataDirectory显示包含此小插图所需的数据对象的目录。

dataDirectory

## [1] "/home/msmith/Applications/R/R-devel/library/EpigeneticsCSAMA/bedfiles"

为了研究这些文件，可以使用文本编辑器或终端模拟器查看它们。

箱子的生成

我们将把基因组平铺成大小为200 bp的不重叠的容器。为此，我们需要小鼠基因组(组装)中染色体大小的信息mm9)．在数据包中，我们提供对象如果(链信息)，用于保存这些数据。对象创建所需的代码如果对象中的获得Si *对象为*mm9**的附录．

数据(si) si

## Seqinfo对象，21个序列来自mm9基因组:## seqnames seqlengthiscircular genome ## chr1 197195432 FALSE mm9 ## chr2 181748087 FALSE mm9 ## chr3 159599783 FALSE mm9 ## chr4 155630120 FALSE mm9 ## chr5 152537259 FALSE mm9 ## ... ... ... ...## chr17 95272651 FALSE mm9 ## chr18 90772031 FALSE mm9 ## chr19 61342430 FALSE mm9 ## chrX 166650296 FALSE mm9 ## chrY 15902555 FALSE mm9

接下来，我们使用tileGenome函数从GenomicRanges包来生成农庄对象，其间隔以200 bp大小的瓦片(箱)覆盖基因组。

binsize = 200 bins = tileGenome(si['chr6']， tilewidth=binsize, cut.last.tile.in.chrom=TRUE) bins

##具有747586范围和0个元数据列的GRanges对象:## seqnames ranges strand ##    ## [1] chr6 [1,200] * ## [2] chr6 [201,400] * ## [3] chr6 [401, 600] * ## [4] chr6 [601, 800] * ## [5] chr6 [801,1000] * ## ... ... ... ...## [747582] chr6 [149516201, 149516400] * ## [747583] chr6 [149516401, 149516600] * ## [747584] chr6 [149516601, 149516800] * ## [747585] chr6 [149516801, 149517000] * ## [747586] chr6 [149517001, 149517037] * ## ------- ## seqinfo:来自mm9基因组的1个序列

装箱

我们现在计算每个bin中有多少读取。为此，我们定义了函数BinChIPseq．它需要两个参数，读取而且垃圾箱这是农庄对象。

BinChIPseq =函数(读取，箱子){mcols(箱子)$score = countOverlaps(箱子，读取)返回(箱子)}

现在我们把它应用到对象上输入，rep1而且rep2．我们获得input.200bins，rep1.200bins而且rep2.200bins，分别是农庄对象，包含输入和ChIP-seq实验的二进制读覆盖。

input.200bins= BinChIPseq( input, bins ) rep1.200bins = BinChIPseq( rep1, bins ) rep2.200bins = BinChIPseq( rep2, bins ) rep1.200bins

##具有747586范围和1个元数据列的GRanges对象:## seqnames ranges strand | score ##    |  ## [1] chr6 [1,200] * | 0 ## [2] chr6 [201,400] * | 0 ## [3] chr6 [401, 600] * | 0 ## [4] chr6 [601, 800] * | 0 ## [5] chr6 [801,1000] * | 0 ## ... ... ... ... . ...## [747582] chr6 [149516201, 149516400] * | 0 ## [747583] chr6 [149516401, 149516600] * | 0 ## [747584] chr6 [149516601, 149516800] * | 0 ## [747585] chr6 [149516801, 149517000] * | 0 ## [747586] chr6 [149517001, 149517037] * | 0 ## ------- ## seqinfo: 1个来自mm9基因组的序列

我们可以画出1000个箱子的覆盖率，从20万个箱子开始。

Plot (200000:201000, rep1.200bins$score[200000:201000]， xlab="chr6"， ylab="每箱计数")

下面，我们将看到更复杂的绘制覆盖率的方法。

导出二进制数据

在分析的这个步骤中，数据可以被可视化和共享了。共享ChIP-seq数据的最常用方法之一是生成.wig、. binwig或. bedgraph文件。它们是存储和大小有效的文件，保存了基因组中有关信号的信息。此外，这些文件可以在IGV和IGB等基因组浏览器中可视化。我们将展示如何将已打包的数据导出为. bedgraph文件。

出口(输入。200箱,= ' input_chr6监狱。bedGraph'， format = "bedGraph") export(rep1.200bins, con='H3K27ac_rep1_chr6. exe ')bedGraph'， format = "bedGraph") export(rep2.200bins, con='H3K27ac_rep2_chr6. bin ')， format = "bedGraph")

在下一节中，我们将看到如何使用可视化床图文件R．

峰的识别

正如我们在本教程的前一节中看到的，H3K27ac标记显示了明确的峰值。在这种情况下，苹果电脑是最常用的寻峰软件之一。ChIP-seq峰值通话也可以在R与BayesPeak包中。但是，我们在这里只讨论最常用的方法和用法苹果电脑．我们跑苹果电脑并在数据包中提供结果。中可以找到获取峰值所需的代码附录小插图。

峰的基本分析R

接下来，我们从数据包中导入隔离峰值的.bed文件。

峰值。rep1= import.bed(file.path(dataDirectory,'Rep1_peaks_ucsc_chr6.bed')) peaks.rep2 = import.bed(file.path(dataDirectory,'Rep2_peaks_ucsc_chr6.bed'))

分析已识别的峰值的第一步是简单地将它们与ChIP-seq和输入轨道一起显示在浏览器中。为此，我们使用AnnotationTrack函数。我们将峰值显示为蓝色的方框。

peaks1。track = AnnotationTrack(峰值。rep1，genome="mm9", name='Peaks Rep. 1', chromosome='chr6', shape='box',fill='blue3',size=2) peaks2.track = AnnotationTrack(peaks.rep2, genome="mm9", name='Peaks Rep. 2', chromosome='chr6', shape='box',fill='blue3',size=2)

我们想象Nanog轨迹。

plotTracks (c(输入。跟踪rep1。追踪,peaks1。追踪,rep2。追踪,peaks2。track, bm, AT)， from=122630000, to=122700000, transcriptAnnotation="symbol"， window="auto"， type="histogram"， cex.title=0.7, fontsize=10)

我们可以看到苹果电脑成功识别出H3K27ac信号高的区域。我们看到两种生物复制都很一致，然而，在某些情况下，只有一个样本才称为峰值。在下一节中，我们将分析如何经常看到峰之间的重叠和孤立的可重复的峰。

维恩图

我们首先找到两个重复的峰值集之间的重叠。

ovlp = findOverlaps(峰值。rep1峰值。Rep2) ovlp

##命中2528次，0元数据列的对象:## queryHits subjectHits ##   ## [1] 1 1 ## [2] 2 2 ## [3] 3 3 ## [4] 4 4 ## [5] 5 5 ## ... ... ...## [2524] 3025 3025 ## [2525] 3026 3026 ## [2526] 3029 3027 ## [2527] 3030 3028 ## [2528] 3031 3030 ## ------- ## queryLength: 3032 / subjectLength: 3032

如果一个复制中的峰值与另一个复制中的多个峰值重叠，则它将在中出现多次ovlp．为了查看有多少峰与另一个副本中的某物重叠，我们计算的两列中的每一列中的唯一峰的数量ovlp取这两个数中较小的数作为共有峰数。

ov = min(长度(唯一的(queryHits(ovlp)))，长度(唯一的(subjectHits(ovlp))))

我们把它画成维恩图，用draw.pairwise.venn函数从维恩图包中。

库(文氏图)

##加载所需的包

draw.pairwise.venn( area1=length(peaks.rep1), area2=length(peaks.rep2), cross.area=ov, category=c("rep1", "rep2"), fill=c("steelblue", "blue3"), cat.cex=0.7)

# #(多边形[GRID.polygon。100年],[GRID.polygon多边形。101年],[GRID.polygon多边形。102年],[GRID.polygon多边形。103年],[GRID.text文本。104年],[GRID.text文本。105年],[GRID.text文本。106年],[GRID.text文本。107),文本[GRID.text.108])

我们将只关注在两个重复中确定的峰值(以下称为富集区)。丰富的区域用绿色表示。

充实。区域=减少(subsetByOverlaps, list(峰值。rep1峰值。rep2)) enr.reg.track = AnnotationTrack(enriched.regions, genome="mm9", name='Enriched regions', chromosome='chr6', shape='box',fill='green3',size=2) plotTracks(c(input.track, rep1.track, peaks1.track, rep2.track, peaks2.track, enr.reg.track, bm, AT), from=122630000, to=122700000, transcriptAnnotation="symbol", window="auto", type="histogram", cex.title=0.5, fontsize=10 )

H3K27ac峰重叠启动子的分离

ChIP序列分析的问题之一是ChIP富集区域与选定的特征类型(如启动子或富集其他修饰的区域)重叠的范围。为此，本文分析了ChIP-seq信号峰值与感兴趣特征之间的重叠。

我们通过测试有多少H3K27ac富集区域重叠启动子区域来举例说明这种分析。

数据(egs)头(egs)

# # 1 # # ensembl_gene_id external_gene_id chromosome_name start_position ENSMUSG00000030270 Cpne9 6 113232301 # # 2 ENSMUSG00000001632 Brpf1 6 113257175 # # 3 ENSMUSG00000030271 Ogg1 6 113276966 # # 4 ENSMUSG00000030272 Camk1 6 113284118 # # 5 ENSMUSG00000048930 Tada3 6 113316019 # # 6 ENSMUSG00000079426 Arpc4 6 113328107 # # end_position链# # 1 113328107 1 # # 2 113274851 1 # # 3 113285062 1 # # 4 113285062 1 # # 5 113327877 1 # # 6 113327877 1

egs$TSS = ifelse(egs$strand == "1"， egs$start_position, egs$end_position) head(egs)

# # 1 # # ensembl_gene_id external_gene_id chromosome_name start_position ENSMUSG00000030270 Cpne9 6 113232301 # # 2 ENSMUSG00000001632 Brpf1 6 113257175 # # 3 ENSMUSG00000030271 Ogg1 6 113276966 # # 4 ENSMUSG00000030272 Camk1 6 113284118 # # 5 ENSMUSG00000048930 Tada3 6 113316019 # # 6 ENSMUSG00000079426 Arpc4 6 113328107 # # end_position链TSS # # 1 113328107 1 113232301 # # 2 113232301 1 113232301 # # 3 113285062 113276966 # # 4 113276966 1 113285062 # # 5 113327877 1 113327877 # # 6 113327877 1113328107

promoter_regions = GRanges(seqnames = Rle(paste0('chr'， egs$chromosome_name))， ranges = IRanges(start = egs$TSS - 200, end = egs$TSS + 200)， strand = Rle(rep("*"， nrow(egs)))， gene = egs$external_gene_id) promoter_regions

|基因##    |  ## [1] chr6 [113232101, 113232501] * | Cpne9 ## [2] chr6 [113256975, 113257375] * | Brpf1 ## [3] chr6 [113276766, 113277166] * | Ogg1 ## [4] chr6 [113293778, 113294178] * | Camk1 ## [5] chr6 [113327677, 113328077] * | Tada3 ## ... ... ... ... . ...## [1969] chr6 [30119537,30119937] * | Mir183 ## [1970] chr6 [93743566, 93743966] * | U6 ## [1971] chr6 [116321886, 116322286] * | SNORA17 ## [1972] chr6 [76150072,76150472] * | U1 ## [1973] chr6 [106951246, 106951646] * | U6 ## ------- # seqinfo: 1个来自未知基因组的序列;没有seqlengths

ovlp2 = finoverlaps(丰富的。区域，promoter_regions) cat(sprintf("%d的%d启动子被一个丰富的区域重叠。"，length(unique(subjectHits(ovlp2)))， length(promoter_regions)))))

1973启动子中的## 634被富集区域重叠。

ovlp2b = findOverlaps(promoter_regions，。区域)cat(sprintf(“%d的%d富集区域重叠启动子。”，长度(唯一的(subjectHits(ovlp2b)))，长度(富集。地区)))

2301个富集区域中的546个与一个启动子重叠。

Promotor_total_length = sum(width(reduce(promoter_regions)

## [1] 766386

Promotor_fraction_of_chromosome_6 = promotor_total_length / seqlength (si)["chr6"]

binom。test(长度(唯一的(subjectHits(ovlp2b))))，长度(充实。(区域)，promotor_fraction_of_chromosome_6)

## ##精确的二项式测试## ##数据:长度(唯一的(subjectHits(ovlp2b)))和长度(.regions) ##成功次数= 546，试验次数= 2301,p值< ## 0.00000000000000022 ##替代假设:真正成功的概率不等于0.005125744 ## 95%置信区间:## 0.2200317 0.2552159 ##样本估计:##成功的概率## 0.2372881

pos.TSS = egs[unique(queryHits(finoverllaps (promoter_regions，。)，] post . tss [1:3，]

## ENSMUSG00000030271 Ogg1 6 113276966 ## 4 ENSMUSG00000030272 Camk1 6 113284118 ## end_position strand TSS ## 2 113274851 1 113257175 ## 3 113285062 1 113276966 ## 4 113293978 -1 113293978

H3K27ac基因启动子的分布分析

在这一部分的分析中，我们展示了如何生成显示ChIP-seq信号在某些基因组位置周围分布的图，这里是一组启动子区域。这些包括热图表示和H3K27ac信号在启动子重叠的H3K27ac峰的平均剖面苹果电脑．这些是ChIP-seq实验中最常执行的分析步骤之一。

在前一节中，我们已经确定了重叠H3K27ac峰的启动子pos.TSS对象)。为了全面了解H3K27ac在相应TSS周围的分布情况，我们将分析区域扩展到\ 1000下午(\ \)bp围绕TSS。我们将每个2000 bp的区域分成20个长度为100 bp的箱子，并对箱子进行排序，“+”链上的基因的位置增加，“-”链上的基因的位置降低。

接下来，我们对启动子区域进行连续的100bp瓷砖拼接。对于每个区域，我们根据基因取向对瓦片进行排序。我们在每个启动子区域创建20个磁贴。

tiles = sapply(1:nrow(pos.TSS)，函数(i) if(pos.TSS$strand[i] == "1") pos.TSS$TSS[i] + seq(-1000, 900，长度。out=20) else post .TSS$TSS[i] + seq(900， -1000，长度。out=20)) tiles = GRanges(tilename = paste(rep(post . tss $ensembl_gene_id, each=20)， 1:20, sep="_")， seqnames = Rle(rep(paste0('chr'， post . tss $chromosome_name)， each=20))， ranges = IRanges(start = as.vector(tiles)， width = 100)， strand = Rle(rep("*"， length(as.vector(tiles))))， seqinfo=si) tiles

##具有12680范围和1个元数据列的GRanges对象:## seqnames ranges strand | tilename ##    |  ## [1] chr6 [113256175, 113256274] * | ENSMUSG00000001632_1 ## [2] chr6 [113256275, 113256374] * | ENSMUSG00000001632_2 ## [3] chr6 [113256375, 113256474] * | ENSMUSG00000001632_3 ## [4] chr6 [113256475, 113256574] * | ENSMUSG00000001632_4 ## [5] chr6 [113256575, 113256674] * | ENSMUSG00000001632_5 ## ... ... ... ... . ...## [12676] chr6 [30119137, 30119236] * | ENSMUSG00000065619_16 ## [12677] chr6 [30118937, 30119036] * | ENSMUSG00000065619_18 ## [12679] chr6 [30118837, 30118936] * | ENSMUSG00000065619_19 ## [12680] chr6 [30118737, 30118836] * | ENSMUSG00000065619_20 ## ------- ## seqinfo:来自mm9基因组的21个序列

接下来，我们计算映射到每个tile的读取数。得到的向量H3K27ac.p下一个用于创建矩阵(H3K27ac.p.matrix)，其中每一行都是h3k27ac富集启动子。每个柱对应一个连续的100bp瓦，在TSS周围的2000bp区域重叠一个H3K27ac峰值。由于我们将每个启动子区域划分为21个瓦，我们得到了一个21列634行的矩阵(H3K27ac峰重叠启动子的数量)。

H3K27ac.p= countOverlaps( tiles, rep1) + countOverlaps( tiles, rep2 ) H3K27ac.p.matrix = matrix( H3K27ac.p, nrow=nrow(pos.TSS), ncol=20, byrow=TRUE )

最后，我们将结果绘制为热图，并将所有包含的启动子的每个瓦的平均值绘制为图。

颜色= colorRampPalette (c(“白”、“红”、“灰色”,“黑色”))(100)布局(垫=矩阵(c(1 2 0 3), 2, 2),宽度= c(2 2 2),高度= c(0.5, 5, 0.5, 5),真的)标准(mar = c(4 4、1.5、1)图像(seq(0,最大值(H3K27ac.p.matrix) length.out = 100), 1,矩阵(seq(0,最大值(H3K27ac.p.matrix) length.out = 100), 100年,1),坳=颜色,xlab =距离TSS, ylab = ",主要=数量的读取,yaxt = ' n ', lwd = 3,轴= TRUE)框(col =‘黑’,lwd = 2)图像(x = seq(-1000、1000、length.out = 20), y = 1: nrow (H3K27ac.p.matrix),z=t(H3K27ac.p.matrix[order(rowsum (H3K27ac.p.matrix))，])， col=colors, xlab='距离TSS (bp)'， ylab='Promoters'， lwd=2) box(col='black'， lwd=2) abline(v=0, lwd=1, col='gray') plot(x=seq(- 1000,1000, length.out=20)， y=colMeans(H3K27ac.p.matrix)， ty='b'， pch=19, col='red4'，lwd=2, ylab='Mean tag count'， xlab='距离TSS (bp)') abline(h=seq(1100,by=5)， v=seq(- 1000,1000, length.out=20)， lwd=0.25, col='gray') box(col='black'， lwd=2)

我们观察到在TSS之后H3K27ac修饰的强烈富集，在TSS的上游区域有一个较弱的H3K27ac修饰峰。

从欧洲核苷酸档案获取数据

欧洲核苷酸档案(http://www.ebi.ac.uk/ena)提供多种类型的原始测序数据、序列组装信息和功能注释。我们下载了小鼠胚胎干细胞(mES细胞)中H3K27ac组蛋白修饰的ChIP-seq实验对应的数据以及该研究的输入对照样本组蛋白H3K27ac分离活性增强子和稳定增强子并预测发育状态由Creyghton等．

wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR066/SRR066787/SRR066787.fastq.gz。wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR066/SRR066766/SRR066766.fastq.gz。wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR066/SRR066767/SRR066767.fastq.gz。

阅读质量

读取质量是所有序列读取分析的第一步。这个包ShortRead提供了一个函数，将从ENA数据库下载的.fastq文件作为输入。我们首先用fastq文件名生成一个向量。

FLS =列表。文件(dataDirectory,”。fastq$"， full=TRUE) names(fls) = sub("。Fastq "， ""， basename(fls))

我们读取这些文件并应用qas函数评估每个文件中的读取质量。最后，我们生成一个超文本标记语言质量报告。

library(ShortRead) qas = lapply(seq_along(fls)， function(i, fls) qa(readFastq(fls[i])， names(fls)[i])， fls) qa = do。call(rbind, qas) rpt = report(qa,dest = 'QA_report.html')

外部文件准备

下一步是将reads与mm9小鼠基因组组装对齐。这是用Bowtie2工具。生成的.sam文件接下来转换为.bam文件，并使用samtools．去除PCR重复序列。BAM文件接下来转换为床文件。为了与其他工具保持一致，在数据预处理的最后一步中，我们在使用的染色体名称前添加了一个' chr '前缀awk．

gunzip SRR066767.fastq.gz gunzip SRR066767.fastq.gz

对齐

bowtie2 -p 8 -q NCBIM37.67 SRR066787。fastq -S ES_input。sam bowtie2 -p 8 -q NCBIM37.67 SRR066766。fastq -S H3K27ac_rep1。sam bowtie2 -p 8 -q NCBIM37.67 SRR066767。fastq -S H3K27ac_rep2.sam

只保留最佳对齐

samtools view -bS -q 40 ES_input. xmlsam > ES_input_bestAlignment。bam samtools view -bS -q 40 H3K27ac_rep1。sam > H3K27ac_rep1_bestAlignment。bam samtools view -bS -q 40 H3K27ac_rep2。sam > h3k27ac_rep2_bestalignments .bam

PCR重复去除

samtools rmdup -s ES_input_bestAlignment。bam ES_input_filtered。bam samtools rmdup -s H3K27ac_rep1_bestAlignment。bam H3K27ac_rep1_filtered。bam samtools rmdup -s H3K27ac_rep2_bestAlignment。bam H3K27ac_rep2_filtered.bam

将读取转换为.bed格式

bedtools bamtobed -i ES_input_filteredbam > ES_input_filtered. .bed bedtools bamtobed -i H3K27ac_rep1_filtered。bam > H3K27ac_rep1_filtered。bed bedtools bamtobed -i H3K27ac_rep2_filtered。bam > H3K27ac_rep2_filtered.bed

额外的准备工作

awk '$0="chr"$0' ES_input_filtered. .bed > ES_input_filtered_ucsc. bedbed awk '$0="chr"$0' H3K27ac_rep1_filtered。bed > H3K27ac_rep1_filtered_ucsc。bed awk '$0="chr"$0' H3K27ac_rep2_filtered。床> H3K27ac_rep2_filtered_ucsc.bed

最后，为了本实验室的目的，我们仅分离了一条染色体(chr6)的数据。

awk '{if($1=="chr6") print $0}' ES_input_filtered_ucsc. awk '{if($1=="chr6") print $0}'bed > ES_input_filtered_ucsc_chr6。' H3K27ac_rep1_filtered_ucsc '{if($1=="chr6") print $0}'bed > H3K27ac_rep1_filtered_ucsc_chr6。bed awk '{if($1=="chr6") print $0}' H3K27ac_rep2_filtered_ucsc。床> H3K27ac_rep2_filtered_ucsc_chr6.bed

library(BSgenome.Mmusculus.UCSC.mm9) genome = BSgenome.Mmusculus.UCSC. BSgenome.Mmusculus.UCSC.mm9)mm9 si = seqinfo(基因组)si = si[paste0('chr'， c(1:19， 'X'， 'Y'))]

library(biomaRt) mart = useMart(biomaRt =" ENSEMBL_MART_ENSEMBL"， dataset =" mmusculus_gene_ensembl"， host="may2012.archive.ensembl.org") fm = Gviz:::.getBMFeatureMap() fm["symbol"] =" external_gene_id"

bm = BiomartGeneRegionTrack(染色体='chr6'，基因组="mm9"， start=122530000, end = 122900000, biomart=mart,filter=list("with_ox_refseq_mrna"=TRUE)， size=4, name="RefSeq"， utr5="red3"， utr3="red3"， protein_coding="black"， col.line=NULL, cex=7, collapseTranscripts="longest"， featureMap=fm)

macs14 -t H3K27ac_rep1_filtered。-c ES_input_filtered_ucsc.使用实例bed -f bed -g mm——nommodel -n Rep1 macs14 -t H3K27ac_rep2_filtered。-c ES_input_filtered_ucsc.使用实例bed -f bed -g mm——nommodel -n Rep2 awk '$0="chr"$0' Rep1_peaks。bed > Rep1_peaks_ucsc。bed awk '$0="chr"$0' Rep2_peaks。bed > Rep2_peaks_ucsc。{if($1=="chr6") print $0}' Rep1_peaks_ucsc. bed awk '{if($1=="chr6")bed > Rep1_peaks_ucsc_chr6。* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *。Rep2_peaks_ucsc_chr6.bed >

listAttributes(mart)[1:3，] ds = useDataset('mmusculus_gene_ensembl'， mart=mart) chroms = 6 egs = getBM(attributes =c ('ensembl_gene_id'，'external_gene_id'， 'chromosome_name'，'start_position'， 'end_position'，'strand')， filters='chromosome_name'， values=chroms, mart=ds)