**Curso Análisis de datos de secuenciación masiva en genómica de poblaciones** Dra. Paulina Mejía Ruíz #### Usaremos vcftools para filtrar los datos por individuos y loci ##### https://vcftools.sourceforge.net/man_latest.html ## Missing data por individuos --missing-indv genera un archivo que reporta la falta de datos por individuo con terminación .imiss ``` vcftools --vcf populations.snps.vcf --missing-indv cat out.imiss #Para ver el archivo ``` Para ver el histograma de frecuencias: ``` mawk '!/IN/' out.imiss | cut -f5 > totalmissing gnuplot << \EOF set terminal dumb size 120, 30 set autoscale unset label set title "Histogram of % missing data per individual" set ylabel "Number of Occurrences" set xlabel "% of missing data" #set yr [0:100000] binwidth=0.01 bin(x,width)=width*floor(x/width) + binwidth/2.0 plot 'totalmissing' using (bin($1,binwidth)):(1.0) smooth freq with boxes pause -1 EOF ``` ## Missing data por loci --max-missing Excluye sitios en proporción del porcentaje de datos faltantes, va de 0 a 1, en donde cero permite sitios con puros datos faltantes y 1 sitios sin datos faltantes ``` vcftools --vcf populations.snps.vcf --max-missing 0.8 --recode --out populations_08.vcf ``` ### ¿Con qué porcentaje de datos faltantes por loci nos quedamos? ### ¿Qué individuos eliminamos de nuestro archivo .vcf? ## Eliminar individuos de un archivo .vcf Nuevamente utilizaremos vcftools: Crea un archivo de texto que contenga los nombres de los individuos que deseas eliminar (cada individuo en una línea). ej. individuals_to_remove.txt: Sample_1 Sample_2 Sample_3 ``` vcftools --vcf positions.snps.vcf --remove individuals_to_remove.txt --recode --out popultaions_filtered.vcf ``` ## Extraer loci en un nuevo .vcf Crear un archivo de texto con las posiciones de los loci que se desean extraer. ej. loci.txt #CHROM POS chr1 1000 chr2 500 chr3 2000 ``` vcftools --vcf position.snps.vcf --positions loci.txt --recode --out output.vcf ``` ## Filtrar tu archivo por profundidad de genotipo El código para aplicar el filtro en vcftools es el siguiente: ``` vcftools --geno-depth --vcf populations.snps.vcf --out GD_Genotypes ``` #### Este comando extrae la información de la profundidad de cobertura genotípica de un archivo .vcf y la guarda en un archivo llamado "GD_Genotypes_global.gdepth". El flag --geno-depth indica que deseas obtener la información de profundidad de cobertura genotípica. ``` cat GD_Genotypes.gdepth | head -10 ``` --minDP <float> # incluye sólo genotipos iguales o mayores al valor --maxDP <float> # incluye sólo genotipos iguales o menores al valor ``` vcftools --vcf populations.snps.vcf --minDP 3 --out populations_DP3.vcf ``` ## Filtrado de profundidad utilizando R Parece que el archivo de salida de Stacks no es totalmente compatible con vcftools para aplicar este comando. Así que lo podemos hacer con R utilizando las siguientes paqueterías: https://devonderaad.github.io/SNPfiltR/index.html https://knausb.github.io/vcfR_documentation/depth_plot.html ``` #instalar librerías install.packages("SNPfiltR") install.packages("vcfR") install.packages("ggplot2") # llamar librerías library(SNPfiltR) library(vcfR) library(ggplot2) ``` ``` # importar archivo vcf vcf <- read.vcfR(populations_10MDloci.vcf) # ver objeto vcf vcf ``` > vcf ***** Object of Class vcfR ***** 108 samples 2812 CHROMs 2,812 variants Object size: 21.1 Mb 0 percent missing data ***** ***** ***** ### Evaluación de profundidad ``` hard_filter(vcf) ``` ![](https://hackmd.io/_uploads/HJ-JoJ3u3.png) ![](https://hackmd.io/_uploads/rJ3Dok3O2.png) ![](https://hackmd.io/_uploads/SkOYo12On.png) ![](https://hackmd.io/_uploads/Hya5jknun.png) ``` hard_filter(vcf=vcf, depth = 3) hard_filter(vcf=vcf, depth = 4) hard_filter(vcf=vcf, depth = 5) ``` ## Estimación de la profundidad por genotipo ``` dp <- extract.gt(vcf, element = "DP", as.numeric=TRUE) df <- as.data.frame(table(dp)) total_count <- sum(df$Freq) # Recuento total de todos los valores df$Percentage <- df$Freq / total_count * 100 # Calcula el porcentaje ``` ## Gráfica de profundidades (DP) ``` df_plot <- df[c(1:50),] ggplot(df_plot, aes(dp, Percentage)) + geom_point() + labs(title='Curso_datos_DP3 (DP)', x='Genotypes depth coverage genotypes (DP)', y='Percentage') ``` ![](https://hackmd.io/_uploads/B1hcCJ3On.png) ### DP promedio ``` mean(dp, na.rm = TRUE) # mean dp= 29.36998 ``` ### ¿Filtratían su archivo por profundidad? ### ¿si, no, por qué y a que valor? ## Si desean filtrar la base a X valor (ej. DP>6) ``` vcf2 <- hard_filter(vcf=vcf, depth = 6) dp2 <- extract.gt(vcf2, element = "DP", as.numeric=TRUE) df2<- as.data.frame(table(dp2)) total_count2 <- sum(df2$Freq) df2$Percentage <- df2$Freq / total_count2 * 100 df_plot2 <- df2[c(1:50),] ggplot(df_plot2, aes(dp2, Percentage)) + geom_point() + labs(title='Global_DP (DP>3)', x='Genotypes depth coverage genotypes (DP)', y='Percentage') ``` ![](https://hackmd.io/_uploads/r1NARE6uh.png) Escribir un .vcf nuevo con los genotipos filtrados (>6) ``` vcfR::write.vcf(vcf2, file = "Genotypes_DP6.vcf.gz") ```