### Importing libraries
library(DESeq2)
## Loading required package: S4Vectors
## Loading required package: stats4
## Loading required package: BiocGenerics
##
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
##
## anyDuplicated, aperm, append, as.data.frame, basename, cbind,
## colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
## get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
## match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
## Position, rank, rbind, Reduce, rownames, sapply, saveRDS, setdiff,
## table, tapply, union, unique, unsplit, which.max, which.min
##
## Attaching package: 'S4Vectors'
## The following object is masked from 'package:utils':
##
## findMatches
## The following objects are masked from 'package:base':
##
## expand.grid, I, unname
## Loading required package: IRanges
## Loading required package: GenomicRanges
## Loading required package: GenomeInfoDb
## Loading required package: SummarizedExperiment
## Loading required package: MatrixGenerics
## Loading required package: matrixStats
## Warning: package 'matrixStats' was built under R version 4.4.3
##
## Attaching package: 'MatrixGenerics'
## The following objects are masked from '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, colVars, colWeightedMads,
## colWeightedMeans, colWeightedMedians, colWeightedSds,
## 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, rowRanks,
## rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
## rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
## rowWeightedSds, rowWeightedVars
## Loading required package: Biobase
## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.
##
## Attaching package: 'Biobase'
## The following object is masked from 'package:MatrixGenerics':
##
## rowMedians
## The following objects are masked from 'package:matrixStats':
##
## anyMissing, rowMedians
library(org.Mm.eg.db)
## Loading required package: AnnotationDbi
##
library(TxDb.Mmusculus.UCSC.mm39.knownGene)
## Loading required package: GenomicFeatures
library(ChIPseeker)
##
## ChIPseeker v1.42.1 Learn more at https://yulab-smu.top/contribution-knowledge-mining/
##
## Please cite:
##
## Qianwen Wang, Ming Li, Tianzhi Wu, Li Zhan, Lin Li, Meijun Chen, Wenqin
## Xie, Zijing Xie, Erqiang Hu, Shuangbin Xu, Guangchuang Yu. Exploring
## epigenomic datasets by ChIPseeker. Current Protocols. 2022, 2(10): e585
library(GenomicRanges)
library(VennDiagram)
## Loading required package: grid
## Loading required package: futile.logger
## Warning: package 'futile.logger' was built under R version 4.4.3
library(ggplot2)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:AnnotationDbi':
##
## select
## The following object is masked from 'package:Biobase':
##
## combine
## The following object is masked from 'package:matrixStats':
##
## count
## The following objects are masked from 'package:GenomicRanges':
##
## intersect, setdiff, union
## The following object is masked from 'package:GenomeInfoDb':
##
## intersect
## The following objects are masked from 'package:IRanges':
##
## collapse, desc, intersect, setdiff, slice, union
## The following objects are masked from 'package:S4Vectors':
##
## first, intersect, rename, setdiff, setequal, union
## The following objects are masked from 'package:BiocGenerics':
##
## combine, intersect, setdiff, union
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tibble)
## Warning: package 'tibble' was built under R version 4.4.3
library(tidyr)
## Warning: package 'tidyr' was built under R version 4.4.3
##
## Attaching package: 'tidyr'
## The following object is masked from 'package:S4Vectors':
##
## expand
library(rtracklayer)
## Warning: package 'rtracklayer' was built under R version 4.4.3
### Reading count table
counts <- read.table("/home/drp46/project_pi_nd392/drp46/YK_2026_ATAC/Main_2026/data/bam/counttable_raw_downsampled_v2.txt",header = TRUE)
count_matrix <- as.matrix(counts[, 4:ncol(counts)])
### Seeting row_name as the chr_start_end
rownames(count_matrix) <- paste(counts$chrom,counts$start,counts$end,sep="_")
keep_peaks <- !grepl("^X_|Y_", rownames(count_matrix))
count_matrix <- count_matrix[keep_peaks, ]
### Metadata
coldata <- data.frame(
row.names = colnames(count_matrix),
genotype = c(
"WT","WT",
"KO","KO",
"KO","KO",
"WT","WT"
),
treatment = c(
"RO","RO",
"Asy","Asy",
"RO","RO",
"Asy","Asy"
)
)
coldata$genotype <- factor(coldata$genotype,levels = c("WT", "KO"))
coldata$treatment <- factor(coldata$treatment,levels = c("Asy", "RO"))
head(count_matrix)
## WTRO.2_S5_L002 WTRO.3_S6_L002 KOAsy.2_S8_L002 KOAsy.3_S9_L002
## 1_3319344_3319798 48 54 41 26
## 1_3515187_3515509 74 90 39 34
## 1_3563678_3564390 108 226 153 170
## 1_3584858_3585607 459 543 456 378
## 1_3681429_3681752 90 128 54 47
## 1_3741642_3742858 323 420 211 179
## KORO.2_S11_L002 KORO.3_S12_L002 WTAsy.2_S2_L002
## 1_3319344_3319798 34 24 168
## 1_3515187_3515509 59 29 89
## 1_3563678_3564390 126 133 315
## 1_3584858_3585607 436 446 736
## 1_3681429_3681752 65 48 145
## 1_3741642_3742858 281 241 342
## WTAsy.3_S3_L002
## 1_3319344_3319798 72
## 1_3515187_3515509 72
## 1_3563678_3564390 208
## 1_3584858_3585607 658
## 1_3681429_3681752 99
## 1_3741642_3742858 302
### Running DESeq2
dds <- DESeqDataSetFromMatrix(countData = count_matrix,colData = coldata,design = ~ genotype + treatment + genotype:treatment)
dds <- DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
sizeFactors(dds)
## WTRO.2_S5_L002 WTRO.3_S6_L002 KOAsy.2_S8_L002 KOAsy.3_S9_L002 KORO.2_S11_L002
## 0.7955230 0.9497626 1.0937403 0.9896128 1.0693573
## KORO.3_S12_L002 WTAsy.2_S2_L002 WTAsy.3_S3_L002
## 1.2187636 1.0669932 0.9353547
### PCA plot
vsd <- vst(dds, blind = TRUE)
mat <- assay(vsd)
pca <- prcomp(t(mat))
percentVar <- (pca$sdev^2) / sum(pca$sdev^2)
pcaData <- data.frame(pca$x,genotype = colData(vsd)$genotype,treatment = colData(vsd)$treatment)
ggplot(
pcaData,
aes(
x = PC1,
y = PC2,
color = genotype,
shape = treatment
)
) +
geom_point(size = 4) +
xlab(
paste0(
"PC1: ",
round(percentVar[1] * 100, 2),
"% variance"
)
) +
ylab(
paste0(
"PC2: ",
round(percentVar[2] * 100, 2),
"% variance"
)
) +
theme_minimal()
### First step for annotating the peaks
peak_coordinates <- data.frame(
peak_id = rownames(count_matrix)
) %>%
separate(
peak_id,
into = c("chr", "start", "end"),
sep = "_",
remove = FALSE,
convert = TRUE
)
head(peak_coordinates)
## peak_id chr start end
## 1 1_3319344_3319798 1 3319344 3319798
## 2 1_3515187_3515509 1 3515187 3515509
## 3 1_3563678_3564390 1 3563678 3564390
## 4 1_3584858_3585607 1 3584858 3585607
## 5 1_3681429_3681752 1 3681429 3681752
## 6 1_3741642_3742858 1 3741642 3742858
### Second and final step to annotate the peaks
peak_coordinates$chr <- paste0("chr", peak_coordinates$chr)
peaks_gr <- GRanges(seqnames = peak_coordinates$chr,ranges = IRanges(start = peak_coordinates$start,end = peak_coordinates$end))
mcols(peaks_gr)$peak_id <- peak_coordinates$peak_id
peakAnno <- annotatePeak(
peaks_gr,
TxDb = TxDb.Mmusculus.UCSC.mm39.knownGene,
tssRegion = c(-3000, 3000),
annoDb = "org.Mm.eg.db"
)
## >> preparing features information... 2026-08-08 23:53:04
## >> identifying nearest features... 2026-08-08 23:53:04
## Warning in .merge_two_Seqinfo_objects(x, y): Each of the 2 combined objects has sequence levels not in the other:
## - in 'x': chrGL456233.2, chrGL456392.1, chrJH584295.1, chrJH584304.1, chrMT, chrMU069434.1, chrMU069435.1
## - in 'y': chrX, chrY, chrM, chr1_GL456210v1_random, chr1_GL456211v1_random, chr1_GL456212v1_random, chr1_GL456221v1_random, chr1_GL456239v1_random, chr1_MU069434v1_random, chr4_JH584295v1_random, chr5_GL456354v1_random, chr5_JH584296v1_random, chr5_JH584297v1_random, chr5_JH584298v1_random, chr5_JH584299v1_random, chr7_GL456219v1_random, chrX_GL456233v2_random, chrY_JH584300v1_random, chrY_JH584301v1_random, chrY_JH584302v1_random, chrY_JH584303v1_random, chrUn_GL456359v1, chrUn_GL456360v1, chrUn_GL456366v1, chrUn_GL456367v1, chrUn_GL456368v1, chrUn_GL456370v1, chrUn_GL456372v1, chrUn_GL456378v1, chrUn_GL456379v1, chrUn_GL456381v1, chrUn_GL456382v1, chrUn_GL456383v1, chrUn_GL456385v1, chrUn_GL456387v1, chrUn_GL456389v1, chrUn_GL456390v1, chrUn_GL456392v1, chrUn_GL456394v1, chrUn_GL456396v1, chrUn_JH584304v1, chrUn_MU069435v1
## Make sure to always combine/compare objects based on the same reference
## genome (use suppressWarnings() to suppress this warning).
## >> calculating distance from peak to TSS... 2026-08-08 23:53:05
## >> assigning genomic annotation... 2026-08-08 23:53:05
## Warning in .merge_two_Seqinfo_objects(x, y): Each of the 2 combined objects has sequence levels not in the other:
## - in 'x': chrGL456233.2, chrGL456392.1, chrJH584295.1, chrJH584304.1, chrMT, chrMU069434.1, chrMU069435.1
## - in 'y': chrX, chrY, chrM, chr1_GL456210v1_random, chr1_GL456211v1_random, chr1_GL456212v1_random, chr1_GL456221v1_random, chr1_GL456239v1_random, chr1_MU069434v1_random, chr4_JH584295v1_random, chr5_GL456354v1_random, chr5_JH584296v1_random, chr5_JH584297v1_random, chr5_JH584298v1_random, chr5_JH584299v1_random, chr7_GL456219v1_random, chrX_GL456233v2_random, chrY_JH584300v1_random, chrY_JH584301v1_random, chrY_JH584302v1_random, chrY_JH584303v1_random, chrUn_GL456359v1, chrUn_GL456360v1, chrUn_GL456366v1, chrUn_GL456367v1, chrUn_GL456368v1, chrUn_GL456370v1, chrUn_GL456372v1, chrUn_GL456378v1, chrUn_GL456379v1, chrUn_GL456381v1, chrUn_GL456382v1, chrUn_GL456383v1, chrUn_GL456385v1, chrUn_GL456387v1, chrUn_GL456389v1, chrUn_GL456390v1, chrUn_GL456392v1, chrUn_GL456394v1, chrUn_GL456396v1, chrUn_JH584304v1, chrUn_MU069435v1
## Make sure to always combine/compare objects based on the same reference
## genome (use suppressWarnings() to suppress this warning).
## Warning in .merge_two_Seqinfo_objects(x, y): Each of the 2 combined objects has sequence levels not in the other:
## - in 'x': chrGL456233.2, chrGL456392.1, chrJH584295.1, chrJH584304.1, chrMT, chrMU069434.1, chrMU069435.1
## - in 'y': chrX, chrY, chrM, chr1_GL456210v1_random, chr1_GL456211v1_random, chr1_GL456212v1_random, chr1_GL456221v1_random, chr1_GL456239v1_random, chr1_MU069434v1_random, chr4_JH584295v1_random, chr5_GL456354v1_random, chr5_JH584296v1_random, chr5_JH584297v1_random, chr5_JH584298v1_random, chr5_JH584299v1_random, chr7_GL456219v1_random, chrX_GL456233v2_random, chrY_JH584300v1_random, chrY_JH584301v1_random, chrY_JH584302v1_random, chrY_JH584303v1_random, chrUn_GL456359v1, chrUn_GL456360v1, chrUn_GL456366v1, chrUn_GL456367v1, chrUn_GL456368v1, chrUn_GL456370v1, chrUn_GL456372v1, chrUn_GL456378v1, chrUn_GL456379v1, chrUn_GL456381v1, chrUn_GL456382v1, chrUn_GL456383v1, chrUn_GL456385v1, chrUn_GL456387v1, chrUn_GL456389v1, chrUn_GL456390v1, chrUn_GL456392v1, chrUn_GL456394v1, chrUn_GL456396v1, chrUn_JH584304v1, chrUn_MU069435v1
## Make sure to always combine/compare objects based on the same reference
## genome (use suppressWarnings() to suppress this warning).
## Warning in .merge_two_Seqinfo_objects(x, y): Each of the 2 combined objects has sequence levels not in the other:
## - in 'x': chrGL456233.2, chrGL456392.1, chrJH584295.1, chrJH584304.1, chrMT, chrMU069434.1, chrMU069435.1
## - in 'y': chrX, chrY, chrM, chr1_GL456210v1_random, chr1_GL456211v1_random, chr1_GL456212v1_random, chr1_GL456221v1_random, chr1_GL456239v1_random, chr1_MU069434v1_random, chr4_JH584295v1_random, chr5_GL456354v1_random, chr5_JH584296v1_random, chr5_JH584297v1_random, chr5_JH584298v1_random, chr5_JH584299v1_random, chr7_GL456219v1_random, chrX_GL456233v2_random, chrY_JH584300v1_random, chrY_JH584301v1_random, chrY_JH584302v1_random, chrY_JH584303v1_random, chrUn_GL456359v1, chrUn_GL456360v1, chrUn_GL456366v1, chrUn_GL456367v1, chrUn_GL456368v1, chrUn_GL456370v1, chrUn_GL456372v1, chrUn_GL456378v1, chrUn_GL456379v1, chrUn_GL456381v1, chrUn_GL456382v1, chrUn_GL456383v1, chrUn_GL456385v1, chrUn_GL456387v1, chrUn_GL456389v1, chrUn_GL456390v1, chrUn_GL456392v1, chrUn_GL456394v1, chrUn_GL456396v1, chrUn_JH584304v1, chrUn_MU069435v1
## Make sure to always combine/compare objects based on the same reference
## genome (use suppressWarnings() to suppress this warning).
## Warning in .merge_two_Seqinfo_objects(x, y): Each of the 2 combined objects has sequence levels not in the other:
## - in 'x': chrGL456233.2, chrGL456392.1, chrJH584295.1, chrJH584304.1, chrMT, chrMU069434.1, chrMU069435.1
## - in 'y': chrX, chrY, chrM, chr1_GL456210v1_random, chr1_GL456211v1_random, chr1_GL456212v1_random, chr1_GL456221v1_random, chr1_GL456239v1_random, chr1_MU069434v1_random, chr4_JH584295v1_random, chr5_GL456354v1_random, chr5_JH584296v1_random, chr5_JH584297v1_random, chr5_JH584298v1_random, chr5_JH584299v1_random, chr7_GL456219v1_random, chrX_GL456233v2_random, chrY_JH584300v1_random, chrY_JH584301v1_random, chrY_JH584302v1_random, chrY_JH584303v1_random, chrUn_GL456359v1, chrUn_GL456360v1, chrUn_GL456366v1, chrUn_GL456367v1, chrUn_GL456368v1, chrUn_GL456370v1, chrUn_GL456372v1, chrUn_GL456378v1, chrUn_GL456379v1, chrUn_GL456381v1, chrUn_GL456382v1, chrUn_GL456383v1, chrUn_GL456385v1, chrUn_GL456387v1, chrUn_GL456389v1, chrUn_GL456390v1, chrUn_GL456392v1, chrUn_GL456394v1, chrUn_GL456396v1, chrUn_JH584304v1, chrUn_MU069435v1
## Make sure to always combine/compare objects based on the same reference
## genome (use suppressWarnings() to suppress this warning).
## >> adding gene annotation... 2026-08-08 23:53:17
## 'select()' returned 1:many mapping between keys and columns
## >> assigning chromosome lengths 2026-08-08 23:53:17
## >> done... 2026-08-08 23:53:17
annotation_df <- as.data.frame(peakAnno)
### Storing the raw counts for master_dataframe generation
raw_counts_df <- as.data.frame(counts(dds, normalized = FALSE)) %>%
rownames_to_column("peak_id")
head(raw_counts_df)
## peak_id WTRO.2_S5_L002 WTRO.3_S6_L002 KOAsy.2_S8_L002
## 1 1_3319344_3319798 48 54 41
## 2 1_3515187_3515509 74 90 39
## 3 1_3563678_3564390 108 226 153
## 4 1_3584858_3585607 459 543 456
## 5 1_3681429_3681752 90 128 54
## 6 1_3741642_3742858 323 420 211
## KOAsy.3_S9_L002 KORO.2_S11_L002 KORO.3_S12_L002 WTAsy.2_S2_L002
## 1 26 34 24 168
## 2 34 59 29 89
## 3 170 126 133 315
## 4 378 436 446 736
## 5 47 65 48 145
## 6 179 281 241 342
## WTAsy.3_S3_L002
## 1 72
## 2 72
## 3 208
## 4 658
## 5 99
## 6 302
### Adding prefix raw_
colnames(raw_counts_df)[-1] <- paste0("raw_",colnames(raw_counts_df)[-1])
head(raw_counts_df)
## peak_id raw_WTRO.2_S5_L002 raw_WTRO.3_S6_L002 raw_KOAsy.2_S8_L002
## 1 1_3319344_3319798 48 54 41
## 2 1_3515187_3515509 74 90 39
## 3 1_3563678_3564390 108 226 153
## 4 1_3584858_3585607 459 543 456
## 5 1_3681429_3681752 90 128 54
## 6 1_3741642_3742858 323 420 211
## raw_KOAsy.3_S9_L002 raw_KORO.2_S11_L002 raw_KORO.3_S12_L002
## 1 26 34 24
## 2 34 59 29
## 3 170 126 133
## 4 378 436 446
## 5 47 65 48
## 6 179 281 241
## raw_WTAsy.2_S2_L002 raw_WTAsy.3_S3_L002
## 1 168 72
## 2 89 72
## 3 315 208
## 4 736 658
## 5 145 99
## 6 342 302
### Storing the normalized values for masterdataframe generation
normalized_counts_df <- as.data.frame(counts(dds, normalized = TRUE)) %>%
rownames_to_column("peak_id")
colnames(normalized_counts_df)[-1] <- paste0("norm_",colnames(normalized_counts_df)[-1])
peak_loadings_df <- as.data.frame(pca$rotation) %>%
rownames_to_column("peak_id")
peak_PCA_annotation <- annotation_df %>%
left_join(peak_loadings_df,by = "peak_id")
head(peak_PCA_annotation)
## seqnames start end width strand peak_id
## 1 chr1 3319344 3319798 455 * 1_3319344_3319798
## 2 chr1 3515187 3515509 323 * 1_3515187_3515509
## 3 chr1 3563678 3564390 713 * 1_3563678_3564390
## 4 chr1 3584858 3585607 750 * 1_3584858_3585607
## 5 chr1 3681429 3681752 324 * 1_3681429_3681752
## 6 chr1 3741642 3742858 1217 * 1_3741642_3742858
## annotation geneChr
## 1 Intron (ENSMUST00000070533.5/497097, intron 2 of 2) 1
## 2 Intron (ENSMUST00000070533.5/497097, intron 1 of 2) 1
## 3 Intron (ENSMUST00000161581.2/ENSMUST00000161581.2, intron 1 of 1) 1
## 4 Intron (ENSMUST00000070533.5/497097, intron 1 of 2) 1
## 5 Intron (ENSMUST00000070533.5/497097, intron 1 of 2) 1
## 6 Promoter (<=1kb) 1
## geneStart geneEnd geneLength geneStrand geneId transcriptId
## 1 3276124 3286567 10444 2 497097 ENSMUST00000162897.2
## 2 3284705 3741721 457017 2 497097 ENSMUST00000070533.5
## 3 3284705 3741721 457017 2 497097 ENSMUST00000070533.5
## 4 3284705 3741721 457017 2 497097 ENSMUST00000070533.5
## 5 3284705 3741721 457017 2 497097 ENSMUST00000070533.5
## 6 3284705 3741721 457017 2 497097 ENSMUST00000070533.5
## distanceToTSS ENSEMBL SYMBOL GENENAME
## 1 -32777 ENSMUSG00000051951 Xkr4 X-linked Kx blood group related 4
## 2 226212 ENSMUSG00000051951 Xkr4 X-linked Kx blood group related 4
## 3 177331 ENSMUSG00000051951 Xkr4 X-linked Kx blood group related 4
## 4 156114 ENSMUSG00000051951 Xkr4 X-linked Kx blood group related 4
## 5 59969 ENSMUSG00000051951 Xkr4 X-linked Kx blood group related 4
## 6 0 ENSMUSG00000051951 Xkr4 X-linked Kx blood group related 4
## PC1 PC2 PC3 PC4 PC5
## 1 -0.003687284 -0.0060195456 -0.0021031129 0.007282426 0.005847236
## 2 -0.004005324 -0.0013736069 -0.0026820003 0.001273850 -0.002440695
## 3 -0.002876285 -0.0060804715 0.0030551848 0.004653895 -0.000967440
## 4 -0.005316011 -0.0049121443 -0.0002945318 0.005568899 0.001439713
## 5 -0.004677738 -0.0027211150 -0.0014958028 0.001708438 -0.001213980
## 6 -0.006043373 0.0001870505 -0.0027713596 -0.001101333 -0.004912735
## PC6 PC7 PC8
## 1 -0.001061368 0.0057394575 0.030080570
## 2 -0.003808964 0.0003800777 -0.084052629
## 3 -0.006166469 0.0071418264 -0.008519383
## 4 0.001543501 -0.0025874704 -0.126459944
## 5 -0.002973707 0.0044076768 -0.202620845
## 6 -0.001949480 0.0036481728 -0.289245982
### Generating Master sheet
master_peak_df <- annotation_df %>%
left_join(raw_counts_df,by = "peak_id") %>%
left_join(normalized_counts_df,by = "peak_id") %>%
left_join(peak_loadings_df,by = "peak_id")
head(master_peak_df)
## seqnames start end width strand peak_id
## 1 chr1 3319344 3319798 455 * 1_3319344_3319798
## 2 chr1 3515187 3515509 323 * 1_3515187_3515509
## 3 chr1 3563678 3564390 713 * 1_3563678_3564390
## 4 chr1 3584858 3585607 750 * 1_3584858_3585607
## 5 chr1 3681429 3681752 324 * 1_3681429_3681752
## 6 chr1 3741642 3742858 1217 * 1_3741642_3742858
## annotation geneChr
## 1 Intron (ENSMUST00000070533.5/497097, intron 2 of 2) 1
## 2 Intron (ENSMUST00000070533.5/497097, intron 1 of 2) 1
## 3 Intron (ENSMUST00000161581.2/ENSMUST00000161581.2, intron 1 of 1) 1
## 4 Intron (ENSMUST00000070533.5/497097, intron 1 of 2) 1
## 5 Intron (ENSMUST00000070533.5/497097, intron 1 of 2) 1
## 6 Promoter (<=1kb) 1
## geneStart geneEnd geneLength geneStrand geneId transcriptId
## 1 3276124 3286567 10444 2 497097 ENSMUST00000162897.2
## 2 3284705 3741721 457017 2 497097 ENSMUST00000070533.5
## 3 3284705 3741721 457017 2 497097 ENSMUST00000070533.5
## 4 3284705 3741721 457017 2 497097 ENSMUST00000070533.5
## 5 3284705 3741721 457017 2 497097 ENSMUST00000070533.5
## 6 3284705 3741721 457017 2 497097 ENSMUST00000070533.5
## distanceToTSS ENSEMBL SYMBOL GENENAME
## 1 -32777 ENSMUSG00000051951 Xkr4 X-linked Kx blood group related 4
## 2 226212 ENSMUSG00000051951 Xkr4 X-linked Kx blood group related 4
## 3 177331 ENSMUSG00000051951 Xkr4 X-linked Kx blood group related 4
## 4 156114 ENSMUSG00000051951 Xkr4 X-linked Kx blood group related 4
## 5 59969 ENSMUSG00000051951 Xkr4 X-linked Kx blood group related 4
## 6 0 ENSMUSG00000051951 Xkr4 X-linked Kx blood group related 4
## raw_WTRO.2_S5_L002 raw_WTRO.3_S6_L002 raw_KOAsy.2_S8_L002 raw_KOAsy.3_S9_L002
## 1 48 54 41 26
## 2 74 90 39 34
## 3 108 226 153 170
## 4 459 543 456 378
## 5 90 128 54 47
## 6 323 420 211 179
## raw_KORO.2_S11_L002 raw_KORO.3_S12_L002 raw_WTAsy.2_S2_L002
## 1 34 24 168
## 2 59 29 89
## 3 126 133 315
## 4 436 446 736
## 5 65 48 145
## 6 281 241 342
## raw_WTAsy.3_S3_L002 norm_WTRO.2_S5_L002 norm_WTRO.3_S6_L002
## 1 72 60.33767 56.85631
## 2 72 93.02057 94.76052
## 3 208 135.75975 237.95419
## 4 658 576.97894 571.72179
## 5 99 113.13313 134.77051
## 6 302 406.02222 442.21575
## norm_KOAsy.2_S8_L002 norm_KOAsy.3_S9_L002 norm_KORO.2_S11_L002
## 1 37.48605 26.27290 31.79480
## 2 35.65746 34.35687 55.17332
## 3 139.88695 171.78436 117.82778
## 4 416.91798 381.96759 407.72152
## 5 49.37187 47.49332 60.78417
## 6 192.91599 180.87883 262.77465
## norm_KORO.3_S12_L002 norm_WTAsy.2_S2_L002 norm_WTAsy.3_S3_L002 PC1
## 1 19.69209 157.45180 76.97615 -0.003687284
## 2 23.79461 83.41196 76.97615 -0.004005324
## 3 109.12698 295.22212 222.37554 -0.002876285
## 4 365.94462 689.78882 703.47647 -0.005316011
## 5 39.38417 135.89590 105.84220 -0.004677738
## 6 197.74138 320.52687 322.87218 -0.006043373
## PC2 PC3 PC4 PC5 PC6
## 1 -0.0060195456 -0.0021031129 0.007282426 0.005847236 -0.001061368
## 2 -0.0013736069 -0.0026820003 0.001273850 -0.002440695 -0.003808964
## 3 -0.0060804715 0.0030551848 0.004653895 -0.000967440 -0.006166469
## 4 -0.0049121443 -0.0002945318 0.005568899 0.001439713 0.001543501
## 5 -0.0027211150 -0.0014958028 0.001708438 -0.001213980 -0.002973707
## 6 0.0001870505 -0.0027713596 -0.001101333 -0.004912735 -0.001949480
## PC7 PC8
## 1 0.0057394575 0.030080570
## 2 0.0003800777 -0.084052629
## 3 0.0071418264 -0.008519383
## 4 -0.0025874704 -0.126459944
## 5 0.0044076768 -0.202620845
## 6 0.0036481728 -0.289245982
### Treatment effect in WT
res_WTRO_vs_WTAsy <- results(dds,name = "treatment_RO_vs_Asy")
### Treatment effect in KO == treatment_effect + interaction
res_KORO_vs_KOAsy <- results(dds,contrast = list(c("treatment_RO_vs_Asy","genotypeKO.treatmentRO")))
### Genotype effect in RO == genotype_effect + interaction
res_KORO_vs_WTRO <- results(dds,contrast = list(c("genotype_KO_vs_WT","genotypeKO.treatmentRO")))
### Genotype effect in Asy
res_KOAsy_vs_WTAsy <- results(dds,name = "genotype_KO_vs_WT")
make_volcano <- function(res, title) {
df <- as.data.frame(res) %>%
rownames_to_column("peak_id") %>%
mutate(
category = case_when(
!is.na(padj) & padj < 0.05 & log2FoldChange > 1 ~ "Up",
!is.na(padj) & padj < 0.05 & log2FoldChange < -1 ~ "Down",
TRUE ~ "Not significant"
),
minus_log10_padj = -log10(
pmax(padj, .Machine$double.xmin)
)
)
n_up <- sum(df$category == "Up", na.rm = TRUE)
n_down <- sum(df$category == "Down", na.rm = TRUE)
ggplot(
df,
aes(
x = log2FoldChange,
y = minus_log10_padj,
color = category
)
) +
geom_point(
alpha = 0.6,
size = 1.5
) +
scale_color_manual(
values = c(
"Up" = "red",
"Down" = "blue",
"Not significant" = "grey"
)
) +
geom_vline(
xintercept = c(-1, 1),
linetype = "dashed"
) +
geom_hline(
yintercept = -log10(0.05),
linetype = "dashed"
) +
annotate(
"text",
x = Inf,
y = Inf,
label = paste0("Up: ", n_up),
hjust = 1.2,
vjust = 2,
color = "red",
size = 5
) +
annotate(
"text",
x = -Inf,
y = Inf,
label = paste0("Down: ", n_down),
hjust = -0.2,
vjust = 2,
color = "blue",
size = 5
) +
labs(
title = title,
x = "log2 Fold Change",
y = "-log10 adjusted p-value"
) +
theme_minimal() +
theme(
plot.title = element_text(
face = "bold",
size = 16
),
legend.position = "bottom"
)
}
p1 <- make_volcano(res_WTRO_vs_WTAsy,"WTRO vs WTAsy")
p2 <- make_volcano(res_KORO_vs_KOAsy,"KORO vs KOAsy")
p3 <- make_volcano(res_KORO_vs_WTRO,"KORO vs WTRO")
p4 <- make_volcano(res_KOAsy_vs_WTAsy,"KOAsy vs WTAsy")
p1
## Warning: Removed 70014 rows containing missing values or values outside the scale range
## (`geom_point()`).
p2
## Warning: Removed 5252 rows containing missing values or values outside the scale range
## (`geom_point()`).
p3
p4
### Generating Master sheet
master_sheet <- master_peak_df
prepare_deseq_result <- function(res, comparison_name) {
res_df <- as.data.frame(res) %>%
rownames_to_column("peak_id")
# Rename all DESeq2 columns except peak_id
colnames(res_df)[-1] <- paste0(comparison_name,"_",colnames(res_df)[-1])
return(res_df)
}
df_WTRO_vs_WTAsy <- prepare_deseq_result(res_WTRO_vs_WTAsy,"WTRO_vs_WTAsy")
df_KORO_vs_KOAsy <- prepare_deseq_result(res_KORO_vs_KOAsy,"KORO_vs_KOAsy")
df_KORO_vs_WTRO <- prepare_deseq_result(res_KORO_vs_WTRO,"KORO_vs_WTRO")
df_KOAsy_vs_WTAsy <- prepare_deseq_result(res_KOAsy_vs_WTAsy,"KOAsy_vs_WTAsy")
master_sheet <- master_sheet %>%
left_join(df_WTRO_vs_WTAsy,by = "peak_id") %>%
left_join(df_KORO_vs_KOAsy,by = "peak_id") %>%
left_join(df_KORO_vs_WTRO,by = "peak_id") %>%
left_join(df_KOAsy_vs_WTAsy,by = "peak_id")
write.table(master_sheet,"/home/drp46/project_pi_nd392/drp46/YK_2026_ATAC/Main_2026/Downsampled_ATAC_Master_Sheet_All_Data_v2.tsv",row.names = FALSE, sep = "\t",quote = FALSE)