Bug in call-ataqc module

[nicolerg]

The call-ataqc module is failing with the following python variable assignment error:

Traceback (most recent call last):
  File "/users/nicolerg/anaconda2/envs/encode-atac-seq-pipeline/bin/encode_ataqc.py", line 355, in <module>
    ataqc()
  File "/users/nicolerg/anaconda2/envs/encode-atac-seq-pipeline/bin/encode_ataqc.py", line 122, in ataqc
    chr_m_reads, fraction_chr_m = get_chr_m(COORDSORT_BAM)
  File "/users/nicolerg/anaconda2/envs/encode-atac-seq-pipeline/bin/run_ataqc.py", line 215, in get_chr_m
    fract_chr_m = float(chr_m_reads) / tot_reads
UnboundLocalError: local variable 'chr_m_reads' referenced before assignment

OS/Platform and dependencies

• OS or Platform: Ubuntu 16.04.4 (durga)
• Cromwell/dxWDL version: cromwell 34-unknown-SNAP
• Conda version: conda 4.5.11

I'm happy to provide error logs if needed, but this seems pretty cut and dry.

[leepc12]

It looks like your BAM file is empty. Please provide error logs (.tar.gz).

[nicolerg]

I'm not sure if either of these are the BAM file you're referring to:
call-ataqc/shard-0/inputs/789850838/${SAMPLE}.trim.merged.bam is 1.5G
call-ataqc/shard-0/inputs/-623782124/${SAMPLE}.trim.merged.nodup.bam is 667M.

debug_64.tar.gz

[leepc12]

This will be fixed in v1.1.4. But if you want to fix it yourself then modify function get_chr_m() in line 203-221 of /users/nicolerg/anaconda2/envs/encode-atac-seq-pipeline/bin/encode_ataqc.py like the following.

def get_chr_m(sorted_bam_file):
    '''
    Get fraction of reads that are mitochondrial (chr M).
    '''
    logging.info('Getting mitochondrial chromosome fraction...')
    chrom_list = pysam.idxstats(sorted_bam_file, split_lines=True)
    tot_reads = 0
    chr_m_reads = 0
    for chrom in chrom_list:
        chrom_stats = chrom.split('\t')
        if chrom_stats[0] == 'chrM':
            chr_m_reads = int(chrom_stats[2])
        tot_reads += int(chrom_stats[2])
    if tot_reads==0:
        fract_chr_m = 0
    else:
        fract_chr_m = float(chr_m_reads) / tot_reads

    return chr_m_reads, fract_chr_m

[nicolerg]

That function is not in encode_ataqc.py, but I modified it in /users/nicolerg/anaconda2/envs/encode-atac-seq-pipeline/bin/run_ataqc.py.