DSL2: genotyping

.github/workflows/ci.yml

@@ -28,11 +28,11 @@ jobs:
           - "latest-everything"
         PARAMS:
           - "-profile test,docker --preprocessing_tool fastp --preprocessing_adapterlist 'https://github.com/nf-core/test-datasets/raw/modules/data/delete_me/fastp/adapters.fasta'"
-          - "-profile test,docker --preprocessing_tool adapterremoval --preprocessing_adapterlist 'https://github.com/nf-core/test-datasets/raw/modules/data/delete_me/adapterremoval/adapterremoval_adapterlist.txt' --sequencing_qc_tool falco"
-          - "-profile test,docker --mapping_tool bwamem --run_mapdamage_rescaling --run_pmd_filtering --run_trim_bam"
-          - "-profile test,docker --mapping_tool bowtie2 --damagecalculation_tool mapdamage --damagecalculation_mapdamage_downsample 100"
+          - "-profile test,docker --preprocessing_tool adapterremoval --preprocessing_adapterlist 'https://github.com/nf-core/test-datasets/raw/modules/data/delete_me/adapterremoval/adapterremoval_adapterlist.txt' --sequencing_qc_tool falco --run_genotyping --genotyping_tool 'freebayes' --genotyping_source 'raw'"
+          - "-profile test,docker --mapping_tool bwamem --run_mapdamage_rescaling --run_pmd_filtering --run_trim_bam --run_genotyping --genotyping_tool 'ug' --genotyping_source 'trimmed'"
+          - "-profile test,docker --mapping_tool bowtie2 --damagecalculation_tool mapdamage --damagecalculation_mapdamage_downsample 100 --run_genotyping --genotyping_tool 'hc' --genotyping_source 'raw'"
           - "-profile test,docker --skip_preprocessing"
-          - "-profile test_humanbam,docker --run_mtnucratio --run_contamination_estimation_angsd --snpcapture_bed 'https://raw.githubusercontent.com/nf-core/test-datasets/eager/reference/Human/1240K.pos.list_hs37d5.0based.bed.gz'"
+          - "-profile test_humanbam,docker --run_mtnucratio --run_contamination_estimation_angsd --snpcapture_bed 'https://raw.githubusercontent.com/nf-core/test-datasets/eager/reference/Human/1240K.pos.list_hs37d5.0based.bed.gz' --run_genotyping --genotyping_tool 'pileupcaller' --genotyping_source 'raw'"
           - "-profile test_multiref,docker" ## TODO add damage manipulation here instead once it goes multiref
     steps:
       - name: Check out pipeline code

CITATIONS.md

