Какие брать концентрации нуклеотидов?

import numpy as np
import matplotlib.pyplot as plt

ddNTP = 0.001 # Вероятность того, что встроится ddNTP
num_molecules = 100000 # Число молекул в реакционной смеси
length = 1000 # Длинна таргетного участка

lens = np.random.geometric(p = ddNTP, size=num_molecules)
lens = lens[lens < length]
lens = np.append(lens, [length] * (num_molecules - len(lens)))

fig, axes = plt.subplots(ncols=2, dpi=150, figsize=(13, 3))
axes[0].hist(lens, edgecolor="k", bins=100)
axes[0].set_xlabel("Molecule length")
axes[1].hist(lens[lens != length], edgecolor="k", bins=100)
axes[1].set_xlabel("Molecule length")

print(f"Коротких отрезкв в {sum(lens < 10) / sum(np.logical_and(lens > length - 10, lens != length)):.2f} раз больше, чем длинных")

Коротких отрезкв в 2.57 раз больше, чем длинных

Определение качества определения позиций

Для начала мы загрузим подгрузим пакеты. По умолчанию их нет, но их можно установить при помощи команды:

pip3 install biopython seaborn matplotlib

from Bio import SeqIO
import seaborn as sns
import pandas as pd

Далее мы посмотрим список файлов (чтобы не прописывать их руками) и подгрузим их все с список records.

record_names = !ls ab1_files
records = []
for record_name in record_names:
    records.append(SeqIO.read("ab1_files/" + record_name, "abi"))

На примере одной последовательности нарисуем распределение качества прочтения.

plt.rcParams["figure.figsize"] = (20, 5)
ax = plt.plot(records[0].letter_annotations["phred_quality"])

Посмотрим на то, что творится с качеством у какой-нибудь другой последовательности.

plt.plot(records[45].letter_annotations["phred_quality"])

[<matplotlib.lines.Line2D at 0x7f88b7313d68>]

У неё тут всё гораздо ровнее в центре. Давайте попытаемся как-то визуализировать распределение качества в двух координатных осях: (а) в зависимости от того, какой нуклеотид у нас по порядку и (б) в зависимости от того, в какой части хроматограммы он оказался в процентах (т. е. отнормируем длины).

data_quality = {
    "Position" : [],
    "Quality" : []
}

for record in records:
    data_quality["Position"] += list(range(1, len(record) + 1))
    data_quality["Quality"] += record.letter_annotations["phred_quality"]
    
data_quality = pd.DataFrame(data_quality)
sns.relplot(data=data_quality, x="Position", y="Quality", kind="line", ci="sd",
            height=5, aspect=4)

<seaborn.axisgrid.FacetGrid at 0x7f88b7340860>

data_quality = {
    "Position (%)" : [],
    "Quality" : []
}

for record in records:
    data_quality["Position (%)"] += list(map(round, list(np.array(list(range(1, len(record) + 1))) / len(record) * 100)))
    data_quality["Quality"] += record.letter_annotations["phred_quality"]
    
data_quality = pd.DataFrame(data_quality)
sns.relplot(data=data_quality, x="Position (%)", y="Quality", kind="line", ci="sd",
            height=5, aspect=4)

<seaborn.axisgrid.FacetGrid at 0x7f88b71b0b38>

Сделаем дополнительный контроль — разобьём на диапозоны длина последовательности // 200.

data_quality = {
    "Position" : [],
    "Quality" : [],
    "Length range" : []
}

step = 150
for record in records:
    data_quality["Position"] += list(range(1, len(record) + 1))
    data_quality["Quality"] += record.letter_annotations["phred_quality"]
    data_quality["Length range"] += len(record) * [str(len(record) // step * step) + " - " + str((len(record) + step) // step * step)]
    
data_quality = pd.DataFrame(data_quality)
data_quality["Length range"].apply(str)
sns.relplot(data=data_quality, x="Position", y="Quality", kind="line", ci="sd",
            height=5, aspect=4, hue="Length range")

<seaborn.axisgrid.FacetGrid at 0x7f88b719f6d8>

sns.relplot(data=data_quality, x="Position", y="Quality", kind="line", ci=None,
            height=5, aspect=4, hue="Length range")

<seaborn.axisgrid.FacetGrid at 0x7f8868881e80>

Вместо качества используем вероятность ошибки

def phred_score(p):
    return (-10 * np.log10(p))

def error_prob(phred):
    return 10 ** (float(phred) / -10)

data_quality = {
    "Position" : [],
    "Error probability" : []
}

for record in records:
    data_quality["Position"] += list(range(1, len(record) + 1))
    data_quality["Error probability"] += record.letter_annotations["phred_quality"]

data_quality = pd.DataFrame(data_quality)
data_quality["Error probability"] = data_quality["Error probability"].apply(error_prob)
ax = sns.relplot(data=data_quality, x="Position", y="Error probability", kind="line", ci="sd",
                 height=5, aspect=4)

y_threshold = 0.1
plt.plot([0, max(data_quality.Position)], [y_threshold, y_threshold], color="red", linestyle="--", linewidth=1)
ax.set(xlim=(0, max(data_quality.Position)))

<seaborn.axisgrid.FacetGrid at 0x7f88688204a8>