Building a custom genome

sxsw-2017.py

#!/opt/venvs/deskgen/bin/python

"""
Script to build genome, call variants, patch, and score
"""
import logging
import multiprocessing
import os
import sys

import genome_fs
from omics_tools import bowtie2, varianttool
LOG = logging.getLogger(__name__)
EXPDIR = os.path.dirname(os.path.abspath(__file__))
logging.basicConfig(filename=os.path.join(EXPDIR, 'error.log'))
USER_EMAIL = 'rileyd@deskgen.com'
GENOME_URL = '/opt/biodata/human/GRCh38.81'
STRAIN = 'GM24385'
SAMPLE = 'SRR5120985'
# load our genome object
GENOME = genome_fs.discover(GENOME_URL)
READ_LEN = '125'
THREADS = multiprocessing.cpu_count()
ALIGNMENT_PARAMS = ' '.join([
    '--phred33',
    '--very-sensitive',
    '-p', str(THREADS)
])

def align_reads(bam_aligned):

    bt2params = ALIGNMENT_PARAMS
    b2index = os.path.join(GENOME_URL, 'formated')
    fq_r1pairedzip = os.path.join(EXPDIR, SAMPLE + '_1.fastq.gz')
    fq_r2pairedzip = os.path.join(EXPDIR, SAMPLE + '_2.fastq.gz')
    bam_aligned = os.path.join(EXPDIR, SAMPLE + '.bam')
    bowtie2.bt2pairedalign(
            bt2params,
            b2index,
            fq_r1pairedzip,
            fq_r2pairedzip,
            bam_aligned
    )


if __name__ == '__main__':
    # Technically best to simply search for .bt2 files
    bam_aligned = os.path.join(EXPDIR, SAMPLE + '.bam')
    if not os.path.exists(bam_aligned):
        align_reads(bam_aligned)
    fasta = os.path.join(GENOME_URL, 'formatted/GRCh38.81.fa')
    chrom_list = GENOME.items
    bam_list = [bam_aligned]
    out_vcf = os.path.join(EXPDIR, SAMPLE+'.vcf')
    keeptemp = True
    varianttool.freebayes_parallel(
        THREADS,
        chrom_list,
        fasta,
        bam_list,
        params,
        out_vcf,
        keeptemp
    )
    sys.exit(0)