diff --git a/src/TP1.v2.R b/src/TP1.v2.R new file mode 100644 index 0000000000000000000000000000000000000000..5bb53cd9367251bb9fae37a49d20b538dc442fa2 --- /dev/null +++ b/src/TP1.v2.R @@ -0,0 +1,184 @@ +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) + diff --git a/src/TP1.v3.R b/src/TP1.v3.R new file mode 100644 index 0000000000000000000000000000000000000000..15b1f1a9b55d579c72a69ff487853acd0a9612a9 --- /dev/null +++ b/src/TP1.v3.R @@ -0,0 +1,31 @@ +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) + diff --git a/src/tp1.R b/src/tp1.R index 72bdb7e2866c2d310e8d5c4303c813cf7a63f777..6d6f7d55f9d82fcb7cc6a66254124b34b9f6189e 100644 --- a/src/tp1.R +++ b/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")) -