Анализ транскриптомов

Сравнение биологических реплик

Рисунок 1. Реплика 1: качество по номеру основания

Рисунок 2. Реплика 2: качество по номеру основания

---

Остальная выдача FastQC имеет незначительные различия; упоминания стоит K-мерное распределение и анализ дубликатов последовательностей:

Рисунок 3. Реплика 1: анализ K-меров

Рисунок 4. Реплика 2: анализ K-меров

---

Из этой выдачи видно, что в репликах встречается несколько отдельных сильно экспрессированных последовательностей – на это указывают отдельные острые пики на графиках.

Рисунок 5. Реплика 1: дубликаты последовательностей

Рисунок 6. Реплика 2: дубликаты последовательностей

---

Видно, что в ридах присутствует очень много повторяющихся последовательностей (особенно в реплике 1) Это весьма характерно для RNA-Seq: для того, чтобы было возможно видеть слабоэкспрессированные последовательности, приходится многократно прочитывать хорошо экспрессированные участки, что создает достаточно большое количество дубликатов чтений.

---

Картирование чтений и анализ выравнивания

С помошью программы HISAT2 риды были откартированы, а затем с помощью пакета SAMTools проанализированы. Команды:
hisat2-build ref-chr19.fasta ref-chr19 > hisat2-build.log
python3 hisat2_extract_splice_sites.py gencode.v19.chr_patch_hapl_scaff.annotation.gtf > splice_sites.tab
grep chr19 splice_sites.tab > chr19.splice.tab
hisat2 --no-softclip -x ref-chr19 -U chr19.1.fastq -S chr19.1.sam --known-splicesite-infile chr19.splice.tab --novel-splicesite-outfile chr19.1.newsplices.tab 2> hisat2.1.log 
hisat2 --no-softclip -x ref-chr19 -U chr19.2.fastq -S chr19.2.sam --known-splicesite-infile chr19.splice.tab --novel-splicesite-outfile chr19.2.newsplices.tab 2> hisat2.2.log 
samtools view chr19.1.sam -o chr19.1.bam
samtools view chr19.2.sam -o chr19.2.bam
samtools sort chr19.1.bam -o chr19.1.sorted.bam
samtools sort chr19.2.bam -o chr19.2.sorted.bam
samtools index chr19.1.sorted.bam
samtools index chr19.2.sorted.bam
samtools depth --reference ref-chr19.fasta chr19.1.sorted.bam > chr19.1.sorted.depth.tab
samtools depth --reference ref-chr19.fasta chr19.2.sorted.bam > chr19.2.sorted.depth.tab

В файлах hisat2.1.log, hisat2.2.log указано, сколько чтений картировано на хромосому. Оказалось, что:
Реплика 1:
Картированы: 34927 ридов;
Не картированы: 9 ридов;
Реплика 2:
Картированы: 24017 ридов;
Не картированы: 6 ридов;

Подсчёт чтений

С помощью набора программ BEDTools был проведён подсчёт чтений в ридах. Команды:
bedtools bamtobed -i chr19.1.sorted.bam > chr19.1.sorted.bed
bedtools bamtobed -i chr19.2.sorted.bam > chr19.2.sorted.bed
bedtools intersect -c      -a gencode.genes.bed -b chr19.1.sorted.bed > chr19.1.c.bed
bedtools intersect -c      -a gencode.genes.bed -b chr19.2.sorted.bed > chr19.2.c.bed
bedtools intersect -u      -a gencode.genes.bed -b chr19.1.sorted.bed > chr19.1.u.bed
bedtools intersect -u      -a gencode.genes.bed -b chr19.2.sorted.bed > chr19.2.u.bed
bedtools intersect -wa -wb -a gencode.genes.bed -b chr19.1.sorted.bed > chr19.1.wab.bed
bedtools intersect -wa -wb -a gencode.genes.bed -b chr19.2.sorted.bed > chr19.2.wab.bed

