SNP

[psilentium]

1. SRA to fastq fastqdump -- split-3
2. FastQC

[avershinina]

started 6-mat-14

fastqc /storage1/home/s_alisa/Hydra/data/SRR1032106_1.fastq && fastqc /storage1/home/s_alisa/Hydra/data/SRR1032106_2.fastq && fastqc /storage1/home/s_alisa/Hydra/data/SRR1033637_1.fastq && fastqc /storage1/home/s_alisa/Hydra/data/SRR1033637_2.fastq

[avershinina]

Всё качество ОК, кроме SRR1032106_2 Его подрезать надо, причём конкретно. Картинки тут

/storage1/home/s_alisa/Hydra/data/fastqc

[agbragin]

Выравнивание SRR1033637 закончено. Результаты в /storage1/home/s_alisa/Hydra/analysis/bowtie2 Суммарный процент выровненных ридов 69.41%.

[agbragin]

Данные с 454 (SRX000112) мы решили не брать?

[avershinina]

Решили не брать, ага.

[agbragin]

Запустил trimmomatic для SRR1032106. Вырезаем адаптеры TruSeq2, обрезаем с конца с порогом 28, выбрасываем риды короче 20. Скрипт для запуска добавлен в репозиторий.

Предлагаю все команды запуска тулов добавлять в виде sh в папку git/scripts. Потом делать: git add -A git commit -m "Description of the repository changes" git push origin master

Так мы сохраним все знания, с таким трудом выковырянные из мануалов=)

[myope]

Продолжили variant-calling pipeline на двух .sam файлах (Hydra/analyses/bowtie2) -- без отсеивания GC-богатых ридов.

samtools view -bS .sam > .bam

Запустили то же самое с опцией, которая делает .bam файл без незамапленных ридов

samtools view -F 4 -bS .sam > _mapped.bam

Сортировка .bam файлов

samtools sort _mapped.bam _mapped.sorted

MPILEUP (reference == indexed reference genome)

s_alisa@dcdell:~/Hydra/analysis/bowtie2$ samtools mpileup -ugf ../../data/hma_ref_Hydra_RP_1.0_chrUn.fa.fai ./SRR1033637_mapped.sorted.bam | bcftools view -bvcg -> raw.bcf

[agbragin]

В образце SRR1033637 контаминации не было, из данных SRR1032106 она была убрана перед выравниванием.

[psilentium]

а почему для mpileup использовался файл fai, а не fa? я решила снова запустить с fasta samtools mpileup -ugf ../../data/hma_ref_Hydra_RP_1.0_chrUn.fa ./SRR1033637_mapped.sorted.bam | bcftools view -bvcg -> raw.SRR1033637.bcf samtools mpileup -ugf ../../data/hma_ref_Hydra_RP_1.0_chrUn.fa ./SRR1032106_mapped.sorted.bam | bcftools view -bvcg -> raw.SRR1032106.bcf

[psilentium]

Convert to vcf bcftools view raw.SRR1033637.bcf| vcfutils.pl varFilter -D100 > raw.SRR1033637.flt.vcf bcftools view raw.SRR1032106.bcf| vcfutils.pl varFilter -D100 > raw.SRR1032106.flt.vcf

[avershinina]

у нас же СНПями вроде Аня занимается? А что дальше с этими vcf-ами сделать можно? Проаннотировать? Извитие, я тут уже всё забыла.

[psilentium]

vcf можно аннотировать с помощью snpeff, но непонятно, что брать в качестве референса.

[psilentium]

Нужно будет еще сделать график плотности snp

[avershinina]

В качестве референса у snpEff нет гидры, поэтому я качнула нематостеллу и на ней прогнала snpEff с двумя файлами, они сейчас перекидываются в папку /Hydra/analysis/snpEff Что это за два файла? Чем оони отличаются? С одного эффектов на гиг набралось, с другого на 80 метров.

