Análisis de cuatro variedades mexicanas de Vitis vinifera y su relación con variedades de vid cultivadas en Europa, Asia y Norte América

Integrantes:

• Jimena Murillo Montero
• Diana Margarita Mojica Muñoz

Descripción y justificación

El cultivo de la vid en América data desde la conquista de la Nueva España. Gracias a la demanda de vino para la evangelización, los viñedos se dispersaron con la formación de misiones en todo el territorio. En Baja California, la región vitivinícola más importante de México, aun prevalecen viñedos históricos, en los que se ha observado buen rendimiento a pesar de las condiciones del cambio climático. Sin embargo, se desconoce la procedencia de las variedades remanentes en estos viñedos. A fin de (aportar en) dilucidar en el origen de vides introducidas por los españoles durante la conquista, se recopilaron los datos de 4 muestras colectadas de viñedos históricos de Baja California, México, así como datos de variedades de vid cultivadas de Europa, Asia y Norte América y se desea comparar sus secuencias de genoma completo anteriormente publicadas, para evaluar su estructura poblacional, y determinar si, las secuencias genómicas de muestras históricas de México, muestran menor variación con respecto a alguno de los grupos antes mencionados. La identificación de las variedades cultivadas en viñedos históricos es importante para el mantenimiento del acervo genético, pues en estas se podría encontrar características importantes para la adaptación.

Antecedentes

La vid (Vitis vinifera L.) es un cultivo de relevancia económica, gracias a su globalización para la producción de vino, y otros usos de la uva. Actualmente existen miles de variedades cultivadas de V. vinifera, sin embargo, el mercado del vino es dominado solo por unos pocos cultivares comerciales (This et. al., 2006).

El origen de su cultivo radica en los países del Caucásico, a partir de su domesticación hace aproximadamente 8 mil años. Desde entonces, la producción de vino guió la distribución del cultivo en países de Europa y Ásia . No obstante, fue hasta la colonización de los países del Nuevo Mundo, que el cultivo de la vid llegó a América (Zhou et. al., 2017). Gracias a los historiadores, se sabe que la evangelización y la formación de misiones en América, permitieron la distribución del cultivo de vid, pues era necesario para la elaboración de vino para consagrar. Para ello, de acuerdo con las ordenanzas de Hernan Cortés, "todo barco que zarpara hacia la Nueva España debía traer consigo vides para plantar". Estas vides debían ser variedades resistentes a las condiciones del viaje. Sin embargo, no se tienen registros confiables sobre el origen ni qué variedades de vid fueron introducidas por los conquistadores. Además, muchos viñedos fueron abandonados con la perdida de misiones (Díaz, 2020). La importancia de poder caracterizar las vides históricas remanentes, radica en el mantenimiento de estos recursos genéticos, pues podrían proporcionar información relevante de resistencia a enfermedades.

Estudios previos, donde analizaron los datos de marcadores moleculares microsatélites, determinaron el posible origen de variedades históricas de Argentina, Chile, Perú y Estados Unidos, encontrando que sus cultivos endémicos corresponden a dos cultivares ancestrales que aun se cultivan en España. Determinaron también que, otros cultivares americanos ancestrales, pueden ser resultado de una mezcla entre variedades cuando se encuentran en simpatría (Tapia et. al., 2007). Otros estudios de datos de resecuenciación de genoma completo, que fueron utilizados para ver la estructura poblacional de variedades silvestres y domesticadas, mediante la realización de un análisis de componentes principales, revelaron un agrupamiento de variedades de vino que reflejan su procedencia estimada (Freitas et. al., 2021).

Descripción del problema y objetivos

Actualmente, las variedades de vid cultivadas en México, corresponden a variedades comerciales. Sin embargo, algunos cultivos de vid históricos introducidos en América durante la epoca de la conquista, aun permanecen remanentes en misiones históricas. No obstante, se desconoce si su origen es de países Europeos, Asiáticos, o si pueden ser resultado de alguna hibridación entre vides cultivadas en América.

Objetivo principal: Dilucidar el posible origen de los cultivos históricos de vid de Baja California, México.

Objetivos específicos:

• Procesar un set de datos, que incluye información sobre el muestreo de variedades cultivadas de vid
• Realizar un análisis de componentes principales (PCA) para ver la estructura y agrupación de las variedades históricas mexicanas, con un set de datos que incluyen variedades de Europa, Asia y Norte América

Descripción del conjunto de datos

El panel de datos a utilizar se recuperó de los datos publicos de secuenciación de genoma completo de Liang et. al. (2019), así como la información correspondiente de cada cultivo. Los datos de secuenciación fueron previamente procesados, haciendo una limpieza de lecturas crudas, alineamiento a genoma de referencia, llamado de variantes, y obtención de eigenvalues para realizar un análisis por componentes principales.

Incluye datos de muestreo de 224 variedades domesticadas de vid. Los cultivares se clasifican dependiendo de su producto final: uva de vino, uva de mesa y uva pasas.

