OSCR

Allele-specific chromatin architecture shapes imprinted domains and coordinates a distal enhancer and antisense transcription at the mouse Mest-Copg2 domain.

Code ↔ Paper

8 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 8 matches · 3 of them tie a paragraph to a whole file, not to given lines: weak matches, whose lines are not tinted
  1. [1] § Methods › Capture Hi-C analysis ↔ Figure1/BB_CHiC_GRCm39_Figure1.sh, lines 49–89 · score 0.85 · hicup2juicer, SNPsplit, usemid, digested, v0, g2
  2. [2] § Results › Methylation-sensitive CTCF at ICRs shapes allelic architecture ↔ Figure4/BB_CHiC_GRCm39_contactanalysisChromHMM_FIgure4.R, lines 1127–1186 · score 0.61 · Cdkn1c, Kcnq1ot1, imprinted genes, Peg10, Peg3, Snrpn
  3. [3] § Results › Methylation-sensitive CTCF at ICRs shapes allelic architecture ↔ Figure2/BB_CHiC_GRCm39_Figure2.sh, the whole file · a weak match · score 0.60 · sDMRs, gDMRs, Methylated alleles, unmethylated, domain
  4. [4] § Results › Methylation-sensitive CTCF at ICRs shapes allelic architecture ↔ Figure2/BB_CHiC_GRCm39_Figure2.sh, the whole file · a weak match · score 0.54 · sDMRs, gDMRs, Coolpuppy, Density, unmethylated, cortex
  5. [5] § Results › Methylation-sensitive CTCF at ICRs shapes allelic architecture ↔ Figure4/BB_CHiC_GRCm39_contactanalysisChromHMM_FIgure4.R, lines 1348–1395 · score 0.54 · Kcnq1ot1, Peg12, Usp29, Kcnk9, Peg13, Sgce
  6. [6] § Results › Imprinted domains exhibit allele-specific chromatin architectures ↔ Figure4/BB_CHiC_GRCm39_contactanalysisChromHMM_FIgure4.R, lines 1127–1186 · score 0.54 · Kcnq1ot1, Peg12, Usp29, Kcnk9, Peg13, Sgce
  7. [7] § Methods › Capture Hi-C analysis ↔ FigureSupp10/BB_CHiC_GRCm39_FigureSupp10.sh, the whole file · a weak match · score 0.52 · Juicer pre, strain, g2, g1, tool, hic
  8. [8] § Methods › Capture Hi-C analysis ↔ Figure1/BB_CHiC_GRCm39_Figure1.sh, lines 49–89 · score 0.51 · Juicer pre, g2, g1, tool, filtered, hic

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

R · 1,400 lines · 53 KB · no license · 3 matches

  1. library(dplyr)
  2. library(ggplot2)
  3. library(scales)
  4. library(tidyverse)
  5. library(valr)
  6. library(palr)
  7. library(GenomicRanges)
  10. #personalized theme
  11. my_theme = theme(
  12. text = element_text(size = 10),
  13. axis.title.x = element_text(size = 12),
  14. axis.title.y = element_text(size = 12),
  15. axis.text = element_text(size = 10, colour = 'black'),
  16. axis.text.x = element_text(vjust = 0.5, size = 10, angle=90),
  17. legend.title=element_text(size=0),
  18. legend.text=element_text(size=12),
  19. legend.position = "top",
  20. #legend.position = "none",
  21. plot.title = element_text(lineheight=.8, face="bold", size = 16),
  22. panel.border = element_blank(),
  23. panel.background = element_rect(fill = 'white'),
  24. panel.grid.minor = element_blank(),
  25. panel.grid.major = element_blank(),
  26. axis.line.x = element_line(colour = 'black', size=0.5, linetype='solid'),
  27. axis.line.y = element_line(colour = 'black', size=0.5, linetype='solid'))
  30. #==============================
  31. func.paircounts.mat <- function(viewpoint) {
  32. #mat positive region
  33. mat_chr_pos <- mat_chr_reg %>%
  34. filter(mat_chr_reg$V3 > viewpoint-2500 & mat_chr_reg$V3 < viewpoint+2500)
  35. #mat negative region
  36. mat_chr_neg <- mat_chr_reg %>%
  37. filter(mat_chr_reg$V5 > viewpoint-2500 & mat_chr_reg$V5 < viewpoint+2500)
  38. #mat merge both regions
  39. mat_chr_pos2 <- as.data.frame(mat_chr_pos$V1)
  40. mat_chr_pos2$V2 <- mat_chr_pos$V5
  41. colnames(mat_chr_pos2) <- c("V1","V2")
  42. mat_chr_neg2 <- as.data.frame(mat_chr_neg$V1)
  43. mat_chr_neg2$V2 <- mat_chr_neg$V3
  44. colnames(mat_chr_neg2) <- c("V1","V2")
  45. mat_chr_plot <- rbind(mat_chr_neg2,mat_chr_pos2)
  46. mat_chr_plot$chr <- chrnumb
  47. mat_chr_plot$allele <- "mat"
  48. #remove proximal contacts
  49. mat_chr_plot <- mat_chr_plot %>% filter(V2 < viewpoint-50000 | V2 > viewpoint+50000)
  50. return(mat_chr_plot)
  51. }
  53. func.paircounts.pat <- function(viewpoint) {
  54. #pat positive region
  55. pat_chr_pos <- pat_chr_reg %>%
  56. filter(pat_chr_reg$V3 > viewpoint-2500 & pat_chr_reg$V3 < viewpoint+2500)
  57. #pat negative region
  58. pat_chr_neg <- pat_chr_reg %>%
  59. filter(pat_chr_reg$V5 > viewpoint-2500 & pat_chr_reg$V5 < viewpoint+2500)
  60. #pat merge both regions
  61. pat_chr_pos2 <- as.data.frame(pat_chr_pos$V1)
  62. pat_chr_pos2$V2 <- pat_chr_pos$V5
  63. colnames(pat_chr_pos2) <- c("V1","V2")
  64. pat_chr_neg2 <- as.data.frame(pat_chr_neg$V1)
  65. pat_chr_neg2$V2 <- pat_chr_neg$V3
  66. colnames(pat_chr_neg2) <- c("V1","V2")
  67. pat_chr_plot <- rbind(pat_chr_neg2,pat_chr_pos2)
  68. pat_chr_plot$chr <- chrnumb
  69. pat_chr_plot$allele <- "pat"
  70. #remove proximal contacts
  71. pat_chr_plot <- pat_chr_plot %>% filter(V2 < viewpoint-50000 | V2 > viewpoint+50000)
  72. return(pat_chr_plot)
  73. }
  75. #==============================
  77. func.activebed <- function(active) {
  78. #contacts from active TSS
  79. bed1 <- as.data.frame(active$chr)
  80. bed1 <- cbind(bed1, active$start)
  81. bed1 <- cbind(bed1, active$V2)
  82. bed1 <- cbind(bed1, active$V1)
  83. colnames(bed1) <- c("chrom","start","end","name")
  84. bed1$chrom <- sub("^", "chr", bed1$chrom)
  85. return(bed1)
  86. }
  88. func.inactivebed <- function(inactive) {
  89. #contacts from inactive TSS
  90. bed2 <- as.data.frame(inactive$chr)
  91. bed2 <- cbind(bed2, inactive$start)
  92. bed2 <- cbind(bed2, inactive$V2)
  93. bed2 <- cbind(bed2, inactive$V1)
  94. colnames(bed2) <- c("chrom","start","end","name")
  95. bed2$chrom <- sub("^", "chr", bed2$chrom)
  96. return(bed2)
  97. }
  99. func.biallele1bed <- function(biallele1) {
  100. #contacts from biallelic TSS
  101. bed3 <- as.data.frame(biallele1$chr)
  102. bed3 <- cbind(bed3, biallele1$start)
  103. bed3 <- cbind(bed3, biallele1$V2)
  104. bed3 <- cbind(bed3, biallele1$V1)
  105. colnames(bed3) <- c("chrom","start","end","name")
  106. bed3$chrom <- sub("^", "chr", bed3$chrom)
  107. return(bed3)
  108. }
  110. func.biallele2bed <- function(biallele2) {
  111. #contacts from biallelic TSS
  112. bed4 <- as.data.frame(biallele2$chr)
  113. bed4 <- cbind(bed4, biallele2$start)
  114. bed4 <- cbind(bed4, biallele2$V2)
  115. bed4 <- cbind(bed4, biallele2$V1)
  116. colnames(bed4) <- c("chrom","start","end","name")
  117. bed4$chrom <- sub("^", "chr", bed4$chrom)
  118. return(bed4)
  119. }
  120. #==============================
  121. cumm.biallele1bed <- NULL
  122. cumm.biallele2bed <- NULL
  123. #==============================
  128. ## this script may not be very straightforward
  129. ## basically, you are grabbing contact information from each expressed gene allele and grouping them according to their expression pattern
  130. ## note: total contact numbers in medium files are normalized to equalize mat and pat numbers
  132. #==============================
  133. #===Region1 file analysis
  134. mat_chr_reg <- read.table("Prejuicer/EXP012.CHiC_Mat_prejuicer_Q30_reg1.medium", header = FALSE, sep = '\t')
  135. mat_chr_reg <- mat_chr_reg[c('V1','V3','V4','V7','V8','V2','V6')]
  136. colnames(mat_chr_reg) <- c('V1','V2','V3','V4','V5','V6','V7')
  137. pat_chr_reg <- read.table("Prejuicer/EXP012.CHiC_Pat_prejuicer_Q30_reg1.medium", header = FALSE, sep = '\t')
  138. pat_chr_reg <- pat_chr_reg[c('V1','V3','V4','V7','V8','V2','V6')]
  139. colnames(pat_chr_reg) <- c('V1','V2','V3','V4','V5','V6','V7')
  140. #==============================
  141. #Sgce pat TSS
  142. gene <- "Sgce"
  143. viewpoint <- 4747180
  144. chrnumb <- 6
  145. exp <- "pat"
  146. noexp <- "mat"
  147. #counting all the contacts from the viewpoint 5kb bin
  148. mat_chr_plot <- func.paircounts.mat(viewpoint)
  149. pat_chr_plot <- func.paircounts.pat(viewpoint)
  150. #which allele is transcriptionally active?
  151. if (exp == "pat"){
  152. active <- pat_chr_plot
  153. inactive <- mat_chr_plot
  154. } else if (exp == "mat"){
  155. active <- mat_chr_plot
  156. inactive <- pat_chr_plot
  157. }
  158. #add viewpoint
  159. active$start <- viewpoint
  160. inactive$start <- viewpoint
  161. #generating bed
  162. activebed <- func.activebed(active)
  163. inactivebed <- func.inactivebed(inactive)
  164. activebed$IG <- gene
  165. inactivebed$IG <- gene
  166. #cummulative
  167. cumm.activebed <- activebed
  168. cumm.inactivebed <- inactivebed
  169. #==============================
  170. #Peg10 pat TSS
  171. gene <- "Peg10"
  172. viewpoint <- 4747306
  173. chrnumb <- 6
  174. exp <- "pat"
  175. noexp <- "mat"
  176. #counting all the contacts from the viewpoint 5kb bin
  177. mat_chr_plot <- func.paircounts.mat(viewpoint)
  178. pat_chr_plot <- func.paircounts.pat(viewpoint)
  179. #which allele is transcriptionally active?
  180. if (exp == "pat"){
  181. active <- pat_chr_plot
  182. inactive <- mat_chr_plot
  183. } else if (exp == "mat"){
  184. active <- mat_chr_plot
  185. inactive <- pat_chr_plot
  186. }
  187. #add viewpoint
  188. active$start <- viewpoint
  189. inactive$start <- viewpoint
  190. #generating bed
  191. activebed <- func.activebed(active)
  192. inactivebed <- func.inactivebed(inactive)
  193. activebed$IG <- gene
  194. inactivebed$IG <- gene
  195. #cummulative
  196. cumm.activebed <- rbind(unique(cumm.activebed),activebed)
  197. cumm.inactivebed <- rbind(unique(cumm.inactivebed),inactivebed)
  198. #==============================
  199. biallelicbed <- read.table("Features/EXP022_CN_BiallelicRNAseq_capture.bed", header=FALSE, sep="\t", stringsAsFactors=FALSE)
  200. filter_gr <- GRanges(seqnames=6, ranges=IRanges(start=3330467 + 1, end=8138330))
  201. query_gr <- GRanges(seqnames=biallelicbed[[1]],
  202. ranges=IRanges(start=biallelicbed[[2]] + 1, end=biallelicbed[[3]]))
  203. hits <- findOverlaps(query_gr, filter_gr, type = "within")
  204. filtered_query <- biallelicbed[queryHits(hits), ]
  205. biallelicregion <- filtered_query
  206. biallelicregion$viewpoint[biallelicregion$V6 == "-"] <- biallelicregion$V3[biallelicregion$V6 == "-"]
  207. biallelicregion$viewpoint[biallelicregion$V6 == "+"] <- biallelicregion$V2[biallelicregion$V6 == "+"]
  208. biallelicregion <- biallelicregion[,c("V4","V1","viewpoint")]
  209. colnames(biallelicregion) <- c("gene","chrnumb","viewpoint")
  210. i=1
  211. for (i in seq_len(nrow(biallelicregion))) {
  212. gene <- biallelicregion$gene[i]
  213. viewpoint <- biallelicregion$viewpoint[i]
  214. chrnumb <- biallelicregion$chrnumb[i]
  215. exp <- "bi"
  216. mat_chr_plot <- func.paircounts.mat(viewpoint)
  217. pat_chr_plot <- func.paircounts.pat(viewpoint)
  218. if (exp == "bi"){
  219. biallele1 <- pat_chr_plot
  220. biallele2 <- mat_chr_plot
  221. }
  222. biallele1$start <- viewpoint
  223. biallele2$start <- viewpoint
  224. biallele1bed <- func.biallele1bed(biallele1)
  225. biallele2bed <- func.biallele2bed(biallele2)
  226. biallele1bed$IG <- gene
  227. biallele2bed$IG <- gene
  228. cumm.biallele1bed <- rbind(unique(cumm.biallele1bed),biallele1bed)
  229. cumm.biallele2bed <- rbind(unique(cumm.biallele2bed),biallele2bed)
  230. }
  231. #==============================
  237. #==============================
  238. #===Region2 file analysis
  239. mat_chr_reg <- read.table("Prejuicer/EXP012.CHiC_Mat_prejuicer_Q30_reg2.medium", header = FALSE, sep = '\t')
  240. mat_chr_reg <- mat_chr_reg[c('V1','V3','V4','V7','V8','V2','V6')]
  241. colnames(mat_chr_reg) <- c('V1','V2','V3','V4','V5','V6','V7')
  242. pat_chr_reg <- read.table("Prejuicer/EXP012.CHiC_Pat_prejuicer_Q30_reg2.medium", header = FALSE, sep = '\t')
  243. pat_chr_reg <- pat_chr_reg[c('V1','V3','V4','V7','V8','V2','V6')]
  244. colnames(pat_chr_reg) <- c('V1','V2','V3','V4','V5','V6','V7')
  245. #==============================
  246. #Mest pat TSS
  247. gene <- "Mest"
  248. viewpoint <- 30737996
  249. chrnumb <- 6
  250. exp <- "pat"
  251. noexp <- "mat"
  252. #counting all the contacts from the viewpoint 5kb bin
  253. mat_chr_plot <- func.paircounts.mat(viewpoint)
  254. pat_chr_plot <- func.paircounts.pat(viewpoint)
  255. #which allele is transcriptionally active?
  256. if (exp == "pat"){
  257. active <- pat_chr_plot
  258. inactive <- mat_chr_plot
  259. } else if (exp == "mat"){
  260. active <- mat_chr_plot
  261. inactive <- pat_chr_plot
  262. }
  263. #add viewpoint
  264. active$start <- viewpoint
  265. inactive$start <- viewpoint
  266. #generating bed
  267. activebed <- func.activebed(active)
  268. inactivebed <- func.inactivebed(inactive)
  269. activebed$IG <- gene
  270. inactivebed$IG <- gene
  271. #cummulative
  272. cumm.activebed <- rbind(unique(cumm.activebed),activebed)
  273. cumm.inactivebed <- rbind(unique(cumm.inactivebed),inactivebed)
  274. #==============================
  275. #Copg2 TSS
  276. gene <- "Copg2"
  277. viewpoint <- 30873712
  278. chrnumb <- 6
  279. exp <- "mat"
  280. noexp <- "pat"
  281. #counting all the contacts from the viewpoint 5kb bin
  282. mat_chr_plot <- func.paircounts.mat(viewpoint)
  283. pat_chr_plot <- func.paircounts.pat(viewpoint)
  284. #which allele is transcriptionally active?
  285. if (exp == "pat"){
  286. active <- pat_chr_plot
  287. inactive <- mat_chr_plot
  288. } else if (exp == "mat"){
  289. active <- mat_chr_plot
  290. inactive <- pat_chr_plot
  291. }
  292. #add viewpoint
  293. active$start <- viewpoint
  294. inactive$start <- viewpoint
  295. #generating bed
  296. activebed <- func.activebed(active)
  297. inactivebed <- func.inactivebed(inactive)
  298. activebed$IG <- gene
  299. inactivebed$IG <- gene
  300. #cummulative
  301. cumm.activebed <- rbind(unique(cumm.activebed),activebed)
  302. cumm.inactivebed <- rbind(unique(cumm.inactivebed),inactivebed)
  303. #==============================
  304. biallelicbed <- read.table("Features/EXP022_CN_BiallelicRNAseq_capture.bed", header=FALSE, sep="\t", stringsAsFactors=FALSE)
  305. filter_gr <- GRanges(seqnames=6, ranges=IRanges(start=30068177+1, end=31232042))
  306. query_gr <- GRanges(seqnames=biallelicbed[[1]],
  307. ranges=IRanges(start=biallelicbed[[2]] + 1, end=biallelicbed[[3]]))
  308. hits <- findOverlaps(query_gr, filter_gr, type = "within")
  309. filtered_query <- biallelicbed[queryHits(hits), ]
  310. biallelicregion <- filtered_query
  311. biallelicregion$viewpoint[biallelicregion$V6 == "-"] <- biallelicregion$V3[biallelicregion$V6 == "-"]
  312. biallelicregion$viewpoint[biallelicregion$V6 == "+"] <- biallelicregion$V2[biallelicregion$V6 == "+"]
  313. biallelicregion <- biallelicregion[,c("V4","V1","viewpoint")]
  314. colnames(biallelicregion) <- c("gene","chrnumb","viewpoint")
  315. i=1
  316. for (i in seq_len(nrow(biallelicregion))) {
  317. gene <- biallelicregion$gene[i]
  318. viewpoint <- biallelicregion$viewpoint[i]
  319. chrnumb <- biallelicregion$chrnumb[i]
  320. exp <- "bi"
  321. mat_chr_plot <- func.paircounts.mat(viewpoint)
  322. pat_chr_plot <- func.paircounts.pat(viewpoint)
  323. if (exp == "bi"){
  324. biallele1 <- pat_chr_plot
  325. biallele2 <- mat_chr_plot
  326. }
  327. biallele1$start <- viewpoint
  328. biallele2$start <- viewpoint
  329. biallele1bed <- func.biallele1bed(biallele1)
  330. biallele2bed <- func.biallele2bed(biallele2)
  331. biallele1bed$IG <- gene
  332. biallele2bed$IG <- gene
  333. cumm.biallele1bed <- rbind(unique(cumm.biallele1bed),biallele1bed)
  334. cumm.biallele2bed <- rbind(unique(cumm.biallele2bed),biallele2bed)
  335. }
  336. #==============================
  342. #==============================
  343. #===Region3 pairs file analysis
  344. mat_chr_reg <- read.table("Prejuicer/EXP012.CHiC_Mat_prejuicer_Q30_reg3.medium", header = FALSE, sep = '\t')
  345. mat_chr_reg <- mat_chr_reg[c('V1','V3','V4','V7','V8','V2','V6')]
  346. colnames(mat_chr_reg) <- c('V1','V2','V3','V4','V5','V6','V7')
  347. pat_chr_reg <- read.table("Prejuicer/EXP012.CHiC_Pat_prejuicer_Q30_reg3.medium", header = FALSE, sep = '\t')
  348. pat_chr_reg <- pat_chr_reg[c('V1','V3','V4','V7','V8','V2','V6')]
  349. colnames(pat_chr_reg) <- c('V1','V2','V3','V4','V5','V6','V7')
  350. #==============================
  351. #==============================
  352. #Peg3 TSS
  353. gene <- "Peg3"
  354. viewpoint <- 6733443
  355. chrnumb <- 7
  356. exp <- "pat"
  357. noexp <- "mat"
  358. #counting all the contacts from the viewpoint 5kb bin
  359. mat_chr_plot <- func.paircounts.mat(viewpoint)
  360. pat_chr_plot <- func.paircounts.pat(viewpoint)
  361. #which allele is transcriptionally active?
  362. if (exp == "pat"){
  363. active <- pat_chr_plot
  364. inactive <- mat_chr_plot
  365. } else if (exp == "mat"){
  366. active <- mat_chr_plot
  367. inactive <- pat_chr_plot
  368. }
  369. #add viewpoint
  370. active$start <- viewpoint
  371. inactive$start <- viewpoint
  372. #generating bed
  373. activebed <- func.activebed(active)
  374. inactivebed <- func.inactivebed(inactive)
  375. activebed$IG <- gene
  376. inactivebed$IG <- gene
  377. #cummulative
  378. cumm.activebed <- rbind(unique(cumm.activebed),activebed)
  379. cumm.inactivebed <- rbind(unique(cumm.inactivebed),inactivebed)
  380. #==============================
  381. #Usp29 TSS
  382. gene <- "Usp29"
  383. viewpoint <- 6733560
  384. chrnumb <- 7
  385. exp <- "pat"
  386. noexp <- "mat"
  387. #counting all the contacts from the viewpoint 5kb bin
  388. mat_chr_plot <- func.paircounts.mat(viewpoint)
  389. pat_chr_plot <- func.paircounts.pat(viewpoint)
  390. #which allele is transcriptionally active?
  391. if (exp == "pat"){
  392. active <- pat_chr_plot
  393. inactive <- mat_chr_plot
  394. } else if (exp == "mat"){
  395. active <- mat_chr_plot
  396. inactive <- pat_chr_plot
  397. }
  398. #add viewpoint
  399. active$start <- viewpoint
  400. inactive$start <- viewpoint
  401. #generating bed
  402. activebed <- func.activebed(active)
  403. inactivebed <- func.inactivebed(inactive)
  404. activebed$IG <- gene
  405. inactivebed$IG <- gene
  406. #cummulative
  407. cumm.activebed <- rbind(unique(cumm.activebed),activebed)
  408. cumm.inactivebed <- rbind(unique(cumm.inactivebed),inactivebed)
  409. #==============================
  410. biallelicbed <- read.table("Features/EXP022_CN_BiallelicRNAseq_capture.bed", header=FALSE, sep="\t", stringsAsFactors=FALSE)
  411. filter_gr <- GRanges(seqnames=7, ranges=IRanges(start=6118501 + 1, end=7334692))
  412. query_gr <- GRanges(seqnames=biallelicbed[[1]],
  413. ranges=IRanges(start=biallelicbed[[2]] + 1, end=biallelicbed[[3]]))
  414. hits <- findOverlaps(query_gr, filter_gr, type = "within")
  415. filtered_query <- biallelicbed[queryHits(hits), ]
  416. biallelicregion <- filtered_query
  417. biallelicregion$viewpoint[biallelicregion$V6 == "-"] <- biallelicregion$V3[biallelicregion$V6 == "-"]
  418. biallelicregion$viewpoint[biallelicregion$V6 == "+"] <- biallelicregion$V2[biallelicregion$V6 == "+"]
  419. biallelicregion <- biallelicregion[,c("V4","V1","viewpoint")]
  420. colnames(biallelicregion) <- c("gene","chrnumb","viewpoint")
  421. i=1
  422. for (i in seq_len(nrow(biallelicregion))) {
  423. gene <- biallelicregion$gene[i]
  424. viewpoint <- biallelicregion$viewpoint[i]
  425. chrnumb <- biallelicregion$chrnumb[i]
  426. exp <- "bi"
  427. mat_chr_plot <- func.paircounts.mat(viewpoint)
  428. pat_chr_plot <- func.paircounts.pat(viewpoint)
  429. if (exp == "bi"){
  430. biallele1 <- pat_chr_plot
  431. biallele2 <- mat_chr_plot
  432. }
  433. biallele1$start <- viewpoint
  434. biallele2$start <- viewpoint
  435. biallele1bed <- func.biallele1bed(biallele1)
  436. biallele2bed <- func.biallele2bed(biallele2)
  437. biallele1bed$IG <- gene
  438. biallele2bed$IG <- gene
  439. cumm.biallele1bed <- rbind(unique(cumm.biallele1bed),biallele1bed)
  440. cumm.biallele2bed <- rbind(unique(cumm.biallele2bed),biallele2bed)
  441. }
  442. #==============================
  447. #==============================
  448. #===Region4 file analysis
  449. mat_chr_reg <- read.table("Prejuicer/EXP012.CHiC_Mat_prejuicer_Q30_reg4.medium", header = FALSE, sep = '\t')
  450. mat_chr_reg <- mat_chr_reg[c('V1','V3','V4','V7','V8','V2','V6')]
  451. colnames(mat_chr_reg) <- c('V1','V2','V3','V4','V5','V6','V7')
  452. pat_chr_reg <- read.table("Prejuicer/EXP012.CHiC_Pat_prejuicer_Q30_reg4.medium", header = FALSE, sep = '\t')
  453. pat_chr_reg <- pat_chr_reg[c('V1','V3','V4','V7','V8','V2','V6')]
  454. colnames(pat_chr_reg) <- c('V1','V2','V3','V4','V5','V6','V7')
  455. #==============================
  456. #==============================
  457. #Ube3a TSS
  458. gene <- "Ube3a"
  459. viewpoint <- 58878498
  460. chrnumb <- 7
  461. exp <- "mat"
  462. noexp <- "pat"
  463. #counting all the contacts from the viewpoint 5kb bin
  464. mat_chr_plot <- func.paircounts.mat(viewpoint)
  465. pat_chr_plot <- func.paircounts.pat(viewpoint)
  466. #which allele is transcriptionally active?
  467. if (exp == "pat"){
  468. active <- pat_chr_plot
  469. inactive <- mat_chr_plot
  470. } else if (exp == "mat"){
  471. active <- mat_chr_plot
  472. inactive <- pat_chr_plot
  473. }
  474. #add viewpoint
  475. active$start <- viewpoint
  476. inactive$start <- viewpoint
  477. #generating bed
  478. activebed <- func.activebed(active)
  479. inactivebed <- func.inactivebed(inactive)
  480. activebed$IG <- gene
  481. inactivebed$IG <- gene
  482. #cummulative
  483. cumm.activebed <- rbind(unique(cumm.activebed),activebed)
  484. cumm.inactivebed <- rbind(unique(cumm.inactivebed),inactivebed)
  485. #==============================
  486. #Snrpn TSS
  487. gene <- "Snrpn"
  488. viewpoint <- 59789967
  489. chrnumb <- 7
  490. exp <- "pat"
  491. noexp <- "mat"
  492. #counting all the contacts from the viewpoint 5kb bin
  493. mat_chr_plot <- func.paircounts.mat(viewpoint)
  494. pat_chr_plot <- func.paircounts.pat(viewpoint)
  495. #which allele is transcriptionally active?
  496. if (exp == "pat"){
  497. active <- pat_chr_plot
  498. inactive <- mat_chr_plot
  499. } else if (exp == "mat"){
  500. active <- mat_chr_plot
  501. inactive <- pat_chr_plot
  502. }
  503. #add viewpoint
  504. active$start <- viewpoint
  505. inactive$start <- viewpoint
  506. #generating bed
  507. activebed <- func.activebed(active)
  508. inactivebed <- func.inactivebed(inactive)
  509. activebed$IG <- gene
  510. inactivebed$IG <- gene
  511. #cummulative
  512. cumm.activebed <- rbind(unique(cumm.activebed),activebed)
  513. cumm.inactivebed <- rbind(unique(cumm.inactivebed),inactivebed)
  514. #==============================
  515. #A230006K03Rik TSS
  516. gene <- "A230006K03Rik"
  517. viewpoint <- 60961397
  518. chrnumb <- 7
  519. exp <- "pat"
  520. noexp <- "mat"
  521. #counting all the contacts from the viewpoint 5kb bin
  522. mat_chr_plot <- func.paircounts.mat(viewpoint)
  523. pat_chr_plot <- func.paircounts.pat(viewpoint)
  524. #which allele is transcriptionally active?
  525. if (exp == "pat"){
  526. active <- pat_chr_plot
  527. inactive <- mat_chr_plot
  528. } else if (exp == "mat"){
  529. active <- mat_chr_plot
  530. inactive <- pat_chr_plot
  531. }
  532. #add viewpoint
  533. active$start <- viewpoint
  534. inactive$start <- viewpoint
  535. #generating bed
  536. activebed <- func.activebed(active)
  537. inactivebed <- func.inactivebed(inactive)
  538. activebed$IG <- gene
  539. inactivebed$IG <- gene
  540. #cummulative
  541. cumm.activebed <- rbind(unique(cumm.activebed),activebed)
  542. cumm.inactivebed <- rbind(unique(cumm.inactivebed),inactivebed)
  543. #==============================
  544. #B230209E15Rik TSS
  545. gene <- "B230209E15Rik"
  546. viewpoint <- 61265075
  547. chrnumb <- 7
  548. exp <- "pat"
  549. noexp <- "mat"
  550. #counting all the contacts from the viewpoint 5kb bin
  551. mat_chr_plot <- func.paircounts.mat(viewpoint)
  552. pat_chr_plot <- func.paircounts.pat(viewpoint)
  553. #which allele is transcriptionally active?
  554. if (exp == "pat"){
  555. active <- pat_chr_plot
  556. inactive <- mat_chr_plot
  557. } else if (exp == "mat"){
  558. active <- mat_chr_plot
  559. inactive <- pat_chr_plot
  560. }
  561. #add viewpoint
  562. active$start <- viewpoint
  563. inactive$start <- viewpoint
  564. #generating bed
  565. activebed <- func.activebed(active)
  566. inactivebed <- func.inactivebed(inactive)
  567. activebed$IG <- gene
  568. inactivebed$IG <- gene
  569. #cummulative
  570. cumm.activebed <- rbind(unique(cumm.activebed),activebed)
  571. cumm.inactivebed <- rbind(unique(cumm.inactivebed),inactivebed)
  572. #==============================
  573. #A330076H08Rik TSS
  574. gene <- "A330076H08Rik"
  575. viewpoint <- 61632103
  576. chrnumb <- 7
  577. exp <- "pat"
  578. noexp <- "mat"
  579. #counting all the contacts from the viewpoint 5kb bin
  580. mat_chr_plot <- func.paircounts.mat(viewpoint)
  581. pat_chr_plot <- func.paircounts.pat(viewpoint)
  582. #which allele is transcriptionally active?
  583. if (exp == "pat"){
  584. active <- pat_chr_plot
  585. inactive <- mat_chr_plot
  586. } else if (exp == "mat"){
  587. active <- mat_chr_plot
  588. inactive <- pat_chr_plot
  589. }
  590. #add viewpoint
  591. active$start <- viewpoint
  592. inactive$start <- viewpoint
  593. #generating bed
  594. activebed <- func.activebed(active)
  595. inactivebed <- func.inactivebed(inactive)
  596. activebed$IG <- gene
  597. inactivebed$IG <- gene
  598. #cummulative
  599. cumm.activebed <- rbind(unique(cumm.activebed),activebed)
  600. cumm.inactivebed <- rbind(unique(cumm.inactivebed),inactivebed)
  601. #==============================
  602. #Ndn TSS
  603. gene <- "Ndn"
  604. viewpoint <- 61998025
  605. chrnumb <- 7
  606. exp <- "pat"
  607. noexp <- "mat"
  608. #counting all the contacts from the viewpoint 5kb bin
  609. mat_chr_plot <- func.paircounts.mat(viewpoint)
  610. pat_chr_plot <- func.paircounts.pat(viewpoint)
  611. #which allele is transcriptionally active?
  612. if (exp == "pat"){
  613. active <- pat_chr_plot
  614. inactive <- mat_chr_plot
  615. } else if (exp == "mat"){
  616. active <- mat_chr_plot
  617. inactive <- pat_chr_plot
  618. }
  619. #add viewpoint
  620. active$start <- viewpoint
  621. inactive$start <- viewpoint
  622. #generating bed
  623. activebed <- func.activebed(active)
  624. inactivebed <- func.inactivebed(inactive)
  625. activebed$IG <- gene
  626. inactivebed$IG <- gene
  627. #cummulative
  628. cumm.activebed <- rbind(unique(cumm.activebed),activebed)
  629. cumm.inactivebed <- rbind(unique(cumm.inactivebed),inactivebed)
  630. #==============================
  631. #Magel2 TSS
  632. gene <- "Magel2"
  633. viewpoint <- 62026727
  634. chrnumb <- 7
  635. exp <- "pat"
  636. noexp <- "mat"
  637. #counting all the contacts from the viewpoint 5kb bin
  638. mat_chr_plot <- func.paircounts.mat(viewpoint)
  639. pat_chr_plot <- func.paircounts.pat(viewpoint)
  640. #which allele is transcriptionally active?
  641. if (exp == "pat"){
  642. active <- pat_chr_plot
  643. inactive <- mat_chr_plot
  644. } else if (exp == "mat"){
  645. active <- mat_chr_plot
  646. inactive <- pat_chr_plot
  647. }
  648. #add viewpoint
  649. active$start <- viewpoint
  650. inactive$start <- viewpoint
  651. #generating bed
  652. activebed <- func.activebed(active)
  653. inactivebed <- func.inactivebed(inactive)
  654. activebed$IG <- gene
  655. inactivebed$IG <- gene
  656. #cummulative
  657. cumm.activebed <- rbind(unique(cumm.activebed),activebed)
  658. cumm.inactivebed <- rbind(unique(cumm.inactivebed),inactivebed)
  659. #==============================
  660. #Mkrn3 TSS
  661. gene <- "Mkrn3"
  662. viewpoint <- 62069901
  663. chrnumb <- 7
  664. exp <- "pat"
  665. noexp <- "mat"
  666. #counting all the contacts from the viewpoint 5kb bin
  667. mat_chr_plot <- func.paircounts.mat(viewpoint)
  668. pat_chr_plot <- func.paircounts.pat(viewpoint)
  669. #which allele is transcriptionally active?
  670. if (exp == "pat"){
  671. active <- pat_chr_plot
  672. inactive <- mat_chr_plot
  673. } else if (exp == "mat"){
  674. active <- mat_chr_plot
  675. inactive <- pat_chr_plot
  676. }
  677. #add viewpoint
  678. active$start <- viewpoint
  679. inactive$start <- viewpoint
  680. #generating bed
  681. activebed <- func.activebed(active)
  682. inactivebed <- func.inactivebed(inactive)
  683. activebed$IG <- gene
  684. inactivebed$IG <- gene
  685. #cummulative
  686. cumm.activebed <- rbind(unique(cumm.activebed),activebed)
  687. cumm.inactivebed <- rbind(unique(cumm.inactivebed),inactivebed)
  688. #==============================
  689. #Peg12 TSS
  690. gene <- "Peg12"
  691. viewpoint <- 62114258
  692. chrnumb <- 7
  693. exp <- "pat"
  694. noexp <- "mat"
  695. #counting all the contacts from the viewpoint 5kb bin
  696. mat_chr_plot <- func.paircounts.mat(viewpoint)
  697. pat_chr_plot <- func.paircounts.pat(viewpoint)
  698. #which allele is transcriptionally active?
  699. if (exp == "pat"){
  700. active <- pat_chr_plot
  701. inactive <- mat_chr_plot
  702. } else if (exp == "mat"){
  703. active <- mat_chr_plot
  704. inactive <- pat_chr_plot
  705. }
  706. #add viewpoint
  707. active$start <- viewpoint
  708. inactive$start <- viewpoint
  709. #generating bed
  710. activebed <- func.activebed(active)
  711. inactivebed <- func.inactivebed(inactive)
  712. activebed$IG <- gene
  713. inactivebed$IG <- gene
  714. #cummulative
  715. cumm.activebed <- rbind(unique(cumm.activebed),activebed)
  716. cumm.inactivebed <- rbind(unique(cumm.inactivebed),inactivebed)
  717. #==============================
  718. biallelicbed <- read.table("Features/EXP022_CN_BiallelicRNAseq_capture.bed", header=FALSE, sep="\t", stringsAsFactors=FALSE)
  719. filter_gr <- GRanges(seqnames=7, ranges=IRanges(start=55345629 + 1, end=62207321))
  720. query_gr <- GRanges(seqnames=biallelicbed[[1]],
  721. ranges=IRanges(start=biallelicbed[[2]] + 1, end=biallelicbed[[3]]))
  722. hits <- findOverlaps(query_gr, filter_gr, type = "within")
  723. filtered_query <- biallelicbed[queryHits(hits), ]
  724. biallelicregion <- filtered_query
  725. biallelicregion$viewpoint[biallelicregion$V6 == "-"] <- biallelicregion$V3[biallelicregion$V6 == "-"]
  726. biallelicregion$viewpoint[biallelicregion$V6 == "+"] <- biallelicregion$V2[biallelicregion$V6 == "+"]
  727. biallelicregion <- biallelicregion[,c("V4","V1","viewpoint")]
  728. colnames(biallelicregion) <- c("gene","chrnumb","viewpoint")
  729. i=1
  730. for (i in seq_len(nrow(biallelicregion))) {
  731. gene <- biallelicregion$gene[i]
  732. viewpoint <- biallelicregion$viewpoint[i]
  733. chrnumb <- biallelicregion$chrnumb[i]
  734. exp <- "bi"
  735. mat_chr_plot <- func.paircounts.mat(viewpoint)
  736. pat_chr_plot <- func.paircounts.pat(viewpoint)
  737. if (exp == "bi"){
  738. biallele1 <- pat_chr_plot
  739. biallele2 <- mat_chr_plot
  740. }
  741. biallele1$start <- viewpoint
  742. biallele2$start <- viewpoint
  743. biallele1bed <- func.biallele1bed(biallele1)
  744. biallele2bed <- func.biallele2bed(biallele2)
  745. biallele1bed$IG <- gene
  746. biallele2bed$IG <- gene
  747. cumm.biallele1bed <- rbind(unique(cumm.biallele1bed),biallele1bed)
  748. cumm.biallele2bed <- rbind(unique(cumm.biallele2bed),biallele2bed)
  749. }
  750. #==============================
  756. #==============================
  757. #===Region5 file analysis
  758. mat_chr_reg <- read.table("Prejuicer/EXP012.CHiC_Mat_prejuicer_Q30_reg5a.medium", header = FALSE, sep = '\t')
  759. mat_chr_reg <- mat_chr_reg[c('V1','V3','V4','V7','V8','V2','V6')]
  760. colnames(mat_chr_reg) <- c('V1','V2','V3','V4','V5','V6','V7')
  761. pat_chr_reg <- read.table("Prejuicer/EXP012.CHiC_Pat_prejuicer_Q30_reg5a.medium", header = FALSE, sep = '\t')
  762. pat_chr_reg <- pat_chr_reg[c('V1','V3','V4','V7','V8','V2','V6')]
  763. colnames(pat_chr_reg) <- c('V1','V2','V3','V4','V5','V6','V7')
  764. #==============================
  765. #H19 TSS not expressed in neurons
  766. #==============================
  767. #Igf2 TSS not expressed in neurons
  768. #==============================
  769. biallelicbed <- read.table("Features/EXP022_CN_BiallelicRNAseq_capture.bed", header=FALSE, sep="\t", stringsAsFactors=FALSE)
  770. filter_gr <- GRanges(seqnames=7, ranges=IRanges(start=141663846 + 1, end=142264758))
  771. query_gr <- GRanges(seqnames=biallelicbed[[1]],
  772. ranges=IRanges(start=biallelicbed[[2]] + 1, end=biallelicbed[[3]]))
  773. hits <- findOverlaps(query_gr, filter_gr, type = "within")
  774. filtered_query <- biallelicbed[queryHits(hits), ]
  775. biallelicregion <- filtered_query
  776. biallelicregion$viewpoint[biallelicregion$V6 == "-"] <- biallelicregion$V3[biallelicregion$V6 == "-"]
  777. biallelicregion$viewpoint[biallelicregion$V6 == "+"] <- biallelicregion$V2[biallelicregion$V6 == "+"]
  778. biallelicregion <- biallelicregion[,c("V4","V1","viewpoint")]
  779. colnames(biallelicregion) <- c("gene","chrnumb","viewpoint")
  780. i=1
  781. for (i in seq_len(nrow(biallelicregion))) {
  782. gene <- biallelicregion$gene[i]
  783. viewpoint <- biallelicregion$viewpoint[i]
  784. chrnumb <- biallelicregion$chrnumb[i]
  785. exp <- "bi"
  786. mat_chr_plot <- func.paircounts.mat(viewpoint)
  787. pat_chr_plot <- func.paircounts.pat(viewpoint)
  788. if (exp == "bi"){
  789. biallele1 <- pat_chr_plot
  790. biallele2 <- mat_chr_plot
  791. }
  792. biallele1$start <- viewpoint
  793. biallele2$start <- viewpoint
  794. biallele1bed <- func.biallele1bed(biallele1)
  795. biallele2bed <- func.biallele2bed(biallele2)
  796. biallele1bed$IG <- gene
  797. biallele2bed$IG <- gene
  798. cumm.biallele1bed <- rbind(unique(cumm.biallele1bed),biallele1bed)
  799. cumm.biallele2bed <- rbind(unique(cumm.biallele2bed),biallele2bed)
  800. }
  801. #==============================
  802. #===Region5-2
  803. mat_chr_reg <- read.table("Prejuicer/EXP012.CHiC_Mat_prejuicer_Q30_reg5b.medium", header = FALSE, sep = '\t')
  804. mat_chr_reg <- mat_chr_reg[c('V1','V3','V4','V7','V8','V2','V6')]
  805. colnames(mat_chr_reg) <- c('V1','V2','V3','V4','V5','V6','V7')
  806. pat_chr_reg <- read.table("Prejuicer/EXP012.CHiC_Pat_prejuicer_Q30_reg5b.medium", header = FALSE, sep = '\t')
  807. pat_chr_reg <- pat_chr_reg[c('V1','V3','V4','V7','V8','V2','V6')]
  808. colnames(pat_chr_reg) <- c('V1','V2','V3','V4','V5','V6','V7')
  809. #=============================
  810. #Kcnq1ot1 TSS
  811. gene <- "Kcnq1ot1"
  812. viewpoint <- 142850286
  813. chrnumb <- 7
  814. exp <- "pat"
  815. noexp <- "mat"
  816. #counting all the contacts from the viewpoint 5kb bin
  817. mat_chr_plot <- func.paircounts.mat(viewpoint)
  818. pat_chr_plot <- func.paircounts.pat(viewpoint)
  819. #which allele is transcriptionally active?
  820. if (exp == "pat"){
  821. active <- pat_chr_plot
  822. inactive <- mat_chr_plot
  823. } else if (exp == "mat"){
  824. active <- mat_chr_plot
  825. inactive <- pat_chr_plot
  826. }
  827. #add viewpoint
  828. active$start <- viewpoint
  829. inactive$start <- viewpoint
  830. #generating bed
  831. activebed <- func.activebed(active)
  832. inactivebed <- func.inactivebed(inactive)
  833. activebed$IG <- gene
  834. inactivebed$IG <- gene
  835. #cummulative
  836. cumm.activebed <- rbind(unique(cumm.activebed),activebed)
  837. cumm.inactivebed <- rbind(unique(cumm.inactivebed),inactivebed)
  838. #=============================
  839. #Cdkn1c TSS
  840. gene <- "Cdkn1c"
  841. viewpoint <- 143014735
  842. chrnumb <- 7
  843. exp <- "mat"
  844. noexp <- "pat"
  845. #counting all the contacts from the viewpoint 5kb bin
  846. mat_chr_plot <- func.paircounts.mat(viewpoint)
  847. pat_chr_plot <- func.paircounts.pat(viewpoint)
  848. #which allele is transcriptionally active?
  849. if (exp == "pat"){
  850. active <- pat_chr_plot
  851. inactive <- mat_chr_plot
  852. } else if (exp == "mat"){
  853. active <- mat_chr_plot
  854. inactive <- pat_chr_plot
  855. }
  856. #add viewpoint
  857. active$start <- viewpoint
  858. inactive$start <- viewpoint
  859. #generating bed
  860. activebed <- func.activebed(active)
  861. inactivebed <- func.inactivebed(inactive)
  862. activebed$IG <- gene
  863. inactivebed$IG <- gene
  864. #cummulative
  865. cumm.activebed <- rbind(unique(cumm.activebed),activebed)
  866. cumm.inactivebed <- rbind(unique(cumm.inactivebed),inactivebed)
  867. #=============================
  868. biallelicbed <- read.table("Features/EXP022_CN_BiallelicRNAseq_capture.bed", header=FALSE, sep="\t", stringsAsFactors=FALSE)
  869. filter_gr <- GRanges(seqnames=7, ranges=IRanges(start=142266803 + 1, end=144176724))
  870. query_gr <- GRanges(seqnames=biallelicbed[[1]],
  871. ranges=IRanges(start=biallelicbed[[2]] + 1, end=biallelicbed[[3]]))
  872. hits <- findOverlaps(query_gr, filter_gr, type = "within")
  873. filtered_query <- biallelicbed[queryHits(hits), ]
  874. biallelicregion <- filtered_query
  875. biallelicregion$viewpoint[biallelicregion$V6 == "-"] <- biallelicregion$V3[biallelicregion$V6 == "-"]
  876. biallelicregion$viewpoint[biallelicregion$V6 == "+"] <- biallelicregion$V2[biallelicregion$V6 == "+"]
  877. biallelicregion <- biallelicregion[,c("V4","V1","viewpoint")]
  878. colnames(biallelicregion) <- c("gene","chrnumb","viewpoint")
  879. i=1
  880. for (i in seq_len(nrow(biallelicregion))) {
  881. gene <- biallelicregion$gene[i]
  882. viewpoint <- biallelicregion$viewpoint[i]
  883. chrnumb <- biallelicregion$chrnumb[i]
  884. exp <- "bi"
  885. mat_chr_plot <- func.paircounts.mat(viewpoint)
  886. pat_chr_plot <- func.paircounts.pat(viewpoint)
  887. if (exp == "bi"){
  888. biallele1 <- pat_chr_plot
  889. biallele2 <- mat_chr_plot
  890. }
  891. biallele1$start <- viewpoint
  892. biallele2$start <- viewpoint
  893. biallele1bed <- func.biallele1bed(biallele1)
  894. biallele2bed <- func.biallele2bed(biallele2)
  895. biallele1bed$IG <- gene
  896. biallele2bed$IG <- gene
  897. cumm.biallele1bed <- rbind(unique(cumm.biallele1bed),biallele1bed)
  898. cumm.biallele2bed <- rbind(unique(cumm.biallele2bed),biallele2bed)
  899. }
  900. #==============================
  905. #==============================
  906. #===Region6 pairs file analysis
  907. mat_chr_reg <- read.table("Prejuicer/EXP012.CHiC_Mat_prejuicer_Q30_reg6.medium", header = FALSE, sep = '\t')
  908. mat_chr_reg <- mat_chr_reg[c('V1','V3','V4','V7','V8','V2','V6')]
  909. colnames(mat_chr_reg) <- c('V1','V2','V3','V4','V5','V6','V7')
  910. pat_chr_reg <- read.table("Prejuicer/EXP012.CHiC_Pat_prejuicer_Q30_reg6.medium", header = FALSE, sep = '\t')
  911. pat_chr_reg <- pat_chr_reg[c('V1','V3','V4','V7','V8','V2','V6')]
  912. colnames(pat_chr_reg) <- c('V1','V2','V3','V4','V5','V6','V7')
  913. #==============================
  914. #Dlk1 TSS
  915. #==============================
  916. #Meg3 TSS1
  917. gene <- "Meg3"
  918. viewpoint <- 109506879
  919. chrnumb <- 12
  920. exp <- "mat"
  921. noexp <- "pat"
  922. #counting all the contacts from the viewpoint 5kb bin
  923. mat_chr_plot <- func.paircounts.mat(viewpoint)
  924. pat_chr_plot <- func.paircounts.pat(viewpoint)
  925. #which allele is transcriptionally active?
  926. if (exp == "pat"){
  927. active <- pat_chr_plot
  928. inactive <- mat_chr_plot
  929. } else if (exp == "mat"){
  930. active <- mat_chr_plot
  931. inactive <- pat_chr_plot
  932. }
  933. #add viewpoint
  934. active$start <- viewpoint
  935. inactive$start <- viewpoint
  936. #generating bed
  937. activebed <- func.activebed(active)
  938. inactivebed <- func.inactivebed(inactive)
  939. activebed$IG <- gene
  940. inactivebed$IG <- gene
  941. #cummulative
  942. cumm.activebed <- rbind(unique(cumm.activebed),activebed)
  943. cumm.inactivebed <- rbind(unique(cumm.inactivebed),inactivebed)
  944. #==============================
  945. #Rtl1 TSS
  946. #==============================
  947. #Rian TSS polycistronic
  948. #==============================
  949. #Mirg TSS polycistronic
  950. #==============================
  951. #=============================
  952. biallelicbed <- read.table("Features/EXP022_CN_BiallelicRNAseq_capture.bed", header=FALSE, sep="\t", stringsAsFactors=FALSE)
  953. filter_gr <- GRanges(seqnames=12, ranges=IRanges(start=108178545 + 1, end=111651244))
  954. query_gr <- GRanges(seqnames=biallelicbed[[1]],
  955. ranges=IRanges(start=biallelicbed[[2]] + 1, end=biallelicbed[[3]]))
  956. hits <- findOverlaps(query_gr, filter_gr, type = "within")
  957. filtered_query <- biallelicbed[queryHits(hits), ]
  958. biallelicregion <- filtered_query
  959. biallelicregion$viewpoint[biallelicregion$V6 == "-"] <- biallelicregion$V3[biallelicregion$V6 == "-"]
  960. biallelicregion$viewpoint[biallelicregion$V6 == "+"] <- biallelicregion$V2[biallelicregion$V6 == "+"]
  961. biallelicregion <- biallelicregion[,c("V4","V1","viewpoint")]
  962. colnames(biallelicregion) <- c("gene","chrnumb","viewpoint")
  963. i=1
  964. for (i in seq_len(nrow(biallelicregion))) {
  965. gene <- biallelicregion$gene[i]
  966. viewpoint <- biallelicregion$viewpoint[i]
  967. chrnumb <- biallelicregion$chrnumb[i]
  968. exp <- "bi"
  969. mat_chr_plot <- func.paircounts.mat(viewpoint)
  970. pat_chr_plot <- func.paircounts.pat(viewpoint)
  971. if (exp == "bi"){
  972. biallele1 <- pat_chr_plot
  973. biallele2 <- mat_chr_plot
  974. }
  975. biallele1$start <- viewpoint
  976. biallele2$start <- viewpoint
  977. biallele1bed <- func.biallele1bed(biallele1)
  978. biallele2bed <- func.biallele2bed(biallele2)
  979. biallele1bed$IG <- gene
  980. biallele2bed$IG <- gene
  981. cumm.biallele1bed <- rbind(unique(cumm.biallele1bed),biallele1bed)
  982. cumm.biallele2bed <- rbind(unique(cumm.biallele2bed),biallele2bed)
  983. }
  984. #==============================
  990. #==============================
  991. #===Region7 pairs file analysis
  992. mat_chr_reg <- read.table("Prejuicer/EXP012.CHiC_Mat_prejuicer_Q30_reg7.medium", header = FALSE, sep = '\t')
  993. mat_chr_reg <- mat_chr_reg[c('V1','V3','V4','V7','V8','V2','V6')]
  994. colnames(mat_chr_reg) <- c('V1','V2','V3','V4','V5','V6','V7')
  995. pat_chr_reg <- read.table("Prejuicer/EXP012.CHiC_Pat_prejuicer_Q30_reg7.medium", header = FALSE, sep = '\t')
  996. pat_chr_reg <- pat_chr_reg[c('V1','V3','V4','V7','V8','V2','V6')]
  997. colnames(pat_chr_reg) <- c('V1','V2','V3','V4','V5','V6','V7')
  998. #==============================
  999. #Kcnk9 TSS
  1000. gene <- "Kcnk9"
  1001. viewpoint <- 72422415
  1002. chrnumb <- 15
  1003. exp <- "mat"
  1004. noexp <- "pat"
  1005. #counting all the contacts from the viewpoint 5kb bin
  1006. mat_chr_plot <- func.paircounts.mat(viewpoint)
  1007. pat_chr_plot <- func.paircounts.pat(viewpoint)
  1008. #which allele is transcriptionally active?
  1009. if (exp == "pat"){
  1010. active <- pat_chr_plot
  1011. inactive <- mat_chr_plot
  1012. } else if (exp == "mat"){
  1013. active <- mat_chr_plot
  1014. inactive <- pat_chr_plot
  1015. }
  1016. #add viewpoint
  1017. active$start <- viewpoint
  1018. inactive$start <- viewpoint
  1019. #generating bed
  1020. activebed <- func.activebed(active)
  1021. inactivebed <- func.inactivebed(inactive)
  1022. activebed$IG <- gene
  1023. inactivebed$IG <- gene
  1024. #cummulative
  1025. cumm.activebed <- rbind(unique(cumm.activebed),activebed)
  1026. cumm.inactivebed <- rbind(unique(cumm.inactivebed),inactivebed)
  1027. #==============================
  1028. #Peg13 TSS
  1029. gene <- "Peg13"
  1030. viewpoint <- 72682173
  1031. chrnumb <- 15
  1032. exp <- "pat"
  1033. noexp <- "mat"
  1034. #counting all the contacts from the viewpoint 5kb bin
  1035. mat_chr_plot <- func.paircounts.mat(viewpoint)
  1036. pat_chr_plot <- func.paircounts.pat(viewpoint)
  1037. #which allele is transcriptionally active?
  1038. if (exp == "pat"){
  1039. active <- pat_chr_plot
  1040. inactive <- mat_chr_plot
  1041. } else if (exp == "mat"){
  1042. active <- mat_chr_plot
  1043. inactive <- pat_chr_plot
  1044. }
  1045. #add viewpoint
  1046. active$start <- viewpoint
  1047. inactive$start <- viewpoint
  1048. #generating bed
  1049. activebed <- func.activebed(active)
  1050. inactivebed <- func.inactivebed(inactive)
  1051. activebed$IG <- gene
  1052. inactivebed$IG <- gene
  1053. #cummulative
  1054. cumm.activebed <- rbind(unique(cumm.activebed),activebed)
  1055. cumm.inactivebed <- rbind(unique(cumm.inactivebed),inactivebed)
  1056. #==============================
  1057. #Trappc9 TSS
  1058. gene <- "Trappc9"
  1059. viewpoint <- 72933053
  1060. chrnumb <- 15
  1061. exp <- "mat"
  1062. noexp <- "pat"
  1063. #counting all the contacts from the viewpoint 5kb bin
  1064. mat_chr_plot <- func.paircounts.mat(viewpoint)
  1065. pat_chr_plot <- func.paircounts.pat(viewpoint)
  1066. #which allele is transcriptionally active?
  1067. if (exp == "pat"){
  1068. active <- pat_chr_plot
  1069. inactive <- mat_chr_plot
  1070. } else if (exp == "mat"){
  1071. active <- mat_chr_plot
  1072. inactive <- pat_chr_plot
  1073. }
  1074. #add viewpoint
  1075. active$start <- viewpoint
  1076. inactive$start <- viewpoint
  1077. #generating bed
  1078. activebed <- func.activebed(active)
  1079. inactivebed <- func.inactivebed(inactive)
  1080. activebed$IG <- gene
  1081. inactivebed$IG <- gene
  1082. #cummulative
  1083. cumm.activebed <- rbind(unique(cumm.activebed),activebed)
  1084. cumm.inactivebed <- rbind(unique(cumm.inactivebed),inactivebed)
  1085. #==============================
  1086. biallelicbed <- read.table("Features/EXP022_CN_BiallelicRNAseq_capture.bed", header=FALSE, sep="\t", stringsAsFactors=FALSE)
  1087. filter_gr <- GRanges(seqnames=15, ranges=IRanges(start=71271582 + 1, end=73517436))
  1088. query_gr <- GRanges(seqnames=biallelicbed[[1]],
  1089. ranges=IRanges(start=biallelicbed[[2]] + 1, end=biallelicbed[[3]]))
  1090. hits <- findOverlaps(query_gr, filter_gr, type = "within")
  1091. filtered_query <- biallelicbed[queryHits(hits), ]
  1092. biallelicregion <- filtered_query
  1093. biallelicregion$viewpoint[biallelicregion$V6 == "-"] <- biallelicregion$V3[biallelicregion$V6 == "-"]
  1094. biallelicregion$viewpoint[biallelicregion$V6 == "+"] <- biallelicregion$V2[biallelicregion$V6 == "+"]
  1095. biallelicregion <- biallelicregion[,c("V4","V1","viewpoint")]
  1096. colnames(biallelicregion) <- c("gene","chrnumb","viewpoint")
  1097. i=1
  1098. for (i in seq_len(nrow(biallelicregion))) {
  1099. gene <- biallelicregion$gene[i]
  1100. viewpoint <- biallelicregion$viewpoint[i]
  1101. chrnumb <- biallelicregion$chrnumb[i]
  1102. exp <- "bi"
  1103. mat_chr_plot <- func.paircounts.mat(viewpoint)
  1104. pat_chr_plot <- func.paircounts.pat(viewpoint)
  1105. if (exp == "bi"){
  1106. biallele1 <- pat_chr_plot
  1107. biallele2 <- mat_chr_plot
  1108. }
  1109. biallele1$start <- viewpoint
  1110. biallele2$start <- viewpoint
  1111. biallele1bed <- func.biallele1bed(biallele1)
  1112. biallele2bed <- func.biallele2bed(biallele2)
  1113. biallele1bed$IG <- gene
  1114. biallele2bed$IG <- gene
  1115. cumm.biallele1bed <- rbind(unique(cumm.biallele1bed),biallele1bed)
  1116. cumm.biallele2bed <- rbind(unique(cumm.biallele2bed),biallele2bed)
  1117. }
  1118. #==============================
  1127. #==============================
  1128. #summary
  1129. cumm.activebed$TSS <- "active allele"
  1130. cumm.inactivebed$TSS <- "inactive allele"
  1131. cumm.biallele1bed$TSS <- "biallele1"
  1132. cumm.biallele2bed$TSS <- "biallele2"
  1133. cumm.all <- rbind(cumm.activebed, cumm.inactivebed,cumm.biallele1bed,cumm.biallele2bed)
  1134. #==============================
  1135. #==============================
  1136. unique(cumm.all$IG)
  1137. summary2 <- cumm.all %>% dplyr::count(TSS)
  1138. summary2$group <- c("mono","bi","bi","mono")
  1139. #==============================
  1140. pdf(file = "EXP012_NumberOfContactsFromTSS_Norm_all_Q30_12_bi-loop.pdf", width = 3.5, height = 3.5)
  1141. p <- ggplot(data=summary2, aes(x=group, y=n, fill=TSS)) +
  1142. geom_bar(stat="identity", position = "fill") +
  1143. scale_y_continuous(expand = c(0, 0)) +
  1144. labs(x = "Imprinted genes", y = "total # contacts from TSS") +
  1145. scale_fill_manual(values=c("chocolate3", "lightblue4","rosybrown3", "gray70")) +
  1146. my_theme + theme(legend.position = "right")
  1147. p
  1148. dev.off()
  1149. #==============================
  1150. cumm.all.filtered3 <- cumm.all %>% filter(TSS %in% c("biallele1","biallele2"))
  1151. summary3 <- cumm.all.filtered3 %>% dplyr::count(IG, TSS)
  1152. #==============================
  1153. pdf(file = "EXP012_NumberOfContactsFromTSS_Norm_individual_fill_Q30_12_bi-loop.pdf", width = 14, height = 3.5)
  1154. p <- ggplot(data=summary3, aes(x=IG, y=n, fill=TSS)) +
  1155. geom_bar(stat="identity", position = "fill") +
  1156. scale_y_continuous(expand = c(0, 0)) +
  1157. labs(x = "Imprinted genes", y = "total # contacts from TSS") +
  1158. scale_fill_manual(values=c("chocolate3", "lightblue4","rosybrown3", "gray70")) +
  1159. my_theme
  1160. p
  1161. dev.off()
  1162. #==============================
  1163. level_order <- c('Sgce','Peg10','Mest','Peg3','Usp29','Snrpn','A230006K03Rik','B230209E15Rik','A330076H08Rik','Ndn','Magel2','Mkrn3','Peg12','Kcnq1ot1','Cdkn1c','Meg3','Peg13',
  1164. 'Copg2','Ube3a','Kcnk9','Trappc9',
  1165. 'Ppp1r9a', 'Casd1', 'Tmem209', 'Zfp583', 'Clcn4','Gabrb3','Ctsd', 'Nap1l4','Cars1')
  1166. cumm.all.filtered <- cumm.all %>% filter(IG %in% level_order)
  1167. cumm.all.filtered <- unique(cumm.all.filtered)
  1168. unique(cumm.all.filtered$IG)
  1169. #==============================
  1170. summary <- cumm.all.filtered %>% dplyr::count(IG, TSS)
  1171. colnames(summary) <- c("IG","TSS","allcounts")
  1172. summary$type[summary$IG %in% c('Sgce','Peg10','Mest','Peg3','Usp29','Snrpn','A230006K03Rik','B230209E15Rik','A330076H08Rik','Ndn','Magel2','Mkrn3','Peg12', 'Kcnq1ot1','Cdkn1c','Meg3','Peg13')] <- "DMR"
  1173. summary$type[summary$IG %in% c('Copg2','Ube3a','Kcnk9','Trappc9')] <- "distal"
  1174. summary$type[summary$IG %in% c('Ppp1r9a', 'Casd1', 'Tmem209', 'Zfp583', 'Clcn4','Gabrb3', 'Ctsd','Nap1l4','Cars1')] <- "nonIG"
  1175. #==============================
  1176. # Basic barplot
  1177. pdf(file = "EXP012_NumberOfContactsFromTSS_Norm_individual_stack_Q30_12_bi-loop.pdf", width = 8, height = 3.5)
  1178. p <- ggplot(data=summary, aes(x=IG, y=allcounts, fill=TSS)) +
  1179. geom_bar(stat="identity") +
  1180. scale_y_continuous(expand = c(0, 0), limit = c(0,3500)) +
  1181. labs(x = "Imprinted genes", y = "total # contacts from TSS") +
  1182. scale_x_discrete(limits = level_order) +
  1183. scale_fill_manual(values=c("chocolate3", "lightblue4","rosybrown3", "gray70")) +
  1184. my_theme
  1185. p
  1186. dev.off()
  1188. pdf(file = "EXP012_NumberOfContactsFromTSS_Norm_individual_fill_Q30_12_bi-loop.pdf", width = 8, height = 3.5)
  1189. p <- ggplot(data=summary, aes(x=IG, y=allcounts, fill=TSS)) +
  1190. geom_bar(stat="identity", position = "fill") +
  1191. scale_y_continuous(expand = c(0, 0)) +
  1192. labs(x = "Imprinted genes", y = "total # contacts from TSS") +
  1193. scale_x_discrete(limits = level_order) +
  1194. scale_fill_manual(values=c("chocolate3", "lightblue4","rosybrown3", "gray70")) +
  1195. my_theme
  1196. p
  1197. dev.off()
  1198. #==============================
  1199. #==============================
  1200. cumm.all.filtered4 <- cumm.all
  1201. cumm.all.filtered4$type <- "nonIG"
  1202. cumm.all.filtered4$type[cumm.all.filtered4$IG %in% c('Sgce','Peg10','Mest','Peg3','Usp29','Snrpn','Peg12', 'Kcnq1ot1','Peg13')] <- "DMR"
  1203. cumm.all.filtered4$type[cumm.all.filtered4$IG %in% c('A230006K03Rik','B230209E15Rik','A330076H08Rik','Ndn','Magel2','Mkrn3','Cdkn1c','Meg3')] <- "DMR"
  1204. cumm.all.filtered4$type[cumm.all.filtered4$IG %in% c('Copg2','Ube3a','Kcnk9','Trappc9')] <- "distal"
  1205. summary4 <- cumm.all.filtered4 %>% dplyr::count(type,TSS)
  1207. pdf(file = "EXP012_NumberOfContactsFromTSS_Norm_groups_stack_Q30_12_bi-loop.pdf", width = 8, height = 3.5)
  1208. p <- ggplot(data=summary4, aes(x=type, y=n, fill=TSS)) +
  1209. geom_bar(stat="identity", position = "fill") +
  1210. scale_y_continuous(expand = c(0, 0)) +
  1211. labs(x = "Imprinted genes", y = "total # contacts from TSS") +
  1212. scale_fill_manual(values=c("chocolate3", "lightblue4","rosybrown3", "gray70")) +
  1213. my_theme + theme(legend.position = "right")
  1214. p
  1215. dev.off()
  1216. #==============================
  1217. #==============================
  1228. ### contact types classified by ChromHMM
  1230. #==============================
  1231. #==============================
  1232. unique(cumm.all$IG)
  1233. cumm.all.filtered <- cumm.all
  1234. cumm.all.filtered <- unique(cumm.all.filtered)
  1235. unique(cumm.all.filtered$IG)
  1236. #==============================
  1237. cumm.all.filtered$active[cumm.all.filtered$IG %in% c('Copg2','Ube3a','Kcnk9','Trappc9','Meg3','Cdkn1c') & cumm.all.filtered$TSS == "active allele"] <- "mat"
  1238. cumm.all.filtered$active[cumm.all.filtered$IG %in% c('Copg2','Ube3a','Kcnk9','Trappc9','Meg3','Cdkn1c') & cumm.all.filtered$TSS == "inactive allele"] <- "pat"
  1239. cumm.all.filtered$active[cumm.all.filtered$IG %in% c('Sgce','Peg10','Mest','Peg3','Usp29','Snrpn','A230006K03Rik','B230209E15Rik','A330076H08Rik','Ndn','Magel2','Mkrn3','Peg12', 'Kcnq1ot1','Peg13') & cumm.all.filtered$TSS == "active allele"] <- "pat"
  1240. cumm.all.filtered$active[cumm.all.filtered$IG %in% c('Sgce','Peg10','Mest','Peg3','Usp29','Snrpn','A230006K03Rik','B230209E15Rik','A330076H08Rik','Ndn','Magel2','Mkrn3','Peg12', 'Kcnq1ot1','Peg13') & cumm.all.filtered$TSS == "inactive allele"] <- "mat"
  1241. cumm.all.filtered$active[cumm.all.filtered$TSS == "biallele1"] <- "mat"
  1242. cumm.all.filtered$active[cumm.all.filtered$TSS == "biallele2"] <- "pat"
  1243. cumm.all.filtered <- unique(cumm.all.filtered)
  1244. cumm.all.filtered$start <- cumm.all.filtered$end
  1245. summary <- cumm.all.filtered %>% dplyr::count(IG, TSS, active)
  1246. colnames(summary) <- c("IG","TSS", "active","allcounts")
  1247. summary$type[summary$IG %in% c('Sgce','Peg10','Mest','Peg3','Usp29','Snrpn','A230006K03Rik','B230209E15Rik','A330076H08Rik','Ndn','Magel2','Mkrn3','Peg12', 'Kcnq1ot1','Cdkn1c','Meg3','Peg13')] <- "DMR"
  1248. summary$type[summary$IG %in% c('Copg2','Ube3a','Kcnk9','Trappc9')] <- "distal"
  1249. summary$type[summary$TSS %in% c("biallele1", "biallele2")] <- "nonIG"
  1250. #==============================
  1251. groups <- c("E1","E2","E3","E4","E5","E6","E7","E8","E9","E10","E11","E12")
  1252. #==============================
  1253. mat.ChromHMM <- read.table("/ChromHMM/Mat_12_segments.bed", header = FALSE, sep = '\t')
  1254. colnames(mat.ChromHMM) <- c("chrom","start","end","group")
  1255. mat.ChromHMM <- mat.ChromHMM[- grep("chr", mat.ChromHMM$chrom),]
  1256. pat.ChromHMM <- read.table("/ChromHMM/Pat_12_segments.bed", header = FALSE, sep = '\t')
  1257. colnames(pat.ChromHMM) <- c("chrom","start","end","group")
  1258. pat.ChromHMM <- pat.ChromHMM[- grep("chr", pat.ChromHMM$chrom),]
  1259. #==============================
  1260. #==============================
  1261. cumm.all.filtered$chrom <- str_replace_all(cumm.all.filtered$chrom,"chr","")
  1262. cumm.all.filtered.mat <- cumm.all.filtered[cumm.all.filtered$active == "mat",]
  1263. cumm.all.filtered.pat <- cumm.all.filtered[cumm.all.filtered$active == "pat",]
  1264. #==============================
  1265. for (group in groups){
  1266. assign(gsub(" ","",paste("summary.mat.counts.",group)), (bed_intersect(cumm.all.filtered.mat, mat.ChromHMM[mat.ChromHMM$group == group,]) %>% count(IG.x, TSS.x)))
  1267. }
  1268. for (group in groups){
  1269. assign(gsub(" ","",paste("summary.pat.counts.",group)), (bed_intersect(cumm.all.filtered.pat, pat.ChromHMM[pat.ChromHMM$group == group,]) %>% count(IG.x, TSS.x)))
  1270. }
  1271. #==============================
  1272. colnames(summary.mat.counts.E1) <- c("IG","TSS","E1")
  1273. colnames(summary.mat.counts.E2) <- c("IG","TSS","E2")
  1274. colnames(summary.mat.counts.E3) <- c("IG","TSS","E3")
  1275. colnames(summary.mat.counts.E4) <- c("IG","TSS","E4")
  1276. colnames(summary.mat.counts.E5) <- c("IG","TSS","E5")
  1277. colnames(summary.mat.counts.E6) <- c("IG","TSS","E6")
  1278. colnames(summary.mat.counts.E7) <- c("IG","TSS","E7")
  1279. colnames(summary.mat.counts.E8) <- c("IG","TSS","E8")
  1280. colnames(summary.mat.counts.E9) <- c("IG","TSS","E9")
  1281. colnames(summary.mat.counts.E10) <- c("IG","TSS","E10")
  1282. colnames(summary.mat.counts.E11) <- c("IG","TSS","E11")
  1283. colnames(summary.mat.counts.E12) <- c("IG","TSS","E12")
  1284. #==============================
  1285. colnames(summary.pat.counts.E1) <- c("IG","TSS","E1")
  1286. colnames(summary.pat.counts.E2) <- c("IG","TSS","E2")
  1287. colnames(summary.pat.counts.E3) <- c("IG","TSS","E3")
  1288. colnames(summary.pat.counts.E4) <- c("IG","TSS","E4")
  1289. colnames(summary.pat.counts.E5) <- c("IG","TSS","E5")
  1290. colnames(summary.pat.counts.E6) <- c("IG","TSS","E6")
  1291. colnames(summary.pat.counts.E7) <- c("IG","TSS","E7")
  1292. colnames(summary.pat.counts.E8) <- c("IG","TSS","E8")
  1293. colnames(summary.pat.counts.E9) <- c("IG","TSS","E9")
  1294. colnames(summary.pat.counts.E10) <- c("IG","TSS","E10")
  1295. colnames(summary.pat.counts.E11) <- c("IG","TSS","E11")
  1296. colnames(summary.pat.counts.E12) <- c("IG","TSS","E12")
  1297. #==============================
  1298. # normalize counts by genome % per each group
  1299. summary.mat.counts.E1$E1 <- summary.mat.counts.E1$E1/0.036
  1300. summary.mat.counts.E2$E2 <- summary.mat.counts.E2$E2/0.317
  1301. summary.mat.counts.E3$E3 <- summary.mat.counts.E3$E3/0.480
  1302. summary.mat.counts.E4$E4 <- summary.mat.counts.E4$E4/0.009
  1303. summary.mat.counts.E5$E5 <- summary.mat.counts.E5$E5/0.038
  1304. summary.mat.counts.E6$E6 <- summary.mat.counts.E6$E6/0.010
  1305. summary.mat.counts.E7$E7 <- summary.mat.counts.E7$E7/0.004
  1306. summary.mat.counts.E8$E8 <- summary.mat.counts.E8$E8/0.010
  1307. summary.mat.counts.E9$E9 <- summary.mat.counts.E9$E9/0.006
  1308. summary.mat.counts.E10$E10 <- summary.mat.counts.E10$E10/0.024
  1309. summary.mat.counts.E11$E11 <- summary.mat.counts.E11$E11/0.015
  1310. summary.mat.counts.E12$E12 <- summary.mat.counts.E12$E12/0.052
  1311. # normalize counts by genome % per each group
  1312. summary.pat.counts.E1$E1 <- summary.pat.counts.E1$E1/0.035
  1313. summary.pat.counts.E2$E2 <- summary.pat.counts.E2$E2/0.310
  1314. summary.pat.counts.E3$E3 <- summary.pat.counts.E3$E3/0.488
  1315. summary.pat.counts.E4$E4 <- summary.pat.counts.E4$E4/0.009
  1316. summary.pat.counts.E5$E5 <- summary.pat.counts.E5$E5/0.038
  1317. summary.pat.counts.E6$E6 <- summary.pat.counts.E6$E6/0.010
  1318. summary.pat.counts.E7$E7 <- summary.pat.counts.E7$E7/0.004
  1319. summary.pat.counts.E8$E8 <- summary.pat.counts.E8$E8/0.010
  1320. summary.pat.counts.E9$E9 <- summary.pat.counts.E9$E9/0.006
  1321. summary.pat.counts.E10$E10 <- summary.pat.counts.E10$E10/0.025
  1322. summary.pat.counts.E11$E11 <- summary.pat.counts.E11$E11/0.015
  1323. summary.pat.counts.E12$E12 <- summary.pat.counts.E12$E12/0.051
  1324. #==============================
  1325. summary.counts.E1 <- rbind(summary.pat.counts.E1, summary.mat.counts.E1)
  1326. summary.counts.E2 <- rbind(summary.pat.counts.E2, summary.mat.counts.E2)
  1327. summary.counts.E3 <- rbind(summary.pat.counts.E3, summary.mat.counts.E3)
  1328. summary.counts.E4 <- rbind(summary.pat.counts.E4, summary.mat.counts.E4)
  1329. summary.counts.E5 <- rbind(summary.pat.counts.E5, summary.mat.counts.E5)
  1330. summary.counts.E6 <- rbind(summary.pat.counts.E6, summary.mat.counts.E6)
  1331. summary.counts.E7 <- rbind(summary.pat.counts.E7, summary.mat.counts.E7)
  1332. summary.counts.E8 <- rbind(summary.pat.counts.E8, summary.mat.counts.E8)
  1333. summary.counts.E9 <- rbind(summary.pat.counts.E9, summary.mat.counts.E9)
  1334. summary.counts.E10 <- rbind(summary.pat.counts.E10, summary.mat.counts.E10)
  1335. summary.counts.E11 <- rbind(summary.pat.counts.E11, summary.mat.counts.E11)
  1336. summary.counts.E12 <- rbind(summary.pat.counts.E12, summary.mat.counts.E12)
  1337. #==============================
  1338. #==============================
  1339. allcounts <- summary
  1340. allcounts <- merge(allcounts, summary.counts.E1, by=c("IG","TSS"),all=TRUE)
  1341. allcounts <- merge(allcounts, summary.counts.E2, by=c("IG","TSS"),all=TRUE)
  1342. allcounts <- merge(allcounts, summary.counts.E3, by=c("IG","TSS"),all=TRUE)
  1343. allcounts <- merge(allcounts, summary.counts.E4, by=c("IG","TSS"),all=TRUE)
  1344. allcounts <- merge(allcounts, summary.counts.E5, by=c("IG","TSS"),all=TRUE)
  1345. allcounts <- merge(allcounts, summary.counts.E6, by=c("IG","TSS"),all=TRUE)
  1346. allcounts <- merge(allcounts, summary.counts.E7, by=c("IG","TSS"),all=TRUE)
  1347. allcounts <- merge(allcounts, summary.counts.E8, by=c("IG","TSS"),all=TRUE)
  1348. allcounts <- merge(allcounts, summary.counts.E9, by=c("IG","TSS"),all=TRUE)
  1349. allcounts <- merge(allcounts, summary.counts.E10, by=c("IG","TSS"),all=TRUE)
  1350. allcounts <- merge(allcounts, summary.counts.E11, by=c("IG","TSS"),all=TRUE)
  1351. allcounts <- merge(allcounts, summary.counts.E12, by=c("IG","TSS"),all=TRUE)
  1352. allcounts[is.na(allcounts)] <- 0
  1353. allcounts <- allcounts[,-4]
  1354. #==============================
  1355. allcounts <- allcounts %>% pivot_longer(cols = 5:16, names_to = "countstype", values_to = "n")
  1356. allcounts$group <- paste(allcounts$IG, allcounts$TSS)
  1357. allcounts$group2 <- paste(allcounts$TSS, allcounts$type)
  1358. sum(allcounts$n)
  1359. levelorder <- c("E2","E3","E8",
  1360. "E1","E4",
  1361. "E9","E11","E12",
  1362. "E5", "E6","E7","E10")
  1363. palette <- c("#000000","#666666","#333333",
  1364. "#993333","#FF6633",
  1365. "#FFFF99","#FFFFCC","#FFFF99",
  1366. "#99CC66","#66CC66","#339966","#336633")
  1367. unique(allcounts$group)
  1368. #==============================
  1369. pdf(file = "EXP012_NumberOfContactsFromTSS_Norm_diffcontacttypes_fill_Q30_12_bi-loop.pdf", width = 6.5, height = 5.5)
  1370. p <- ggplot(data=allcounts, aes(x=factor(group2, level=c("active allele DMR","inactive allele DMR","active allele distal","inactive allele distal","biallele1 nonIG","biallele2 nonIG")), y=n, fill=factor(countstype, level=levelorder))) +
  1371. geom_bar(stat="identity", position = "fill") +
  1372. scale_y_continuous(expand = c(0, 0)) +
  1373. labs(x = "Imprinted genes", y = "Normalized # contacts from TSS") +
  1374. scale_fill_manual(values=palette) +
  1375. my_theme + theme(legend.position = "right", legend.text=element_text(size=8))
  1376. p
  1377. dev.off()
  1378. #==============================
  1379. length(unique(allcounts$IG))
  1380. allcounts$n.ave[allcounts$IG %in% c('Copg2','Ube3a','Kcnk9','Trappc9')] <- 4
  1381. allcounts$n.ave[allcounts$IG %in% c('Sgce','Peg10','Mest','Peg3','Usp29','Snrpn','A230006K03Rik','B230209E15Rik','A330076H08Rik',
  1382. 'Ndn','Magel2','Mkrn3','Peg12', 'Kcnq1ot1','Cdkn1c','Meg3','Peg13')] <- 17
  1383. allcounts$n.ave[allcounts$type == "nonIG"] <- length(unique(allcounts$IG)) - 4 - 17
  1384. allcounts$n <- allcounts$n/allcounts$n.ave
  1385. #==============================
  1386. pdf(file = "EXP012_NumberOfContactsFromTSS_Norm_diffcontacttypes_Q30_12_bi-loop.pdf", width = 6.5, height = 5.5)
  1387. p <- ggplot(data=allcounts, aes(x=factor(group2, level=c("active allele DMR","inactive allele DMR","active allele distal","inactive allele distal","biallele1 nonIG","biallele2 nonIG")), y=n, fill=factor(countstype, level=levelorder))) +
  1388. geom_bar(stat="identity") +
  1389. scale_y_continuous(expand = c(0, 0)) +
  1390. labs(x = "Imprinted genes", y = "Normalized # contacts from TSS") +
  1391. scale_fill_manual(values=palette) +
  1392. my_theme + theme(legend.position = "right", legend.text=element_text(size=8))
  1393. p
  1394. dev.off()
  1395. #==============================

