\documentclass[11pt, oneside]{article} % use "amsart" instead of "article" for AMSLaTeX format %\usepackage{geometry} % See geometry.pdf to learn the layout options. There are lots. %\geometry{letterpaper} % ... or a4paper or a5paper or ... %\geometry{landscape} % Activate for for rotated page geometry %\usepackage[parfill]{parskip} % Activate to begin paragraphs with an empty line rather than an indent %\usepackage{graphicx} % Use pdf, png, jpg, or eps with pdflatex; use eps in DVI mode % TeX will automatically convert eps --> pdf in pdflatex %\usepackage{amssymb} \usepackage[utf8]{inputenc} %\usepackage[cyr]{aeguill} %\usepackage[francais]{babel} %\usepackage{hyperref} \title{Positive selection on genes interacting with SARS-Cov2, comparison of different analysis} \author{Marie Cariou} \date{February 2021} % Activate to display a given date or no date \begin{document} \maketitle \tableofcontents \newpage \section{Files manipulations} \subsection{Complete table} <<>>= workdir<-"/home/adminmarie/Documents/CIRI_BIBS_projects/2020_05_Etienne_covid/" tab<-read.delim(paste0(workdir, "covid_comp/covid_comp_complete.txt"), h=T, sep="\t") dim(tab) tab$Gene.name<-as.character(tab$Gene.name.x) tab$Gene.name[tab$PreyGene=="MARC1"]<-"MARC1" @ \subsection{Read DGINN Young table} DGINN-Young-primate table correspond to DGINN results, on the SAME alignment as Young-primate. I will merge the 2 tables. <<>>= dginnY<-read.delim(paste0(workdir, "data/summary_primate_young.res"), fill=T, h=T) dim(dginnY) names(dginnY) @ <<>>= add_col<-function(method="PamlM1M2"){ tmp<-dginnY[dginnY$Method==method, c("Gene", "Omega", "PosSel", "PValue", "NbSites", "PSS")] names(tmp)<-c("Gene.name", paste0("Omega_", method), paste0("PosSel_", method), paste0("PValue_", method), paste0("NbSites_", method), paste0("PSS_", method)) tab<-merge(tab, tmp, by="Gene.name") return(tab) } tab<-add_col("PamlM1M2") tab<-add_col("PamlM7M8") tab<-add_col("BppM1M2") tab<-add_col("BppM7M8") # Manip pour la colonne BUSTED tmp<-dginnY[dginnY$Method=="BUSTED",c("Gene", "Omega", "PosSel", "PValue")] names(tmp)<-c("Gene.name", "Omega_BUSTED", "PosSel_BUSTED", "PValue_BUSTED") tab<-merge(tab, tmp, by="Gene.name") tmp<-dginnY[dginnY$Method=="MEME",c("Gene", "NbSites", "PSS")] names(tmp)<-c("Gene.name", "NbSites_MEME", "PSS_MEME") tab<-merge(tab, tmp, by="Gene.name") dim(tab) @ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Comparisons Primates} \subsection{DGINN results on Janet Young's alignments (DGINN-Young-primate) VS Janet Young's results} Comparaison des Omega: colonne L "whole.gene.dN.dS.model.0" VS colonne "omega" dans la sortie de dginn. <<omegaM7M8>>= tab$whole.gene.dN.dS.model.0<-as.numeric(as.character(tab$whole.gene.dN.dS.model.0)) plot(tab$whole.gene.dN.dS.model.0, tab$Omega_PamlM7M8, xlab="Omega Young-primate", ylab="Omega DGINN-Young-primate") abline(0,1) outlier<-tab[tab$whole.gene.dN.dS.model.0<0.2 & tab$Omega_PamlM7M8>0.4,] text(x=outlier$whole.gene.dN.dS.model.0, y=(outlier$Omega_PamlM7M8+0.01), outlier$Gene.name) outlier<-tab[tab$whole.gene.dN.dS.model.0<0.6 & tab$Omega_PamlM7M8>0.7,] text(x=outlier$whole.gene.dN.dS.model.0, y=(outlier$Omega_PamlM7M8+0.01), outlier$Gene.name) @ \subsection{DGINN results on Janet Young's alignments (DGINN-Young-primate) VS DGINN-full's results} Comparaison des Omega: colonne L "whole.gene.dN.dS.model.0" VS colonne "omega" dans la sortie de dginn. <<omegaM7M8_2>>= tab$'dginn.primate_omegaM0Bpp'<-as.numeric(as.character(tab$'dginn.primate_omegaM0Bpp')) plot(tab$'dginn.primate_omegaM0Bpp', tab$Omega_PamlM7M8, xlab="DGINN-full's", ylab="Omega DGINN-Young-primate") abline(0,1) outlier<-tab[tab$'dginn.primate_omegaM0Bpp'>0.4 & tab$Omega_PamlM7M8<0.2,] text(x=outlier$'dginn.primate_omegaM0Bpp', y=(outlier$Omega_PamlM7M8+0.01), outlier$Gene.name) outlier<-tab[tab$'dginn.primate_omegaM0Bpp'>0.5 & tab$Omega_PamlM7M8<0.4,] text(x=outlier$'dginn.primate_omegaM0Bpp', y=(outlier$Omega_PamlM7M8+0.01), outlier$Gene.name) outlier<-tab[tab$'dginn.primate_omegaM0Bpp'>0.2 & tab$Omega_PamlM7M8>0.6,] text(x=outlier$'dginn.primate_omegaM0Bpp', y=(outlier$Omega_PamlM7M8+0.01), outlier$Gene.name) @ \subsection{Janet Young's results (Young-primate) VS DGINN-full's results} Comparaison des Omega: colonne L "whole.gene.dN.dS.model.0" VS colonne "omega" dans la sortie de dginn. <<omegaM7M8_3>>= plot(tab$whole.gene.dN.dS.model.0, as.numeric(as.character(tab$'dginn.primate_omegaM0Bpp')), xlab="Omega Young-primate", ylab="DGINN-full's") abline(0,1) outlier<-tab[tab$whole.gene.dN.dS.model.0<0.4 & as.numeric(as.character(tab$'dginn.primate_omegaM0Bpp'))>0.5,] text(x=outlier$whole.gene.dN.dS.model.0, y=outlier$'dginn.primate_omegaM0Bpp', outlier$Gene.name) outlier<-tab[tab$whole.gene.dN.dS.model.0>0.7 & as.numeric(as.character(tab$'dginn.primate_omegaM0Bpp'))>0,] text(x=outlier$whole.gene.dN.dS.model.0, y=outlier$'dginn.primate_omegaM0Bpp', outlier$Gene.name) outlier<-tab[tab$whole.gene.dN.dS.model.0<0.1 & as.numeric(as.character(tab$'dginn.primate_omegaM0Bpp'))>0.3,] text(x=outlier$whole.gene.dN.dS.model.0+0.03, y=outlier$'dginn.primate_omegaM0Bpp', outlier$Gene.name) @ \section{Overlap} \subsection{Mondrian} <<mondrianprimates>>= library(Mondrian) monddata<-as.data.frame(tab$Gene.name) dim(monddata) dginnyoungtmp<-rowSums(cbind(tab$PosSel_PamlM1M2=="Y", tab$PosSel_PamlM7M8=="Y", tab$PosSel_BppM1M2=="Y", tab$PosSel_BppM7M8=="Y", tab$PosSel_BUSTED=="Y")) #monddata$primates_dginn_young<-ifelse(tmp$PosSel_PamlM7M8=="Y", 1,0) dginnfulltmp<-rowSums(cbind(tab$'dginn.primate_BUSTED'=="Y", tab$'dginn.primate_BppM1M2'=="Y", tab$'dginn.primate_BppM7M8'=="Y", tab$'dginn.primate_codemlM1M2'=="Y", tab$'dginn.primate_codemlM7M8'=="Y")) monddata$primates_young<-ifelse(tab$pVal.M8vsM7<0.05, 1, 0) #monddata$primates_cooper<-ifelse(tab$cooper.primates.M7.M8_p_val<0.05, 1, 0) monddata$primates_dginn_young<-ifelse(dginnyoungtmp>=3, 1,0) monddata$primates_dginn_full<-ifelse(dginnfulltmp>=3, 1,0) mondrian(na.omit(monddata[,2:4]), labels=c("Young", "DGINN-Young >=3", "DGINN-full >=3" )) ##### monddata$primates_dginn_young<-ifelse(dginnyoungtmp>=4, 1,0) monddata$primates_dginn_full<-ifelse(dginnfulltmp>=4, 1,0) mondrian(na.omit(monddata[,2:4]), labels=c("Young", "DGINN-Young >=4", "DGINN-full >=4")) @ Comparison of results with the same method. <<>>= monddata$primates_dginn_young<-tab$PosSel_BppM7M8=="Y" monddata$primates_dginn_full<-tab$'dginn.primate_codemlM7M8'=="Y" mondrian(na.omit(monddata[,2:4]), labels=c("Young", "DGINN-Young", "DGINN-full"), main="posel codeml M7M8") @ \subsection{subsetR} Just another representation of the same result, for now, I focuse on the gene positive in 3 methodes for DGINN analysis. <<subsetprimates>>= library(UpSetR) upsetdata<-as.data.frame(tab$Gene.name) upsetdata$primates_young<-ifelse(tab$pVal.M8vsM7<0.05, 1, 0) ### upsetdata$primates_dginn_young<-ifelse(dginnyoungtmp>=3, 1,0) upsetdata$primates_dginn_full<-ifelse(dginnfulltmp>=3, 1,0) upset(na.omit(upsetdata), nsets = 3, matrix.color = "#DC267F", main.bar.color = "#648FFF", sets.bar.color = "#FE6100") @ \section{Gene List} <<setup, include=FALSE, cache=FALSE, tidy=TRUE>>= options(tidy=TRUE, width=70) @ List of the 34 genes found under positive selection in all analysis. <<>>= upsetdata$`tab$Gene.name`[(upsetdata$primates_young==TRUE & upsetdata$primates_dginn_young==TRUE & upsetdata$primates_dginn_full==TRUE)] @ List of the 13 genes found under positive selection in both Young analysis and DGINN-Young alignments (but not full-DGINN). <<>>= upsetdata$`tab$Gene.name`[(upsetdata$primates_young==TRUE & upsetdata$primates_dginn_young==TRUE & upsetdata$primates_dginn_full==FALSE)] @ List of the 1 gene found under positive selection in both DGINN analysis, but not Young. <<>>= upsetdata$`tab$Gene.name`[(upsetdata$primates_young==FALSE & upsetdata$primates_dginn_young==TRUE & upsetdata$primates_dginn_full==TRUE)] @ List of the 8 genes found under positive selection in both Young analysis and full-DGINN, but not DGINN-young. <<>>= upsetdata$`tab$Gene.name`[(upsetdata$primates_young==TRUE & upsetdata$primates_dginn_young==FALSE & upsetdata$primates_dginn_full==TRUE)] @ List of the 18 genes found under positive selection ONLY in Young analysis. <<>>= upsetdata$`tab$Gene.name`[(upsetdata$primates_young==TRUE & upsetdata$primates_dginn_young==FALSE & upsetdata$primates_dginn_full==FALSE)] @ List of the 1 genes found under positive selection ONLY in DGINN-Young. <<>>= upsetdata$`tab$Gene.name`[(upsetdata$primates_young==FALSE & upsetdata$primates_dginn_young==TRUE & upsetdata$primates_dginn_full==FALSE)] @ List of the 44 genes found under positive selection ONLY in full-DGINN. <<>>= upsetdata$`tab$Gene.name`[(upsetdata$primates_young==FALSE & upsetdata$primates_dginn_young==FALSE & upsetdata$primates_dginn_full==TRUE)] @ <<echo=FALSE, results="hide">>= dginnfulltmp<-rowSums(cbind(tab$'dginn.primate_BUSTED'=="Y", tab$'dginn.primate_BppM1M2'=="Y", tab$'dginn.primate_BppM7M8'=="Y", tab$'dginn.primate_codemlM1M2'=="Y", tab$'dginn.primate_codemlM7M8'=="Y")) tab$Gene.name[dginnfulltmp>=4 & is.na(dginnfulltmp)==F] tab$Gene.name[dginnfulltmp>=3 & is.na(dginnfulltmp)==F] tmp<-tab[dginnfulltmp>=4 & is.na(dginnfulltmp)==F, c("Gene.name","dginn.primate_BUSTED", "dginn.primate_BppM1M2", "dginn.primate_BppM7M8","dginn.primate_codemlM1M2","dginn.primate_codemlM7M8")] write.table(tmp, "geneList_DGINN_full_primate_pos4.txt", row.names=F, quote=F) @ <<shiny, fig.height=11, echo=FALSE, results="hide", fig="hide">>= makeFig1 <- function(df){ # prepare data for colors etc colMethods <- c("deepskyblue4", "darkorange" , "deepskyblue3" , "mediumseagreen" , "yellow3" , "black") nameMethods <- c("BUSTED", "BppM1M2", "BppM7M8", "codemlM1M2", "codemlM7M8", "MEME") metColor <- data.frame(Name = nameMethods , Col = colMethods , stringsAsFactors = FALSE) # subset for this specific figure #df <- df[df$nbY >= 1, ] # to drop genes found by 0 methods (big datasets) xt <- df[, c("BUSTED", "BppM1M2", "BppM7M8", "codemlM1M2", "codemlM7M8")] xt$Gene <- df$Gene nbrMeth <- 5 # reverse order of dataframe so that genes with the most Y are at the bottom (to be on top of the barplot) xt[,1:5] <- ifelse(xt[,1:5] == "Y", 1, 0) # sort and Filter the 0 lines xt<-xt[order(rowSums(xt[,1:5])),] xt<-xt[rowSums(xt[,1:5])>2,] row.names(xt)<-xt$Gene xt<-xt[,1:5] colFig1 <- metColor[which(metColor$Name %in% colnames(xt)) , ] ##### PART 1 : NUMBER OF METHODS par(xpd = NA , mar=c(2,7,4,0) , oma = c(0,0,0,0) , mgp = c(3,0.3,0)) h = barplot( t(xt), border = NA , axes = F , col = adjustcolor(colFig1$Col, alpha.f = 1), horiz = T , las = 2 , main = "Methods detecting positive selection" , cex.main = 0.85, cex.names = min(50/nrow(xt), 1.5) ) axis(3, line = 0, at = c(0:nbrMeth), label = c("0", rep("", nbrMeth -1), nbrMeth), tck = 0.02) legend("bottomleft", horiz = T, border = colFig1$Col, legend = colFig1$Name, fill = colFig1$Col, cex = 0.8, bty = "n", xpd = NA ) } df<-read.delim(paste0(workdir, "/data/DGINN_202005281649summary_cleaned.csv"), fill=T, h=T, sep=",") #makeFig1(df) @ \end{document}