La información de cada cultivo está cubierta en 29 columnas: Type,AvgSpotLen,Bases,BioProject,BioSample,BioSampleModel,Bytes,Center Name,Consent,Cultivar,DATASTORE filetype,DATASTORE provider,DATASTORE region,dev_stage,Experiment,geo_loc_name_country,geo_loc_name_country_continent,geo_loc_name,Instrument,Library Name,LibraryLayout,LibrarySelection,LibrarySource,Number,Organism,Platform,ReleaseDate,create_date,version,Sample Name,SRA Study,TISSUE.

Para el Análisis de Componentes Principales, los datos de secuenciación fueron previamente procesados con el software PLINK 1.9. PLINK convierte los datos genéticos, representados por las letras A, T, C, G, a un formato numérico. Comúnmente, los genotipos son codificados como 0, 1, o 2, donde:

• 0 representa homocigoto para el alelo mayoritario.
• 1 representa heterocigoto.
• 2 representa homocigoto para el alelo minoritario.

A partir de esto, el programa estandariza los datos y construye la matriz de covarianza, que utiliza para descomponer en vectores propios. Estos últimos valores son los datos contenidos en el dataset para construir el gráfico PCA.

Procesamiento y visualización de los datos

Instalación de librerías

!pip install numpy
!pip install pandas
!pip install seaborn
!pip install matplotlib

Requirement already satisfied: numpy in /usr/local/lib/python3.10/dist-packages (1.25.2)
Requirement already satisfied: pandas in /usr/local/lib/python3.10/dist-packages (2.0.3)
Requirement already satisfied: python-dateutil>=2.8.2 in /usr/local/lib/python3.10/dist-packages (from pandas) (2.8.2)
Requirement already satisfied: pytz>=2020.1 in /usr/local/lib/python3.10/dist-packages (from pandas) (2023.4)
Requirement already satisfied: tzdata>=2022.1 in /usr/local/lib/python3.10/dist-packages (from pandas) (2024.1)
Requirement already satisfied: numpy>=1.21.0 in /usr/local/lib/python3.10/dist-packages (from pandas) (1.25.2)
Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.10/dist-packages (from python-dateutil>=2.8.2->pandas) (1.16.0)
Requirement already satisfied: seaborn in /usr/local/lib/python3.10/dist-packages (0.13.1)
Requirement already satisfied: numpy!=1.24.0,>=1.20 in /usr/local/lib/python3.10/dist-packages (from seaborn) (1.25.2)
Requirement already satisfied: pandas>=1.2 in /usr/local/lib/python3.10/dist-packages (from seaborn) (2.0.3)
Requirement already satisfied: matplotlib!=3.6.1,>=3.4 in /usr/local/lib/python3.10/dist-packages (from seaborn) (3.7.1)
Requirement already satisfied: contourpy>=1.0.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib!=3.6.1,>=3.4->seaborn) (1.2.1)
Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.10/dist-packages (from matplotlib!=3.6.1,>=3.4->seaborn) (0.12.1)
Requirement already satisfied: fonttools>=4.22.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib!=3.6.1,>=3.4->seaborn) (4.53.0)
Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib!=3.6.1,>=3.4->seaborn) (1.4.5)
Requirement already satisfied: packaging>=20.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib!=3.6.1,>=3.4->seaborn) (24.0)
Requirement already satisfied: pillow>=6.2.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib!=3.6.1,>=3.4->seaborn) (9.4.0)
Requirement already satisfied: pyparsing>=2.3.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib!=3.6.1,>=3.4->seaborn) (3.1.2)
Requirement already satisfied: python-dateutil>=2.7 in /usr/local/lib/python3.10/dist-packages (from matplotlib!=3.6.1,>=3.4->seaborn) (2.8.2)
Requirement already satisfied: pytz>=2020.1 in /usr/local/lib/python3.10/dist-packages (from pandas>=1.2->seaborn) (2023.4)
Requirement already satisfied: tzdata>=2022.1 in /usr/local/lib/python3.10/dist-packages (from pandas>=1.2->seaborn) (2024.1)
Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.10/dist-packages (from python-dateutil>=2.7->matplotlib!=3.6.1,>=3.4->seaborn) (1.16.0)
Requirement already satisfied: matplotlib in /usr/local/lib/python3.10/dist-packages (3.7.1)
Requirement already satisfied: contourpy>=1.0.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib) (1.2.1)
Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.10/dist-packages (from matplotlib) (0.12.1)
Requirement already satisfied: fonttools>=4.22.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib) (4.53.0)
Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib) (1.4.5)
Requirement already satisfied: numpy>=1.20 in /usr/local/lib/python3.10/dist-packages (from matplotlib) (1.25.2)
Requirement already satisfied: packaging>=20.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib) (24.0)
Requirement already satisfied: pillow>=6.2.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib) (9.4.0)
Requirement already satisfied: pyparsing>=2.3.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib) (3.1.2)
Requirement already satisfied: python-dateutil>=2.7 in /usr/local/lib/python3.10/dist-packages (from matplotlib) (2.8.2)
Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.10/dist-packages (from python-dateutil>=2.7->matplotlib) (1.16.0)

import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
%matplotlib inline

Importar los datos

## Montar Drive
from google.colab import drive
drive.mount('/content/drive')

Mounted at /content/drive

### Carga de datos ###

## Set de datos
dataset = "/content/drive/MyDrive/Proyecto curso python/SraRunTable_CEU.csv"