BB_CHiC_GRCm39_contactanalysisChromHMM_FIgure4.R at commit 016d4b2, no license · at the source

Overview

Authors: Bongmin Bae1, Katherine Gu1, Daniel Loftus1, Amanda J Whipple1
ORCID iDs: Bongmin Bae, Katherine Gu, Amanda J Whipple
  1. Department of Molecular and Cellular Biology, Harvard University, Cambridge, MA USA
Institutions: Harvard University (United States)
Journal: Nature communications, volume 17, issue 1, article 9543
Dates: received 30 April 2026; accepted 31 July 2026; published online 8 August 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1038/s41467-026-76506-3 · PMID 42702594 · PMCID PMC13547231 · OpenAlex W7201952788
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: genetics / omics (modality), mouse (organism), cellular / molecular (subfield)
Methods: Spectral & time-frequency, Statistics, Connectivity, Evoked potentials
Keywords: Imprinting, DNA methylation
MeSH: Chromatin*, Enhancer Elements, Genetic*, Genomic Imprinting*, Alleles, Animals, CCCTC-Binding Factor, Cerebral Cortex, DNA Methylation, Female, Male, Mice, Mice, Inbred C57BL, Neurons, Promoter Regions, Genetic, Transcription, Genetic (* major topic)
Topic: Genetic Syndromes and Imprinting (Genetics, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: U.S. Department of Health & Human Services | NIH | National Institute of General Medical Sciences (NIGMS) (F32GM156080, R35GM146921); American Heart Association (American Heart Association, Inc.) (26POST1543266); NEI NIH HHS (P30 EY012196); NIGMS NIH HHS (F32 GM156080, R35 GM146921)
Citations: not cited yet (Europe PMC); 54 references in the paper
Research resources: RRID:SCR_018673

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repositories

Its files are read in the Code ↔ Paper reader above, with 8 matches between paragraphs and lines of code.

Whipple-Lab/capture-hic-imprinting

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 016d4b249d40783a8163af765a2a63d1bfbd1506, 1 June 2026
Languages: R (19), Shell (3)
Size: 23 files, 22 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: data.table (16 files), ggplot2 (11 files), reshape2 (8 files), tidyverse (3 files)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
23 files
At the source: github.com/Whipple-Lab/capture-hic-imprinting

Zenodo 21224497

License: CC-BY-4.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: the references
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: data.table (16 files), ggplot2 (11 files), reshape2 (8 files), tidyverse (3 files)
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
23 files
At the source:

Code availability statement

The paper has a code availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

Read it in the paper: doi.org/10.1038/s41467-026-76506-3.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 44 scripts, each with its path and the digest of its content;
  • 8 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

Data availability statement

The paper has a data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

Read it in the paper: doi.org/10.1038/s41467-026-76506-3.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 4 authors, 2 keywords, 15 MeSH terms, 4 funders, 54 references, 1 RRID.

Cite

This paper

Bae, B., Gu, K., Loftus, D., & Whipple, A. J. (2026). Allele-specific chromatin architecture shapes imprinted domains and coordinates a distal enhancer and antisense transcription at the mouse Mest-Copg2 domain. Nature communications, 17(1), 9543. https://doi.org/10.1038/s41467-026-76506-3

BibTeX

@article{bae2026allele,
author = {Bae, Bongmin and Gu, Katherine and Loftus, Daniel and Whipple, Amanda J},
title = {{Allele-specific chromatin architecture shapes imprinted domains and coordinates a distal enhancer and antisense transcription at the mouse Mest-Copg2 domain}},
journal = {Nature communications},
year = {2026},
month = aug,
volume = {17},
number = {1},
pages = {9543},
publisher = {Nature Publishing Group},
issn = {2041-1723},
doi = {10.1038/s41467-026-76506-3},
url = {https://doi.org/10.1038/s41467-026-76506-3},
pmid = {42702594},
pmcid = {PMC13547231}
}

RIS

TY - JOUR
AU - Bae, Bongmin
AU - Gu, Katherine
AU - Loftus, Daniel
AU - Whipple, Amanda J
TI - Allele-specific chromatin architecture shapes imprinted domains and coordinates a distal enhancer and antisense transcription at the mouse Mest-Copg2 domain
T2 - Nature communications
J2 - Nat Commun
PY - 2026
DA - 2026/08/08
VL - 17
IS - 1
SP - 9543
SN - 2041-1723
PB - Nature Publishing Group
DO - 10.1038/s41467-026-76506-3
UR - https://doi.org/10.1038/s41467-026-76506-3
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41467-026-76506-3",
"type": "article-journal",
"title": "Allele-specific chromatin architecture shapes imprinted domains and coordinates a distal enhancer and antisense transcription at the mouse Mest-Copg2 domain",
"container-title": "Nature communications",
"author": [
{
"family": "Bae",
"given": "Bongmin"
},
{
"family": "Gu",
"given": "Katherine"
},
{
"family": "Loftus",
"given": "Daniel"
},
{
"family": "Whipple",
"given": "Amanda J"
}
],
"container-title-short": "Nat Commun",
"volume": "17",
"issue": "1",
"page": "9543",
"DOI": "10.1038/s41467-026-76506-3",
"PMID": "42702594",
"PMCID": "PMC13547231",
"ISSN": "2041-1723",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41467-026-76506-3",
"language": "en",
"issued": {
"date-parts": [
[
2026,
8,
8
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.21203/rs.3.rs-9927928/v1 [code]
Genome-wide and allele-resolved maps of the radial architecture of the mouse genome
Journal: Research Square (preprint)
In common: reshape2, data.table, ggplot2, 1 other tool, genetics / omics, mouse, 2 references
[2] doi:10.1371/journal.pgen.1012081
ADNP regulates chromatin architecture and lineage fidelity during neural differentiation.
Journal: PLoS genetics
In common: mouse, 5 references
[3] doi:10.1093/bib/bbag096 [code]
scDIAGRAM: detecting chromatin compartments from individual single-cell Hi-C matrix without imputation or reference features.
Journal: Briefings in bioinformatics
In common: mouse, 5 references
[4] doi:10.1038/s41467-026-73325-4 [code]
A scalable Tn5-based method for genome-wide DNA methylation profiling in development and disease.
Journal: Nature communications
In common: reshape2, ggplot2, tidyverse, genetics / omics, 2 references
[5] doi:10.1016/j.stemcr.2026.102930 [code]
ZFHX4 is necessary for dopaminergic neuron differentiation and controls cell cycle by regulating LIN28A.
Journal: Stem cell reports
In common: reshape2, ggplot2, tidyverse, genetics / omics, cellular / molecular, 1 reference
[6] doi:10.1002/imt2.70163 [code]
Spatial multi-omics unveils sphingolipid metabolic reprogramming within the retinal pathological niche.
Journal: iMeta
In common: reshape2, data.table, ggplot2, 1 other tool, genetics / omics, mouse, cellular / molecular
[7] doi:10.1038/s12276-026-01801-4 [code]
Putative glioblastoma origin-like cells in the subventricular zone: isolation and characterization.
Journal: Experimental & molecular medicine
In common: reshape2, data.table, ggplot2, 1 other tool, genetics / omics, mouse, cellular / molecular
[8] doi:10.1371/journal.pcbi.1014573 [code]
Cell-type-specific m1A dynamics are associated with microglial phenotypic transition and neuronal metabolic adaptation during spinal cord injury.
Journal: PLoS computational biology
In common: reshape2, data.table, ggplot2, 1 other tool, genetics / omics, mouse, cellular / molecular
[9] doi:10.1172/jci.insight.207270 [code]
Progressive hypothalamic neuroinflammation in ovariectomized mice parallels aging-related transcriptomic changes in the female human hypothalamus.
Journal: JCI insight
In common: reshape2, data.table, ggplot2, 1 other tool, genetics / omics, mouse, cellular / molecular
[10] doi:10.1186/s13073-026-01704-z [code]
Gene expression profiling enables refined parcellation of cortical layers in the heterogeneous human cerebral cortex.
Journal: Genome medicine
In common: reshape2, data.table, ggplot2, 1 other tool, genetics / omics, mouse, cellular / molecular

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.