From 08c57764807609aac73cbd21a09be5952c095cd1 Mon Sep 17 00:00:00 2001
From: OTT Oceane <o22025448@V-PP-47-L-084.salsa.univ-amu.fr>
Date: Wed, 23 Oct 2024 17:00:54 +0200
Subject: [PATCH] k
---
src/TP1.v2.R | 190 ---------------------------------------------------
src/TP1.v3.R | 39 -----------
2 files changed, 229 deletions(-)
delete mode 100644 src/TP1.v2.R
delete mode 100644 src/TP1.v3.R
diff --git a/src/TP1.v2.R b/src/TP1.v2.R
deleted file mode 100644
index 0db5778..0000000
--- a/src/TP1.v2.R
+++ /dev/null
@@ -1,190 +0,0 @@
-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)
-
-# 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")
-
-clinical_csv_path = "data/GWAStutorial_clinical.csv"
-
-# 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) }
-
-
-
-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_csv_path, 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
deleted file mode 100644
index abefe18..0000000
--- a/src/TP1.v3.R
+++ /dev/null
@@ -1,39 +0,0 @@
-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)
-
-# On charge le fichier clinique
-clinical_csv_path = "data/GWAStutorial_clinical.csv"
-clinical<- read.csv(clinical_csv_path, colClasses = c("character", "factor", "factor", rep("numeric", 4)))
-rownames(clinical)<-clinical$FamID
-#print(head(clinical))
-
-protein.coding.coords.fname<-"data/ProCodgene_coords.csv"
-
-
-rdata_path = file.path(output_dir, "TP1_asbvg.RData")
-save.image(rdata_path)
-
--
GitLab