---

Ген	Реплики	Белок/РНК	Ридов/ген	Направление	Начало, конец (chr19)	Функция
DAPK3	2	Protein coding	85	→	3958451 → 3971121	Death-associated protein kinase 3 – induces morphological changes in apoptosis when overexpressed in mammalian cells.
MIR637	2	miRNA	3	→	3961412 → 3961510	MicroRNA 637 – is involved in post-transcriptional regulation of gene expression.
EEF2	1, 2	Protein coding	400131	→	3976054 → 3985461	Eukaryotic elongation factor 2 – encodes a member of the GTP-binding translation elongation factor family.
SNORD37	1, 2	snoRNA	906	→	3980505 → 3984570	Small nucleolar RNA – functions in the modification of other small nuclear RNAs.

---

Расширенный анализ

1. Получить из файла в выравниванием файл с чтениями в формате fastq:
   bedtools bamtofastq -i chr19.1.bam -fq chr19.1.bam.fastq

   ---

2. Получить файл с нуклеотидной последовательностью для одного из покрытых чтениями генов:
   bedtools getfasta -fi ref-chr19.fasta -bed stdin -fo EEF2.fasta <<< "`head -n 1 chr19.1.u.bed`"

   ---

3. Разбить хромосому на фрагменты по 1 млн нуклеотидов:
   printf "chr19\t0\t`infoseq -only -length -noheading -sequence ref-chr19.fasta`\n" | bedtools makewindows -b stdin -w 1000000 > chr19.windows.tab
   Длина 19-й хромосомы – 59`128`983 п.н. В результате получилось 60 интервалов.
    
4. Объединить чтения в кластеры:
   bedtools cluster -i chr19.1.sorted.bed -s > chr19.1.cluster.bed

   ---

5. Набрать из хромосомы 1000 случайных фрагментов по 200 нуклеотидов:
   printf "chr19\t`infoseq -only -length -noheading -sequence ref-chr19.fasta`\n" > gf.txt
   bedtools random -g gf.txt -n 1000 -l 200 -seed 561 > chr19.random.bed # Parameter -g seems to fail to accept info from stdin 
   rm gf.txt

   ---

6. Получить координаты 3`-области одного из покрытых чтениями генов длиной в 1000 нуклеотидов:
   

   ---

7. Получить координаты одного из покрытых чтениями генов, расширенные на 1000 нуклеотидов в обе стороны:
   printf "chr19\t`infoseq -only -length -noheading -sequence ref-chr19.fasta`\n" > gf.txt
   echo "`head -n 1 chr19.1.u.bed`" | bedtools slop -i stdin -g gf.txt -b 1000 > EEF2.slop.bed
   rm gf.txt

   ---

8. Получить координаты одного из покрытых чтениями генов, сдвинутые на 500 нуклеотидов ближе к началу хромосомы:
   printf "chr19\t`infoseq -only -length -noheading -sequence ref-chr19.fasta`\n" > gf.txt
   echo "`head -n 1 chr19.1.u.bed`" | bedtools shift -i stdin -g gf.txt -s -500 > EEF2.shift.bed
   rm gf.txt

   ---

9. Получить непересекающиеся фрагменты, соответствующие области, покрытой чтениями:
   printf "chr19\t`infoseq -only -length -noheading -sequence ref-chr19.fasta`\n" > gf.txt
   bedtools genomecov -i chr19.1.sorted.bed -g gf.txt -bg | grep -w '1$' > chr19.1.single.bed
   rm gf.txt

   ---

10. Получите файл с координатами интервалов, покрытых чтениями, с информацией о покрытии:
    printf "chr19\t`infoseq -only -length -noheading -sequence ref-chr19.fasta`\n" > gf.txt
    bedtools genomecov -i chr19.1.sorted.bed -g gf.txt -bg > chr19.1.coverage.bed
    rm gf.txt

    ---