Нуклеотидный BLAST

Часть I. Анализ качества чтений
У нас был файл chr20.1.fastq. Мы произвели контроль качества чтений с помощью программы FastqC с помощью следующей команды:

fastqc chr20.1.fastq

Далее мы убрали все чтения (как в практикуме 11) качеством менее 20 и длиной менее 50 с помощью команд:

java -jar /nfs/srv/databases/ngs/suvorova/trimmomatic/trimmomatic-0.30.jar SE -phred33 chr20.1.fastq chr20.1.v2.fastq TRAILING:20

java -jar /nfs/srv/databases/ngs/suvorova/trimmomatic/trimmomatic-0.30.jar SE -phred33 chr20.1.v2.fastq chr20.1.v3.fastq MINLEN:50

Далее был снова произведён контроль качества чтений с помощью программы FastqC:

fastqc chr20.1.v3.fastq

Как можно видеть, наша обработка ридов программой trimmomatic не имела особого смысла.

Часть II и III. Картирование чтений и анализ выравнивания.
Мы использовали индексированную хромосому из практикума 11.
Затем мы строили выравнивания прочтений и референса с помощью комманды:

 hisat2 -x hisat2_index_base -U chr20.1.v3.fastq --no-softclip -S chr20pr12.sam

Мы убирали параметр --no-spliced-alignment, потому что это анализ транскриптома, и разрывы допустимы, так как в зрелой мРНК вырезаются интроны.
Далее переводим выравнивания из формата .sam в бинарный формат .bam с помощью команды

samtools view -b chr20pr12.sam >> chr20pr12.bam

Сортируем выравнивания по координате в референсе начала чтения:

 samtools sort chr20pr12.bam pr12out

Индексируем отсортированный .bam файл:

samtools index pr12out.bam

Ридов откартировано на хромосому: 2498(все чтения откартированы).Все выровнены 1 раз.

Часть IV и V. Подсчёт чтений.
Мы пользовались программой htseq-count со следующей командой:

htseq-count -f bam -s yes -i gene_id -m intersection-nonempty pr12out.bam /P/y14/term3/block4/SNP/rnaseq_reads/gencode.v19.chr_patch_hapl_scaff.annotation.gtf >> count

-f необходимо указать формат входного файла(sam или bam)
-s цепь
-i атрибут формата GFF, используемый как feature ID
-m обработка ридов, перекрывающихся более одного раза.(Берём intersection-nonempty чтобы посмотреть количество непустых перекрываний)
script.py показал нам вот такую выдачу:

ENSG00000125835.13      2601
ENSG00000251806.1       4

2605 ридов легли в 2 гена: ENSG00000125835.13(2601 рид) и ENSG00000251806.1(4 рида). Не выровняно 42 последовательности:


__no_feature    274
__ambiguous     609
__too_low_aQual 0
__not_aligned   42
__alignment_not_unique  0

Геномные браузеры