Маппирование чтений на референс
1. Картирование чтений на референсный геном
В этом практикуме мы должны были картировать почищенные парные (прямые и обратные) чтения на референс и подготовить их к поиску вариантов.
Выполняем с помощью команды:
bwa mem -t 10 chr9.fna ../trimming_dna/dna_f_p.fastq.gz ../trimming_dna/dna_r_p.fastq.gz > dna.sam 2> bwa_mem_log.txt (
лог)
В начале указали референс, затем два файла с парными чтениями, выход перенаправили в
.sam файл и, через
перенаправление из потока ошибок, сохраняем лог-файл, на всякий случай. Число ядер для ускорения процесса задано с помощью
-t.
2. Анализ sam файла с выравниванием
a. SAM-формат (Sequence Alignment Map) - текстовый формат файла, содержащий
биологические последовательности, выровненные по референсу.
Формат SAM состоит из заголовка и секции выравнивания.
Заголовки начинаются с символа «@», который отличает их от раздела выравнивания.
Разделы выравнивания содержат 11 обязательных полей, а также переменное количество дополнительных полей. Попробуем посмотреть, как выглядят первые 10 строк нашего sam файла
с помощью команды
head dna.sam (
файл):
@SQ SN:NC_000009.12 LN:138394717
@PG ID:bwa PN:bwa VN:0.7.17-r1188 CL:bwa mem chr9.fna ../trimming_dna/dna_f_p.fastq.gz ../trimming_dna/dna_r_p.fastq.gz
SRR10720412.2 77 * 0 0 * * 0 0 GCTGAGTGATGGAGGGGGATAAGTCTTCTCGATAATAAACCTTCCCCAGGCCATCGGTCGTGTCCAAGTCCGCCA GGGD@F>EEBDDFFDDEEDEGGEEGGBGGFGEGBBGDGDDGGGGGGEDGGGEF?EGIGGEGIIGIEGBEGEEEEEGFFFFICC#CCB?=AB AS:i:0 XS:i:0
SRR10720412.3 141 * 0 0 * * 0 0 TGTGTTCCTTTTGCTTTGTNGTGTCTCCTGAGCAACCGCCNNTNNNNACAGTAAACATATTTTCTTTAATTATTA GEGEGIGIHHDIHIGFFFB#C=CAAIIIGIBHIIHEFCEF##?####498A=?AAIIIIIE774#####8>##7889#4##<<7><#######0<>::@<=GEGIFIGHFBDGGDDIHIHI AS:i:0 XS:i:0
SRR10720412.6 77 * 0 0 * * 0 0 CATTTTTATTGCATTTTTGCAGTTNNTNCCATTTATTGAAAATNACAATTGTAACAATTTTCTTA HHHHHHEHHHHHDHHHHHHH>@???BBGHHEHEFB#EDBBBBHEHHHHFHHHHHHHH AS:i:0 XS:i:0
SRR10720412.6 141 * 0 0 * * 0 0 ATGGTTTCTTTACAGAANNNNNNTCNNCACCNTNNNTTCGNNNNNNNNNCTGTTCGAGGTGGAGTGACAGGACGT GGGGFIIIIIIIIHI88######59##?<:<#<###<8;<#########657559DDGG@GFED@DDEBE@;FA< AS:i:0 XS:i:0
b.
Таблица 1. Описание первых 12 столбцов тела файла в формате SAM
| № |
Поле |
Тип данных |
Описание |
| 1 |
QNAME |
String |
Название рида (парным будет присвоено одинаковое) |
| 2 |
FLAG |
Int |
Представляет собой комбинацию нескольких бинарных флагов, обозначающих парность, картированность и т. д. |
| 3 |
RNAME |
String |
Название референса |
| 4 |
POS |
Int |
Один минус самая левая позиция выравнивания |
| 5 |
MAPQ |
Int |
Качество маппинга. MAPQ = −10 log10P, где P - вероятность того, что выравнивание неверно |
| 6 |
CIGAR |
String |
Описание выравнивания, в записи которого используется набор операций (совпадение, инсерция и др.) |
| 7 |
RNEXT |
String |
Название референса пары |
| 8 |
PNEXT |
Int |
Позиция в референсе, которой соответствует начало пары. Аналогично четвёртому полю. |
| 9 |
TLEN |
Int |
Расстояние между крайними точками выравнивания |
| 10 |
SEQ |
String |
Последовательность рида (регистр не имеет значения) |
| 11 |
QUAL |
String |
Качество в кодировке ASCII по знакомой нам шкале Phred+33 |
| 12 |
TAGs |
Mixed |
Специальные тэги, предназначенные для хранения информации о прочтении или выравнивании |
c. С помощью команды
ls -sh выяснили, что наш
.sam файл весит 16 Гб.
3. Получение и индексирование bam файла
a. SAM файл очень тяжелый.
Двоичным эквивалентом формата SAM является формат BAM (Binary Alignment Map), в котором те же данные хранятся в сжатом двоичном представлении
Чтобы получить сортированный по самой левой координате
.bam файл, переконвертируем его с помощью команды:
samtools sort -O bam -@ 10 dna.sam -o dna.bam
где
-O bam (outformat) задает формат выходного файла, а
-o (output) - файл выхода,
а
-@ задаем число ядер на побольше, чтобы считалось быстрее.
b. С помощью команды
ls -sh выяснили, что наш
.bam файл весит 4,7 Гб.
c. BAM-файл тоже можно использовать для работы, но для этого его нужно проиндексировать с помощью команды:
samtools index -b -@ 10 dna.bam dna.bai
где
-b указывает на формат
.bai, хотя и этот формат в любом случае установлен по умолчанию.
4. Анализ bam файла
a. Заглянуть в bam файл просто так не получится - он бинарный.
Попробуем получить статистику о файле с помощью команды:
samtools flagstat -@ 10 dna.bam > dna_flag.txt (
файл)
b. В файле мы видим 12,6% картированных чтений, но не все из них корректно картированы.
Мы ожидаем, что чтения из одной пары должны картироваться недалеко друг от друга и быть направлены друг к другу - таких обнаруживаем 9,67%.
Число небольшое, но объясняется оно тем, что мы картируем данные полногеномного секвенирования, картированные на одну из хромосом.
5. Получение картированных чтений
a. Попробуем вытащить из выравнивания чтения, картированные только на нашу хромосому.
Сначала, проиндексируем fasta-файл с хромосомой с помощью команды (из файла
.fai узнаем название хромосомы):
samtools faidx chr9.fna -o chr9.fai
Теперь используем:
samtools view -h dna.bam NC_000009.12 > dna.chr9.sam
samtools view -bS dna.chr9.sam > dna.chr9.bam - конвертируем в bam
b. Посмотрим же теперь статистику полученных ридов из нашей хромосомы:
samtools flagstat -@ 10 dna.chr9.bam > dna_chr9_flag.txt (
файл)
Теперь у нас 89,42% картированных чтений, и 72,23% корректно картированных
c. Дела стали лучше, теперь мы имеем большую часть ридов, относящихся к нашей хромосоме, хотя и не все.
d. Почему мы получили чтения, картированные на нужную нам хромосому, но при этом количество картированных чтений не 100%?
По всей видимости, samtools плохо распознает некорректно картированные чтения из-за дупликаций и наложений при маппинге.
e. Теперь получим только правильно картированные пары чтений с помощью команды:
samtools view -@ 10 -f 0x2 -bS dna.chr9.bam > dna_chr9_correct.bam
f. Снова посмотрим статистику:
samtools flagstat -@ 10 dna_chr9_correct.bam > dna_chr9_correct_flag.txt (
файл)
Теперь у нас 100% картированных, и 100% корректно картированных.
g. Теперь у нас есть файл с корректно картированными чтениями, samtools видит их правильно, поэтому логично, что их теперь 100%.
h. Проиндексируем наш файлик с ридами для дальнейшей работы:
samtools index -b -@ 10 dna_chr9_correct.bam dna_chr9_correct.bai
6. Подготовка выравнивания к поиску вариантов
a. Маркируем дублированные чтения с помощью команды:
picard MarkDuplicates -M metrix.file.txt -I dna_chr9_correct.bam -O markdup.bam
b. Посмотрим статистику дублированных чтений:
samtools flagstat -@ 10 markdup.bam > markdup.txt (
файл)
Дублированных чтений оказалось 257 010.
Variant calling
В этом практикуме из картированных на референс прочетний нужно было получить, отфильтровать и минимально исследовать варианты, а также поработать с покрытием.
1. Получение вариантов
a. Из bam файла с маркированными дублированными чтениями получим файл с вариантами с помощью bcftools:
bcftools mpileup --threads 10 -f ../bwa/chr9.fna ../bwa/markdup.bam | bcftools call -mv -o vc_chr9.vcf
Команда mpileup создает
.vcf файл с частотами встречаемости вариантов во множественном выравнивании.
Число ядер задаем как -threads, индексированный референс хромосомы как -f, после чего передаем выход
команде call, которая осуществляет "вызов" вариантов, с параметрами
-m - модель вызова для многоаллельных и редких вариантов,
-v - запись в файл только сайтов с вариантами,
и -o - выходной файл.
b. Рассмотрим первые 35 строчек файла с помощью
head -35 vc_chr9.vcf (некоторые строки я убрал для лакончиности):
##fileformat=VCFv4.2
##FILTER=
##bcftoolsVersion=1.9+htslib-1.9
##bcftoolsCommand=mpileup --threads 10 -f ../bwa/chr9.fna ../bwa/markdup.bam
##reference=file://../bwa/chr9.fna
##contig=
##ALT=
##INFO=
...
##FORMAT=
...
##INFO=
...
##bcftools_callVersion=1.9+htslib-1.9
##bcftools_callCommand=call -mv -o vc_chr9.vcf; Date=Wed Dec 2 11:08:14 2020
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT ../bwa/markdup.bam
NC_000009.12 10111 . C A 13.7431 . DP=15;VDB=0.327477;SGB=-0.511536;RPB=0.400803;MQB=1;MQSB=0.914584;BQB=0.914584;MQ0F=0.133333;ICB=1;HOB=0.5;AC=1;AN=2;DP4=7,3,3,0;MQ=30 GT:PL 0/1:47,0,142
NC_000009.12 10465 . C A 10.7923 . DP=2;SGB=-0.379885;MQ0F=0;AC=2;AN=2;DP4=0,0,1,0;MQ=60 GT:PL 1/1:40,3,0
NC_000009.12 10469 . G C 10.7923 . DP=1;SGB=-0.379885;MQ0F=0;AC=2;AN=2;DP4=0,0,1,0;MQ=60 GT:PL 1/1:40,3,0
NC_000009.12 10473 . G C 10.7923 . DP=1;SGB=-0.379885;MQ0F=0;AC=2;AN=2;DP4=0,0,1,0;MQ=60 GT:PL 1/1:40,3,0
NC_000009.12 10482 . G C 3.18683 . DP=2;SGB=-0.379885;RPB=1;MQB=1;BQB=1;MQ0F=0;ICB=1;HOB=0.5;AC=1;AN=2;DP4=1,0,1,0;MQ=60 GT:PL 0/1:33,0,19
NC_000009.12 10506 . A G 27.4023 . DP=5;VDB=0.18;SGB=-0.453602;RPB=0;MQB=0;BQB=0.75;MQ0F=0;ICB=1;HOB=0.5;AC=1;AN=2;DP4=2,0,2,0;MQ=56 GT:PL 0/1:60,0,54
VCF формат - формат записи вариантов последовательностей. Заголовок начинается с названия и содержит метаданные, описывающие тело файла.
Строки заголовков начинаются с символа #. Специальные ключевые слова в заголовке обозначаются ##.
Немного упомянём и поля, присутствующие в файле.
Таблица 2. Описание первых 9 столбцов файла в формате VCF
| № |
Название |
Описание |
| 1 |
CHROM |
Имя последовательности (обычно хромосомы), в которой вызывается вариант |
| 2 |
POS |
Один минус позиция вариант |
| 3 |
ID |
Индивидуальный номер варианта |
| 4 |
REF |
Референсный нуклеотид (-ы, в случае инделя) |
| 5 |
ALT |
Список аллельных вариантов |
| 6 |
QUAL |
Качество варианта, в зависимости от частот аллелей |
| 7 |
FILTER |
Флаг, указывающий, какой из заданного набора фильтров прошёл вариант |
| 8 |
INFO |
Ключевые значения для описания варианта |
| 9 |
FORMAT |
Описание образца |
c. Посмотрим статистику vcf файла с помощью:
bcftools stats vc_chr9.vcf > vc_chr9_stats.txt (
файл)
d. Исходя из файла, мы видим, что всего обнаружено 151 464 варианта, из которых 148 959 являются однонуклеотидными заменами (SNP), а 2505 - вставкам и делециями.
2. Фильтрация вариантов
a. Отфильтруем варианты по качеству больше 30 и по покрытию больше 50:
bcftools filter -i'%QUAL>30 && DP>50' --threads 10 vc_chr9.vcf -o vc_filtered.vcf
b. Посмотрим статистику вариантов после фильтрации:
bcftools stats vc_filtered.vcf > vc_filtered_stats.txt (
файл)
Исходя из файла, мы видим, что всего обнаружено 10 316 вариантов, из которых 10 172 являются однонуклеотидными заменами (SNP), а 144 - вставкам и делециями.
Логично, что после фильтрации их стало сильно меньше, зато среди них точно есть функционально значимые.
3. Изучение покрытия
Однако, поскольку перед нами не стояло задачи исследовать конкретные участки, просто посмотрим на хорошо покрытые.
a. Получим участки генома, покрытые нашими чтениями хотя бы 1 раз. В качестве входного bam файла взяли файл с маркированными дублированными чтениями,
с опцией -bg, чтобы получить выход в формате BedGraph:
bedtools genomecov -bg -ibam ../bwa/markdup.bam > gencov.bed (
файл)
b. Формат BedGraph представлен в данном случае в виде таблички из 4 колонок:
название референса (в нашем случае хромосома), начало и конец варианта (обычно с разницей в 1), и покрытие
c. Отберите только участки генома, покрытые более 50 раз, а затем для удобства отсортируем их по убыванию покрытия:
awk '{if ($4 > 50) {print $0}}' gencov.bed | sort -k4,4nr > gencov_geq50.txt (
файл)
d. Выясним, сколько таких вариантов получилось с помощью подсчёта строк:
wc -l gencov_geq50.txt - их получилось 593 932
Для сравнения, у нас было
wc -l wc -l gencov.bed - 2 436 694 покрытых хотя бы единожды вариантов, то есть примерно пятая часть покрыта более чем 50 раз.
e. Теперь поищем варианты, которые попали в хорошо покрытые участки, подав на вход отфильтрованный vcf файл из пункта 2а и bed файл из пункта 3с:
bedtools intersect -a vc_filtered.vcf -b gencov.bed > intersect.txt (
файл)
wc -l intersect.txt - всего таких вариантов 10 723.
Помогите, я устал писать команды.