Учебный сайт Левина Ильи, 3-й семестр

Подготовка данных NGS для картирования

Задание 1: Описание файла с референсом

Заглянув в свой файл с референсом, я с удивлением обнаружил, что он очень похож на обычный fasta-файл с нуклеотидной последовательностью. После символа ">" в первой строке файла идёт, видимо, RefSeq AC (или АС какой-нибудь другой базы данных), после чего идёт описание последовательности: у меня оказалась первичная сборка 3-й хромосомы человека. Далее, после этой строки идёт ОЧЕНЬ МНОГО строк самой последовательности. Стоит отметить, что первые строк 80-90 состоят исключительно из букв "N":

>NC_000003.12 Homo sapiens chromosome 3, GRCh38.p13 Primary Assembly
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
. . .
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
CTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTCACCCTACCCTAACCCTAACCCTAACCCTC
ACCCTCACCCTCACTCAACCCTAACCCTCACCCTAACCTCCACCCTAACCCTAACCCCTAACCCTAACCCTAACCCTAAC
. . .

Задание 2: Индексация референса

Делается это с помощью инструмента bwa - Burrows-Wheeler Alignment Tool, то есть инструмента выравнивания. В нём мы обращаемся к команде index с параметром -a и опцией bwtsw, то есть мы хотим проиндексировать наш референс, используя алгоритм BWT-SW для построения индекса, так как именно он [алгоритм] работает по всему геному человека, что нам и нужно. Так выглядит сессия индексации в командной строке:

$ bwa index -a bwtsw chr3.fna
[bwa_index] Pack FASTA... 1.81 sec
[bwa_index] Construct BWT for the packed sequence...
[BWTIncCreate] textLength=396591118, availableWord=39905300
[BWTIncConstructFromPacked] 10 iterations done. 65826158 characters processed.
[BWTIncConstructFromPacked] 20 iterations done. 121609790 characters processed.
[BWTIncConstructFromPacked] 30 iterations done. 171186446 characters processed.
[BWTIncConstructFromPacked] 40 iterations done. 215246238 characters processed.
[BWTIncConstructFromPacked] 50 iterations done. 254402686 characters processed.
[BWTIncConstructFromPacked] 60 iterations done. 289201022 characters processed.
[BWTIncConstructFromPacked] 70 iterations done. 320125822 characters processed.
[BWTIncConstructFromPacked] 80 iterations done. 347607806 characters processed.
[BWTIncConstructFromPacked] 90 iterations done. 372029854 characters processed.
[BWTIncConstructFromPacked] 100 iterations done. 393732190 characters processed.
[bwt_gen] Finished constructing BWT in 102 iterations.
[bwa_index] 115.72 seconds elapse.
[bwa_index] Update BWT... 1.15 sec
[bwa_index] Pack forward-only FASTA... 1.11 sec
[bwa_index] Construct SA from BWT and Occ... 44.91 sec
[main] Version: 0.7.17-r1188
[main] CMD: bwa index -a bwtsw chr3.fna
[main] Real time: 170.674 sec; CPU: 164.697 sec

Получается, на вход мы подали наш референс, а на выход получили его проиндексированную версию в виде 5-и файлов разного назначения. Ниже будет приведён состав директории после использования bwa index:

$ ls -la
total 534984
drwxr-xr-x 2 lewis@kodomo 1020      4096 Dec 26 16:42 .
drwxr-xr-x 3 lewis@kodomo 1020      4096 Dec 26 16:11 ..
-rw-r--r-- 1 lewis@kodomo 1020 200774323 Dec 26 16:14 chr3.fna
-rw-r--r-- 1 lewis@kodomo 1020       435 Dec 26 16:42 chr3.fna.amb
-rw-r--r-- 1 lewis@kodomo 1020       100 Dec 26 16:42 chr3.fna.ann
-rw-r--r-- 1 lewis@kodomo 1020 198295660 Dec 26 16:42 chr3.fna.bwt
-rw-r--r-- 1 lewis@kodomo 1020  49573891 Dec 26 16:42 chr3.fna.pac
-rw-r--r-- 1 lewis@kodomo 1020  99147832 Dec 26 16:42 chr3.fna.sa

Задание 3: Описание образца

Таблица 1. Описание образца SRA
Критерий	Описание
Ссылка	Тут
Прибор	Illumina Genome Analyzer IIx
Организм	Homo sapiens
Стратегия секвенирования	OTHER
Парно-/непарноконцевые чтения	Парноконцевые
Сколько чтений ожидается	41,398,792

Задание 4: Проверка качества исходных чтений

Для проверки качества чтений я использовал программу fastQC на kodomo, на узле calc. Применял я её таким образом:

$ fastqc SRR10720412_1.fastq.gz
    
$ fastqc SRR10720412_2.fastq.gz

Одним из выходных файлов этой программы был html-файл, который я перенёс себе на компьютер и смотрел его уже у себя.

