Differences

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

--- genetica:wgs [2020/05/19 09:57] – [Paralelizacion] osotolongo
+++ genetica:wgs [2020/08/04 10:58] (current) – external edit 127.0.0.1
@@ Line 3: / Line 3: @@
 El objetivo es definir un procedimiento que procese un numero grande de secuencias WGS en el menor tiempo posible. Para ello,una vez definido el pipeline deberemos automatizar las tareas e integrarlas en el schedule manager del cluster ([[:cluster|Como usar el cluster sin morir en el intento]]).
+** Nota: Aunque el script de ejecución final está integrado en el cluster de Fundació ACE, puede ser modificado fácilmente para integrarlo en cualquier otro cluster que use [[https://slurm.schedmd.com/|SLURM]] y probablemente pueda modificarse para otros //schedule managers//. **
 ===== tl;dr =====
 <code>
@@ Line 52: / Line 53: @@
 {{ :genetica:flujo_wgs.png |}}
+Entonces, una vez definido el pipeline, la secuencia de ejecucion de cada trozo y el modo de paralelizacion seria sencillo definir el flujo de ejecucion si supieramos que archivos iniciales debe cargar cada ejecucion inicial (los //bwa//).
 ===== Parsing =====
+Lo primero que ddebemos saber es que el directorio de //input// esta compuesto por una lista de subdirectorios cada uno perteneciente a un sujeto distinto.
+<code bash>
+[osotolongo@detritus HUMehbE]$ pwd
+/the_dysk/BGI_exome/F18FTSEUET0180/HUMehbE
+[osotolongo@detritus HUMehbE]$ ls
+BGI_Data_List_F18FTSEUET0180_filled.pdf                    md5sum.check  seq-1   seq-11  seq-13  seq-15  seq-17  seq-19  seq-20  seq-22  seq-24  seq-26  seq-28  seq-4  seq-6  seq-8
+BGI_Sequencing_Analysis_Report_F18FTSEUET0180_HUMehbE.pdf  md5sum.txt    seq-10  seq-12  seq-14  seq-16  seq-18  seq-2   seq-21  seq-23  seq-25  seq-27  seq-3   seq-5  seq-7  seq-9
+</code>
+Para cada sujeto hay una lista de ocho archivos que deben parearse segun el nombre de archivo.
+<code>
+[osotolongo@detritus HUMehbE]$ ls seq-5/
+V300016291_L01_549_1.fq.gz  V300016291_L01_550_1.fq.gz  V300016291_L01_551_1.fq.gz  V300016291_L01_552_1.fq.gz
+V300016291_L01_549_2.fq.gz  V300016291_L01_550_2.fq.gz  V300016291_L01_551_2.fq.gz  V300016291_L01_552_2.fq.gz
+</code>
+En este ejemplo deben parearse,
+  * V300016291_L01_549_1.fq.gz y V300016291_L01_549_2.fq.gz
+  * V300016291_L01_550_1.fq.gz y V300016291_L01_550_2.fq.gz
+  * V300016291_L01_551_1.fq.gz y V300016291_L01_551_2.fq.gz
+  * V300016291_L01_552_1.fq.gz y V300016291_L01_552_2.fq.gz
+asi que la estructura a seguir es bastante clara. No obstante a que los nombres de archivo varian entre sujetos, siguen la misma estructura.
+<code>
+[osotolongo@detritus HUMehbE]$ ls seq-10
+V300016281_L01_553_1.fq.gz  V300016281_L01_554_1.fq.gz  V300016281_L01_555_1.fq.gz  V300016281_L01_556_1.fq.gz
+V300016281_L01_553_2.fq.gz  V300016281_L01_554_2.fq.gz  V300016281_L01_555_2.fq.gz  V300016281_L01_556_2.fq.gz
+</code>
+Lo que voy hacer entonces es definir un //hash// donde se guarda la informacion relativa a cada sujeto (incluidos los patrones en los archivos originales) y a partir de ahi construir los scripts.
+Veamos, primero leo el //input dir//
+<code perl>
+my $src_dir = shift;
+$src_dir =~ s/\/$//;
+my $find_rule = File::Find::Rule->new;
+my @ppaths = $find_rule->maxdepth(1)->directory->in($src_dir);
+@ppaths = grep {!/$src_dir$/} @ppaths;
+my %ltpaths = map { /.*\/(.*)?$/ => $_ } @ppaths;
+my %finfo;
+my %lpaths;
+</code>
+luego, en dependencia de las opciones puedo reducir la lista o no,
+<code perl>
+if($cfile && -e $cfile && -f $cfile){
+        my @cuts = read_file $cfile;
+        chomp @cuts;
+        foreach my $cut (@cuts){
+                if(grep {/$cut/} %ltpaths){
+                        $lpaths{$cut} = $ltpaths{$cut};
+                }
+        }
+}else{
+        %lpaths = %ltpaths;
+}
+</code>
+y ahora lleno el //hash// de datos,
+<code perl>
+foreach my $pollo (sort keys %lpaths){
+        opendir PD, $lpaths{$pollo} or next;
+        my @fqs = grep { !/^\./} readdir PD;
+        #$finfo{$pollo}{'name'} = $pollo;
+        $finfo{$pollo}{'path'} = $lpaths{$pollo};
+        foreach my $fq (@fqs) {
+                $_=$fq;
+                my @ids  = /^V(\d*?)\_L(\d*?)\_(\d*?)\_\d\.fq\.gz$/;
+                $finfo{$pollo}{'fq_name'} = 'V'.$ids[0].'_L'.$ids[1];
+                push @{ $finfo{$pollo}{'fq_list'} }, $ids[2];
+        }
+}
+</code>
+===== Programatic tree =====
+++++ Ya casi esta, ahora por cada sujeto guardado pueden definirse los primeros 4 procesos,|
+<code perl>
+        for (my $i=0; $i<4; $i++){
+                my $fqid = $finfo{$pollo}{'fq_list'}[2*$i];
+                my $orderfile = $outdir.'/bwa_'.$pollo.'_'.$i.'.sh';
+                open ORD, ">$orderfile" or die "Couldnt create file";
+                print ORD '#!/bin/bash'."\n";
+                print ORD '#SBATCH -J sam_'.$pollo."\n";
+                print ORD '#SBATCH --time=24:0:0'."\n";
+                print ORD '#SBATCH -o '.$outdir.'/bwa_'.$pollo.'-%j'."\n";
+                print ORD '#SBATCH -c 8'."\n";
+                print ORD '#SBATCH --mem-per-cpu=4G'."\n";
+                print ORD '#SBATCH -p fast'."\n";
+                print ORD '#SBATCH --mail-type=FAIL,TIME_LIMIT,STAGE_OUT'."\n"; #no quieres que te mande email de todo
+                print ORD '#SBATCH --mail-user='."$ENV{'USER'}\n";
+                print ORD $bwa.' mem -t 4 -R "@RG\tID:'.$finfo{$pollo}{'fq_name'}.'_'.$fqid.'\tSM:'.$pname.'\tPL:BGI\tPI:380" -M '.$ref_dir.'/Homo_sapiens_assembly38 '.$finfo{$pollo}{'path'}.'/'.$finfo{$pollo}{'fq_name'}.'_'.$fqid.'_1.fq.gz '.$finfo{$pollo}{'path'}.'/'.$finfo{$pollo}{'fq_name'}.'_'.$fqid.'_2.fq.gz > '.$tmpdir.'/'.$pname.'_'.$i.'.sam';
+                print ORD "\n";
+                close ORD;
+                $gsconv.= 'I='.$tmpdir.'/'.$pname.'_'.$i.'.sam ';
+                system("sbatch $orderfile");
+        }
+</code>
+++++
+++++ El numero 5 se hace que dependa de estos 4,|
+<code perl>
+        my $orderfile = $outdir.'/merge_'.$pollo.'.sh';
+        open ORD, ">$orderfile";
+        print ORD '#!/bin/bash'."\n";
+        print ORD '#SBATCH -J sam_'.$pollo."\n";
+        print ORD '#SBATCH -p fast'."\n";
+        print ORD '#SBATCH -o '.$outdir.'/merge_'.$pollo.'-%j'."\n";
+        print ORD '#SBATCH -c 8'."\n";
+        print ORD '#SBATCH --mem-per-cpu=4G'."\n";
+        print ORD '#SBATCH -p fast'."\n";
+        print ORD '#SBATCH --mail-type=FAIL,TIME_LIMIT,STAGE_OUT'."\n"; #no quieres que te mande email de todo
+        print ORD '#SBATCH --mail-user='."$ENV{'USER'}\n";
+        print ORD $picard.' MergeSamFiles '.$gsconv.' O='.$tmpdir.'/'.$pname.'.sam'."\n";
+        print ORD $picard.' SortSam I='.$tmpdir.'/'.$pname.'.sam O='.$tmpdir.'/'.$pname.'_sorted.bam SORT_ORDER=coordinate'."\n";
+        print ORD $samtools.' index '.$tmpdir.'/'.$pname.'_sorted.bam'."\n";
+        close ORD;
+        my $order = 'sbatch --dependency=singleton '.$orderfile;
+        my $ujobid = `$order`;
+        $ujobid = ( split ' ', $ujobid )[ -1 ];
+</code>
+++++
+++++ El 6 y el 7 dependen del 5,|
+<code perl>
+        $orderfile = $outdir.'/verify_'.$pollo.'.sh';
+        open ORD, ">$orderfile";
+        print ORD '#!/bin/bash'."\n";
+        print ORD '#SBATCH -J verify_'.$pollo."\n";
+        print ORD '#SBATCH -p fast'."\n";
+        print ORD '#SBATCH -o '.$outdir.'/verify_'.$pollo.'-%j'."\n";
+        print ORD '#SBATCH -c 4'."\n";
+        print ORD '#SBATCH --mem-per-cpu=4G'."\n";
+        print ORD '#SBATCH -p fast'."\n";
+        print ORD '#SBATCH --mail-type=FAIL,TIME_LIMIT,STAGE_OUT'."\n"; #no quieres que te mande email de todo
+        print ORD '#SBATCH --mail-user='."$ENV{'USER'}\n";
+        print ORD $verifybamib. ' --vcf '.$ref_dir.'/hapmap_3.3.hg38.vcf.gz --bam '.$tmpdir.'/'.$pname.'_sorted.bam --chip-none --maxDepth 1000 --precise --verbose --ignoreRG --out '.$resdir.'/'.$pname.'_verifybam |& grep -v "Skipping marker"'."\n";
+        close ORD;
+        $order = 'sbatch --dependency=afterok:'.$ujobid.' '.$orderfile;
+        my $jobid = `$order`;
+        $jobid = ( split ' ', $jobid )[ -1 ];
+        push @jobids, $jobid;
+</code>
+<code perl>
+        $orderfile = $outdir.'/validate_'.$pollo.'.sh';
+        open ORD, ">$orderfile";
+        print ORD '#!/bin/bash'."\n";
+        print ORD '#SBATCH -J validate_'.$pollo."\n";
+        print ORD '#SBATCH -p fast'."\n";
+        print ORD '#SBATCH -o '.$outdir.'/validate_'.$pollo.'-%j'."\n";
+        print ORD '#SBATCH -c 12'."\n";
+        print ORD '#SBATCH --mem-per-cpu=4G'."\n";
+        print ORD '#SBATCH -p fast'."\n";
+        print ORD '#SBATCH --mail-type=FAIL,TIME_LIMIT,STAGE_OUT'."\n"; #no quieres que te mande email de todo
+        print ORD '#SBATCH --mail-user='."$ENV{'USER'}\n";
+        print ORD $picard.' ValidateSamFile IGNORE=MATE_NOT_FOUND IGNORE=MISSING_READ_GROUP IGNORE=RECORD_MISSING_READ_GROUP I='.$tmpdir.'/'.$pname.'_sorted.bam'."\n";
+        print ORD $picard.' MarkDuplicates I='.$tmpdir.'/'.$pname.'_sorted.bam O='.$tmpdir.'/'.$pname.'_GATKready.bam METRICS_FILE='.$resdir.'/'.$pname.'_metrics.txt QUIET=true MAX_RECORDS_IN_RAM=2000000 ASSUME_SORTED=TRUE CREATE_INDEX=TRUE'."\n";
+        print ORD $gatk3.' -T DepthOfCoverage -R '.$ref_dir.'/Homo_sapiens_assembly38.fasta -nt 1 -ct 10 -ct 15 -ct 20 -ct 30 --omitDepthOutputAtEachBase --omitIntervalStatistics --omitLocusTable -L '.$ref_dir.'/MGI_Exome_Capture_V5_bis2.bed -I '.$tmpdir.'/'.$pname.'_GATKready.bam -o '.$resdir.'/'.$pname.'_exome_coverage'."\n";
+        print ORD $gatk4.' BaseRecalibrator -I '.$tmpdir.'/'.$pname.'_GATKready.bam -R '.$ref_dir.'/Homo_sapiens_assembly38.fasta --known-sites '.$ref_dir.'/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz --known-sites '.$ref_dir.'/dbsnp_146.hg38.vcf.gz --known-sites '.$ref_dir.'/1000G_phase1.snps.high_confidence.hg38.vcf.gz -O '.$resdir.'/'.$pname.'_recal_data.table1'."\n";
+        print ORD $gatk4.' ApplyBQSR -R '.$ref_dir.'/Homo_sapiens_assembly38.fasta -I '.$tmpdir.'/'.$pname.'_GATKready.bam -bqsr-recal-file '.$resdir.'/'.$pname.'_recal_data.table1 -O '.$resdir.'/'.$pname.'_recal.bam'."\n";
+        print ORD $gatk4.' AnalyzeCovariates -bqsr '.$resdir.'/'.$pname.'_recal_data.table1 --plots '.$resdir.'/'.$pname.'_AnalyzeCovariates.pdf'."\n";
+        print ORD $gatk4.' BaseRecalibrator -I '.$resdir.'/'.$pname.'_recal.bam -R '.$ref_dir.'/Homo_sapiens_assembly38.fasta --known-sites '.$ref_dir.'/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz --known-sites '.$ref_dir.'/dbsnp_146.hg38.vcf.gz --known-sites '.$ref_dir.'/1000G_phase1.snps.high_confidence.hg38.vcf.gz -O '.$resdir.'/'.$pname.'_recal_data.table2'."\n";
+        print ORD $gatk4.' AnalyzeCovariates -before '.$resdir.'/'.$pname.'_recal_data.table1 -after '.$resdir.'/'.$pname.'_recal_data.table2 -plots '.$resdir.'/'.$pname.'_before-after-plots.pdf'."\n";
+        print ORD $gatk4.' HaplotypeCaller -R '.$ref_dir.'/Homo_sapiens_assembly38.fasta -I '.$resdir.'/'.$pname.'_recal.bam -ERC GVCF --dbsnp '.$ref_dir.'/dbsnp_146.hg38.vcf.gz -O '.$resdir.'/'.$pname.'_raw.snps.indels.g.vcf.gz'."\n";
+        print ORD $gatk4_l.' VariantEval -R '.$ref_dir.'/Homo_sapiens_assembly38.fasta -L '.$ref_dir.'/MGI_Exome_Capture_V5_bis2.bed -D '.$ref_dir.'/dbsnp_146.hg38.vcf.gz -O '.$resdir.'/'.$pname.'_eval.gatkreport --eval '.$resdir.'/'.$pname.'_raw.snps.indels.g.vcf.gz'."\n";
+        close ORD;
+        $order = 'sbatch --dependency=afterok:'.$ujobid.' '.$orderfile;
+        $jobid = `$order`;
+        $jobid = ( split ' ', $jobid )[ -1 ];
+        push @jobids, $jobid;
+</code>
+++++
+++++ y por ultimo, hay un proceso que depende de que terminen todos los demas (TODOS) y lo que hace es limpiar los archivos temporales y enviar un email de finalizacion, |
+<code perl>
+my $orderfile = $outdir.'/wgs_end.sh';
+open ORD, ">$orderfile";
+print ORD '#!/bin/bash'."\n";
+print ORD '#SBATCH -J wgs_end'."\n";
+print ORD '#SBATCH --mail-type=END'."\n"; #email cuando termine o falle
+print ORD '#SBATCH --mail-user='."$ENV{'USER'}\n";
+print ORD '#SBATCH -o '.$outdir.'/wgs_end-%j'."\n";
+unless ($debug){
+        print ORD "rm -rf $w_dir/*/tmp\n";
+}else{
+        print ORD ":\n";
+}
+close ORD;
+my $sjobs = join(',',@jobids);
+my $order = 'sbatch --depend=afterok:'.$sjobs.' '.$orderfile;
+print "$order\n";
+exec($order);
+</code>
+++++
+++++ La magia completa aqui, en menos de 200 lineas |
+<code perl wgs.pl>
+#!/usr/bin/perl
+# Copyright 2020 O. Sotolongo <asqwerty@gmail.com>
+#
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+use strict; use warnings;
+use File::Find::Rule;
+use File::Slurp qw(read_file);
+use Data::Dump qw(dump);
+#################################################################
+################### rellenar los paths sistema ##################
+#################################################################
+my $ref_dir = '/the_dysk/BGI_exome/reference';
+my $bwa ='/nas/usr/local/bin/bwa';
+#my $picard = 'java -jar /nas/usr/local/bin/picard.jar';
+my $picard = 'java -Djava.io.tmpdir=/nas/'.$ENV{'USER'}.'/tmp/ -Xmx8g -jar /nas/usr/local/bin/picard.jar';
+my $samtools = '/nas/software/samtools/bin/samtools';
+my $verifybamib = '/nas/usr/local/bin/verifyBamID';
+my $gatk3 = 'java -jar /nas/usr/local/opt/gatk3/GenomeAnalysisTK.jar';
+#my $gatk4 = '/nas/usr/local/opt/gatk4/gatk --java-options "-DGATK_STACKTRACE_ON_USER_EXCEPTION=true -Xmx16G"';
+#my $gatk4_l = '/nas/usr/local/opt/gatk4/gatk --java-options "-DGATK_STACKTRACE_ON_USER_EXCEPTION=true -Xmx16G -Xmx16G -XX:+UseConcMarkSweepGC"';
+my $gatk4 = 'singularity run --cleanenv -B /nas:/nas -B /the_dysk:/the_dysk /usr/local/bin/gatk4.simg gatk --java-options "-DGATK_STACKTRACE_ON_USER_EXCEPTION=true -Xmx16G"';
+my $gatk4_l = 'singularity run --cleanenv -B /nas:/nas -B /the_dysk:/the_dysk /usr/local/bin/gatk4.simg gatk --java-options "-DGATK_STACKTRACE_ON_USER_EXCEPTION=true -Xmx16G -XX:+UseConcMarkSweepGC"';
+#################################################################
+#################################################################
+#################################################################
+my $cfile;
+my $outdatadir;
+my $debug = 0;
+@ARGV = ("-h") unless @ARGV;
+while (@ARGV and $ARGV[0] =~ /^-/) {
+    $_ = shift;
+    last if /^--$/;
+    if (/^-cut/) { $cfile = shift; chomp($cfile);}
+    if (/^-o/) { $outdatadir = shift;}
+    if (/^-g/) { $debug = 1;}
+}
+my $src_dir = shift;
+$src_dir =~ s/\/$//;
+my $find_rule = File::Find::Rule->new;
+my @ppaths = $find_rule->maxdepth(1)->directory->in($src_dir);
+@ppaths = grep {!/$src_dir$/} @ppaths;
+my %ltpaths = map { /.*\/(.*)?$/ => $_ } @ppaths;
+my %finfo;
+my %lpaths;
+if($cfile && -e $cfile && -f $cfile){
+        my @cuts = read_file $cfile;
+        chomp @cuts;
+        foreach my $cut (@cuts){
+                if(grep {/$cut/} %ltpaths){
+                        $lpaths{$cut} = $ltpaths{$cut};
+                }
+        }
+}else{
+        %lpaths = %ltpaths;
+}
+# Creo el entorno de ejecucion
+my $tmp_shit = '/nas/'.$ENV{'USER'}.'/tmp/';
+mkdir $tmp_shit;
+my $wdir;
+if($outdatadir){
+        mkdir $outdatadir;
+        $wdir = $outdatadir;
+}else{
+        $wdir = $ENV{'PWD'};
+}
+my $outdir = $wdir.'/slurm';
+mkdir $outdir;
+#Recopilo la informacion de los archivos de entrada
+foreach my $pollo (sort keys %lpaths){
+        opendir PD, $lpaths{$pollo} or next;
+        my @fqs = grep { !/^\./} readdir PD;
+        #$finfo{$pollo}{'name'} = $pollo;
+        $finfo{$pollo}{'path'} = $lpaths{$pollo};
+        foreach my $fq (@fqs) {
+                $_=$fq;
+                my @ids  = /^V(\d*?)\_L(\d*?)\_(\d*?)\_\d\.fq\.gz$/;
+                $finfo{$pollo}{'fq_name'} = 'V'.$ids[0].'_L'.$ids[1];
+                push @{ $finfo{$pollo}{'fq_list'} }, $ids[2];
+        }
+}
+#Creo y ejecuto los procesos
+my @jobids;
+foreach my $pollo (sort keys %finfo){
+        my $udir = $wdir.'/'.$pollo;
+        mkdir $udir;
+        my $tmpdir = $udir.'/tmp';
+        mkdir $tmpdir;
+        my $resdir = $udir.'/results';
+        mkdir $resdir;
+        my $gsconv = "";
+        my @ujobids;
+        (my $pname = $pollo) =~ s/-//g;
+        for (my $i=0; $i<4; $i++){
+                my $fqid = $finfo{$pollo}{'fq_list'}[2*$i];
+                my $orderfile = $outdir.'/bwa_'.$pollo.'_'.$i.'.sh';
+                open ORD, ">$orderfile" or die "Couldnt create file";
+                print ORD '#!/bin/bash'."\n";
+                print ORD '#SBATCH -J sam_'.$pollo."\n";
+                print ORD '#SBATCH --time=24:0:0'."\n";
+                print ORD '#SBATCH -o '.$outdir.'/bwa_'.$pollo.'-%j'."\n";
+                print ORD '#SBATCH -c 8'."\n";
+                print ORD '#SBATCH --mem-per-cpu=4G'."\n";
+                print ORD '#SBATCH -p fast'."\n";
+                print ORD '#SBATCH --mail-type=FAIL,TIME_LIMIT,STAGE_OUT'."\n"; #no quieres que te mande email de todo
+                print ORD '#SBATCH --mail-user='."$ENV{'USER'}\n";
+                print ORD $bwa.' mem -t 4 -R "@RG\tID:'.$finfo{$pollo}{'fq_name'}.'_'.$fqid.'\tSM:'.$pname.'\tPL:BGI\tPI:380" -M '.$ref_dir.'/Homo_sapiens_assembly38 '.$finfo{$pollo}{'path'}.'/'.$finfo{$pollo}{'fq_name'}.'_'.$fqid.'_1.fq.gz '.$finfo{$pollo}{'path'}.'/'.$finfo{$pollo}{'fq_name'}.'_'.$fqid.'_2.fq.gz > '.$tmpdir.'/'.$pname.'_'.$i.'.sam';
+                print ORD "\n";
+                close ORD;
+                $gsconv.= 'I='.$tmpdir.'/'.$pname.'_'.$i.'.sam ';
+                system("sbatch $orderfile");
+        }
+        my $orderfile = $outdir.'/merge_'.$pollo.'.sh';
+        open ORD, ">$orderfile";
+        print ORD '#!/bin/bash'."\n";
+        print ORD '#SBATCH -J sam_'.$pollo."\n";
+        print ORD '#SBATCH -p fast'."\n";
+        print ORD '#SBATCH -o '.$outdir.'/merge_'.$pollo.'-%j'."\n";
+        print ORD '#SBATCH -c 8'."\n";
+        print ORD '#SBATCH --mem-per-cpu=4G'."\n";
+        print ORD '#SBATCH -p fast'."\n";
+        print ORD '#SBATCH --mail-type=FAIL,TIME_LIMIT,STAGE_OUT'."\n"; #no quieres que te mande email de todo
+        print ORD '#SBATCH --mail-user='."$ENV{'USER'}\n";
+        print ORD $picard.' MergeSamFiles '.$gsconv.' O='.$tmpdir.'/'.$pname.'.sam'."\n";
+        print ORD $picard.' SortSam I='.$tmpdir.'/'.$pname.'.sam O='.$tmpdir.'/'.$pname.'_sorted.bam SORT_ORDER=coordinate'."\n";
+        print ORD $samtools.' index '.$tmpdir.'/'.$pname.'_sorted.bam'."\n";
+        close ORD;
+        my $order = 'sbatch --dependency=singleton '.$orderfile;
+        my $ujobid = `$order`;
+        $ujobid = ( split ' ', $ujobid )[ -1 ];
+        $orderfile = $outdir.'/verify_'.$pollo.'.sh';
+        open ORD, ">$orderfile";
+        print ORD '#!/bin/bash'."\n";
+        print ORD '#SBATCH -J verify_'.$pollo."\n";
+        print ORD '#SBATCH -p fast'."\n";
+        print ORD '#SBATCH -o '.$outdir.'/verify_'.$pollo.'-%j'."\n";
+        print ORD '#SBATCH -c 4'."\n";
+        print ORD '#SBATCH --mem-per-cpu=4G'."\n";
+        print ORD '#SBATCH -p fast'."\n";
+        print ORD '#SBATCH --mail-type=FAIL,TIME_LIMIT,STAGE_OUT'."\n"; #no quieres que te mande email de todo
+        print ORD '#SBATCH --mail-user='."$ENV{'USER'}\n";
+        print ORD $verifybamib. ' --vcf '.$ref_dir.'/hapmap_3.3.hg38.vcf.gz --bam '.$tmpdir.'/'.$pname.'_sorted.bam --chip-none --maxDepth 1000 --precise --verbose --ignoreRG --out '.$resdir.'/'.$pname.'_verifybam |& grep -v "Skipping marker"'."\n";
+        close ORD;
+        $order = 'sbatch --dependency=afterok:'.$ujobid.' '.$orderfile;
+        my $jobid = `$order`;
+        $jobid = ( split ' ', $jobid )[ -1 ];
+        push @jobids, $jobid;
+        $orderfile = $outdir.'/validate_'.$pollo.'.sh';
+        open ORD, ">$orderfile";
+        print ORD '#!/bin/bash'."\n";
+        print ORD '#SBATCH -J validate_'.$pollo."\n";
+        print ORD '#SBATCH -p fast'."\n";
+        print ORD '#SBATCH -o '.$outdir.'/validate_'.$pollo.'-%j'."\n";
+        print ORD '#SBATCH -c 12'."\n";
+        print ORD '#SBATCH --mem-per-cpu=4G'."\n";
+        print ORD '#SBATCH -p fast'."\n";
+        print ORD '#SBATCH --mail-type=FAIL,TIME_LIMIT,STAGE_OUT'."\n"; #no quieres que te mande email de todo
+        print ORD '#SBATCH --mail-user='."$ENV{'USER'}\n";
+        print ORD $picard.' ValidateSamFile IGNORE=MATE_NOT_FOUND IGNORE=MISSING_READ_GROUP IGNORE=RECORD_MISSING_READ_GROUP I='.$tmpdir.'/'.$pname.'_sorted.bam'."\n";
+        print ORD $picard.' MarkDuplicates I='.$tmpdir.'/'.$pname.'_sorted.bam O='.$tmpdir.'/'.$pname.'_GATKready.bam METRICS_FILE='.$resdir.'/'.$pname.'_metrics.txt QUIET=true MAX_RECORDS_IN_RAM=2000000 ASSUME_SORTED=TRUE CREATE_INDEX=TRUE'."\n";
+        print ORD $gatk3.' -T DepthOfCoverage -R '.$ref_dir.'/Homo_sapiens_assembly38.fasta -nt 1 -ct 10 -ct 15 -ct 20 -ct 30 --omitDepthOutputAtEachBase --omitIntervalStatistics --omitLocusTable -L '.$ref_dir.'/MGI_Exome_Capture_V5_bis2.bed -I '.$tmpdir.'/'.$pname.'_GATKready.bam -o '.$resdir.'/'.$pname.'_exome_coverage'."\n";
+        print ORD $gatk4.' BaseRecalibrator -I '.$tmpdir.'/'.$pname.'_GATKready.bam -R '.$ref_dir.'/Homo_sapiens_assembly38.fasta --known-sites '.$ref_dir.'/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz --known-sites '.$ref_dir.'/dbsnp_146.hg38.vcf.gz --known-sites '.$ref_dir.'/1000G_phase1.snps.high_confidence.hg38.vcf.gz -O '.$resdir.'/'.$pname.'_recal_data.table1'."\n";
+        print ORD $gatk4.' ApplyBQSR -R '.$ref_dir.'/Homo_sapiens_assembly38.fasta -I '.$tmpdir.'/'.$pname.'_GATKready.bam -bqsr-recal-file '.$resdir.'/'.$pname.'_recal_data.table1 -O '.$resdir.'/'.$pname.'_recal.bam'."\n";
+        print ORD $gatk4.' AnalyzeCovariates -bqsr '.$resdir.'/'.$pname.'_recal_data.table1 --plots '.$resdir.'/'.$pname.'_AnalyzeCovariates.pdf'."\n";
+        print ORD $gatk4.' BaseRecalibrator -I '.$resdir.'/'.$pname.'_recal.bam -R '.$ref_dir.'/Homo_sapiens_assembly38.fasta --known-sites '.$ref_dir.'/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz --known-sites '.$ref_dir.'/dbsnp_146.hg38.vcf.gz --known-sites '.$ref_dir.'/1000G_phase1.snps.high_confidence.hg38.vcf.gz -O '.$resdir.'/'.$pname.'_recal_data.table2'."\n";
+        print ORD $gatk4.' AnalyzeCovariates -before '.$resdir.'/'.$pname.'_recal_data.table1 -after '.$resdir.'/'.$pname.'_recal_data.table2 -plots '.$resdir.'/'.$pname.'_before-after-plots.pdf'."\n";
+        print ORD $gatk4.' HaplotypeCaller -R '.$ref_dir.'/Homo_sapiens_assembly38.fasta -I '.$resdir.'/'.$pname.'_recal.bam -ERC GVCF --dbsnp '.$ref_dir.'/dbsnp_146.hg38.vcf.gz -O '.$resdir.'/'.$pname.'_raw.snps.indels.g.vcf.gz'."\n";
+        print ORD $gatk4_l.' VariantEval -R '.$ref_dir.'/Homo_sapiens_assembly38.fasta -L '.$ref_dir.'/MGI_Exome_Capture_V5_bis2.bed -D '.$ref_dir.'/dbsnp_146.hg38.vcf.gz -O '.$resdir.'/'.$pname.'_eval.gatkreport --eval '.$resdir.'/'.$pname.'_raw.snps.indels.g.vcf.gz'."\n";
+        close ORD;
+        $order = 'sbatch --dependency=afterok:'.$ujobid.' '.$orderfile;
+        $jobid = `$order`;
+        $jobid = ( split ' ', $jobid )[ -1 ];
+        push @jobids, $jobid;
+}
+my $orderfile = $outdir.'/wgs_end.sh';
+open ORD, ">$orderfile";
+print ORD '#!/bin/bash'."\n";
+print ORD '#SBATCH -J wgs_end'."\n";
+print ORD '#SBATCH --mail-type=END'."\n"; #email cuando termine o falle
+print ORD '#SBATCH --mail-user='."$ENV{'USER'}\n";
+print ORD '#SBATCH -o '.$outdir.'/wgs_end-%j'."\n";
+unless ($debug){
+        print ORD "rm -rf $w_dir/*/tmp\n";
+}else{
+        print ORD ":\n";
+}
+close ORD;
+my $sjobs = join(',',@jobids);
+my $order = 'sbatch --depend=afterok:'.$sjobs.' '.$orderfile;
+print "$order\n";
+exec($order);
+</code>
+++++
+{{:genetica:its_fucking_cool.gif?300|
+-)}}
+===== Ejecucion =====
+La forma correcta de ejecutar el script es,
+<code>
+$ ./wgs.pl -o /the_dysk/BGI_exome/F18FTSEUET0180/WGS /the_dysk/BGI_exome/F18FTSEUET0180/HUMehbE/
+</code>
+La opción //-o// indica al script donde guardar los resultados, en este caso en // /the_dysk/BGI_exome/F18FTSEUET0180/WGS //. Los datos se leeran desde ///the_dysk/BGI_exome/F18FTSEUET0180/HUMehbE/ //,
+La opcion //-o// no es obligatoria. En caso de no suministrar un directorio de salida, esta se escribira en el directorio actual desde donde se ejeuta el script.
+Los resultados se escriben en el directorio de salida, haciendo un subdirectorio para cada sujeto analizado y siguiendo las mismas convenciones del //input//,
+<code>
+[osotolongo@brick03 WGS]$ ls /the_dysk/BGI_exome/F18FTSEUET0180/WGS
+seq-1	seq-11	seq-13	seq-15	seq-17	seq-19	seq-20	seq-22	seq-24	seq-26	seq-28	seq-4  seq-6  seq-8  slurm
+seq-10	seq-12	seq-14	seq-16	seq-18	seq-2	seq-21	seq-23	seq-25	seq-27	seq-3	seq-5  seq-7  seq-9
+</code>
+Dentro de cada sujeto se crea u n directorio //results// con los resultados del analisis,
+<code bash>
+[osotolongo@brick03 WGS]$ tree seq-1
+seq-1
+└── results
+    ├── seq1_AnalyzeCovariates.pdf
+    ├── seq1_before-after-plots.pdf
+    ├── seq1_eval.gatkreport
+    ├── seq1_exome_coverage.sample_statistics
+    ├── seq1_exome_coverage.sample_summary
+    ├── seq1_metrics.txt
+    ├── seq1_raw.snps.indels.g.vcf.gz
+    ├── seq1_raw.snps.indels.g.vcf.gz.tbi
+    ├── seq1_recal.bai
+    ├── seq1_recal.bam
+    ├── seq1_recal_data.table1
+    ├── seq1_recal_data.table2
+    ├── seq1_verifybam.depthSM
+    ├── seq1_verifybam.log
+    └── seq1_verifybam.selfSM
+directory, 15 files
+</code>
+----
+=== trick ===
+Esta estructura permitiria ejecutar,
+<code>
+$ ./wgs.pl -o /the_dysk/BGI_exome/F18FTSEUET0180/HUMehbE /the_dysk/BGI_exome/F18FTSEUET0180/HUMehbE
+</code>
+ y **en caso de no tener problemas con los permisos**, los resultados quedarian en un directorio dentro del mismo sujeto de proyecto. No se ha hecho asi por defecto previendo precisamente los problemas de permisos.
+----
 ===== Troubleshooting =====
+Pueden ocurrir multitud de problemas en la ejecucion del script. Afortunadamente un fallo en uno de las tareas afecta solamente al sujeto individual sobre el que se realiza la tarea. El resto de sujetos se procesara correctamente.
+Al ocurrir un fallo en algun script, el sistema enviara un email informando del fallo y el //jobid// asociado. Toda la informacion ira en el //subject// del email.
+==== ejemplo 1 ====
+Aqui se informa que la tarea con //jobid = 154540// ha fallado. Tambien sabemos que esta relacionada con el sujeto //seq-27//.
+<code>
+SLURM Job_id=154540 Name=sam_seq-27 Failed, Run time 00:28:28, FAILED, ExitCode 1
+</code>
+Ahora, dentro del directorio de output se crea un directorio llamado //slurm// donde se guardan todos los logs de ejecucion de las tareas. Hemos de buscar este fallo en particular por el //jobid//,
+<code>
+[osotolongo@brick03 WGS]$ ls slurm/ | grep 154540
+bwa_seq-27-154540
+</code>
+una vez localizado el //log// correcto podeos buscar el error manualmente,
+<code>
+[osotolongo@brick03 WGS]$ less slurm/bwa_seq-27-154540
+</code>
+o simplemente hacer algun grep,
+<code>
+[osotolongo@brick03 WGS]$ grep -i error slurm/bwa_seq-27-154540
+[gzread] Input/output error
+</code>
+Como informacion adicional, en el directorio slurm tambien se guardan los scripts a ejecutar para cada parte del proceso,
+<code bash>
+[osotolongo@brick03 WGS]$ ls slurm/ | grep seq-27 | grep sh
+bwa_seq-27_0.sh
+bwa_seq-27_1.sh
+bwa_seq-27_2.sh
+bwa_seq-27_3.sh
+merge_seq-27.sh
+validate_seq-27.sh
+verify_seq-27.sh
+</code>
+El nombre del log nos dice que ha fallado alguno de los //bwa//. En caso de querer depurar mas lo que ha pasado podemos leer el log y encontrar cual de los comandos ha fallado exactamente.
+==== ejemplo 2 ====
+Algo menos basico, recibimos un email con subject,
+<code>
+SLURM Job_id=154595 Name=validate_seq-8 Failed, Run time 01:01:18, FAILED, ExitCode 2
+</code>
+y buscamos el log correspondiente,
+<code>
+[osotolongo@brick03 WGS]$ ls slurm/ | grep 154595
+validate_seq-8-154595
+</code>
+Aqui si hacemos,
+<code>
+[osotolongo@brick03 WGS]$ grep -i error slurm/validate_seq-8-154595
+</code>
+obtenemos un monton de informacion ya que el error ocurrio a un nivel y todos los niveles posteriores arrojan errores. Asi que no hay mas remediio que leerse el log con calma, y vemos que mientras se hace **MarkDuplicates** tenemos un error de //Java//,
+<code>
+Exception in thread "main" picard.PicardException: Exception writing ReadEnds to file.
+        at picard.sam.markduplicates.util.ReadEndsForMarkDuplicatesCodec.decode(ReadEndsForMarkDuplicatesCodec.java:110)
+        at picard.sam.markduplicates.util.ReadEndsForMarkDuplicatesCodec.decode(ReadEndsForMarkDuplicatesCodec.java:32)
+        at htsjdk.samtools.util.SortingCollection$FileRecordIterator.advance(SortingCollection.java:630)
+        at htsjdk.samtools.util.SortingCollection$FileRecordIterator.next(SortingCollection.java:620)
+        at htsjdk.samtools.util.PeekIterator.peek(PeekIterator.java:69)
+        at htsjdk.samtools.util.SortingCollection$PeekFileRecordIteratorComparator.compare(SortingCollection.java:655)
+        at htsjdk.samtools.util.SortingCollection$PeekFileRecordIteratorComparator.compare(SortingCollection.java:652)
+        at java.util.TreeMap.compare(TreeMap.java:1295)
+        at java.util.TreeMap.put(TreeMap.java:538)
+        at java.util.TreeSet.add(TreeSet.java:255)
+        at htsjdk.samtools.util.SortingCollection$MergingIterator.next(SortingCollection.java:566)
+        at picard.sam.markduplicates.MarkDuplicates.generateDuplicateIndexes(MarkDuplicates.java:729)
+        at picard.sam.markduplicates.MarkDuplicates.doWork(MarkDuplicates.java:259)
+        at picard.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:305)
+        at picard.cmdline.PicardCommandLine.instanceMain(PicardCommandLine.java:103)
+        at picard.cmdline.PicardCommandLine.main(PicardCommandLine.java:113)
+Caused by: java.io.IOException: failed to read chunk
+        at org.xerial.snappy.SnappyInputStream.hasNextChunk(SnappyInputStream.java:433)
+        at org.xerial.snappy.SnappyInputStream.read(SnappyInputStream.java:167)
+        at java.io.DataInputStream.readFully(DataInputStream.java:195)
+        at java.io.DataInputStream.readLong(DataInputStream.java:416)
+        at picard.sam.markduplicates.util.ReadEndsForMarkDuplicatesCodec.decode(ReadEndsForMarkDuplicatesCodec.java:92)
+        ... 15 more
+</code>
+Este error desencadena todos los demas errores.
+Dado que sabemos que este log corresponde al script //validate_seq-8.sh// podemos encontrar rapidamente el comando culpable del error,
+<code>
+[osotolongo@brick03 WGS]$ grep MarkDuplicates slurm/validate_seq-8.sh
+java -Djava.io.tmpdir=/nas/osotolongo/tmp/ -Xmx8g -jar /nas/usr/local/bin/picard.jar MarkDuplicates I=/the_dysk/BGI_exome/F18FTSEUET0180/WGS/seq-8/tmp/seq8_sorted.bam O=/the_dysk/BGI_exome/F18FTSEUET0180/WGS/seq-8/tmp/seq8_GATKready.bam METRICS_FILE=/the_dysk/BGI_exome/F18FTSEUET0180/WGS/seq-8/tmp/seq8_metrics.txt QUIET=true MAX_RECORDS_IN_RAM=2000000 ASSUME_SORTED=TRUE CREATE_INDEX=TRUE
+</code>
+e intentar reparar lo que este mal.
+==== Opciones extra ====
+=== lista de inclusion ===
+Digamos que de un pool de sujetos han fallado 2 o 3. No queremos tener que ejecutar el analisis para todos los sujetos, solo para los que han fallado. No queremos mover los sujetos que han fallado a un nuevo directorio pues atenta contra la organizacion del trabajo y el uso correcto del tiempo.
+Para esto existe la opcion //-cut// que indica al script que se analice exclusivamente los sujetos contenidos en un archivo de input.
+La forma correcta de utilizarlo es la siguiente. Ejemplo, si solo queremos que se ejecute el analisis en los sujetos //seq-8// y //seq-27//, hacemos el siguiente archivo,
+<code>
+[osotolongo@brick03 WGS]$ cat /nas/osotolongo/wgs/only.txt
+seq-8
+seq-27
+</code>
+y tras esto ejecutamos algo como,
+<code>
+$ ./wgs.pl -cut /nas/osotolongo/wgs/only.txt -o /the_dysk/BGI_exome/F18FTSEUET0180/WGS /the_dysk/BGI_exome/F18FTSEUET0180/HUMehbE
+</code>
+que hara todo el procedimiento pero solo para estos sujetos.
+=== debug ===
+Para ayudar en la localizacion de errores existe la opcion //-g//, para indicar que no se borren los archivos temporales. El comando,
+<code>
+$ ./wgs.pl -g -o /the_dysk/BGI_exome/F18FTSEUET0180/WGS /the_dysk/BGI_exome/F18FTSEUET0180/HUMehbE
+</code>
+deja un subdirectorio //tmp// por cada sujeto con sus archivos temporales correspondientes. Ejemplo para seq-5, estos archivos estarian en // /the_dysk/BGI_exome/F18FTSEUET0180/WGS/seq-5/tmp/ //.