Esempi di statistica descrittiva e inferenziale/Analisi del proprio genoma

Indice del libro

Installazione librerieModifica

install.packages("BiocManager")
BiocManager::install("gwascat")

Caricamento librerieModifica

 library(dplyr)
 library(ggplot2)
 library(stringr)
 library(gwascat)


Parte 1: DatiModifica

La 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 datiModifica

genome_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.