df = pd.read_csv(dataset)
df

## eigenvectors para PCA
eigenvec = "/content/drive/MyDrive/Proyecto curso python/CEU-2m-pruned_topca.eigenvec"
pca_data = pd.read_table(eigenvec,sep=" ",header=None)

## Set de datos dos (incluye 'uso')
dataset2 = "/content/drive/MyDrive/Proyecto curso python/41467_2019_9135_MOESM4_ESM.xlsx"
df2 = pd.read_excel(dataset2)

## Ver los dataframes
print(df.head())
print(pca_data.head())
print(df2.head())

          Run Assay Type  AvgSpotLen          Bases   BioProject  \
0          S1        WGS         300  not collected          NaN   
1          S2        WGS         300  not collected          NaN   
2          S3        WGS         300  not collected          NaN   
3          S4        WGS         300  not collected          NaN   
4  SRR5891597        WGS         300     7284698100  PRJNA393611   

      BioSample BioSampleModel         Bytes                     Center Name  \
0           NaN          Plant           NaN                       CINVESTAV   
1           NaN          Plant           NaN                       CINVESTAV   
2           NaN          Plant           NaN                       CINVESTAV   
3           NaN          Plant           NaN                       CINVESTAV   
4  SAMN07341438          Plant  2.808682e+09  YUNNAN AGRICULTURAL UNIVERSITY   

  Consent  ... LibrarySource Number        Organism  Platform  \
0  public  ...       GENOMIC  10231  Vitis vinifera  ILLUMINA   
1  public  ...       GENOMIC  10232  Vitis vinifera  ILLUMINA   
2  public  ...       GENOMIC  10233  Vitis vinifera  ILLUMINA   
3  public  ...       GENOMIC  10234  Vitis vinifera  ILLUMINA   
4  public  ...       GENOMIC  10235  Vitis vinifera  ILLUMINA   

            ReleaseDate Sample Name  SRA Study Tissue           create_date  \
0  2020-03-12T00:00:00Z         E12        NaN   leaf  2023-03-05T09:00:00Z   
1  2020-03-12T00:00:00Z        E11A        NaN   leaf  2023-03-05T09:00:00Z   
2  2020-03-12T00:00:00Z          E1        NaN   leaf  2023-03-05T09:00:00Z   
3  2020-03-12T00:00:00Z         E15        NaN   leaf  2023-03-05T09:00:00Z   
4  2019-03-22T00:00:00Z      TA-298  SRP114740   leaf  2017-08-03T05:56:00Z   

  version  
0       1  
1       1  
2       1  
3       1  
4       1  

[5 rows x 33 columns]
      0      1         2         3         4         5         6         7   \
0    E12    E12 -0.014686 -0.016190  0.053925 -0.074761 -0.014066 -0.070570   
1   E11A   E11A -0.014524 -0.015733  0.054707 -0.076237 -0.014130 -0.071402   
2     E1     E1 -0.011885 -0.027222  0.073832 -0.020151  0.062249 -0.057693   
3    E15    E15 -0.014885 -0.023335  0.060087 -0.004371  0.019015 -0.029426   
4  TA-29  TA-29  0.044947 -0.065460  0.070841 -0.055787 -0.104528  0.265947   

         8         9   ...        12        13        14        15        16  \
0 -0.001969 -0.043521  ...  0.003277 -0.092974 -0.006528 -0.002712 -0.068238   
1 -0.002164 -0.043265  ...  0.003843 -0.094429 -0.009198 -0.004783 -0.068545   
2 -0.002218 -0.077083  ...  0.059356 -0.056765 -0.021118  0.022149 -0.038081   
3  0.001844 -0.074204  ...  0.013128 -0.005396 -0.008955  0.007440 -0.020285   
4  0.002019 -0.030907  ... -0.141452  0.167367 -0.139908  0.324391  0.077352   

         17        18        19        20        21  
0  0.074481  0.059498 -0.001976 -0.061126 -0.166349  
1  0.074812  0.060030 -0.001830 -0.060826 -0.167435  
2  0.049547  0.063412 -0.034344 -0.080176 -0.112781  
3  0.060897  0.003371 -0.046752 -0.030797 -0.011192  
4  0.100560  0.076298 -0.233575  0.158972  0.071994  

[5 rows x 22 columns]
  Database ID\n(U: US-National Plant Germplasm System\nV: VitisInternational Variety Catalogue)  \
0                                                NaN                                              
1                                                NaN                                              
2                                                NaN                                              
3                                                NaN                                              
4                                            V. 2624                                              

       Latin Name Cultivar Name\n(Common Name) Sample ID  \
0  Vitis vinifera                Listan Prieto        E1   
1  Vitis vinifera                Rosa del Peru      E11A   
2  Vitis vinifera                Rosa del Peru       E12   
3  Vitis vinifera                     Palomino       E15   
4  Vitis vinifera               Christmas Rose     TA-10   

  Taxonomic Comment \nAccording to \nPhylogenetic Tree  Country of Origin  \
0                                                NaN               Mexico   
1                                                NaN               Mexico   
2                                                NaN               Mexico   
3                                                NaN               Mexico   
4                                                NaN        United States   

  Cultivation\n(W: Wild;  C: Cultivated)  Usage Genetic Classification  \
