Permalink
Cannot retrieve contributors at this time
Name already in use
A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?
BAP_trophoblast_45X/AE_analysis_v01_pass.R
Go to fileThis commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
315 lines (269 sloc)
23.2 KB
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| fx.gg_col <- function(n) { | |
| hues = seq(15, 375, length = n + 1) | |
| hcl(h = hues, l = 65, c = 100)[1:n] | |
| } | |
| # vtf.ESC <- integrate allelic analysis | |
| # ae & df are dataframes, mn are numeric matrices, vc are character vectors, vtf are boolean vectors/matrices, nF and nM are numbers of female and male samples (not doing much with male) | |
| # pulling in the data, adding a pseudocount (pc) of 1 for log-transformation later - automatically subtracted for escapee binomial & ratio | |
| pc<-1 | |
| ae.maleA<-as.data.frame(read.csv("male_tbl_gene_ae_counts_df_aCount.csv",header = T,row.names = 4)) | |
| ae.maleB<-as.data.frame(read.csv("male_tbl_gene_ae_counts_df_bCount.csv",header = T,row.names = 4)) | |
| min.male.AB<-ae.maleA | |
| nM<-length(min.male.AB)-3 | |
| min.male.AB[,c(4:(3+nM))]<-matrix(as.numeric(apply(cbind(unlist(ae.maleA[,c(4:(3+nM))]+pc),unlist(ae.maleB[,c(4:(3+nM))]+pc)), MARGIN = 1, min)), ncol= nM) | |
| ae.male.AB<-ae.maleA | |
| ae.maleA[,c(4:(3+nM))]<-ae.maleA[,c(4:(3+nM))]+pc | |
| ae.maleB[,c(4:(3+nM))]<-ae.maleB[,c(4:(3+nM))]+pc | |
| ae.male.AB[,c(4:(3+nM))]<-ae.maleA[,c(4:(3+nM))]+ae.maleB[,c(4:(3+nM))] | |
| ae.femA<-as.data.frame(read.csv("clean_female_tbl_gene_ae_counts_df_aCount.csv",header = T,row.names = 4)) | |
| ae.femB<-as.data.frame(read.csv("clean_female_tbl_gene_ae_counts_df_bCount.csv",header = T,row.names = 4)) | |
| min.fem.AB<-ae.femA | |
| nF<-length(min.fem.AB)-3 | |
| min.fem.AB[,c(4:(3+nF))]<-matrix(as.numeric(apply(cbind(unlist(ae.femA[,c(4:(3+nF))]+pc),unlist(ae.femB[,c(4:(3+nF))]+pc)), MARGIN = 1, min)), ncol= nF) | |
| ae.fem.AB<-ae.femA | |
| ae.femA[,c(4:(3+nF))]<-ae.femA[,c(4:(3+nF))]+pc | |
| ae.femB[,c(4:(3+nF))]<-ae.femB[,c(4:(3+nF))]+pc | |
| ae.fem.AB[,c(4:(3+nF))]<-ae.femA[,c(4:(3+nF))]+ae.femB[,c(4:(3+nF))] | |
| ae.fem.AB<-ae.fem.AB[order(ae.fem.AB$chr, ae.fem.AB$start),] | |
| min.fem.AB<-min.fem.AB[order(min.fem.AB$chr, min.fem.AB$start),] | |
| #combine replicates into by line | |
| ae.fem.sum<-cbind(rowSums(ae.fem.AB[,which(str_detect(colnames(ae.fem.AB), regex("C6XX", ignore_case = TRUE)))]-pc), | |
| rowSums(ae.fem.AB[,which(str_detect(colnames(ae.fem.AB), regex("C19XX", ignore_case = TRUE)))]-pc), | |
| rowSums(ae.fem.AB[,which(str_detect(colnames(ae.fem.AB), regex("C23XX", ignore_case = TRUE)))]-pc), | |
| rowSums(ae.fem.AB[,which(str_detect(colnames(ae.fem.AB), regex("C2XO", ignore_case = TRUE)))]-pc), | |
| rowSums(ae.fem.AB[,which(str_detect(colnames(ae.fem.AB), regex("C8XO", ignore_case = TRUE)))]-pc), | |
| rowSums(ae.fem.AB[,which(str_detect(colnames(ae.fem.AB), regex("C9XO", ignore_case = TRUE)))])-pc) | |
| colnames(ae.fem.sum)<-c("C6XX", "C19XX", "C23XX", "C2XO", "C8XO", "C9XO") | |
| ae.fem.min<-cbind(rowSums(min.fem.AB[,which(str_detect(colnames(min.fem.AB), regex("C6XX", ignore_case = TRUE)))]-pc), | |
| rowSums(min.fem.AB[,which(str_detect(colnames(min.fem.AB), regex("C19XX", ignore_case = TRUE)))]-pc), | |
| rowSums(min.fem.AB[,which(str_detect(colnames(min.fem.AB), regex("C23XX", ignore_case = TRUE)))]-pc), | |
| rowSums(min.fem.AB[,which(str_detect(colnames(min.fem.AB), regex("C2XO", ignore_case = TRUE)))]-pc), | |
| rowSums(min.fem.AB[,which(str_detect(colnames(min.fem.AB), regex("C8XO", ignore_case = TRUE)))]-pc), | |
| rowSums(min.fem.AB[,which(str_detect(colnames(min.fem.AB), regex("C9XO", ignore_case = TRUE)))]-pc)) | |
| colnames(ae.fem.min)<-c("C6XX", "C19XX", "C23XX", "C2XO", "C8XO", "C9XO") | |
| #individual-sample depth-estimated number of mis-aligned reads based on biallelic rate of XO samples, then summed across replicates (works for chrX only of course) | |
| ae.fem.subXO<-round(ae.fem.sum * rowMeans(ae.fem.min[,which(str_detect(colnames(ae.fem.min), regex("XO", ignore_case = TRUE)))]) | |
| /rowMeans(ae.fem.sum[,which(str_detect(colnames(ae.fem.sum), regex("XO", ignore_case = TRUE)))])) | |
| # subtract from ae.fem.min df or both min and sum dfs | |
| #ae.fem.sum<-ae.fem.sum - ae.fem.subXO | |
| ae.fem.min<-ae.fem.min - ae.fem.subXO | |
| # set negative read counts (rounded) after subtraction to 0 | |
| ae.fem.min[ae.fem.min<0]<-0 | |
| ae.fem.sum[ae.fem.sum<0]<-0 | |
| #Power calculation: http://www.stat.wisc.edu/~st571-1/10-power-4.pdf | |
| sPcut<-0.05 | |
| POW<-0.9 | |
| P0<-0.2 #double of 0.1 cutoff used below, ie. estimating 90% power of identifying all escapees with an allelic ratio of 0.2 | |
| P1<-0.01 # read error estimate | |
| N<-c(1:1000) | |
| K<-qbinom(sPcut, N, P0)-1 | |
| #ALPHA<-pbinom(K, N, P0) | |
| POWER<-pbinom(q = K, size = N,prob =P1) | |
| minREAD<-min(which(POWER>=POW)) | |
| pBi.fem.AB<-matrix(as.numeric(pbinom(q = ae.fem.min, size = ae.fem.sum,prob = (P0/2),lower.tail = F)), ncol=length(colnames(ae.fem.sum))) | |
| pBi.fem.AB[ae.fem.sum<minREAD]<-NA | |
| pBi.fem.AB<-cbind(ae.fem.AB[,c(1:3)], pBi.fem.AB) | |
| colnames(pBi.fem.AB)<-c(colnames(ae.fem.AB[,c(1:3)]), colnames(ae.fem.min)) | |
| ratio.fem.AB<-ae.fem.min/ae.fem.sum | |
| ratio.fem.AB[ae.fem.sum<minREAD]<-NA | |
| ratio.fem.AB<-cbind(ae.fem.AB[,c(1:3)], ratio.fem.AB) | |
| colnames(ratio.fem.AB)<-c(colnames(ae.fem.AB[,c(1:3)]), colnames(ae.fem.min)) | |
| # can choose to omit na, leaving only genes for which calls were made across all XX samples | |
| #pBi.fem.AB<-na.omit(pBi.fem.AB) | |
| #ratio.fem.AB<-na.omit(ratio.fem.AB) | |
| # check pBi and ratio cutoffs | |
| mtf.ratio<-ratio.fem.AB >= P0/2 #0.1 for powercalc 0.2 ratio | |
| mtf.pBi<-pBi.fem.AB<=sPcut | |
| mtf.esc<-mtf.ratio & mtf.pBi | |
| #output escapee gene counts | |
| summary(mtf.esc[which(ratio.fem.AB$chr == "chrX"),]) | |
| # output escapee gene counts after removing rare XO calls | |
| mtf.esc[which(mtf.esc[,7] | mtf.esc[,8] | mtf.esc[,9]),4:9]<-NA | |
| summary(mtf.esc[which(ratio.fem.AB$chr == "chrX"),]) | |
| esc6<-names(which(mtf.esc[which(ratio.fem.AB$chr == "chrX"),4])) | |
| esc19<-names(which(mtf.esc[which(ratio.fem.AB$chr == "chrX"),5])) | |
| esc23<-names(which(mtf.esc[which(ratio.fem.AB$chr == "chrX"),6])) | |
| esc0<-names(which(mtf.esc[which(ratio.fem.AB$chr == "chrX"),7] | mtf.esc[which(ratio.fem.AB$chr == "chrX"),8] | mtf.esc[which(ratio.fem.AB$chr == "chrX"),9])) | |
| vc.esc.all<-names(which(mtf.esc[which(ratio.fem.AB$chr == "chrX"),4] | mtf.esc[which(ratio.fem.AB$chr == "chrX"),5] | mtf.esc[which(ratio.fem.AB$chr == "chrX"),6])) | |
| require(stringr) | |
| mn.ae.fem<-as.matrix(ae.fem.AB[,c(4:(3+nF))]-pc) | |
| mn.min.fem<-as.matrix(min.fem.AB[,c(4:(3+nF))]-pc) | |
| mn.fem.subXO<-round(mn.ae.fem * rowMeans(mn.min.fem[,which(str_detect(colnames(mn.min.fem), regex("XO", ignore_case = TRUE)))]) | |
| /rowMeans(mn.ae.fem[,which(str_detect(colnames(mn.ae.fem), regex("XO", ignore_case = TRUE)))])) | |
| #individual-sample depth-estimated number of mis-aligned reads based on biallelic rate of XO samples (works for chrX only of course) | |
| # subtract from mn.min.fem or both min and ae mns | |
| #mn.ae.fem<-mn.ae.fem - mn.fem.subXO | |
| mn.min.fem<-mn.min.fem - mn.fem.subXO | |
| mn.ae.fem[mn.ae.fem<0]<-0 | |
| mn.min.fem[mn.min.fem<0]<-0 | |
| # R is on a scale from -1 to 1 and captures the alternate Xi in C6 - not used here yet | |
| R.fem.AB<-((ae.femA[,c(4:(nF+3))]-ae.femB[,c(4:(nF+3))])/ae.fem.AB[,c(4:(nF+3))]) | |
| #can multiply by sign of R if want to capture allele identity in allelic ratio analysis | |
| mn.test<-(mn.min.fem/mn.ae.fem) # * sign(R.fem.AB) *sign(rowMeans(R.fem.AB[,c(grep("XO", colnames(R.fem.AB)))])) | |
| # label low counts in individual samples | |
| maskLowAEcounts<-F | |
| if(maskLowAEcounts) {mn.test[which(mn.ae.fem<round(minREAD/4), arr.ind = T)]<-NA} | |
| require(pheatmap) | |
| require(gplots) | |
| mn.test2<-mn.test | |
| rownames(mn.test2)<-t(matrix(unlist(strsplit(rownames(mn.test), "[.]")),nrow=2))[,1] | |
| vc.esc.all2<-t(matrix(unlist(strsplit(vc.esc.all, "[.]")),nrow=2))[,1] | |
| mtf.esc2<-mtf.esc | |
| rownames(mtf.esc2)<-t(matrix(unlist(strsplit(rownames(mtf.esc2), "[.]")),nrow=2))[,1] | |
| vc.ae.Tuk2017<-rownames(mtf.esc2)[which(rownames(mtf.esc2) %in% rownames(df.Tuk2017))] | |
| class(mtf.esc2)<-"character" | |
| mtf.esc2[which(mtf.esc2 == "TRUE", arr.ind = T)]<-"escape" | |
| mtf.esc2[which(mtf.esc2 == "FALSE", arr.ind = T)]<-"inactive" | |
| mtf.esc2[which(mtf.esc2 == "NA", arr.ind = T)]<-"na" | |
| # re-run only from top of block (after powercalc) | |
| colnames(ratio.fem.AB)<-paste("ae.ratio", colnames(ratio.fem.AB), sep = "_") | |
| colnames(pBi.fem.AB)<-paste("ae.pBi", colnames(pBi.fem.AB), sep = "_") | |
| colnames(mtf.esc)<-paste("ae.esc", colnames(mtf.esc), sep = "_") | |
| df.pheatanno_row<-cbind(df.Tuk2017[vc.ae.Tuk2017, c(5:6)],mtf.esc2[vc.ae.Tuk2017,c(4:6)]) | |
| df.pheatanno_col<-data.frame(t(matrix(unlist(strsplit(colnames(mn.test2), "_")),nrow=4))[,3]) | |
| colnames(df.pheatanno_col)<-c("line") | |
| rownames(df.pheatanno_col)<-colnames(mn.test2) | |
| ls.pheatcolor = list(Region = c(PAR = "black", nonPAR = "lightgrey"), | |
| Reported_XCI_status = c(Escape = "yellow", Inactive = "blue", Unknown = "lightgrey", Variable = "lightgreen"), | |
| C6XX = c(escape = "yellow", inactive = "blue", na = "white" ), | |
| C19XX = c(escape = "yellow", inactive = "blue", na = "white" ), | |
| C23XX = c(escape = "yellow", inactive = "blue", na = "white" ), | |
| line = c(C19XX = "darkorange2", C23XX = "indianred", C6XX = "darkolivegreen", C2XO = "grey", C8XO = "grey", C9XO = "grey") | |
| ) | |
| pheatmap(mn.test2[vc.esc.all2,], legend =T, labels_row = df[vc.esc.all2, "gene_name"], | |
| breaks = seq(0,0.5,length.out = 51), colorpanel(50, low = "white", high = "red"), | |
| annotation_row = df.pheatanno_row[vc.esc.all2,], annotation_col = df.pheatanno_col, cluster_rows = F, fontsize_row = 5, fontsize_col = 5, annotation_colors = ls.pheatcolor, | |
| border_color = NA, filename = "allelicHM_escapees.pdf", cellwidth = 5, cellheight = 5) #, width=7.5, height = 10) | |
| pheatmap(mn.test2[rownames(na.omit(df.pheatanno_row)),], legend =T, labels_row = All.df.ENSG[rownames(na.omit(df.pheatanno_row)), "SYMBOL"], | |
| breaks = seq(0,0.5,length.out = 51), colorpanel(50, low = "white", high = "red"), | |
| annotation_row = df.pheatanno_row[rownames(na.omit(df.pheatanno_row)),], annotation_col = df.pheatanno_col, cluster_rows = F , fontsize_row = 4, fontsize_col = 5, annotation_colors = ls.pheatcolor, | |
| border_color = NA, filename = "allelicHM_allX.pdf", width =5, height = 10) #cellwidth = 4, cellheight = 4) | |
| vc.PAIR<-c("ZFX","TXLNG","MXRA5","TBL1X", "USP9X", "KDM5C", "NLGN4X", "KDM6A" , "EIF1AX", "PRKX", "RPS4X", "TMSB4X", "DDX3X") #"TSPYL2") | |
| vc.PAIRY<-c("ZFY","TXLNGY","MXRA5Y","TBL1Y", "USP9Y", "KDM5D", "NLGN4Y", "UTY" , "EIF1AY", "PRKY", "RPS4Y1", "TMSB4Y", "DDX3Y") #"TSPY1") | |
| vtf.PAIR<-(df$gene_name %in% vc.PAIR) | |
| vtf.PAIRY<-(df$gene_name %in% vc.PAIRY) | |
| vtf.ESC<-(rownames(df) %in% vc.esc.all) | |
| vtf.oESC<-(vtf.ESC & !vtf.PAR & !vtf.PAIR) | |
| vtf.FAIL<-((rownames(df) %in% esc19) | (rownames(df) %in% esc23)) & !(rownames(df) %in% esc6) | |
| pdf("Escapee_Counts.pdf", width = 6, height= 6, onefile=T, paper="letter") | |
| par(mfrow=c(2,2), mai = rep(0.6, 4), cex=0.6, cex.lab = 0.7, cex.axis = 0.7, omi = rep(0, 4)) | |
| pch.xo=1 | |
| pch.c6=20 | |
| pch.c19=20 | |
| pch.c23=20 | |
| col.xo=adjustcolor( "grey", alpha.f = 0.7) | |
| col.c6=adjustcolor( "darkolivegreen", alpha.f = 0.7) | |
| col.c19=adjustcolor( "darkorange2", alpha.f = 0.7) | |
| col.c23=adjustcolor( "indianred", alpha.f = 0.7) | |
| fcex=0.5 | |
| pcex=0.8 | |
| col.line=rgb(0,0,0,0.5) | |
| vc.GENES<-df$gene_name[which(vtf.PAIR | vtf.PAR)]#& !(vtf.PAR | vtf.PAIR))] | |
| plot(rowMeans(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),which(str_detect(colnames(ae.femA), regex("XO", ignore_case = TRUE)))]), | |
| rowMeans(ae.femB[rownames(df[which((df$gene_name %in% vc.GENES)),]), which(str_detect(colnames(ae.femB), regex("XO", ignore_case = TRUE)))]), | |
| log="xy",pch=pch.xo, col=col.xo,xlab="mean allele A count", ylab = "mean allele B count", main = "PAR/PAIR genes", cex=pcex) | |
| points(rowMeans(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),which(str_detect(colnames(ae.femA), regex("C19XX", ignore_case = TRUE)))]), | |
| rowMeans(ae.femB[rownames(df[which((df$gene_name %in% vc.GENES)),]), which(str_detect(colnames(ae.femB), regex("C19XX", ignore_case = TRUE)))]), | |
| pch=pch.c19, col=col.c19, cex=pcex) | |
| points(rowMeans(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),which(str_detect(colnames(ae.femA), regex("C23XX", ignore_case = TRUE)))]), | |
| rowMeans(ae.femB[rownames(df[which((df$gene_name %in% vc.GENES)),]), which(str_detect(colnames(ae.femB), regex("C23XX", ignore_case = TRUE)))]), | |
| pch=pch.c23, col=col.c23, cex=pcex) | |
| points(rowMeans(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),which(str_detect(colnames(ae.femA), regex("C6XX", ignore_case = TRUE)))]), | |
| rowMeans(ae.femB[rownames(df[which((df$gene_name %in% vc.GENES)),]), which(str_detect(colnames(ae.femB), regex("C6XX", ignore_case = TRUE)))]), | |
| pch=pch.c6, col=col.c6, cex=pcex) | |
| text(rowMeans(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),which(str_detect(colnames(min.fem.AB), regex("XO", ignore_case = TRUE)))]), | |
| rowMeans(ae.femB[rownames(df[which((df$gene_name %in% vc.GENES)),]), which(str_detect(colnames(ae.fem.AB), regex("XO", ignore_case = TRUE)))]), | |
| labels=df[rownames(na.omit(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),])),"gene_name"],cex=fcex, pos=1, offset = 0.3, col=col.xo) | |
| text(rowMeans(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),which(str_detect(colnames(min.fem.AB), regex("C19XX", ignore_case = TRUE)))]), | |
| rowMeans(ae.femB[rownames(df[which((df$gene_name %in% vc.GENES)),]), which(str_detect(colnames(ae.fem.AB), regex("C19XX", ignore_case = TRUE)))]), | |
| labels=df[rownames(na.omit(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),])),"gene_name"],cex=fcex, pos=1, offset = 0.3, col=col.c19) | |
| text(rowMeans(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),which(str_detect(colnames(min.fem.AB), regex("C23XX", ignore_case = TRUE)))]), | |
| rowMeans(ae.femB[rownames(df[which((df$gene_name %in% vc.GENES)),]), which(str_detect(colnames(ae.fem.AB), regex("C23XX", ignore_case = TRUE)))]), | |
| labels=df[rownames(na.omit(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),])),"gene_name"],cex=fcex, pos=1, offset = 0.3, col=col.c23) | |
| text(rowMeans(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),which(str_detect(colnames(min.fem.AB), regex("C6XX", ignore_case = TRUE)))]), | |
| rowMeans(ae.femB[rownames(df[which((df$gene_name %in% vc.GENES)),]), which(str_detect(colnames(ae.fem.AB), regex("C6XX", ignore_case = TRUE)))]), | |
| labels=df[rownames(na.omit(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),])),"gene_name"],cex=fcex, pos=1, offset = 0.3, col=col.c6) | |
| abline(coef = c(0,1), untf = T, lty=1, col=col.line) | |
| abline(coef = c(0,1/3), untf = T, lty=2, col=col.line) | |
| abline(coef = c(0,3), untf = T, lty=2, col=col.line) | |
| abline(coef = c(0,1/10), untf = T, lty=3, col=col.line) | |
| abline(coef = c(0,9), untf = T, lty=3, col=col.line) | |
| vc.GENES<-df$gene_name[which(vtf.ESC & !(vtf.PAR | vtf.PAIR))] | |
| plot(rowMeans(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),which(str_detect(colnames(ae.femA), regex("XO", ignore_case = TRUE)))]), | |
| rowMeans(ae.femB[rownames(df[which((df$gene_name %in% vc.GENES)),]), which(str_detect(colnames(ae.femB), regex("XO", ignore_case = TRUE)))]), | |
| log="xy",pch=pch.xo, col=col.xo,xlab="mean allele A count", ylab = "mean allele B count", main = "Non-PAR/PAIR escapees", cex=pcex) | |
| points(rowMeans(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),which(str_detect(colnames(ae.femA), regex("C19XX", ignore_case = TRUE)))]), | |
| rowMeans(ae.femB[rownames(df[which((df$gene_name %in% vc.GENES)),]), which(str_detect(colnames(ae.femB), regex("C19XX", ignore_case = TRUE)))]), | |
| pch=pch.c19, col=col.c19, cex=pcex) | |
| points(rowMeans(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),which(str_detect(colnames(ae.femA), regex("C23XX", ignore_case = TRUE)))]), | |
| rowMeans(ae.femB[rownames(df[which((df$gene_name %in% vc.GENES)),]), which(str_detect(colnames(ae.femB), regex("C23XX", ignore_case = TRUE)))]), | |
| pch=pch.c23, col=col.c23, cex=pcex) | |
| points(rowMeans(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),which(str_detect(colnames(ae.femA), regex("C6XX", ignore_case = TRUE)))]), | |
| rowMeans(ae.femB[rownames(df[which((df$gene_name %in% vc.GENES)),]), which(str_detect(colnames(ae.femB), regex("C6XX", ignore_case = TRUE)))]), | |
| pch=pch.c6, col=col.c6, cex=pcex) | |
| text(rowMeans(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),which(str_detect(colnames(min.fem.AB), regex("XO", ignore_case = TRUE)))]), | |
| rowMeans(ae.femB[rownames(df[which((df$gene_name %in% vc.GENES)),]), which(str_detect(colnames(ae.fem.AB), regex("XO", ignore_case = TRUE)))]), | |
| labels=df[rownames(na.omit(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),])),"gene_name"],cex=fcex, pos=1, offset = 0.3, col=col.xo) | |
| text(rowMeans(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),which(str_detect(colnames(min.fem.AB), regex("C19XX", ignore_case = TRUE)))]), | |
| rowMeans(ae.femB[rownames(df[which((df$gene_name %in% vc.GENES)),]), which(str_detect(colnames(ae.fem.AB), regex("C19XX", ignore_case = TRUE)))]), | |
| labels=df[rownames(na.omit(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),])),"gene_name"],cex=fcex, pos=1, offset = 0.3, col=col.c19) | |
| text(rowMeans(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),which(str_detect(colnames(min.fem.AB), regex("C23XX", ignore_case = TRUE)))]), | |
| rowMeans(ae.femB[rownames(df[which((df$gene_name %in% vc.GENES)),]), which(str_detect(colnames(ae.fem.AB), regex("C23XX", ignore_case = TRUE)))]), | |
| labels=df[rownames(na.omit(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),])),"gene_name"],cex=fcex, pos=1, offset = 0.3, col=col.c23) | |
| text(rowMeans(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),which(str_detect(colnames(min.fem.AB), regex("C6XX", ignore_case = TRUE)))]), | |
| rowMeans(ae.femB[rownames(df[which((df$gene_name %in% vc.GENES)),]), which(str_detect(colnames(ae.fem.AB), regex("C6XX", ignore_case = TRUE)))]), | |
| labels=df[rownames(na.omit(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),])),"gene_name"],cex=fcex, pos=1, offset = 0.3, col=col.c6) | |
| abline(coef = c(0,1), untf = T, lty=1, col=col.line) | |
| abline(coef = c(0,1/3), untf = T, lty=2, col=col.line) | |
| abline(coef = c(0,3), untf = T, lty=2, col=col.line) | |
| abline(coef = c(0,1/10), untf = T, lty=3, col=col.line) | |
| abline(coef = c(0,9), untf = T, lty=3, col=col.line) | |
| vc.GENES<-df$gene_name[which(vtf.FAIL)] | |
| plot(rowMeans(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),which(str_detect(colnames(ae.femA), regex("XO", ignore_case = TRUE)))]), | |
| rowMeans(ae.femB[rownames(df[which((df$gene_name %in% vc.GENES)),]), which(str_detect(colnames(ae.femB), regex("XO", ignore_case = TRUE)))]), | |
| log="xy",pch=pch.xo, col=col.xo,xlab="mean allele A count", ylab = "mean allele B count", main = "XX19/23 escapees repressed in XX6", cex=pcex) | |
| points(rowMeans(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),which(str_detect(colnames(ae.femA), regex("C19XX", ignore_case = TRUE)))]), | |
| rowMeans(ae.femB[rownames(df[which((df$gene_name %in% vc.GENES)),]), which(str_detect(colnames(ae.femB), regex("C19XX", ignore_case = TRUE)))]), | |
| pch=pch.c19, col=col.c19, cex=pcex) | |
| points(rowMeans(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),which(str_detect(colnames(ae.femA), regex("C23XX", ignore_case = TRUE)))]), | |
| rowMeans(ae.femB[rownames(df[which((df$gene_name %in% vc.GENES)),]), which(str_detect(colnames(ae.femB), regex("C23XX", ignore_case = TRUE)))]), | |
| pch=pch.c23, col=col.c23, cex=pcex) | |
| points(rowMeans(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),which(str_detect(colnames(ae.femA), regex("C6XX", ignore_case = TRUE)))]), | |
| rowMeans(ae.femB[rownames(df[which((df$gene_name %in% vc.GENES)),]), which(str_detect(colnames(ae.femB), regex("C6XX", ignore_case = TRUE)))]), | |
| pch=pch.c6, col=col.c6, cex=pcex) | |
| text(rowMeans(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),which(str_detect(colnames(min.fem.AB), regex("XO", ignore_case = TRUE)))]), | |
| rowMeans(ae.femB[rownames(df[which((df$gene_name %in% vc.GENES)),]), which(str_detect(colnames(ae.fem.AB), regex("XO", ignore_case = TRUE)))]), | |
| labels=df[rownames(na.omit(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),])),"gene_name"],cex=fcex, pos=1, offset = 0.3, col=col.xo) | |
| text(rowMeans(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),which(str_detect(colnames(min.fem.AB), regex("C19XX", ignore_case = TRUE)))]), | |
| rowMeans(ae.femB[rownames(df[which((df$gene_name %in% vc.GENES)),]), which(str_detect(colnames(ae.fem.AB), regex("C19XX", ignore_case = TRUE)))]), | |
| labels=df[rownames(na.omit(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),])),"gene_name"],cex=fcex, pos=1, offset = 0.3, col=col.c19) | |
| text(rowMeans(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),which(str_detect(colnames(min.fem.AB), regex("C23XX", ignore_case = TRUE)))]), | |
| rowMeans(ae.femB[rownames(df[which((df$gene_name %in% vc.GENES)),]), which(str_detect(colnames(ae.fem.AB), regex("C23XX", ignore_case = TRUE)))]), | |
| labels=df[rownames(na.omit(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),])),"gene_name"],cex=fcex, pos=1, offset = 0.3, col=col.c23) | |
| text(rowMeans(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),which(str_detect(colnames(min.fem.AB), regex("C6XX", ignore_case = TRUE)))]), | |
| rowMeans(ae.femB[rownames(df[which((df$gene_name %in% vc.GENES)),]), which(str_detect(colnames(ae.fem.AB), regex("C6XX", ignore_case = TRUE)))]), | |
| labels=df[rownames(na.omit(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),])),"gene_name"],cex=fcex, pos=1, offset = 0.3, col=col.c6) | |
| abline(coef = c(0,1), untf = T, lty=1, col=col.line) | |
| abline(coef = c(0,1/3), untf = T, lty=2, col=col.line) | |
| abline(coef = c(0,3), untf = T, lty=2, col=col.line) | |
| abline(coef = c(0,1/10), untf = T, lty=3, col=col.line) | |
| abline(coef = c(0,9), untf = T, lty=3, col=col.line) | |
| vc.GENES<-df$gene_name[which((vtf.X | vtf.PAR) & !vtf.ESC)] | |
| plot(rowMeans(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),which(str_detect(colnames(ae.femA), regex("XO", ignore_case = TRUE)))]), | |
| rowMeans(ae.femB[rownames(df[which((df$gene_name %in% vc.GENES)),]), which(str_detect(colnames(ae.femB), regex("XO", ignore_case = TRUE)))]), | |
| log="xy",pch=pch.xo, col=col.xo,xlab="mean allele A count", ylab = "mean allele B count", main = "All silenced genes", cex=pcex) | |
| points(rowMeans(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),which(str_detect(colnames(ae.femA), regex("C19XX", ignore_case = TRUE)))]), | |
| rowMeans(ae.femB[rownames(df[which((df$gene_name %in% vc.GENES)),]), which(str_detect(colnames(ae.femB), regex("C19XX", ignore_case = TRUE)))]), | |
| pch=pch.c19, col=col.c19, cex=pcex) | |
| points(rowMeans(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),which(str_detect(colnames(ae.femA), regex("C23XX", ignore_case = TRUE)))]), | |
| rowMeans(ae.femB[rownames(df[which((df$gene_name %in% vc.GENES)),]), which(str_detect(colnames(ae.femB), regex("C23XX", ignore_case = TRUE)))]), | |
| pch=pch.c23, col=col.c23, cex=pcex) | |
| points(rowMeans(ae.femA[rownames(df[which((df$gene_name %in% vc.GENES)),]),which(str_detect(colnames(ae.femA), regex("C6XX", ignore_case = TRUE)))]), | |
| rowMeans(ae.femB[rownames(df[which((df$gene_name %in% vc.GENES)),]), which(str_detect(colnames(ae.femB), regex("C6XX", ignore_case = TRUE)))]), | |
| pch=pch.c6, col=col.c6, cex=pcex) | |
| abline(coef = c(0,1), untf = T, lty=1, col=col.line) | |
| abline(coef = c(0,1/3), untf = T, lty=2, col=col.line) | |
| abline(coef = c(0,3), untf = T, lty=2, col=col.line) | |
| abline(coef = c(0,1/10), untf = T, lty=3, col=col.line) | |
| abline(coef = c(0,9), untf = T, lty=3, col=col.line) | |
| dev.off() |