Количество пар чтений	41,398,792
Совпадает ли количество "прямых" и "обратных" чтений	ДА
Совпадает ли количество чтений с ожидаемым	ДА

Per sequence quality base

Прямые чтения:

Обратные чтения:

Ну, собственно, как мы можем заметить, чтения нам попались с довольно хорошим качеством, все в "зелёной" зоне, хоть и есть такие, качество которых доходит и до середины желтой зоны. При этом среденее качество прочтения нуклеотида от 1-го до последнего убывает от 40 до 35. Такой предварительный взгляд на качество чтений говорит о том, что изначально у нас чтения очень хорошего качество, но всё ещё нуждаются в дополнительной обработке.

Per sequence quality scores

Прямые чтения:

Обратные чтения:

Как мы можем судить по этим графикам распределения, что у "прямых", что у "обратных" чтений больше всего штук с качеством прочтения 39. При этом чтений с качеством прочтения меньше 22 уже пренебрежительно мало, что нам в очередной раз даёт понять, что чтения у нас получились качественные.

Sequence Length Distribution

Прямые чтения:

Обратные чтения:

Как мы можем заметить, распределение ридов по длине у "прямых" и "обратных" одинаково - 75 нуклеотидов. Таким образом мы можем сказать, что все риды у нас получились по 75 нуклеотидов.

Задание 5: Фильтрация чтений

Для фильтрации чтений я использовал программу trimmomatic:

$ java -jar /usr/share/java/trimmomatic.jar PE -phred33 SRR10720412_1.fastq.gz SRR10720412_2.fastq.gz 
    SRR10720412_1_paired.fq.gz SRR10720412_1_unpaired.fq.gz SRR10720412_2_paired.fq.gz SRR10720412_2_unpaired.fq.gz 
    ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 TRAILING:20 MINLEN:50

Что значит каждый параметр:

PE - параметр, поясняющий программе, что у нас парноконцевые чтения (pair-ended);
-phred33 - этим параметром мы выбираем одну из форм записи качества прочтения нуклеотида, у наших чтений она [форма записи] соответствует Phred+33;
Далее 2-мя аргументами у нас идут имена файлов с прямыми и обратными чтениями, а следующие за ними 4 - имена выходных файлов;
ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 - параметр, задающий удаление адаптеров Illumina (если они ещё вдруг остались);
TRAILING:20 - удаляет с конца чтений нуклеотиды с качеством ниже 20;
MINLEN:50 - "оставляет в живых" только чтения, длина которых получилась больше 50 нуклеотидов.

Задание 6: Проверка качества триммированных чтений

Для проверки качества триммированных чтений я вновь использовал программу fastQC. Полученные в HTML результаты я скопировал себе на компьютер и дальше рассматривал на компьютере. Использовал я программу на узле calc таким образом:

$ fastqc SRR10720412_1_paired.fq.gz
    
$ fastqc SRR10720412_1_unpaired.fq.gz

$ fastqc SRR10720412_2_paired.fq.gz

$ fastqc SRR10720412_2_unpaired.fq.gz

Пар чтений осталось: 39844329 (96.25%);

Per sequence quality base

Paired

Прямые чтения:

Обратные чтения:

Unpaired

Прямые чтения:

Обратные чтения:

Как мы можем заметить по этим картинкам, paired чтения прочитались секвенатором лучше и точнее, нежели чем unpaired чтения, так как качества прочтения нуклеотидов из unpaired чтений распределены шире, некоторые отдельные нуклеотиды прочитаны даже хуже 20. При этом же довольно много прочтений находятся в "жёлтой" зоне, но всё же основная масса прочтений находится в "зелёной" зоне. В это же время АБСОЛЮТНО ВСЕ нуклеотиды из paired чтений прочитаны просто идеально, все нуклеотиды в "зелёной" зоне.

В это же время качество paired чтений лучше, чем качество изначальных чтений. Это видно по тому, что ширина блоков уменьшилась, не осталось ни одного нуклеотида в "жёлтой" зоне, у всех качество прочтения 30 и более, в то время как у изначальных чтений самое низкое качество прочтения нуклеотида было 24.

Per sequence quality scores

Paired

Прямые чтения:

Обратные чтения:

От изначальных распределений эти графики отличаются лишь тем, что ось абсцисс сокращена: на изначальных распределениях она начиналась с качества 2, а на этих - с 18 у прямых и 13 у обратных. Это нам говорит, видимо, о том, что trimmomatic как раз поудалял эти чтения, оставив более качественные.

Sequence Length Distribution

Paired

Прямые чтения:

Обратные чтения:

Наибольшее количество чтений всё ещё имеет длину 75 нуклеотидов, но при этом у нас появилось достаточно много более коротких чтений с длинами от 49 до 74 нуклеотидов. Это всё последствия работы trimmomatic: он удалил достаточно много нуклеотидов с плохими прочтениями, что и стало причиной появления более коротких чтений. Изначально ВСЕ чтения были длиной 75 нуклеотидов.