L3 GeBAP - TP Analyse de QTL
From silico.biotoul.fr
Line 1: | Line 1: | ||
+ | <span style='color: #990000'> !! Si vous travaillez avec Rstudio Cloud, il suffit de vous connecter avec vos identifiants sur https://rstudio.cloud, et de créer un "New Project" !! </span> | ||
+ | |||
== Lecture et exploration du tableau de données == | == Lecture et exploration du tableau de données == | ||
Revision as of 08:12, 1 February 2021
!! Si vous travaillez avec Rstudio Cloud, il suffit de vous connecter avec vos identifiants sur https://rstudio.cloud, et de créer un "New Project" !!
Contents |
Lecture et exploration du tableau de données
Lecture:
pheno=read.table("http://silico.biotoul.fr/site/images/8/89/L3P_phenotypes.txt", sep="\t", header=TRUE) # nécessite une connexion internet
Accéder directement aux variables simplement en donnant leurs noms:
attach(pheno) names(pheno)
Résumé des données et histogramme
summary(prop_brunissement) hist(prop_brunissement,col="orange",main="distribution de la proportion de brunissement") abline(v=1,lwd=3,lty=2) text(0.95,40,"F83") abline(v=0.5,lwd=3,lty=2) text(0.45,40,"A17")
Analyse QTl avec le package qtl de R
Installer et Charger le package "qtl"
install.packages("qtl") library(qtl)
Lecture d'un fichier de données pour R/qtl
data <- read.cross("csv", file="http://silico.biotoul.fr/site/images/8/84/L3P_data.csv")
Séparer légèrement les marqueurs confondus
data <- jittermap(data); data
Résumé des données et représentation graphique (commandes spécifiques du package "qtl")
plot(data) # info données manquantes, carte génétique et phénotypes markernames(data) # noms des marqueurs nchr(data) # nombre de chromosomes nmar(data) # nombre de marqueurs par chromosomes
Analyse QTL avec méthode ANOVA
out.mr <- scanone(data, method="mr",pheno.col=1) # analyse ANOVA pour chaque marqueur out.mr # résultats interval mapping plot(out.mr, ylab="LOD score") # graphe du LOD score de l'ANOVA pour chaque marqueur plot(out.mr[out.mr$chr==3,], ylab="LOD score") # zoom sur chromosome 3 max(out.mr) # marqueur avec LOD score le plus fort summary(out.mr,threshold =3) # LOD score le plus important pour chaque chromosome, avec un seuil à 3 plotPXG(data, "mtic742",pheno.col=1) # boxplots phénotypiques pour les génotypes au marqueur le plus significatif avec l'ANOVA
Analyse QTL avec méthode Interval Mapping (IM)
data <- calc.genoprob(data, step=1, error.prob=0.001) # calculer la probabilité de chaque génotype possible à une position donnée, en fonction des marqueurs bordants out.em <- scanone(data, method="em",pheno.col=1) # analyse interval mapping pour chaque marqueur vrai (plus marqueurs entre 2 marqueurs bordants) out.em # résultats interval mapping plot(out.em,ylab="LOD score",col="blue",ylim=c(0,16)) # graphe du LOD de l'interval mapping pour chaque marqueur plot(out.em,out.mr, col=c(4,1), lty=1:2, ylab="LOD score",ylim=c(0,16)) # comparaison résultats interval mapping et ANOVA plot(out.em,out.mr, chr=3,col=c(4,1), lty=1:2, ylab="LOD score",ylim=c(0,16)) # zoom sur chromosome 3 max(out.em) # marqueur avec LOD score le plus fort summary(out.em,threshold=3) # LOD score le plus important pour chaque chromosome, avec un seuil à 3
Estimation de l'intervalle de confiance du QTL sur le chromosome 3 (on peut aussi le faire pour le QTL du chromosome 4)
out.em[out.em$chr==3,] lodint(out.em, 3, 1.5) # 1.5 LOD support (chute de 1.5 du LOD de part et d'autre du QTL) lodint(out.em, 3, 1.5, expandtomarkers=T) # étend l'intervalle aux marqueurs les plus proches # à faire pour le qtl du chromosome 4
Boxplots phénotypiques pour les génotypes aux marqueurs aux bornes de l'intervalle de confiance [mtic879,mtic1049] => englobe mtic742
par(mfrow=c(1,3)) plotPXG(data, "mtic879",pheno.col=1) plotPXG(data, "mtic742",pheno.col=1) plotPXG(data, "mtic1049",pheno.col=1)
Héritabilité due à chaque QTL (% de variabilité phénotypique expliquée par le marqueur)
# heritabilité du premier QTL (chromosome 3) 1-10^(-2*summary(out.em,threshold=3)[1,]$lod/nind(data)) # heritabilité du deuxième QTL (chromosome 4) 1-10^(-2*summary(out.em,threshold=3)[2,]$lod/nind(data))