Анализ результатов секвенирования по Сэнгеру

Гуков, 201

Задание 1

Результаты работы:

• Консенсусная последовательность
• Выравнивание прямой и обратной последовательнойстей до корректировки ошибок
• Референс
• Выравнивание прямой и обратной последовательнойстей после корректировки ошибок
• 44_F.ab1
• 44_R.ab1

Для начала оценим хроматограммы для прямой и обратной последовательностей метрикой оценки качества сигнала - Q, которая может быть расчитана по формуле: $$Q = {-10}\cdot\log_{10}{p}$$ Где Q - качество сигнала для данной вероятности p.

Удобно построить графики зависимостей качества сигнала от позиции конкретного нуклеотида. Данные для построения графиков уже находятся в бинарных файлах расширения ab1. В коде ниже файлы 44_F.ab1 и 44_R.ab1 хранятся в списке records. Далее создаётся pandas.DataFrame с именем data_quality, в котором первый столбец - имя последовательности из 44_*.ab1, второй столбец - рассматриваемая позиция, а третий столбец - мера качества для данной позиции. Данные в data_quality находятся в формате long-form data для удобной работы с модулем визуализации данных seaborn.

По графику можно сказать, что хроматограммы достаточно качественные (более точная оценка ниже), видны начальные и конечные нечитаемые участки.

from Bio import SeqIO
import seaborn as sns
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

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

data_quality = {
    'name': [],
    'position': [],
    'quality': []
}

for k in range(len(records)):
    for i, n in enumerate(records[k].letter_annotations['phred_quality'], start = 1):
        data_quality['position'].append(i)
        data_quality['quality'].append(n)
        data_quality['name'].append(records[k].name)


data_quality = pd.DataFrame(data_quality)
f, ax = plt.subplots(2, 1, figsize=(14, 8), dpi = 150)
sns.set_context('paper')
sns.lineplot(ax = ax[0], data = data_quality, x = 'position', y = 'quality', linewidth = .7, hue = 'name', ci = 'sd', palette = 'BuPu')
ax[0].set_xlabel('')
sns.lineplot(ax = ax[1], data = data_quality, x = 'position', y = 'quality', linewidth = .7, ci = 'sd', color = '#815DA6')
plt.show()

Можно более точно описать "насколько графики качественные", используя возможности pandas, к примеру метод describe():

round(data_quality.loc[data_quality['name'] == records[0].name]['quality'].describe(), 2) # Для прямой цепи

count    719.00
mean      32.50
std       12.07
min        1.00
25%       26.00
50%       33.00
75%       39.00
max       58.00
Name: quality, dtype: float64

round(data_quality.loc[data_quality['name'] == records[1].name]['quality'].describe(), 2) # Для обратной цепи

count    718.00
mean      29.04
std       10.75
min        1.00
25%       24.00
50%       28.00
75%       35.00
max       58.00
Name: quality, dtype: float64

Среднее значение Q будет обозначаться как Q_mean.

Q_mean прямой цепи равняется 32.5, а Q_mean обратной цепи - 29.04. Данные значения являются показателями того, что хроматограммы действительно хорошие.

Таже с помощью библиотеки seaborn можно построить box plot ("ящик с усами") и визуально оценить как себя ведёт мера качества взависимости от хроматограммы. В принципе, график только подтверждает качественность хроматограмм, потому что значения равные Q = 1 (очень высокая вероятность ошибки в данной позиции) даже не входят в сам box plot, а являются выбросами.

plt.figure(figsize=(6, 6), dpi = 100)
sns.boxplot(data = data_quality, x = 'name', y = 'quality', palette = 'BuPu')
plt.show()

Попробуем сделать интересенее и посмотрим график зависимости вероятности ошибки от позиции потенциальной ошибки. Для расчета вероятности ошибки воспользуемся формулой качества сигнала: $$Q = {-10}\cdot\log_{10}{p}$$ Так как Q известно изначально, можно выразить p: $$p = {10}^{\frac{Q}{-10}}$$

В дальнейшей работе будем придерживаться идеи, что если вероятность ошибки больше 20%, то лучше проверить это по хроматограмме.

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

data_quality = {
    'name': [],
    'position': [],
    'error probability': []
}

for k in range(len(records)):
    for i, n in enumerate(records[k].letter_annotations['phred_quality'], start = 1):
        data_quality['position'].append(i)
        data_quality['error probability'].append(n)
        data_quality['name'].append(records[k].name)
        
plt.figure(figsize=(14, 8), dpi = 100)
sns.set_context('paper')
data_quality = pd.DataFrame(data_quality)
data_quality['error probability'] = round(data_quality['error probability'].apply(error_prob), 2)
ax = sns.relplot(data = data_quality, x = 'position', y = 'error probability', hue = 'name', kind = 'line', ci = 'sd', height=5, aspect=4, palette = 'BuPu')

y_threshold = .2
plt.plot([0, max(data_quality.position)], [y_threshold, y_threshold], color = 'red', linestyle = '--', linewidth = 1)
ax.set(xlim = (0, max(data_quality.position)))
plt.show()

<Figure size 1400x800 with 0 Axes>

В основном позиции с высокой вероятностью ошибки находятся либо в начале, либо в конце. Это, так называемые, начальный и конечный нечитаемые учатстики.

Ниже представлены таблицы, в которую помещены позиции, вероятность ошибки в которых больше, чем 20%. Первая таблица - это таблица с матричной цепью (forward). Вторая таблица описывает не хорошие позиции обратной (reverse) цепи. Данные таблицы были сделаны для того, чтобы был ориентир на те позиции, которые требуют к себе повышенное внимание.

По первой таблице, можно сразу же предположить, что 5'-1:20-3' - это начальный нечитаемый участок, а 5'-675:719-3' - это конечный нечитаемый участок, однако утверждения необходимо проверить непосредственно на хроматограмме.

data_quality.loc[(data_quality['error probability'] > .2) & (data_quality['name'] == '44_F')]

	name	position	error probability
1	44_F	2	0.79
2	44_F	3	0.79
3	44_F	4	0.79
4	44_F	5	0.32
5	44_F	6	0.79
7	44_F	8	0.79
8	44_F	9	0.79
9	44_F	10	0.32
10	44_F	11	0.32
11	44_F	12	0.32
12	44_F	13	0.32
13	44_F	14	0.32
14	44_F	15	0.79
16	44_F	17	0.32
17	44_F	18	0.79
18	44_F	19	0.79
19	44_F	20	0.32
131	44_F	132	0.79
236	44_F	237	0.79
363	44_F	364	0.79
385	44_F	386	0.50
470	44_F	471	0.79
667	44_F	668	0.79
674	44_F	675	0.79
677	44_F	678	0.32
680	44_F	681	0.79
686	44_F	687	0.63
694	44_F	695	0.79
704	44_F	705	0.79
707	44_F	708	0.79
708	44_F	709	0.79
711	44_F	712	0.79
712	44_F	713	0.79
714	44_F	715	0.79
715	44_F	716	0.79
716	44_F	717	0.32

Точно также как и о матричной цепи, о комплиментарной можно сказать, что начальным нечитаемым участком является последовательность 5'-1:21-3', а конечным - 5'-717-718-3'. Напоминаю, что эти значения - всего лишь ориентир. Необходимо проверять значения по хроматограммам.

data_quality.loc[(data_quality['error probability'] > .2) & (data_quality['name'] == '44_R')]

	name	position	error probability
721	44_R	3	0.50
722	44_R	4	0.50
723	44_R	5	0.50
724	44_R	6	0.79
725	44_R	7	0.32
726	44_R	8	0.32
727	44_R	9	0.32
730	44_R	12	0.79
731	44_R	13	0.25
733	44_R	15	0.32
734	44_R	16	0.32
735	44_R	17	0.32
746	44_R	28	0.32
747	44_R	29	0.32
748	44_R	30	0.32
749	44_R	31	0.32
760	44_R	42	0.79
944	44_R	226	0.79
986	44_R	268	0.25
1059	44_R	341	0.40
1112	44_R	394	0.79
1249	44_R	531	0.79
1288	44_R	570	0.79
1321	44_R	603	0.63
1395	44_R	677	0.79
1396	44_R	678	0.79
1399	44_R	681	0.32
1435	44_R	717	0.79
1436	44_R	718	0.32

Непосредственно посмотрев на хроматограмы были приянты решения? которые совпадают с решением программы Chromas:

• У прямой цепи начальный нечитаемый участок 5'_1:23_3', а конечный - 5'_675:719_3'
• У обратной цепи начальный нечитаемый учатсок 5'_1:31_3', а конечный - 5'_714:718_3'

Переменные, которые объявлены ниже не содержат начальных и конечных нечитаемых участков.

from Bio.Seq import Seq
seqF = Seq('TTTTATATTTGGCGCCTGAGCAGGTACAGTAGGGACTGCCATGAGAAAAATTATACGAGTTGAACTTTCTCAGCCAGGCTCTTTAATACAAGATGATCAAGTATATAANGTTATGGTAACGGCCCACGCCTTCGTCATGATATTTTTTATGGTAATGCCCATAATGATAGGGGGGTTTGGCAAATGACTTGTCCCACTAATGTTAGGAGCGCCNGATATGGCTTTCCCCCGAATGAAAAAAATGAGATTTTGGCTACTACCCCCAGCTTTTATACTTCTTCTAGCTTCAGCTGCAAACGAAGGAGGAGTAGGCACTGGATGAACTATTTATCCCCCTTTGNCAGGCCCTACCGCACATGCAGNAGGCTGCGTAGACCTCGCAATTTTTTCTCTCCACCTAGCAGGTGCGTCTTCAATTATGGCCTCAATAAAATTTATTACAACTATNATAAATATGCGTAGGCCCGGCATGACCATGGATCGACTTCCACTTTTTGCTTGATCTATTTTCTTAACAACTATATTACTACTCCTTTCTCTGCCTGTTTTAGCAGGAGCTATTACAATGCTATTAACTGATCGTAACATAAAAACAACGTTTTTTGATCCTACAGGAGGAGGAGACCCAATACTTTTCCAACATTNATTTTG')
seqR = Seq('NTATTGGGTCNCCTCCTCCTGTAGGATCAAAAAACGTTGTTTTTATGTTANGATCAGTTAATAGCATTGTAATAGCTCCTGCTAAAACAGGCAGAGAAAGGAGTAGTAATATAGTTGTTAAGAAAATAGATCAAGCAAAAAGTGGAAGTCGATCCATGGTCATGCCGGGCCTACGCATATTTATAATAGTTGTANTAAATTTTATTGAGGCCATAATTGAAGACGCACCTGCTAGGNGGAGAGAAAAAATTGCGAGGTCTACGCAGCCTCCTGCATGTGCGGTAGGGCCTGACAAAGGGGGATAAATAGNTCATCCAGTGCCTACTCCTCCTTCGTTTGCAGCTGAAGCTAGAAGAAGTATANAAGCTGGGGGTAGTAGCCAAAATCTCATTTTTTTCATTCGGGGGAAAGCCATATCAGGCGCTCCTAACATTAGTGGGACAAGTCATTTGCCAAACCCCCCTATCATTATGGGCATTACCATAAAAAATATCANGACNAAGGCGTGGGCCGTTACCATAACTTTATATACTTGATCNTCTTGTATTAAAGAGCCTGGCTGAGAAAGTTCNACTCGTATAATTTTTCTCATGGCAGTCCCTACTGTACCTGCTCAGGCGCCAAATATAAAATATAGTGTTCCAANNTCNTTGTGATTCGTCGACAACTGGCCGTCGTTTTA')

seqRrevcom = seqR.reverse_complement()
print(seqRrevcom)

TAAAACGACGGCCAGTTGTCGACGAATCACAANGANNTTGGAACACTATATTTTATATTTGGCGCCTGAGCAGGTACAGTAGGGACTGCCATGAGAAAAATTATACGAGTNGAACTTTCTCAGCCAGGCTCTTTAATACAAGANGATCAAGTATATAAAGTTATGGTAACGGCCCACGCCTTNGTCNTGATATTTTTTATGGTAATGCCCATAATGATAGGGGGGTTTGGCAAATGACTTGTCCCACTAATGTTAGGAGCGCCTGATATGGCTTTCCCCCGAATGAAAAAAATGAGATTTTGGCTACTACCCCCAGCTTNTATACTTCTTCTAGCTTCAGCTGCAAACGAAGGAGGAGTAGGCACTGGATGANCTATTTATCCCCCTTTGTCAGGCCCTACCGCACATGCAGGAGGCTGCGTAGACCTCGCAATTTTTTCTCTCCNCCTAGCAGGTGCGTCTTCAATTATGGCCTCAATAAAATTTANTACAACTATTATAAATATGCGTAGGCCCGGCATGACCATGGATCGACTTCCACTTTTTGCTTGATCTATTTTCTTAACAACTATATTACTACTCCTTTCTCTGCCTGTTTTAGCAGGAGCTATTACAATGCTATTAACTGATCNTAACATAAAAACAACGTTTTTTGATCCTACAGGAGGAGGNGACCCAATAN

from Bio import pairwise2
alignments = pairwise2.align.globalxx(seqF, seqRrevcom)
for alig in alignments:
    print(pairwise2.format_alignment(*alig))

T--------------TT-T--A----T-A--------TT-------T-------------GGCGCCTGAGCAGGTACAGTAGGGACTGCCATGAGAAAAATTATACGAGTT-GAACTTTCTCAGCCAGGCTCTTTAATACAAGAT-GATCAAGTATATAAN-GTTATGGTAACGGCCCACGCCTTC-GTCA-TGATATTTTTTATGGTAATGCCCATAATGATAGGGGGGTTTGGCAAATGACTTGTCCCACTAATGTTAGGAGCGCCN-GATATGGCTTTCCCCCGAATGAAAAAAATGAGATTTTGGCTACTACCCCCAGCTTT-TATACTTCTTCTAGCTTCAGCTGCAAACGAAGGAGGAGTAGGCACTGGATGAA-CTATTTATCCCCCTTTGN-CAGGCCCTACCGCACATGCAGN-AGGCTGCGTAGACCTCGCAATTTTTTCTCTCCA-CCTAGCAGGTGCGTCTTCAATTATGGCCTCAATAAAATTTAT-TACAACTATN-ATAAATATGCGTAGGCCCGGCATGACCATGGATCGACTTCCACTTTTTGCTTGATCTATTTTCTTAACAACTATATTACTACTCCTTTCTCTGCCTGTTTTAGCAGGAGCTATTACAATGCTATTAACTGATCG-TAACATAAAAACAACGTTTTTTGATCCTACAGGAGGAGGA-GACCCAATACTTTTCCAACATTNATTTTG
|              || |  |    | |        ||       |             ||||||||||||||||||||||||||||||||||||||||||||||||| | ||||||||||||||||||||||||||||||||  ||||||||||||||  |||||||||||||||||||||||  |||  ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||  ||||||||||||||||||||||||||||||||||||||||||||||||||||| || ||||||||||||||||||||||||||||||||||||||||||||||||||| | |||||||||||||||||  |||||||||||||||||||||  ||||||||||||||||||||||||||||||||  |||||||||||||||||||||||||||||||||||||||||  |||||||||  |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||  |||||||||||||||||||||||||||||||||||||||  ||||| | |    |     |  |      
TAAAACGACGGCCAGTTGTCGACGAATCACAANGANNTTGGAACACTATATTTTATATTTGGCGCCTGAGCAGGTACAGTAGGGACTGCCATGAGAAAAATTATACGAG-TNGAACTTTCTCAGCCAGGCTCTTTAATACAAGA-NGATCAAGTATATAA-AGTTATGGTAACGGCCCACGCCTT-NGTC-NTGATATTTTTTATGGTAATGCCCATAATGATAGGGGGGTTTGGCAAATGACTTGTCCCACTAATGTTAGGAGCGCC-TGATATGGCTTTCCCCCGAATGAAAAAAATGAGATTTTGGCTACTACCCCCAGC-TTNTATACTTCTTCTAGCTTCAGCTGCAAACGAAGGAGGAGTAGGCACTGGATG-ANCTATTTATCCCCCTTTG-TCAGGCCCTACCGCACATGCAG-GAGGCTGCGTAGACCTCGCAATTTTTTCTCTCC-NCCTAGCAGGTGCGTCTTCAATTATGGCCTCAATAAAATTTA-NTACAACTAT-TATAAATATGCGTAGGCCCGGCATGACCATGGATCGACTTCCACTTTTTGCTTGATCTATTTTCTTAACAACTATATTACTACTCCTTTCTCTGCCTGTTTTAGCAGGAGCTATTACAATGCTATTAACTGATC-NTAACATAAAAACAACGTTTTTTGATCCTACAGGAGGAGG-NGACCC-A-A----T-----A--N------
  Score=617

Я не знаю как решить проблему отображения, поэтому взял последовательности из переменных seqF и seqRrevcov и пошёл делать выравнивание в Jalview. В Jalview было пройзведено выравнивание программой Muscle. Далее, с использованием Chromas, непосредственно в Jalview производился анализ последовательностей для установления консенсусной цепи.

Примеры ошибок, с которыми произошла встреча во время просмотра выравнивания и хроматограмм:

Сверху 44_F, снизу 44_R_rev_com

• Наложение пиков T и C на 144 позиции 44_R_rev_com. Было принято решение поставить T, так как на 44_F красный пик на данной позиции.

• Наложение пиков T и G на 214 позиции 44_F. Было принято решение поставить T, так как на 44_R_rev_com красный пик на данной позиции.

• Выравнивание показало, что на 507 позиции 44_F находися пропуск. Однако, можно увидеть, что там точно находится T.

• На 33 позиции 44_rev_com находится очень интересная ситуация: A и G пересекаются, однако, моё решение такого, что там нет нуклеотида вообще.

Задание 2

Скрин был сделан из хроматограммы NN_G10.ab1. Вероятно, там несколько молекул ДНК, поэтому много пиков накладываются друг на друга в одной и той же позиции.