Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • mcariou/2020_dginn_covid19
  • ciri/ps_sars-cov-2/2021_dginn_covid19
2 results
Show changes
Commits on Source (24)
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 1175 additions and 4188 deletions
COVID_PAMLresults_332hits_plusBatScreens_plusDGINN_20200506.txt deleted 100644 → 0
View file @ 77f559a1
source diff could not be displayed: it is too large. Options to address this: view the blob.
COVID_PAMLresults_332hits_plusBatScreens_plusDGINN_20201014.txt deleted 100644 → 0
View file @ 77f559a1
source diff could not be displayed: it is too large. Options to address this: view the blob.
README.md 0 → 100644
View file @ 2f6095f9
# Evolutionary history of SARS-CoV-2 interactome in bats and primates identifies key virus-host interfaces and conflicts
## Introduction
The current COVID-19 pandemic is caused by a novel coronavirus strain, SARS-CoV-2. It originated from the cross-species transmission of a coronavirus from the bat reservoir, directly or through an intermediate host to humans. This catastrophic spillover underlines the necessity to better understand how viruses and hosts have shaped one another over evolutionary time.
Pathogenic viruses put a selective pressure on the host-viral interacting proteins. Identifying which host genes bear signatures of such evolutionary conflict (e.g. positive selection) can lead to the identification of the proteins that have been the most relevant in the response to a virus family. Here, we have used this evolutionary framework to decipher which interactions between the SARS-CoV-2-like viruses and our cells have been important in vivo. In addition, identifying traces of positive selection in different hosts phylogenetic lineages also sheds lights on ancient epidemics and how virus-host determinants may be species specific. This may help to understand differences in susceptibility and pathogenicity to SARS-CoV-like viruses between hosts.
To achieve this, we characterized the evolutionary history of the SARS-CoV-2 interactome identified in in vitro studies: 332 host proteins identified by mass-spectrometry by [Gordon and collaborators](https://www.nature.com/articles/s41586-020-2286-9), as well as two essential SARS-CoV-2 entry factors, the angiotensin converting enzyme 2 (ACE2) and the transmembrane serine protease 2 (TMPRSS2) genes. We characterized their evolution in primates (tracing the human history) and in bats (the natural viral reservoir). To do so, we used [DGINN](https://academic.oup.com/nar/article/48/18/e103/5907962?login=true), a novel computational pipeline to Detect Genetic INNovations in protein-coding genes, which embeds gold-standard methods to perform phylogenetic and positive selection analyses in a high-throughput manner.
## Data formating
Requisite R packages: formatR, tinytex
~
Script to merge DGINN outputs from different batch of analysis and included or correct rows corresponding to genes ran on corrected alignmenents.
```
rnw_scripts/covid_comp_script0_table.pdf
```
Input tables in **data/**.
Output tables in **out_tab/**
The tables output from this script will be used for the following analysis steps.
## Primates Results
**Requisite R packages**: *Mondrian, UpSetR*
~
Script to compare primates screen with Gordon et al.'s positive selection analysis.
```
rnw_scripts/covid_comp_dataset.pdf
```
Main input tables in **out_tab/** and Young's result table in **data/**.
Output tables in **figure/1_xxx**
## Comparison between datasets primates and bats
**Requisite R packages**: *Mondrian, UpSetR, dendextend, ggraph, igraph, tidyverse, viridis.*
~
Script to compare bats and primates screen.
```
rnw_scripts/covid_comp_dataset.pdf
```
Input tables in **out_tab/**.
Output tables in **figure/2_xxx**
## Comparison with MAIC score and pancorona analysis
Script to compare the DGINN screen results to [MAIC](https://www.nature.com/articles/s41598-020-79033-3) score and [pancorona data](https://translational-medicine.biomedcentral.com/articles/10.1186/s12967-020-02480-z).
```
rnw_scripts/covid_comp_maic_pancorona.pdf
```
Input tables in **out_tab/**.
Output tables in **figure/3_xxx**
## GWAS analysis
Get sites under positive selection in FYCO1.
```
rnw_scripts/covid_comp_gwas.pdf
```
Input tables in **out_tab/**.
Output tables in **figure/4_xxx**
covid_comp_bats.Rnw deleted 100644 → 0
View file @ 77f559a1
\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{Mars 2021} % Activate to display a given date or no date
\begin{document}
\maketitle
\tableofcontents
\newpage
\section{Data}
Analysis were formatted by the script covid\_comp\_script0\_table.Rnw.
<<>>=
home<-"/home/adminmarie/Documents/"
workdir<-paste0(home,"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=="MTARC1"]<-"MTARC1"
@
\section{Comparison Bats}
\subsection{Cooper-bats results VS DGINN-bats results}
<<omegaM7M8bats>>=
tab$bats_omegaM0codeml[tab$bats_omegaM0codeml=="na"]<-NA
plot(tab$cooper.batsAverage_dNdS,
as.numeric(as.character(tab$bats_omegaM0codeml)),
xlab="Omega Cooper-bats",
ylab="Omega DGINN-bats")
abline(0,1)
abline(lm(as.numeric(as.character(tab$bats_omegaM0codeml))~
tab$cooper.batsAverage_dNdS),
col="red")
outlier<-tab[tab$cooper.batsAverage_dNdS>0.35 &
as.numeric(as.character(tab$bats_omegaM0codeml))<0.3,]
text(x=outlier$cooper.batsAverage_dNdS,
y=as.numeric(as.character(outlier$bats_omegaM0codeml)),
outlier$Gene.name)
@
\section{Overlap}
\subsection{Data}
<<subbats>>=
tmp<-na.omit(tab[,c("Gene.name", "bats_codemlM7M8_p.value",
"hawkins_Positive.Selection..M8vM8a.p.value",
"cooper.batsM7.M8_p_value", "bats_BUSTED",
"bats_BppM1M2", "bats_BppM7M8", "bats_codemlM1M2",
"bats_codemlM7M8")])
tmp$bats_codemlM7M8_p.value[tmp$bats_codemlM7M8_p.value=="na"]<-NA
tmp$bats_codemlM7M8_p.value<-as.numeric(
as.character(tmp$bats_codemlM7M8_p.value))
dim(tmp)
@
170 genes (present in the 3 experiments)
\subsection{Mondrian}
<<mondrianbats>>=
library(Mondrian)
monddata<-as.data.frame(tmp$Gene.name)
monddata$bats_hawkins<-ifelse(
tmp$hawkins_Positive.Selection..M8vM8a.p.value<0.05, 1, 0)
monddata$bats_cooper<-ifelse(
tmp$cooper.batsM7.M8_p_value<0.05, 1, 0)
dginntmp<-rowSums(cbind(tmp$bats_codemlM1M2=="Y",
tmp$bats_codemlM7M8=="Y",
tmp$bats_BppM1M2=="Y",
tmp$bats_BppM7M8=="Y",
tmp$bats_BUSTED=="Y"))
monddata$bats_dginn<-ifelse(dginntmp>=3, 1,0)
mondrian(monddata[,2:4],
labels=c("DGINN >=3", "hawkins", "Cooper"))
monddata$bats_dginn<-ifelse(dginntmp>=4, 1,0)
mondrian(monddata[,2:4],
labels=c("DGINN >=4", "hawkins", "Cooper"))
@
\subsection{subsetR}
<<subsetbats>>=
library(UpSetR)
upsetdata<-as.data.frame(tmp$Gene.name)
upsetdata$bats_hawkins<-ifelse(
tmp$hawkins_Positive.Selection..M8vM8a.p.value<0.05, 1, 0)
upsetdata$bats_cooper<-ifelse(
tmp$cooper.batsM7.M8_p_value<0.05, 1, 0)
upsetdata$bats_dginn<-ifelse(dginntmp>=3, 1,0)
upset(upsetdata, nsets = 3, matrix.color = "#DC267F",
main.bar.color = "#648FFF", sets.bar.color = "#FE6100")
upsetdata$bats_dginn<-ifelse(dginntmp>=4, 1,0)
upset(upsetdata, nsets = 3, matrix.color = "#DC267F",
main.bar.color = "#648FFF", sets.bar.color = "#FE6100")
@
<<>>=
source("covid_comp_shiny.R")
df<-read.delim(paste0(workdir,
"/data/DGINN_202005281649summary_cleaned.csv"),
fill=T, h=T, sep=",")
names(df)
dftmp<-tab[,c("bats_File", "bats_Name",
"Gene.name", "bats_GeneSize",
"bats_NbSpecies", "bats_omegaM0Bpp",
"bats_omegaM0codeml", "bats_BUSTED",
"bats_BUSTED_p.value", "bats_MEME_NbSites",
"bats_MEME_PSS", "bats_BppM1M2",
"bats_BppM1M2_p.value", "bats_BppM1M2_NbSites",
"bats_BppM1M2_PSS", "bats_BppM7M8",
"bats_BppM7M8_p.value", "bats_BppM7M8_NbSites",
"bats_BppM7M8_PSS", "bats_codemlM1M2",
"bats_codemlM1M2_p.value", "bats_codemlM1M2_NbSites",
"bats_codemlM1M2_PSS", "bats_codemlM7M8",
"bats_codemlM7M8_p.value", "bats_codemlM7M8_NbSites" ,
"bats_codemlM7M8_PSS")]
names(dftmp)<-names(df)
makeFig1(dftmp)
@
\end{document}
covid_comp_bats.pdf deleted 100644 → 0
View file @ 77f559a1
File deleted
covid_comp_complete_old.txt deleted 100644 → 0
View file @ 77f559a1
source diff could not be displayed: it is too large. Options to address this: view the blob.
covid_comp_maic.pdf deleted 100644 → 0
View file @ 77f559a1
File deleted
covid_comp_maic.tex deleted 100644 → 0
View file @ 77f559a1
\documentclass[11pt, oneside]{article}\usepackage[]{graphicx}\usepackage[]{color}
% maxwidth is the original width if it is less than linewidth
% otherwise use linewidth (to make sure the graphics do not exceed the margin)
\makeatletter
\def\maxwidth{ %
\ifdim\Gin@nat@width>\linewidth
\linewidth
\else
\Gin@nat@width
\fi
}
\makeatother
\definecolor{fgcolor}{rgb}{0.345, 0.345, 0.345}
\newcommand{\hlnum}[1]{\textcolor[rgb]{0.686,0.059,0.569}{#1}}%
\newcommand{\hlstr}[1]{\textcolor[rgb]{0.192,0.494,0.8}{#1}}%
\newcommand{\hlcom}[1]{\textcolor[rgb]{0.678,0.584,0.686}{\textit{#1}}}%
\newcommand{\hlopt}[1]{\textcolor[rgb]{0,0,0}{#1}}%
\newcommand{\hlstd}[1]{\textcolor[rgb]{0.345,0.345,0.345}{#1}}%
\newcommand{\hlkwa}[1]{\textcolor[rgb]{0.161,0.373,0.58}{\textbf{#1}}}%
\newcommand{\hlkwb}[1]{\textcolor[rgb]{0.69,0.353,0.396}{#1}}%
\newcommand{\hlkwc}[1]{\textcolor[rgb]{0.333,0.667,0.333}{#1}}%
\newcommand{\hlkwd}[1]{\textcolor[rgb]{0.737,0.353,0.396}{\textbf{#1}}}%
\let\hlipl\hlkwb
\usepackage{framed}
\makeatletter
\newenvironment{kframe}{%
\def\at@end@of@kframe{}%
\ifinner\ifhmode%
\def\at@end@of@kframe{\end{minipage}}%
\begin{minipage}{\columnwidth}%
\fi\fi%
\def\FrameCommand##1{\hskip\@totalleftmargin \hskip-\fboxsep
\colorbox{shadecolor}{##1}\hskip-\fboxsep
% There is no \\@totalrightmargin, so:
\hskip-\linewidth \hskip-\@totalleftmargin \hskip\columnwidth}%
\MakeFramed {\advance\hsize-\width
\@totalleftmargin\z@ \linewidth\hsize
\@setminipage}}%
{\par\unskip\endMakeFramed%
\at@end@of@kframe}
\makeatother
\definecolor{shadecolor}{rgb}{.97, .97, .97}
\definecolor{messagecolor}{rgb}{0, 0, 0}
\definecolor{warningcolor}{rgb}{1, 0, 1}
\definecolor{errorcolor}{rgb}{1, 0, 0}
\newenvironment{knitrout}{}{} % an empty environment to be redefined in TeX
\usepackage{alltt} % 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, maic}
\author{Marie Cariou}
\date{March 2021} % Activate to display a given date or no date
\IfFileExists{upquote.sty}{\usepackage{upquote}}{}
\begin{document}
\maketitle
\tableofcontents
\newpage
\section{Data}
output from covid\_comp\_dataset.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{tablo}\hlkwb{<-}\hlkwd{read.table}\hlstd{(}\hlstr{"primatesVbats.csv"}\hlstd{,}
\hlkwc{h}\hlstd{=T,} \hlkwc{sep}\hlstd{=}\hlstr{","}\hlstd{)}
\end{alltt}
\end{kframe}
\end{knitrout}
Output MAIC formatted by Léa. This table includes the DGINN "score".
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{home}\hlkwb{<-}\hlstr{"/home/adminmarie/Documents/"}
\hlstd{workdir}\hlkwb{<-}\hlkwd{paste0}\hlstd{(home,} \hlstr{"CIRI_BIBS_projects/2020_05_Etienne_covid/"}\hlstd{)}
\hlstd{maic}\hlkwb{<-}\hlkwd{read.table}\hlstd{(}\hlkwd{paste0}\hlstd{(workdir,} \hlstr{"data/covid_comp_maic.txt"}\hlstd{),}
\hlkwc{h}\hlstd{=T)}
\end{alltt}
\end{kframe}
\end{knitrout}
\section{MAIC}
\subsection{Boxplot}
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{par}\hlstd{(}\hlkwc{mfrow}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{2}\hlstd{,}\hlnum{1}\hlstd{))}
\hlkwd{boxplot}\hlstd{(maic}\hlopt{$}\hlstd{rank}\hlopt{~}\hlstd{maic}\hlopt{$}\hlstd{nbats,} \hlkwc{notch}\hlstd{=}\hlnum{TRUE}\hlstd{,} \hlkwc{varwidth}\hlstd{=}\hlnum{TRUE}\hlstd{,} \hlkwc{xlab}\hlstd{=}\hlstr{"score DGINN"}\hlstd{,} \hlkwc{ylab}\hlstd{=}\hlstr{"rank MAIC"}\hlstd{,} \hlkwc{main}\hlstd{=}\hlstr{"Bats"}\hlstd{)}
\end{alltt}
{\ttfamily\noindent\color{warningcolor}{\#\# Warning in bxp(list(stats = structure(c(21, 825, 1664, 2860, 5392, 15, 625.5, : some notches went outside hinges ('box'): maybe set notch=FALSE}}\begin{alltt}
\hlkwd{stripchart}\hlstd{(maic}\hlopt{$}\hlstd{rank}\hlopt{~}\hlstd{maic}\hlopt{$}\hlstd{nbats,} \hlkwc{method}\hlstd{=}\hlstr{"jitter"}\hlstd{,} \hlkwc{vertical}\hlstd{=}\hlnum{TRUE}\hlstd{,} \hlkwc{pch}\hlstd{=}\hlnum{1}\hlstd{,} \hlkwc{cex}\hlstd{=}\hlnum{0.3}\hlstd{,} \hlkwc{add}\hlstd{=}\hlnum{TRUE}\hlstd{)}
\hlkwd{boxplot}\hlstd{(maic}\hlopt{$}\hlstd{rank}\hlopt{~}\hlstd{maic}\hlopt{$}\hlstd{nprimates,} \hlkwc{notch}\hlstd{=}\hlnum{TRUE}\hlstd{,} \hlkwc{xlab}\hlstd{=}\hlstr{"score DGINN"}\hlstd{,} \hlkwc{ylab}\hlstd{=}\hlstr{"rank MAIC"}\hlstd{,} \hlkwc{main}\hlstd{=}\hlstr{"Primates"}\hlstd{)}
\hlkwd{stripchart}\hlstd{(maic}\hlopt{$}\hlstd{rank}\hlopt{~}\hlstd{maic}\hlopt{$}\hlstd{nprimates,} \hlkwc{method}\hlstd{=}\hlstr{"jitter"}\hlstd{,} \hlkwc{vertical}\hlstd{=}\hlnum{TRUE}\hlstd{,} \hlkwc{pch}\hlstd{=}\hlnum{1}\hlstd{,} \hlkwc{cex}\hlstd{=}\hlnum{0.3}\hlstd{,} \hlkwc{add}\hlstd{=}\hlnum{TRUE}\hlstd{)}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/boxplot-1}
\end{knitrout}
\subsection{Dotchart}
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{tmp}\hlkwb{<-}\hlstd{maic[maic}\hlopt{$}\hlstd{nbats}\hlopt{>=}\hlnum{3}\hlstd{,} \hlkwd{c}\hlstd{(}\hlstr{"gene"}\hlstd{,} \hlstr{"rank"}\hlstd{,} \hlstr{"nbats"}\hlstd{)]}
\hlstd{tmp}\hlkwb{<-}\hlstd{tmp[}\hlkwd{order}\hlstd{(tmp}\hlopt{$}\hlstd{rank,} \hlkwc{decreasing} \hlstd{=} \hlnum{TRUE}\hlstd{),]}
\hlstd{tmp}\hlopt{$}\hlstd{col}\hlkwb{<-}\hlstr{"black"}
\hlstd{tmp}\hlopt{$}\hlstd{col[tmp}\hlopt{$}\hlstd{gene}\hlopt{==}\hlstr{"ACE2"}\hlstd{]}\hlkwb{<-}\hlstr{"red"}
\hlstd{tmp}\hlopt{$}\hlstd{col[tmp}\hlopt{$}\hlstd{gene}\hlopt{==}\hlstr{"TMPRSS2"}\hlstd{]}\hlkwb{<-}\hlstr{"red"}
\hlstd{tmp}\hlopt{$}\hlstd{pch[tmp}\hlopt{$}\hlstd{nbats}\hlopt{==}\hlnum{5}\hlstd{]}\hlkwb{<-}\hlnum{1}
\hlstd{tmp}\hlopt{$}\hlstd{pch[tmp}\hlopt{$}\hlstd{nbats}\hlopt{==}\hlnum{4}\hlstd{]}\hlkwb{<-}\hlnum{20}
\hlstd{tmp}\hlopt{$}\hlstd{pch[tmp}\hlopt{$}\hlstd{nbats}\hlopt{==}\hlnum{3}\hlstd{]}\hlkwb{<-}\hlnum{4}
\hlkwd{dotchart}\hlstd{(tmp}\hlopt{$}\hlstd{rank,} \hlkwc{main}\hlstd{=}\hlstr{"Bats DGINN >=3"}\hlstd{,} \hlkwc{xlab}\hlstd{=}\hlstr{"rank MAIC"}\hlstd{,} \hlkwc{labels}\hlstd{=tmp}\hlopt{$}\hlstd{gene,} \hlkwc{pch}\hlstd{=tmp}\hlopt{$}\hlstd{pch,} \hlkwc{col}\hlstd{=tmp}\hlopt{$}\hlstd{col)}
\hlkwd{legend}\hlstd{(}\hlstr{"topright"}\hlstd{,} \hlkwd{c}\hlstd{(}\hlstr{"5 (score DGINN)"}\hlstd{,} \hlstr{"4"}\hlstd{,} \hlstr{"3"}\hlstd{),} \hlkwc{pch}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{1}\hlstd{,}\hlnum{20}\hlstd{,}\hlnum{4}\hlstd{))}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/dotbats-1}
\end{knitrout}
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{tmp}\hlkwb{<-}\hlstd{maic[maic}\hlopt{$}\hlstd{nprimates}\hlopt{>=}\hlnum{3}\hlstd{,} \hlkwd{c}\hlstd{(}\hlstr{"gene"}\hlstd{,} \hlstr{"rank"}\hlstd{,} \hlstr{"nprimates"}\hlstd{)]}
\hlstd{tmp}\hlkwb{<-}\hlstd{tmp[}\hlkwd{order}\hlstd{(tmp}\hlopt{$}\hlstd{rank,} \hlkwc{decreasing} \hlstd{=} \hlnum{TRUE}\hlstd{),]}
\hlstd{tmp}\hlopt{$}\hlstd{pch[tmp}\hlopt{$}\hlstd{nprimates}\hlopt{==}\hlnum{5}\hlstd{]}\hlkwb{<-}\hlnum{1}
\hlstd{tmp}\hlopt{$}\hlstd{pch[tmp}\hlopt{$}\hlstd{nprimates}\hlopt{==}\hlnum{4}\hlstd{]}\hlkwb{<-}\hlnum{20}
\hlstd{tmp}\hlopt{$}\hlstd{pch[tmp}\hlopt{$}\hlstd{nprimates}\hlopt{==}\hlnum{3}\hlstd{]}\hlkwb{<-}\hlnum{4}
\hlstd{tmp}\hlopt{$}\hlstd{col}\hlkwb{<-}\hlstr{"black"}
\hlstd{tmp}\hlopt{$}\hlstd{col[tmp}\hlopt{$}\hlstd{gene}\hlopt{==}\hlstr{"ACE2"}\hlstd{]}\hlkwb{<-}\hlstr{"red"}
\hlstd{tmp}\hlopt{$}\hlstd{col[tmp}\hlopt{$}\hlstd{gene}\hlopt{==}\hlstr{"TMPRSS2"}\hlstd{]}\hlkwb{<-}\hlstr{"red"}
\hlkwd{dotchart}\hlstd{(tmp}\hlopt{$}\hlstd{rank,} \hlkwc{main}\hlstd{=}\hlstr{"Primates DGINN >=3"}\hlstd{,} \hlkwc{xlab}\hlstd{=}\hlstr{"rank MAIC"}\hlstd{,} \hlkwc{labels}\hlstd{=tmp}\hlopt{$}\hlstd{gene,} \hlkwc{pch}\hlstd{=tmp}\hlopt{$}\hlstd{pch,} \hlkwc{cex}\hlstd{=}\hlnum{0.8}\hlstd{,} \hlkwc{col}\hlstd{=tmp}\hlopt{$}\hlstd{col)}
\hlkwd{legend}\hlstd{(}\hlstr{"topright"}\hlstd{,} \hlkwd{c}\hlstd{(}\hlstr{"5 (score DGINN)"}\hlstd{,} \hlstr{"4"}\hlstd{,} \hlstr{"3"}\hlstd{),} \hlkwc{pch}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{1}\hlstd{,}\hlnum{20}\hlstd{,}\hlnum{4}\hlstd{))}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/dotprimates-1}
\end{knitrout}
\section{Pan Corona}
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{pancorona}\hlkwb{<-}\hlkwd{read.table}\hlstd{(}\hlkwd{paste0}\hlstd{(workdir,} \hlstr{"data/pancorona_S5.csv"}\hlstd{),}
\hlkwc{h}\hlstd{=T,} \hlkwc{fill} \hlstd{=} \hlnum{TRUE}\hlstd{,} \hlkwc{sep}\hlstd{=}\hlstr{"\textbackslash{}t"}\hlstd{)}
\hlkwd{names}\hlstd{(pancorona)}\hlkwb{<-}\hlkwd{c}\hlstd{(}\hlstr{"tmp.Gene.name"}\hlstd{,} \hlkwd{names}\hlstd{(pancorona)[}\hlopt{-}\hlnum{1}\hlstd{])}
\hlcom{# Genes en commun}
\hlstd{pancorona}\hlopt{$}\hlstd{tmp.Gene.name[pancorona}\hlopt{$}\hlstd{tmp.Gene.name} \hlopt{%in%} \hlstd{tablo}\hlopt{$}\hlstd{tmp.Gene.name]}
\end{alltt}
\begin{verbatim}
## [1] TBK1 MARK3 GIGYF2 MARK2 G3BP1 LARP1 ACE2 PABPC1
## [9] TMPRSS2 AP3B1 CLCC1 CSDE1 HECTD1 MARK1 MEPCE PDE4DIP
## [17] POR PRKAR2B RAB5C RTN4 SRP54 UBAP2 UBAP2L UBXN8
## [25] SPART BZW2 EIF4E2 SMOC1 STOML2 DDX21 FAM98A G3BP2
## [33] MOV10 PABPC4 UPF1
## 105 Levels: ACE2 ANPEP AP3B1 ATXN2L BTF3 BZW2 CKAP5 CLCC1 ... YTHDF2
\end{verbatim}
\begin{alltt}
\hlcom{# Uniquement dans le tableau pancorona}
\hlkwd{sort}\hlstd{(pancorona}\hlopt{$}\hlstd{tmp.Gene.name[(pancorona}\hlopt{$}\hlstd{tmp.Gene.name} \hlopt{%in%} \hlstd{tablo}\hlopt{$}\hlstd{tmp.Gene.name)}\hlopt{==}\hlnum{FALSE}\hlstd{])}
\end{alltt}
\begin{verbatim}
## [1] ANPEP ATXN2L BTF3 CKAP5 CTSB CTSL
## [7] CYB5R3 DDX1 DDX5 DDX58 DHX9 DNM1L
## [13] EEF1A1 EIF2A EIF3F EIF4B EZR FLNA
## [19] FURIN FUS GSK3A GSK3B HDLBP HNRNPA1
## [25] HNRNPD HNRNPF HNRNPU IFIH1 IGF2BP1 IKBKB
## [31] IKBKE IRF3 ISG15 KPNA3 KPNB1 MYH9
## [37] NCL POLD1 POLR2B PRKRA RBM14 RCHY1
## [43] RPL13A RPL26 RPS13 RPS17 RPS19 RPS9
## [49] SDCBP SERBP1 SGTA SLC1A5 SNAP47 SSB
## [55] STING1 SYNCRIP TANC1 TBCB TMPRSS11D TRAF3
## [61] TUBA4A TUBB2A TUBB4A TUBB6 USP10 VPS36
## [67] XRCC5 XRCC6 YBX1 YTHDF2
## 105 Levels: ACE2 ANPEP AP3B1 ATXN2L BTF3 BZW2 CKAP5 CLCC1 ... YTHDF2
\end{verbatim}
\begin{alltt}
\hlcom{## Uniquement dans tableau }
\hlkwd{sort}\hlstd{(tablo}\hlopt{$}\hlstd{tmp.Gene.name[(tablo}\hlopt{$}\hlstd{tmp.Gene.name} \hlopt{%in%} \hlstd{pancorona}\hlopt{$}\hlstd{tmp.Gene.name)}\hlopt{==}\hlnum{FALSE}\hlstd{])}
\end{alltt}
\begin{verbatim}
## [1] AAR2 AASS AATF ABCC1 ACAD9 ACADM
## [7] ACSL3 ADAM9 ADAMTS1 AGPS AKAP8 AKAP8L
## [13] AKAP9 ALG11 ALG5 ALG8 ANO6 AP2A2
## [19] AP2M1 ARF6 ATE1 ATP13A3 ATP1B1 ATP6AP1
## [25] ATP6V1A BAG5 BCKDK BRD2 BRD4 CCDC86
## [31] CDK5RAP2 CENPF CEP112 CEP135 CEP250 CEP350
## [37] CEP68 CHMP2A CHPF CHPF2 CISD3 CIT
## [43] CLIP4 CNTRL COL6A1 COLGALT1 COMT COQ8B
## [49] CRTC3 CSNK2A2 CSNK2B CUL2 CWC27 CYB5B
## [55] DCAF7 DCAKD DCTPP1 DDX10 DNAJC11 DNAJC19
## [61] DNMT1 DPH5 DPY19L1 ECSIT EDEM3 EIF4H
## [67] ELOC EMC1 ERC1 ERGIC1 ERLEC1 ERMP1
## [73] ERO1B ERP44 ETFA EXOSC2 EXOSC3 EXOSC5
## [79] EXOSC8 F2RL1 FAM162A FAM8A1 FAR2 FASTKD5
## [85] FBLN5 FBN1 FBN2 FBXL12 FKBP10 FKBP15
## [91] FKBP7 FOXRED2 FYCO1 GCC1 GCC2 GDF15
## [97] GFER GGCX GGH GHITM GLA GNB1
## [103] GNG5 GOLGA2 GOLGA3 GOLGA7 GOLGB1 GORASP1
## [109] GPAA1 GPX1 GRIPAP1 GRPEL1 GTF2F2 HDAC2
## [115] HEATR3 HMOX1 HOOK1 HS2ST1 HS6ST2 HSBP1
## [121] HYOU1 IDE IL17RA IMPDH2 INHBE INTS4
## [127] ITGB1 JAKMIP1 KDELC1 KDELC2 LARP4B LARP7
## [133] LMAN2 LOX MAP7D1 MARC1 MAT2B MDN1
## [139] MIB1 MIPOL1 MOGS MPHOSPH10 MRPS2 MRPS25
## [145] MRPS27 MRPS5 MTCH1 MYCBP2 NARS2 NAT14
## [151] NDFIP2 NDUFAF1 NDUFAF2 NDUFB9 NEK9 NEU1
## [157] NGDN NGLY1 NIN NINL NLRX1 NOL10
## [163] NPC2 NPTX1 NSD2 NUP210 NUP214 NUP54
## [169] NUP58 NUP62 NUP88 NUP98 NUTF2 OS9
## [175] PCNT PCSK5 PCSK6 PDZD11 PIGO PIGS
## [181] PITRM1 PKP2 PLAT PLD3 PLEKHA5 PLEKHF2
## [187] PLOD2 PMPCA PMPCB POFUT1 POLA1 POLA2
## [193] PPIL3 PPT1 PRIM1 PRIM2 PRKACA PRKAR2A
## [199] PRRC2B PSMD8 PTBP2 PTGES2 PUSL1 PVR
## [205] QSOX2 RAB10 RAB14 RAB18 RAB1A RAB2A
## [211] RAB7A RAB8A RAE1 RALA RAP1GDS1 RBM28
## [217] RBM41 RBX1 RDX REEP5 REEP6 RETREG3
## [223] RHOA RIPK1 RNF41 RPL36 RRP9 SAAL1
## [229] SBNO1 SCAP SCARB1 SCCPDH SDF2 SELENOS
## [235] SEPSECS SIL1 SIRT5 SLC25A21 SLC27A2 SLC30A6
## [241] SLC30A7 SLC30A9 SLC44A2 SLC9A3R1 SLU7 SNIP1
## [247] SRP19 SRP72 STC2 STOM SUN2 TAPT1
## [253] TARS2 TBCA TBKBP1 TCF12 THTPA TIMM10
## [259] TIMM10B TIMM29 TIMM8B TIMM9 TLE1 TLE3
## [265] TM2D3 TMED5 TMEM39B TMEM97 TOMM70 TOR1A
## [271] TOR1AIP1 TRIM59 TRMT1 TUBGCP2 TUBGCP3 TYSND1
## [277] UGGT2 USP54 VPS11 VPS39 WASHC4 WFS1
## [283] YIF1A ZC3H18 ZC3H7A ZDHHC5 ZNF318 ZNF503
## [289] ZYG11B
## 324 Levels: AAR2 AASS AATF ABCC1 ACAD9 ACADM ACE2 ACSL3 ... ZYG11B
\end{verbatim}
\end{kframe}
\end{knitrout}
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{pancorona}\hlkwb{<-}\hlstd{pancorona[,}\hlkwd{c}\hlstd{(}\hlstr{"tmp.Gene.name"}\hlstd{,} \hlstr{"TOTAL"}\hlstd{)]}
\hlstd{pandginn}\hlkwb{<-}\hlkwd{na.omit}\hlstd{(}\hlkwd{merge}\hlstd{(pancorona, tablo,} \hlkwc{by}\hlstd{=}\hlstr{"tmp.Gene.name"}\hlstd{,} \hlkwc{all.x}\hlstd{=}\hlnum{TRUE}\hlstd{))}
\hlstd{pandginn}\hlkwb{<-}\hlstd{pandginn[}\hlkwd{order}\hlstd{(pandginn}\hlopt{$}\hlstd{nprimates),]}
\hlstd{pandginn}\hlkwb{<-}\hlstd{pandginn[}\hlkwd{order}\hlstd{(pandginn}\hlopt{$}\hlstd{TOTAL),]}
\hlkwd{dotchart}\hlstd{(}\hlkwd{as.matrix}\hlstd{(pandginn[,}\hlnum{2}\hlstd{]),} \hlkwc{labels} \hlstd{= pandginn}\hlopt{$}\hlstd{tmp.Gene.name,} \hlkwc{xlim}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{0}\hlstd{,}\hlnum{5}\hlstd{))}
\hlkwd{points}\hlstd{(pandginn[,}\hlnum{4}\hlstd{],} \hlnum{1}\hlopt{:}\hlkwd{nrow}\hlstd{(pandginn),} \hlkwc{col}\hlstd{=}\hlstr{"blue"}\hlstd{,} \hlkwc{pch}\hlstd{=}\hlnum{20}\hlstd{,} \hlkwc{cex}\hlstd{=}\hlnum{0.7}\hlstd{)}
\hlkwd{points}\hlstd{(pandginn[,}\hlnum{3}\hlstd{],} \hlnum{1}\hlopt{:}\hlkwd{nrow}\hlstd{(pandginn),} \hlkwc{col}\hlstd{=}\hlstr{"blue"}\hlstd{,} \hlkwc{pch}\hlstd{=}\hlnum{4}\hlstd{)}
\hlkwd{legend}\hlstd{(}\hlstr{"bottomright"}\hlstd{,} \hlkwd{c}\hlstd{(}\hlstr{"pancorona score"}\hlstd{,} \hlstr{"dginn primate score"}\hlstd{,} \hlstr{"dginn bats score"}\hlstd{),} \hlkwc{pch}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{1}\hlstd{,}\hlnum{20}\hlstd{,}\hlnum{4}\hlstd{),} \hlkwc{col}\hlstd{=}\hlkwd{c}\hlstd{(}\hlstr{"black"}\hlstd{,} \hlstr{"blue"}\hlstd{,} \hlstr{"blue"}\hlstd{))}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/pancorona-1}
\end{knitrout}
\end{document}
covid_comp_primate-Young-focus.Rnw deleted 100644 → 0
View file @ 77f559a1
\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{March 2021} % Activate to display a given date or no date
\begin{document}
\maketitle
\tableofcontents
\newpage
\section{Files manipulations}
\subsection{Complete table}
<<>>=
home<-"/home/adminmarie/Documents/"
workdir<-paste0(home, "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"))
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_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}
covid_comp_primate-Young-focus.pdf deleted 100644 → 0
View file @ 77f559a1
File deleted
covid_comp_primate-Young-focus_backup.Rnw deleted 100644 → 0
View file @ 77f559a1
\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{October 2020} % Activate to display a given date or no date
\begin{document}
\maketitle
\tableofcontents
\newpage
\section{Files manipulations}
\subsection{Read Janet Young's table}
<<>>=
workdir<-"/home/adminmarie/Documents/CIRI_BIBS_projects/2020_05_Etienne_covid/"
tab<-read.delim(paste0(workdir,
"data/COVID_PAMLresults_332hits_plusBatScreens_2020_Apr14.csv"),
fill=T, h=T, dec=",")
dim(tab)
#names(tab)
@
\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)
@
\subsection{Joining Young and DGINN Young table}
\textit{I hide some code corresponding to verifications of gene names coherence between tables}
<<results="hide", echo=FALSE>>=
head(tab)[,1:5]
# gene avec un nom bizar dans certaines colomne
tab[158,1:10]
#
length(unique(dginnY$Gene))
length(unique(tab$PreyGene))
length(unique(tab$Gene.name))
#quelle paire de colonne contient le plus de noms identiques
sum(unique(dginnY$Gene) %in% unique(tab$PreyGene))
sum(unique(dginnY$Gene) %in% unique(tab$Gene.name))
# dginn$Gene et tab$Gene.name presque identiques sauf 1 ligne.
# Je soupçonne que c'est celle là:
tab[158,1:10]
# Verif:
tab[,1:10][(tab$Gene.name %in% unique(dginnY$Gene))==F,]
# yep
# Remplacement manuel par
as.character(unique(dginnY$Gene)[(unique(dginnY$Gene) %in% tab$Gene.name)==F])
# dans le tableau de Janet
val_remp=as.character(unique(dginnY$Gene)[(unique(dginnY$Gene) %in% tab$Gene.name)==F])
tab$Gene.name<-as.character(tab$Gene.name)
tab$Gene.name[158]<-val_remp
sum(unique(dginnY$Gene) %in% unique(tab$Gene.name))
@
<<>>=
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")
@
\subsection{Read DGINN Table}
<<>>=
dginnT<-read.delim(paste0(workdir,
"/data/DGINN_202005281649summary_cleaned.csv"),
fill=T, h=T, sep=",")
dim(dginnT)
names(dginnT)
# Number of genes in dginn-primate output not present in the original table
dginnT[(dginnT$Gene %in% tab$Gene.name)==F,"Gene"]
# This includes paralogs, recombinations found by DGINN
# and additionnal genes included on purpose
# Number of genes from the original list not present in DGINN output
tab[(tab$Gene.name %in% dginnT$Gene)==F,"Gene.name"]
names(dginnT)<-c("File", "Name", "Gene.name", "GeneSize", "dginn-primate_NbSpecies", "dginn-primate_omegaM0Bpp",
"dginn-primate_omegaM0codeml", "dginn-primate_BUSTED", "dginn-primate_BUSTED.p.value",
"dginn-primate_MEME.NbSites", "dginn-primate_MEME.PSS", "dginn-primate_BppM1M2",
"dginn-primate_BppM1M2.p.value", "dginn-primate_BppM1M2.NbSites", "dginn-primate_BppM1M2.PSS",
"dginn-primate_BppM7M8", "dginn-primate_BppM7M8.p.value", "dginn-primate_BppM7M8.NbSites",
"dginn-primate_BppM7M8.PSS", "dginn-primate_codemlM1M2", "dginn-primate_codemlM1M2.p.value",
"dginn-primate_codemlM1M2.NbSites", "dginn-primate_codemlM1M2.PSS", "dginn-primate_codemlM7M8",
"dginn-primate_codemlM7M8.p.value", "dginn-primate_codemlM7M8.NbSites", "dginn-primate_codemlM7M8.PSS")
@
\subsection{Join Table and DGINN table}
<<>>=
tab<-merge(tab,dginnT, by="Gene.name", all.x=T)
@
\subsection{Write new table}
<<>>=
write.table(tab,
"COVID_PAMLresults_332hits_plusBatScreens_plusDGINN_20201014.txt",
row.names=F, quote=F, sep="\t")
@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\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>>=
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)
@
\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 & tab$'dginn-primate_omegaM0Bpp'>0.5,]
text(x=outlier$whole.gene.dN.dS.model.0,
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.
<<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")
###
upsetdata$primates_dginn_young<-ifelse(dginnyoungtmp>=4, 1,0)
upsetdata$primates_dginn_full<-ifelse(dginnfulltmp>=4, 1,0)
upset(na.omit(upsetdata), nsets = 3, matrix.color = "#DC267F",
main.bar.color = "#648FFF", sets.bar.color = "#FE6100")
@
\section{Gene List}
Genes under positive selection for at least 4 methods.
<<>>=
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)
@
\section{Shiny like}
<<shiny, fig.height=11>>=
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}
covid_comp_primate-Young-focus_backup.pdf deleted 100644 → 0
View file @ 77f559a1
File deleted
covid_comp_primate-Young-focus_backup.tex deleted 100644 → 0
View file @ 77f559a1
This diff is collapsed.
covid_comp_primate.Rnw deleted 100644 → 0
View file @ 77f559a1
\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{October 2020} % Activate to display a given date or no date
\begin{document}
\maketitle
\tableofcontents
\newpage
\section{Data}
Analysis were formatted by the script covid\_comp\_script0\_table.Rnw.
<<>>=
home<-"/home/adminmarie/Documents/"
workdir<-paste0(home, "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=="MTARC1"]<-"MTARC1"
@
\section{Comparisons Primates}
\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_1>>=
tab$dginn.primate_omegaM0Bpp[tab$dginn.primate_omegaM0Bpp=="na"]<-NA
tab$dginn.primate_omegaM0Bpp<-as.numeric(as.character(
tab$dginn.primate_omegaM0Bpp))
plot(tab$whole.gene.dN.dS.model.0,
tab$dginn.primate_omegaM0Bpp,
xlab="Omega Young-primate",
ylab="DGINN-full's",
cex=0.3)
abline(0,1)
abline(lm(tab$dginn.primate_omegaM0Bpp~tab$whole.gene.dN.dS.model.0),
col="red")
outlier<-tab[tab$whole.gene.dN.dS.model.0<0.4 &
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.1 &
tab$dginn.primate_omegaM0Bpp>0.3,]
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.33 &
tab$dginn.primate_omegaM0Bpp<0.2,]
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.6 &
tab$dginn.primate_omegaM0Bpp<0.6,]
text(x=outlier$whole.gene.dN.dS.model.0,
y=outlier$dginn.primate_omegaM0Bpp,
outlier$Gene.name)
@
\subsection{Janet Young's results (Young-primate) VS Cooper's result}
Comparaison des Omega: colonne L "whole.gene.dN.dS.model.0" VS colonne "cooper.primates.Average\_dNdS".
<<omegaM7M8_2>>=
tab$cooper.primates.Average_dNdS<-as.numeric(as.character(
tab$cooper.primates.Average_dNdS))
plot(tab$whole.gene.dN.dS.model.0,
tab$cooper.primates.Average_dNdS,
xlab="Omega Young-primate",
ylab="Omega Cooper-primate",
cex=0.3)
abline(0,1)
abline(lm(tab$cooper.primates.Average_dNdS~tab$whole.gene.dN.dS.model.0),
col="red")
outlier<-tab[tab$whole.gene.dN.dS.model.0<0.15 &
tab$cooper.primates.Average_dNdS>0.4,]
text(x=outlier$whole.gene.dN.dS.model.0,
y=outlier$cooper.primates.Average_dNdS,
outlier$Gene.name, cex=0.5)
outlier<-tab[tab$whole.gene.dN.dS.model.0<0.3 &
tab$cooper.primates.Average_dNdS>0.5,]
text(x=outlier$whole.gene.dN.dS.model.0,
y=outlier$cooper.primates.Average_dNdS,
outlier$Gene.name, cex=0.5)
outlier<-tab[tab$whole.gene.dN.dS.model.0>0.3 &
tab$cooper.primates.Average_dNdS<0.1,]
text(x=outlier$whole.gene.dN.dS.model.0,
y=outlier$cooper.primates.Average_dNdS,
outlier$Gene.name, cex=0.5)
@
\subsection{Cooper's results (Cooper-primate) VS DGINN-full's results}
Comparaison des Omega: colonne "cooper.primates.Average\_dNdS" VS colonne "omega" dans la sortie de dginn.
<<omegaM7M8_3>>=
plot(tab$cooper.primates.Average_dNd,
tab$dginn.primate_omegaM0Bpp,
xlab="Omega Cooper-primate",
ylab="DGINN-full's",
cex=0.3)
abline(0,1)
abline(lm(tab$dginn.primate_omegaM0Bpp~tab$cooper.primates.Average_dNd), col="red")
outlier<-tab[tab$cooper.primates.Average_dNd<0.4 &
tab$dginn.primate_omegaM0Bpp>0.5,]
text(x=outlier$cooper.primates.Average_dNd,
y=outlier$dginn.primate_omegaM0Bpp,
outlier$Gene.name, cex=0.5)
outlier<-tab[tab$cooper.primates.Average_dNd<0.1 &
tab$dginn.primate_omegaM0Bpp>0.3,]
text(x=outlier$cooper.primates.Average_dNd,
y=outlier$dginn.primate_omegaM0Bpp,
outlier$Gene.name, cex=0.5)
outlier<-tab[tab$cooper.primates.Average_dNd>0.7 &
tab$dginn.primate_omegaM0Bpp<0.3,]
text(x=outlier$cooper.primates.Average_dNd,
y=outlier$dginn.primate_omegaM0Bpp,
outlier$Gene.name, cex=0.5)
outlier<-tab[tab$cooper.primates.Average_dNd>0.45 &
tab$dginn.primate_omegaM0Bpp<0.2,]
text(x=outlier$cooper.primates.Average_dNd,
y=outlier$dginn.primate_omegaM0Bpp,
outlier$Gene.name, cex=0.5)
@
\section{Overlap}
\subsection{Mondrian}
<<mondrianprimates>>=
library(Mondrian)
monddata<-as.data.frame(tab$Gene.name)
dim(monddata)
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$primate_cooper<-ifelse(
tab$cooper.primates.M7.M8_p_value<0.05, 1, 0)
monddata$primates_dginn_full<-ifelse(
dginnfulltmp>=3, 1,0)
mondrian(na.omit(monddata[,2:4]),
labels=c("Young", "Cooper", "DGINN-full >=3" ))
monddata$primates_dginn_full<-ifelse(
dginnfulltmp>=4, 1,0)
mondrian(na.omit(monddata[,2:4]),
labels=c("Young", "Cooper", "DGINN-full >=4"))
@
\subsection{subsetR}
Just another representation of the same result.
<<subsetprimates>>=
library(UpSetR)
upsetdata<-as.data.frame(tab$Gene.name)
upsetdata$primates_young<-ifelse(tab$pVal.M8vsM7<0.05, 1, 0)
upsetdata$primate_cooper<-ifelse(
tab$cooper.primates.M7.M8_p_value<0.05, 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")
###
upsetdata$primates_dginn_full<-ifelse(dginnfulltmp>=4, 1,0)
upset(na.omit(upsetdata), nsets = 3, matrix.color = "#DC267F",
main.bar.color = "#648FFF", sets.bar.color = "#FE6100")
@
\section{Gene List}
Genes under positive selection for at least 4 methods.
<<>>=
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)
@
\section{Shiny like}
<<shiny, fig.height=11>>=
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
)
}
@
<<>>=
source("covid_comp_shiny.R")
df<-read.delim(paste0(workdir,
"/data/DGINN_202005281649summary_cleaned.csv"),
fill=T, h=T, sep=",")
names(df)
dftmp<-tab[,c("File", "Name", "Gene.name",
"GeneSize", "dginn.primate_NbSpecies", "dginn.primate_omegaM0Bpp",
"dginn.primate_omegaM0codeml", "dginn.primate_BUSTED", "dginn.primate_BUSTED.p.value",
"dginn.primate_MEME.NbSites", "dginn.primate_MEME.PSS", "dginn.primate_BppM1M2",
"dginn.primate_BppM1M2.p.value", "dginn.primate_BppM1M2.NbSites", "dginn.primate_BppM1M2.PSS",
"dginn.primate_BppM7M8", "dginn.primate_BppM7M8.p.value", "dginn.primate_BppM7M8.NbSites",
"dginn.primate_BppM7M8.PSS", "dginn.primate_codemlM1M2", "dginn.primate_codemlM1M2.p.value",
"dginn.primate_codemlM1M2.NbSites","dginn.primate_codemlM1M2.PSS", "dginn.primate_codemlM7M8",
"dginn.primate_codemlM7M8.p.value", "dginn.primate_codemlM7M8.NbSites" , "dginn.primate_codemlM7M8.PSS")]
names(dftmp)<-names(df)
makeFig1(dftmp)
@
\end{document}
covid_comp_primate.pdf deleted 100644 → 0
View file @ 77f559a1
File deleted
covid_comp_primate.tex deleted 100644 → 0
View file @ 77f559a1
This diff is collapsed.
covid_comp_script0_table_backup.Rnw deleted 100644 → 0
View file @ 77f559a1
\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{Janvier 2021} % Activate to display a given date or no date
\begin{document}
\maketitle
\tableofcontents
\newpage
\section{Files manipulations}
\subsection{Read Janet Young's table}
<<>>=
workdir<-"/home/adminmarie/Documents/CIRI_BIBS_projects/2020_05_Etienne_covid/"
tab<-read.delim(paste0(workdir,
"data/COVID_PAMLresults_332hits_plusBatScreens_2020_Apr14.csv"),
fill=T, h=T, dec=",")
dim(tab)
@
\subsection{Read DGINN Young table}
<<>>=
dginnY<-read.delim(paste0(workdir,
"data/summary_primate_young.res"),
fill=T, h=T)
dim(dginnY)
@
\subsection{Joining Young and DGINN Young table}
<<>>=
# correct gene names (MARC1)
val_remp=as.character(unique(dginnY$Gene)[(unique(dginnY$Gene) %in%
tab$Gene.name)==F])
tab$Gene.name<-as.character(tab$Gene.name)
tab$Gene.name[158]<-val_remp
sum(unique(dginnY$Gene) %in% unique(tab$Gene.name))
@
<<>>=
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")
@
\subsection{Read DGINN Table}
<<>>=
dginnT<-read.delim(paste0(workdir,
"data/DGINN_202005281649summary_cleaned.csv"),
fill=T, h=T, sep=",")
dim(dginnT)
names(dginnT)
# Number of genes in dginn-primate output not present in the original table
dginnT[(dginnT$Gene %in% tab$Gene.name)==F,"Gene"]
# This includes paralogs, recombinations found by DGINN and additionnal genes
# included on purpose
# Number of genes from the original list not present in DGINN output
tab[(tab$Gene.name %in% dginnT$Gene)==F,"Gene.name"]
names(dginnT)<-c("File", "Name", "Gene.name", "GeneSize",
"dginn-primate_NbSpecies", "dginn-primate_omegaM0Bpp",
"dginn-primate_omegaM0codeml", "dginn-primate_BUSTED",
"dginn-primate_BUSTED.p.value", "dginn-primate_MEME.NbSites",
"dginn-primate_MEME.PSS", "dginn-primate_BppM1M2",
"dginn-primate_BppM1M2.p.value", "dginn-primate_BppM1M2.NbSites",
"dginn-primate_BppM1M2.PSS", "dginn-primate_BppM7M8",
"dginn-primate_BppM7M8.p.value", "dginn-primate_BppM7M8.NbSites",
"dginn-primate_BppM7M8.PSS", "dginn-primate_codemlM1M2",
"dginn-primate_codemlM1M2.p.value", "dginn-primate_codemlM1M2.NbSites",
"dginn-primate_codemlM1M2.PSS", "dginn-primate_codemlM7M8",
"dginn-primate_codemlM7M8.p.value", "dginn-primate_codemlM7M8.NbSites",
"dginn-primate_codemlM7M8.PSS")
@
<<eval=FALSE>>=
table(dginnT$`dginn-primate_BUSTED`)
table(dginnT$`dginn-primate_codemlM1M2`)
table(dginnT$`dginn-primate_codemlM7M8`)
table(dginnT$`dginn-primate_BppM1M2`)
table(dginnT$`dginn-primate_BppM7M8`)
table(dginnT$`dginn-primate_BUSTED`=="na",dginnT$`dginn-primate_codemlM1M2`=="na", dginnT$`dginn-primate_codemlM7M8`=="na",
dginnT$`dginn-primate_BppM1M2`=="na", dginnT$`dginn-primate_BppM7M8`=="na" )
@
\subsection{Join Table and DGINN table}
<<>>=
tab<-merge(tab,dginnT, by="Gene.name", all.x=T)
table(tab$`dginn-primate_BUSTED`)
table(tab$`dginn-primate_codemlM1M2`)
table(tab$`dginn-primate_codemlM7M8`)
table(tab$`dginn-primate_BppM1M2`)
table(tab$`dginn-primate_BppM7M8`)
table(tab$`dginn-primate_BUSTED`=="na" | tab$`dginn-primate_codemlM1M2`=="na" | tab$`dginn-primate_codemlM7M8`=="na" |
tab$`dginn-primate_BppM1M2`=="na"| tab$`dginn-primate_BppM7M8`=="na" )
@
\subsection{Add DGINN results on bat dataset}
DGINN results from different analysis.
<<>>=
# original table
dginnbats<-read.delim(paste0(workdir,
"data/DGINN_202005281339summary_cleaned.tab"),
fill=T, h=T)
# rerun on corrected alignment
dginnbatsnew1<-read.delim(paste0(workdir,
"data/DGINN_202011262248_summary.tab"),
fill=T, h=T)
dginnbatsnew2<-read.delim(paste0(workdir,
"data/DGINN_202012192053_summary.tab"),
fill=T, h=T)
# colomne choice, BUSTED and Bppml form first file, codeml from the other one
dginnbatsnew<-dginnbatsnew1
dginnbatsnew$omegaM0codeml<-dginnbatsnew2$omegaM0codeml
dginnbatsnew$codemlM1M2<-dginnbatsnew2$codemlM1M2
dginnbatsnew$codemlM1M2_p.value<-dginnbatsnew2$codemlM1M2_p.value
dginnbatsnew$codemlM1M2_NbSites<-dginnbatsnew2$codemlM1M2_NbSites
dginnbatsnew$codemlM1M2_PSS<-dginnbatsnew2$codemlM1M2_PSS
dginnbatsnew$codemlM7M8<-dginnbatsnew2$codemlM7M8
dginnbatsnew$codemlM7M8_p.value<-dginnbatsnew2$codemlM7M8_p.value
dginnbatsnew$codemlM7M8_NbSites<-dginnbatsnew2$codemlM7M8_NbSites
dginnbatsnew$codemlM7M8_PSS<-dginnbatsnew2$codemlM7M8_PSS
####
## RIPK1 is actually a primat results
## 1. Take it and put it at the right place
ripk1<-as.vector(dginnbatsnew[dginnbatsnew$Gene=="RIPK1",])
tab$`dginn-primate_omegaM0Bpp`<-as.numeric(as.character(tab$`dginn-primate_omegaM0Bpp`))
tab$`dginn-primate_BUSTED.p.value`<-as.numeric(as.character(tab$`dginn-primate_BUSTED.p.value`))
tab$`dginn-primate_BppM1M2.p.value`<-as.numeric(as.character(tab$`dginn-primate_BppM1M2.p.value`))
tab$`dginn-primate_BppM7M8.p.value`<-as.numeric(as.character(tab$`dginn-primate_BppM7M8.p.value`))
tab$`dginn-primate_BppM7M8.PSS`<-as.numeric(as.character(tab$`dginn-primate_BppM7M8.PSS`))
tab$`dginn-primate_codemlM1M2.p.value`<-as.numeric(as.character(tab$`dginn-primate_codemlM1M2.p.value`))
tab$`dginn-primate_codemlM1M2.PSS`<-as.numeric(as.character(tab$`dginn-primate_codemlM1M2.PSS`))
tab$`dginn-primate_codemlM7M8.p.value`<-as.numeric(as.character(tab$`dginn-primate_codemlM7M8.p.value`))
tab$`dginn-primate_codemlM7M8.PSS`<-as.numeric(as.character(tab$`dginn-primate_codemlM7M8.PSS`))
tab[tab$Gene.name=="RIPK1","GeneSize"]<-ripk1$GeneSize
tab[tab$Gene.name=="RIPK1","dginn-primate_NbSpecies"]<-ripk1$NbSpecies
tab[tab$Gene.name=="RIPK1","dginn-primate_omegaM0Bpp"]<-ripk1$omegaM0Bpp
tab[tab$Gene.name=="RIPK1","dginn-primate_omegaM0codeml"]<-ripk1$omegaM0codeml
tab[tab$Gene.name=="RIPK1","dginn-primate_BUSTED"]<-ripk1$BUSTED
tab[tab$Gene.name=="RIPK1","dginn-primate_BUSTED.p.value"]<-ripk1$BUSTED_p.value
tab[tab$Gene.name=="RIPK1","dginn-primate_MEME.NbSites"]<-ripk1$MEME_NbSites
tab[tab$Gene.name=="RIPK1","dginn-primate_MEME.PSS"]<-as.numeric(as.character(ripk1$MEME_PSS))
tab[tab$Gene.name=="RIPK1","dginn-primate_BppM1M2"]<-ripk1$BppM1M2
tab[tab$Gene.name=="RIPK1","dginn-primate_BppM1M2.p.value"]<-ripk1$BppM1M2_p.value
tab[tab$Gene.name=="RIPK1","dginn-primate_BppM1M2.NbSites"]<-ripk1$BppM1M2_NbSites
tab[tab$Gene.name=="RIPK1","dginn-primate_BppM1M2.PSS"]<-ripk1$BppM1M2_PSS
tab[tab$Gene.name=="RIPK1","dginn-primate_BppM7M8"]<-ripk1$BppM7M8
tab[tab$Gene.name=="RIPK1","dginn-primate_BppM7M8.p.value"]<-ripk1$BppM7M8_p.value
tab[tab$Gene.name=="RIPK1","dginn-primate_BppM7M8.NbSites"]<-ripk1$BppM7M8_NbSites
tab[tab$Gene.name=="RIPK1","dginn-primate_BppM7M8.PSS"]<-ripk1$BppM7M8_PSS
tab[tab$Gene.name=="RIPK1","dginn-primate_codemlM1M2"]<-ripk1$codemlM1M2
tab[tab$Gene.name=="RIPK1","dginn-primate_codemlM1M2.p.value"]<-ripk1$codemlM1M2_p.value
tab[tab$Gene.name=="RIPK1","dginn-primate_codemlM1M2.NbSites"]<-ripk1$codemlM1M2_NbSites
tab[tab$Gene.name=="RIPK1","dginn-primate_codemlM1M2.PSS"]<-ripk1$codemlM1M2_PSS
tab[tab$Gene.name=="RIPK1","dginn-primate_codemlM7M8"]<-ripk1$codemlM7M8
tab[tab$Gene.name=="RIPK1","dginn-primate_codemlM7M8.p.value"]<-ripk1$codemlM7M8_p.value
tab[tab$Gene.name=="RIPK1","dginn-primate_codemlM7M8.NbSites"]<-ripk1$codemlM7M8_NbSites
tab[tab$Gene.name=="RIPK1","dginn-primate_codemlM7M8.PSS"]<-ripk1$codemlM7M8_PSS
## 2. Remove it
dginnbatsnew<-dginnbatsnew[dginnbatsnew$Gene!="RIPK1",]
## suppress redundant lines
dginnbats<-dginnbats[(dginnbats$Gene %in% dginnbatsnew$Gene)==FALSE,]
names(dginnbatsnew)<-names(dginnbats)
##############"
dginnbatsnew[,4]<-as.numeric(dginnbatsnew[,4])
dginnbats[,6]<-as.numeric(as.character(dginnbats[,6]))
dginnbats[,8]<-as.character(dginnbats[,8])
dginnbats[,12]<-as.character(dginnbats[,12])
dginnbats[,13]<-as.numeric(as.character(dginnbats[,13]))
dginnbats[,16]<-as.character(dginnbats[,16])
dginnbats[,17]<-as.numeric(as.character(dginnbats[,17]))
## replace by new data
dginnbats<-rbind(dginnbats, dginnbatsnew)
names(dginnbats)<-c("File", "bats_Name", "cooper.batsGene", paste0("bats_",
names(dginnbats)[-(1:3)]))
names(dginnbats)
tab<-merge(tab,dginnbats, by="cooper.batsGene", all.x=T)
@
\subsection{Write the new table}
<<>>=
write.table(tab, "covid_comp_complete_old.txt", row.names=FALSE, quote=FALSE, sep="\t")
@
\section{Second Table}
Table containing the DGINN results for both Primates and bats. Conserve all genes.
\subsection{Primates}
<<>>=
dginnT<-read.delim(paste0(workdir,
"data/DGINN_202005281649summary_cleaned.csv"),
fill=T, h=T, sep=",")
dim(dginnT)
names(dginnT)
# Rename the columns to include primate
names(dginnT)<-c("File", "Name", "Gene.name", "GeneSize",
"dginn-primate_NbSpecies", "dginn-primate_omegaM0Bpp",
"dginn-primate_omegaM0codeml", "dginn-primate_BUSTED",
"dginn-primate_BUSTED.p.value", "dginn-primate_MEME.NbSites",
"dginn-primate_MEME.PSS", "dginn-primate_BppM1M2",
"dginn-primate_BppM1M2.p.value", "dginn-primate_BppM1M2.NbSites",
"dginn-primate_BppM1M2.PSS", "dginn-primate_BppM7M8",
"dginn-primate_BppM7M8.p.value", "dginn-primate_BppM7M8.NbSites",
"dginn-primate_BppM7M8.PSS", "dginn-primate_codemlM1M2",
"dginn-primate_codemlM1M2.p.value", "dginn-primate_codemlM1M2.NbSites",
"dginn-primate_codemlM1M2.PSS", "dginn-primate_codemlM7M8",
"dginn-primate_codemlM7M8.p.value", "dginn-primate_codemlM7M8.NbSites",
"dginn-primate_codemlM7M8.PSS")
@
\subsection{Bats}
<<>>=
# original table
dginnbats<-read.delim(paste0(workdir,
"data/DGINN_202005281339summary_cleaned-LE201108.txt"),
fill=T, h=T)
# rerun on corrected alignment
dginnbatsnew<-read.delim(paste0(workdir,
"data/DGINN_202011262248_hyphybpp-202012192053_codeml-summary.txt"),
fill=T, h=T)
@
<<>>=
# Add both columns
dginnbatsnew$Lucie.s.comments<-""
dginnbatsnew$Action.taken<-""
# Homogenize column names
dginnbats$BUSTED_p.value<-dginnbats$BUSTED.p.value
dginnbats$MEME_NbSites<-dginnbats$MEME.NbSites
dginnbats$MEME_PSS<-dginnbats$MEME.PSS
dginnbats$BppM1M2_p.value<-dginnbats$BppM1M2.p.value
dginnbats$BppM1M2_NbSites<-dginnbats$BppM1M2.NbSites
dginnbats$BppM1M2_PSS<-dginnbats$BppM1M2.PSS
dginnbats$BppM7M8_p.value<-dginnbats$BppM7M8.p.value
dginnbats$BppM7M8_NbSites<-dginnbats$BppM7M8.NbSites
dginnbats$BppM7M8_PSS<-dginnbats$BppM7M8.PSS
dginnbats$codemlM1M2_p.value<-dginnbats$codemlM1M2.p.value
dginnbats$codemlM1M2_NbSites<-dginnbats$codemlM1M2.NbSites
dginnbats$codemlM1M2_PSS<-dginnbats$codemlM1M2.PSS
dginnbats$codemlM7M8_p.value<-dginnbats$codemlM7M8.p.value
dginnbats$codemlM7M8_NbSites<-dginnbats$codemlM7M8.NbSites
dginnbats$codemlM7M8_PSS<-dginnbats$codemlM7M8.PSS
@
<<>>=
# Order columns in the same order in both tables
dginnbats<-dginnbats[,names(dginnbatsnew)]
names(dginnbatsnew) %in% names(dginnbats)
names(dginnbats)==names(dginnbatsnew)
# Put RIPK aside
ripk1<-dginnbatsnew[dginnbatsnew$Gene=="RIPK1",1:27]
# Add it to primate table
names(ripk1)<-names(dginnT)
ripk1$`dginn-primate_omegaM0Bpp`<-as.factor(ripk1$`dginn-primate_omegaM0Bpp`)
ripk1$`dginn-primate_BUSTED.p.value`<-as.factor(ripk1$`dginn-primate_BUSTED.p.value`)
ripk1$`dginn-primate_BppM1M2.p.value`<-as.factor(ripk1$`dginn-primate_BppM1M2.p.value`)
ripk1$`dginn-primate_BppM7M8.p.value`<-as.factor(ripk1$`dginn-primate_BppM7M8.p.value`)
dginnT<-rbind(dginnT, ripk1)
## Remove it Ripk1 from bats
dginnbatsnew<-dginnbatsnew[dginnbatsnew$Gene!="RIPK1",]
## suppress redundant lines
dginnbats<-dginnbats[(dginnbats$Gene %in% dginnbatsnew$Gene)==FALSE,]
names(dginnbatsnew)<-names(dginnbats)
## replace by new data
dginnbatsnew$omegaM0Bpp<-as.factor(dginnbatsnew$omegaM0Bpp)
dginnbatsnew$BppM1M2_p.value<-as.factor(dginnbatsnew$BppM1M2_p.value)
dginnbatsnew$BppM7M8_p.value<-as.factor(dginnbatsnew$BppM7M8_p.value)
dginnbats<-rbind(dginnbats, dginnbatsnew)
names(dginnbats)<-c("bats_File", "bats_Name", "Gene.name", paste0("bats_",
names(dginnbats)[-(1:3)]))
names(dginnbats)
@
\subsection{Merged table}
<<setup, include=FALSE, cache=FALSE, tidy=TRUE>>=
options(tidy=TRUE, width=70)
@
<<>>=
#tidy.opts = list(width.cutoff = 60)
dim(dginnT)
dginnT$Gene.name
dim(dginnbats)
dginnbats$Gene.name
@
Manual corrections:
TMPRSS2 in bats
<<>>=
dginnbats[dginnbats$Gene.name=="TMPRSS2",]
# keeping the uncut one
# renaming the other one TMPRSS2_cut
dginnbats$Gene.name<-as.character(dginnbats$Gene.name)
dginnbats[dginnbats$bats_File=="TMPRSS2_bat_select_cut_mafft_prank","Gene.name"]<-"TMPRSS2_cut"
@
RIPK1: ANcestral version kept, suppress it "RIPK1\_sequences\_filtered\_longestORFs\_mafft\_mincov\_prank"
<<>>=
dginnT<-dginnT[dginnT$File!="RIPK1_sequences_filtered_longestORFs_mafft_mincov_prank",]
@
REEP6 eA et B
<<>>=
dginnbats$Gene.name<-as.character(dginnbats$Gene.name)
dginnbats[dginnbats$bats_File=="REEP6_sequences_filtered_longestORFs_D210gp1_prank", "Gene.name"]<-"REEP6_old"
dginnbats[dginnbats$bats_File=="REEP6_LA_bat_select_mafft_prank", "Gene.name"]<-"REEP6"
dginnbats[dginnbats$bats_File=="REEP6_LB_bat_select_mafft_prank", "Gene.name"]<-"REEP6_like"
@
GNG5
<<>>=
dginnT$Gene.name<-as.character(dginnT$Gene.name)
dginnT[dginnT$File=="GNG5_sequences_filtered_longestORFs_D189gp2_prank", "Gene.name"]<-"GNG5_like"
@
<<>>=
dim(dginnbats)
dim(dginnT)
# genes in common
common<-dginnT$Gene.name[dginnT$Gene.name %in% dginnbats$Gene.name]
common
length(dginnT$Gene.name[dginnT$Gene.name %in% dginnbats$Gene.name])
# genes only in primates
onlyprimates<-dginnT$Gene.name[(dginnT$Gene.name %in% dginnbats$Gene.name)==FALSE]
onlyprimates
length(dginnT$Gene.name[(dginnT$Gene.name %in% dginnbats$Gene.name)==FALSE])
# genes only in bats
onlybats<-dginnbats$Gene.name[(dginnbats$Gene.name %in% dginnT$Gene.name)==FALSE]
onlybats
length(dginnbats$Gene.name[(dginnbats$Gene.name %in% dginnT$Gene.name)==FALSE])
@
<<>>=
tab<-merge(dginnT, dginnbats, by="Gene.name", all.x=T, all.y=T)
dim(tab)
# add column "shared"/"only bats"/"only primates"
tab$status<-""
tab$status[tab$Gene.name %in% common]<-"shared"
tab$status[tab$Gene.name %in% onlyprimates]<-"onlyprimates"
tab$status[tab$Gene.name %in% onlybats]<-"onlybats"
table(tab$status)
write.table(tab, "covid_comp_alldginn.txt", sep="\t")
@
\section{Complete data}
Merge the previous tab with J Young's original table. \textbf{Will replace the 1st part of this script}
<<>>=
young<-read.delim(paste0(workdir,
"data/COVID_PAMLresults_332hits_plusBatScreens_2020_Apr14.csv"),
fill=T, h=T, dec=",")
dim(young)
young$PreyGene<-as.character(young$PreyGene)
young$PreyGene[young$PreyGene=="MTARC1"]<-"MARC1"
@
How many genes in the Young table are not in the DGINN table. And who are they?
<<>>=
table(young$PreyGene %in% tab$Gene.name)
young[(young$PreyGene %in% tab$Gene.name)==FALSE, "PreyGene"]
tab[(tab$Gene.name %in% young$PreyGene)==FALSE, "Gene.name"]
@
Merge them and keep only the krogan genes
<<>>=
tablo<-merge(young, tab, by="Gene.name", all.x=TRUE)
write.table(tablo, "covid_comp_complete.txt", row.names=FALSE, quote=TRUE, sep="\t")
@
\end{document}
data/COVID_PAMLresults_332hits_plusBatScreens_2020_Apr14.csv 0 → 100755
View file @ 2f6095f9
This diff is collapsed.
data/DGINN_202005281339summary_cleaned-LE201108.txt 0 → 100755
View file @ 2f6095f9
This diff is collapsed.
data/DGINN_202005281649summary_cleaned.csv 0 → 100755
View file @ 2f6095f9
This diff is collapsed.