[agbragin]

Чтобы аннотировать с использованием snpEff, нужен сделать базу данных разметки. Для этого нам нужны гены в формате GFF. Промежуточный итог относительно поиска SNP: у нас есть VCF для двух геномов (raw.SRR1032106.flt.vcf и raw.SRR1033637.flt.vcf), нужно как-то проанализировать варианты, которые в них содержатся.

Вопрос 1: почему один файл в 10 раз больше другого?)

[psilentium]

вообще и сами sra архивы отличались в 5 раз. однако я сейчас посмотрела, и меня смутил файл SRR1032106.log, потому что в нем нет результатов выравнивания, будто бы не до конца оно прошло. а в SSRR1033637.log все ок

[psilentium]

статистика и графики http://vcftools.sourceforge.net/htslib.html#stats

s_alisa@dcdell:~/Hydra/analysis/bowtie2/bcftools$ ./bcftools stats ~/Hydra/analysis/bowtie2/raw.SRR1032106.flt.vcf.gz > SRR1032106.stat.vchk

./plot-vcfstats SRR1033637.stat.vchk -p plots1/

[agbragin]

VCF-файлы скопированы в отдельную папку:

~/Hydra/analysis/vcf

Фильтр по качеству с порогом 30.

vcftools --vcf SRR1033637.vcf --minQ 30 --recode --out SRR1033637_qual_filtered.vcf > SRR1033637.vcf.log
vcftools --vcf SRR1032106.vcf --minQ 30 --recode --out SRR1032106_qual_filtered.vcf > SRR1032106.vcf.log

[agbragin]

Варианты отфильтрованы по повторам.

vcftools --vcf SRR1032106_qual_filtered.vcf --exclude-bed ../repeats/Hydra_RP_1.0_repeats_all.bed --recode --out SRR1032106_qual_filtered_masked >> SRR1032106.vcf.log

Информация о количестве отфильтрованных на каждом шаге вариантов в соответствующих лог-файлах.

[avershinina]

А получилось так, что лог minQ перезаписался логом повторов? Т. е. инфа After filtering, kept 2258558 out of a possible 3866781 Sites

• это про повторы?

[agbragin]

Нет, логи не переписывались, логи дописывались. Первый запуск - фильтр по качеству (SRR1033637 192190 out of a possible 397654 Sites, SRR1032106 3866781 out of a possible 5172032 Sites) второй - по повторам (SRR1033637 118363 out of a possible 192190 Sites, SRR1032106 2258558 out of a possible 3866781 Sites).

[avershinina]

Для отчёта, что б не забыть:

• 106: было 4.6 млн SNP + 500к indels стало 2 млн SNP + 200k indels 1 change every 344 bases, что составляет ~25% - как-то очень много ts\tv не изменилось (как такое возможно лол)
• 637 было 200к SNP + 200k indels стало 100к SNP 1 change every 6,256 bases (как такое возможно лол)

В статье: "We estimate the SNP rate to be 0.69 +/- 0.04% (or 1 SNP per 144 bp on average)"

[avershinina]

Если я правильно понимаю аутпут про гены, то: найдено ~ 90k "генов" ~ 30k из них - с длиной экзона больше чем 900 п.н. Большая часть сдлиной экзона а-ля 3-5, их наверно считать не надо.

[avershinina]

Статистика по гетерозиготности:

• 106

Type Total Homo Hetero SNP 2,073,918 1,947,442 126,476 INS 107,724 103,891 3,833 DEL 103,745 97,419 6,326 Total 2,285,412 2,148,777 136,635

• 637

Type Total Homo Hetero SNP 46337 8882 37455 INS 27162 7989 19173 DEL 45628 11684 33944 Total 119176 28604 90572

[avershinina]

В 637 гетеро ненормально больше, но он сам какой-то кривой, поэтому предлагаю ориентироваться на 106.

Что ещё надо сделать, кто помнит?