# -*- coding: utf-8 -*-
"""Quality_reads_fastq
# Checking quality of the reads - Working with FASTQ files
### Getting the data

!rm -f SRR003265.filt.fastq.gz 2>/dev/null
!wget -nd ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/NA18489/sequence_read/SRR003265.filt.fastq.gz

"""### Reading FASTQ files"""
from collections import defaultdict
import gzip
# import seaborn as sns
import matplotlib.pyplot as plt
from Bio import SeqIO

recs = SeqIO.parse('./SRR494102.fastq', 'fastq')
print(rec)
print(rec.id, rec.description, rec.seq)
print(rec.letter_annotations) # quality scores of our reads, per letter

"""### Distribution of nucleotide reads"""
def get_distribution_of_nuc_reads():
    count = defaultdict(int)
    for rec in recs:
        for letter in rec.seq:
            count[letter] += 1
    tot = sum(count.values())
    for letter, count in count.items():
        print('%s: %.2f %d' % (letter, 100. * count / tot, count))

"""N - uncalled base"""

recs = SeqIO.parse(gzip.open('SRR003265.filt.fastq.gz', 'rt', encoding='UTF-8'), 'fastq')
n_cnt = defaultdict(int)
for rec in recs:
    for i, letter in enumerate(rec.seq):
        pos = i + 1
        if letter == 'N':
            n_cnt[pos] += 1
seq_len = max(n_cnt.keys())
positions = range(1, seq_len + 1)
fig, ax = plt.subplots(figsize=(16,9))
ax.plot(positions, [n_cnt[x] for x in positions])
ax.set_xlim(1, seq_len)
pass

"""### Quality of reads"""

recs = SeqIO.parse(gzip.open('SRR003265.filt.fastq.gz', 'rt', encoding='utf-8'), 'fastq')
cnt_qual = defaultdict(int)
for rec in recs:
    for i, qual in enumerate(rec.letter_annotations['phred_quality']):
        if i < 25:
            continue
        cnt_qual[qual] += 1
tot = sum(cnt_qual.values())
for qual, cnt in cnt_qual.items():
    print('%d: %.2f %d' % (qual, 100. * cnt / tot, cnt))

recs = SeqIO.parse(gzip.open('SRR003265.filt.fastq.gz', 'rt', encoding='utf-8'), 'fastq')
qual_pos = defaultdict(list)
for rec in recs:
    for i, qual in enumerate(rec.letter_annotations['phred_quality']):
        if i < 25 or qual == 40:
            continue
        pos = i + 1
        qual_pos[pos].append(qual)

vps = []

poses = list(qual_pos.keys())
poses.sort()

for pos in poses:
    vps.append(qual_pos[pos])

fig, ax = plt.subplots(figsize=(16,9))
plt.boxplot(data=vps, ax=ax)
ax.set_xticklabels([str(x) for x in range(26, max(qual_pos.keys()) + 1)])

