Практикум 15 - сборка de novo
Код доступа проекта по секвенированию Buchnera aphidicola str. Tuc7: SRR4240380
Скачиваю файл с чтениями:
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR424/001/SRR4240381/SRR4240381.fastq.gz
Смотрю на качество чтений:
fastqc SRR4240381.fastq.gz
Всего чтений до триммирования: 13 710 994
Качество чтений не очень, но, может быть, после триммирования станет лучше.
Меня немного смутил график, отражающий сожержание адаптеров в чтениях. Судя по всему, мы ожидаем их увидеть, но не видим. Может быть, FastQC просто не видит те адаптеры, которые есть в данных чтениях.
- Подготовка чтений программой trimmomatic
- Запуск velveth
- data – название папки, в которую будет выдача
- 31 – длина k-меров
- -short – чтения короткие и не парные
- -fastq.gz – тип файла вывода
- output2.fastq.gz – файл с чтениями
- Запуск velvetg
- lgth – длина контига, выраженная в количестве узлов, которые объединены в контиг
- in и out – количество рёбер, присоединённых к 5'- и 3'-концам контига соответственно
- short1_cov и short1_Ocov – покрытие, посчитанное двумя разными способами; мне показалось, что в данном случае полученное число в обоих случаях получилось одном и тем же
- Анализ
- ID 6
- ID 9
- ID 49
Так как чтения непаноконцевые, используется TrimmomaticSE. ILLUMINACLIP убирает адаптеры из чтений, выравнивнивая их с
адаптерами.
Синтаксис таков: ILLUMINACLIP:<файл с последовательностями адаптеров>:<максимальное допустимое количество мисматчей>:<какой-то
параметр для PE, видимо>:<точность соответствия последовательности адаптера и последовательности чтения, что бы это ни значило>
Конкатенирую файлы с адаптерами в один:
cat *.fa > ../cringebutfree/pr15/adapters.fasta
Триммирую адаптеры:
TrimmomaticSE -phred33 SRR4240381.fastq.gz output1.fastq.gz
ILLUMINACLIP:adapters.fasta:2:7:7
Посмотрим, что из этого получилось. Выдача FastQC.
Количество "выживших" последовательностей: 13 704 864, то есть 0,04% последовательностей отсеялись полностью. Некоторые последовательности порезались в той или иной степени, и самые короткие уйдут на следующем шаге.
Теперь нужно убрать с правого конца нуклеотиды с качеством ниже 20 и оставить те чтения, у которых длина больше 32.
TrimmomaticSE -phred33 output1.fastq.gz output2.fastq.gz TRAILING:20 MINLEN:32
"Выжило" 11 219 816 последовательностей, это 81.87%.
Можно заметить, что качество чтений всё ещё довольно грустное, но на самом конце оно стало немного лучше. При этом нуклеотидный состав чтений на конце стал сильно неравномерным. По-хорошему, надо сделать триммирование со скользящим окном, но это позже.
Прежде чем запускать velvetg, нужно подготовить данные:
velveth data 31 -short -fastq.gz output2.fastq.gz
где
velvetg data
N50: 5987
Структура файла stats.txt следующая:
ID lgth out in long_cov short1_cov short1_Ocov short2_cov short2_Ocov long_nb short1_nb short2_nb
1 11333 0 0 0.000000 30.694873 30.694873 0.000000 0.000000 0 0 0
2 3903 0 0 0.000000 26.503203 26.503203 0.000000 0.000000 0 0 0
3 9973 0 0 0.000000 34.071593 34.071593 0.000000 0.000000 0 0 0
Bash сложный и непонятный, поэтому я немного схитрила:
contigs = {}
with open(path) as file:
file.readline() #Убираю заголовок
#Тут длина будет уже в нуклеотидах
for line in file.readlines():
line = line.split()
contigs[line[0]] = [int(line[1]) + 30, float(line[5])]
#Самые длинные контиги
longest = sorted(contigs, key=lambda x: contigs[x][0], reverse=True)[:3]
cols = ['ID', 'length', 'coverage']
for elem in cols:
print(elem.ljust(10), end='')
print()
for contig in longest:
print(contig.ljust(10), end='')
for elem in contigs[contig]:
print(str(elem).ljust(10), end='')
print()
Выдача выглядит следующим образом:
ID length coverage 6 69701 33.069843 9 27504 38.781939 49 26092 35.108817
Это ID, длина (в нуклеотидах) и покрытие трёх самых длинных контигов.
Теперь найду среднее покрытие и "выбивающиеся" по покрытию контиги:
num = 0
key in contigs.keys():
num += contigs[key][1]
num /= len(contigs)
print(f'Среднее значение: {num:.03f}')
print()
for elem in cols:
print(elem.ljust(10), end='')
print()
for key in sorted(contigs, key=lambda x: contigs[x][1]):
if contigs[key][1] < num / 5 or contigs[key][1] > num * 5:
print(key.ljust(10), end='')
for elem in contigs[key]:
print(str(elem).ljust(10), end='')
print()
Выдача выглядит так:
Среднее значение: 31.557
ID length coverage
2840 34 1.0
2933 33 1.0
2935 33 1.666667
2855 31 2.0
...
178 61 282.612903
105 338 283.016234
152 56 363.115385
287 61 542.258065
507 31 53018.0
Можно обратить внимание, что контиг с наибольшим покрытием состоит всего из одного узла. Вообще, есть очень много последовательностей с покрытием сильно меньше среднего, потому что этот выброс сильно повлиял на среднее.
Получаю последовательности самых длинных последовательностей:
path = 'contigs.fa'
flag = False
seq = ''
seqs = []
count = 1
with open(path) as file:
for line in file:
if line.startswith('>'):
if seq:
with open(f'contig{count}.fasta', 'w') as out:
out.write(seq)
count += 1
seq = ''
if line.split('_')[1] in longest:
seq += line
flag = True
else:
flag = False
elif flag:
seq += line
В итоге получаю три файла: contig1.fasta, contig2.fasta и contig3.fasta, содержащте три самых длинных контига.
Теперь нужно выровнять с хромосомой Buchnera aphidicola (GenBank/EMBL AC — CP009253). Я использовала megsblast.
В собранном контиге и банковском геноме были выбраны разные направления в качестве положительных.
| Координаты хромосомы | Координаты контига | % идентичности | Количество мисматчей | Количество инделей |
|---|---|---|---|---|
| 467412-474667 | 62810-55499 | 77.027 | 1491 | 171 |
| 500370-508806 | 29444-20999 | 75.609 | 1756 | 261 |
| 510441-516539 | 19305-13130 | 78.569 | 1143 | 147 |
| 523117-528679 | 6125-550 | 76.873 | 1105 | 159 |
| 462496-467424 | 67781-62837 | 76.987 | 991 | 134 |
| 481997-488128 | 48162-42085 | 74.074 | 1306 | 249 |
| 474844-480660 | 55396-49521 | 74.189 | 1280 | 209 |
| 517766-521500 | 11869-8140 | 77.208 | 763 | 78 |
| 496111-500325 | 33838-29560 | 75.306 | 912 | 123 |
| 493487-494864 | 36632-35256 | 80.144 | 260 | 15 |
| 480874-481548 | 49315-48633 | 82.096 | 107 | 15 |
| 528794-529211 | 480-75 | 84.038 | 40 | 22 |
| 495033-495148 | 35122-35004 | 90.000 | 7 | 4 |
Контиг покрыл участок, который выбран как начало координат в банковском геноме (он кольцевой).
| Координаты хромосомы | Координаты контига | % идентичности | Количество мисматчей | Количество инделей |
|---|---|---|---|---|
| 2004-11103 | 14875-23964 | 78.388 | 1748 | 190 |
| 615561-620926 | 49-5431 | 77.860 | 1077 | 104 |
| 621055-627104 | 5661-11710 | 75.830 | 1246 | 208 |
| 13994-14465 | 26974-27448 | 82.008 | 77 | 8 |
Снова были выбраны разные направления в качестве положительных.
| Координаты хромосомы | Координаты контига | % идентичности | Количество мисматчей | Количество инделей |
|---|---|---|---|---|
| 101718-108876 | 23620-16459 | 76.599 | 1484 | 186 |