# ANALYSE RRBS # Donnees : population EPP, individus de terrain (F0) et individus du labo (F5) # Traiment en amont : trimming, QC, bismark et samtools # Fichiers de depart : sortie samtools en sorted.bam # ech 1 à 24 : EPP F5 (23 ech) # ech 69 à 110 : EPP FO (40 ech) - ces fichiers sont déjà en .txt # Ech supprimés : 9 et 21 ; 70, 72 #------------------------------ 1. Environnement de travail -------------------------------------------------- getwd() ls() #voir tous les objets R qui y sont crees rm(list =ls()) #effacer les objets R creer pendant la session R #------------------------------ 2. Charger les packages necessaires ------------------------------------------ BiocManager::install(c("Rqc","methylKit","ShortRead","QuasR")) library(methylKit) library(Rqc) library(ShortRead) library(QuasR) packageVersion("methylKit") #------------------------------ 3. Methylation calling ------------------------------------------------------- #The methylation call files are basically text files that contain percent methylation score per base. # STEP 1. From sorted Bismark alignment files (.sorted.bam) to methylRawList objects (.txt) using processBismarkAln function ?processBismarkAln getwd() analysisFiles <- list("E:/rrbs/Bam_sorted/RRBS/1.sorted.bam", "E:/rrbs/Bam_sorted/RRBS/2.sorted.bam", "E:/rrbs/Bam_sorted/RRBS/3.sorted.bam", "E:/rrbs/Bam_sorted/RRBS/4.sorted.bam", "E:/rrbs/Bam_sorted/RRBS/5.sorted.bam", "E:/rrbs/Bam_sorted/RRBS/6.sorted.bam", "E:/rrbs/Bam_sorted/RRBS/7.sorted.bam", "E:/rrbs/Bam_sorted/RRBS/8.sorted.bam", "E:/rrbs/Bam_sorted/RRBS/9.sorted.bam", "E:/rrbs/Bam_sorted/RRBS/11.sorted.bam", "E:/rrbs/Bam_sorted/RRBS/12.sorted.bam", "E:/rrbs/Bam_sorted/RRBS/13.sorted.bam", "E:/rrbs/Bam_sorted/RRBS/14.sorted.bam", "E:/rrbs/Bam_sorted/RRBS/15.sorted.bam", "E:/rrbs/Bam_sorted/RRBS/16.sorted.bam", "E:/rrbs/Bam_sorted/RRBS/17.sorted.bam", "E:/rrbs/Bam_sorted/RRBS/18.sorted.bam", "E:/rrbs/Bam_sorted/RRBS/19.sorted.bam", "E:/rrbs/Bam_sorted/RRBS/20.sorted.bam", "E:/rrbs/Bam_sorted/RRBS/21.sorted.bam", "E:/rrbs/Bam_sorted/RRBS/22.sorted.bam", "E:/rrbs/Bam_sorted/RRBS/23.sorted.bam", "E:/rrbs/Bam_sorted/RRBS/24.sorted.bam") sampleIDs <- list("1","2","3","4","5","6","7","8","9","11","12","13","14","15","16","17","18","19","20","21","22","23","24") treatmentSpecification <- c(rep(1, times = 23)) my.methRaw <- processBismarkAln(location = analysisFiles, sample.id = sampleIDs, assembly="ASM164957v2", read.context="CpG", treatment = treatmentSpecification, save.folder=getwd()) my.methRaw #------------------------------ 1. Selection des fichiers .txt ------------------------- sampleIDs <- list("1","2","3","4","5","6","7","8","11","12","13","14","15","16","17","18","19","20","22","23","24", "69","71","73","74","75","76","77","78","79","80","81","82","83","84","85","86","87","88","89","90", "91","92","93","94","95","96","97","98","99","100","101","102","103","104","105","106","107","108","109","110") treatmentSpecification <- c(rep(1, times = 21),rep(2, times = 40)) #F5 = Groupe 1 ; F0 = Groupe 2 #Covariate = length, pour representer l'age covar=data.frame(length=c(2.858, 2.732, 2.837, 2.956, 2.752, 2.84, 2.77, 2.87, 2.672, 2.772, 2.795, 2.687, 2.674, 2.507, 2.756, 2.77, 2.744, 2.706, 2.799, 2.599, 2.833, 3.322, 3.552, 2.841, 3.903, 3.184, 2.885, 3.367, 2.603, 2.882, 3.154, 2.763, 3.12, 3.033, 3.358, 2.967, 2.792, 2.829, 3.227, 3.228, 3.005, 2.764, 2.766, 3.033, 2.917, 2.712, 2.616, 2.733, 2.845, 2.599, 3.032, 2.725, 2.798, 2.934, 2.834, 2.607, 2.894, 2.597, 3.024, 2.637, 3.298)) View(covar) txt.files <-list("E:/rrbs/Bam_sorted/RRBS/1_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/2_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/3_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/4_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/5_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/6_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/7_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/8_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/11_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/12_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/13_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/14_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/15_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/16_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/17_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/18_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/19_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/20_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/22_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/23_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/24_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/69_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/71_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/73_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/74_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/75_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/76_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/77_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/78_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/79_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/80_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/81_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/82_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/83_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/84_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/85_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/86_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/87_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/88_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/89_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/90_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/91_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/92_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/93_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/94_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/95_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/96_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/97_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/98_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/99_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/100_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/101_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/102_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/103_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/104_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/105_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/106_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/107_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/108_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/109_CpG.txt", "E:/rrbs/Bam_sorted/RRBS/110_CpG.txt") txt.files length(txt.files) # STEP 2. Reading the methylation call files with methRead function -> create a methylation object from the coverage files, and include sample ID and treatment information. ?methRead myobj= methRead(txt.files, sample.id = sampleIDs, assembly = "ASM164957v2", treatment = treatmentSpecification, context = "CpG", dbdir = getwd(), mincov=10) #nbr de fois que ca a ete seq par un meme endroit # --> infos sur la localisation genomique des C, infos sur les C methyles/non methyles length(myobj) class(myobj) head(myobj) myobj #Sauvergader l'objet cree pour ne pas repasser par le long process de creation de myobj saveRDS(myobj, file = "myobj.RDS") rm(list = ls()) getwd() myobj=readRDS("myobj.RDS") myobj length(myobj) class(myobj) head(myobj) #------------------------------ 2. Descriptive statistics on samples ------------------------------------------------------- #Controle de qualite supplementaire #les valeurs de methylation doivent avoir une distribution bimodale #il faut que la plupart des CpG ont soit une methylation elevee, soit une methylation faible #deviation par rapport a la bimodalite peuvent etre due a une mauvaise qualite experimentales (prob traitement bisulfite) #1) Distribution of %methylation : #---------------------------------- getMethylationStats(myobj[[1]],plot=F,both.strands=FALSE) getMethylationStats(myobj[[1]],plot=TRUE,both.strands=FALSE) getMethylationStats(myobj[[5]],plot=TRUE,both.strands=FALSE) getMethylationStats(myobj[[10]],plot=TRUE,both.strands=FALSE) getMethylationStats(myobj[[15]],plot=TRUE,both.strands=FALSE) getMethylationStats(myobj[[20]],plot=TRUE,both.strands=FALSE) getMethylationStats(myobj[[25]],plot=TRUE,both.strands=FALSE) getMethylationStats(myobj[[30]],plot=TRUE,both.strands=FALSE) getMethylationStats(myobj[[35]],plot=TRUE,both.strands=FALSE) getMethylationStats(myobj[[40]],plot=TRUE,both.strands=FALSE) getMethylationStats(myobj[[45]],plot=TRUE,both.strands=FALSE) getMethylationStats(myobj[[50]],plot=TRUE,both.strands=FALSE) getMethylationStats(myobj[[55]],plot=TRUE,both.strands=FALSE) getMethylationStats(myobj[[60]],plot=TRUE,both.strands=FALSE) getMethylationStats(myobj[[61]],plot=TRUE,both.strands=FALSE) #2) Analyse des coverages et filtre de myobj : #--------------------------------------------- #verifier les valeurs de couverture #les bases avec une couverture tres elevee sont alarmantes --> prob de biais PCR dans la couverture experimentale getCoverageStats(myobj[[1]], plot = T, both.strands = F) getCoverageStats(myobj[[1]], plot = F, both.strands = F) # Donne les valeurs. Valeur max = 8164 getCoverageStats(myobj[[5]], plot = T, both.strands = F) getCoverageStats(myobj[[5]], plot = F, both.strands = F) # Donne les valeurs. Valeur max = 45529 getCoverageStats(myobj[[10]], plot = T, both.strands = F) getCoverageStats(myobj[[10]], plot = F, both.strands = F) # Donne les valeurs. Valeur max = 18279 getCoverageStats(myobj[[15]], plot = T, both.strands = F) getCoverageStats(myobj[[15]], plot = F, both.strands = F) # Donne les valeurs. Valeur max = 18221 getCoverageStats(myobj[[20]], plot = T, both.strands = F) getCoverageStats(myobj[[20]], plot = F, both.strands = F) # Donne les valeurs. Valeur max = 17109 getCoverageStats(myobj[[25]], plot = T, both.strands = F) getCoverageStats(myobj[[25]], plot = F, both.strands = F) # Donne les valeurs. Valeur max = 217081 getCoverageStats(myobj[[30]], plot = T, both.strands = F) getCoverageStats(myobj[[30]], plot = F, both.strands = F) # Donne les valeurs. Valeur max = 355442 getCoverageStats(myobj[[40]], plot = T, both.strands = F) getCoverageStats(myobj[[40]], plot = F, both.strands = F) # Donne les valeurs. Valeur max = 355237 getCoverageStats(myobj[[50]], plot = T, both.strands = F) getCoverageStats(myobj[[50]], plot = F, both.strands = F) # Donne les valeurs. Valeur max = 226422 getCoverageStats(myobj[[55]], plot = T, both.strands = F) getCoverageStats(myobj[[55]], plot = F, both.strands = F) # Donne les valeurs. Valeur max = 388719 getCoverageStats(myobj[[60]], plot = T, both.strands = F) getCoverageStats(myobj[[60]], plot = F, both.strands = F) # Donne les valeurs. Valeur max = 42978 filtered.myobj=filterByCoverage(myobj, lo.count =10, lo.perc = NULL, hi.count = NULL, hi.perc = 99.9 ) filtered.myobj getCoverageStats(filtered.myobj[[1]], plot = T, both.strands = F) getCoverageStats(filtered.myobj[[1]], plot = F, both.strands = F) # Donne les valeurs. Valeur max = 814 getCoverageStats(filtered.myobj[[5]], plot = T, both.strands = F) getCoverageStats(filtered.myobj[[5]], plot = F, both.strands = F) # Donne les valeurs. Valeur max = 773 getCoverageStats(filtered.myobj[[10]], plot = T, both.strands = F) getCoverageStats(filtered.myobj[[10]], plot = F, both.strands = F) # Donne les valeurs. Valeur max = 681 getCoverageStats(filtered.myobj[[20]], plot = T, both.strands = F) getCoverageStats(filtered.myobj[[20]], plot = F, both.strands = F) # Donne les valeurs. Valeur max = 739 getCoverageStats(filtered.myobj[[25]], plot = T, both.strands = F) getCoverageStats(filtered.myobj[[25]], plot = F, both.strands = F) # Donne les valeurs. Valeur max = 652 getCoverageStats(filtered.myobj[[30]], plot = T, both.strands = F) getCoverageStats(filtered.myobj[[30]], plot = F, both.strands = F) # Donne les valeurs. Valeur max = 667 getCoverageStats(filtered.myobj[[40]], plot = T, both.strands = F) getCoverageStats(filtered.myobj[[40]], plot = F, both.strands = F) # Donne les valeurs. Valeur max = 622 getCoverageStats(filtered.myobj[[50]], plot = T, both.strands = F) getCoverageStats(filtered.myobj[[50]], plot = F, both.strands = F) # Donne les valeurs. Valeur max = 642 getCoverageStats(filtered.myobj[[60]], plot = T, both.strands = F) getCoverageStats(filtered.myobj[[60]], plot = F, both.strands = F) # Donne les valeurs. Valeur max = 594 saveRDS(filtered.myobj, file = "my_filtered_obj.RDS") rm(list = ls()) filtered.myobj=readRDS("my_filtered_obj.RDS") filtered.myobj length(filtered.myobj) class(filtered.myobj) head(filtered.myobj) #CpG #Fusion d'echantillons dans un seul tableau #3) Rassembler en un objet les 4 sous-objets: Objet Meth #---------------------------------------------------------- ?unite methall=unite(filtered.myobj, destrand = T,min.per.group = NULL) methall # 25717 rows = nbr CpGs etudies qu'on retrouve chez tous head(methall) class(methall) head(methall) #filter les CpG en fonction de la variation #supprimer les CpG qui ne sont pas variables #supprimer les C qui sont potentiellement des mutations C-->T #4) Coefficient of variablity of object Meth: #-------------------------------------------- pm=percMethylation(methall) #Creer la matrice de pourcentage de methylation par CpG head(pm, 10) #Donne la matrice de pourcentage de methylation par CpG summary(pm) View(pm) write.csv2(x=pm, file = "pm.csv") rm(list =ls()) #effacer les objets R crees pendant la session R pm <- read.csv2("pm.csv") #!!!! Ajouter le titre "Position" a la 1ère colonne (80803 CpGs) View(pm) CV=function(x){sd(x)/mean(x)*100} #On cree une formule pour calculer les coefficients de correlation F5=apply(pm[,1:21], MARGIN = 1, FUN = CV) #MARGIN : 1 pour les lignes, 2 pour les colonnes. Ici, on veut comparer les lignes. summary(F5) #65.515 hist(F5) F0=apply(pm[,22:61], MARGIN = 1, FUN = CV) #MARGIN : 1 pour les lignes, 2 pour les colonnes. Ici, on veut comparer les lignes. summary(F0) #47.8599 hist(F0) #------------------------------ 5. Analyse multivariee de meth------------------------------------------------------- #1) PCA and clustering : #------------------------ #Coefficient de correlation getCorrelation(methall, plot=F) #Affichage des coef de correlation entre les conditions getCorrelation(methall, plot=T) #Clustering #deux clusters sont fusionnés si la variance est minimisée par rapport aux autres opérations de fusion possibles #construire un dendogramme qui montre la relation entre les clusters, utilise la distance de corrélation et la méthode de Ward pour le regroupement hiérarchique ?clusterSamples clusterSamples(methall, dist = "correlation", method = "ward.D", plot = T) #Regroupement en clusters #PCA ?PCASamples #screeplot pour la variance expliquee des composantes principales PCASamples(methall, screeplot=TRUE) #nuage de points de nos echantillons sur les axes PC1 et PC2 #voir comment ils se regroupent dans ces dimensions #on s'attend a voir des echantillons similaires proches les uns des autres pc=PCASamples(methall,obj.return = T,comp = c(1,2)) prcomp(methall) summary(pc) pc$x write.csv2(pc$x, file="pca_epigenet") DMC=calculateDiffMeth(methall, overdispersion = "MN", test="Chisq", covariates = covar, adjust = "SLIM") DMC #25717 all.DMC=getMethylDiff(DMC, difference=10, qvalue=0.01, type="all") #4232 rows all.DMC=getMethylDiff(DMC, difference=30, qvalue=0.01, type="all") #2350 rows (2172 et 178) all.DMC=getMethylDiff(DMC, difference=50, qvalue=0.01, type="all") #131 all.DMC=getMethylDiff(DMC, difference=70, qvalue=0.01, type="all") #13 (10 et 3) all.DMC hyper.DMC=getMethylDiff(DMC, difference=10, qvalue=0.01, type="hyper") hyper.DMC hypo.DMC=getMethylDiff(DMC, difference=10, qvalue=0.01, type="hypo") hypo.DMC #REGIONS #------------------------------ 3. TILES creation and descriptive analysis ------------------------- #1) Création des regions (ou tiles): #--------------------------------- tiles300 <- tileMethylCounts(filtered.myobj, win.size = 300, step.size = 300, cov.bases = 10) #region de 300pb avec un saut de 300 #toujours le meme coverage tiles300[[1]] # On a bien des tiles de 300 pb head(tiles300[[1]]) length(tiles300) class(tiles300) #2) Rassembler en un objet les 4 sous-objets: function "unite" #------------------------------------------------------------- tiles300_30=methylKit::unite(tiles300, destrand = T, min.per.group = 21L) #30 = min de poisson de population tiles300_30 #5799 rows = nombre de regions de 300 bp etudiees #3) Coefficient of variablity of object tiles: #---------------------------------------------- CV=function(x){sd(x)/mean(x)*100} #On cree une formule pour calculer les coefficients de correlation pm300=percMethylation(tiles300_30) head(pm300, 10) summary(pm300) F5_300=apply(pm300[,1:21], MARGIN = 1, FUN = CV) summary(F5_300) #mean=43.6540 hist(F5_300) F0_300=apply(pm300[,22:61], MARGIN = 1, FUN = CV) summary(F0_300) #mean=17.8 hist(F0_300) #4) PCA and clustering : #------------------------ #Coefficient de correlation getCorrelation(tiles300_30, plot=F) #Affichage des coef de correlation entre les conditions getCorrelation(tiles300_30, plot=T) #Clustering clusterSamples(tiles300_30, dist = "correlation", method = "ward", plot = T) #Regroupement en clusters #PCA pc=PCASamples(tiles300_30, obj.return = T, adj.lim = c(1,1), comp = c(1,2)) pc=PCASamples(tiles300_30) #------------------------------ 4. Finding differentially methylated regions between populations ------------ #1) Regression logisitique #------------------------- #Fonction calculateDiffMeth : calcule les differences de methylation #Fonction get.methylDiff -> Pour extraire ce qui est significatif ?calculateDiffMeth dm.lr.DMR=calculateDiffMeth(tiles300_30, overdispersion = "MN", test="Chisq", covariates = covar, adjust = "SLIM") dm.lr.DMR all.diff.DMR=getMethylDiff(dm.lr.DMR, difference=10, qvalue=0.01, type="all") #232 rows (215 et 17) all.diff.DMR=getMethylDiff(dm.lr.DMR, difference=32, qvalue=0.01, type="all") #2 rows all.diff.DMR=getMethylDiff(dm.lr.DMR, difference=25, qvalue=0.01, type="all") #6 rows all.diff.DMR=getMethylDiff(dm.lr.DMR, difference=20, qvalue=0.01, type="all") #13 rows (12 et 1) all.diff.DMR hyper.DMR=getMethylDiff(dm.lr.DMR, difference=25, qvalue=0.01, type="hyper") hyper.DMR hypo.DMR=getMethylDiff(dm.lr.DMR, difference=25, qvalue=0.01, type="hypo") hypo.DMR #We can also visualize the distribution of hypo/hyper-methylated bases/regions per chromosome using the following function. #In this case, the example set includes only one chromosome diffMethPerChr(dm.lr.DMR,plot=TRUE,qvalue.cutoff=0.01, meth.cutoff=25) getContext(dm.lr.DMR) #------------------------------ 5. Annotating differentially methylated bases or regions ----------------------------------- #Exporter en format csv les objets "all.diff.DMR", "hyper.DMR" et "hypo.DMR" pour les importer dans SeqMonk write.csv(all.diff.DMR.bis, file = "DMR_all") write.csv(hyper.DMR, file = "DMR_hyper") write.csv(hypo.DMR, file = "DMR_hypo") write.csv2(dm.lr.DMR.bis, file = "analyzed_regions")