Skip to content
Snippets Groups Projects
Commit 08c57764 authored by OTT Oceane's avatar OTT Oceane
Browse files

k

parent 30fa1a81
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/TP1.v2.R deleted 100644 → 0
+ 0
− 190
View file @ 30fa1a81
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)
src/TP1.v3.R deleted 100644 → 0
+ 0
− 39
View file @ 30fa1a81
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)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment