V1

2023/01/27阅读：52主题：科技蓝

# ATAC-seq分析：差异分析（10）

## 1. 识别非冗余峰

``peaks <- dir("~/Downloads/ATAC_Workshop/ATAC_Data/ATAC_Peaks_forCounting/", pattern = "*.narrowPeak",    full.names = TRUE)myPeaks <- lapply(peaks, ChIPQC:::GetGRanges, simple = TRUE)allPeaksSet_nR <- reduce(unlist(GRangesList(myPeaks)))overlap <- list()for (i in 1:length(myPeaks)) {    overlap[[i]] <- allPeaksSet_nR %over% myPeaks[[i]]}overlapMatrix <- do.call(cbind, overlap)colnames(overlapMatrix) <- basename(peaks)mcols(allPeaksSet_nR) <- overlapMatrix``
``allPeaksSet_nR[1:2, ]``

``blklist <- import.bed("data/ENCFF547MET.bed.gz")nrToCount <- allPeaksSet_nR[!allPeaksSet_nR %over% blklist & !seqnames(allPeaksSet_nR) %in%    "chrM"]nrToCount``

## 2. 差异计数

``library(Rsubread)occurrences <- rowSums(as.data.frame(elementMetadata(nrToCount)))nrToCount <- nrToCount[occurrences >= 2, ]nrToCount``

``library(GenomicAlignments)bamsToCount <- dir("~/Downloads/ATAC_Workshop/ATAC_Data/ATAC_BAM_forCounting/", full.names = TRUE,    pattern = "*.\\.bam\$")myCounts <- summarizeOverlaps(consensusToCount, bamsToCount, singleEnd = FALSE)colnames(myCounts) <- c("HindBrain_1", "HindBrain_2", "Kidney_1", "Kidney_2", "Liver_1",    "Liver_2")``

## 3. DESeq2

``library(DESeq2)load("data/myCounts.RData")Group <- factor(c("HindBrain", "HindBrain", "Kidney", "Kidney", "Liver", "Liver"))metaData <- data.frame(Group, row.names = colnames(myCounts))atacDDS <- DESeqDataSetFromMatrix(assay(myCounts), metaData, ~Group, rowRanges = rowRanges(myCounts))atacDDS <- DESeq(atacDDS)``

``KidneyMinusHindbrain <- results(atacDDS, c("Group", "Kidney", "HindBrain"), format = "GRanges")KidneyMinusHindbrain <- KidneyMinusHindbrain[order(KidneyMinusHindbrain\$pvalue)]KidneyMinusHindbrain``

``library(TxDb.Mmusculus.UCSC.mm10.knownGene)toOverLap <- promoters(TxDb.Mmusculus.UCSC.mm10.knownGene, 500, 500)KidneyMinusHindbrain <- KidneyMinusHindbrain[(!is.na(KidneyMinusHindbrain\$padj) &    KidneyMinusHindbrain\$padj < 0.05) & KidneyMinusHindbrain %over% toOverLap, ]myReport <- makebedtable(KidneyMinusHindbrain, "KidneyMinusHindbrain.html", getwd())browseURL(myReport)``

## 4. 差异注释

``anno_KidneyMinusHindbrain <- annotatePeak(KidneyMinusHindbrain, TxDb = TxDb.Mmusculus.UCSC.mm10.knownGene,    verbose = FALSE)DB_ATAC <- as.data.frame(anno_KidneyMinusHindbrain)DB_ATAC[1, ]``

``library(clusterProfiler)go <- enrichGO(DB_ATAC\$geneId, OrgDb = "org.Mm.eg.db", ont = "BP", maxGSSize = 5000)go[1:2, 1:6]``

V1