fixed fai index

README.md

@@ -31,6 +31,10 @@
     SampleName : your sample ID
     fq1: fastq file for a given sample
     fq2: read 2 for paired-end reads
+    (sp1   ra_R1.fastq.gz   ra_R2.fastq.gz)
+
+    a sample my have multiple fastq files separated by a ','
+    (sp1   ra_R1.fastq.gz,rb_R1.fastq.gz   ra_R2.fastq.gz,rb_R2.fastq.gz)
 
     - samples.csv
 
@@ -51,6 +55,13 @@
 
 
 ### 5) optionaly run post vcf filtering (example use for GWAS)
+    1) filter by read depth
+    2) keep only biallelic allels
+    3) max missing sample for a site
+    4) filter by MAF
+    5) keep only polymorph sites
+
+
 `sbatch run_post_filters_and_stats.sh`
 
 #### Documentation being written (still)

Snakefile_para_gatk_PEmerged.smk

@@ -384,9 +384,9 @@ rule multiqc_bam:
 # 3-1) create a reference dictionnary the ref must not be a sl and a ref.fa => ref.dict (not ref.fa.dict)
 rule gatk4_ref_dict:
     input:
-        ref = config["REFPATH"] + "/" + config["GENOME"]
+        ref  = "{refp}/{ref}".format(refp=config["REFPATH"],ref=config["GENOME"]),
+        fai  = "{refp}/{ref}".format(refp=config["REFPATH"],ref=config["GENOME_FAI"]),
     output:
-        fai  = "{refp}/{ref}.fai".format(refp=config["REFPATH"],ref=config["GENOME"]),
         dic  = "{refp}/{ref}.dict".format(refp=config["REFPATH"],ref=config["GENOME"]) #config["REFPATH"] + "/" + config["GENOME"] +".dict"
     params:
         bind = config["BIND"],
@@ -395,7 +395,7 @@ rule gatk4_ref_dict:
     shell:
         """
         singularity exec {params.bind} {params.gatk4_bin} gatk CreateSequenceDictionary -R {input.ref}
-        singularity exec {params.bind} {params.samtools_bin} samtools faidx {input.ref}
+        #singularity exec {params.bind} {params.samtools_bin} samtools faidx {input.ref}
         singularity exec {params.bind} {params.samtools_bin} samtools dict -o {output.dic} {input.ref}
         """
 
@@ -404,8 +404,8 @@ rule gatk4_ref_dict:
 rule gatk4_hc:
     input:
         #"{outdir}/variant/chrslist.OK".format(outdir=config["outdir"]),
-        bam = "{outdir}/mapped/{{sample}}_sorted.bam".format(outdir=config["outdir"]),
-        refdict = config["REFPATH"] + "/" + config["GENOME"] +".dict",
+        bam     = "{outdir}/mapped/{{sample}}_sorted.bam".format(outdir=config["outdir"]),
+        refdict = "{refp}/{ref}.dict".format(refp=config["REFPATH"],ref=config["GENOME"]),
     output:
         gvcf="{outdir}/variant/gatk_gvcf/{{sample}}-{{mychr}}.g.vcf.gz".format(outdir=config["outdir"])
     params:
@@ -430,7 +430,7 @@ rule gatk4_hc:
 rule gatk4_gvcf_map_file:
     input:
         gvcfs = expand("{outdir}/variant/gatk_gvcf/{sample}-{{mychr}}.g.vcf.gz", outdir=config["outdir"], sample=samples['SampleName']),
-        fai       = "{refp}/{ref}.fai".format(refp=config["REFPATH"],ref=config["GENOME"])
+        fai  = "{refp}/{ref}".format(refp=config["REFPATH"],ref=config["GENOME_FAI"]),
     output:
         gvcfmap      = "{outdir}/variant/gvcf_{{mychr}}_list.map".format(outdir=config["outdir"]),
     params: