script

src/TP1.v2.R

+options(repos = c(CRAN = "https://cloud.r-project.org"))
+
+# Définir le répertoire de travail de manière dynamique
+m2bsgreprod <- normalizePath(".")
+data.dir <- file.path(m2bsgreprod, "data")
+output_dir <- file.path(m2bsgreprod, "results", "TP1")
+
+# Créer les répertoires si nécessaire
+dir.create(output_dir, showWarnings = FALSE, recursive = TRUE)
+
+# Fichiers d'entrées
+gwas.fn<-lapply(c(bed='bed', fam='fam', bim ='bim', gds='gds'), function(n) sprintf("%s/GWAStutorial.%s", data.dir, n))
+gwas.fn.2<-lapply(c(bed='bed', fam='fam', bim ='bim', gds='gds'), function(n) sprintf("%s/data.%s", data.dir, n))
+clinical.fn<-sprintf("%s/GWAStutorial_clinical.csv", data.dir)
+onethou.fn<-lapply(c(info='info', ped='ped'), function(n) sprintf("%s/chr16_1000g_CEU.%s", data.dir, n))
+protein.coding.coords.fname<-sprintf("%s/ProCodgene_coords.csv", data.dir)
+
+# Sauvegarde des "modules"
+working.data.fname <- function(num) { sprintf("%s/working_%s.Rdata", output_dir, num) }
+
+# Les données analysées nécessitant beaucoup de RAM, nous allons sélectionner aléatoirement 250000 SNPs et réecrire des fichiers bed, bim, fam
+penncath_bed_path = file.path(data.dir, "penncath.bed")
+penncath_bim_path = file.path(data.dir, "penncath.bim")
+penncath_fam_path = file.path(data.dir, "penncath.fam")
+
+geno <- snpStats::read.plink(penncath_bed_path, penncath_bim_path, penncath_fam_path, select.snps=sample(1:861473, 25000, replace = FALSE ), na.strings = ("-9"))
+
+plink_base=file.path(output_dir, "plink_base")
+snpStats::write.plink(plink_base, snps=geno$genotypes, pedigree=geno$fam[,1], id=geno$fam[,1], mother=geno$fam[,4], sex=geno$fam[,5], phenotype=geno$fam[,6], chromosome = geno$map[,1], genetic.distance = geno$map[,3], position = geno$map[,4], allele.1 = geno$map[,5], allele.2 = geno$map[,6], na.code = ("-9"))
+
+genoBim<-geno$map
+colnames(genoBim)<-c("chr", "SNP", "gen.dist", "position", "A1", "A2")
+
+genotype<-geno$genotype
+
+genoFam<-geno$fam
+
+# On commence par libérer de l'espace
+rm(geno)
+
+# On charge le fichier clinique
+clinical<- read.csv(clinical.fn, colClasses = c("character", "factor", "factor", rep("numeric", 4)))
+rownames(clinical)<-clinical$FamID
+print(head(clinical))
+
+# On conserve uniquement les individus présents dans les 2 fichiers : génotype et clinique
+genotype <- genotype[clinical$FamID, ]
+print(genotype)
+
+save(genotype, genoBim, clinical, file = working.data.fname(1))
+
+# Creation d'un résumé statistique lié aux SNPs
+snpsum.col<-col.summary(genotype)
+print(head(snpsum.col))
+
+# Création des seuils
+call<-0.95
+minor<-0.01
+
+# Application des filtres sur la MAF et le call rate
+use<- with(snpsum.col, (!is.na(MAF) & MAF > minor) & Call.rate >= call)
+use[is.na(use)]<-FALSE
+cat(ncol(genotype)-sum(use), "SNPs will be removed due to low MAF or call rate. \n")
+
+# Creation d'un jeu de données et de résumé stat excluant les SNPs n'ayant pas passé les seuils
+genotype<- genotype[,use]
+snpsum.col<- snpsum.col[use,]
+print(genotype)
+
+save(genotype, snpsum.col, genoBim, clinical, file=working.data.fname(2))
+
+# Création des statistiques liés aux échantillons (Call rate, hétérozygotie)
+snpsum.row<-row.summary(genotype)
+
+# Ajout du coefficient de consanguinité|F| à snpsum.row
+MAF<-snpsum.col$MAF
+callmatrix<-!is.na(genotype) # matrice de valeurs logiques les génotypes non NA prennent la valeur logique TRUE 
+hetExp<-callmatrix %*% (2*MAF*(1-MAF)) # produit matriciel matrice logique par le vecteur numérique de la fréquence de génotype hétérozygote => en gros on additionne les 2pq des SNPs génotypés (non NA, valeur logique TRUE) pour chaque individu (donc la proportion de SNPs hétérozygotes attendus) 
+hetObs<-with(snpsum.row, Heterozygosity*(ncol(genotype))*Call.rate) # calcule la proportion de SNP hétérozygotes observés
+snpsum.row$hetF<-1-(hetObs/hetExp)
+
+# Création des seuils
+sampcall<-0.95
+hetcutoff<-0.1
+
+sampleuse<-with(snpsum.row, !is.na(Call.rate) & Call.rate >sampcall & abs(hetF)<=hetcutoff)
+sampleuse[is.na(sampleuse)]<-FALSE # élimine les données manquantes
+cat(nrow(genotype)-sum(sampleuse), "Subjects will beremoved due to low sample call rate or inbreeding coefficient")
+
+# Creation d'un jeu de données (incluant données cliniques et génotypes) excluant les échantillons (individus) n'ayant pas passé les seuils
+genotype<-genotype[sampleuse,]
+clinical<-clinical[rownames(genotype),]
+
+# Vérification des liens de parenté
+ld.thresh<-0.2 # seuil de déséquilibre de liaison
+kin.thresh<-0.1 # seuil de parenté
+
+# Création d'un fichier gds nécessaire pour les fonctions de la library SNPRelate
+snpgdsBED2GDS(gwas.fn.2$bed, gwas.fn.2$fam, gwas.fn.2$bim, gwas.fn.2$gds)
+genofile<-openfn.gds(gwas.fn.2$gds, readonly=FALSE)
+
+# Suppression des suffixes -1 automatiquement ajoutés aux échantillons
+gds.ids<-read.gdsn(index.gdsn(genofile, "sample.id"))
+gds.ids<-sub("-1", "", gds.ids)
+add.gdsn(genofile, "sample.id", gds.ids, replace=TRUE)
+
+# Elagage des SNPs en DL pour l'analyse IBD
+set.seed(1000) # la position de début de l'analyse est aléatoire. set seed permet de fixer cette position aléatoire pour chaque itération
+geno.sample.ids<- rownames(genotype)
+snpSUB<-snpgdsLDpruning(genofile, ld.threshold = ld.thresh, 
+                        sample.id = geno.sample.ids, # Analyse seulement les ech. filtr?s
+                        snp.id = colnames(genotype)) # Analyse seulement les SNPs filtr?s
+snpset.ibd<-unlist(snpSUB, use.names=FALSE)
+cat(length(snpset.ibd), "SNPs will be used in IBD analysis\n")
+
+# Recherche des coefficients IBD par paire à l'aide de la méthode des moments.  Le résultat est une table indiquant le lien de parenté entre chaque paire d'échantillon
+ibd<-snpgdsIBDMoM(genofile, kinship = TRUE, sample.id=geno.sample.ids, snp.id = snpset.ibd, num.thread = 1)
+ibdcoeff<-snpgdsIBDSelection(ibd) # comparaison par paire d'individus
+head(ibdcoeff)
+#k0: probabilité de partager ZERO alleles, k1   probabilité de partager 1 allele, kinship: coefficient de parenté
+
+# Recherche des individus présentant un lien de parenté potentiel
+ibdcoeff <- ibdcoeff[ ibdcoeff$kinship >= kin.thresh, ]
+
+# Suppression itérative des échantillons/individus avec un lien de parenté élevé en débutant par l'échantillon ayant le plus d'appariement
+related.samples <- NULL
+while ( nrow(ibdcoeff) > 0 ) {
+    sample.counts <- arrange(count(c(ibdcoeff$ID1, ibdcoeff$ID2)), -freq)
+    rm.sample <- sample.counts[1, 'x']
+    cat("Removing sample", as.character(rm.sample), 'too closely related to', 
+        sample.counts[1, 'freq'],'other samples.\n')
+
+    ibdcoeff <- ibdcoeff[ibdcoeff$ID1 != rm.sample & ibdcoeff$ID2 != rm.sample,]
+    related.samples <- c(as.character(rm.sample), related.samples)
+}
+
+# Filtre de la matrice genotype et du data frame clinical afin de ne conserver que les individus sans lien de parenté
+genotype <- genotype[ !(rownames(genotype) %in% related.samples), ]
+clinical <- clinical[ !(clinical$FamID %in% related.samples), ]
+
+geno.sample.ids <- rownames(genotype)
+cat(length(related.samples), "similar samples removed due to correlation coefficient >=", kin.thresh,"\n") 
+print(genotype) 
+
+pca <- snpgdsPCA(genofile, sample.id = geno.sample.ids,  snp.id = snpset.ibd, num.thread=1)
+
+# Création d'un data frame des 2 composantes principales (PC)
+pop_pheno <- read.gdsn(index.gdsn(genofile, "sample.annot/phenotype"))
+samp.id <- read.gdsn(index.gdsn(genofile, "sample.id"))
+pctab <- data.frame(sample.id = pca$sample.id,
+                  CAD = factor(pop_pheno)[match(pca$sample.id, samp.id)],
+                  PC1 = pca$eigenvect[,1],    # 1ère composante principale
+                  PC2 = pca$eigenvect[,2],    # 2nde composante principale
+                  stringsAsFactors = FALSE)
+head(pctab)
+
+# Graphe des 2 PCs
+plot(pctab$PC2, pctab$PC1, col=as.integer(pctab$CAD), xlab="Principal Component 2", ylab="Principal Component 1", main = "Ancestry Plot")
+legend("bottomright", legend=levels(pctab$CAD), pch="o", col=1:4)
+
+
+closefn.gds(genofile)
+save(genotype, genoBim, clinical, file=working.data.fname(3))
+
+# Exclusions des SNPs déviant de l'HWE (uniquement pour les témoins)
+hardy <- 10^-6 # Seuil d'HWE
+
+CADcontrols <- as.character(clinical[ clinical$CAD==0, 'FamID' ])
+snpsum.colCont <- col.summary(genotype[CADcontrols,])
+HWEuse <- with(snpsum.colCont, !is.na(z.HWE) & ( abs(z.HWE) < abs( qnorm(hardy/2) ) ) )
+rm(snpsum.colCont)
+
+HWEuse[is.na(HWEuse)] <- FALSE          # Suppression des valeurs manquantes
+cat(ncol(genotype)-sum(HWEuse),"SNPs will be removed due to high HWE.\n") 
+
+# Suppression des SNPs déviant de l'HWE de la matrice genotype  
+genotype <- genotype[,HWEuse]
+print(genotype)           
+
+save(genotype, genoBim, clinical, file=working.data.fname(4))
+
+rdata_path = file.path(output_dir, "TP1_asbvg.RData")
+save.image(rdata_path)
+

