Практикум 13. Анализ данных NGS: поиск вариантов, изучение покрытия
Этот практикум посвящён поиску вариантов (SNP, инделей) в анализируемых прочтениях, а также их фильтрации по качеству и покрытию.
Получение вариантов
Для этого задания было предложено воспользоваться средствами bcftools. Далее приведена использованная команда:
bcftools mpileup --threads 10 -f chr4.fna correct_markd.bam | bcftools call -mv -o variants_chr4.vcf
Здесь bcftools mpileup - алгоритм, который на выходе даёт vcf-файл с вероятностями генотипов (?), -f chr4.fna - имя референсной последовательности, --threads 10 - параметр, указывающий использовать 10 потоков, correct_markd.bam - файл с корректно картированными прочтениями и маркированными дублями, bcftools call - вторая команда, принимающая на вход выдачу первой и вызывающая варианты, -mv - опции, позволяющие работать с редкими и многоаллельными вариантами и записывать в файл только варианты соответственно, -o variants_chr4.vcf - выходной файл.
Далее приведён фрагмент полученного vcf-файла.
##fileformat=VCFv4.2 ##FILTER=<ID=PASS,Description="All filters passed"> ##bcftoolsVersion=1.9+htslib-1.9 ##bcftoolsCommand=mpileup --threads 10 -f chr4.fna correct_markd.bam ##reference=file://chr4.fna ##contig=<ID=NC_000004.12,length=190214555> ##ALT=<ID=*,Description="Represents allele(s) other than observed."> ##INFO=<ID=INDEL,Number=0,Type=Flag,Description="Indicates that the variant is an INDEL."> ##INFO=<ID=IDV,Number=1,Type=Integer,Description="Maximum number of reads supporting an indel"> ##INFO=<ID=IMF,Number=1,Type=Float,Description="Maximum fraction of reads supporting an indel"> ##INFO=<ID=DP,Number=1,Type=Integer,Description="Raw read depth"> ##INFO=<ID=VDB,Number=1,Type=Float,Description="Variant Distance Bias for filtering splice-site artefacts in RNA-seq data (bigger is better)",Version="3"> ##INFO=<ID=RPB,Number=1,Type=Float,Description="Mann-Whitney U test of Read Position Bias (bigger is better)"> ##INFO=<ID=MQB,Number=1,Type=Float,Description="Mann-Whitney U test of Mapping Quality Bias (bigger is better)"> ##INFO=<ID=BQB,Number=1,Type=Float,Description="Mann-Whitney U test of Base Quality Bias (bigger is better)"> ##INFO=<ID=MQSB,Number=1,Type=Float,Description="Mann-Whitney U test of Mapping Quality vs Strand Bias (bigger is better)"> ##INFO=<ID=SGB,Number=1,Type=Float,Description="Segregation based metric."> ##INFO=<ID=MQ0F,Number=1,Type=Float,Description="Fraction of MQ0 reads (smaller is better)"> ##FORMAT=<ID=PL,Number=G,Type=Integer,Description="List of Phred-scaled genotype likelihoods"> ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype"> ##INFO=<ID=ICB,Number=1,Type=Float,Description="Inbreeding Coefficient Binomial test (bigger is better)"> ##INFO=<ID=HOB,Number=1,Type=Float,Description="Bias in the number of HOMs number (smaller is better)"> ##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes for each ALT allele, in the same order as listed"> ##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes"> ##INFO=<ID=DP4,Number=4,Type=Integer,Description="Number of high-quality ref-forward , ref-reverse, alt-forward and alt-reverse bases"> ##INFO=<ID=MQ,Number=1,Type=Integer,Description="Average mapping quality"> ##bcftools_callVersion=1.9+htslib-1.9 ##bcftools_callCommand=call -mv -o variants_chr4.vcf; Date=Sat Dec 19 22:38:39 2020 #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT correct_markd.bam NC_000004.12 10047 . A G 3.77461 . DP=7;VDB=0.58;SGB=-0.453602;RPB=0.333333;MQB=0.333333;MQSB=1;BQB=0.833333;MQ0F=0.142857;ICB=1;HOB=0.5;AC=1;AN=2;DP4=3,0,1,1;MQ=18 GT:PL 0/1:34,0,29 NC_000004.12 12278 . G C 13.6078 . DP=3;VDB=0.909987;SGB=-0.511536;MQSB=1;MQ0F=0.333333;AC=2;AN=2;DP4=0,0,1,2;MQ=14 GT:PL 1/1:43,9,0 NC_000004.12 12295 . A G 9.88514 . DP=3;VDB=0.878371;SGB=-0.511536;MQ0F=0.333333;AC=2;AN=2;DP4=0,0,0,3;MQ=14 GT:PL 1/1:39,9,0 NC_000004.12 12318 . G A 8.99921 . DP=2;VDB=0.14;SGB=-0.453602;MQ0F=0;AC=2;AN=2;DP4=0,0,0,2;MQ=21 GT:PL 1/1:38,6,0 NC_000004.12 15036 . G A 21.3048 . DP=4;VDB=0.06;SGB=-0.453602;RPB=1;MQB=1;MQSB=1;BQB=1;MQ0F=0.25;AC=2;AN=2;DP4=0,1,1,1;MQ=19 GT:PL 1/1:48,3,0 NC_000004.12 16698 . T C 4.36621 . DP=6;VDB=0.66;SGB=-0.453602;RPB=0.333333;MQB=0.333333;MQSB=0.666667;BQB=0.333333;MQ0F=0.166667;ICB=1;HOB=0.5;AC=1;AN=2;DP4=1,2,2,0;MQ=19 GT:PL 0/1:35,0,30 NC_000004.12 16754 . G A 4.34871 . DP=4;VDB=0.06;SGB=-0.453602;RPB=0;MQB=0.5;MQSB=1;BQB=0.5;MQ0F=0.25;ICB=1;HOB=0.5;AC=1;AN=2;DP4=1,1,2,0;MQ=18 GT:PL 0/1:35,0,17 NC_000004.12 16843 . A G 5.01626 . DP=4;VDB=0.9;SGB=-0.453602;RPB=1;MQB=0.25;BQB=0.5;MQ0F=0;ICB=1;HOB=0.5;AC=1;AN=2;DP4=0,2,0,2;MQ=32 GT:PL 0/1:36,0,28
Файл vcf состоит из строк заголовка, начинающихся с «#», и расположенных ниже записей в 9 полях:
• CHROM - последовательность, в которой был вызван вариант
• POS - координата варианта в последовательности
• ID - идентификатор варианта
• REF - нуклеотид в референсе
• ALT - нуклеотид в варианте
• QUAL - качество
• FILTER - пройденные вариантом фильтры
• INFO - дополнительная информация о варианте
• FORMAT - дополнительные поля для описания образца
Затем данный файл был проанализирован при помощи команды:
bcftools stats variants_chr4.vcf > stats_var_chr4.txt
Из файла stats_var_chr4.txt была получена информация о том, что всего было обнаружено 507 069 вариантов, из них 496 003 - однонуклеотидные полиморфизмы (SNP), 11 066 - индели.
Фильтрация вариантов
Варианты были отфильтрованы по качеству и покрытию при помощи bcftools:
bcftools filter -i'%QUAL>30 && DP>50' variants_chr4.vcf -o filt_var_chr4.vcf
Файл с фильтрованными вариантами был также проанализирован:
bcftools stats filt_var_chr4.vcf > stats_filt_var.txt
В файле filt_var_chr4.vcf осталось только 23 091 вариантов: 22 754 SNP и 337 инделей; общее количество вариантов уменьшилось почти в 22 раза, при этом доля SNP возросла в сравнении с нефильтроваными вариантами.
Изучение покрытия
Для этого с помощью программы bedtools были получены участки генома, покрытые прочтениями хотя бы один раз:
bedtools genomecov -bg -ibam correct_markd.bam > covered_1.bed
Далее приведён фрагмент файла covered_1.bed:
NC_000004.12 3255 3274 2 NC_000004.12 9998 10007 1 NC_000004.12 10007 10009 2 NC_000004.12 10009 10010 3 NC_000004.12 10010 10016 4 NC_000004.12 10016 10017 3 NC_000004.12 10017 10032 4 NC_000004.12 10032 10035 3 NC_000004.12 10035 10036 2 NC_000004.12 10036 10037 3 NC_000004.12 10037 10041 4 NC_000004.12 10041 10045 6 NC_000004.12 10045 10049 7 NC_000004.12 10049 10052 6 NC_000004.12 10052 10062 7 NC_000004.12 10062 10066 8 NC_000004.12 10066 10069 10 NC_000004.12 10069 10079 9 NC_000004.12 10079 10080 8 NC_000004.12 10080 10081 11
Файл состоит из четырёх полей: имя референса, начальная координата прочтения, конечная координата прочтения, покрытие соответственно.
Затем требовалось отобрать только участки с покрытием более 50. Для этого была использована команда awk:
awk '$4 > 50' covered_1.bed > covered_50.bed
Далее число таких участков (число заптсей в файле covered_50.bed) было подсчитано при помощи команды:
wc -l covered_50.bed
В результате получилось 1 309 173 участков.
Наконец, предстояло отделить варианты, которые попали на эти участки с покрытием более 50. Сделано это было средствами bedtools:
bedtools intersect -a filt_var_chr4.vcf -b covered_50.bed > variants_covered_50.vcf