diff --git a/src/extp2.R b/src/extp2.R deleted file mode 100644 index ed23afb7b725a459958e2815d06f9353bdbbbabe..0000000000000000000000000000000000000000 --- a/src/extp2.R +++ /dev/null @@ -1,28 +0,0 @@ - -wdir ="/amuhome/s23014817/m2bsgreprod/sources" -dir.create(wdir, showWarnings = F, recursive = T) -setwd(wdir) - -library(devtools) - -if (!require("BiocManager", quietly = TRUE)) #boucle qui permet de dire que s'il n'est pas installé il faut l'installer mais s'il est déjà installé pas besoin - install.packages("BiocManager") - -if (!require("snpStats", quietly = TRUE)) - BiocManager::install("snpStats") - -if (!require("SNPRelate", quietly = TRUE)) - BiocManager::install("SNPRelate") - -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, penncath -penncath_bed_path = "/amuhome/s23014817/m2bsgreprod/results/data/penncath.bed" -penncath_bim_path = "/amuhome/s23014817/m2bsgreprod/results/data/penncath.bim" -penncath_fam_path = "/amuhome/s23014817/m2bsgreprod/results/data/penncath.fam" -geno <- read.plink(penncath_bed_path, penncath_bim_path, penncath_fam_path, select.snps=sample(1:861473, 25000, replace = FALSE ), na.strings = ("-9")) - - - - diff --git a/src/tp2rmd.R b/src/tp2rmd.R deleted file mode 100644 index dbf0a3a0b69ebe9522be5cfc0d0d1d7154cfde3f..0000000000000000000000000000000000000000 --- a/src/tp2rmd.R +++ /dev/null @@ -1,123 +0,0 @@ -#on télécharge le fichier de données - -load("results/tp1/TP1_asbvg.RData") - -#on crée un dossier dans lequel se trouveront les fichiers de sortie - -output_dir ="results/tp2" -dir.create(output2_dir, showWarnings = F, recursive = T) - -#on installe les librairies nécessaires - -if (!require("BiocManager", quietly = TRUE)) # - install.packages("BiocManager") - -if (!require("SNPRelate", quietly = TRUE)) - BiocManager::install("SNPRelate") - -#on charge la librairie - -#library(SNPRelate) - -#on ouvre le fichier contenant les données SNPs - -#genofile <- SNPRelate::snpgdsOpen(fichier$gds, readonly = TRUE) - -# De nouveau on effectue une étape d'élagage -# Placer le seuil de LD à 0.2 - -ld.thresh <- 0.2 - -set.seed(1000) -geno.sample.ids <- rownames(genotype) -snpSUB <- SNPRelate::snpgdsLDpruning(genofile, ld.threshold = ld.thresh, - sample.id = geno.sample.ids, # On analyse uniquement les individus filtrés - snp.id = colnames(genotype)) # On analyse uniquement les SNPs filtrés - -snpset.pca <- unlist(snpSUB, use.names=FALSE) -cat(length(snpset.pca),"SNPs will be used in PCA analysis. \n") - -pca <- SNPRelate::snpgdsPCA(genofile, sample.id = geno.sample.ids, snp.id = snpset.pca, num.thread=1) - -# Calcul et sauvegarde des 10 composantes principales en un data frame où chaque colonne est une PC - -pcs <- data.frame(FamID = pca$sample.id, pca$eigenvect[,1 : 10], stringsAsFactors = FALSE) -colnames(pcs)[2:11]<-paste("pc", 1:10, sep = "") - -print(head(pcs)) # Exercice 1 -closefn.gds(genofile) -save(genotype, genoBim, clinical, pcs, file=working.data.fname(5)) - -#Exercice 1 - -# Lecture des données de 1000G pour le chromosome 16 -pedfile.bed <- read.pedfile(plink_base.bed$ped, snps = plink_base.bed$info, which=1) -pedfile.fam <- read.pedfile(plink_base.fam$ped, snps = plink_base.fam$info , which=1) -pedffile.bim <- read.pedfile(plink_base.bim$ped, snps = plink_base.bim$info, which=1) - - -# Obtenir/stocker les données génotypiques pour le chromosome 16 -genoMatrix <- thougeno$genotypes - -# Obtenir/stocker les positions chromosomiques pour chaque SNP -support <- thougeno$map -colnames(support)<-c("SNP", "position", "A1", "A2") -head(support) - -# Imputation des SNPs non génotypés issus des données 1000G -# Générer une "sous"-matrice comportant les SNPs génotypés et filtrés du chromosome 16 -presSnps <- colnames(genotype) -presDatChr <- genoBim[genoBim$SNP %in% presSnps & genoBim$chr==16, ] -targetSnps <- presDatChr$SNP - -# Générer deux "sous"-matrice de SNPs issus de 1000G : "present" (matrice des SNPs 1000G typés dans nos données) et "missing" (matrice des SNPs 1000G non-typés dans nos données). Ces matrices sont nécessaires à l'établissement des règles d'imputation -is.present <- colnames(genoMatrix) %in% targetSnps -missing <- genoMatrix[,!is.present] -present <- genoMatrix[,is.present] -print(missing) -print(present) - -#Exercice 2 - -# Obtenir et stocker les positions des SNPs utilisés dans les règles d'imputation -pos.pres <- support$position[is.present] -pos.miss <- support$position[!is.present] - -# Calcul et sauvegarde des règles d'imputation à l'aide de snp.imputation() -rules <- snp.imputation(present, missing, pos.pres, pos.miss) - -#Exercice 3 - -# Suppression des imputations ayant échouées -rules <- rules[can.impute(rules)] -cat("Imputation rules for", length(rules), "SNPs were estimated\n") - -#Exercice 4 - -# Contrôles-qualité : incertitude et faible MAF estimée -# Fixation des seuils -r2threshold <- 0.7 -minor <- 0.01 - -# Filtre sur l'incertitude et la faible MAF -rules <- rules[imputation.r2(rules) >= r2threshold] -cat(length(rules),"imputation rules remain after imputations with low certainty were removed\n") - -rules <- rules[imputation.maf(rules) >= minor] -cat(length(rules),"imputation rules remain after MAF filtering\n") - -# Obtenir le génotype des SNPs imputés -target <- genotype[,targetSnps] -imputed <- impute.snps(rules, target, as.numeric=FALSE) -print(imputed) - - -#on crée le fichier de sortie - -rdata_path=file.path(output_dir, "TP2_asbvg.RData") -save.image(rdata2_path) - - - - -