# 目录
1.Module 1 - Introduction to RNA sequencing
2.Module 2 - RNA-seq Alignment and Visualization
3.Module 3 - Expression and Differential Expression
4.Module 4 - Isoform Discovery and Alternative Expression
- Reference Guided Transcript Assembly
- de novo Transcript Assembly
- Transcript Assembly Merge
- Differential Splicing
- Splicing Visualization
5.Module 5 - De novo transcript reconstruction
6.Module 6 - Functional Annotation of Transcripts
# 3.4 Alignment Free Expression Estimation (Kallisto)
# 获取转录本的 fasta 序列
请注意，我们已经在 RNA-seq 教程的前面有了参考基因组序列的 fasta 序列。然而，Kallisto 直接作用于目标 cDNA / 转录本序列。记住，我们有 22 号染色体上基因的转录模型。这些转录模型是以 GTF 格式从 Ensembl 下载的。GTF 包含组成每个转录本的外显子的坐标描述，但不包含转录本序列本身。所以目前我们没有 Kallisto 需要的转录序列。我们有很多地方可以得到这些转录序列。
为了将 Kallisto 结果与 StringTie 的表达式结果进行比较，我们将创建一个定制的 Fasta 文件，该文件对应于用于 StringTie 分析的文本。我们如何以 Fasta 格式获得这些转录序列？
我们可以下载完整的人类 fasta 转录本数据库并只提取出 22 号染色体上的基因。或者我们可以使用来自 tophat 的名为 gtf_to_fasta 的工具从我们的 GTF 文件中生成 fasta 序列。这种方法很方便，因为它还包括控制中的 ERCC 峰值序列，允许我们为这些特征生成 Kallisto 丰度估计值。
gtf_to_fasta chr22_with_ERCC92.gtf chr22_with_ERCC92.fa chr22_ERCC92_transcripts.fa |
使用 less 查看文件 chr22_ERCC92_transcripts.fa。请注意，该文件有混乱的文字记录名称。使用下面的 hairball perl 一行代码来整理每个 fasta 序列的标题行
cat chr22_ERCC92_transcripts.fa | perl -ne 'if ($_ =~/^\>\d+\s+\w+\s+(ERCC\S+)[\+\-]/){print ">$1\n"} elsif ($_ =~ /\d+\s+(ENST\d+)/){print ">$1\n"} else {print $_}' > chr22_ERCC92_transcripts.clean.fa | |
wc -l chr22_ERCC92_transcripts*.fa | |
126025 chr22_ERCC92_transcripts.clean.fa | |
126025 chr22_ERCC92_transcripts.fa | |
252050 总用量 |
使用 less chr22_ERCC92_transcripts.clean.fa 查看生成的 “clean” 文件。使用 tail 查看文件末尾 chr22_ERCC92_transcripts.clean.fa。请注意，我们对 22 号染色体上的每个集合转录本都有一个 fasta 记录，我们对每个 ERCC spike-in 序列都有一个额外的 fasta 记录。
cat chr22_ERCC92_transcripts.clean.fa | grep ">" | perl -ne '$_ =~ s/\>//; print $_' | sort | uniq > transcript_id_list.txt |
# Build a Kallisto transcriptome index
请记住，Kallisto 并不进行比对或使用参考基因组序列。相反，它执行伪比对以确定读段与目标 (本例中为转录序列) 的兼容性。然而，与 Tophat 或 STAR 等对齐算法类似，Kallisto 需要一个索引来高效、快速地评估这种兼容性。
kallisto index --index=chr22_ERCC92_transcripts_kallisto_index ../chr22_ERCC92_transcripts.clean.fa |
# Generate abundance estimates for all samples using Kallisto
正如使用 StringTie 和 HT-Seq 所做的那样，使用 Kallisto 为每个演示样本生成转录本丰度。
kallisto quant --index=chr22_ERCC92_transcripts_kallisto_index --output-dir=UHR_Rep1_ERCC-Mix1 --threads=4 --plaintext UHR_Rep1_ERCC-Mix1_Build37-ErccTranscripts-chr22.read1.fastq.gz UHR_Rep1_ERCC-Mix1_Build37-ErccTranscripts-chr22.read2.fastq.gz | |
kallisto quant --index=chr22_ERCC92_transcripts_kallisto_index --output-dir=UHR_Rep2_ERCC-Mix1 --threads=4 --plaintext UHR_Rep2_ERCC-Mix1_Build37-ErccTranscripts-chr22.read1.fastq.gz UHR_Rep2_ERCC-Mix1_Build37-ErccTranscripts-chr22.read2.fastq.gz | |
kallisto quant --index=chr22_ERCC92_transcripts_kallisto_index --output-dir=UHR_Rep3_ERCC-Mix1 --threads=4 --plaintext UHR_Rep3_ERCC-Mix1_Build37-ErccTranscripts-chr22.read1.fastq.gz UHR_Rep3_ERCC-Mix1_Build37-ErccTranscripts-chr22.read2.fastq.gz | |
kallisto quant --index=chr22_ERCC92_transcripts_kallisto_index --output-dir=HBR_Rep1_ERCC-Mix2 --threads=4 --plaintext HBR_Rep1_ERCC-Mix2_Build37-ErccTranscripts-chr22.read1.fastq.gz HBR_Rep1_ERCC-Mix2_Build37-ErccTranscripts-chr22.read2.fastq.gz | |
kallisto quant --index=chr22_ERCC92_transcripts_kallisto_index --output-dir=HBR_Rep2_ERCC-Mix2 --threads=4 --plaintext HBR_Rep2_ERCC-Mix2_Build37-ErccTranscripts-chr22.read1.fastq.gz HBR_Rep2_ERCC-Mix2_Build37-ErccTranscripts-chr22.read2.fastq.gz | |
kallisto quant --index=chr22_ERCC92_transcripts_kallisto_index --output-dir=HBR_Rep3_ERCC-Mix2 --threads=4 --plaintext HBR_Rep3_ERCC-Mix2_Build37-ErccTranscripts-chr22.read1.fastq.gz HBR_Rep3_ERCC-Mix2_Build37-ErccTranscripts-chr22.read2.fastq.gz |
创建一个 TSV 文件，其中包含所有六个样本的 TPM 丰度估计值。
paste */abundance.tsv | cut -f 1,2,5,10,15,20,25,30 > transcript_tpms_all_samples.tsv | |
ls -1 */abundance.tsv | perl -ne 'chomp $_; if ($_ =~ /(\S+)\/abundance\.tsv/){print "\t$1"}' | perl -ne 'print "target_id\tlength$_\n"' > header.tsv | |
cat header.tsv transcript_tpms_all_samples.tsv | grep -v "tpm" > transcript_tpms_all_samples.tsv2 | |
mv transcript_tpms_all_samples.tsv2 transcript_tpms_all_samples.tsv | |
rm -f header.tsv |
head transcript_tpms_all_samples.tsv | |
tail transcript_tpms_all_samples.tsv |
