admin管理员组文章数量:1037775
在R语言中的 ATACseq 数据分析全流程实战（二）
续上一篇实战：在R语言中的 ATACseq 数据分析全流程实战（一）
（续）
比对结果：运行了一个小时多一点
代码语言：javascript代码运行次数：0运行复制//================================ Summary =================================\\
|| ||
|| Total fragments : 56,598,621 ||
|| Mapped : 49,993,415 (88.3%) ||
|| Uniquely mapped : 49,993,415 ||
|| Multi-mapping : 0 ||
|| ||
|| Unmapped : 6,605,206 ||
|| ||
|| Properly paired : 46,024,907 ||
|| Not properly paired : 3,968,508 ||
|| Singleton : 1,592,254 ||
|| Chimeric : 312,410 ||
|| Unexpected strandness : 26,237 ||
|| Unexpected fragment length : 222,426 ||
|| Unexpected read order : 1,815,181 ||
|| ||
|| Indels : 56,370 ||
|| ||
|| Running time : 69.4 minutes ||
|| ||
\\============================================================================//
ATAC_50K_2.bam
Total_fragments 56598621
Mapped_fragments 49993415
Uniquely_mapped_fragments 49993415
Multi_mapping_fragments 0
Unmapped_fragments 6605206
Properly_paired_fragments 46024907
Singleton_fragments 1592254
More_than_one_chr_fragments 312410
Unexpected_strandness_fragments 26237
Unexpected_template_length 222426
Inversed_mapping 1815181
Indels 56370
4.4 sam转bam并构建索引
首先，也是需要创建索引：
代码语言：javascript代码运行次数：0运行复制# 对bam进行排序，并构建索引
library(Rsamtools)
outBAM <- "ATAC_50K_2.bam"
sortedBAM <- "ATAC_50K_2_sorted.bam"
sortBam(outBAM, gsub("\\.bam","",sortedBAM) )
indexBam(sortedBAM)
在ATAC-seq中，可以检查比对上的reads在染色体上的分布情况，使用idxstatsBam()函数来检查每个染色体上比对上的reads数量。
ATAC-seq以其在线粒体染色体上的高信号而闻名，我们可以在这里检查这一点。
代码语言：javascript代码运行次数：0运行复制library(Rsamtools)
mappedReads <- idxstatsBam(sortedBAM)
mappedReads
# seqnames seqlength mapped unmapped
# 1 chr1 249250621 3712101 0
# 2 chr2 243199373 2880468 0
# 3 chr3 198022430 2575872 0
# 4 chr4 191154276 2289209 0
# 5 chr5 180915260 2652191 0
# 6 chr6 171115067 2122969 0
# 7 chr7 159138663 1827870 0
# 8 chr8 146364022 1887888 0
# 9 chr9 141213431 1284648 0
# 10 chr10 135534747 1614971 0
# 11 chr11 135006516 1712776 0
# 12 chr12 133851895 1555460 0
# 13 chr13 115169878 1024664 0
# 14 chr14 107349540 900443 0
# 15 chr15 102531392 728044 0
# 16 chr16 90354753 951870 0
# 17 chr17 81195210 894934 0
# 18 chr18 78077248 882604 0
# 19 chr19 59128983 626703 0
# 20 chr20 63025520 666760 0
# 21 chr21 48129895 369333 0
# 22 chrX 155270560 1472796 0
# 23 chrY 59373566 81610 0
# 24 chrM 16571 63678392 0
# 25 * 0 0 14802666
# 绘图展示
library(ggplot2)
ggplot(mappedReads, aes(seqnames, mapped, fill = seqnames)) +
geom_bar(stat = "identity") +
coord_flip()
在这个例子中，可以看到比对到线粒体基因组的比率很高的情况：
5.比对后处理
比对完了之后，我们就可以开始处理比对文件了。
5.1 Proper pairs
首先，可以查看 ATAC-seq 数据预期的片段长度分布，使用GenomicAlignments包读取新比对的数据。使用ScanBamParam()和scanBamFlag()函数来控制只读取 Proper pairs的 reads 到R中。
what：指定哪些信息被读取进来which：可以指定读取的染色体号与范围
################### 比对后处理
library(GenomicAlignments)
flags = scanBamFlag(isProperPair = TRUE)
myParam = ScanBamParam(flag = flags,
what = c("qname", "mapq", "isize"),
which = GRanges("chr20",IRanges(1, 63025520)))
myParam
atacReads <- readGAlignmentPairs(sortedBAM, param = myParam)
class(atacReads)
返回的结果是一个GAlignmentPairs对象。
5.2 MapQ打分
代码语言：javascript代码运行次数：0运行复制############## 提取Q打分并绘制分布图
# 使用first()和second()来访问GAlignments对象，以分别获取第一条或第二条读取的信息
read1 <- GenomicAlignments::first(atacReads)
read2 <- GenomicAlignments::second(atacReads)
read1[1, ]
read2[1, ]
# 提取 MapQ 打分
read1MapQ <- mcols(read1)$mapq
read2MapQ <- mcols(read2)$mapq
read1MapQ[1:2]
# 计算MapQ 打分 分布
read1MapQFreqs <- table(read1MapQ)
read2MapQFreqs <- table(read2MapQ)
read1MapQFreqs
read2MapQFreqs
# 绘制
library(ggplot2)
toPlot <- data.frame(MapQ = c(names(read1MapQFreqs), names(read2MapQFreqs)),
Frequency = c(read1MapQFreqs,read2MapQFreqs),
Read = c(rep("Read1", length(read1MapQFreqs)), rep("Read2",length(read2MapQFreqs)))
)
toPlot$MapQ <- factor(toPlot$MapQ, levels = unique(sort(as.numeric(toPlot$MapQ))))
head(toPlot)
ggplot(toPlot, aes(x = MapQ, y = Frequency, fill = MapQ)) +
geom_bar(stat = "identity") +
facet_grid(~Read)
比对Q值分布如下：
5.3 提取insert sizes
可以从GAlignments对象的elementMetadata()中检索插入片段的大小：
############### 提取插入片段大小
atacReads_read1 <- GenomicAlignments::first(atacReads)
insertSizes <- abs(elementMetadata(atacReads_read1)$isize)
head(insertSizes)
fragLenSizes <- table(insertSizes)
fragLenSizes[1:5]
library(ggplot2)
toPlot <- data.frame(InsertSize = as.numeric(names(fragLenSizes)),
Count = as.numeric(fragLenSizes))
fragLenPlot <- ggplot(toPlot, aes(x = InsertSize, y = Count)) +
geom_line() +
theme_bw()
fragLenPlot
ATAC-seq应该代表对应于无核小体、单核小体和多核小体部分的各种片段长度的混合。来看下获取的20号染色体上插入片段的长度分布：
对计数应用log2转换，以更清晰地展示核小体模式，并根据原文中的方法，对无核小体（小于100bp）、单核小体（180bp-247bp）和双核小体（315bp-437bp）区域进行注释：
代码语言：javascript代码运行次数：0运行复制fragLenPlot <- ggplot(toPlot, aes(x = InsertSize, y = Count)) +
geom_line() +
scale_y_continuous(trans = "log2") +
geom_vline(xintercept = c(100), colour = "darkgreen") + # 无核小体（小于100bp）
geom_vline(xintercept = c(180,247), colour = "red") + # 单核小体（180bp-247bp）
geom_vline(xintercept = c(315, 437), colour = "darkblue") + # 双核小体（315bp-437bp）
theme_bw()
fragLenPlot
结果如下：
可以与原文中的对比，这两者的结果非常一致：
;gid=article-figures&pid=figure-2-uid-1
此外，还可以通过使用插入片段的大小来将比对reads分为代表无核小体和有核小体占据的reads。
在这里，通过使用插入片段的大小来过滤read，创建代表无核小体、单核小体和双核小体的读取的BAM文件。生成的reads可以被写回到BAM文件中，以便在分析的其他部分使用，或者通过rtracklayer包中的函数在IGV等程序中进行可视化。
代码语言：javascript代码运行次数：0运行复制################ 创建代表无核小体、单核小体和双核小体的读取的BAM文件
atacReads_NucFree <- atacReads[insertSizes < 100, ]
atacReads_MonoNuc <- atacReads[insertSizes > 180 & insertSizes < 240, ]
atacReads_diNuc <- atacReads[insertSizes > 315 & insertSizes < 437, ]
nucFreeRegionBam <- gsub("\\.bam", "_nucFreeRegions\\.bam", sortedBAM)
monoNucBam <- gsub("\\.bam", "_monoNuc\\.bam", sortedBAM)
diNucBam <- gsub("\\.bam", "_diNuc\\.bam", sortedBAM)
library(rtracklayer)
export(atacReads_NucFree, nucFreeRegionBam, format = "bam")
export(atacReads_MonoNuc, monoNucBam, format = "bam")
export(atacReads_diNuc, diNucBam, format = "bam")
将生成的文件ATAC_50K_2_sorted_openRegionRPM.bw导入到IGV中：
6.ATACseqQC
ATACseqQC这个包可以计算ATACseq数据的两个非常重要的QC指标：PCR Bottleneck Coefficients (PBC1 and PBC2)。
与ChIPQC类似，ATACseqQC函数包含了一个工作流函数，该函数仅通过一个BAM文件路径参数就能获取大部分所需的质量控制（QC）信息。
代码语言：javascript代码运行次数：0运行复制############# ATACseqQC
# 安装
# BiocManager::install("ATACseqQC")
library(ATACseqQC)
ATACQC <- bamQC("ATAC_50K_2_sorted.bam")
# ATACQC object
names(ATACQC)
## [1] "totalQNAMEs" "duplicateRate"
## [3] "mitochondriaRate" "properPairRate"
## [5] "unmappedRate" "hasUnmappedMateRate"
## [7] "notPassingQualityControlsRate" "nonRedundantFraction"
## [9] "PCRbottleneckCoefficient_1" "PCRbottleneckCoefficient_2"
## [11] "MAPQ" "idxstats"
返回的ATACQC对象中包含了非常多的质控指标，如duplicateRate, non-redundant fraction, distribution of signal across chromosomes, mitochondrial fraction等，还有刚刚说的**PCRbottleneckCoefficient_1 和 PCRbottleneckCoefficient_2。
6.1 指标1：PCR Bottleneck Coefficient
PCR瓶颈系数（PCR Bottleneck Coefficient，PBC）是衡量文库复杂性的一项指标，用于识别在ATAC样本制备过程中可能发生的PCR偏差/过度扩增。 PCRbottleneckCoefficient_1的计算公式为PBC1 = N1/Nd，其中N1是基因组中只有一个唯一比read的位置数，Nd是至少有一个比read的位置数。 PCRbottleneckCoefficient_2的计算公式为PBC2 = N1/N2，其中N1是基因组中只有一个唯一比read的位置数，N2是只有两个唯一比read的位置数。
示例：有20个reads，其中16个reads是唯一比对read，4个不唯一对比（比对到两个位置，一个位置上有两个read），那么：
PBC1 = 16/18=0.889，PBC2=16/2=8
取值范围：
PBC1：<0.7表示严重瓶颈，0.7-0.9表示中度瓶颈，>0.9表示没有瓶颈
PBC2：<1表示严重瓶颈，1-3表示中度瓶颈，>3 表示没有瓶颈
代码语言：javascript代码运行次数：0运行复制############# ATACseqQC
library(ATACseqQC)
ATACQC <- bamQC("ATAC_50K_2_sorted.bam")
# ATACQC object
names(ATACQC)
## [1] "totalQNAMEs" "duplicateRate"
## [3] "mitochondriaRate" "properPairRate"
## [5] "unmappedRate" "hasUnmappedMateRate"
## [7] "notPassingQualityControlsRate" "nonRedundantFraction"
## [9] "PCRbottleneckCoefficient_1" "PCRbottleneckCoefficient_2"
## [11] "MAPQ" "idxstats"
ATACQC$PCRbottleneckCoefficient_1
# [1] 0.6652245
ATACQC$PCRbottleneckCoefficient_2
# [1] 3.413779
6.2 课后习题
题目：.html
答案：.html
代码：.R
（未完待续~）
本文参与 腾讯云自媒体同步曝光计划，分享自微信公众号。原始发表：2025-03-13，如有侵权请联系 cloudcommunity@tencent 删除索引数据分析对象函数数据在R语言中的 ATACseq 数据分析全流程实战（二）
续上一篇实战：在R语言中的 ATACseq 数据分析全流程实战（一）
（续）
比对结果：运行了一个小时多一点
代码语言：javascript代码运行次数：0运行复制//================================ Summary =================================\\
|| ||
|| Total fragments : 56,598,621 ||
|| Mapped : 49,993,415 (88.3%) ||
|| Uniquely mapped : 49,993,415 ||
|| Multi-mapping : 0 ||
|| ||
|| Unmapped : 6,605,206 ||
|| ||
|| Properly paired : 46,024,907 ||
|| Not properly paired : 3,968,508 ||
|| Singleton : 1,592,254 ||
|| Chimeric : 312,410 ||
|| Unexpected strandness : 26,237 ||
|| Unexpected fragment length : 222,426 ||
|| Unexpected read order : 1,815,181 ||
|| ||
|| Indels : 56,370 ||
|| ||
|| Running time : 69.4 minutes ||
|| ||
\\============================================================================//
ATAC_50K_2.bam
Total_fragments 56598621
Mapped_fragments 49993415
Uniquely_mapped_fragments 49993415
Multi_mapping_fragments 0
Unmapped_fragments 6605206
Properly_paired_fragments 46024907
Singleton_fragments 1592254
More_than_one_chr_fragments 312410
Unexpected_strandness_fragments 26237
Unexpected_template_length 222426
Inversed_mapping 1815181
Indels 56370
4.4 sam转bam并构建索引
首先，也是需要创建索引：
代码语言：javascript代码运行次数：0运行复制# 对bam进行排序，并构建索引
library(Rsamtools)
outBAM <- "ATAC_50K_2.bam"
sortedBAM <- "ATAC_50K_2_sorted.bam"
sortBam(outBAM, gsub("\\.bam","",sortedBAM) )
indexBam(sortedBAM)
在ATAC-seq中，可以检查比对上的reads在染色体上的分布情况，使用idxstatsBam()函数来检查每个染色体上比对上的reads数量。
ATAC-seq以其在线粒体染色体上的高信号而闻名，我们可以在这里检查这一点。
代码语言：javascript代码运行次数：0运行复制library(Rsamtools)
mappedReads <- idxstatsBam(sortedBAM)
mappedReads
# seqnames seqlength mapped unmapped
# 1 chr1 249250621 3712101 0
# 2 chr2 243199373 2880468 0
# 3 chr3 198022430 2575872 0
# 4 chr4 191154276 2289209 0
# 5 chr5 180915260 2652191 0
# 6 chr6 171115067 2122969 0
# 7 chr7 159138663 1827870 0
# 8 chr8 146364022 1887888 0
# 9 chr9 141213431 1284648 0
# 10 chr10 135534747 1614971 0
# 11 chr11 135006516 1712776 0
# 12 chr12 133851895 1555460 0
# 13 chr13 115169878 1024664 0
# 14 chr14 107349540 900443 0
# 15 chr15 102531392 728044 0
# 16 chr16 90354753 951870 0
# 17 chr17 81195210 894934 0
# 18 chr18 78077248 882604 0
# 19 chr19 59128983 626703 0
# 20 chr20 63025520 666760 0
# 21 chr21 48129895 369333 0
# 22 chrX 155270560 1472796 0
# 23 chrY 59373566 81610 0
# 24 chrM 16571 63678392 0
# 25 * 0 0 14802666
# 绘图展示
library(ggplot2)
ggplot(mappedReads, aes(seqnames, mapped, fill = seqnames)) +
geom_bar(stat = "identity") +
coord_flip()
在这个例子中，可以看到比对到线粒体基因组的比率很高的情况：
5.比对后处理
比对完了之后，我们就可以开始处理比对文件了。
5.1 Proper pairs
首先，可以查看 ATAC-seq 数据预期的片段长度分布，使用GenomicAlignments包读取新比对的数据。使用ScanBamParam()和scanBamFlag()函数来控制只读取 Proper pairs的 reads 到R中。
what：指定哪些信息被读取进来which：可以指定读取的染色体号与范围
################### 比对后处理
library(GenomicAlignments)
flags = scanBamFlag(isProperPair = TRUE)
myParam = ScanBamParam(flag = flags,
what = c("qname", "mapq", "isize"),
which = GRanges("chr20",IRanges(1, 63025520)))
myParam
atacReads <- readGAlignmentPairs(sortedBAM, param = myParam)
class(atacReads)
返回的结果是一个GAlignmentPairs对象。
5.2 MapQ打分
代码语言：javascript代码运行次数：0运行复制############## 提取Q打分并绘制分布图
# 使用first()和second()来访问GAlignments对象，以分别获取第一条或第二条读取的信息
read1 <- GenomicAlignments::first(atacReads)
read2 <- GenomicAlignments::second(atacReads)
read1[1, ]
read2[1, ]
# 提取 MapQ 打分
read1MapQ <- mcols(read1)$mapq
read2MapQ <- mcols(read2)$mapq
read1MapQ[1:2]
# 计算MapQ 打分 分布
read1MapQFreqs <- table(read1MapQ)
read2MapQFreqs <- table(read2MapQ)
read1MapQFreqs
read2MapQFreqs
# 绘制
library(ggplot2)
toPlot <- data.frame(MapQ = c(names(read1MapQFreqs), names(read2MapQFreqs)),
Frequency = c(read1MapQFreqs,read2MapQFreqs),
Read = c(rep("Read1", length(read1MapQFreqs)), rep("Read2",length(read2MapQFreqs)))
)
toPlot$MapQ <- factor(toPlot$MapQ, levels = unique(sort(as.numeric(toPlot$MapQ))))
head(toPlot)
ggplot(toPlot, aes(x = MapQ, y = Frequency, fill = MapQ)) +
geom_bar(stat = "identity") +
facet_grid(~Read)
比对Q值分布如下：
5.3 提取insert sizes
可以从GAlignments对象的elementMetadata()中检索插入片段的大小：
############### 提取插入片段大小
atacReads_read1 <- GenomicAlignments::first(atacReads)
insertSizes <- abs(elementMetadata(atacReads_read1)$isize)
head(insertSizes)
fragLenSizes <- table(insertSizes)
fragLenSizes[1:5]
library(ggplot2)
toPlot <- data.frame(InsertSize = as.numeric(names(fragLenSizes)),
Count = as.numeric(fragLenSizes))
fragLenPlot <- ggplot(toPlot, aes(x = InsertSize, y = Count)) +
geom_line() +
theme_bw()
fragLenPlot
ATAC-seq应该代表对应于无核小体、单核小体和多核小体部分的各种片段长度的混合。来看下获取的20号染色体上插入片段的长度分布：
对计数应用log2转换，以更清晰地展示核小体模式，并根据原文中的方法，对无核小体（小于100bp）、单核小体（180bp-247bp）和双核小体（315bp-437bp）区域进行注释：
代码语言：javascript代码运行次数：0运行复制fragLenPlot <- ggplot(toPlot, aes(x = InsertSize, y = Count)) +
geom_line() +
scale_y_continuous(trans = "log2") +
geom_vline(xintercept = c(100), colour = "darkgreen") + # 无核小体（小于100bp）
geom_vline(xintercept = c(180,247), colour = "red") + # 单核小体（180bp-247bp）
geom_vline(xintercept = c(315, 437), colour = "darkblue") + # 双核小体（315bp-437bp）
theme_bw()
fragLenPlot
结果如下：
可以与原文中的对比，这两者的结果非常一致：
;gid=article-figures&pid=figure-2-uid-1
此外，还可以通过使用插入片段的大小来将比对reads分为代表无核小体和有核小体占据的reads。
在这里，通过使用插入片段的大小来过滤read，创建代表无核小体、单核小体和双核小体的读取的BAM文件。生成的reads可以被写回到BAM文件中，以便在分析的其他部分使用，或者通过rtracklayer包中的函数在IGV等程序中进行可视化。
代码语言：javascript代码运行次数：0运行复制################ 创建代表无核小体、单核小体和双核小体的读取的BAM文件
atacReads_NucFree <- atacReads[insertSizes < 100, ]
atacReads_MonoNuc <- atacReads[insertSizes > 180 & insertSizes < 240, ]
atacReads_diNuc <- atacReads[insertSizes > 315 & insertSizes < 437, ]
nucFreeRegionBam <- gsub("\\.bam", "_nucFreeRegions\\.bam", sortedBAM)
monoNucBam <- gsub("\\.bam", "_monoNuc\\.bam", sortedBAM)
diNucBam <- gsub("\\.bam", "_diNuc\\.bam", sortedBAM)
library(rtracklayer)
export(atacReads_NucFree, nucFreeRegionBam, format = "bam")
export(atacReads_MonoNuc, monoNucBam, format = "bam")
export(atacReads_diNuc, diNucBam, format = "bam")
将生成的文件ATAC_50K_2_sorted_openRegionRPM.bw导入到IGV中：
6.ATACseqQC
ATACseqQC这个包可以计算ATACseq数据的两个非常重要的QC指标：PCR Bottleneck Coefficients (PBC1 and PBC2)。
与ChIPQC类似，ATACseqQC函数包含了一个工作流函数，该函数仅通过一个BAM文件路径参数就能获取大部分所需的质量控制（QC）信息。
代码语言：javascript代码运行次数：0运行复制############# ATACseqQC
# 安装
# BiocManager::install("ATACseqQC")
library(ATACseqQC)
ATACQC <- bamQC("ATAC_50K_2_sorted.bam")
# ATACQC object
names(ATACQC)
## [1] "totalQNAMEs" "duplicateRate"
## [3] "mitochondriaRate" "properPairRate"
## [5] "unmappedRate" "hasUnmappedMateRate"
## [7] "notPassingQualityControlsRate" "nonRedundantFraction"
## [9] "PCRbottleneckCoefficient_1" "PCRbottleneckCoefficient_2"
## [11] "MAPQ" "idxstats"
返回的ATACQC对象中包含了非常多的质控指标，如duplicateRate, non-redundant fraction, distribution of signal across chromosomes, mitochondrial fraction等，还有刚刚说的**PCRbottleneckCoefficient_1 和 PCRbottleneckCoefficient_2。
6.1 指标1：PCR Bottleneck Coefficient
PCR瓶颈系数（PCR Bottleneck Coefficient，PBC）是衡量文库复杂性的一项指标，用于识别在ATAC样本制备过程中可能发生的PCR偏差/过度扩增。 PCRbottleneckCoefficient_1的计算公式为PBC1 = N1/Nd，其中N1是基因组中只有一个唯一比read的位置数，Nd是至少有一个比read的位置数。 PCRbottleneckCoefficient_2的计算公式为PBC2 = N1/N2，其中N1是基因组中只有一个唯一比read的位置数，N2是只有两个唯一比read的位置数。
示例：有20个reads，其中16个reads是唯一比对read，4个不唯一对比（比对到两个位置，一个位置上有两个read），那么：
PBC1 = 16/18=0.889，PBC2=16/2=8
取值范围：
PBC1：<0.7表示严重瓶颈，0.7-0.9表示中度瓶颈，>0.9表示没有瓶颈
PBC2：<1表示严重瓶颈，1-3表示中度瓶颈，>3 表示没有瓶颈
代码语言：javascript代码运行次数：0运行复制############# ATACseqQC
library(ATACseqQC)
ATACQC <- bamQC("ATAC_50K_2_sorted.bam")
# ATACQC object
names(ATACQC)
## [1] "totalQNAMEs" "duplicateRate"
## [3] "mitochondriaRate" "properPairRate"
## [5] "unmappedRate" "hasUnmappedMateRate"
## [7] "notPassingQualityControlsRate" "nonRedundantFraction"
## [9] "PCRbottleneckCoefficient_1" "PCRbottleneckCoefficient_2"
## [11] "MAPQ" "idxstats"
ATACQC$PCRbottleneckCoefficient_1
# [1] 0.6652245
ATACQC$PCRbottleneckCoefficient_2
# [1] 3.413779
6.2 课后习题
题目：.html
答案：.html
代码：.R
（未完待续~）
本文参与 腾讯云自媒体同步曝光计划，分享自微信公众号。原始发表：2025-03-13，如有侵权请联系 cloudcommunity@tencent 删除索引数据分析对象函数数据本文标签： 在R语言中的 ATACseq 数据分析全流程实战（二）
版权声明：本文标题：在R语言中的 ATACseq 数据分析全流程实战（二） 内容由热心网友自发贡献，该文观点仅代表作者本人， 转载请联系作者并注明出处：http://it.en369.cn/jiaocheng/1748290116a2280652.html， 本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌抄袭侵权/违法违规的内容，一经查实，本站将立刻删除。
发表评论