0                                      C   Wine                    NaN   
1                                      C   Wine                    NaN   
2                                      C   Wine                    NaN   
3                                      C   Wine                    NaN   
4                                      C  Table   Intraspecific Hybrid   

                      Parent 1                          Parent 2 Ploidy\n(N)  \
0                          NaN                               NaN           2   
1                          NaN                               NaN           2   
2                          NaN                               NaN           2   
3                          NaN                               NaN           2   
4  (Hunisa X Emperor) X Nocera  (Hunisa X Emperor) X Pirovano 75           2   

  Fruit Skin Color  Sex of Flower  Berry Weight\n(g)  Sample Source Note  
0            Rouge  Hermaphrodite                NaN         Mexico  NaN  
1            Rouge  Hermaphrodite                NaN         Mexico  NaN  
2            Rouge  Hermaphrodite                NaN         Mexico  NaN  
3            Blanc  Hermaphrodite                NaN         Mexico  NaN  
4            Rouge  Hermaphrodite              8.648  United States  NaN  

Procesamiento set de datos

Seleccionar solo las columnas SampleID, BioSample, Center Name, Cultivar, geo_loc_name_country, geo_loc_name_country_continent, Number, Organism, Tissue, ya que estas contienen la información más relevante del dataset.

## Columnas de interés
columnas = ['Sample Name', 'BioSample', 'Center Name', 'Cultivar', 'geo_loc_name_country',
            'geo_loc_name_country_continent', 'Number', 'Organism', 'Tissue']
## Columnas del df antes
columnas_antes = df.columns
print(columnas_antes)

## Columnas del df después
df = df[columnas]
print(df.columns)

Index(['Run', 'Assay Type', 'AvgSpotLen', 'Bases', 'BioProject', 'BioSample',
       'BioSampleModel', 'Bytes', 'Center Name', 'Consent', 'Cultivar',
       'DATASTORE filetype', 'DATASTORE provider', 'DATASTORE region',
       'dev_stage', 'Experiment', 'geo_loc_name_country',
       'geo_loc_name_country_continent', 'geo_loc_name', 'Instrument',
       'Library Name', 'LibraryLayout', 'LibrarySelection', 'LibrarySource',
       'Number', 'Organism', 'Platform', 'ReleaseDate', 'Sample Name',
       'SRA Study', 'Tissue', 'create_date', 'version'],
      dtype='object')
Index(['Sample Name', 'BioSample', 'Center Name', 'Cultivar',
       'geo_loc_name_country', 'geo_loc_name_country_continent', 'Number',
       'Organism', 'Tissue'],
      dtype='object')

Añadir columna "Usage"

Dependiendo del valor en la columna "Sample Name" en el df, agregar una nueva columna, que indique si tiene uso de uva de mesa (Table), uva de vino (Wine), utilizando los datos del df2

for i in range(len(df)):
    condition = df2["Sample ID"] == df.loc[i, "Sample Name"]
    try:
        value = df2.loc[condition, "Usage"].astype(str).values[0]
        if pd.isna(value) or value == 'nan':  # Check for NaN or 'nan' string
            df.loc[i, "Usage"] = "Unknown"
        else:
            df.loc[i, "Usage"] = value
    except (IndexError, TypeError):
        df.loc[i, "Usage"] = "Unknown"

fields = ["Sample Name", "Usage"]
df.loc[:,fields]

df["Usage"].head(20)

0                 Unknown
1                 Unknown
2                 Unknown
3                 Unknown
4                   Table
5                    Wine
6                   Table
7                    Wine
8                   Table
9                   Table
10                  Table
11            Wine, Table
12            Wine, Table
13                   Wine
14                  Table
15          Table, Raisin
16                  Table
17                  Table
18          Table, Raisin
19    Wine, Table, Raisin
Name: Usage, dtype: object

Eliminar filas donde la columna 'country' tiene valores "uncalculated".

Contar cuántas filas tiene valor "uncalculated" en la columna country, y mostrar el data frame antes y después de eliminar las filas

print(df["geo_loc_name_country"].value_counts)
condicion = df["geo_loc_name_country"] == "uncalculated"


df = df[-condicion]

print(df["geo_loc_name_country"].value_counts)

<bound method IndexOpsMixin.value_counts of 0            Mexico
1            Mexico
2            Mexico
3            Mexico
4             Japan
           ...     
219           Japan
220    uncalculated
221    uncalculated
222    uncalculated
223           Japan
Name: geo_loc_name_country, Length: 224, dtype: object>
<bound method IndexOpsMixin.value_counts of 0       Mexico
1       Mexico
2       Mexico
3       Mexico
4        Japan
        ...   
216    Lebanon
217    Moldova
218    Moldova
219      Japan
223      Japan
Name: geo_loc_name_country, Length: 129, dtype: object>

Procesamiento datos para PCA

Dar nombre a las columnas del dataframe de eigenvectors.

• Primeras dos columnas como Sample Name y Sample ID
• Columnas restantes como PC1, PC2, ... , PCn