src/TP1.v3.R

+options(repos = c(CRAN = "https://cloud.r-project.org"))
+output_dir = "results/tp1"
+dir.create(output_dir, showWarnings = F, recursive = T)
+
+# Les données analysées nécessitant beaucoup de RAM, nous allons sélectionner aléatoirement 250000 SNPs et réecrire des fichiers bed, bim, fam
+penncath_bed_path = "results/data/penncath.bed"
+penncath_bim_path = "results/data/penncath.bim"
+penncath_fam_path = "results/data/penncath.fam"
+
+geno <- snpStats::read.plink(penncath_bed_path, penncath_bim_path, penncath_fam_path, select.snps=sample(1:861473, 25000, replace = FALSE ), na.strings = ("-9"))
+
+
+plink_base=file.path(output_dir, "plink_base")
+snpStats::write.plink(plink_base, snps=geno$genotypes, pedigree=geno$fam[,1], id=geno$fam[,1], mother=geno$fam[,4], sex=geno$fam[,5], phenotype=geno$fam[,6], chromosome = geno$map[,1], genetic.distance = geno$map[,3], position = geno$map[,4], allele.1 = geno$map[,5], allele.2 = geno$map[,6], na.code = ("-9"))
+
+genoBim<-geno$map
+colnames(genoBim)<-c("chr", "SNP", "gen.dist", "position", "A1", "A2")
+#head(genoBim)
+
+genotype<-geno$genotype
+#print(genotype)
+
+genoFam<-geno$fam
+#head(genoFam)
+
+# On commence par libérer de l'espace
+rm(geno)
+
+rdata_path = file.path(output_dir, "TP1_asbvg.RData")
+save.image(rdata_path)
+

src/tp1.R

@@ -3,6 +3,7 @@ dir.create(wdir, showWarnings = F, recursive = T)
 setwd(wdir)
 
 #library(devtools)
+#library(plyr)
 
 if (!require("BiocManager", quietly = TRUE))
   install.packages("BiocManager")
@@ -13,10 +14,13 @@ if (!require("snpStats", quietly = TRUE))
 if (!require("SNPRelate", quietly = TRUE))
 BiocManager::install("SNPRelate")
 
+# Charger les bibliothèques
+library(snpStats)
+library(SNPRelate)
+
 # Les données analysées nécessitant beaucoup de RAM, nous allons sélectionner aléatoirement 250000 SNPs et réecrire des fichiers bed, bim, fam
 penncath_bed_path = "results/data/penncath.bed"
 penncath_bim_path = "results/data/penncath.bim"
 penncath_fam_path = "results/data/penncath.fam"
 geno <- snpStats::read.plink(penncath_bed_path, penncath_bim_path, penncath_fam_path, select.snps=sample(1:861473, 25000, replace = FALSE ), na.strings = ("-9"))
 
-