User Tools

Site Tools


genetica:pre_1kg

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

genetica:pre_1kg [2013/05/03 08:30]
osotolongo [DUPLICATE MARKERS FOUND]
genetica:pre_1kg [2020/08/04 10:58]
Line 1: Line 1:
-{{ ../dna_orbit_animated.gif?100}} 
-====== Convertir las bases de datos de 1000Genome a formato plink ====== 
  
-===== Bajando los archivos ===== 
- 
-Primeramente hay que obtener las bases de datos del proyecto //1000Genome//. Estos archivos estan en formato VCF (Variant Call Format) en el ftp de //1000Genome//, 
- 
-ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20110521/ 
- 
-Observese que esta es la última //release// de los archivos pero hay otras. La forma más fácil de bajar todos los archivos es hacer un **wget** para todos los ficheros. Hacemos la lista primero, 
- 
-<code txt files2download.txt> 
-ALL.chr1.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz 
-ALL.chr1.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz.tbi 
-ALL.chr10.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz 
-ALL.chr10.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz.tbi 
-ALL.chr11.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz 
-ALL.chr11.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz.tbi 
-ALL.chr12.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz 
-ALL.chr12.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz.tbi 
-ALL.chr13.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz 
-ALL.chr13.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz.tbi 
-ALL.chr14.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz 
-ALL.chr14.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz.tbi 
-ALL.chr15.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz 
-ALL.chr15.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz.tbi 
-ALL.chr16.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz 
-ALL.chr16.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz.tbi 
-ALL.chr17.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz 
-ALL.chr17.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz.tbi 
-ALL.chr18.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz 
-ALL.chr18.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz.tbi 
-ALL.chr19.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz 
-ALL.chr19.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz.tbi 
-ALL.chr2.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz 
-ALL.chr2.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz.tbi 
-ALL.chr20.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz 
-ALL.chr20.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz.tbi 
-ALL.chr21.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz 
-ALL.chr21.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz.tbi 
-ALL.chr22.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz 
-ALL.chr22.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz.tbi 
-ALL.chr3.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz 
-ALL.chr3.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz.tbi 
-ALL.chr4.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz 
-ALL.chr4.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz.tbi 
-ALL.chr5.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz 
-ALL.chr5.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz.tbi 
-ALL.chr6.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz 
-ALL.chr6.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz.tbi 
-ALL.chr7.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz 
-ALL.chr7.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz.tbi 
-ALL.chr8.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz 
-ALL.chr8.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz.tbi 
-ALL.chr9.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz 
-ALL.chr9.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz.tbi 
-</code> 
-**Nota:** Los archivos .tbi son los indices que usan las herramientas de procesamiento de VCF para poder leer mas rapidos las DBs. Los he puesto aqui para que no sorprenda su existencia. 
- 
-y despues la bajamos, 
- 
-<code bash> 
-alias fget='wget -U "Mozilla/5.0 (X11; Linux x86_64; rv:10.0.11) Gecko/20121121 Firefox/10.0.11"'; 
-for x in `cat files2download.txt`; do fget ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20110521/$x ; done  
-</code> 
- 
-===== Convirtiendo con vcf_to_ped_converter.pl ===== 
- 
-Una vez que tengamos los archivos vamos a buscar los conversores necesarios. La descripcion de como se convierte puede verse en http://www.1000genomes.org/faq/can-i-convert-vcf-files-plinkped-format 
- 
-Hay dos scripts que son necesarios //vcf_to_ped_converter.pl// y //tabix//. El primero es un script de perl cuya ultima version puede obtenerse de, ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/browser/vcf_to_ped_converter/ 
- 
-El segundo es un paquete que forma parte de las [[http://samtools.sourceforge.net/|SAM Tools]] y que puede bajarse de http://sourceforge.net/projects/samtools/files/tabix/. La pagina de manual esa en [[http://samtools.sourceforge.net/tabix.shtml]]. Una vez bajado y **compilado** se debe configurar el //PATH// del usuario para que se tenga acceso a el. Todo esto está hecho en **detritus** y lo que faltaria por usuario es añadir el //PATH// al archivo //~/.bash_profile//. Algo así, 
- 
-<code bash> 
-PATH=$PATH:/opt/tabix 
-</code> 
- 
-Cuando se intenta convertir los archivos que se han bajado da un error, 
-<code bash> 
-$ /opt/gntics/vcf_to_ped_convert.pl -vcf /media/1000Genome/ALL.chr1.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz -sample_panel_file /media/1000Genome/phase1_integrated_calls.20101123.ALL.panel -population CEU -region 1:10583-249239465 
-[tabix] the index file either does not exist or is older than the vcf file. Please reindex. 
-tabix exited with status 1 at /opt/gntics/vcf_to_ped_convert.pl line 198. 
-</code> 
-Lo que esta diciendo es que hay que reindexar el archivo. Esto se hace con, 
-<code bash> 
-tabix -p vcf /media/1000Genome/ALL.chr1.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz 
-</code> 
-Despues de esto ya es posible convertir el archivo a formato //ped//, 
-<code bash> 
-$ /opt/gntics/vcf_to_ped_convert.pl -vcf /media/1000Genome/ALL.chr1.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz -sample_panel_file /media/1000Genome/phase1_integrated_calls.20101123.ALL.panel -population CEU -population TSI -population IBS -region 1:10583-249239465 
-Created 1_10583-249239465.info and 1_10583-249239465.ped 
-</code> 
-Aqui ya casi hemos terminado. Falta sólo convertir el archivo //.info// a un archivo //.map//. Esto se hace, 
-<code bash> 
-awk {'print "1 "$1" 0 "$2'} 1_10583-249239465.info > 1_10583-249239465.map 
-</code> 
-==== batching ==== 
-Reindex: 
-<code bash> 
-for x in `cat files2download.txt | grep -v tbi`; do tabix -f -p vcf $x; done 
-</code> 
-Convert: 
-<code bash> 
-for x in `cat files2download.txt | grep -v tbi`;  
-do 
-chr=$(echo $x | sed 's/.*chr\(.*\)\.phase.*/\1/'); 
-p0=$(zcat $x | grep "^$chr" | head -n 1 | awk {'print $2'}); 
-pf=$(zcat $x | tail -n 1 | awk {'print $2'}); 
-/opt/gntics/vcf_to_ped_convert.pl -vcf $x -sample_panel_file phase1_integrated_calls.20101123.ALL.panel -population CEU -population TSI -population IBS -region $chr:$p0-$pf; 
-done; 
-</code> 
-info2map: 
-<code bash> 
-for x in *.info; 
-do 
-chr=$(echo $x | awk -F"_" {'print $1'}); 
-awk -v chr="$chr" {'print chr" "$1" 0 "$2'} $x > ${x%.info}.map; 
-done; 
-</code> 
-===== Convirtiendo con VCFtools ===== 
- 
-[[http://vcftools.sourceforge.net/index.html|VCFtools]] es una herramienta para manipular los archivos VCF. Entre otras cosas permite convertir a otros formatos tal y como se indica en la [[http://vcftools.sourceforge.net/options.html|documentacion]]. Tras bajar y **compilar** este paquete ha de añadirse al //PATH// del usuario en //~/.bash_profile//, 
- 
-<code bash> 
-PATH=$PATH:/opt/vcftools_0.1.10/bin 
-</code> 
- 
-Despues ejecutamos la conversión a //tped//, (la conversion a //ped// da problemas) 
-<code bash> 
-vcftools --gzvcf ALL.chr9.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz --plink-tped --out chr9 
-</code> 
- 
-y hacemos la traspuesta al mismo tiempo que el binario 
- 
-<code bash> 
-plink --tped chr9.tped --tfam chr9.tfam --alleleACGT --make-bed --out chr9_ok 
-</code> 
- 
-**Ventajas** 
-   * Puede convertirse **todo** el archivo vcf, no hay que especificar la poblacion. 
-   * Tampoco hay que especificar las regiones a convertir. 
-**Desventajas** 
-   * No hay una forma clara de especificar la poblacion //target// por lo que en principio hay que procesar todo junto (La cosa va a tardar bastante mas) 
- 
-==== batching ==== 
- 
-Para correr esto por todos los archivos, 
- 
-<code bash> 
-for x in `cat files2download.txt | grep -v tbi`;  
-do 
-chr=$(echo $x | sed 's/.*chr\(.*\)\.phase.*/\1/'); 
-vcftools --gzvcf $x --plink-tped --out chr${chr}_tmp; 
-plink --tped chr${chr}_tmp.tped --tfam chr${chr}_tmp.tfam --alleleACGT --make-bed --out all_chr${chr}; 
-rm -rf chr${chr}_tmp.*; 
-done; 
-</code> 
- 
-===== Preprocesando con plink ===== 
- 
-Primero hay que recodificar los alelos a ACGT y de paso los pasamos a binario (solo es necesario si hemos converitdo con el primer metodo) 
- 
-<code bash> 
-for x in *.ped; 
-do 
-plink --file ${x%.ped} --out ${x%.ped} --alleleACGT --make-bed; 
-done; 
-</code> 
- 
-Vamos a intentar juntarlos todos 
- 
-<code bash> 
-afr=(`ls *.bed`); 
-for (( i=1; i<${#afr[@]}; i++ )); 
-do 
-x=${afr[$i]}; 
-echo "${x} ${x%.bed}.bim ${x%.bed}.fam" >> allfiles.txt; 
-done; 
-x=${afr[0]}; 
-plink --bfile ${x%.bed} --merge-list allfiles.txt --make-bed --out 1000genome_CEU_merged 
-</code> 
- 
-===== y ya ta ===== 
- 
-Ahora hay que seguir fundiendo la DB que tenemos con nuestros datos: [[plink_1kg_impute| para poder imputar]] 
- 
-===== porqueria que puede pasar ====== 
- 
-==== DUPLICATE MARKERS FOUND ==== 
-<code> 
-$ vcftools --gzvcf /media/1000Genome/ALL.chr22.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz --plink-tped --out chr22_tmp && plink --tped chr22_tmp.tped --tfam chr22_tmp.tfam --alleleACGT --make-bed --out all_chr22 && rm -rf chr22_tmp.* 
- 
-VCFtools - v0.1.10 
-(C) Adam Auton 2009 
- 
-Parameters as interpreted: 
- --gzvcf /media/1000Genome/ALL.chr22.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz 
- --out chr22_tmp 
- --plink-tped 
- 
-Using zlib version: 1.2.3 
-Versions of zlib >= 1.2.4 will be *much* faster when reading zipped VCF files. 
-Reading Index file. 
-File contains 494328 entries and 1092 individuals. 
-Applying Required Filters. 
-After filtering, kept 1092 out of 1092 Individuals 
-After filtering, kept 494328 out of a possible 494328 Sites 
-Writing PLINK TPED file ... Writing PLINK TFAM file ... Done. 
-Run Time = 1231.00 seconds 
- 
-@----------------------------------------------------------@ 
-|        PLINK!           v1.07      |   10/Aug/2009     | 
-|----------------------------------------------------------| 
-|  (C) 2009 Shaun Purcell, GNU General Public License, v2  | 
-|----------------------------------------------------------| 
-|  For documentation, citation & bug-report instructions:  | 
-|        http://pngu.mgh.harvard.edu/purcell/plink/        | 
-@----------------------------------------------------------@ 
- 
-Web-based version check ( --noweb to skip ) 
-Connecting to web...  OK, v1.07 is current 
- 
-Writing this text to log file [ all_chr22.log ] 
-Analysis started: Thu May  2 09:37:16 2013 
- 
-Options in effect: 
- --tped chr22_tmp.tped 
- --tfam chr22_tmp.tfam 
- --alleleACGT 
- --make-bed 
- --out all_chr22 
- 
-Reading pedigree information from [ chr22_tmp.tfam ]  
-1092 individuals read from [ chr22_tmp.tfam ]  
-0 individuals with nonmissing phenotypes 
-Assuming a disease phenotype (1=unaff, 2=aff, 0=miss) 
-Missing phenotype value is also -9 
-0 cases, 0 controls and 1092 missing 
-0 males, 0 females, and 1092 of unspecified sex 
-Warning, found 1092 individuals with ambiguous sex codes 
-Writing list of these individuals to [ all_chr22.nosex ] 
-494328 (of 494328) markers to be included from [ chr22_tmp.tped ] 
- *** WARNING *** DUPLICATE MARKERS FOUND ***  
-Duplicate marker name found: [ rs11457237 ] 
-Duplicate marker name found: [ rs113940759 ] 
-Duplicate marker name found: [ rs71904485 ] 
- 
-Before frequency and genotyping pruning, there are 494328 SNPs 
-1092 founders and 0 non-founders found 
-Total genotyping rate in remaining individuals is 1 
-0 SNPs failed missingness test ( GENO > 1 ) 
-0 SNPs failed frequency test ( MAF < 0 ) 
-After frequency and genotyping pruning, there are 494328 SNPs 
-After filtering, 0 cases, 0 controls and 1092 missing 
-After filtering, 0 males, 0 females, and 1092 of unspecified sex 
-Writing pedigree information to [ all_chr22.fam ]  
-Writing map (extended format) information to [ all_chr22.bim ]  
-Writing genotype bitfile to [ all_chr22.bed ]  
-Using (default) SNP-major mode 
- 
-Analysis finished: Thu May  2 09:41:48 2013 
- 
-</code> 
- 
-asi que voy y lo convierto a ascii. 
- 
-<code bash> 
-$ plink --bfile all_chr22 --recode --out all_chr22 
-</code> 
- 
-y despues edito a mano los marcadores dobles y pongo a uno //a// y al otro //b//. 
-Ejemplo: 
-<code> 
-22 rs11457237a 0 34030843 
-22 rs11457237b 0 34030846 
-</code> 
-<del> Ahora deja ver si son el mismo, 
- 
-<code> 
-$ plink --file all_chr22 --twolocus rs11457237a rs11457237b --allow-no-sex 
-</code></del> 
- 
-Bueno esta parte no la entendi asi que voy a borrar uno de los duplicados utilizando el criterio loreal (//Because I'm worth it//). 
- 
-<code> 
-$ cat rmsnps.txt 
-rs11457237b 
-rs113940759b 
-rs71904485b 
-</code> 
- 
-<code> 
-$ plink --file all_chr22 --exclude rmsnps.txt --out all_chr22_tmp 
-</code> 
genetica/pre_1kg.txt · Last modified: 2020/08/04 10:58 (external edit)