# Esto se hizo para generar el nombre de las columnas PC
PC_columns = []
for i in range(1,len(pca_data.columns)-1): # aunque son dos columnas las que no son PC, se quita solo 1, por el valor que no considera el len (que iria de 0 al final)
  caracter = str(i)
  name = 'PC'+ caracter
  PC_columns.append(name) # Esto va añadiendo a la lista cada salida del ciclo
PC_columns

['PC1',
 'PC2',
 'PC3',
 'PC4',
 'PC5',
 'PC6',
 'PC7',
 'PC8',
 'PC9',
 'PC10',
 'PC11',
 'PC12',
 'PC13',
 'PC14',
 'PC15',
 'PC16',
 'PC17',
 'PC18',
 'PC19',
 'PC20']

pca_data.columns = ['Sample Name','Sample ID'] + PC_columns
pca_data.columns
pca_data

	Sample Name	Sample ID	PC1	PC2	PC3	PC4	PC5	PC6	PC7	PC8	...	PC11	PC12	PC13	PC14	PC15	PC16	PC17	PC18	PC19	PC20
0	E12	E12	-0.014686	-0.016190	0.053925	-0.074761	-0.014066	-0.070570	-0.001969	-0.043521	...	0.003277	-0.092974	-0.006528	-0.002712	-0.068238	0.074481	0.059498	-0.001976	-0.061126	-0.166349
1	E11A	E11A	-0.014524	-0.015733	0.054707	-0.076237	-0.014130	-0.071402	-0.002164	-0.043265	...	0.003843	-0.094429	-0.009198	-0.004783	-0.068545	0.074812	0.060030	-0.001830	-0.060826	-0.167435
2	E1	E1	-0.011885	-0.027222	0.073832	-0.020151	0.062249	-0.057693	-0.002218	-0.077083	...	0.059356	-0.056765	-0.021118	0.022149	-0.038081	0.049547	0.063412	-0.034344	-0.080176	-0.112781
3	E15	E15	-0.014885	-0.023335	0.060087	-0.004371	0.019015	-0.029426	0.001844	-0.074204	...	0.013128	-0.005396	-0.008955	0.007440	-0.020285	0.060897	0.003371	-0.046752	-0.030797	-0.011192
4	TA-29	TA-29	0.044947	-0.065460	0.070841	-0.055787	-0.104528	0.265947	0.002019	-0.030907	...	-0.141452	0.167367	-0.139908	0.324391	0.077352	0.100560	0.076298	-0.233575	0.158972	0.071994
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
219	TA-283	TA-283	-0.017452	-0.016879	0.036911	-0.088469	-0.053740	0.047101	-0.000239	-0.009757	...	-0.017475	-0.022053	0.013440	0.009233	-0.030485	0.052002	-0.022209	0.049497	0.012616	-0.075352
220	TA-209	TA-209	-0.025009	0.012676	-0.037408	-0.002282	-0.065865	-0.003752	-0.014067	0.075457	...	-0.011797	-0.045304	-0.083485	0.083388	0.011899	-0.047824	-0.032725	0.019612	0.007137	-0.003353
221	TA-204	TA-204	-0.020803	0.007552	-0.035712	0.015591	-0.085044	-0.007878	-0.012688	0.084355	...	-0.025399	-0.033390	-0.060646	0.097038	0.014038	-0.081245	-0.044973	0.037149	-0.005514	0.022293
222	TA-205	TA-205	-0.019320	0.004735	-0.035951	0.006743	-0.088148	0.001385	-0.007459	0.077632	...	-0.018559	-0.034754	-0.052400	0.100654	0.013964	-0.063865	-0.048938	0.044306	-0.003964	0.005114
223	TA-88	TA-88	-0.015556	-0.012067	0.029640	-0.100086	-0.015574	0.063558	0.002726	0.001665	...	-0.064565	-0.032204	0.044884	-0.019625	-0.070008	-0.024771	-0.028793	0.079466	0.078044	-0.031704

224 rows × 22 columns

## Se pegó la columna de país del df, utilizando la columna en común Sample Name
pca_data = pd.merge(pca_data,df[['Sample Name','geo_loc_name_country',]],on='Sample Name')
pca_data

## Se pegó la columna de Usage del df2, utilizando la columna en común Sample ID
pca_data = pd.merge(pca_data,df2[['Sample ID','Usage',]],on='Sample ID')
pca_data

	Sample Name	Sample ID	PC1	PC2	PC3	PC4	PC5	PC6	PC7	PC8	...	PC13	PC14	PC15	PC16	PC17	PC18	PC19	PC20	geo_loc_name_country	Usage
