Отчёт по практикуму 14 третьего семестра. Задания по сборке de novo
Подготовка чтений программой Trimmomatic
Для начала я скачал одноконцевые чтения РНК из хлоропластов резуховидки Таля (Arabidopsis thaliana) при помощи команды:
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR172/008/SRR1724088/SRR1724088.fastq.gz
Адаптеры для одноконцевых чтений были объединены в один файл:
cat TruSeq2-SE.fa TruSeq3-SE.fa > adapters.fa
Сперва программой fastqc были получены данные о чтениях без удаления остатка адаптеров:
fastqc SRR1724088.fastq.gz
Затем программой Trimmotatic были удалены адаптеры одноконцевых чтений:
TrimmomaticSE -phred33 SRR1724088.fastq.gz output_trimmed.fastq.gz ILLUMINACLIP:adapters.fa:2:7:7
Таким образом было удалено ~ 3,4% чтений от изначального файла. Размер исходного файла - 733M, а полностью обработанного - 691M.
Работа программы Velveth
Подготовка чтений для дальнейшего анализа осуществлялась программой velveth, k-меры = 31. Чтения - короткие и непарные.
velveth . 31 -short -fastq output_trimmed_trailed.fastq.gz
Получение контигов
Команда velvetg . -exp_cov auto -cov_cutoff 70 собственно создаёт контиги на основе генома.
При помощи скрипта:
#!/usr/bin/env python3
import statistics
contigs = []
with open("contigs.fa") as f:
for line in f:
if line.startswith(">"):
header = line.strip()
parts = header.split("_")
length = int(parts[3])
cov = float(parts[5])
contigs.append((length, cov, header))
# сортировка по длине (по убыванию)
contigs.sort(reverse=True)
# 3 самых длинных
top3 = contigs[:3]
print("3 самых длинных:")
for length, cov, header in top3:
print(header)
print("length:", length, "coverage:", cov)
print()
# Среднее
mean_top3 = sum(cov for length, cov, header in top3) / 3
print("Среднее покрытие топ-3:", mean_top3)
# N50
lengths = [length for length, cov, header in contigs]
total = sum(lengths)
half = total / 2
running = 0
n50 = None
for length in lengths:
running += length
if running >= half:
n50 = length
break
print("N50:", n50)
# Пороги аномальных контигов
low = mean_top3 / 5
high = mean_top3 * 5
print("Порог низкого покрытия:", low)
print("Порог высокого покрытия:", high)
# Аномальные
anomalous = [
(length, cov, header)
for length, cov, header in contigs
if cov < low or cov > high
]
print("\nВсего аномальных:", len(anomalous))
# 2 самых высоких покртия
print("\n2 контига с самым большим покрытием (аномальные):")
for length, cov, header in sorted(anomalous, key=lambda x: x[1], reverse=True)[:2]:
print(header)
print("length:", length, "coverage:", cov)
print()
# 2 самых низких покрытия
print("2 контига с самым маленьким покрытием (аномальные):")
for length, cov, header in sorted(anomalous, key=lambda x: x[1])[:2]:
print(header)
print("length:", length, "coverage:", cov)
print()
на языке Python были выведены характеристики контигов
3 самых длинных:
>NODE_308_length_5386_cov_576.270691
length: 5386 coverage: 576.270691
>NODE_387_length_2939_cov_302.083710
length: 2939 coverage: 302.08371
>NODE_376_length_2315_cov_713.514038
length: 2315 coverage: 713.514038
Медиана: 530.62
N50: 663
Порог низкого покрытия: 106
Порог высокого покрытия: 2653
Всего аномальных: 77
2 контига с самым большим покрытием (аномальные):
>NODE_236_length_51_cov_137320.718750
length: 51 coverage: 137320.71875
>NODE_50_length_80_cov_40150.261719
length: 80 coverage: 40150.261719
2 контига с самым маленьким покрытием (аномальные):
>NODE_369_length_268_cov_71.156715
length: 268 coverage: 71.156715
>NODE_638_length_86_cov_72.802322
length: 86 coverage: 72.802322
Контиги с аномально малым покрытием возможно попали в чтения из ядра клетки. С аномально большим покрытием связаны чтения хлоропластного генома.
Анализ контигов
Для megablast была взята хлоропластная хромосома Arabidopsis thaliana из RefSeq Genomes.
Самый длинный контиг >NODE_308 не выровнялся на геном. Возможно, он является частью ядерного генома или что-то пошло не так при сборке. Но результат megablast получился нулевым.
Контиг >NODE_387 выровнялся на один участок, результаты выравнивания:
