scripts R

src/extp2.R

-
-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"))
-
-
-
-

src/tp2rmd.R

-#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) 
-
-
-
-
-