0	E12	E12	-0.014686	-0.016190	0.053925	-0.074761	-0.014066	-0.070570	-0.001969	-0.043521	...	-0.006528	-0.002712	-0.068238	0.074481	0.059498	-0.001976	-0.061126	-0.166349	Mexico	Wine
1	E11A	E11A	-0.014524	-0.015733	0.054707	-0.076237	-0.014130	-0.071402	-0.002164	-0.043265	...	-0.009198	-0.004783	-0.068545	0.074812	0.060030	-0.001830	-0.060826	-0.167435	Mexico	Wine
2	E1	E1	-0.011885	-0.027222	0.073832	-0.020151	0.062249	-0.057693	-0.002218	-0.077083	...	-0.021118	0.022149	-0.038081	0.049547	0.063412	-0.034344	-0.080176	-0.112781	Mexico	Wine
3	E15	E15	-0.014885	-0.023335	0.060087	-0.004371	0.019015	-0.029426	0.001844	-0.074204	...	-0.008955	0.007440	-0.020285	0.060897	0.003371	-0.046752	-0.030797	-0.011192	Mexico	Wine
4	TA-29	TA-29	0.044947	-0.065460	0.070841	-0.055787	-0.104528	0.265947	0.002019	-0.030907	...	-0.139908	0.324391	0.077352	0.100560	0.076298	-0.233575	0.158972	0.071994	Japan	Table
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
219	TA-283	TA-283	-0.017452	-0.016879	0.036911	-0.088469	-0.053740	0.047101	-0.000239	-0.009757	...	0.013440	0.009233	-0.030485	0.052002	-0.022209	0.049497	0.012616	-0.075352	Japan	Wine
220	TA-209	TA-209	-0.025009	0.012676	-0.037408	-0.002282	-0.065865	-0.003752	-0.014067	0.075457	...	-0.083485	0.083388	0.011899	-0.047824	-0.032725	0.019612	0.007137	-0.003353	uncalculated	Table
221	TA-204	TA-204	-0.020803	0.007552	-0.035712	0.015591	-0.085044	-0.007878	-0.012688	0.084355	...	-0.060646	0.097038	0.014038	-0.081245	-0.044973	0.037149	-0.005514	0.022293	uncalculated	Table
222	TA-205	TA-205	-0.019320	0.004735	-0.035951	0.006743	-0.088148	0.001385	-0.007459	0.077632	...	-0.052400	0.100654	0.013964	-0.063865	-0.048938	0.044306	-0.003964	0.005114	uncalculated	Table
223	TA-88	TA-88	-0.015556	-0.012067	0.029640	-0.100086	-0.015574	0.063558	0.002726	0.001665	...	0.044884	-0.019625	-0.070008	-0.024771	-0.028793	0.079466	0.078044	-0.031704	Japan	Table

224 rows × 24 columns

Como se observó previamente, el set de datos general incluye el nombre de las muestras y el país al que corresponden. Dependiendo del país de procedencia, crea una nueva columna "region", de manera que agrupe a las muestras según corresponda a el país de esa fila.

Se agrupan de la siguiente forma, según el país de procedencia:

• Norte América: USA
• Sudamérica: Brazil
• Europa occidental: United Kingdom, Spain, France, Italy, Switzerland, Germany.
• Europa central: Albania, Austria, Bulgaria, Georgia, Greece, Hungary, Moldova, Romania, Russia, Turkey.
• Africa: Algeria
• Asia oriental: Japan, Korea, China, Uzbekistan, Azerbaijan, Armenia, Lebanon.
• Not collected

# Agrupaciones
mex = ['Mexico']
nor_am = ['USA']
sou_am = ['Brazil']
eur_occ = ['United Kingdom', 'Spain', 'France', 'Italy', 'Switzerland', 'Germany']
eur_cen = ['Albania', 'Austria', 'Bulgaria', 'Georgia', 'Greece', 'Hungary', 'Moldova', 'Romania', 'Russia', 'Turkey']
africa = ['Algeria']
asia_or = ['Japan', 'Korea', 'China', 'Uzbekistan', 'Azerbaijan', 'Armenia', 'Lebanon']

# Función para categorizar
def asignar_region(fila):
  for i in range(0,len(eur_occ)-1):
    if fila['geo_loc_name_country'] == eur_occ[i] :
      return "eur_occ"
  for i in range(0,len(eur_cen)-1):
    if fila['geo_loc_name_country'] == eur_cen[i]:
      return "eur_cen"
  for i in range(0,len(asia_or)):
    if fila['geo_loc_name_country'] == asia_or[i]:
      return "asia_or"
  if fila['geo_loc_name_country'] == nor_am[0]:
    return "nor_am"
  elif fila['geo_loc_name_country'] == sou_am[0]:
    return "sou_am"
  elif fila['geo_loc_name_country'] == africa[0]:
    return "africa"
  elif fila['geo_loc_name_country'] == mex[0]:
    return "mex"

# Aplicar la función a cada fila y agregar una nueva columna "region"
pca_data["region"] = pca_data.apply(asignar_region, axis=1)

# Resultado
pca_data[['Sample Name','region']]

	Sample Name	region
0	E12	mex
1	E11A	mex
2	E1	mex
3	E15	mex
4	TA-29	asia_or
...	...	...
219	TA-283	asia_or
220	TA-209	None
221	TA-204	None
222	TA-205	None
223	TA-88	asia_or

224 rows × 2 columns

## Para eliminar las filas que no tienen valores en región, ya que no tenían información de país
region_NAindex = pca_data[pd.isna(pca_data['region'])].index
pca_data = pca_data.drop(region_NAindex,axis=0) # axis=1 es para columnas
pca_data[['Sample Name','region']]

	Sample Name	region
0	E12	mex
1	E11A	mex
2	E1	mex
3	E15	mex
4	TA-29	asia_or
...	...	...
216	TA-337	asia_or
217	TA-244	eur_cen
218	TA-248	eur_cen
219	TA-283	asia_or
223	TA-88	asia_or