@@ -100,10 +100,27 @@
 
 - [QualiMap](https://doi.org/10.1093/bioinformatics/btv566)
 
-  > QualiMap Okonechnikov, K., Conesa, A., & García-Alcalde, F. (2016). Qualimap 2: advanced multi-sample quality control for high-throughput sequencing data. Bioinformatics , 32(2), 292–294. Download: http://qualimap.bioinfo.cipf.es/
+  > QualiMap Okonechnikov, K., Conesa, A., & García-Alcalde, F. (2016). Qualimap 2: advanced multi-sample quality control for high-throughput sequencing data. Bioinformatics , 32(2), 292–294. doi: [10.1093/bioinformatics/btv566](https://doi.org/10.1093/bioinformatics/btv566).
 
 - [DamageProfiler](https://doi.org/10.1093/bioinformatics/btab190)
-  > DamageProfiler Neukamm, J., Peltzer, A., & Nieselt, K. (2020). DamageProfiler: Fast damage pattern calculation for ancient DNA. In Bioinformatics (btab190). doi: [10.1093/bioinformatics/btab190](https://doi.org/10.1093/bioinformatics/btab190). Download: https://github.com/Integrative-Transcriptomics/DamageProfiler
+
+  > DamageProfiler Neukamm, J., Peltzer, A., & Nieselt, K. (2020). DamageProfiler: Fast damage pattern calculation for ancient DNA. In Bioinformatics (btab190). doi: [10.1093/bioinformatics/btab190](https://doi.org/10.1093/bioinformatics/btab190).
+
+- [GATK 3.5](https://console.cloud.google.com/storage/browser/gatk)
+
+  > DePristo M, Banks E, Poplin R, Garimella K, Maguire J, Hartl C, Philippakis A, del Angel G, Rivas MA, Hanna M, McKenna A, Fennell T, Kernytsky A, Sivachenko A, Cibulskis K, Gabriel S, Altshuler D, Daly M. (2011). A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nature Genetics, 43(5), 491–498. doi: [10.1038/ng.806](https://doi.org/10.1038/ng.806).
+
+- [GATK 4.X](https://github.com/broadinstitute/gatk/releases)
+
+  > Poplin R, Ruano-Rubio V, DePristo MA, Fennell TJ, Carneiro MO, Van der Auwera GA, Kling DE, Gauthier LD, Levy-Moonshine A, Roazen D, Shakir K, Thibault J, Chandran S, Whelan C, Lek M, Gabriel S, Daly MJ, Neale B, MacArthur DG, Banks E. (2017). Scaling accurate genetic variant discovery to tens of thousands of samples bioRxiv, 201178. doi: [10.1101/201178](https://doi.org/10.1101/201178).
+
+- [FreeBayes](https://github.com/freebayes/freebayes)
+
+  > Garrison E, Marth G. Haplotype-based variant detection from short-read sequencing. arXiv preprint arXiv:1207.3907 \[q-bio.GN] 2012. doi: [10.48550/arXiv.1207.3907](https://doi.org/10.48550/arXiv.1207.3907).
+
+- [BCFtools](https://github.com/samtools/bcftools)
+
+  > Li H. A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data. Bioinformatics (2011) 27(21) 2987-93.doi: [10.1093/bioinformatics/btr509](https://doi.org/10.1093/bioinformatics/btr509).
 
 ## Software packaging/containerisation tools
 

bin/collect_genotypes.py

@@ -0,0 +1,102 @@
+#!/usr/bin/env python
+
+# MIT License (c) Thiseas C. Lamnidis (@TCLamnidis)
+
+import argparse
+import filecmp
+
+def file_len(fname):
+    with open(fname) as f:
+        for i, l in enumerate(f):
+            pass
+    return i + 1
+
+## A function to return the number of genotypes per line in a .geno file.
+def file_width(fname):
+    with open(fname) as f:
+        for i in f:
+            return(len(i.strip()))
+            break
+
+## A function to check that there are no duplicate individual IDs across ind files.
+def check_for_duplicate_ids(indf1, indf2):
+    with open(indf1) as f:
+        inds1 = [x.strip().split()[0] for x in f.readlines()]
+    with open(indf2) as f:
+        inds2 = [x.strip().split()[0] for x in f.readlines()]
+    intersection = set(inds1).intersection(inds2)
+    if len(intersection) > 0:
+        raise IOError("Input .ind files contain duplicate individual IDs. Duplicates: {}".format(intersection))
+
+## Function to check that the snp files are identical
+def check_snp_files(snpf1, snpf2):
+    if not filecmp.cmp(snpf1, snpf2):
+        raise IOError("Input .snp files are not identical.")
+
+## Function to check the consistency of an eigenstrat database
+def validate_eigenstrat(genof, snpf, indf):
+    dimsGeno = [file_len(genof), file_width(genof)]
+    linesSnp = file_len(snpf)
+    linesInd = file_len(indf)
+
+    # print(dimsGeno,linesSnp,linesInd)
+    ##Check geno and snp compatibility
+    if dimsGeno[0] !=    linesSnp:
+        raise IOError("Input .snp and .geno files do not match.")
+
+    ##Check geno and ind compatibility
+    if dimsGeno[1] !=    linesInd:
+        raise IOError("Input .ind and .geno files do not match.")
+
+VERSION = "1.0.0"
+
+parser = argparse.ArgumentParser(usage="%(prog)s (-i <Input file prefix>) (-c <input ind file> | -R | -E) [-L <SAMPLE LIST> | -S Ind [-S Ind2]] [-o <OUTPUT FILE PREFIX>]" , description="A tool to check two different EingenStrat databses for shared individuals, and extract or remove individuals from an EigenStrat database.")
+parser._optionals.title = "Available options"
+parser.add_argument("-g1", "--genoFn1", type = str, metavar = "<GENO FILE 1 NAME>",    required = True, help = "The path to the input geno file of the first dataset.")
+parser.add_argument("-s1", "--snpFn1",  type = str, metavar = "<SNP FILE 1 NAME>",     required = True, help = "The path to the input snp file of the first dataset.")
+parser.add_argument("-i1", "--indFn1",  type = str, metavar = "<IND FILE 1 NAME>",     required = True, help = "The path to the input ind file of the first dataset.")
+parser.add_argument("-g2", "--genoFn2", type = str, metavar = "<GENO FILE 2 NAME>",    required = True, help = "The path to the input geno file of the second dataset.")
+parser.add_argument("-s2", "--snpFn2",  type = str, metavar = "<SNP FILE 2 NAME>",     required = True, help = "The path to the input snp file of the second dataset.")
+parser.add_argument("-i2", "--indFn2",  type = str, metavar = "<IND FILE 2 NAME>",     required = True, help = "The path to the input ind file of the second dataset.")
+parser.add_argument("-o", "--output",   type = str, metavar = "<OUTPUT FILES PREFIX>", required = True, help = "The desired output file prefix. Three output files are created, <OUTPUT FILES PREFIX>.geno , <OUTPUT FILES PREFIX>.snp and <OUTPUT FILES PREFIX>.ind .")
+parser.add_argument("-v", "--version", action='version', version="{}".format(VERSION), help="Print the version and exit.")
+args = parser.parse_args()
+
+## Open input files
+GenoFile1 = open(args.genoFn1, "r")
+SnpFile1  = open(args.snpFn1,  "r")
+IndFile1  = open(args.indFn1,  "r")
+
+GenoFile2 = open(args.genoFn2, "r")
+# SnpFile2  = open(args.snpFn2,  "r") ## Never actually read in line by line
+IndFile2  = open(args.indFn2,  "r")
+
+## open output files
+GenoFileOut = open(args.output + ".geno", "w")
+SnpFileOut  = open(args.output + ".snp",  "w")
+IndFileOut  = open(args.output + ".ind",  "w")
+
+## Perform basic validation on inputs
+validate_eigenstrat(args.genoFn1, args.snpFn1, args.indFn1)
+validate_eigenstrat(args.genoFn2, args.snpFn2, args.indFn2)
+check_for_duplicate_ids(args.indFn1, args.indFn2)
+check_snp_files(args.snpFn1, args.snpFn2)
+
+## Now actually merge the data
+## Geno
+for line1, line2 in zip(GenoFile1, GenoFile2):
+    geno_line="{}{}".format(line1.strip(),line2.strip())
+    print(geno_line, file=GenoFileOut)
+
+## Snp
+##  Copying the file would be faster, but this way we do not rely on the os or external packages.
+##  We already checked that the snp files are byte-identical, so we can just copy one of them.
+for line in SnpFile1:
+    print(line.strip(), file=SnpFileOut)
+
+## Ind
+##  The indfiles are simply concatenated in the same order as the geno file.
+for line in IndFile1:
+    print(line.strip(), file=IndFileOut)
+for line in IndFile2:
+    print(line.strip(), file=IndFileOut)

conf/modules.config

@@ -961,4 +961,172 @@ process {
             pattern: '*.flagstat'
         ]
     }
+
+    //
+    // GENOTYPING
+    //
+
+    withName: SAMTOOLS_MPILEUP_PILEUPCALLER {
+        tag = { "${meta.reference}|${meta.strandedness}" }
+        ext.args = [
+            "-B",
+            "-q ${params.genotyping_pileupcaller_min_base_quality}",
+            "-Q ${params.genotyping_pileupcaller_min_map_quality}",
+        ].join(' ').trim()
+        ext.prefix = { "${meta.strandedness}_${meta.reference}" }
+        publishDir = [
+            enabled: false
+        ]
+    }
+
+    withName: SEQUENCETOOLS_PILEUPCALLER {
+        tag = { "${meta.reference}|${meta.strandedness}" }
+        ext.args = {[
+            "--${params.genotyping_pileupcaller_method}",
+            params.genotyping_pileupcaller_transitions_mode == "SkipTransitions" ? "--skipTransitions" : params.genotyping_pileupcaller_transitions_mode == "TransitionsMissing" ? "--transitionsMissing" : "",
+            "${meta.strandedness}" == 'single' ? "--singleStrandMode" : "" ,
+            "--sampleNames", meta.sample_id.join(","),
+            "-e pileupcaller.${meta.strandedness}.${meta.reference}"
+        ].join(' ').trim() }
+        ext.prefix = { "${meta.strandedness}_${meta.reference}" }
+        publishDir = [
+            enabled: false // Not published because the output goes through COLLECT_GENOTYPES
+        ]
+    }
+
+    withName: COLLECT_GENOTYPES {
+        tag = { "${meta.reference}" }
+        ext.prefix = { "pileupcaller_genotypes_${meta.reference}" }
+        publishDir = [
+            path: { "${params.outdir}/genotyping/" },
+            mode: params.publish_dir_mode,
+            enabled: true,
+            pattern: '*.{geno,snp,ind}'
+        ]
+    }
+
+    withName: EIGENSTRATDATABASETOOLS_EIGENSTRATSNPCOVERAGE {
+        tag = { "${meta.reference}" }
+        ext.args = { "-j ${prefix}.json" }
+        ext.prefix = { "pileupcaller_genotypes_${meta.reference}_coverage" }
+        publishDir = [
+            path: { "${params.outdir}/genotyping/" },
+            mode: params.publish_dir_mode,
+            enabled: true,
+            pattern: '*.{tsv}'
+        ]
+    }
+
+    withName: GATK_REALIGNERTARGETCREATOR {
+        tag = { "${meta.reference}|${meta.sample_id}" }
+        ext.args = [
+            params.genotyping_gatk_ug_defaultbasequalities > 0 ? "--defaultBaseQualities ${params.genotyping_gatk_ug_defaultbasequalities}" : "", // Empty string since GATK complains if its default of -1 is provided.
+        ].join(' ').trim()
+        ext.prefix = { "${meta.sample_id}_${meta.reference}_realigntarget" }
+        publishDir = [
+            enabled: false
+        ]
+    }
+
+    withName: GATK_INDELREALIGNER {
+        tag = { "${meta.reference}|${meta.sample_id}" }
+        ext.args = [
+            params.genotyping_gatk_ug_defaultbasequalities > 0 ? "--defaultBaseQualities ${params.genotyping_gatk_ug_defaultbasequalities}" : "", // Empty string since GATK complains if its default of -1 is provided.
+        ].join(' ').trim()
+        ext.prefix = { "${meta.sample_id}_${meta.reference}_realigned" }
+        publishDir = [
+            path: { "${params.outdir}/genotyping/IndelRealigner" },
+            mode: params.publish_dir_mode,
+            enabled: params.genotyping_gatk_ug_keeprealignbam,
+            pattern: '*.{bam,bai}'
+        ]
+    }
+
+    withName: GATK_UNIFIEDGENOTYPER {
+        tag = { "${meta.reference}|${meta.sample_id}" }
+        ext.args = {[
+            "--sample_ploidy ${meta2.ploidy}",
+            "-stand_call_conf ${params.genotyping_gatk_call_conf}",
+            "-dcov ${params.genotyping_gatk_ug_downsample}",
+            "--output_mode ${params.genotyping_gatk_ug_out_mode}",
+            "--genotype_likelihoods_model ${params.genotyping_gatk_ug_genotype_mode}",
+            params.genotyping_gatk_ug_defaultbasequalities > 0 ? "--defaultBaseQualities ${params.genotyping_gatk_ug_defaultbasequalities}" : "", // Empty string since GATK complains if its default of -1 is provided.
+        ].join(' ').trim() }
+        ext.prefix = { "${meta.sample_id}_${meta.reference}" }
+        publishDir = [
+            path: { "${params.outdir}/genotyping/" },
+            mode: params.publish_dir_mode,
+            enabled: true,
+            pattern: '*.vcf.gz'
+        ]
+    }
+
+    withName: BCFTOOLS_INDEX_UG {
+        tag = { "${meta.reference}|${meta.sample_id}" }
+        ext.args = "--tbi" //tbi indices for consistency with GATK HC
+        ext.prefix = { "${meta.sample_id}_${meta.reference}" }
+        publishDir = [
+            path: { "${params.outdir}/genotyping/" },
+            mode: params.publish_dir_mode,
+            enabled: true,
+            pattern: '*.vcf.gz.tbi'
+        ]
+    }
+
+    withName: GATK4_HAPLOTYPECALLER {
+        tag = { "${meta.reference}|${meta.sample_id}" }
+        ext.args = {[
+            // Option names have changed from underscore_separated to hyphen-separated in GATK4
+            "--sample-ploidy ${meta2.ploidy}",
+            "-stand-call-conf ${params.genotyping_gatk_call_conf}",
+            "--output-mode ${params.genotyping_gatk_hc_out_mode}",
+            "--emit-ref-confidence ${params.genotyping_gatk_hc_emitrefconf}",
+        ].join(' ').trim() }
+        ext.prefix = { "${meta.sample_id}_${meta.reference}" }
+        publishDir = [
+            path: { "${params.outdir}/genotyping/" },
+            mode: params.publish_dir_mode,
+            enabled: true,
+            pattern: '*.{vcf.gz,vcf.gz.tbi}'
+        ]
+    }
+
+    withName: FREEBAYES {
+        tag = { "${meta.reference}|${meta.sample_id}" }
+        ext.args = {[
+            "-p ${ref_meta.ploidy}",
+            "-C ${params.genotyping_freebayes_min_alternate_count}",
+            params.genotyping_freebayes_skip_coverage != 0 ? "-g ${params.genotyping_freebayes_skip_coverage}" : "",
+        ].join(' ').trim() }
+        ext.prefix = { "${meta.sample_id}_${meta.reference}" }
+        publishDir = [
+            path: { "${params.outdir}/genotyping/" },
+            mode: params.publish_dir_mode,
+            enabled: true,
+            pattern: '*.vcf.gz'
+        ]
+    }
+
+    withName: BCFTOOLS_INDEX_FREEBAYES {
+        tag = { "${meta.reference}|${meta.sample_id}" }
+        ext.args = "--tbi" //tbi indices for consistency with GATK HC
+        ext.prefix = { "${meta.sample_id}_${meta.reference}" }
+        publishDir = [
+            path: { "${params.outdir}/genotyping/" },
+            mode: params.publish_dir_mode,
+            enabled: true,
+            pattern: '*.vcf.gz.tbi'
+        ]
+    }
+
+    withName: BCFTOOLS_STATS_GENOTYPING {
+        tag = { "${meta.reference}|${meta.sample_id}" }
+        ext.prefix = { "${meta.sample_id}_${meta.reference}" }
+        publishDir = [
+            path: { "${params.outdir}/genotyping/" },
+            mode: params.publish_dir_mode,
+            enabled: true,
+            pattern: '*.txt'
+        ]
+    }
 }

conf/test_humanbam.config

@@ -38,8 +38,8 @@ params {
     // //Sex Determination
     // sexdeterrmine_bedfile = 'https://raw.githubusercontent.com/nf-core/test-datasets/eager/reference/Human/1240K.pos.list_hs37d5.0based.bed.gz'
     // // Genotyping
-    // pileupcaller_bedfile = 'https://raw.githubusercontent.com/nf-core/test-datasets/eager/reference/Human/1240K.pos.list_hs37d5.0based.bed.gz'
-    // pileupcaller_snpfile = 'https://raw.githubusercontent.com/nf-core/test-datasets/eager/reference/Human/1240K_covered_in_JK2067_downsampled_s0.1.numeric_chromosomes.snp'
+    genotyping_pileupcaller_bedfile = 'https://raw.githubusercontent.com/nf-core/test-datasets/eager/reference/Human/1240K.pos.list_hs37d5.0based.bed.gz'
+    genotyping_pileupcaller_snpfile = 'https://raw.githubusercontent.com/nf-core/test-datasets/eager/reference/Human/1240K_covered_in_JK2067_downsampled_s0.1.numeric_chromosomes.snp'
 
 
     // BAM filtering