Esempi di statistica descrittiva e inferenziale/Analisi del proprio genoma
Installazione librerie
modificainstall.packages("BiocManager")
BiocManager::install("gwascat")
Caricamento librerie
modifica library(dplyr)
library(ggplot2)
library(stringr)
library(gwascat)
Parte 1: Dati
modificaLa società californiana 23andMe fornisce a pagamento in tutto il mondo un dataset contenente i dati relativi al proprio genotipo, ottenuti analizzando un campione di saliva inviato dal cliente per posta. Questi dati possono essere utilizzati per scopi di ricerca, educativi ed informativi ma non per uso medico. Ciascuna riga del dataset corrisponde ad un singolo SNP contenente 4 colonne :
- l'identificatore univoco rsid dell'SNP
- il cromosoma contenente l'SNP
- la posizione numerica all'interno del DNA nel cromosoma
- il genotipo formato da 2 lettere tra le seguenti: adenina (A), citosina (C), guanina (G), e timina (T).
Un signore ha reso pubblico il proprio dataset di SNP fornito da 23andME ed è possibile scaricarlo da qui : https://www.kaggle.com/zusmani/mygenome
Caricamento del dataset:
genome_zeeshan_usmani <- read.csv("genome_zeeshan_usmani.csv")
I primi 6 SNP contenuti nel dataset:
head(genome_zeeshan_usmani)
rsid chromosome position genotype 1 rs12564807 1 734462 AA 2 rs3131972 1 752721 AG 3 rs148828841 1 760998 AC 4 rs12124819 1 776546 AA 5 rs115093905 1 787173 GG 6 rs11240777 1 798959 GG
E' possibile scaricare un dataset di links a pubblicazioni scientifiche relative ad SNP contenenti l' allele di rischio collegato ad una certa malattia o predisposizione. Ad esempio nell'SNP rs533123 l'allele di rischio potrebbe essere A per cui bisogna indagare se il proprio genotipo relativo a quell'SNP lo contiene. Quindi unendo i 2 datasets si potrà sapere se il proprio genotipo contiene quel particolare allele di rischio relativo ad uno specifico tratto o malattia.
Caricamento dataset:
options(timeout = 300)
updated_gwas_data <- as.data.frame(makeCurrentGwascat())
Le variabili contenute nel dataset di pubblicazioni scientifiche sono le seguenti:
names(updated_gwas_data)
[1] "seqnames" "start" [3] "end" "width" [5] "strand" "DATE.ADDED.TO.CATALOG" [7] "PUBMEDID" "FIRST.AUTHOR" [9] "DATE" "JOURNAL" [11] "LINK" "STUDY" [13] "DISEASE.TRAIT" "INITIAL.SAMPLE.SIZE" [15] "REPLICATION.SAMPLE.SIZE" "REGION" [17] "CHR_ID" "CHR_POS" [19] "REPORTED.GENE.S." "MAPPED_GENE" [21] "UPSTREAM_GENE_ID" "DOWNSTREAM_GENE_ID" [23] "SNP_GENE_IDS" "UPSTREAM_GENE_DISTANCE" [25] "DOWNSTREAM_GENE_DISTANCE" "STRONGEST.SNP.RISK.ALLELE" [27] "SNPS" "MERGED" [29] "SNP_ID_CURRENT" "CONTEXT" [31] "INTERGENIC" "RISK.ALLELE.FREQUENCY" [33] "P.VALUE" "PVALUE_MLOG" [35] "P.VALUE..TEXT." "OR.or.BETA" [37] "X95..CI..TEXT." "PLATFORM..SNPS.PASSING.QC." [39] "CNV" "MAPPED_TRAIT" [41] "MAPPED_TRAIT_URI" "STUDY.ACCESSION" [43] "GENOTYPING.TECHNOLOGY"
Parte 2 : Esplorazione dati
modificagenome_zeeshan_usmani %>%
group_by(chromosome) %>%
summarise(total=n()) %>%
mutate(chromosome=reorder(chromosome,total)) %>%
ggplot(aes(chromosome,total, fill=chromosome))+
geom_bar(stat="identity")+
coord_flip()+
geom_text(aes(label=total), hjust=0, size=2)+
guides(fill=FALSE)+
ylab("Totale SNP") +
xlab("Cromosoma") +
ggtitle("Numero totale di SNP per cromosoma ",subtitle = "nel genoma umano fornito")
E' possibile cercare nel sito https://www.snpedia.com/ gli SNP relativi ad una certa malattia o tratto e cercarli nel genoma del signore . Ad esempio relativamente al rischio alcolico questo è l'SNP del signore:
genome_zeeshan_usmani %>%
filter(rsid=='rs1799971')
rsid chromosome position genotype 1 rs1799971 6 154360797 AG
Si uniscono i 2 datasets: quello relativo agli SNP del signore e quello relativo alle pubblicazioni scientifiche . In tal modo si avranno a disposizione gli alleli di rischio e le pubblicazioni scientifiche relativi agli SNP del signore.
output_data <- inner_join(genome_zeeshan_usmani, updated_gwas_data, by = c("rsid" = "SNPS"))
Si visualizzano le variabili di interesse nelle prime 6 righe:
output_data %>%
select(DISEASE.TRAIT, STRONGEST.SNP.RISK.ALLELE, genotype,STUDY)%>%
head()
DISEASE.TRAIT STRONGEST.SNP.RISK.ALLELE genotype 1 Skin reflectance (Melanin index) rs2272756-A AG 2 Pancreatic cancer rs13303010-G AA 3 Body mass index rs3934834-G CT 4 Urate levels rs9442380-T CC 5 Epithelial ovarian cancer rs9442387-? CT 6 Blood protein levels rs3766186-C AC
STUDY 1 A genome-wide association study of skin and iris pigmentation among individuals of South Asian ancestry. 2 Genome-wide meta-analysis identifies five new susceptibility loci for pancreatic cancer. 3 Linkage and genome-wide association analysis of obesity-related phenotypes: association of weight with the MGAT1 gene. 4 Urate, Blood Pressure, and Cardiovascular Disease: Evidence From Mendelian Randomization and Meta-Analysis of Clinical Trials. 5 Genome-wide association studies identify susceptibility loci for epithelial ovarian cancer in east Asian women. 6 Co-regulatory networks of human serum proteins link genetics to disease.
Si estraggono dal dataset le righe relative agli alleli di rischio contenuti nel genotipo del signore.Ad esempio se l'allele è T e nel genotipo si trova AT.
allele <- str_sub(output_data$STRONGEST.SNP.RISK.ALLELE, -1)
df <- data.frame()
for (r in 1:nrow(output_data) ) {
if (allele[r]=="?" | allele[r]=="*") next
if (str_detect(output_data[r,4], allele[r]))
df <- rbind(df,output_data[r,])
if(r%%1000==0) print(r)
}
df %>%
select(DISEASE.TRAIT, STRONGEST.SNP.RISK.ALLELE, genotype,STUDY) %>%
head()
DISEASE.TRAIT STRONGEST.SNP.RISK.ALLELE genotype 1 Skin reflectance (Melanin index) rs2272756-A AG 6 Blood protein levels rs3766186-C AC 11 Serum alkaline phosphatase levels rs1123571-A AG 16 Body mass index rs7535528-G GG 18 Autoimmune thyroid disease rs2234167-A AG 19 Medication use (thyroid preparations) rs2234167-A AG
Si estraggono dal dataset le righe relative agli alleli di rischio contenuti nel genotipo del signore 2 volte. Ad esempio se l'allele è A e nel genotipo si trova AA.
df1 <- data.frame()
allele <- str_sub(df$STRONGEST.SNP.RISK.ALLELE, -1)
for (r in 1:nrow(df) ) {
if (str_count(df[r,4], allele[r])==2)
df1 <- rbind(df1,df[r,])
if(r%%1000==0) print(r)
}
df1 %>%
select(DISEASE.TRAIT, STRONGEST.SNP.RISK.ALLELE, genotype,STUDY)
DISEASE.TRAIT STRONGEST.SNP.RISK.ALLELE genotype 16 Body mass index rs7535528-G GG 23 Body mass index rs6667605-T TT 35 Serum total cholesterol levels rs1456465-G GG 36 Mean corpuscular volume rs1569419-C CC 37 Plateletcrit rs1569419-C CC 38 Platelet distribution width rs1569419-C CC
df %>%
group_by(DISEASE.TRAIT) %>%
summarise(total =n())%>%
mutate(DISEASE.TRAIT=reorder(DISEASE.TRAIT,total)) %>%
filter(total>100) %>%
ggplot(aes(DISEASE.TRAIT,total, fill=DISEASE.TRAIT))+
geom_bar(stat="identity")+
coord_flip()+
geom_text(aes(label=total), hjust=0, size=2)+
guides(fill=FALSE)+
ylab("Numero di pubblicazioni") +
xlab("Tratto di malattia") +
ggtitle("Numero di pubblicazioni con singolo allele di rischio",subtitle = "per tratto di malattia maggiore di 100 nel genoma umano fornito")
df1 %>%
group_by(DISEASE.TRAIT) %>%
summarise(total =n())%>%
mutate(DISEASE.TRAIT=reorder(DISEASE.TRAIT,total)) %>%
filter(total>50) %>%
ggplot(aes(DISEASE.TRAIT,total, fill=DISEASE.TRAIT))+
geom_bar(stat="identity")+
coord_flip()+
geom_text(aes(label=total), hjust=0, size=2)+
guides(fill=FALSE)+
ylab("Numero di pubblicazioni") +
xlab("Tratto di malattia") +
ggtitle("Numero di pubblicazioni con DOPPIO allele di rischio",subtitle = "per tratto di malattia maggiore di 100 nel genoma umano fornito")
Per ogni SNP si può leggere lo studio fornito tramite il suo link o cercarlo su SNPedia ...
Ad esempio per la schizofrenia :
df1 %>%
filter(DISEASE.TRAIT=="Schizophrenia") %>%
select(DISEASE.TRAIT, STRONGEST.SNP.RISK.ALLELE, genotype,STUDY) %>%
head()
DISEASE.TRAIT STRONGEST.SNP.RISK.ALLELE genotype
1 Schizophrenia rs533123-A AA
2 Schizophrenia rs3001723-A AA
3 Schizophrenia rs1198588-T TT
4 Schizophrenia rs1389724-T TT
5 Schizophrenia rs7521837-T TT
6 Schizophrenia rs7527939-C CC
STUDY
1 Genome-wide association study of schizophrenia in Ashkenazi Jews.
2 Genome-wide association study of schizophrenia in Ashkenazi Jews.
3 Genome-wide association analysis identifies 13 new risk loci for schizophrenia.
4 Genome-wide association study of schizophrenia in Ashkenazi Jews.
5 Genome-wide association study of schizophrenia in Ashkenazi Jews.
6 Whole-genome-wide association study in the Bulgarian population reveals HHAT as schizophrenia susceptibility gene.