Практикум 14. Сборка de novo
Код доступа: SRR1722713
Команда для скачивания файла
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR172/003/SRR1722713/SRR1722713.fastq.gz
Исходный файл имеет вес 821 МБ
Рабочая директория: /mnt/scratch/NGS/sevartem/pr14
1. Подготовка чтений программой trimmomatic
Я объединил файлы всех адаптеров Illumina в один файл adapters_all.fasta:
cat /mnt/scratch/NGS/adapters/*.fa > adapters_all.fasta
Первый запуск TrimmomaticSE:
TrimmomaticSE -threads 8 -phred33 SRR1722713.fastq.gz SRR1722713_no_adapters.fastq.gz ILLUMINACLIP:adapters_all.fasta:2:7:7 -summary summary1.txt
Полученный файл SRR1722713_no_adapters.fastq.gz имеет вес 805 МБ
Из 13241431 осталось 13235572, то есть было отброшено 5859 чтений (0.04%)
Второй запуск TrimmomaticSE:
TrimmomaticSE -threads 8 -phred33 SRR1722713_no_adapters.fastq.gz SRR1722713_trimmed.fastq.gz TRAILING:20 MINLEN:32 -summary summary2.txt
Полученный файл SRR1722713_trimmed.fastq.gz имеет вес 773 МБ
Из 13235572 чтений осталось 12767826, то есть было отброшено 467746 чтений (3.53%)й
Логи запусков сохранены в файлах log1.txt и log2.txt. Summary cохранены в файлах summary1.txt и summary2.txt
Параметры trimmomatic:
- -threads 8 – работа в 8 потоках для большей скорости
- -phred33 – кодировка
- ILLUMINACLIP:adapters_all.fasta:2:7:7 – здесь:
- 2 – максимальное число несовпадений при поиске адаптера
- 7 – порог для палиндромного выравнивания (здесь он игнорируется, так как выравнивание не парное)
- 7 – порог для простого выравнивания: чем он выше, тем строже поиск адаптера. Значение 7 является умеренным
- TRAILING:20 – удаляет нуклеотиды с качеством ниже 20 с правого конца чтения
- MINLEN:32 – чтения короче 32 нуклеотидов отбрасываются
Стоит отметить, что в логе присутствуют сообщения Skipping duplicate Clipping Sequence. Это значит, что некоторые адаптеры были одинаковыми, однако это не влияет на результат
Статистика:
- Входные чтения: 13241431
- Оставшиеся: 12767826 (96.42%)
- Отброшенные: 473605 (3.58%)
2. Сборка транскриптома de novo (Velvet)
2.1. Построение графа де Брёйна (velveth)
Применена команда
velveth velveth_output 31 -short -fastq.gz SRR1722713_trimmed.fastq.gz
В ходе ее работы была создана директория velveth_output с файлами Log, Roadmaps и Sequences
Пояснение параметров:
- velveth_output – имя выходной директории
- 31 – длина k-мера (hash_length)
- -short – указание, что чтения короткие и одиночные (single-end), а не парные. (для выдачи коротких парных чтений используется отдельная опция shortPaired)
- -fastq.gz – формат входного файла (сжатый fastq)
- SRR1722713_trimmed.fastq.gz – имя входного файла
2.2. Сборка контигов (velvethg)
Применена команда
velvetg velveth_output
В данном случае дополнительные опции не требуются, поскольку
- Чтения одиночные, и опции для парных концов (-ins_length и связанные) не требуются
- Цель этапа – первичная сборка для оценки характеристик без предварительной фильтрации
Результат сборки – файл velveth_output/contigs.fa.
N50 = 126. Найдено с помощью скрипта Python. Команда:
python get_N50.py contigs.fa
С помощью другого скрипта top_contigs.py я нашел топ-3 самых длинных контига. Команда для запуска:
python top_contigs.py contigs.fa 3
Самые длинные контиги:
- NODE_34306: длина 2210, покрытие 9.159728
- NODE_57331: длина 1899, покрытие 7.360716
- NODE_17779: длина 1673, покрытие 8.439928
Среднее покрытие по ним: 8.320124
С помощью третьего скрипта anomal_contigs.py я нашел 680 контигов с аномально высокими и 16959 контигов аномально низкими покрытиями (общее число контигов: 53145, то есть доля аномальных составляет 33.19%)
Один из контигов с аномально большим покрытием NODE_44186_length_39_cov_1280.051270 находится посередине гена psaA (CDS: NP_051059.1) с функцией photosystem I P700 chlorophyll a apoprotein A1
Один из контигов с аномально малым покрытием >NODE_138268_length_62_cov_1.048387 находится в начале гена BGLU16 (CDS: NP_191572.1) с функцией beta glucosidase 16
3. Анализ полученных контигов
Я провел BLAST для трех самых длинных контигов по базе данных RefSeq Genome database по таксону Arabidopsis thaliana (taxid:3702). Параметры алгоритмов оставались стандарнтыми
3.1. BLAST контига NODE_34306
Длина 2210, покрытие 9.159728
При стандартной длине слова (28) не было получено результатов. При уменьшении длины слова до 16 был получен один результат с довольно большим E-value и очень низким покрытием
Таблица 1. Выдача BLAST для контига NODE_34306
| Max Score | Total Score | Query Cover | E value | Per. Ident | Acc. Len | Accession |
|---|---|---|---|---|---|---|
| 45.4 | 45.4 | 2% | 0.006 | 91.18% | 23459830 | NC_003074.8 |
Вывод: контиг NODE_34306 не показывает значимого сходства с хлоропластным геномом. Найденный фрагмент длиной 39 п.н. с 91% идентичности может быть либо результатом случайного совпадения, либо указывать на наличие консервативного мотива в ядерном гене. Основная часть контига (98%) не имеет гомологии в референсном геноме A. thaliana при использованных параметрах BLAST.
Рис. 1. Карта локального сходства для BLAST контига NODE_34306
3.2. BLAST контига NODE_57331
Длина 1899, покрытие 7.360716
При стандартной длине слова (28) не было получено результатов. При уменьшении длины слова до 16 и увеличении expect trashhold до 1.00 было получено 5 результатов с довольно большими E-value и очень низким покрытием
Таблица 2. Выдача BLAST для контига NODE_57331
| Max Score | Total Score | Query Cover | E value | Per. Ident | Acc. Len | Accession |
|---|---|---|---|---|---|---|
| 39.9 | 39.9 | 2% | 0.22 | 86.49% | 23459830 | NC_003074.8 |
| 39.9 | 192 | 3% | 0.22 | 85.37% | 30427671 | NC_003070.9 |
| 38.1 | 76.1 | 3% | 0.80 | 87.88% | 26975502 | NC_003076.8 |
| 38.1 | 38.1 | 1% | 0.80 | 100.00% | 18585056 | NC_003075.7 |
| 38.1 | 38.1 | 1% | 0.80 | 90.00% | 19698289 | NC_003071.7 |
Выбранная находка для DotPlot: NC_003070.9:12619675-12619637
Вывод: контиг NODE_57331 демонстрирует слабое сходство с межгенным участком хромосомы 1, а также с другими хромосомами (2, 3, 5). Это типично для повторяющихся или малоконсервативных последовательностей. Хлоропластное происхождение не подтверждается.
Рис. 3. Карта локального сходства для BLAST контига NODE_57331
3.3. BLAST контига NODE_17779
Длина 1673, покрытие 8.439928
ри стандартной длине слова (28) не было получено результатов. При уменьшении длины слова до 16 и увеличении expect trashhold до 1.00 было получено 2 результата с довольно большими E-value и очень низким покрытием
Таблица 3. Выдача BLAST для контига NODE_17779
| Max Score | Total Score | Query Cover | E value | Per. Ident | Acc. Len | Accession |
|---|---|---|---|---|---|---|
| 39.9 | 78.0 | 2% | 0.20 | 90.62% | 18585056 | NC_003075.7 |
| 38.1 | 38.1 | 1% | 0.70 | 100.00% | 23459830 | NC_003074.8 |
Выбранная находка для DotPlot: NC_003075.7:3802699-3802668
Вывод: Как и в предыдущих случаях, контиг NODE_17779 имеет лишь очень короткий участок сходства с ядерным геномом, не относящимся к хлоропласту. Основная часть контига (98%) остаётся неидентифицированной.
Рис. 3. Карта локального сходства для BLAST контига NODE_17779
Вывод
В ходе практикума выполнена сборка de novo транскриптома хлоропластов Arabidopsis thaliana на основе одиночных прочтений Illumina (SRR1722713). После фильтрации Trimmomatic (удалено 3,58% чтений) и сборки Velvet (k=31) получено 53 145 контигов с N50=126. Три самых длинных контига (длиной 2210, 1899 и 1673 п.н.) при BLAST-анализе не проявили значимого сходства с хлоропластным геномом, демонстрируя лишь короткие выравнивания на ядерные хромосомы, что, вероятно, связано с фрагментированностью сборки и неоптимальными параметрами. Тем не менее, анализ контигов с экстремальным покрытием позволил идентифицировать короткий фрагмент (39 п.н.) гена psaA фотосистемы I с покрытием 1280×, что подтверждает присутствие высокоэкспрессируемых хлоропластных транскриптов в исходных данных. Таким образом, основные этапы обработки данных NGS и сборки de novo освоены, а полученные результаты иллюстрируют как потенциальные возможности метода, так и его ограничения при работе с реальными транскриптомными данными.