125 rows × 2 columns

Eliminar una de las columnas repetidas de nombre de muestras

pca_data = pca_data.drop("Sample Name",axis=1) # axis=1 es para columnas
pca_data.columns

Index(['Sample ID', 'PC1', 'PC2', 'PC3', 'PC4', 'PC5', 'PC6', 'PC7', 'PC8',
       'PC9', 'PC10', 'PC11', 'PC12', 'PC13', 'PC14', 'PC15', 'PC16', 'PC17',
       'PC18', 'PC19', 'PC20', 'geo_loc_name_country', 'Usage', 'region'],
      dtype='object')

Procesamiento datos para K-means a partir de variables del PCA

# Clustering
from sklearn.cluster import AffinityPropagation
from sklearn.cluster import KMeans
from sklearn.cluster import AgglomerativeClustering
from sklearn.cluster import DBSCAN

from sklearn import metrics
from sklearn.preprocessing import OneHotEncoder
import pandas as pd
from scipy.cluster.hierarchy import dendrogram, linkage, fcluster
import matplotlib.pyplot as plt

#Se calcula el error cuadrado con k diferentes (distorción)
pca_data.dropna(inplace=True)
data2 = pca_data.drop(columns =["Sample ID","region","geo_loc_name_country","Usage"])
data2 = data2.to_numpy()

distortions = []
for i in range(1, 11):
    km = KMeans(
        n_clusters=i, init='random',
        n_init=10, max_iter=300,
        tol=1e-04, random_state=0)
    km.fit(data2)
    distortions.append(km.inertia_)

# plot
plt.plot(range(1, 11), distortions, marker='o')
plt.xlabel('Número de agrupamientos')
plt.ylabel('Distorción')
plt.show()

Se decide seleccionar 6 clústers distintos, que es la zona donde se ve una disminución en la baja de la distorción.

Resultados

Gráficos descriptivos de los datasets utilizados

Uso de las Uvas analizadas

frequency = df["Usage"].value_counts()

# Create bar plot
plt.figure(figsize=(15, 5))
plt.bar(frequency.index, frequency.values, color='lightblue')
plt.xlabel('Usage')
plt.ylabel('Frecuencia')
plt.title('Uso de la Uva')
plt.show()

Origen de la Uva Analizada

frequency = pca_data["region"].value_counts()

# Create bar plot
plt.figure(figsize=(15, 5))
plt.bar(["Asia Oriental", "Europa occidental", "Europa central", "Norte América", "Sur América", "África", "México"], frequency.values, color='green')
plt.xlabel('Continente')
plt.ylabel('Frecuencia')
plt.title('Origen de la Uva')
plt.show()

La mayoría de las uvas analizadas provienen de Asia Orienta y hay una cantidad reducida que provienen de Sur América y África.

Mapa con la frecuencias de las uvas según el país de origen

con = df["geo_loc_name_country"] == "USA"
df.loc[con,"geo_loc_name_country"] = "United States of America"
frequency = df["geo_loc_name_country"].value_counts()
frequency

geo_loc_name_country
China                       29
Japan                       24
United States of America    15
France                      15
Italy                        5
Romania                      4
Hungary                      3
Greece                       3
Uzbekistan                   3
Germany                      3
Georgia                      3
Lebanon                      2
Moldova                      2
Armenia                      2
Azerbaijan                   2
Spain                        1
Switzerland                  1
Brazil                       1
Russia                       1
United Kingdom               1
Bulgaria                     1
Algeria                      1
Albania                      1
Austria                      1
Turkey                       1
Name: count, dtype: int64

import geopandas as gpd

data = {
    'Country': frequency.index,
    'Frequency': frequency.values
}
data = pd.DataFrame(data)

# Load world shapefile from geopandas datasets
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))

# Merge your data with the world GeoDataFrame
world = world.merge(data, how='left', left_on='name', right_on='Country')

# Aggregate data by continent
continent_freq = world.groupby('continent').agg({'Frequency': 'sum'}).reset_index()

base = world.boundary.plot(edgecolor='gray', figsize=(20, 15))

# Merge the continent frequency data back to the world dataframe for mapping
world = world.merge(continent_freq, on='continent', suffixes=('', '_continent'))

world = world[world['continent'] != "Antarctica"]

# Plot
world.plot(column='Frequency_continent', ax=base, legend=True,
           legend_kwds={'label': "Frequency by Continent"},
           cmap='OrRd')  # Color map can be changed according to preference

plt.title('Continent Frequency Map')
plt.show()

<ipython-input-117-17dd33b5ead8>:10: FutureWarning: The geopandas.dataset module is deprecated and will be removed in GeoPandas 1.0. You can get the original 'naturalearth_lowres' data from https://www.naturalearthdata.com/downloads/110m-cultural-vectors/.
  world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))

Principal Component Analysis

Regiones

El análisis de componentes principales no mostró mucha varibilidad entre las muestras, pues todas son variedades cultivadas sin tanta distancia genética. Sin embargo, es posible distinguir los grupos de regiones agrupados por colores. Las muestras que más se alejan, posiblemente sean variedades silvestes.

colors_PC23 = ["#3366CC","#DC3912","#FF9900","#109618","#990099","#0099C6","#DD4477"]
plt.figure(figsize=(12, 6))
sns.scatterplot(data=pca_data, x="PC1", y="PC2", hue="region",palette=colors_PC23)

<Axes: xlabel='PC1', ylabel='PC2'>

Para el caso de las variedades mexicanas de interés, es posible notar que estas se agrupan cerca de las variedades de Europa Occidental, lo cual cumple con la hipótesis de que las muestras provienen de España, y no de regiones de Asia o Norte América.

colors_PC23 = ["#990099","#DC3912","#DE9ED6","#FF9900","#0099C6","#3366CC","#109618"]
plt.figure(figsize=(12, 6))
pca_plot = sns.scatterplot(data=pca_data, x="PC2", y="PC3", hue="region",palette=colors_PC23)

Uso

El PCA también nos permitió ver las agrupaciones dependiendo del uso que se le da a los cultivos, lo cuál indica que sí hay diferencias genéticas dependiendo si es uva de mesa o uva de vino.

plt.figure(figsize=(12, 6))
sns.scatterplot(data=pca_data, x="PC1", y="PC2", hue="Usage")

<Axes: xlabel='PC1', ylabel='PC2'>

Así mismo, es clara la distinción entre uvas de vino, uvas de mesa, y aquellas que cumplen ambas funciones, quedan justo entre estos dos grupos, exponiendo claramente el principio de un PCA.

plt.figure(figsize=(12, 6))
sns.scatterplot(data=pca_data, x="PC1", y="PC3", hue="Usage")

<Axes: xlabel='PC1', ylabel='PC3'>

Resultados del K-means

# Generación de clusters

Z = linkage(data2, method='complete')
max_clusters = 6
clusters = fcluster(Z, max_clusters, criterion='maxclust')
plt.figure(figsize=(10, 7))
dendrogram(Z)

# Agregar una línea horizontal para mostrar el corte
max_d = Z[-(max_clusters-1), 2]  # max_clusters-1 porque Z tiene una combinación menos que el número de clusters
plt.axhline(y=max_d, c='k')

plt.title("Dendrograma con corte en {} clusters".format(max_clusters))
plt.xlabel("Muestras")
plt.ylabel("Distancia")
plt.show()

# Imprimir las asignaciones de clusters
pca_data["cluster"] = clusters

cross_tab = pd.crosstab(pca_data['cluster'], pca_data['region'])

# Mostrar la tabla cruzada
print(cross_tab)

region   africa  asia_or  eur_cen  eur_occ  nor_am  sou_am
cluster                                                   
1             0        2        0        0       0       0
2             1       58       18       22      15       1
3             0        1        0        0       0       0
4             0        1        0        0       0       0
5             0        0        0        1       0       0
6             0        0        1        0       0       0

Conclusiones

• Gran parte de las muestras del dataset inicial carece de información relevante, como país de procedencia o nombre del cultivo. Es importante no considerar las filas con información faltante, pues puede afectar el análisis.

• El conjunto de datos utilizado proporciona información relevante para inferir el origen de las muestras de vid históricas de México. Como se observa en el análisis de componentes principales, las muestras mexicanas están más cerca de aquellas vides de Europa Occidental, como era de esperarse, sabiendo que fueron los españoles quienes introdujeron la vid en México. Sin embargo, aun hay mucha similitud con las variedades de Asia y Europa Oriental, pues su tiempo de divergencia no es tan grande, haciendo que la variación no sea tan notoria.

• El PCA que agrupa a las muestras por su uso, permite también observar las diferencias entre aquellas que se usan como uva de mesa, y uva para vino, observando que las vides mexicanas históricas se agrupan mayormente con las uvas para vino, pues no existe evidencia de mezcla con cultivos de uva de mesa.

• Tomando en cuenta las diferencias genéticas a través de K-means podemos observar que la mayoría de las uvas analizadas tienen grandes similitudes y sólo unas cuanta suvas de Asia, Europa Central y Europa Occidental se distinguen del resto.

Referencias

• Díaz, M. (2020). The vine and the wine in the world of Hernán Cortés. Fundación Academia Europea e Iberoamericana de Yuste.
• Freitas, S., Gazda, M., Rebelo, M. (2021). Pervasive hybridization with local wild relatives in Western European grapevine varieties. Science Advances, 8584, 7(47).
• Liang, Z., Duan, S., Sheng, J. (2019). Whole-genome resequencing of 472 Vitis accessions for grapevine diversity and demographic history analyses. Nature Communications, 1-12, 10(1).
• Tapia, A., Cabezas, J., Cabello, F. (2007). Determining the Spanish origin of representative ancient American grapevine varieties. American Journal of Enology and Viticulture, 242-251, 58(2)
• This, P., Lacombe, T., Thomas, M. (2006). Historical origins and genetic diversity of wine grapes. Trends in Genetics, 511-519, 22(9).
• Zhou, Y., Massonnet, M., Sanjak, J. (2017). Evolutionary genomics of grape (Vitis vinifera ssp. vinifera) domestication. Proceedings of the National Academy of Sciences of the United States of America, 11715-11720, 114(44).